diff --git a/activitysim/abm/models/trip_scheduling_choice.py b/activitysim/abm/models/trip_scheduling_choice.py index a5e17eb41..63f19f5a1 100644 --- a/activitysim/abm/models/trip_scheduling_choice.py +++ b/activitysim/abm/models/trip_scheduling_choice.py @@ -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() diff --git a/activitysim/abm/models/util/bias_logsums.py b/activitysim/abm/models/util/bias_logsums.py index 26793dad6..51b0ae250 100644 --- a/activitysim/abm/models/util/bias_logsums.py +++ b/activitysim/abm/models/util/bias_logsums.py @@ -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 diff --git a/activitysim/core/configuration/base.py b/activitysim/core/configuration/base.py index b0cf506fb..f4258fb94 100644 --- a/activitysim/core/configuration/base.py +++ b/activitysim/core/configuration/base.py @@ -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. """ diff --git a/activitysim/core/configuration/top.py b/activitysim/core/configuration/top.py index f74e2a30e..6b1c8a6c5 100644 --- a/activitysim/core/configuration/top.py +++ b/activitysim/core/configuration/top.py @@ -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 @@ -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. diff --git a/activitysim/core/interaction_sample.py b/activitysim/core/interaction_sample.py index 84acadefb..58c522d7e 100644 --- a/activitysim/core/interaction_sample.py +++ b/activitysim/core/interaction_sample.py @@ -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 @@ -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( @@ -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, @@ -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 @@ -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)`. @@ -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. @@ -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}." @@ -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 ] @@ -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." ) @@ -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) @@ -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(): @@ -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}") diff --git a/activitysim/core/interaction_sample_simulate.py b/activitysim/core/interaction_sample_simulate.py index 2083a6530..1b56f9d24 100644 --- a/activitysim/core/interaction_sample_simulate.py +++ b/activitysim/core/interaction_sample_simulate.py @@ -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: @@ -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: @@ -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. diff --git a/activitysim/core/interaction_simulate.py b/activitysim/core/interaction_simulate.py index bb3213498..76c2c6acd 100644 --- a/activitysim/core/interaction_simulate.py +++ b/activitysim/core/interaction_simulate.py @@ -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") @@ -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") diff --git a/activitysim/core/logit.py b/activitysim/core/logit.py index 9c3f668d2..d9ef5d8f4 100644 --- a/activitysim/core/logit.py +++ b/activitysim/core/logit.py @@ -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, @@ -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") diff --git a/activitysim/core/test/test_interaction_sample.py b/activitysim/core/test/test_interaction_sample.py index 2d929053d..e6f1af710 100644 --- a/activitysim/core/test/test_interaction_sample.py +++ b/activitysim/core/test/test_interaction_sample.py @@ -23,7 +23,7 @@ def state() -> workflow.State: def test_interaction_sample_ignores_stable_positions_without_global_eet( state, monkeypatch ): - # Do not support stable alt positions or tracking total alts when running with MC sampling + # Do not support stable alt positions or tracking total alts when running with inverse-CDF sampling # to not introduce any additional changes while adding eet simulation support to ensure no # regressions. We can add these features later if desired. captured = {} @@ -139,7 +139,7 @@ def test_interaction_sample_parity(state): index=pd.Index(["chooser_attr * alt_attr"], name="Expression"), ) - # Run Monte Carlo with replacement. + # Run inverse-CDF sampling with replacement. state.settings.use_explicit_error_terms = False state.rng().set_base_seed(42) state.rng().add_channel("person_id", choosers) @@ -198,14 +198,14 @@ def test_interaction_sample_parity(state): assert choices_eet["alt_id"].isin(alternatives.index).all() shares = { - "monte_carlo": _weighted_shares(choices_mnl), + "inverse_cdf": _weighted_shares(choices_mnl), "poisson": _weighted_shares(choices_poisson), "eet": _weighted_shares(choices_eet), } for left, right in [ - ("monte_carlo", "poisson"), - ("monte_carlo", "eet"), + ("inverse_cdf", "poisson"), + ("inverse_cdf", "eet"), ("poisson", "eet"), ]: all_alts = set(shares[left].index) | set(shares[right].index) @@ -428,7 +428,7 @@ def _shares_for_sample( return choices, _weighted_shares(choices) -def test_interaction_sample_eet_sampling_under_mc_simulation(state): +def test_interaction_sample_eet_sampling_under_inverse_cdf_simulation(state): # use_eet=False + sample_method="eet" was silently ignored before the # sampling/simulation decoupling. The dispatch now keys on sampling_method # only, so this combo must produce shares that match use_eet=True + eet. @@ -483,9 +483,9 @@ def test_interaction_sample_eet_sampling_under_mc_simulation(state): ) -def test_interaction_sample_poisson_sampling_under_mc_simulation(state): - # use_eet=False + sample_method="poisson" used to silently fall through to MC - # sampling and then have pick_count forced to 1, corrupting results. After +def test_interaction_sample_poisson_sampling_under_inverse_cdf_simulation(state): + # use_eet=False + sample_method="poisson" used to silently fall through to + # inverse-CDF sampling and then have pick_count forced to 1, corrupting results. After # decoupling, the combo must run the Poisson path and match use_eet=True + poisson. num_choosers = 100_000 num_alts = 100 @@ -530,7 +530,7 @@ def test_interaction_sample_poisson_sampling_under_mc_simulation(state): # Poisson contract: pick_count must be uniformly 1 assert (choices_mc_sim["pick_count"] == 1).all(), ( - "Poisson sampling under MC simulation must produce pick_count=1; got " + "Poisson sampling under inverse-CDF simulation must produce pick_count=1; got " f"{choices_mc_sim['pick_count'].value_counts().to_dict()}" ) @@ -544,7 +544,7 @@ def test_interaction_sample_poisson_sampling_under_mc_simulation(state): ) -def test_interaction_sample_mc_sampling_under_eet_simulation(state): +def test_interaction_sample_inverse_cdf_sampling_under_eet_simulation(state): num_choosers = 100_000 num_alts = 100 sample_size = 10 @@ -570,7 +570,7 @@ def test_interaction_sample_mc_sampling_under_eet_simulation(state): spec, sample_size, use_eet=False, - sample_method="monte_carlo", + sample_method="inverse_cdf", seed=42, step_name="test_mc_under_mc_sim", ) @@ -581,7 +581,7 @@ def test_interaction_sample_mc_sampling_under_eet_simulation(state): spec, sample_size, use_eet=True, - sample_method="monte_carlo", + sample_method="inverse_cdf", seed=42, step_name="test_mc_under_eet_sim", ) @@ -590,7 +590,7 @@ def test_interaction_sample_mc_sampling_under_eet_simulation(state): for alt in all_alts: diff = abs(shares_mc_sim.get(alt, 0.0) - shares_eet_sim.get(alt, 0.0)) assert diff < 0.01, ( - f"MC sampling shares should not depend on simulation mode at alt {alt}: " + f"Inverse-CDF sampling shares should not depend on simulation mode at alt {alt}: " f"mc_sim={shares_mc_sim.get(alt, 0.0):.4f}, " f"eet_sim={shares_eet_sim.get(alt, 0.0):.4f}, diff={diff:.4f}" ) @@ -657,8 +657,9 @@ def _reference_poisson_sampled_values(probs_np, draws, sample_size): An alternative ends up in the choice set if its Bernoulli draw succeeded, or if the chooser drew nothing at all and the alternative is one of the `sample_size` most - likely. Those events are disjoint, so the probability of an alternative being in the - returned set is `q_i + P0 * 1{i in fallback set}` for every chooser and both branches. + likely *available* (p > 0) alternatives. Those events are disjoint, so the + probability of an alternative being in the returned set is + `q_i + P0 * 1{i in fallback set}` for every chooser and both branches. Returns the sparse chooser-by-alternative array of reported probabilities, with np.nan for alternatives that are not in the choice set. @@ -673,6 +674,8 @@ def _reference_poisson_sampled_values(probs_np, draws, sample_size): k = min(sample_size, probs_np.shape[1]) top_k = np.argsort(-probs_np, axis=1, kind="stable")[:, :k] np.put_along_axis(in_fallback, top_k, True, axis=1) + # unavailable alternatives never enter the choice set + in_fallback &= probs_np > 0 # the implementation skips the P0 term where it cannot matter; mirror that here so # the comparison stays exact (see POISSON_EMPTY_SAMPLE_TOLERANCE) @@ -847,6 +850,56 @@ def test_make_sample_choices_poisson_consumes_no_extra_randoms_on_empty_draw(): pd.testing.assert_frame_equal(choices_df, expected) +def test_make_sample_choices_poisson_fallback_excludes_unavailable_alternatives(): + # a chooser with fewer available (p > 0) alternatives than the fallback window must + # not have its fallback set padded with unavailable alternatives: those would enter + # the final choice set carrying a large positive correction term log(1/P0) + probs = pd.DataFrame( + [[0.60, 0.40, 0.00, 0.00]], + index=pd.Index([11], name="person_id"), + columns=np.arange(4), + ) + sample_size = 3 + alternatives = pd.DataFrame(index=pd.Index([100, 300, 700, 900], name="alt_id")) + # both available alternatives fail their inclusion draw, forcing the fallback + fail_draw = np.array([[0.99, 0.99, 0.99, 0.99]], dtype=np.float64) + state = _DummyState(_SequentialDummyRng([fail_draw])) + + choices_df = interaction_sample.make_sample_choices_poisson( + chunk_sizer=_DummyChunkSizer(), + probs=probs, + alternatives=alternatives, + sample_size=sample_size, + alt_col_name="alt_id", + state=state, + trace_label="test_make_sample_choices_poisson_fallback_excludes_unavailable_alternatives", + ) + + # only the two available alternatives are returned, each reported at q_i + P0, + # even though the fallback window min(sample_size, n_alts) = 3 is wider + inclusion_probs = 1 - np.power(1 - probs.to_numpy(), sample_size) + empty_sample_prob = np.prod(1 - inclusion_probs, axis=1)[0] + expected = pd.DataFrame( + { + "person_id": [11, 11], + "prob": [ + inclusion_probs[0, 0] + empty_sample_prob, + inclusion_probs[0, 1] + empty_sample_prob, + ], + "alt_id": [100, 300], + } + ) + pd.testing.assert_frame_equal(choices_df, expected) + + # the reference implementation agrees + pd.testing.assert_frame_equal( + choices_df, + _reference_poisson_choices_df( + probs, fail_draw, sample_size, alternatives, "alt_id" + ), + ) + + def test_make_sample_choices_poisson_reported_prob_is_total_inclusion_probability(): # Monte Carlo check that the reported `prob` really is the probability of the # alternative ending up in the choice set, counting both the Bernoulli draw and the diff --git a/activitysim/core/test/test_logit.py b/activitysim/core/test/test_logit.py index 3cee651d2..b98d225d2 100644 --- a/activitysim/core/test/test_logit.py +++ b/activitysim/core/test/test_logit.py @@ -683,7 +683,7 @@ def test_make_choices_vs_eet_same_distribution(): utils = pd.DataFrame([utils_values] * n_draws, columns=columns) - # Probability-based (Monte Carlo) path — independent RNG + # Probability-based (inverse-CDF) path — independent RNG mc_rng = np.random.default_rng(42) class MCDummyRNG: diff --git a/docs/conf.py b/docs/conf.py index c148e6b6b..590c976a4 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -53,7 +53,7 @@ "dev-guide/_generated2/*", "users-guide/_generated/*", ] -myst_enable_extensions = ["colon_fence"] +myst_enable_extensions = ["colon_fence", "dollarmath"] myst_heading_anchors = 3 nb_merge_streams = True numpydoc_show_class_members = False diff --git a/docs/dev-guide/explicit-error-terms.md b/docs/dev-guide/explicit-error-terms.md index 5211e9fbd..a55bacc89 100644 --- a/docs/dev-guide/explicit-error-terms.md +++ b/docs/dev-guide/explicit-error-terms.md @@ -54,32 +54,32 @@ to draw error terms of all fundamental alternatives. For EET to reduce simulation noise, it is important that alternatives of a choice situation keep the same unobserved error term in different scenario runs. If unchanged alternatives keep the same unobserved draws, changes to choices between scenarios can only happen when -the observed utility of an alternative increases. This is not the case for the Monte Carlo -simulation method, where the draws are based on probabilities, which necessarily change for -all alternatives if any observed utility changes. This combined with sensitivity to small -differences in the final CDF draw when comparing nearby scenarios means that EET removes -noise from scenario comparisons. - -Note that the both MC and EET are simulating the same model, so individual runs with identical -inputs but varying global seed will lead to the same statistical results for individual -output metrics. EET's properties become apparent when comparing two model runs with different -inputs. Because error terms are aligned, the variance of the estimator of the indicator, e.g., -mode choice shift or VMT difference, is reduced. In other words, difference metrics are more +the observed utility of an alternative increases. This is not the case for the inverse-CDF +(probability-based) simulation method, where the draws are based on probabilities, which +necessarily change for all alternatives if any observed utility changes. This combined with +sensitivity to small differences in the final CDF draw when comparing nearby scenarios means +that EET removes noise from scenario comparisons. + +Note that both inverse-CDF and EET simulation are simulating the same model, so individual runs +with identical inputs but varying global seed will lead to the same statistical results for +individual output metrics. EET's properties become apparent when comparing two model runs with +different inputs. Because error terms are aligned, the variance of the estimator of the indicator, +e.g., mode choice shift or VMT difference, is reduced. In other words, difference metrics are more precise estimators under EET. In mathematical terms, for any two metrics $X$ (baseline) and $Y$ (scenario), the variance of the difference $X - Y$ is -$$\text{Var}(X - Y) = \text{Var}(X) + \text{Var}(Y) - 2,\text{Cov}(X, Y)$$ +$$\text{Var}(X - Y) = \text{Var}(X) + \text{Var}(Y) - 2\,\text{Cov}(X, Y)$$ EET deliberately drives $\text{Cov}(X, Y)$ up by aligning error terms, so $\text{Var}(X-Y)$ collapses even though $\text{Var}(X)$ and $\text{Var}(Y)$ individually are unchanged. -In practice, models are often run once for each scenario. EET is still usefull because the +In practice, models are often run once for each scenario. EET is still useful because the lower the noise of the estimator, the higher the chance that a single run is representative. -In other words, the noise level of comparison metrics is lower. Additionally, under MC small -but real benefits can show up as negative in a single run. Under EET, the sign of the effect -is far more trustworthy. +In other words, the noise level of comparison metrics is lower. Additionally, under inverse-CDF +simulation small but real benefits can show up as negative in a single run. Under EET, the sign +of the effect is far more trustworthy. Independent of any statistical argument, under EET, choice changes between two runs are attributable to utility changes which can be helpful for model development, sensitivity @@ -123,11 +123,11 @@ individual comparison runs. ### Runtime and memory usage EET draws one error term per chooser and alternative, which requires many more random numbers -than MC's one per chooser. For models with many alternatives, this can lead to a large amount -of random numbers being calculated. The implementation of EET avoids materialization of large -chooser-alternative arrays of error terms in memory so that the memory usage is in line with MC -simulation. -Regarding runtimes, EET with default settings currently carries a runtime penalty of about 3-10% +than the inverse-CDF method's one per chooser. For models with many alternatives, this can lead +to a large amount of random numbers being calculated. The implementation of EET avoids +materialization of large chooser-alternative arrays of error terms in memory so that the memory +usage is in line with inverse-CDF simulation. +Regarding runtimes, EET with default settings currently carries a runtime penalty of a few percent per demand model run. However, when run in combination with an assignment model the overall system can converge faster and this can reduce the overall model runtime penalty. @@ -167,12 +167,12 @@ do not have a corresponding EET implementation because there are no utilities to ### Unavailable choices utility convention For EET, only utility differences matter, and therefore the outcome for two utilities that are -very small, say -10000 and -10001, is identical to the outcome for 0 and 1. For MC, utilities -have to be exponentiated and therefore floating point precision dictates the smallest and largest -utility that can be used in practice. ActivitySim models historically often use a utility of --999 to make alternatives practically unavailable. That value is below the utility threshold -used in the probability-based path, which is about -691 because ActivitySim clips -exponentiated utilities at 1e-300. To keep behavior consistent, EET treats alternatives with +very small, say -10001 and -10000, is identical to the outcome for 0 and 1. For inverse-CDF +simulation, utilities have to be exponentiated and therefore floating point precision dictates +the smallest and largest utility that can be used in practice. ActivitySim models historically +often use a utility of -999 to make alternatives practically unavailable. That value is below +the utility threshold used in the probability-based path, which is about -691 because ActivitySim +clips exponentiated utilities at 1e-300. To keep behavior consistent, EET treats alternatives with utilities at or below that threshold as unavailable; see `activitysim.core.logit.validate_utils`. ### Normalization diff --git a/docs/dev-guide/sampling-methods.md b/docs/dev-guide/sampling-methods.md index 700c15c30..b98a7e85b 100644 --- a/docs/dev-guide/sampling-methods.md +++ b/docs/dev-guide/sampling-methods.md @@ -32,7 +32,7 @@ sampling utilities, while the corresponding final-choice specs in ## Available Methods -- `monte_carlo`: importance sampling with replacement using probabilities and uniform draws +- `inverse_cdf`: importance sampling with replacement using probabilities and uniform draws - `eet`: importance sampling with replacement using explicit error-term draws - `poisson`: importance sampling via independent Poisson inclusion sampling based on probabilities @@ -41,7 +41,7 @@ sampling utilities, while the corresponding final-choice specs in At the top level, `sample_method` may be set in `settings.yaml`. When it is omitted, ActivitySim preserves the intended default behavior: -- if `use_explicit_error_terms` is `False`, `interaction_sample` defaults to `monte_carlo` +- if `use_explicit_error_terms` is `False`, `interaction_sample` defaults to `inverse_cdf` - if `use_explicit_error_terms` is `True`, `interaction_sample` defaults to `poisson` Individual models may override this default through: @@ -64,20 +64,21 @@ The sampled-choice workflow is: This is the standard sample-of-alternatives pattern: the sampling stage uses an approximation, and the final stage corrects for it. -### Monte Carlo and EET-with-replacement +### Inverse-CDF and EET-with-replacement -The `monte_carlo` and `eet` sampling methods both draw alternatives with replacement. As a result, +The `inverse_cdf` and `eet` sampling methods both draw alternatives with replacement. As a result, duplicates are possible within a chooser's sampled set, and sampled shares track repeated-draw MNL behavior closely. The difference between them is how each draw is made: -- `monte_carlo` draws from analytical probabilities using uniform random numbers +- `inverse_cdf` draws from analytical probabilities using uniform random numbers against the + cumulative distribution - `eet` draws explicit EV1 error terms and chooses the utility-plus-error argmax `eet` freezes the error terms for each chooser-alternative pair across repeated draws, so that unchanged alternatives can keep the same unobserved draws, which can greatly reduce -scenario-to-scenario sampling noise compared to `monte_carlo`. However, `eet` is more expensive to +scenario-to-scenario sampling noise compared to `inverse_cdf`. However, `eet` is more expensive to run because it requires many more random draws and more complex logic to avoid materializing large chooser-alternative arrays of error terms in memory. @@ -87,14 +88,24 @@ chooser-alternative arrays of error terms in memory. pair is sampled independently with inclusion probability $1 - (1 - p)^s$, where $p$ is the original choice probability and $s$ is the configured sample size. A single inclusion draw is made for each alternative. This is much cheaper than repeated draws for -`eet`, and unlike ``monte_carlo``, it can still benefit from stable alignment of random draws to -alternatives, so it can provide improved noise reduction compared to `monte_carlo` without the full +`eet`, and unlike ``inverse_cdf``, it can still benefit from stable alignment of random draws to +alternatives, so it can provide improved noise reduction compared to `inverse_cdf` without the full cost of `eet` and therefore it is the default when running with explicit error terms, see {ref}`explicit-error-terms-dev`. - +numerical noise. The interaction-sample tests document this explicitly. + +Under `poisson`, the configured sample size $s$ is a rate parameter rather than a count of draws, +and it is deliberately not clamped to the number of alternatives (`inverse_cdf` and `eet` clamp it, +which is statistically harmless for with-replacement draws because the omitted $\log s$ term in the +correction is constant per chooser). When a chooser has fewer available alternatives than $s$, its +probability mass is concentrated and the inclusion probabilities saturate towards 1: the chooser +receives essentially its whole availability set, each alternative with a correction term near +$\log(1/1) = 0$, and the final choice approaches exact MNL over the true availability set. No +special-casing is needed for such choosers; their expected sample size $\sum_i q_i$ is simply +smaller than $s$. A chooser can occasionally receive no sampled alternatives under Poisson sampling, because each alternative is tested independently. The probability of this happening for a given chooser is @@ -107,7 +118,9 @@ Because the probabilities sum to one and $1 - p \le e^{-p}$, this is bounded abo regardless of how the probabilities are distributed. It is therefore negligible at the sample sizes these models use (at most $10^{-13}$ for a sample size of 30), but not negligible at small sample sizes, or for a chooser whose probability mass is spread very thinly. If it happens, that chooser -falls back to its $s$ highest-probability alternatives. +falls back to its $\min(s, n)$ highest-probability *available* alternatives, where $n$ is the +number of alternatives with non-zero probability. Zero-probability alternatives are never included, +so an unavailable alternative cannot enter the choice set through either branch. The fallback is deliberately deterministic and draws no random numbers, so every chooser advances its random number channel by exactly the same amount whether or not the fallback fires. A retry or @@ -137,7 +150,7 @@ positive, and the fallback can only add inclusion mass on top of that, never rem - `prob` - `pick_count` -For `monte_carlo` and `eet`, `pick_count` is the number of times the alternative was selected in +For `inverse_cdf` and `eet`, `pick_count` is the number of times the alternative was selected in the repeated with-replacement draws. For `poisson`, `pick_count` is always `1`, because an alternative is either included or not included. For all methods, `prob` is the quantity used in the correction term, but it means different things for different methods. ActivitySim's final @@ -149,7 +162,7 @@ np.log(df.pick_count/df.prob) This is the sample-of-alternatives correction factor used in the final choice model. -For `monte_carlo` and `eet`, `prob` is the one-draw sampling probability implied by the +For `inverse_cdf` and `eet`, `prob` is the one-draw sampling probability implied by the approximate sampling utility, and `pick_count` is the number of times that alternative appeared in the repeated sample. McFadden's utility correction term for repeated with-replacement sampling is `log(pick_count / (sample_size * prob)) = log(pick_count / prob) - log(sample_size)`. ActivitySim @@ -158,7 +171,7 @@ chooser and therefore does not affect choice probabilities. For `poisson`, `prob` is the inclusion probability of the alternative in the sampled set, not the one-draw choice probability. Specifically, if the original approximate choice probability is $p$ -and the configured sample size is $s$, then the inclusion probably of the Bernoulli trial is +and the configured sample size is $s$, then the inclusion probability of the Bernoulli trial is $$ q_i = 1 - (1 - p_i)^s @@ -185,7 +198,7 @@ not depend on how many times a given chooser was redrawn. Ranking the probabilities to find the fallback set costs about as much as the Bernoulli draw itself, so the implementation evaluates the fallback term only for choosers whose $P_0$ exceeds -`POISSON_EMPTY_SAMPLE_TOLERANCE`, which is set to $1e-12$, plus every chooser that actually drew +`POISSON_EMPTY_SAMPLE_TOLERANCE`, which is set to $10^{-12}$, plus every chooser that actually drew nothing. Since $P_0 \le e^{-s}$, this branch is never evaluated above a sample size of 27. Dropping the term understates `prob` by $P_0$, so the relative error on the correction is $P_0 / q_i$, which is only large for an alternative whose own inclusion probability is far below $P_0$. But such an @@ -205,25 +218,25 @@ than for the with-replacement methods. Runtime and noise characteristics differ across methods. -- `monte_carlo` is the fastest method. It draws one uniform random number per repeated sample for +- `inverse_cdf` is the fastest method. It draws one uniform random number per repeated sample for each chooser, but it also has the most simulation noise because small changes in approximate probabilities can change the sampled set substantially. - `poisson` is also relatively inexpensive. It draws one uniform random number per - chooser-alternative pair. With stable alternative alignment it is much less noisy - than Monte Carlo. + chooser-alternative pair (with stable alternative alignment, one per chooser and + stable-universe alternative, so inactive alternatives also consume draws). With stable + alternative alignment it is much less noisy than inverse-CDF sampling. - `eet` is the slowest sampling method. It draws one EV1 error term per chooser, alternative, and repeated sample draw. In return, it produces the most stable sampled sets across scenarios because unchanged alternatives keep the same unobserved error draws and only observed utility changes can change the sampled set. Note that `eet` does not remove the dependence on the approximate sampling utility itself: if that -utility changes, the sampled set can still change. What it removes is the extra Monte Carlo noise -from the sampling draw. `poisson` also benefits from stable alignment per alternative, but unlike -`eet` it still depends on probability-based inclusion tests. The practical effect on scenario -comparisons is ultimately empirical, but expected to be small. This was found to be the case for -test scenarios with an increase in employment in some zones, and therefore the sampling utility, -for the SANDAG example model. `poisson` is therefore the default sampling method when running in -explicit error term simulation mode. +utility changes, the sampled set can still change. What it removes is the extra noise from the +probability-space sampling draw. `poisson` also benefits from stable alignment per alternative, but +unlike `eet` it still depends on probability-based inclusion tests. The practical effect on +scenario comparisons is expected to be negligible, and empirical tests with an increase in +employment in some zones for the SANDAG example model confirm this. `poisson` is therefore the +default sampling method when running in explicit error term simulation mode. ## References diff --git a/docs/users-guide/sampling-methods.rst b/docs/users-guide/sampling-methods.rst index 22c47b0d5..c1b7324db 100644 --- a/docs/users-guide/sampling-methods.rst +++ b/docs/users-guide/sampling-methods.rst @@ -9,13 +9,13 @@ and location choice. Available methods are: -* ``monte_carlo``: importance sampling with replacement using probabilities and uniform draws +* ``inverse_cdf``: importance sampling with replacement using probabilities and uniform draws * ``eet``: importance sampling with replacement using explicit error-term draws * ``poisson``: independent Poisson inclusion sampling using probabilities Default behavior depends on the global simulation method setting: -* if ``use_explicit_error_terms: False``, the default sampling method is ``monte_carlo`` +* if ``use_explicit_error_terms: False``, the default sampling method is ``inverse_cdf`` * if ``use_explicit_error_terms: True``, the default sampling method is ``poisson`` However, any method can be used with either simulation method and can be set @@ -37,16 +37,16 @@ are simulated elsewhere in ActivitySim. Practical differences: -* ``monte_carlo`` and ``eet`` both sample with replacement, so duplicated sampled alternatives +* ``inverse_cdf`` and ``eet`` both sample with replacement, so duplicated sampled alternatives are possible and their aggregate sampled shares track repeated-draw MNL behavior more closely. * ``poisson`` samples alternatives by inclusion probability, so each sampled alternative appears at most once per chooser. This can change raw sampled shares in highly peaked cases, even though the downstream sampling correction remains well defined. -* ``monte_carlo`` is the fastest method, followed by ``poisson``, with ``eet`` being the slowest. +* ``inverse_cdf`` is the fastest method, followed by ``poisson``, with ``eet`` being the slowest. However, for models like location choice, most runtime comes from logsum calculations and the - total difference between ``monte_carlo`` and ``poisson`` sampling is usually very small. + total difference between ``inverse_cdf`` and ``poisson`` sampling is usually very small. * ``poisson`` is the current default when running with simulation method explicit error terms because it avoids repeated chooser-by-alternative explicit-error draws during sampling while - still providing improved noise reduction compared to Monte Carlo sampling. + still providing improved noise reduction compared to inverse-CDF sampling. For implementation details and runtime considerations, see :doc:`/dev-guide/sampling-methods`. diff --git a/docs/users-guide/ways_to_run.rst b/docs/users-guide/ways_to_run.rst index 18363fa3d..025d323b1 100644 --- a/docs/users-guide/ways_to_run.rst +++ b/docs/users-guide/ways_to_run.rst @@ -291,12 +291,12 @@ ____________________ ActivitySim makes heavy use of micro-simulation. Most model components are discrete choice models with an inherent random component, and for each choice situation a single outcome is generated. -With the default Monte Carlo draw method, ActivitySim first calculates analytical probabilities from the +With the default inverse-CDF (probability-based) draw method, ActivitySim first calculates analytical probabilities from the systematic utilities of a multinomial or nested logit model and then makes one draw from the cumulative distribution for each chooser. Explicit Error Terms (EET) replaces that final draw with a direct random-utility simulation by drawing the unobserved portion of utility (error term) for each chooser-alternative pair, adding it to the systematic utility, and selecting the alternative with the highest -total utility. Both methods simulate the same underlying model, but EET can be less affected by Monte Carlo +total utility. Both methods simulate the same underlying model, but EET can be less affected by simulation noise when comparing scenarios and can make some comparisons easier to interpret. This is because the selected alternative is the one with the highest total utility after adding the explicit error term, and if the explicit error term is consistent between a base and scenario run then @@ -311,6 +311,7 @@ To enable EET for a model run, set the global switch in ``settings.yaml``: Enable or disable this setting consistently across all runs being compared. For more details, including scenario comparison considerations, see :doc:`/dev-guide/explicit-error-terms`. + .. _skip_failed_choices_ways_to_run : Skip Failed Choices