From 258797aad390fa6f57549daa7746513a6ab9110c Mon Sep 17 00:00:00 2001 From: Jan Zill Date: Sun, 26 Jul 2026 17:19:59 +1000 Subject: [PATCH 1/6] removes duplicate call --- activitysim/core/interaction_sample_simulate.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/activitysim/core/interaction_sample_simulate.py b/activitysim/core/interaction_sample_simulate.py index c7d9d93d6..014397238 100644 --- a/activitysim/core/interaction_sample_simulate.py +++ b/activitysim/core/interaction_sample_simulate.py @@ -430,9 +430,6 @@ def _interaction_sample_simulate( # resulting pandas Int64Index has one element per chooser row and is in same order as choosers choices = alternatives[choice_column].take(positions + first_row_offsets) - # resulting pandas Int64Index has one element per chooser row and is in same order as choosers - choices = alternatives[choice_column].take(positions + first_row_offsets) - # create a series with index from choosers and the index of the chosen alternative choices = pd.Series(choices, index=choosers.index) From 0faf29480f0b24957b8fad18f1c8c5099bf3d2bb Mon Sep 17 00:00:00 2001 From: Jan Zill Date: Sun, 26 Jul 2026 17:59:29 +1000 Subject: [PATCH 2/6] addresses pr review comments --- .../abm/models/joint_tour_participation.py | 6 +-- activitysim/core/interaction_sample.py | 54 ++++++++----------- .../core/interaction_sample_simulate.py | 4 +- activitysim/core/logit.py | 7 +++ activitysim/core/random.py | 7 ++- 5 files changed, 39 insertions(+), 39 deletions(-) diff --git a/activitysim/abm/models/joint_tour_participation.py b/activitysim/abm/models/joint_tour_participation.py index 9b6a2a3d9..17013ebc8 100644 --- a/activitysim/abm/models/joint_tour_participation.py +++ b/activitysim/abm/models/joint_tour_participation.py @@ -219,10 +219,10 @@ def participants_chooser( # anybody with probability > 0 is forced to join the joint tour if state.settings.use_explicit_error_terms: # need "is valid choice" such that we certainly choose those with non-zero values, - # and do not choose others. Let's use 3.0 as large value here. + # and do not choose others. probs_or_utils[choice_col] = np.where( probs_or_utils[choice_col] > logit.UTIL_MIN, - 3.0, + logit.UTIL_LARGE_ENOUGH, logit.UTIL_UNAVAILABLE, ) non_choice_col = [ @@ -230,7 +230,7 @@ def participants_chooser( ][0] probs_or_utils[non_choice_col] = np.where( probs_or_utils[choice_col] <= logit.UTIL_MIN, - 3.0, + logit.UTIL_LARGE_ENOUGH, logit.UTIL_UNAVAILABLE, ) else: diff --git a/activitysim/core/interaction_sample.py b/activitysim/core/interaction_sample.py index 62eeababe..f45236b32 100644 --- a/activitysim/core/interaction_sample.py +++ b/activitysim/core/interaction_sample.py @@ -34,52 +34,42 @@ InteractionSampleMethod = typing.Literal["monte_carlo", "eet", "poisson"] -def _resolve_sample_method( - state: workflow.State, - compute_settings: ComputeSettings | None, -) -> InteractionSampleMethod: - sampling_method = None - if compute_settings is not None: - sampling_method = compute_settings.sample_method - if sampling_method is None: - sampling_method = state.settings.sample_method - if sampling_method is None: - sampling_method = ( - "poisson" if state.settings.use_explicit_error_terms else "monte_carlo" - ) - if sampling_method not in typing.get_args(InteractionSampleMethod): - raise ValueError( - f"Unsupported sample_method {sampling_method!r}; expected one of {typing.get_args(InteractionSampleMethod)}" - ) - logger.debug(f"Using sample_method={sampling_method}") - return sampling_method - - def resolve_sample_method( state: workflow.State, - model_settings, + settings: ComputeSettings | None = None, ) -> InteractionSampleMethod: """ - Resolve the sampling method for a model from its `model_settings`. + Resolve the sampling method to use, from most to least specific setting. Parameters ---------- state : workflow.State - model_settings : a pydantic model exposing an optional + settings : ComputeSettings or component settings, optional + Either a `ComputeSettings` directly, or a pydantic model exposing a `compute_settings` attribute (typically a `LogitComponentSettings` - subclass). If absent or None, the method is resolved purely from + subclass). If neither, the method is resolved purely from `state.settings`. Returns ------- sampling_method : InteractionSampleMethod """ - sample_compute_settings = getattr(model_settings, "compute_settings", None) - if sample_compute_settings is not None: - sample_compute_settings = sample_compute_settings.subcomponent_settings( - "sample" + # accept either a ComputeSettings or a component settings object wrapping one + compute_settings = getattr(settings, "compute_settings", settings) + + sampling_method = getattr(compute_settings, "sample_method", None) + if sampling_method is None: + sampling_method = state.settings.sample_method + if sampling_method is None: + sampling_method = ( + "poisson" if state.settings.use_explicit_error_terms else "monte_carlo" ) - return _resolve_sample_method(state, sample_compute_settings) + if sampling_method not in typing.get_args(InteractionSampleMethod): + raise ValueError( + f"Unsupported sample_method {sampling_method!r}; expected one of {typing.get_args(InteractionSampleMethod)}" + ) + logger.debug(f"Using sample_method={sampling_method}") + return sampling_method def _poisson_sample_alternatives_inner( @@ -780,7 +770,7 @@ def _interaction_sample( state.tracing.dump_df(DUMP, utilities, trace_label, "utilities") - sampling_method = _resolve_sample_method(state, compute_settings) + 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": @@ -1091,7 +1081,7 @@ def interaction_sample( if not choosers.index.is_monotonic_increasing: assert choosers.index.is_monotonic_increasing - sampling_method = _resolve_sample_method(state, compute_settings) + sampling_method = resolve_sample_method(state, compute_settings) logger.debug(f" interaction_sample sample method = {sampling_method}") if sampling_method == "monte_carlo": diff --git a/activitysim/core/interaction_sample_simulate.py b/activitysim/core/interaction_sample_simulate.py index 014397238..2083a6530 100644 --- a/activitysim/core/interaction_sample_simulate.py +++ b/activitysim/core/interaction_sample_simulate.py @@ -337,9 +337,7 @@ def _interaction_sample_simulate( if zero_probs.any(): # copied from proabability below, fix when that gets fixed # FIXME this is kind of gnarly, but we force choice of first alt - utilities_df.loc[ - zero_probs, 0 - ] = 3.0 # arbitrary value much larger than UTIL_UNAVAILABLE + utilities_df.loc[zero_probs, 0] = logit.UTIL_LARGE_ENOUGH # positions is series with the chosen alternative represented as a column index in utilities_df # which is an integer between zero and num alternatives in the alternative sample diff --git a/activitysim/core/logit.py b/activitysim/core/logit.py index 78f1ca114..60057ab2c 100644 --- a/activitysim/core/logit.py +++ b/activitysim/core/logit.py @@ -29,6 +29,13 @@ UTIL_MIN = np.log(EXP_UTIL_MIN, dtype=np.float64) UTIL_UNAVAILABLE = 1000.0 * (UTIL_MIN - 1.0) +# Utility assigned to an alternative that must be chosen over UTIL_UNAVAILABLE +# alternatives, e.g. when forcing a choice for a chooser that would otherwise +# have no available alternative. Any finite value works: the gap to +# UTIL_UNAVAILABLE is many orders of magnitude larger than the EV1 error terms +# added by the explicit error term choice methods. +UTIL_LARGE_ENOUGH = 3.0 + PROB_MIN = 0.0 PROB_MAX = 1.0 diff --git a/activitysim/core/random.py b/activitysim/core/random.py index 3524bc155..ccf1cdb55 100644 --- a/activitysim/core/random.py +++ b/activitysim/core/random.py @@ -461,7 +461,12 @@ def gumbel_choice_positions_for_df( positions = np.empty(n_rows, dtype=np.int64) if alt_nrs_df is not None: - assert alt_nrs_df.shape == utilities.shape + assert alt_nrs_df.index.equals( + utilities.index + ), "alt_nrs_df and utilities must share the same index" + assert alt_nrs_df.columns.equals( + utilities.columns + ), "alt_nrs_df and utilities must share the same columns" if n_rands is None: raise ValueError("n_rands is required when alt_nrs_df is provided") alt_nr_values = alt_nrs_df.to_numpy() From 5e5e8f253773e4c34d322d649885fad1f99dc1e6 Mon Sep 17 00:00:00 2001 From: Jan Zill Date: Sun, 26 Jul 2026 22:25:15 +1000 Subject: [PATCH 3/6] changes poisson sampling to be deterministic when no alternative sampled --- activitysim/core/interaction_sample.py | 287 ++++++++------- .../core/test/test_interaction_sample.py | 333 ++++++++++-------- docs/dev-guide/sampling-methods.md | 78 +++- 3 files changed, 391 insertions(+), 307 deletions(-) diff --git a/activitysim/core/interaction_sample.py b/activitysim/core/interaction_sample.py index f45236b32..dca38442b 100644 --- a/activitysim/core/interaction_sample.py +++ b/activitysim/core/interaction_sample.py @@ -33,6 +33,26 @@ InteractionSampleMethod = typing.Literal["monte_carlo", "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 +# that actually draw nothing always get the term regardless of this threshold, so it only +# governs how exact the reported probability is for choosers whose draw succeeded. +# +# Skipping the term avoids ranking the probability array for choosers that will essentially +# never need a fallback set, which costs about as much as the Bernoulli draw itself. The +# error it admits is bounded by the fact that P0 <= exp(-sample_size) whenever the +# probabilities sum to one (from 1 - p <= exp(-p)), so for a sample size of 30 the dropped +# term is at most 1e-13 and this branch is never taken at all for sample sizes above 27. +# +# Dropping the term understates prob by P0, so the relative error on the correction term +# log(1/prob) is P0/q_i, which is only large for an alternative whose own inclusion +# probability q_i is far below P0. Such an alternative can only be affected if it is +# sampled, which happens with probability q_i -- the same small quantity. That coupling +# bounds the expected number of chooser-alternative pairs whose correction is wrong by +# more than delta at n_choosers * sample_size * TOLERANCE / (exp(delta) - 1), i.e. below +# 1e-4 pairs off by more than 1 util in a one-million-chooser run. +POISSON_EMPTY_SAMPLE_TOLERANCE = 1e-12 + def resolve_sample_method( state: workflow.State, @@ -84,9 +104,8 @@ def _poisson_sample_alternatives_inner( """ Draw one Bernoulli inclusion decision per chooser-alternative pair. - Returns a dense 2-D array aligned to `probs` where sampled alternatives - contain their Poisson inclusion probability and unsampled alternatives are - `np.nan`. + Returns a dense 2-D boolean array aligned to `probs` that is True for the + sampled chooser-alternative pairs. """ if stable_alt_positions is None and n_total_alts is None: rands = rng.random_for_df(probs, n=probs.shape[1]) @@ -101,64 +120,32 @@ def _poisson_sample_alternatives_inner( "stable_alt_positions and n_total_alts must both be provided or omitted together" ) chunk_sizer.log_df(trace_label, "rands", rands) - return np.where( - rands < poisson_inclusion_probs_values, poisson_inclusion_probs_values, np.nan - ) + return rands < poisson_inclusion_probs_values -def _poisson_fallback_sample_alternatives( - probs: pd.DataFrame, +def _poisson_fallback_positions( + probs_values: np.ndarray, sample_size: int, - rng: Random, - trace_label: str | None, - chunk_sizer: ChunkSizer, - stable_alt_positions: np.ndarray | None = None, - n_total_alts: int | None = None, ) -> np.ndarray: """ - Fallback sampler used when Poisson retries still leave empty chooser rows. - - This path samples exactly `sample_size` distinct alternatives per chooser - without replacement by ranking one random score per alternative. The - returned array uses the same sparse chooser-by-alternative representation as - the Poisson path: chosen alternatives are `1.0`, unchosen alternatives are - `np.nan`. + Fallback choice set for choosers whose Poisson draw sampled no alternative. + + Returns the column positions of the `sample_size` highest-probability + alternatives for each row of `probs_values` (all of them if there are fewer + than `sample_size` alternatives), with ties broken by column position. + + 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, + which a data-dependent retry or redraw scheme cannot do. Because the fallback + set is a deterministic function of the probabilities, the probability that an + alternative ends up in the returned choice set still has an exact closed form + (see `_poisson_sample_alternatives`). """ - if sample_size > probs.shape[1]: - logger.info( - f"Poisson fallback sampling without replacement with sample_size={sample_size} > number of alternatives=" - + f"{probs.shape[1]}; returning all alternatives for {len(probs)} choosers" - ) - return np.full(probs.shape, 1.0) - - if stable_alt_positions is None and n_total_alts is None: - fallback_rands = rng.random_for_df(probs, n=probs.shape[1]) - elif stable_alt_positions is not None and n_total_alts is not None: - fallback_rands = rng.random_for_df_stable_alt_positions( - probs, - stable_alt_positions=stable_alt_positions, - n_total_alts=n_total_alts, - ) - else: - raise ValueError( - "stable_alt_positions and n_total_alts must both be provided or omitted together" - ) - chunk_sizer.log_df(trace_label, "fallback_rands", fallback_rands) - - chosen_positions = np.argpartition( - fallback_rands, - kth=sample_size - 1, - axis=1, - )[:, :sample_size] - - fallback_sampled_values = np.full(probs.shape, np.nan) - chooser_positions = np.repeat(np.arange(len(probs)), sample_size) - fallback_sampled_values[ - chooser_positions, - chosen_positions.reshape(-1), - ] = 1.0 - - return fallback_sampled_values + k = min(sample_size, probs_values.shape[1]) + # stable sort of the negated probabilities gives descending probability order + # with ties broken by column position + return np.argsort(-probs_values, axis=1, kind="stable")[:, :k] def make_sample_choices_eet( @@ -231,105 +218,109 @@ def _poisson_sample_alternatives( """ Build a Poisson-sampled choice set for each chooser. - The primary path performs independent Poisson inclusion draws for every chooser-alternative pair and retries any - chooser row that sampled no alternatives. Internally the sampler maintains a sparse chooser-by-alternative array - where sampled cells hold the probability to carry forward as `prob` and unsampled cells are np.nan. - - If a chooser still has no sampled alternatives after 10 retries, we fall back to sampling exactly sample_size - distinct alternatives without replacement and force those chosen probabilities to `1.0` so the sampling correction - factor cancels out. In practice we expect this to be very rare with reasonable sample sizes and not too small - choice sets, but it is a known issue with Poisson sampling that we want to guard against. Note that if this - fallback is triggered it can lead to inconsistent random numbers between two scenarios if the number of retries it - takes in each scenario differs, but again we expect this to be very rare and the alternative is potentially - infinite retries or raising an error. - - returns: DataFrame with one row per sampled chooser-alternative pair and columns for chooser index, alt_col_name, - and prob (the Poisson inclusion probability for that pair). - - In the case of Poisson sampling, the inclusion probability for each chooser-alternative pair is the probability - that the alternative was included in the sample at least once across the sample_size draws, which is the - reciprocal of it never being drawn in sample_size draws, so 1-(1-p)^sample_size where p is the - original choice probability. To make Poisson sampling interchangeable with other sampling methods, we return the - inclusion probabilities i.e. the true probability of being sampled. Pick_count will be 1 by definition - (poisson sampling returns a yes/no for each alternative, so if an alternative is included in the sample it is - included once) and the standard sampling correction factor can be recovered as np.log(df.pick_count/df.prob) - = np.log(1/inclusion_prob). + Every chooser-alternative pair gets one independent Bernoulli inclusion draw with + probability + + q_i = 1 - (1 - p_i) ** sample_size + + 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 + 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)`. + + Because the draws are independent, a chooser can end up with no sampled alternatives + at all. That happens with probability + + P0 = prod_j (1 - q_j) + + 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 + 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. + + Determinism also makes the reported probability exact. Alternative `i` ends up in the + returned choice set if it was drawn, or if nothing was drawn and it is in the fallback + set. An alternative that was drawn cannot also have been in an empty draw, so these are + disjoint, and the fallback set is fixed rather than random, giving + + prob_i = q_i + P0 * 1{i in fallback set} + + unconditionally: the same value for every chooser whatever branch it actually took, and + bounded by 1 because P0 <= 1 - q_i. Reporting `q_i` alone would understate the + correction for exactly the choosers most at risk of an empty draw. Ranking the + probabilities to find the fallback set costs about as much as the Bernoulli draw itself, + so the term is only evaluated for choosers whose P0 exceeds + `POISSON_EMPTY_SAMPLE_TOLERANCE`, plus every chooser that actually drew nothing. See + that constant for the error this admits. + + Note this is a different sampling design from retrying until the draw is non-empty, + which would give `q_i / (1 - P0)` instead. Both are valid; this one has an exact closed + form that does not depend on how many times a chooser was redrawn. + + returns: DataFrame with one row per sampled chooser-alternative pair and columns for + chooser index, alt_col_name, and prob. """ - inclusion_probs_values = 1.0 - np.power( - 1.0 - probs.to_numpy(copy=False), sample_size - ) + probs_values = probs.to_numpy(copy=False) - sampled_values = np.full(inclusion_probs_values.shape, np.nan) - - n = 0 - active_row_positions = np.arange(len(probs), dtype=np.int64) - - while active_row_positions.size > 0: - probs_subset = probs.iloc[active_row_positions] - # Each retry call advances the per-row offset by n_total_alts on the - # underlying RNG channel (see SimpleChannel.random_for_df_stable_alt_positions). - # The number of retries is data-dependent, so two scenarios with slightly - # different inclusion probabilities for the same chooser can end up with - # different downstream RNG offsets — defeating cross-scenario stability - # for that chooser. This is rare (only fires when the first Poisson draw - # yields zero samples) but worth keeping in mind for base/project - # comparisons; see the docstring above for context. - sampled_results_subset = _poisson_sample_alternatives_inner( - probs_subset, - inclusion_probs_values[active_row_positions], - state.get_rn_generator(), - trace_label, - chunk_sizer, - stable_alt_positions=stable_alt_positions, - n_total_alts=n_total_alts, + # q_i: probability of alternative i being included at least once in sample_size draws + inclusion_probs = 1.0 - np.power(1.0 - probs_values, sample_size) + + # P0: probability that a chooser's Bernoulli draws include nothing at all. Must be + # computed before inclusion_probs is updated in place below. + empty_sample_probs = np.prod(1.0 - inclusion_probs, axis=1) + + sampled = _poisson_sample_alternatives_inner( + probs, + inclusion_probs, + state.get_rn_generator(), + trace_label, + chunk_sizer, + stable_alt_positions=stable_alt_positions, + n_total_alts=n_total_alts, + ) + chunk_sizer.log_df(trace_label, "sampled", sampled) + + empty_rows = ~sampled.any(axis=1) + + n_empty = int(empty_rows.sum()) + 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"{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}." ) - no_alts_sampled_mask = np.isnan(sampled_results_subset).all(axis=1) - sampled_values[ - active_row_positions[~no_alts_sampled_mask] - ] = sampled_results_subset[~no_alts_sampled_mask] - - if no_alts_sampled_mask.any(): - failed_row_positions = active_row_positions[no_alts_sampled_mask] - extra_per_retry = ( - f"{n_total_alts}" if n_total_alts is not None else "n_total_alts" - ) - logger.warning( - f"Poisson sampling of alternatives failed for {len(failed_row_positions)} " - f"chooser(s) with {n=}, retrying. Note: retried choosers consume an extra " - f"{extra_per_retry} randoms per retry on this RNG channel, which can cause " - f"downstream RNG offset divergence vs scenarios that did not retry." - ) - logger.debug( - f"Sampled size was {sample_size}, poisson method mean expected sample size was" - + f" {inclusion_probs_values[failed_row_positions].sum(axis=1).mean():.1f}, actual sampled mean was" - + f" {np.isfinite(sampled_values[failed_row_positions]).sum(axis=1).mean():.1f} and highest zero" - + f" selection prob was {(1.0 - inclusion_probs_values[failed_row_positions]).prod(axis=1).max():.2g}" - ) - active_row_positions = failed_row_positions - else: # All choosers have at least one alternative in sample set - break + # inclusion_probs is updated in place from q_i to the reported probability + # q_i + P0 * 1{i in fallback set} + fallback_rows = np.nonzero( + empty_rows | (empty_sample_probs > POISSON_EMPTY_SAMPLE_TOLERANCE) + )[0] + if fallback_rows.size > 0: + fallback_cols = _poisson_fallback_positions( + probs_values[fallback_rows], sample_size + ) + row_positions = np.repeat(fallback_rows, fallback_cols.shape[1]) + col_positions = fallback_cols.reshape(-1) + inclusion_probs[row_positions, col_positions] += empty_sample_probs[ + row_positions + ] - n += 1 - if n == 10: - logger.info( - "Poisson choice set sampling exceeded 10 retries; falling back to random sampling for %s choosers", - len(active_row_positions), - ) - fallback_sampled_values = _poisson_fallback_sample_alternatives( - probs.iloc[active_row_positions], - sample_size, - state.get_rn_generator(), - trace_label, - chunk_sizer, - stable_alt_positions=stable_alt_positions, - n_total_alts=n_total_alts, - ) - sampled_values[active_row_positions] = fallback_sampled_values - break + # ...but only the choosers that actually drew nothing take the fallback set + takes_fallback = empty_rows[row_positions] + sampled[row_positions[takes_fallback], col_positions[takes_fallback]] = True - chooser_positions, alt_positions = np.nonzero(~np.isnan(sampled_values)) + chooser_positions, alt_positions = np.nonzero(sampled) chooser_col_name = probs.index.name or "index" if len(chooser_positions) == 0: @@ -338,7 +329,7 @@ def _poisson_sample_alternatives( choices_df = pd.DataFrame( { chooser_col_name: probs.index.to_numpy()[chooser_positions], - "prob": sampled_values[chooser_positions, alt_positions], + "prob": inclusion_probs[chooser_positions, alt_positions], alt_col_name: alternatives.index.to_numpy()[alt_positions], } ) diff --git a/activitysim/core/test/test_interaction_sample.py b/activitysim/core/test/test_interaction_sample.py index a25ebc5f3..b520c264c 100644 --- a/activitysim/core/test/test_interaction_sample.py +++ b/activitysim/core/test/test_interaction_sample.py @@ -650,17 +650,59 @@ def gumbel_max_positions_for_df( ) -def _expected_choices_df(sampled_alternatives, alternatives, alt_col_name): - return ( - sampled_alternatives.rename_axis("alt_idx", axis=1) - .stack() - .reset_index(name="prob") - .assign(**{alt_col_name: lambda df: alternatives.index.values[df["alt_idx"]]}) - .drop(columns=["alt_idx"]) +def _reference_poisson_sampled_values(probs_np, draws, sample_size): + """ + Independent re-derivation of the documented Poisson sampling result, used to check + the implementation against the formula rather than against itself. + + 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. + + Returns the sparse chooser-by-alternative array of reported probabilities, with + np.nan for alternatives that are not in the choice set. + """ + inclusion_probs = 1.0 - np.power(1.0 - probs_np, sample_size) + empty_sample_probs = np.prod(1.0 - inclusion_probs, axis=1) + + sampled = draws < inclusion_probs + empty_rows = ~sampled.any(axis=1) + + in_fallback = np.zeros(probs_np.shape, dtype=bool) + 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) + + # the implementation skips the P0 term where it cannot matter; mirror that here so + # the comparison stays exact (see POISSON_EMPTY_SAMPLE_TOLERANCE) + material = empty_sample_probs > interaction_sample.POISSON_EMPTY_SAMPLE_TOLERANCE + reported = inclusion_probs + empty_sample_probs[:, None] * ( + in_fallback & (material | empty_rows)[:, None] + ) + + sampled = sampled | (empty_rows[:, None] & in_fallback) + return np.where(sampled, reported, np.nan) + + +def _reference_poisson_choices_df( + probs, draws, sample_size, alternatives, alt_col_name +): + """Flatten `_reference_poisson_sampled_values` into the expected choices frame.""" + sampled_values = _reference_poisson_sampled_values( + probs.to_numpy(), draws, sample_size + ) + chooser_idx, alt_idx = np.nonzero(~np.isnan(sampled_values)) + return pd.DataFrame( + { + probs.index.name: probs.index.to_numpy()[chooser_idx], + "prob": sampled_values[chooser_idx, alt_idx], + alt_col_name: alternatives.index.to_numpy()[alt_idx], + } ) -def test_poisson_sample_alternatives_inner_returns_masked_inclusion_probs(): +def test_poisson_sample_alternatives_inner_returns_inclusion_mask(): probs = pd.DataFrame( [[0.2, 0.4, 0.6], [0.1, 0.3, 0.5]], index=pd.Index([11, 17], name="person_id"), @@ -683,50 +725,47 @@ def test_poisson_sample_alternatives_inner_returns_masked_inclusion_probs(): probs, inclusion_probs_values, rng, - trace_label="test_poisson_sample_alternatives_inner_returns_masked_inclusion_probs", + trace_label="test_poisson_sample_alternatives_inner_returns_inclusion_mask", chunk_sizer=_DummyChunkSizer(), ) expected = np.array( - [[0.36, np.nan, 0.84], [np.nan, 0.51, np.nan]], - dtype=np.float64, + [[True, False, True], [False, True, False]], + dtype=bool, ) - np.testing.assert_allclose(sampled, expected, equal_nan=True) + np.testing.assert_array_equal(sampled, expected) -def test_poisson_fallback_sample_alternatives_selects_distinct_positions_with_prob_one(): - probs = pd.DataFrame( +def test_poisson_fallback_positions_selects_highest_probability_alternatives(): + probs_values = np.array( [[0.20, 0.30, 0.50, 0.00], [0.40, 0.10, 0.30, 0.20]], - index=pd.Index([11, 17], name="person_id"), - columns=np.arange(4), - ) - rng = _SequentialDummyRng( - [ - np.array( - [[0.90, 0.10, 0.40, 0.20], [0.05, 0.70, 0.60, 0.10]], - dtype=np.float64, - ) - ] + dtype=np.float64, ) - sampled = interaction_sample._poisson_fallback_sample_alternatives( - probs=probs, - sample_size=2, - rng=rng, - trace_label="test_poisson_fallback_sample_alternatives_selects_distinct_positions_with_prob_one", - chunk_sizer=_DummyChunkSizer(), - ) + positions = interaction_sample._poisson_fallback_positions(probs_values, 2) - expected = np.array( - [[np.nan, 1.0, np.nan, 1.0], [1.0, np.nan, np.nan, 1.0]], - dtype=np.float64, + # highest probability first, so [0.50, 0.30] and [0.40, 0.30] + np.testing.assert_array_equal(positions, np.array([[2, 1], [0, 2]])) + + +def test_poisson_fallback_positions_breaks_ties_by_column_and_caps_at_alt_count(): + probs_values = np.array([[0.25, 0.25, 0.25, 0.25]], dtype=np.float64) + + # ties resolve to the leading columns, deterministically + np.testing.assert_array_equal( + interaction_sample._poisson_fallback_positions(probs_values, 2), + np.array([[0, 1]]), ) - np.testing.assert_allclose(sampled, expected, equal_nan=True) + # asking for more alternatives than exist returns all of them + np.testing.assert_array_equal( + interaction_sample._poisson_fallback_positions(probs_values, 99), + np.array([[0, 1, 2, 3]]), + ) -def test_poisson_sample_alternatives_retries_and_returns_expected_frames(): +def test_poisson_sample_alternatives_returns_expected_frames(): probs = pd.DataFrame( [ [0.20, 0.60, 0.10, 0.05], @@ -738,36 +777,16 @@ def test_poisson_sample_alternatives_retries_and_returns_expected_frames(): ) sample_size = 2 alternatives = pd.DataFrame(index=pd.Index([100, 300, 700, 900], name="alt_id")) - expected_inclusion_probs = 1 - (1 - probs) ** sample_size - expected_sampled_alternatives = pd.DataFrame( + # the middle chooser samples nothing and takes the fallback set + draws = np.array( [ - [expected_inclusion_probs.iloc[0, 0], np.nan, np.nan, np.nan], - [ - expected_inclusion_probs.iloc[1, 0], - expected_inclusion_probs.iloc[1, 1], - np.nan, - np.nan, - ], - [np.nan, np.nan, expected_inclusion_probs.iloc[2, 2], np.nan], + [0.10, 0.90, 0.50, 0.90], + [0.90, 0.90, 0.90, 0.90], + [0.80, 0.90, 0.20, 0.80], ], - index=probs.index, - columns=probs.columns, - ) - state = _DummyState( - _SequentialDummyRng( - [ - np.array( - [ - [0.10, 0.90, 0.50, 0.90], - [0.90, 0.90, 0.90, 0.90], - [0.80, 0.90, 0.20, 0.80], - ], - dtype=np.float64, - ), - np.array([[0.10, 0.05, 0.70, 0.80]], dtype=np.float64), - ] - ) + dtype=np.float64, ) + state = _DummyState(_SequentialDummyRng([draws])) choices_df = interaction_sample._poisson_sample_alternatives( chunk_sizer=_DummyChunkSizer(), @@ -776,19 +795,21 @@ def test_poisson_sample_alternatives_retries_and_returns_expected_frames(): sample_size=sample_size, alt_col_name="alt_id", state=state, - trace_label="test_poisson_sample_alternatives_retries_and_returns_expected_frames", + trace_label="test_poisson_sample_alternatives_returns_expected_frames", ) - expected_choices_df = _expected_choices_df( - expected_sampled_alternatives, - alternatives, - "alt_id", + expected = _reference_poisson_choices_df( + probs, draws, sample_size, alternatives, "alt_id" ) + pd.testing.assert_frame_equal(choices_df, expected) - pd.testing.assert_frame_equal(choices_df, expected_choices_df) + # the fallback chooser gets the two most likely alternatives, 100 and 700 + assert choices_df.loc[choices_df.person_id == 17, "alt_id"].tolist() == [100, 700] -def test_poisson_sample_alternatives_falls_back_to_random_sampling_after_ten_retries(): +def test_poisson_sample_alternatives_consumes_no_extra_randoms_on_empty_draw(): + # the fallback must not draw again: _SequentialDummyRng raises IndexError if the + # sampler asks for a second block, so a single draw array is the assertion here probs = pd.DataFrame( [[0.20, 0.30, 0.50]], index=pd.Index([11], name="person_id"), @@ -797,8 +818,7 @@ def test_poisson_sample_alternatives_falls_back_to_random_sampling_after_ten_ret sample_size = 2 alternatives = pd.DataFrame(index=pd.Index([100, 300, 700], name="alt_id")) fail_draw = np.array([[0.99, 0.99, 0.99]], dtype=np.float64) - fallback_draw = np.array([[0.10, 0.80, 0.20]], dtype=np.float64) - state = _DummyState(_SequentialDummyRng([fail_draw] * 10 + [fallback_draw])) + state = _DummyState(_SequentialDummyRng([fail_draw])) choices_df = interaction_sample._poisson_sample_alternatives( chunk_sizer=_DummyChunkSizer(), @@ -807,21 +827,74 @@ def test_poisson_sample_alternatives_falls_back_to_random_sampling_after_ten_ret sample_size=sample_size, alt_col_name="alt_id", state=state, - trace_label="test_poisson_sample_alternatives_falls_back_to_random_sampling_after_ten_retries", + trace_label="test_poisson_sample_alternatives_consumes_no_extra_randoms_on_empty_draw", ) - expected_sampled_alternatives = pd.DataFrame( - [[1.0, np.nan, 1.0]], - index=probs.index, - columns=probs.columns, + # the two most likely alternatives, reported at q_i + P0 + 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, 1] + empty_sample_prob, + inclusion_probs[0, 2] + empty_sample_prob, + ], + "alt_id": [300, 700], + } ) - expected_choices_df = _expected_choices_df( - expected_sampled_alternatives, - alternatives, - "alt_id", + + pd.testing.assert_frame_equal(choices_df, expected) + + +def test_poisson_sample_alternatives_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 + # fallback. Every chooser is identical, so the empirical inclusion rate across + # choosers estimates that probability directly. sample_size=1 over 6 uniform + # alternatives makes empty draws frequent (P0 = (5/6)^6 ~ 0.33), which is what puts + # the fallback term under test. + n_choosers = 200_000 + n_alts = 6 + sample_size = 1 + + probs = pd.DataFrame( + np.full((n_choosers, n_alts), 1.0 / n_alts), + index=pd.Index(np.arange(n_choosers), name="person_id"), + columns=np.arange(n_alts), + ) + alternatives = pd.DataFrame(index=pd.Index(np.arange(n_alts) * 10, name="alt_id")) + + inclusion_probs = 1 - np.power(1 - probs.to_numpy(), sample_size) + empty_sample_prob = np.prod(1 - inclusion_probs, axis=1)[0] + assert empty_sample_prob > 0.3 + + draws = np.random.default_rng(20260726).random((n_choosers, n_alts)) + state = _DummyState(_SequentialDummyRng([draws])) + + choices_df = interaction_sample._poisson_sample_alternatives( + chunk_sizer=_DummyChunkSizer(), + probs=probs, + alternatives=alternatives, + sample_size=sample_size, + alt_col_name="alt_id", + state=state, + trace_label="test_poisson_sample_alternatives_reported_prob_is_total_inclusion_probability", ) - pd.testing.assert_frame_equal(choices_df, expected_choices_df) + # identical choosers must get an identical reported prob per alternative, whether + # they reached the choice set through the Bernoulli draw or through the fallback + reported = choices_df.groupby("alt_id")["prob"].agg(["min", "max", "first"]) + np.testing.assert_allclose(reported["min"], reported["max"], rtol=1e-12) + + # ties in the fallback resolve to the first column, so only alternative 0 carries + # the extra P0 mass + expected_reported = np.full(n_alts, inclusion_probs[0, 0]) + expected_reported[0] += empty_sample_prob + np.testing.assert_allclose(reported["first"], expected_reported, rtol=1e-12) + + empirical = choices_df.groupby("alt_id").size() / n_choosers + np.testing.assert_allclose(empirical, expected_reported, atol=0.005) def test_poisson_sample_alternatives_repeat_alignment_chooser_dominant_heterogeneity(): @@ -875,23 +948,14 @@ def test_poisson_sample_alternatives_repeat_alignment_chooser_dominant_heterogen trace_label="test_repeat_alignment_chooser_heterogeneity", ) - probs_np = probs.to_numpy() - inclusion_probs = 1 - np.power(1 - probs_np, sample_size) - sampled_values = np.where(poisson_draws < inclusion_probs, inclusion_probs, np.nan) - chooser_idx, alt_idx = np.nonzero(~np.isnan(sampled_values)) - - expected = pd.DataFrame( - { - "person_id": chooser_index.to_numpy()[chooser_idx], - "prob": sampled_values[chooser_idx, alt_idx], - "alt_id": alternatives.index.to_numpy()[alt_idx], - } + expected = _reference_poisson_choices_df( + probs, poisson_draws, sample_size, alternatives, "alt_id" ) pd.testing.assert_frame_equal(out.reset_index(drop=True), expected) -def test_poisson_sample_alternatives_retry_matches_materialized_path(): +def test_poisson_sample_alternatives_matches_materialized_path(): chooser_index = pd.Index([201, 202, 203], name="person_id") choosers = pd.DataFrame(index=chooser_index) alternatives = pd.DataFrame(index=pd.Index([10, 11, 12, 13], name="alt_id")) @@ -908,8 +972,7 @@ def test_poisson_sample_alternatives_retry_matches_materialized_path(): ], dtype=np.float64, ) - retry_draw = np.array([[0.40, 0.10, 0.90, 0.90]], dtype=np.float64) - state = _DummyState(_SequentialDummyRng([poisson_draws, retry_draw])) + state = _DummyState(_SequentialDummyRng([poisson_draws])) probs = interaction_sample.logit.utils_to_probs( state, @@ -930,26 +993,8 @@ def test_poisson_sample_alternatives_retry_matches_materialized_path(): trace_label="test_fused_rng_matches_materialized", ) - probs_np = probs.to_numpy() - inclusion_probs = 1 - np.power(1 - probs_np, sample_size) - sampled_values = np.full(inclusion_probs.shape, np.nan) - first_pass = np.where(poisson_draws < inclusion_probs, inclusion_probs, np.nan) - first_pass_empty = np.isnan(first_pass).all(axis=1) - sampled_values[~first_pass_empty] = first_pass[~first_pass_empty] - retry_pass = np.where( - retry_draw < inclusion_probs[first_pass_empty], - inclusion_probs[first_pass_empty], - np.nan, - ) - sampled_values[first_pass_empty] = retry_pass - chooser_idx, alt_idx = np.nonzero(~np.isnan(sampled_values)) - - expected = pd.DataFrame( - { - "person_id": choosers.index.values[chooser_idx], - "prob": sampled_values[chooser_idx, alt_idx], - "alt_id": alternatives.index.values[alt_idx], - } + expected = _reference_poisson_choices_df( + probs, poisson_draws, sample_size, alternatives, "alt_id" ) pd.testing.assert_frame_equal(out.reset_index(drop=True), expected) @@ -1116,26 +1161,18 @@ def test_poisson_sample_alternatives_stable_alt_mapping_matches_materialized_pat n_total_alts=n_total_alts, ) - probs_np = probs.to_numpy() - inclusion_probs = 1 - np.power(1 - probs_np, sample_size) - active_uniforms = dense_uniforms[:, stable_alt_positions] - sampled_values = np.where( - active_uniforms < inclusion_probs, inclusion_probs, np.nan - ) - chooser_idx, alt_idx = np.nonzero(~np.isnan(sampled_values)) - - expected = pd.DataFrame( - { - "person_id": choosers.index.values[chooser_idx], - "prob": sampled_values[chooser_idx, alt_idx], - "alt_id": alternatives.index.values[alt_idx], - } + expected = _reference_poisson_choices_df( + probs, + dense_uniforms[:, stable_alt_positions], + sample_size, + alternatives, + "alt_id", ) pd.testing.assert_frame_equal(out.reset_index(drop=True), expected) -def test_poisson_sample_alternatives_falls_back_after_retries(): +def test_poisson_sample_alternatives_falls_back_to_most_likely_alternatives(): chooser_index = pd.Index([301, 302], name="person_id") choosers = pd.DataFrame(index=chooser_index) alternatives = pd.DataFrame(index=pd.Index([10, 12, 14], name="alt_id")) @@ -1145,20 +1182,13 @@ def test_poisson_sample_alternatives_falls_back_after_retries(): ) sample_size = 2 fail_draw = np.full((2, 3), 0.99, dtype=np.float64) - fallback_draw = np.array( - [ - [0.40, 0.10, 0.20], - [0.30, 0.20, 0.90], - ], - dtype=np.float64, - ) - state = _DummyState(_SequentialDummyRng([fail_draw] * 10 + [fallback_draw])) + state = _DummyState(_SequentialDummyRng([fail_draw])) probs = interaction_sample.logit.utils_to_probs( state, utilities, allow_zero_probs=False, - trace_label="test_falls_back_after_retries", + trace_label="test_falls_back_to_most_likely_alternatives", overflow_protection=True, trace_choosers=choosers, ) @@ -1170,15 +1200,22 @@ def test_poisson_sample_alternatives_falls_back_after_retries(): sample_size=sample_size, alt_col_name="alt_id", state=state, - trace_label="test_falls_back_after_retries", + trace_label="test_falls_back_to_most_likely_alternatives", ) - expected = pd.DataFrame( - { - "person_id": [301, 301, 302, 302], - "prob": [1.0, 1.0, 1.0, 1.0], - "alt_id": [12, 14, 10, 12], - } - ) + # neither chooser sampled anything, so both take their two most likely + # alternatives: utilities [0.0, 0.3, -0.2] -> alts 12, 10 and + # [1.0, 0.2, 0.4] -> alts 10, 14 + assert out["alt_id"].tolist() == [10, 12, 10, 14] + assert out["person_id"].tolist() == [301, 301, 302, 302] + expected = _reference_poisson_choices_df( + probs, fail_draw, sample_size, alternatives, "alt_id" + ) pd.testing.assert_frame_equal(out.reset_index(drop=True), expected) + + # reported prob is q_i + P0, strictly below 1 and strictly above the bare q_i + inclusion_probs = 1 - np.power(1 - probs.to_numpy(), sample_size) + empty_sample_probs = np.prod(1 - inclusion_probs, axis=1) + assert (out["prob"] < 1.0).all() + assert (empty_sample_probs > 0).all() diff --git a/docs/dev-guide/sampling-methods.md b/docs/dev-guide/sampling-methods.md index 53c098450..700c15c30 100644 --- a/docs/dev-guide/sampling-methods.md +++ b/docs/dev-guide/sampling-methods.md @@ -97,12 +97,36 @@ noticeably from repeated-draw MNL shares in highly peaked cases. This is structu numerical noise. The interaction-sample tests document this explicitly. --> A chooser can occasionally receive no sampled alternatives under Poisson sampling, because each -alternative is tested independently. In the models that use sampling in ActivitySim, this should be -rare. If it happens, the sampler retries that chooser row up to 10 times and then falls back to a -simple without-replacement random sample. - +alternative is tested independently. The probability of this happening for a given chooser is + +$$ +P_0 = \prod_j (1 - p_j)^s +$$ + +Because the probabilities sum to one and $1 - p \le e^{-p}$, this is bounded above by $e^{-s}$ +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. + +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 +redraw scheme cannot do this: the number of retries is data-dependent, so two nearby scenarios would +consume different numbers of randoms for the same chooser and desynchronise every draw after it, +which is undesired when running in explicit error term simulation mode. Determinism also keeps the +reported `prob` exact, see below. + +Taking the $s$ most likely alternatives is only sound as a rare repair, not as a sampling method in +its own right. Used on its own it would give every selected alternative an inclusion probability of +1 and every other alternative an inclusion probability of 0, so the correction term would be the +same constant for all selected alternatives and would cancel out of the choice entirely. The result +is a plain MNL over the top $s$ alternatives, and the choice mass on all remaining alternatives is +lost. Because the sampling utility is a deliberately cheap approximation, the alternatives it ranks +poorly are not the same ones the final utility ranks poorly, so this is a systematic bias rather +than sampling noise. The problem is the deterministic *exclusion*, not the determinism itself: an +inclusion probability of 1 is perfectly valid, but one of 0 cannot be corrected for by any +weighting. Here the Bernoulli draw keeps every alternative's inclusion probability strictly +positive, and the fallback can only add inclusion mass on top of that, never remove it. ### Sampling Correction @@ -134,12 +158,41 @@ 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 returned `prob` is: +and the configured sample size is $s$, then the inclusion probably of the Bernoulli trial is $$ -1 - (1 - p)^s +q_i = 1 - (1 - p_i)^s $$ +An alternative ends up in the returned choice set either because its inclusion draw succeeded, +or because the chooser drew nothing at all and the alternative is in the fallback set. Something +that was drawn cannot also have been part of an empty draw, so these two events are disjoint, and +because the fallback set is deterministic rather than random the returned `prob` is + +$$ +\text{prob}_i = q_i + P_0 \cdot 1\{i \in \text{fallback set}\} +$$ + +Note this is the *unconditional* probability, not either branch on its own. Conditional on the draw +being non-empty the inclusion probability is $q_i / (1 - P_0)$, and conditional on it being empty it +is $1\{i \in \text{fallback set}\}$; mixing those with weights $1 - P_0$ and $P_0$ recovers the +expression above. The unconditional form is the one the correction needs, and it is also what makes +the reported `prob` independent of which branch a given chooser happened to take. + +The conditional form $q_i / (1 - P_0)$ is what a design that retried until the draw was non-empty +would have to report. Both designs are valid samplers; this one has an exact closed form that does +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 +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 +alternative can only be affected if it is sampled, which happens with probability $q_i$, the same +small quantity. That coupling keeps the expected number of materially wrong corrections far below +one for any model size. + Since `pick_count` is always `1` for `poisson`, the correction becomes $\log(1 / \text{prob})$. This means that all three methods use the same correction expression, @@ -156,8 +209,8 @@ Runtime and noise characteristics differ across methods. 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 possible retries for chooser rows that initially sample no - alternatives. With stable alternative alignment it is much less noisy than Monte Carlo. + chooser-alternative pair. With stable alternative alignment it is much less noisy + than Monte Carlo. - `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 @@ -167,7 +220,10 @@ Note that `eet` does not remove the dependence on the approximate sampling utili 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. +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. ## References From 5cd92bcb8bdd55dda6d89d8303d76b689943e084 Mon Sep 17 00:00:00 2001 From: Jan Zill Date: Sun, 26 Jul 2026 22:35:28 +1000 Subject: [PATCH 4/6] do not allow zero probs for nl leafs --- activitysim/core/simulate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/activitysim/core/simulate.py b/activitysim/core/simulate.py index 4d143e119..dd5a22869 100644 --- a/activitysim/core/simulate.py +++ b/activitysim/core/simulate.py @@ -1553,7 +1553,7 @@ def eval_nl( if state.settings.use_explicit_error_terms: raw_utilities = logit.validate_utils( - state, raw_utilities, allow_zero_probs=True, trace_label=trace_label + state, raw_utilities, allow_zero_probs=False, trace_label=trace_label ) if custom_chooser: From ad476ea884be5d4368348146d1183a2908314d70 Mon Sep 17 00:00:00 2001 From: Jan Zill Date: Mon, 27 Jul 2026 00:08:09 +1000 Subject: [PATCH 5/6] handle non-available alternatives in make_choices_utility_based like MC in make_choices --- activitysim/core/logit.py | 38 ++++++++++++-- activitysim/core/random.py | 32 +++++++----- activitysim/core/simulate.py | 2 +- activitysim/core/test/test_logit.py | 76 ++++++++++++++++++++++++++++ activitysim/core/test/test_random.py | 72 ++++++++++++++++++++++++++ 5 files changed, 202 insertions(+), 18 deletions(-) diff --git a/activitysim/core/logit.py b/activitysim/core/logit.py index 60057ab2c..9c3f668d2 100644 --- a/activitysim/core/logit.py +++ b/activitysim/core/logit.py @@ -647,6 +647,7 @@ def make_choices_utility_based( nest_spec=None, alts_context: AltsContext | None = None, alt_nrs_df: pd.DataFrame | None = None, + allow_zero_probs: bool = False, ) -> tuple[pd.Series, pd.Series]: """ Make choices for each chooser from among a set of alternatives based on utilities by adding @@ -673,6 +674,11 @@ def make_choices_utility_based( the alt_nrs for each alternative for each chooser. This is used to index into the random numbers when sampling EET terms for multinomial logit models, and should contain -999 for any alternatives that are not available for a given chooser. Should be provided along with `alts_context`. + allow_zero_probs : bool + If False, report choosers that have no available alternative at all. This is the + utility-space counterpart of `allow_bad_probs` in `make_choices`, and takes the + name `allow_zero_probs` because the condition it suppresses is the same one that + `validate_utils` and `utils_to_probs` describe by that name. Returns ------- @@ -684,14 +690,36 @@ def make_choices_utility_based( Notes ----- - Bad-row reporting (e.g., a chooser whose alternatives are all `UTIL_UNAVAILABLE`) is the - responsibility of `validate_utils()`, which is invoked at every EET call site - (interaction_sample, interaction_sample_simulate, simulate.eval_mnl, simulate.eval_nl) - BEFORE this function is called. EET argmax always returns a valid integer position; - we do not re-check here. + 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") + if not allow_zero_probs: + # An alternative counts as available when its utility exceeds UTIL_MIN, which is the + # same threshold validate_utils uses before clamping to UTIL_UNAVAILABLE. Testing the + # threshold directly rather than testing the row sum against a multiple of + # UTIL_UNAVAILABLE means this also holds for callers that arrive without having run + # validate_utils, where unavailable alternatives may still carry their raw value. + no_available_alts = ~(utilities.to_numpy() > UTIL_MIN).any(axis=1) + if no_available_alts.any(): + report_bad_choices( + state, + no_available_alts, + utilities, + state.settings.skip_failed_choices, + trace_label=tracing.extend_trace_label( + trace_label, "no_available_alts" + ), + msg="no alternative is available", + trace_choosers=trace_choosers, + ) + if nest_spec is None: choices = make_choices_explicit_error_term_mnl( state, diff --git a/activitysim/core/random.py b/activitysim/core/random.py index ccf1cdb55..e3ac6eac3 100644 --- a/activitysim/core/random.py +++ b/activitysim/core/random.py @@ -441,7 +441,9 @@ def gumbel_choice_positions_for_df( alt_nrs_df : pandas.DataFrame, optional DataFrame aligned to `utilities` whose values identify which dense alternative each utility column corresponds to. Use `MASKED_ALT_ID` (-999) for masked or - unavailable positions; any other negative value raises ValueError. + unavailable positions; any other negative value raises ValueError. A row that + is masked in every column is taken to mean the chooser has no alternative to + choose from, and returns position 0. n_rands : int, optional Number of EV1 draws to generate per chooser row. Required when `alt_nrs_df` is provided and may exceed the visible number of utility columns. @@ -479,6 +481,7 @@ def gumbel_choice_positions_for_df( f"{MASKED_ALT_ID} sentinel: {offenders}" ) masked = alt_nr_values == MASKED_ALT_ID + active_mask = ~masked safe_alt_nrs = np.where(masked, 0, alt_nr_values) else: if n_rands is None: @@ -487,7 +490,7 @@ def gumbel_choice_positions_for_df( raise ValueError( "n_rands must equal utilities.shape[1] when alt_nrs_df is omitted" ) - alt_nr_values = masked = safe_alt_nrs = None + alt_nr_values = masked = active_mask = safe_alt_nrs = None generators = self._generators_for_df(utilities) @@ -500,16 +503,21 @@ def gumbel_choice_positions_for_df( utility_row - np.log(-np.log(row_randoms)) ) else: - # Masked positions are set to -inf so they cannot win argmax, - # and the gumbel transform is skipped for them entirely. - row_mask = masked[row_num] - candidate_values = np.full(n_alts, -np.inf, dtype=np.float64) - active = ~row_mask - if active.any(): - candidate_values[active] = utility_row[active] - np.log( - -np.log(row_randoms[safe_alt_nrs[row_num, active]]) - ) - positions[row_num] = np.argmax(candidate_values) + # Masked positions can never be chosen, so apply the gumbel transform and argmax + # only to the active ones. flatnonzero returns ascending indices, so ties resolve + # to the lowest column position, which is what a full-width argmax would have done. + active = np.flatnonzero(active_mask[row_num]) + if active.size == 0: + # The chooser has no alternative available at all. Return the first + # column by convention. Note that logit.make_choices_utility_based filters + # out choosers with no available alternatives, so this case will only occur + # when explicitly allowed by the caller. + positions[row_num] = 0 + continue + gumbel = utility_row[active] - np.log( + -np.log(row_randoms[safe_alt_nrs[row_num, active]]) + ) + positions[row_num] = active[np.argmax(gumbel)] self.row_states.loc[utilities.index, "offset"] += n_rands return positions diff --git a/activitysim/core/simulate.py b/activitysim/core/simulate.py index dd5a22869..9e3f38337 100644 --- a/activitysim/core/simulate.py +++ b/activitysim/core/simulate.py @@ -1553,7 +1553,7 @@ def eval_nl( if state.settings.use_explicit_error_terms: raw_utilities = logit.validate_utils( - state, raw_utilities, allow_zero_probs=False, trace_label=trace_label + state, raw_utilities, trace_label=trace_label ) if custom_chooser: diff --git a/activitysim/core/test/test_logit.py b/activitysim/core/test/test_logit.py index e8a5b5143..3cee651d2 100644 --- a/activitysim/core/test/test_logit.py +++ b/activitysim/core/test/test_logit.py @@ -119,6 +119,82 @@ def test_validate_utils_allows_zero_probs(): assert validated.iloc[0, 1] == logit.UTIL_UNAVAILABLE +def _rng_state_for(choosers): + state = workflow.State().default_settings() + state.settings.skip_failed_choices = False + state.rng().set_base_seed(0) + state.rng().begin_step("test_step") + state.rng().add_channel("persons", choosers) + return state + + +def test_make_choices_utility_based_reports_when_no_alternative_available(): + # the counterpart of make_choices reporting bad probabilities: an argmax over a row + # with nothing available still returns a position, so it has to be reported here + choosers = pd.DataFrame(index=pd.Index([1, 2], name="person_id")) + utils = pd.DataFrame( + [[0.5, 1.0], [logit.UTIL_UNAVAILABLE, logit.UTIL_UNAVAILABLE]], + index=choosers.index, + ) + state = _rng_state_for(choosers) + + with pytest.raises(InvalidTravelError) as excinfo: + logit.make_choices_utility_based(state, utils) + + assert "no alternative is available" in str(excinfo.value) + + +def test_make_choices_utility_based_reports_raw_unavailable_values(): + # callers that reach here without validate_utils still carry raw values rather than + # UTIL_UNAVAILABLE, so the check tests the UTIL_MIN threshold rather than a row sum + choosers = pd.DataFrame(index=pd.Index([1], name="person_id")) + utils = pd.DataFrame( + [[logit.UTIL_MIN - 1.0, logit.UTIL_MIN - 2.0]], index=choosers.index + ) + state = _rng_state_for(choosers) + + with pytest.raises(InvalidTravelError, match="no alternative is available"): + logit.make_choices_utility_based(state, utils) + + +def test_make_choices_utility_based_allows_zero_probs(): + # callers that have already sanctioned the situation opt out, as with + # allow_bad_probs in make_choices + choosers = pd.DataFrame(index=pd.Index([1], name="person_id")) + utils = pd.DataFrame( + [[logit.UTIL_UNAVAILABLE, logit.UTIL_UNAVAILABLE]], index=choosers.index + ) + state = _rng_state_for(choosers) + + choices, _ = logit.make_choices_utility_based(state, utils, allow_zero_probs=True) + + assert len(choices) == 1 + + +def test_make_choices_utility_based_reports_fully_masked_alt_nrs_row(): + # the alt_nrs_df case needs no separate handling: padding columns carry a utility + # below UTIL_MIN, so a row that is masked everywhere is also a row with nothing + # available and trips the same check + choosers = pd.DataFrame(index=pd.Index([1, 2], name="person_id")) + utils = pd.DataFrame( + [[0.5, 1.0], [logit.UTIL_UNAVAILABLE, logit.UTIL_UNAVAILABLE]], + index=choosers.index, + ) + alt_nrs_df = pd.DataFrame( + [[0, 1], [random.MASKED_ALT_ID, random.MASKED_ALT_ID]], + index=choosers.index, + ) + state = _rng_state_for(choosers) + + with pytest.raises(InvalidTravelError, match="no alternative is available"): + logit.make_choices_utility_based( + state, + utils, + alts_context=AltsContext.from_num_alts(2), + alt_nrs_df=alt_nrs_df, + ) + + # # `utils_to_probs` Tests # diff --git a/activitysim/core/test/test_random.py b/activitysim/core/test/test_random.py index d5f84bd12..1cda242f6 100644 --- a/activitysim/core/test/test_random.py +++ b/activitysim/core/test/test_random.py @@ -330,3 +330,75 @@ def test_gumbel_choice_positions_for_df_matches_dense_alt_mapping(): npt.assert_array_equal(observed_positions, expected_positions) npt.assert_allclose(next_random_after_fused, next_random_after_materialized) + + +def test_gumbel_choice_positions_for_df_masked_columns_never_win(): + # padded columns carry a high utility here, so if they were eligible they would + # win every argmax; only the single active column of each row may be returned + persons = pd.DataFrame( + {"household_id": [1, 1, 1]}, + index=pd.Index([41, 42, 43], name="person_id"), + ) + utilities = pd.DataFrame( + [[0.0, 99.0, 99.0], [99.0, 0.0, 99.0], [99.0, 99.0, 0.0]], + index=persons.index, + ) + alt_nrs_df = pd.DataFrame( + [ + [0, random.MASKED_ALT_ID, random.MASKED_ALT_ID], + [random.MASKED_ALT_ID, 1, random.MASKED_ALT_ID], + [random.MASKED_ALT_ID, random.MASKED_ALT_ID, 2], + ], + index=persons.index, + ) + + rng = random.Random() + rng.set_base_seed(0) + rng.begin_step("test_step") + rng.add_channel("persons", persons) + positions = rng.gumbel_choice_positions_for_df( + utilities, alt_nrs_df=alt_nrs_df, n_rands=3 + ) + rng.end_step("test_step") + + npt.assert_array_equal(positions, [0, 1, 2]) + + +def test_gumbel_choice_positions_for_df_fully_masked_row_falls_back_to_first_column(): + # MASKED_ALT_ID marks padded *or unavailable* slots, so an all-masked row means the + # chooser has no alternative available. That returns position 0, mirroring the Monte + # Carlo path's probs.loc[zero_probs, 0] = 1.0, and it must not disturb the choice or + # the random number stream of any other chooser. + persons = pd.DataFrame( + {"household_id": [1, 1, 1]}, + index=pd.Index([51, 52, 53], name="person_id"), + ) + utilities = pd.DataFrame( + [[2.0, 1.0], [0.3, 1.2], [0.7, 0.4]], + index=persons.index, + ) + all_active = pd.DataFrame([[0, 2], [0, 2], [0, 2]], index=persons.index) + with_masked_row = pd.DataFrame( + [[0, 2], [random.MASKED_ALT_ID, random.MASKED_ALT_ID], [0, 2]], + index=persons.index, + ) + + def run(alt_nrs_df): + rng = random.Random() + rng.set_base_seed(0) + rng.begin_step("test_step") + rng.add_channel("persons", persons) + positions = rng.gumbel_choice_positions_for_df( + utilities, alt_nrs_df=alt_nrs_df, n_rands=3 + ) + following = rng.random_for_df(persons) + rng.end_step("test_step") + return positions, following + + baseline_positions, baseline_following = run(all_active) + masked_positions, masked_following = run(with_masked_row) + + assert masked_positions[1] == 0 + npt.assert_array_equal(masked_positions[[0, 2]], baseline_positions[[0, 2]]) + # the masked row still consumes its n_rands draws, so offsets stay aligned + npt.assert_allclose(masked_following, baseline_following) From 035de9f869a3e59880054182d698f134eb8ac64e Mon Sep 17 00:00:00 2001 From: Jan Zill Date: Mon, 27 Jul 2026 08:46:31 +1000 Subject: [PATCH 6/6] harmonize sample choice maker naming --- activitysim/core/interaction_sample.py | 6 +-- .../core/test/test_interaction_sample.py | 42 +++++++++---------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/activitysim/core/interaction_sample.py b/activitysim/core/interaction_sample.py index dca38442b..84acadefb 100644 --- a/activitysim/core/interaction_sample.py +++ b/activitysim/core/interaction_sample.py @@ -140,7 +140,7 @@ def _poisson_fallback_positions( which a data-dependent retry or redraw scheme cannot do. Because the fallback set is a deterministic function of the probabilities, the probability that an alternative ends up in the returned choice set still has an exact closed form - (see `_poisson_sample_alternatives`). + (see `make_sample_choices_poisson`). """ k = min(sample_size, probs_values.shape[1]) # stable sort of the negated probabilities gives descending probability order @@ -204,7 +204,7 @@ def make_sample_choices_eet( return choices_df -def _poisson_sample_alternatives( +def make_sample_choices_poisson( chunk_sizer: ChunkSizer, probs: pd.DataFrame, alternatives: pd.DataFrame, @@ -912,7 +912,7 @@ def _interaction_sample( n_total_alts=n_total_alts, ) else: # sampling_method == "poisson" - choices_df = _poisson_sample_alternatives( + choices_df = make_sample_choices_poisson( chunk_sizer, probs, alternatives, diff --git a/activitysim/core/test/test_interaction_sample.py b/activitysim/core/test/test_interaction_sample.py index b520c264c..2d929053d 100644 --- a/activitysim/core/test/test_interaction_sample.py +++ b/activitysim/core/test/test_interaction_sample.py @@ -765,7 +765,7 @@ def test_poisson_fallback_positions_breaks_ties_by_column_and_caps_at_alt_count( ) -def test_poisson_sample_alternatives_returns_expected_frames(): +def test_make_sample_choices_poisson_returns_expected_frames(): probs = pd.DataFrame( [ [0.20, 0.60, 0.10, 0.05], @@ -788,14 +788,14 @@ def test_poisson_sample_alternatives_returns_expected_frames(): ) state = _DummyState(_SequentialDummyRng([draws])) - choices_df = interaction_sample._poisson_sample_alternatives( + 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_poisson_sample_alternatives_returns_expected_frames", + trace_label="test_make_sample_choices_poisson_returns_expected_frames", ) expected = _reference_poisson_choices_df( @@ -807,7 +807,7 @@ def test_poisson_sample_alternatives_returns_expected_frames(): assert choices_df.loc[choices_df.person_id == 17, "alt_id"].tolist() == [100, 700] -def test_poisson_sample_alternatives_consumes_no_extra_randoms_on_empty_draw(): +def test_make_sample_choices_poisson_consumes_no_extra_randoms_on_empty_draw(): # the fallback must not draw again: _SequentialDummyRng raises IndexError if the # sampler asks for a second block, so a single draw array is the assertion here probs = pd.DataFrame( @@ -820,14 +820,14 @@ def test_poisson_sample_alternatives_consumes_no_extra_randoms_on_empty_draw(): fail_draw = np.array([[0.99, 0.99, 0.99]], dtype=np.float64) state = _DummyState(_SequentialDummyRng([fail_draw])) - choices_df = interaction_sample._poisson_sample_alternatives( + 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_poisson_sample_alternatives_consumes_no_extra_randoms_on_empty_draw", + trace_label="test_make_sample_choices_poisson_consumes_no_extra_randoms_on_empty_draw", ) # the two most likely alternatives, reported at q_i + P0 @@ -847,7 +847,7 @@ def test_poisson_sample_alternatives_consumes_no_extra_randoms_on_empty_draw(): pd.testing.assert_frame_equal(choices_df, expected) -def test_poisson_sample_alternatives_reported_prob_is_total_inclusion_probability(): +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 # fallback. Every chooser is identical, so the empirical inclusion rate across @@ -872,14 +872,14 @@ def test_poisson_sample_alternatives_reported_prob_is_total_inclusion_probabilit draws = np.random.default_rng(20260726).random((n_choosers, n_alts)) state = _DummyState(_SequentialDummyRng([draws])) - choices_df = interaction_sample._poisson_sample_alternatives( + 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_poisson_sample_alternatives_reported_prob_is_total_inclusion_probability", + trace_label="test_make_sample_choices_poisson_reported_prob_is_total_inclusion_probability", ) # identical choosers must get an identical reported prob per alternative, whether @@ -897,7 +897,7 @@ def test_poisson_sample_alternatives_reported_prob_is_total_inclusion_probabilit np.testing.assert_allclose(empirical, expected_reported, atol=0.005) -def test_poisson_sample_alternatives_repeat_alignment_chooser_dominant_heterogeneity(): +def test_repeat_alignment_chooser_heterogeneity(): # Edge case: utilities are close across alternatives but vary strongly by chooser. # This checks that the flattened Poisson result keeps chooser/prob alignment. chooser_index = pd.Index([101, 102, 103, 104, 105, 106], name="person_id") @@ -938,7 +938,7 @@ def test_poisson_sample_alternatives_repeat_alignment_chooser_dominant_heterogen trace_choosers=choosers, ) - out = interaction_sample._poisson_sample_alternatives( + out = interaction_sample.make_sample_choices_poisson( chunk_sizer=_DummyChunkSizer(), probs=probs, alternatives=alternatives, @@ -955,7 +955,7 @@ def test_poisson_sample_alternatives_repeat_alignment_chooser_dominant_heterogen pd.testing.assert_frame_equal(out.reset_index(drop=True), expected) -def test_poisson_sample_alternatives_matches_materialized_path(): +def test_make_sample_choices_poisson_matches_materialized_path(): chooser_index = pd.Index([201, 202, 203], name="person_id") choosers = pd.DataFrame(index=chooser_index) alternatives = pd.DataFrame(index=pd.Index([10, 11, 12, 13], name="alt_id")) @@ -983,7 +983,7 @@ def test_poisson_sample_alternatives_matches_materialized_path(): trace_choosers=choosers, ) - out = interaction_sample._poisson_sample_alternatives( + out = interaction_sample.make_sample_choices_poisson( chunk_sizer=_DummyChunkSizer(), probs=probs, alternatives=alternatives, @@ -1120,7 +1120,7 @@ def test_make_sample_choices_eet_stable_alt_mapping_matches_materialized_path(): pd.testing.assert_frame_equal(out.reset_index(drop=True), expected) -def test_poisson_sample_alternatives_stable_alt_mapping_matches_materialized_path(): +def test_make_sample_choices_poisson_stable_alt_mapping_matches_materialized_path(): chooser_index = pd.Index([311, 312], name="person_id") choosers = pd.DataFrame(index=chooser_index) alternatives = pd.DataFrame(index=pd.Index([10, 12, 14], name="alt_id")) @@ -1144,19 +1144,19 @@ def test_poisson_sample_alternatives_stable_alt_mapping_matches_materialized_pat state, utilities, allow_zero_probs=False, - trace_label="test_poisson_sample_alternatives_stable_alt_mapping_matches_materialized_path", + trace_label="test_make_sample_choices_poisson_stable_alt_mapping_matches_materialized_path", overflow_protection=True, trace_choosers=choosers, ) - out = interaction_sample._poisson_sample_alternatives( + out = 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_poisson_sample_alternatives_stable_alt_mapping_matches_materialized_path", + trace_label="test_make_sample_choices_poisson_stable_alt_mapping_matches_materialized_path", stable_alt_positions=stable_alt_positions, n_total_alts=n_total_alts, ) @@ -1172,7 +1172,7 @@ def test_poisson_sample_alternatives_stable_alt_mapping_matches_materialized_pat pd.testing.assert_frame_equal(out.reset_index(drop=True), expected) -def test_poisson_sample_alternatives_falls_back_to_most_likely_alternatives(): +def test_make_sample_choices_poisson_falls_back_to_most_likely_alternatives(): chooser_index = pd.Index([301, 302], name="person_id") choosers = pd.DataFrame(index=chooser_index) alternatives = pd.DataFrame(index=pd.Index([10, 12, 14], name="alt_id")) @@ -1188,19 +1188,19 @@ def test_poisson_sample_alternatives_falls_back_to_most_likely_alternatives(): state, utilities, allow_zero_probs=False, - trace_label="test_falls_back_to_most_likely_alternatives", + trace_label="test_make_sample_choices_poisson_falls_back_to_most_likely_alternatives", overflow_protection=True, trace_choosers=choosers, ) - out = interaction_sample._poisson_sample_alternatives( + out = 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_falls_back_to_most_likely_alternatives", + trace_label="test_make_sample_choices_poisson_falls_back_to_most_likely_alternatives", ) # neither chooser sampled anything, so both take their two most likely