Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 22 additions & 8 deletions activitysim/abm/models/trip_scheduling_choice.py
Original file line number Diff line number Diff line change
Expand Up @@ -279,14 +279,28 @@ def run_trip_scheduling_choice(
) in chunk.adaptive_chunked_choosers(state, indirect_tours, trace_label):
# Sort the choosers and get the schedule alternatives
choosers = choosers.sort_index()
# FIXME-EET: For explicit error term choices, we need a stable alternative ID. Currently, we use
# SCHEDULE_ID, which justs enumerates all schedule alternatives, of which there are choosers times
# alternative, in the order they are processed, which depends on if there stops on outward/return leg.
# We might want to change SCHEDULE_ID to a fixed pattern of all possible combinations of
# (outbound, main, inbound) duration for the maximum possible tour duration (max time window). For
# 30min intervals, this leads to 1225 alternatives and therefore reasonable memory-wise for random numbers.
# It looks like all that would need to change for this is the generation of the schedule alternatives and
# the lookup of choices as elements in schedule after simulation because choosers are indexed by tour_id.
# FIXME-EET: under use_explicit_error_terms, error terms here are aligned positionally, not keyed
# to a stable alternative ID: no alts_context is passed to _interaction_sample_simulate (this direct
# private call also bypasses the public wrapper's warning about that), so each tour draws one EV1
# error per alternative slot from its own tour_id-keyed channel, and the j-th draw attaches to the
# j-th row of the tour's block in `schedules`. SCHEDULE_ID plays no role in the alignment; it is a
# per-call running enumeration used only to look up the chosen row after the simulation.
#
# The per-tour row order is canonical: get_pattern_index_and_arrays sorts each tour's feasible
# windows lexicographically by (main, outbound, inbound) duration via np.unique, independent of
# chunk composition and processing order. Error terms are therefore stable across scenarios for
# every tour whose stop pattern (outbound/inbound) and duration are unchanged. They are NOT aligned
# for a tour whose duration or stop pattern changes: the lexicographic enumeration shifts, so
# position j maps to a different duration triple.
#
# To keep draws aligned across such changes too, key them to a canonical universe of schedule
# patterns -- e.g. all (outbound, inbound) duration pairs up to the maximum time window, with the
# main leg implied by the tour duration; for 30min intervals that is 1225 stable IDs, reasonable
# memory-wise for random numbers, and a duration change then keeps the error terms of unchanged
# (outbound, inbound) allocations. This would mean making SCHEDULE_ID that stable pattern ID and
# passing an alts_context (see the alt_nrs_df machinery in _interaction_sample_simulate); the
# post-simulation lookup below would then need to match on (tour_id, SCHEDULE_ID) pairs since IDs
# would repeat across tours.

schedules = generate_schedule_alternatives(choosers).sort_index()

Expand Down
4 changes: 2 additions & 2 deletions activitysim/abm/models/util/bias_logsums.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,11 @@ def maybe_bias_logsums(state: workflow.State, choices_df: pd.DataFrame, model_se
else:
logger.warning(
"Using Poisson sampling method for location choice logsum calculations. Currently the logsums results will"
+ " differ from those obtained with monte_carlo or eet sampling by a constant shift of"
+ " differ from those obtained with inverse_cdf or eet sampling by a constant shift of"
+ f" log({model_settings.SAMPLE_SIZE}) if using the common correction factor"
+ " `log(pick_count / prob)` in location choice specs. The results of the Poisson method are unbiased,"
+ " i.e., they agree with the results obtained with a full destination sample, unlike those for"
+ " monte_carlo or eet sampling."
+ " inverse_cdf or eet sampling."
)

return choices_df
4 changes: 2 additions & 2 deletions activitysim/core/configuration/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -135,12 +135,12 @@ class ComputeSettings(PydanticBase):
Sharrow settings for a component.
"""

sample_method: None | Literal["monte_carlo", "eet", "poisson"] = None
sample_method: None | Literal["inverse_cdf", "eet", "poisson"] = None
"""
Override the alternative sampling method used by `interaction_sample`.

When unset, `interaction_sample` preserves legacy behavior: it uses
`monte_carlo` when explicit error terms are off and `poisson` when they
`inverse_cdf` when explicit error terms are off and `poisson` when they
are on.
"""

Expand Down
11 changes: 6 additions & 5 deletions activitysim/core/configuration/top.py
Original file line number Diff line number Diff line change
Expand Up @@ -786,17 +786,18 @@ def _check_store_skims_in_shm(self):
Make choice from random utility model by drawing from distribution of unobserved
part of utility and taking the maximum of total utility.

Defaults to standard Monte Carlo method, i.e., calculating probabilities and then
drawing a single uniform random number to draw from cumulative probabily.
Defaults to the standard inverse-CDF (probability-based) method, i.e., calculating
probabilities and then drawing a single uniform random number against the
cumulative probability.

.. versionadded:: 1.6
"""

sample_method: None | Literal["monte_carlo", "eet", "poisson"] = None
sample_method: None | Literal["inverse_cdf", "eet", "poisson"] = None
"""
Sampling method to use in `activitysim.core.interaction_sample`.

When unset, `monte_carlo` is used when `use_explicit_error_terms` is false and
When unset, `inverse_cdf` is used when `use_explicit_error_terms` is false and
`poisson` is used when it is true.

.. versionadded:: 1.6
Expand All @@ -806,7 +807,7 @@ def _check_store_skims_in_shm(self):
"""
Whether to apply a bias of `log(sample_size)` to the Poisson sampling results.
This is a temporary workaround to align Poisson sampling results with the biased
results of the monte_carlo and eet sampling methods, such that models that were
results of the inverse_cdf and eet sampling methods, such that models that were
estimated with historical biased sampling results can be run with Poisson sampling
without needing to re-estimate the model.

Expand Down
71 changes: 47 additions & 24 deletions activitysim/core/interaction_sample.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@

DUMP = False

InteractionSampleMethod = typing.Literal["monte_carlo", "eet", "poisson"]
InteractionSampleMethod = typing.Literal["inverse_cdf", "eet", "poisson"]

# Threshold on P0, the probability that a chooser's Poisson draw comes up empty, below
# which the fallback term is dropped from the reported inclusion probabilities. Choosers
Expand Down Expand Up @@ -82,7 +82,7 @@ def resolve_sample_method(
sampling_method = state.settings.sample_method
if sampling_method is None:
sampling_method = (
"poisson" if state.settings.use_explicit_error_terms else "monte_carlo"
"poisson" if state.settings.use_explicit_error_terms else "inverse_cdf"
)
if sampling_method not in typing.get_args(InteractionSampleMethod):
raise ValueError(
Expand Down Expand Up @@ -134,6 +134,10 @@ def _poisson_fallback_positions(
alternatives for each row of `probs_values` (all of them if there are fewer
than `sample_size` alternatives), with ties broken by column position.

A row with fewer than `sample_size` positive-probability alternatives gets
zero-probability positions in its trailing columns; the caller drops those
pairs so that an unavailable alternative can never enter the choice set.

This is deliberately *deterministic* and consumes no random numbers, so every
chooser row advances its RNG channel by exactly the same amount whether or not
the fallback fires. That keeps random number streams aligned across scenarios,
Expand Down Expand Up @@ -166,7 +170,7 @@ def make_sample_choices_eet(

Each chooser receives `sample_size` EV1 draw sets and the argmax-over-utility
winner is recorded per draw, so duplicates are possible (same with-replacement
semantics as the Monte Carlo sampling path).
semantics as the inverse-CDF sampling path).

`utilities` drives the Gumbel argmax. `probs` (the MNL choice probabilities
computed from the same utilities by the caller) supplies the `prob` column
Expand Down Expand Up @@ -225,7 +229,7 @@ def make_sample_choices_poisson(

where `p_i` is the chooser's MNL choice probability for alternative `i`. That is the
probability the alternative would have been drawn at least once across `sample_size`
Monte Carlo draws, which is what makes Poisson sampling interchangeable with the other
inverse-CDF draws, which is what makes Poisson sampling interchangeable with the other
sampling methods. `pick_count` is 1 by definition (the draw is a yes/no per
alternative), so the standard sampling correction factor is recoverable in the usual
way as `np.log(df.pick_count / df.prob)`.
Expand All @@ -238,8 +242,11 @@ def make_sample_choices_poisson(
Since the probabilities sum to one and 1 - p <= exp(-p), this is bounded above by
exp(-sample_size): very small at the sample sizes these models use (~1e-13
for `sample_size=30`), but not negligible at small sample sizes or for a chooser whose
probability mass is spread thinly. Those choosers fall back to the `sample_size`
highest-probability alternatives (see `_poisson_fallback_positions`). The fallback is
probability mass is spread thinly. Those choosers fall back to their
`min(sample_size, n_available)` highest-probability available alternatives, where
`n_available` counts the alternatives with non-zero probability (see
`_poisson_fallback_positions`) -- an unavailable alternative can never enter the
choice set through either branch. The fallback is
deterministic and draws no random numbers, so every chooser advances its RNG channel by
exactly the same amount whether or not it fires -- unlike a retry scheme, this cannot
desynchronise random number streams between scenarios.
Expand Down Expand Up @@ -293,9 +300,9 @@ def make_sample_choices_poisson(
if n_empty > 0:
logger.warning(
f"Poisson sampling drew an empty choice set for {n_empty} of {len(probs)} "
f"chooser(s) in {trace_label}; falling back to the "
f"{min(sample_size, probs_values.shape[1])} highest-probability alternatives "
f"for those choosers. Highest empty-sample probability was "
f"chooser(s) in {trace_label}; falling back to (at most) the "
f"{min(sample_size, probs_values.shape[1])} highest-probability available "
f"alternatives for those choosers. Highest empty-sample probability was "
f"{empty_sample_probs[empty_rows].max():.2g} against a requested sample size "
f"of {sample_size} and a mean expected sample size of "
f"{inclusion_probs[empty_rows].sum(axis=1).mean():.1f}."
Expand All @@ -312,6 +319,17 @@ def make_sample_choices_poisson(
)
row_positions = np.repeat(fallback_rows, fallback_cols.shape[1])
col_positions = fallback_cols.reshape(-1)

# drop zero-probability pairs: a chooser with fewer available (p > 0)
# alternatives than the fallback window would otherwise get it padded with
# unavailable alternatives, which would enter the final choice set carrying a
# large positive correction term log(1/P0). The fallback set remains a
# deterministic function of the probabilities, so the closed form for the
# reported prob is unchanged.
available = probs_values[row_positions, col_positions] > 0.0
row_positions = row_positions[available]
col_positions = col_positions[available]

inclusion_probs[row_positions, col_positions] += empty_sample_probs[
row_positions
]
Expand Down Expand Up @@ -763,10 +781,10 @@ def _interaction_sample(

sampling_method = resolve_sample_method(state, compute_settings)

# Estimation requires MC sampling and MC choice for now
if estimation.manager.enabled and sampling_method != "monte_carlo":
# Estimation requires inverse-CDF sampling and choice for now
if estimation.manager.enabled and sampling_method != "inverse_cdf":
raise ValueError(
f"{trace_label}: estimation requires monte_carlo sampling and choice. Set sample_method='monte_carlo'"
f"{trace_label}: estimation requires inverse_cdf sampling and choice. Set sample_method='inverse_cdf'"
+ " (or leave it unset) and use_explicit_error_terms=False for estimation runs."
)

Expand Down Expand Up @@ -817,7 +835,7 @@ def _interaction_sample(
column_labels=["alternative", "probability"],
)

if sampling_method == "monte_carlo":
if sampling_method == "inverse_cdf":
del utilities
chunk_sizer.log_df(trace_label, "utilities", None)

Expand Down Expand Up @@ -872,8 +890,8 @@ def _interaction_sample(
del probs
chunk_sizer.log_df(trace_label, "probs", None)
else:
# eet and poisson: optionally trim choosers with all-zero probs. The MC
# path handles this inside make_sample_choices
# eet and poisson: optionally trim choosers with all-zero probs. The
# inverse-CDF path handles this inside make_sample_choices
if allow_zero_probs:
non_zero = probs.sum(axis=1) != 0
if not non_zero.any():
Expand Down Expand Up @@ -1075,20 +1093,25 @@ def interaction_sample(
sampling_method = resolve_sample_method(state, compute_settings)
logger.debug(f" interaction_sample sample method = {sampling_method}")

if sampling_method == "monte_carlo":
# The MC sampling path (make_sample_choices) does not consume stable_alt_positions
# or n_total_alts. Null them out so callers that conservatively pass values along
# don't accidentally rely on them under MC sampling.
if sampling_method == "inverse_cdf":
# The inverse-CDF sampling path (make_sample_choices) does not consume
# stable_alt_positions or n_total_alts. Null them out so callers that
# conservatively pass values along don't accidentally rely on them under
# inverse-CDF sampling.
stable_alt_positions = None
n_total_alts = None

# FIXME - legacy logic - not sure this is needed or even correct?
if sampling_method != "poisson":
# legacy clamp for the with-replacement methods; statistically harmless because
# the omitted log(sample_size) term in the correction is constant per chooser
sample_size = min(sample_size, len(alternatives.index))
# with poisson sampling, definitely don't want to reduce sample size - it's not a sample size but a number
# of theoretical draws. Another options would be to disable sampling if # alts < sample size to ensure
# all are included (but this wouldn't behave well if there were land use changes in the project case which
# switched regimes)
# with poisson sampling the sample size must not be clamped: it is not a count of
# draws but the rate parameter of the inclusion probabilities. When a chooser has
# fewer available alternatives than sample_size, its inclusion probabilities
# saturate towards 1 and the final choice approaches exact MNL over its full
# availability set, which is the desired behavior. (Disabling sampling entirely
# in that regime would behave badly if a project-case land use change switched
# regimes.)

logger.debug(f" interaction_sample sample size = {sample_size}")

Expand Down
15 changes: 7 additions & 8 deletions activitysim/core/interaction_sample_simulate.py
Original file line number Diff line number Diff line change
Expand Up @@ -91,8 +91,7 @@ def _interaction_sample_simulate(

choices : pandas.Series
A series where index should match the index of the choosers DataFrame
and values will match the index of the alternatives DataFrame -
choices are simulated in the standard Monte Carlo fashion
and values will match the index of the alternatives DataFrame

if want_logsums is True:

Expand Down Expand Up @@ -546,8 +545,7 @@ def interaction_sample_simulate(

choices : pandas.Series
A series where index should match the index of the choosers DataFrame
and values will match the index of the alternatives DataFrame -
choices are simulated in the standard Monte Carlo fashion
and values will match the index of the alternatives DataFrame

if want_logsums is True:

Expand All @@ -565,10 +563,11 @@ def interaction_sample_simulate(
# are NOT guaranteed to be consistent across scenarios that differ in alternative
# availability. We cannot make this a hard error today because two production callers
# rely on the warning-only fallback:
# - trip_scheduling_choice: SCHEDULE_ID is a per-call enumeration that depends on
# chunk composition and tour duration distribution (see FIXME in
# trip_scheduling_choice.py:282-289 for the proposed redesign that would key
# SCHEDULE_ID to a fixed (OB, MAIN, IB) duration tuple).
# - trip_scheduling_choice: draws align positionally to each tour's canonical
# schedule enumeration, which is stable while a tour's stop pattern and duration
# are unchanged but shifts when either changes (see the FIXME-EET in
# trip_scheduling_choice.py for the proposed stable-id redesign; note it calls
# _interaction_sample_simulate directly and so does not pass through this warning).
# - tour_od_choice: OD id is a string concatenation `f"{orig}_{dest}"`; a stable
# integer universe would be O(n_zones^2) error terms per chooser, which is too
# large to allocate.
Expand Down
6 changes: 2 additions & 4 deletions activitysim/core/interaction_simulate.py
Original file line number Diff line number Diff line change
Expand Up @@ -704,8 +704,7 @@ def _interaction_simulate(
-------
ret : pandas.Series
A series where index should match the index of the choosers DataFrame
and values will match the index of the alternatives DataFrame -
choices are simulated in the standard Monte Carlo fashion
and values will match the index of the alternatives DataFrame
"""

trace_label = tracing.extend_trace_label(trace_label, "interaction_simulate")
Expand Down Expand Up @@ -1031,8 +1030,7 @@ def interaction_simulate(
-------
choices : pandas.Series
A series where index should match the index of the choosers DataFrame
and values will match the index of the alternatives DataFrame -
choices are simulated in the standard Monte Carlo fashion
and values will match the index of the alternatives DataFrame
"""

trace_label = tracing.extend_trace_label(trace_label, "interaction_simulate")
Expand Down
13 changes: 3 additions & 10 deletions activitysim/core/logit.py
Original file line number Diff line number Diff line change
Expand Up @@ -584,7 +584,9 @@ def make_choices_explicit_error_term_nl(
pandas.Series
Choice indices aligned to `alt_utilities` columns.
"""
# TODO assert alts_context and alt_nrs_df are both None - no sampling from nested models for now.
assert (
alts_context is None and alt_nrs_df is None
), f"{trace_label} - Sampling from nested models is not implemented, do not pass alts_context or alt_nrs_df."

utilities_incl_unobs = sample_nested_logit_exact_leaf_error_terms(
state,
Expand Down Expand Up @@ -688,15 +690,6 @@ def make_choices_utility_based(
rands : pandas.Series
A series of 0s for compatibility with make_choices. For EET, we do not have per-row random numbers.

Notes
-----
An argmax always returns a position, so a chooser with no available alternative gets a
choice here rather than an error: with every utility at `UTIL_UNAVAILABLE` the alternatives
are tied and the error terms decide, and with every utility at `-inf` the first column wins.
The Monte Carlo path does not go quiet in that situation -- `make_choices` reports it unless
`allow_bad_probs` is set -- so this function reports it too. Most callers reach here having
already run `validate_utils`, which makes the same check; the duplication is deliberate and
mirrors `make_choices` re-checking what `utils_to_probs` has already looked at.
"""
trace_label = tracing.extend_trace_label(trace_label, "make_choices_utility_based")

Expand Down
Loading
Loading