diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 12d970e18..587f957ec 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -14,7 +14,7 @@ repos:
     args: ["--profile", "black", "--filter-files"]
 
 - repo: https://github.com/psf/black
-  rev: 22.6.0
+  rev: 22.12.0
   hooks:
     - id: black
       language_version: python3.13
diff --git a/activitysim/abm/models/location_choice.py b/activitysim/abm/models/location_choice.py
index 7f032a8ae..8a1ea9611 100644
--- a/activitysim/abm/models/location_choice.py
+++ b/activitysim/abm/models/location_choice.py
@@ -17,6 +17,7 @@
 )
 from activitysim.core.interaction_sample import interaction_sample
 from activitysim.core.interaction_sample_simulate import interaction_sample_simulate
+from activitysim.core.logit import AltsContext
 from activitysim.core.util import reindex
 from activitysim.core.exceptions import DuplicateWorkflowTableError
 
@@ -603,6 +604,7 @@ def run_location_simulate(
     chunk_tag,
     trace_label,
     skip_choice=False,
+    alts_context: AltsContext | None = None,
 ):
     """
     run location model on location_sample annotated with mode_choice logsum
@@ -712,6 +714,7 @@ def run_location_simulate(
         compute_settings=model_settings.compute_settings.subcomponent_settings(
             "simulate"
         ),
+        alts_context=alts_context,
     )
 
     if not want_logsums:
@@ -788,6 +791,13 @@ def run_location_choice(
         if choosers.shape[0] == 0:
             logger.info(f"{trace_label} skipping segment {segment_name}: no choosers")
             continue
+        # using land use rather than size terms in case something goes 0 base -> nonzero project, double
+        # check if that would be in dest_size_terms as a zero
+        alts_context = AltsContext.from_series(
+            dest_size_terms.index
+        )  # index zone_id, not ALT_DEST_COL_NAME
+        # assumes that dest_size_terms will always contain zeros for non-attractive zones, i.e. it will have the
+        # same length as land_use
 
         # - location_sample
         location_sample_df = run_location_sample(
@@ -841,6 +851,7 @@ def run_location_choice(
                 trace_label, "simulate.%s" % segment_name
             ),
             skip_choice=skip_choice,
+            alts_context=alts_context,
         )
 
         if estimator:
@@ -1019,6 +1030,10 @@ def iterate_location_choice(
     ) = None  # initialize to None, will be populated in first iteration
 
     for iteration in range(1, max_iterations + 1):
+        # Force reset the setting at the start of each shadow price iteration for consistency
+        logger.info("Resetting random number seeds for iteration {}".format(iteration))
+        state.get_rn_generator().end_step(trace_label)
+        state.get_rn_generator().begin_step(trace_label)
         persons_merged_df_ = persons_merged_df.copy()
 
         if spc.use_shadow_pricing and iteration > 1:
diff --git a/activitysim/abm/models/parking_location_choice.py b/activitysim/abm/models/parking_location_choice.py
index 32f3aabee..743b82261 100644
--- a/activitysim/abm/models/parking_location_choice.py
+++ b/activitysim/abm/models/parking_location_choice.py
@@ -21,6 +21,7 @@
 from activitysim.core.configuration.base import PreprocessorSettings
 from activitysim.core.configuration.logit import LogitComponentSettings
 from activitysim.core.interaction_sample_simulate import interaction_sample_simulate
+from activitysim.core.logit import AltsContext
 from activitysim.core.tracing import print_elapsed_time
 from activitysim.core.util import assign_in_place, drop_unused_columns
 from activitysim.core.exceptions import DuplicateWorkflowTableError
@@ -112,6 +113,7 @@ def parking_destination_simulate(
     chunk_size,
     trace_hh_id,
     trace_label,
+    alts_context: AltsContext | None = None,
 ):
     """
     Chose destination from destination_sample (with od_logsum and dp_logsum columns added)
@@ -150,6 +152,7 @@ def parking_destination_simulate(
         trace_label=trace_label,
         trace_choice_name="parking_loc",
         explicit_chunk_size=model_settings.explicit_chunk,
+        alts_context=alts_context,
     )
 
     # drop any failed zero_prob destinations
@@ -211,6 +214,9 @@ def choose_parking_location(
     )
     destination_sample.index = np.repeat(trips.index.values, len(alternatives))
     destination_sample.index.name = trips.index.name
+    # using destination_sample would also be right because destination_sample isn't a sample here,
+    # but that could change
+    alts_context = AltsContext.from_series(alternatives[alt_dest_col_name])
 
     destinations = parking_destination_simulate(
         state,
@@ -223,6 +229,7 @@ def choose_parking_location(
         chunk_size=chunk_size,
         trace_hh_id=trace_hh_id,
         trace_label=trace_label,
+        alts_context=alts_context,
     )
 
     if want_sample_table:
diff --git a/activitysim/abm/models/trip_destination.py b/activitysim/abm/models/trip_destination.py
index 853cfc35e..3ad35470c 100644
--- a/activitysim/abm/models/trip_destination.py
+++ b/activitysim/abm/models/trip_destination.py
@@ -32,6 +32,7 @@
 from activitysim.core.configuration.logit import LocationComponentSettings
 from activitysim.core.interaction_sample import interaction_sample
 from activitysim.core.interaction_sample_simulate import interaction_sample_simulate
+from activitysim.core.logit import AltsContext
 from activitysim.core.skim_dictionary import DataFrameMatrix
 from activitysim.core.tracing import print_elapsed_time
 from activitysim.core.util import assign_in_place, reindex
@@ -950,6 +951,7 @@ def trip_destination_simulate(
     skim_hotel,
     estimator,
     trace_label,
+    alts_context: AltsContext | None = None,
 ):
     """
     Chose destination from destination_sample (with od_logsum and dp_logsum columns added)
@@ -1036,6 +1038,7 @@ def trip_destination_simulate(
         trace_choice_name="trip_dest",
         estimator=estimator,
         explicit_chunk_size=model_settings.explicit_chunk,
+        alts_context=alts_context,
     )
 
     if not want_logsums:
@@ -1126,7 +1129,10 @@ def choose_trip_destination(
         destination_sample["dp_logsum"] = 0.0
 
     t0 = print_elapsed_time("%s.compute_logsums" % trace_label, t0, debug=True)
-
+    alt_dest_col_name = model_settings.ALT_DEST_COL_NAME
+    alts = alternatives.index
+    assert alts.name == alt_dest_col_name
+    alts_context = AltsContext.from_series(alts)
     destinations = trip_destination_simulate(
         state,
         primary_purpose=primary_purpose,
@@ -1138,6 +1144,7 @@ def choose_trip_destination(
         skim_hotel=skim_hotel,
         estimator=estimator,
         trace_label=trace_label,
+        alts_context=alts_context,
     )
 
     dropped_trips = ~trips.index.isin(destinations.index)
diff --git a/activitysim/abm/models/trip_scheduling_choice.py b/activitysim/abm/models/trip_scheduling_choice.py
index 81d908ef1..3eb695feb 100644
--- a/activitysim/abm/models/trip_scheduling_choice.py
+++ b/activitysim/abm/models/trip_scheduling_choice.py
@@ -20,6 +20,7 @@
 )
 from activitysim.core.configuration.logit import LogitComponentSettings
 from activitysim.core.interaction_sample_simulate import _interaction_sample_simulate
+from activitysim.core.logit import AltsContext
 from activitysim.core.skim_dataset import SkimDataset
 from activitysim.core.skim_dictionary import SkimDict
 
@@ -314,6 +315,9 @@ def run_trip_scheduling_choice(
                 estimator=None,
                 chunk_sizer=chunk_sizer,
                 compute_settings=model_settings.compute_settings,
+                alts_context=AltsContext(
+                    schedules[SCHEDULE_ID].min(), schedules[SCHEDULE_ID].max()
+                ),
             )
 
             assert len(choices.index) == len(choosers.index)
diff --git a/activitysim/abm/tables/shadow_pricing.py b/activitysim/abm/tables/shadow_pricing.py
index 04c5eafc2..5e0e1fc70 100644
--- a/activitysim/abm/tables/shadow_pricing.py
+++ b/activitysim/abm/tables/shadow_pricing.py
@@ -171,6 +171,7 @@ def __init__(
 
         self.num_processes = num_processes
         self.use_shadow_pricing = bool(state.settings.use_shadow_pricing)
+        logger.info("Use shadow pricing: %s", self.use_shadow_pricing)
         self.saved_shadow_price_file_path = (
             None  # set by read_saved_shadow_prices if loaded
         )
diff --git a/activitysim/core/chunk.py b/activitysim/core/chunk.py
index f0074683d..019f6abae 100644
--- a/activitysim/core/chunk.py
+++ b/activitysim/core/chunk.py
@@ -1244,6 +1244,7 @@ def adaptive_chunked_choosers(
             chunk_size = math.ceil(num_choosers * explicit_chunk_size)
         else:
             chunk_size = math.ceil(explicit_chunk_size / num_processes)
+        logger.info(f"Running {trace_label} with {chunk_size} choosers")
     elif chunk_size is None:
         chunk_size = state.settings.chunk_size
 
diff --git a/activitysim/core/interaction_sample.py b/activitysim/core/interaction_sample.py
index 93834c690..e2cc51001 100644
--- a/activitysim/core/interaction_sample.py
+++ b/activitysim/core/interaction_sample.py
@@ -3,10 +3,10 @@
 from __future__ import annotations
 
 import logging
+import typing
 
 import numpy as np
 import pandas as pd
-
 from activitysim.core import (
     chunk,
     estimation,
@@ -17,16 +17,35 @@
     util,
     workflow,
 )
+from activitysim.core.chunk import ChunkSizer
 from activitysim.core.configuration.base import ComputeSettings
 from activitysim.core.exceptions import SegmentedSpecificationError
 from activitysim.core.skim_dataset import DatasetWrapper
 from activitysim.core.skim_dictionary import SkimWrapper
 
+if typing.TYPE_CHECKING:
+    from activitysim.core.random import Random
+
 logger = logging.getLogger(__name__)
 
 DUMP = False
 
 
+def _poisson_sample_alternatives_inner(
+    alternative_count: int,
+    probs: pd.DataFrame,
+    poisson_inclusion_probs: pd.DataFrame,
+    rng: Random,
+    trace_label: str | None,
+    chunk_sizer: ChunkSizer,
+) -> pd.DataFrame:
+    rands = rng.random_for_df(probs, n=alternative_count)
+    chunk_sizer.log_df(trace_label, "rands", rands)
+    sampled_mask = rands < poisson_inclusion_probs
+    sampled_results = poisson_inclusion_probs.where(sampled_mask)
+    return sampled_results
+
+
 def make_sample_choices_utility_based(
     state: workflow.State,
     choosers,
@@ -36,8 +55,8 @@ def make_sample_choices_utility_based(
     alternative_count,
     alt_col_name,
     allow_zero_probs,
-    trace_label,
-    chunk_sizer,
+    trace_label: str,
+    chunk_sizer: ChunkSizer,
 ):
     assert isinstance(utilities, pd.DataFrame)
     assert utilities.shape == (len(choosers), alternative_count)
@@ -60,32 +79,6 @@ def make_sample_choices_utility_based(
 
     utils_array = utilities.to_numpy()
     chunk_sizer.log_df(trace_label, "utils_array", utils_array)
-    chosen_destinations = []
-
-    rands = state.get_rn_generator().gumbel_for_df(utilities, n=alternative_count)
-    chunk_sizer.log_df(trace_label, "rands", rands)
-
-    # TODO-EET [janzill Jun2022]: using for-loop to keep memory usage low, an array of dimension
-    #  (len(choosers), alternative_count, sample_size) can get very large. Probably better to
-    #  use chunking for this.
-    for i in range(sample_size):
-        # created this once for memory logging
-        if i > 0:
-            rands = state.get_rn_generator().gumbel_for_df(
-                utilities, n=alternative_count
-            )
-        chosen_destinations.append(np.argmax(utils_array + rands, axis=1))
-    chosen_destinations = np.concatenate(chosen_destinations, axis=0)
-
-    chunk_sizer.log_df(trace_label, "chosen_destinations", chosen_destinations)
-
-    del utils_array
-    chunk_sizer.log_df(trace_label, "utils_array", None)
-    del rands
-    chunk_sizer.log_df(trace_label, "rands", None)
-
-    chooser_idx = np.tile(np.arange(utilities.shape[0]), sample_size)
-    chunk_sizer.log_df(trace_label, "chooser_idx", chooser_idx)
 
     probs = logit.utils_to_probs(
         state,
@@ -95,28 +88,88 @@ def make_sample_choices_utility_based(
         overflow_protection=not allow_zero_probs,
         trace_choosers=choosers,
     )
-    chunk_sizer.log_df(trace_label, "probs", probs)
+    inclusion_probs, sampled_alternatives = _poisson_sample_alternatives(
+        alternative_count, chunk_sizer, probs, sample_size, state, trace_label
+    )
 
-    choices_df = pd.DataFrame(
-        {
-            alt_col_name: alternatives.index.values[chosen_destinations],
-            "prob": probs.to_numpy()[chooser_idx, chosen_destinations],
-            choosers.index.name: choosers.index.values[chooser_idx],
-        }
+    # Stack removes the NaNs (the ones that weren't sampled)
+    # and gives us a multi-index of (person_id, alt_id)
+    choices_df = (
+        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"])
     )
-    chunk_sizer.log_df(trace_label, "choices_df", choices_df)
 
-    del chooser_idx
-    chunk_sizer.log_df(trace_label, "chooser_idx", None)
-    del chosen_destinations
-    chunk_sizer.log_df(trace_label, "chosen_destinations", None)
-    del probs
-    chunk_sizer.log_df(trace_label, "probs", None)
+    # Here we return the inclusion probabilities i.e. the true probability of being sampled and (ab)use the fact
+    # that pick_count=1 by definition and ln(1)=0 and recover the standard sample correction term.
+    # In non-Poisson sampling, we would return the probs of sampling an alternative once
+    # and the sampling correction factor np.log(df.pick_count/df.prob) is applied to the simulate utilities.
+    # TODO is it safe change the meaning of df.prob, given it's referenced in expression csvs?
+    #   (but the alternative is to update all the expression CSV for sampling?)
+    return choices_df, inclusion_probs
 
-    # handing this off to caller
-    chunk_sizer.log_df(trace_label, "choices_df", None)
 
-    return choices_df
+def _poisson_sample_alternatives(
+    alternative_count,
+    chunk_sizer: ChunkSizer,
+    probs: pd.DataFrame,
+    sample_size,
+    state: workflow.State,
+    trace_label: str,
+) -> tuple[pd.DataFrame, pd.DataFrame]:
+    # compute the inclusion probability as the reciprocal of alt never being drawn
+    #  -- these are common, so compute once upfront
+    exclusion_probs = (1 - probs) ** sample_size
+    inclusion_probs = 1 - exclusion_probs
+
+    n = 0
+    probs_subset = probs
+    inclusion_probs_subset = inclusion_probs
+    sampled_alternatives = pd.DataFrame(
+        0.0, index=inclusion_probs.index, columns=inclusion_probs.columns
+    )
+    while True:
+        sampled_results_subset = _poisson_sample_alternatives_inner(
+            alternative_count,
+            probs_subset,
+            inclusion_probs_subset,
+            state.get_rn_generator(),
+            trace_label,
+            chunk_sizer,
+        )
+        no_alts_sampled_mask = sampled_results_subset.isna().all(axis=1)
+        alts_with_sampled_alternatives = sampled_results_subset[~no_alts_sampled_mask]
+        sampled_alternatives.loc[
+            alts_with_sampled_alternatives.index, :
+        ] = alts_with_sampled_alternatives
+        if no_alts_sampled_mask.any():
+            # TODO if this happens in base but the project case is such that something is picked, random numbers won't
+            #  be consistent - we're asserting that this is very rare models where the sample size is not too small
+            logger.info(f"Poisson sampling of alternatives failed with {n=}, retrying")
+            # TODO put this behind a debug guard, because it will be slow
+            logger.info(
+                f"Sampled size was {sample_size}, poisson method mean expected sample size was {inclusion_probs.sum(axis=1).mean():.1f}, actual sampled mean was {(sampled_alternatives > 0).sum(axis=1).mean():.1f} and highest zero selection prob was {(exclusion_probs).product(axis=1).max():.2g}"
+            )
+            probs_subset = probs[no_alts_sampled_mask]
+            inclusion_probs_subset = inclusion_probs[no_alts_sampled_mask]
+
+        else:  # All alternatives are fine
+            break
+
+        n += 1
+        if n == 10:
+            choosers_no_alts_sampled = sampled_results_subset[no_alts_sampled_mask]
+            msg = (
+                f"Poisson choice set sampling failed after 10 attempts for these cases:\n"
+                f"{choosers_no_alts_sampled}\n{probs_subset}"
+            )
+            raise ValueError(msg)
+
+    chunk_sizer.log_df(trace_label, "sampled_alternatives", sampled_alternatives)
+
+    return inclusion_probs, sampled_alternatives
 
 
 def make_sample_choices(
@@ -227,7 +280,7 @@ def _interaction_sample(
     locals_d=None,
     trace_label=None,
     zone_layer=None,
-    chunk_sizer=None,
+    chunk_sizer: ChunkSizer | None = None,
     compute_settings: ComputeSettings | None = None,
 ):
     """
@@ -292,6 +345,11 @@ def _interaction_sample(
         pick_count : int
             number of duplicate picks for chooser, alt
     """
+    assert (
+        chunk_sizer is not None
+    ), "chunk_sizer cannot be None but old nullable signature is preserved"
+    # TODO it's probably safe to reorder these arguments to make chunk_sizer mandatory since
+    #   _interaction_sample is private?
 
     have_trace_targets = state.tracing.has_trace_targets(choosers)
     trace_ids = None
@@ -812,7 +870,13 @@ def interaction_sample(
         assert choosers.index.is_monotonic_increasing
 
     # FIXME - legacy logic - not sure this is needed or even correct?
-    sample_size = min(sample_size, len(alternatives.index))
+    if not state.settings.use_explicit_error_terms:
+        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)
+
     logger.info(f" --- interaction_sample sample size = {sample_size}")
 
     result_list = []
diff --git a/activitysim/core/interaction_sample_simulate.py b/activitysim/core/interaction_sample_simulate.py
index e73f64f4f..3f1ce1896 100644
--- a/activitysim/core/interaction_sample_simulate.py
+++ b/activitysim/core/interaction_sample_simulate.py
@@ -9,6 +9,7 @@
 
 from activitysim.core import chunk, interaction_simulate, logit, tracing, util, workflow
 from activitysim.core.configuration.base import ComputeSettings
+from activitysim.core.logit import AltsContext
 from activitysim.core.simulate import set_skim_wrapper_targets
 from activitysim.core.exceptions import SegmentedSpecificationError
 
@@ -34,6 +35,7 @@ def _interaction_sample_simulate(
     *,
     chunk_sizer: chunk.ChunkSizer,
     compute_settings: ComputeSettings | None = None,
+    alts_context: AltsContext | None = None,
 ):
     """
     Run a MNL simulation in the situation in which alternatives must
@@ -220,9 +222,6 @@ def _interaction_sample_simulate(
     )
     chunk_sizer.log_df(trace_label, "interaction_utilities", interaction_utilities)
 
-    del interaction_df
-    chunk_sizer.log_df(trace_label, "interaction_df", None)
-
     if have_trace_targets:
         state.tracing.trace_interaction_eval_results(
             trace_eval_results,
@@ -258,26 +257,37 @@ def _interaction_sample_simulate(
     # (we want to insert dummy utilities at the END of the list of alternative utilities)
     # inserts is a list of the indices at which we want to do the insertions
     inserts = np.repeat(last_row_offsets, max_sample_count - sample_counts)
-
+    del last_row_offsets
     del sample_counts
     chunk_sizer.log_df(trace_label, "sample_counts", None)
 
     # insert the zero-prob utilities to pad each alternative set to same size
     padded_utilities = np.insert(interaction_utilities.utility.values, inserts, -999)
+    padded_alt_nrs = np.insert(interaction_df[choice_column], inserts, -999)
     chunk_sizer.log_df(trace_label, "padded_utilities", padded_utilities)
     del inserts
-
+    del interaction_df
     del interaction_utilities
+    # TODO EET: chunk_sizer logging could be more precise
+    chunk_sizer.log_df(trace_label, "interaction_df", None)
     chunk_sizer.log_df(trace_label, "interaction_utilities", None)
 
     # reshape to array with one row per chooser, one column per alternative
     padded_utilities = padded_utilities.reshape(-1, max_sample_count)
+    padded_alt_nrs = padded_alt_nrs.reshape(-1, max_sample_count)
 
     # convert to a dataframe with one row per chooser and one column per alternative
     utilities_df = pd.DataFrame(padded_utilities, index=choosers.index)
+    # alt_nrs_df has columns for each alt in the choice set, with values indicating which alt_id
+    # they correspond to (as opposed to the 0-n index implied by the column number).
+    if alts_context is not None:
+        alt_nrs_df = pd.DataFrame(padded_alt_nrs, index=choosers.index)
+    else:
+        alt_nrs_df = None  # if we don't provide the number of dense alternatives, assume that we'll use the old approach
     chunk_sizer.log_df(trace_label, "utilities_df", utilities_df)
 
     del padded_utilities
+    del padded_alt_nrs
     chunk_sizer.log_df(trace_label, "padded_utilities", None)
 
     if have_trace_targets:
@@ -320,7 +330,12 @@ def _interaction_sample_simulate(
         # 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
         positions, rands = logit.make_choices_utility_based(
-            state, utilities_df, trace_label=trace_label, trace_choosers=choosers
+            state,
+            utilities_df,
+            trace_label=trace_label,
+            trace_choosers=choosers,
+            alts_context=alts_context,
+            alt_nrs_df=alt_nrs_df,
         )
 
         del utilities_df
@@ -451,6 +466,7 @@ def interaction_sample_simulate(
     skip_choice=False,
     explicit_chunk_size=0,
     *,
+    alts_context: AltsContext | None = None,
     compute_settings: ComputeSettings | None = None,
 ):
     """
@@ -496,6 +512,12 @@ def interaction_sample_simulate(
     explicit_chunk_size : float, optional
         If > 0, specifies the chunk size to use when chunking the interaction
         simulation. If < 1, specifies the fraction of the total number of choosers.
+    alts_context: int, optional
+        The number of alternatives available in the choice set in the absense of sampling.
+        This is used with EET simulation to ensure consistent random numbers across the whole alternative set
+        ( as the sampled set may change between base and project). When not provided,
+        the fallback approach is used which may result in frozen error terms being applied to the wrong alternatives
+        if the choice set changes.
 
     Returns
     -------
@@ -551,6 +573,7 @@ def interaction_sample_simulate(
             skip_choice,
             chunk_sizer=chunk_sizer,
             compute_settings=compute_settings,
+            alts_context=alts_context,
         )
 
         result_list.append(choices)
diff --git a/activitysim/core/logit.py b/activitysim/core/logit.py
index 0030168bb..aa1ed5572 100644
--- a/activitysim/core/logit.py
+++ b/activitysim/core/logit.py
@@ -4,6 +4,8 @@
 
 import logging
 import warnings
+from dataclasses import dataclass
+from typing import Union
 
 import numpy as np
 import pandas as pd
@@ -344,12 +346,85 @@ def utils_to_probs(
     return probs
 
 
+FREEZE_RANDOM_NUMBERS_FOR_DENSE_ALTERNATIVE_SET = True
+
+
+@dataclass
+class AltsContext:
+    """Representation of the alternatives without carrying around that full array."""
+
+    min_alt_id: int
+    max_alt_id: int
+
+    def __post_init__(self):
+        # e.g. for zero based zones max_alt_id = n_alts - 1
+        # but for 1 based zones, we don't need to add extra padding
+        self.n_rands_to_sample = max(self.max_alt_id, self.n_alts_to_cover_max_id)
+
+    @classmethod
+    def from_series(cls, ser: Union[pd.Series, pd.Index]) -> "AltsContext":
+        min_alt_id = ser.min()
+        max_alt_id = ser.max()
+        return cls(min_alt_id, max_alt_id)
+
+    @classmethod
+    def from_num_alts(cls, num_alts: int, zero_based: bool = True) -> "AltsContext":
+        if zero_based:
+            offset = -1
+        else:
+            offset = 0
+        return cls(min_alt_id=1 + offset, max_alt_id=num_alts + offset)
+
+    @property
+    def n_alts_to_cover_max_id(self) -> int:
+        """If zones were non-consecutive, this could be a big over-estimate."""
+        return self.max_alt_id + 1
+
+
 # TODO-EET: add doc string, tracing
-def add_ev1_random(state: workflow.State, df: pd.DataFrame):
+def add_ev1_random(
+    state: workflow.State,
+    df: pd.DataFrame,
+    alt_info: AltsContext | None = None,
+    alt_nrs_df: pd.DataFrame | None = None,
+):
+
     nest_utils_for_choice = df.copy()
-    nest_utils_for_choice += state.get_rn_generator().gumbel_for_df(
-        nest_utils_for_choice, n=nest_utils_for_choice.shape[1]
-    )
+    assert (alt_info is None) == (
+        alt_nrs_df is None
+    ), "n_zones and alt_nrs_df must both be provided or omitted together"
+
+    if alt_nrs_df is not None and FREEZE_RANDOM_NUMBERS_FOR_DENSE_ALTERNATIVE_SET:
+        assert alt_info is not None  # narrowing for mypy
+
+        idx_array = alt_nrs_df.values
+        mask = idx_array == -999
+        safe_idx = np.where(
+            mask, 1, idx_array
+        )  # replace -999 with a temp value inbounds
+        # generate random number for all alts - this is wasteful, but ensures that the same zone
+        #  gets the same random number if the sampled choice set changes between base and project
+        # (alternatively, one could seed a channel for (persons x zones) and use the zone seed to ensure consistency.
+        # Trade off is needing to seed (persons x zones) rows and multiindex channels to
+        # avoid extra random numbers generated here. Quick benchmark suggests seeding per row is likely slower
+        rands_dense = state.get_rn_generator().gumbel_for_df(
+            nest_utils_for_choice, n=alt_info.n_alts_to_cover_max_id
+        )
+        # generate n=alt_info.max_alt_id+1 rather than n_alts so that indexing works
+        # (this is drawing a random number for a redundant zeroth zone in 1 based zoning systems)
+        # TODO deal with non 0->n-1 indexed land use more efficiently? ideally do where alt_nrs_df is constructed,
+        #  not on the fly here. Potentially via state.get_injectable('network_los').get_skim_dict('taz').zone_ids
+        rands = np.take_along_axis(rands_dense, safe_idx, axis=1)
+        rands[
+            mask
+        ] = 0  # zero out the masked zones so they don't have the util adjustment of alt 0
+    else:
+        # old behaviour, to remove
+        rands = state.get_rn_generator().gumbel_for_df(
+            nest_utils_for_choice, n=nest_utils_for_choice.shape[1]
+        )
+
+    nest_utils_for_choice += rands
     return nest_utils_for_choice
 
 
@@ -401,8 +476,14 @@ def make_choices_explicit_error_term_nl(
 
 
 # TODO-EET: add doc string, tracing
-def make_choices_explicit_error_term_mnl(state, utilities, trace_label):
-    utilities_incl_unobs = add_ev1_random(state, utilities)
+def make_choices_explicit_error_term_mnl(
+    state,
+    utilities,
+    trace_label,
+    alts_context: AltsContext | None = None,
+    alt_nrs_df: pd.DataFrame | None = None,
+):
+    utilities_incl_unobs = add_ev1_random(state, utilities, alts_context, alt_nrs_df)
     choices = np.argmax(utilities_incl_unobs.to_numpy(), axis=1)
     # TODO-EET: reporting like for zero probs
     assert not np.isnan(choices).any(), f"No choice for {trace_label}"
@@ -411,11 +492,19 @@ def make_choices_explicit_error_term_mnl(state, utilities, trace_label):
 
 
 def make_choices_explicit_error_term(
-    state, utilities, alt_order_array, nest_spec=None, trace_label=None
+    state,
+    utilities,
+    alt_order_array,
+    nest_spec=None,
+    trace_label=None,
+    alts_context: AltsContext | None = None,
+    alt_nrs_df: pd.DataFrame | None = None,
 ):
     trace_label = tracing.extend_trace_label(trace_label, "make_choices_eet")
     if nest_spec is None:
-        choices = make_choices_explicit_error_term_mnl(state, utilities, trace_label)
+        choices = make_choices_explicit_error_term_mnl(
+            state, utilities, trace_label, alts_context, alt_nrs_df
+        )
     else:
         choices = make_choices_explicit_error_term_nl(
             state, utilities, alt_order_array, nest_spec, trace_label
@@ -431,13 +520,15 @@ def make_choices_utility_based(
     trace_label: str = None,
     trace_choosers=None,
     allow_bad_probs=False,
+    alts_context: AltsContext | None = None,
+    alt_nrs_df: pd.DataFrame | None = None,
 ) -> tuple[pd.Series, pd.Series]:
     trace_label = tracing.extend_trace_label(trace_label, "make_choices_utility_based")
 
     # TODO-EET: index of choices for nested utilities is different than unnested - this needs to be consistent for
     #  turning indexes into alternative names to keep code changes to minimum for now
     choices = make_choices_explicit_error_term(
-        state, utilities, name_mapping, nest_spec, trace_label
+        state, utilities, name_mapping, nest_spec, trace_label, alts_context, alt_nrs_df
     )
     # TODO-EET: rands - log all zeros for now
     rands = pd.Series(np.zeros_like(utilities.index.values), index=utilities.index)
diff --git a/activitysim/core/simulate.py b/activitysim/core/simulate.py
index ed0b34452..5a23eb7c8 100644
--- a/activitysim/core/simulate.py
+++ b/activitysim/core/simulate.py
@@ -4,6 +4,7 @@
 
 import logging
 import time
+import typing
 import warnings
 from collections import OrderedDict
 from collections.abc import Callable
@@ -32,7 +33,8 @@
     LogitNestSpec,
     TemplatedLogitComponentSettings,
 )
-from activitysim.core.estimation import Estimator
+
+from activitysim.core.exceptions import ModelConfigurationError
 from activitysim.core.fast_eval import fast_eval
 from activitysim.core.simulate_consts import (
     ALT_LOSER_UTIL,
@@ -40,7 +42,10 @@
     SPEC_EXPRESSION_NAME,
     SPEC_LABEL_NAME,
 )
-from activitysim.core.exceptions import ModelConfigurationError
+
+if typing.TYPE_CHECKING:
+    from activitysim.core.estimation import Estimator
+
 
 logger = logging.getLogger(__name__)
 
diff --git a/activitysim/core/test/test_logit.py b/activitysim/core/test/test_logit.py
index c82606981..436e4cb17 100644
--- a/activitysim/core/test/test_logit.py
+++ b/activitysim/core/test/test_logit.py
@@ -3,13 +3,16 @@
 from __future__ import annotations
 
 import os.path
+import re
 
 import numpy as np
 import pandas as pd
 import pandas.testing as pdt
 import pytest
 
-from activitysim.core import logit, workflow
+
+from activitysim.core import logit, workflow, random
+from activitysim.core.logit import add_ev1_random, AltsContext
 from activitysim.core.simulate import eval_variables
 
 
@@ -131,6 +134,75 @@ def test_make_choices_only_one():
     )
 
 
+def reset_step(state, name="test_step"):
+    state.get_rn_generator().end_step(name)
+    state.get_rn_generator().begin_step(name)
+
+
+def test_make_choices_utility_based_sampled_alts():
+    """Test the situation of making choices from a sampled choice set"""
+    # TODO should these tests go in test_random?
+    state = workflow.State().default_settings()
+    # Make explicit that there's two indexing schemes - the raw alts, and the 0 based internals
+    utils_project_raw = pd.DataFrame(
+        {"a": 10.582999, "b": 10.680792, "c": 10.710443},
+        index=pd.Index([0], name="person_id"),
+    )
+    # zero based indexes
+    utils_project = utils_project_raw.rename(columns={"a": 0, "b": 1, "c": 2})
+    utils_base = utils_project_raw[["a", "c"]].rename(columns={"a": 0, "c": 1})
+
+    assert utils_project.index.name == "person_id"
+    state.get_rn_generator().add_channel("persons", utils_project)
+    state.get_rn_generator().begin_step("test_step")
+    # mock base case, where alt 1 is omitted (it was improved in the project)
+    # this situation is quite common with poisson sampling with a variable choice set size,
+    # but it can also happen in with-replacement EET sampling e.g. if alt 2 had a pick_count of 2 in the base case.
+    # In principle, it can also be problematic for non-sampled choices where there is a base project difference in the
+    # availability of alternatives .e.g a new mode was introduced in the project case
+
+    utils_project_with_rands = add_ev1_random(state, utils_project)
+    rands_project = utils_project_with_rands - utils_project
+    reset_step(state)
+    utils_base_with_rands = add_ev1_random(state, utils_base)
+    rands_base = utils_base_with_rands - utils_base
+    rands_base_labeled = rands_base.rename(columns={0: "a", 1: "c"})
+    rands_project_labeled = rands_project.rename(columns={0: "a", 1: "b", 2: "c"})
+    with pytest.raises(
+        AssertionError, match=re.escape('(column name="c") are different')
+    ):
+        # TODO this should pass
+        pdt.assert_frame_equal(
+            rands_base_labeled, rands_project_labeled.loc[:, rands_base_labeled.columns]
+        )
+    # document incorrect invariant - first two columns have the same random numbers:
+    pdt.assert_frame_equal(rands_base, rands_project.iloc[:, :2])
+
+    # revised approach
+    reset_step(state)
+    alt_nrs_df = pd.DataFrame({0: 0, 1: 1, 2: 2}, index=utils_project_raw.index)
+    alt_info = AltsContext.from_num_alts(3, zero_based=True)
+    utils_project_with_rands = add_ev1_random(
+        state, utils_project, alt_info=alt_info, alt_nrs_df=alt_nrs_df
+    )
+    rands_project = utils_project_with_rands - utils_project
+    reset_step(state)
+
+    # alt "b" is missing from the sampled choice set, alt_nrs_df is set to reflect that
+    alt_nrs_df = pd.DataFrame({0: 0, 1: 2}, index=utils_project_raw.index)
+    utils_base_with_rands = add_ev1_random(
+        state, utils_base, alt_info=alt_info, alt_nrs_df=alt_nrs_df
+    )
+    rands_base = utils_base_with_rands - utils_base
+    rands_base_labeled = rands_base.rename(columns={0: "a", 1: "c"})
+    rands_project_labeled = rands_project.rename(columns={0: "a", 1: "b", 2: "c"})
+
+    # Corrected invariant holds true
+    pdt.assert_frame_equal(
+        rands_base_labeled, rands_project_labeled.loc[:, rands_base_labeled.columns]
+    )
+
+
 def test_make_choices_real_probs(utilities):
     state = workflow.State().default_settings()
     probs = logit.utils_to_probs(state, utilities, trace_label=None)