Stabilization: modernize deps and fix the instabilities that surfaced

.github/workflows/python-package.yml

@@ -11,12 +11,31 @@ on:
 
 jobs:
   build:
-
-    runs-on: ubuntu-latest
+    name: Test (${{ matrix.os }}, Python ${{ matrix.python-version }})
+    runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
+      # Deliberately not a square matrix, so the combinations are listed out
+      # rather than generated. Linux carries the full Python sweep; Windows and
+      # macOS get a single spot-check each. Those two exist to catch platform
+      # divergence rather than version divergence -- word size (`int` is int32
+      # on Windows) and differences in the compiled solver -- which one Python
+      # version exercises just as well as four. A square os x python matrix
+      # would double the job count to buy very little.
       matrix:
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        include:
+        - os: ubuntu-latest
+          python-version: "3.10"
+        - os: ubuntu-latest
+          python-version: "3.11"
+        - os: ubuntu-latest
+          python-version: "3.12"
+        - os: ubuntu-latest
+          python-version: "3.13"
+        - os: windows-latest
+          python-version: "3.12"
+        - os: macos-latest
+          python-version: "3.12"
 
     steps:
     - uses: actions/checkout@v4
@@ -35,3 +54,12 @@ jobs:
 
     - name: Run tests
       run: uv run pytest
+
+    - name: Print determinism fingerprint (cross-platform check)
+      # Unperturbed run of the quantized integerizer. Compare the printed
+      # hash across this matrix's jobs by hand -- an equal hash on every OS
+      # confirms quantization also makes CBC's own platform-dependent tie
+      # choice converge, not just our own pre-solver rounding.
+      env:
+        PYTHONPATH: ${{ github.workspace }}
+      run: uv run python tests/determinism_probe.py 0 1

.gitignore

@@ -1,5 +1,7 @@
 scripts/calm_validation_results/
 
+.venv/
+
 regress/
 .idea
 .ipynb_checkpoints

docs/conf.py

@@ -16,12 +16,14 @@
 import sphinx_rtd_theme
 
 # -- Get Package Version --------------------------------------------------
-with open("../setup.py") as file:
-    lines = file.readlines()
-    for line in lines:
-        if "version" in line:
-            VERSION = line.replace("version='", "").replace("',", "").replace(" ", "")
-            print("package version: " + VERSION)
+# the version lives in pyproject.toml now; setup.py was removed
+try:
+    from importlib.metadata import version as get_version
+
+    VERSION = get_version("populationsim")
+except Exception:
+    VERSION = "unknown"
+print("package version: " + VERSION)
 
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the

populationsim/balancing/simul_balancer.py

@@ -174,7 +174,7 @@ def balance(self):
             # max_delta=
         )
 
-        status = dict(zip(("converged", "iter", "delta", "max_gamma_dif"), status))
+        status = dict(zip(("converged", "iter", "delta", "max_gamma_dif"), status, strict=True))
 
         # dataframe with sub_zone_weights in columns, and zero weight rows restored
         self.sub_zone_weights = pd.DataFrame(

populationsim/balancing/single_balancer.py

@@ -161,7 +161,7 @@ def balance(self):
         )
 
         # Label the status
-        status = dict(zip(("converged", "iter", "delta", "max_gamma_dif"), status))
+        status = dict(zip(("converged", "iter", "delta", "max_gamma_dif"), status, strict=True))
 
         # weights dataframe
         weights = pd.DataFrame(index=self.incidence_table.index)

populationsim/core/input.py

@@ -153,7 +153,7 @@ def read_from_table_info(table_info):
                     map_col = parent_table[f"_original_{lookup_col}"]
                 except KeyError:
                     map_col = parent_table[lookup_col]
-                remapper = dict(zip(map_col, parent_table.index))
+                remapper = dict(zip(map_col, parent_table.index, strict=True))
                 df[colname] = df[colname].apply(remapper.get)
 
     # set index
@@ -239,6 +239,15 @@ def _read_csv_with_fallback_encoding(filepath, dtypes=None):
     if dtypes:
         # although the dtype argument suppresses the DtypeWarning, it does not coerce recognized types (e.g. int)
         for c, dtype in dtypes.items():
-            df[c] = df[c].astype(dtype)
+            missing = df[c].isna()
+            converted = df[c].astype(dtype)
+            # pandas 2 casts missing values to the literal string "nan" under
+            # astype(str), while pandas 3's string dtype preserves NA. Control
+            # expressions such as `persons.PComm.isna()` then silently evaluate
+            # to all-False on pandas 2, dropping a whole control column, so
+            # restore NA explicitly for string casts.
+            if dtype in ("str", "string", str) and missing.any():
+                converted = converted.where(~missing)
+            df[c] = converted
 
     return df

populationsim/core/mp_tasks.py

@@ -1030,7 +1030,7 @@ def run_sub_simulations(
     """
 
     def log_queued_messages():
-        for process, queue in zip(procs, queues):
+        for process, queue in zip(procs, queues, strict=True):
             while not queue.empty():
                 msg = queue.get(block=False)
                 model_name = msg["model"]
@@ -1140,7 +1140,7 @@ def check_proc_status():
         queues.append(q)
 
     # - start processes
-    for _, p in zip(list(range(num_simulations)), procs):
+    for _, p in zip(list(range(num_simulations)), procs, strict=True):
         info(f"start process {p.name}")
         p.start()
 

populationsim/integerizing/reproducibility.py

@@ -0,0 +1,33 @@
+import numpy as np
+
+
+def quantize_weights(weights, quantum):
+    """Snap weights to an absolute grid when quantization is enabled."""
+    values = np.asarray(weights, dtype=np.float64)
+
+    if quantum in (None, False, 0, 0.0):
+        return values
+    if isinstance(quantum, bool):
+        raise ValueError("INTEGERIZER_QUANTUM must be a positive finite number")
+
+    try:
+        quantum = float(quantum)
+    except (TypeError, ValueError) as err:
+        raise ValueError(
+            "INTEGERIZER_QUANTUM must be a positive finite number"
+        ) from err
+    if not np.isfinite(quantum) or quantum <= 0:
+        raise ValueError("INTEGERIZER_QUANTUM must be a positive finite number")
+
+    quantized = np.rint(values / quantum) * quantum
+    # Zero has structural meaning in the integerizers: it makes a household
+    # ineligible for selection. Quantization must not turn an eligible,
+    # positive household into an ineligible one. Tiny positive values share
+    # one canonical sentinel; the existing log overflow guard maps it to the
+    # same finite objective coefficient in every case.
+    positive_to_zero = (values > 0) & (quantized == 0)
+    quantized[positive_to_zero] = np.nextafter(0.0, 1.0)
+    if not np.isfinite(quantized).all():
+        raise ValueError("INTEGERIZER_QUANTUM produced non-finite weights")
+
+    return quantized

populationsim/integerizing/simul_integerizer.py

@@ -8,6 +8,7 @@
 
 from populationsim.core import config
 from populationsim.integerizing.smart_round import smart_round
+from populationsim.integerizing.reproducibility import quantize_weights
 from populationsim.integerizing import lp_ortools, lp_cvx
 
 logger = logging.getLogger(__name__)
@@ -64,7 +65,11 @@ def __init__(
         else:
             self.integerizer_func = lp_ortools.np_simul_integerizer_ortools
 
-        self.timeout_in_seconds = config.setting("INTEGIZER_TIMEOUT", 60)
+        self.timeout_in_seconds = config.setting("INTEGERIZER_TIMEOUT", 60)
+        # Snapping weights to a grid before rounding makes household
+        # selection immune to sub-ULP cross-platform float noise, so this
+        # defaults to on -- see tests/test_reproducibility.py.
+        self.quantum = config.setting("INTEGERIZER_QUANTUM", 1e-6)
 
     def integerize(self):
 
@@ -80,7 +85,9 @@ def integerize(self):
         sub_incidence = self.incidence_df[self.sub_controls_df.columns]
         sub_incidence = sub_incidence.values.astype(np.float64)
 
-        sub_float_weights = self.sub_weights.values.transpose().astype(np.float64)
+        sub_float_weights = quantize_weights(
+            self.sub_weights.values.transpose(), self.quantum
+        )
         sub_int_weights = sub_float_weights.astype(int)
         sub_resid_weights = sub_float_weights % 1.0
 
@@ -202,7 +209,10 @@ def integerize(self):
         sub_zone_count = len(self.sub_weights.columns)
         for i in range(sub_zone_count):
             integerized_weights[i] = smart_round(
-                sub_int_weights[i], resid_weights_out[i], total_household_controls[i]
+                sub_int_weights[i],
+                resid_weights_out[i],
+                total_household_controls[i],
+                tie_break_by_position=bool(self.quantum),
             )
 
         # integerized_weights df: one column of integerized weights per sub_zone

populationsim/integerizing/single_integerizer.py

@@ -9,6 +9,7 @@
 from populationsim.core import config
 from populationsim.integerizing.constants import STATUS_OPTIMAL
 from populationsim.integerizing.smart_round import smart_round
+from populationsim.integerizing.reproducibility import quantize_weights
 from populationsim.integerizing import lp_cvx, lp_ortools
 
 logger = logging.getLogger(__name__)
@@ -60,15 +61,19 @@ def __init__(
         else:
             self.integerizer_func = lp_ortools.np_integerizer_ortools
 
-        self.timeout_in_seconds = config.setting("INTEGIZER_TIMEOUT", 60)
+        self.timeout_in_seconds = config.setting("INTEGERIZER_TIMEOUT", 60)
+        # Snapping weights to a grid before rounding makes household
+        # selection immune to sub-ULP cross-platform float noise, so this
+        # defaults to on -- see tests/test_reproducibility.py.
+        self.quantum = config.setting("INTEGERIZER_QUANTUM", 1e-6)
 
     def integerize(self):
 
         sample_count = len(self.incidence_table.index)
         control_count = len(self.incidence_table.columns)
 
         incidence = self.incidence_table.values.transpose().astype(np.float64)
-        float_weights = np.asanyarray(self.float_weights).astype(np.float64)
+        float_weights = quantize_weights(self.float_weights, self.quantum)
         relaxed_control_totals = np.asanyarray(self.relaxed_control_totals).astype(
             np.float64
         )
@@ -165,7 +170,10 @@ def integerize(self):
             )
 
             integerized_weights = smart_round(
-                int_weights, resid_weights, self.total_hh_control_value
+                int_weights,
+                resid_weights,
+                self.total_hh_control_value,
+                tie_break_by_position=bool(self.quantum),
             )
 
         self.weights = pd.DataFrame(index=self.incidence_table.index)

populationsim/integerizing/smart_round.py

@@ -1,7 +1,9 @@
 import numpy as np
 
 
-def smart_round(int_weights, resid_weights, target_sum):
+def smart_round(
+    int_weights, resid_weights, target_sum, tie_break_by_position=False
+):
     """
     Round weights while ensuring (as far as possible that result sums to target_sum)
 
@@ -32,8 +34,14 @@ def smart_round(int_weights, resid_weights, target_sum):
 
     # Order the residual weights and round at the tipping point where target_sum is achieved
     if int_shortfall > 0:
-        # indices of the int_shortfall highest resid_weights
-        i = np.argsort(resid_weights)[-int_shortfall:]
+        if tie_break_by_position:
+            # Sort by descending residual, breaking exact ties by original
+            # position. Unlike argsort, the tie behavior is explicit and stable.
+            positions = np.arange(len(resid_weights))
+            i = np.lexsort((positions, -resid_weights))[:int_shortfall]
+        else:
+            # Preserve historical output when reproducibility mode is disabled.
+            i = np.argsort(resid_weights)[-int_shortfall:]
 
         # add 1 to the integer weights that we want to round upwards
         rounded_weights[i] += 1

populationsim/integerizing/wrappers.py

@@ -129,7 +129,7 @@ def reshape_result(
         zone_weights_df[sub_geography] = zone_id
         zone_weights_df["balanced_weight"] = float_weights[zone_name].values
         zone_weights_df["integer_weight"] = (
-            integerized_weights[zone_name].astype(int).values
+            integerized_weights[zone_name].astype(np.int64).values
         )
 
         integer_weights_list.append(zone_weights_df)
@@ -538,7 +538,7 @@ def do_sequential_integerizing(
         zone_weights_df[weights.index.name] = weights.index
         zone_weights_df[sub_geography] = zone_id
         zone_weights_df["balanced_weight"] = weights.values
-        zone_weights_df["integer_weight"] = integer_weights.astype(int).values
+        zone_weights_df["integer_weight"] = integer_weights.astype(np.int64).values
 
         if status in STATUS_SUCCESS:
             integerized_weights_list.append(zone_weights_df)

populationsim/steps/expand_households.py

@@ -70,7 +70,7 @@ def expand_households():
         grouper = household_groups.groupby("group_id")
         group_hh_probs = [0] * len(grouper)
         for group_id, df in grouper:
-            hh_ids = list(df[household_id_col])
+            hh_ids = df[household_id_col].values  # preserves int64 dtype as numpy array
             probs = list(df.sample_weight / df.sample_weight.sum())
             group_hh_probs[group_id] = [hh_ids, probs]
 

populationsim/steps/setup_data_structures.py

@@ -182,7 +182,7 @@ def build_control_table(geo, control_spec, crosswalk_df):
     controls = pd.concat(controls_list, axis=1)
 
     # rename columns from seed_col to target
-    columns = {c: t for c, t in zip(control_spec.control_field, control_spec.target)}
+    columns = {c: t for c, t in zip(control_spec.control_field, control_spec.target, strict=True)}
     controls.rename(columns=columns, inplace=True)
 
     # reorder columns to match order of control_spec rows

pyproject.toml

@@ -6,14 +6,14 @@ authors = [
     { name = "Ben Stabler", email = "ben.stabler@rsginc.com" }
     ]
 readme = "README.md"
-requires-python = ">=3.9,<3.13"
+requires-python = ">=3.10,<3.14"
 dependencies = [
-    "ortools>=5.1.4045",
+    "numpy>=2.0.0,<3",
+    "ortools>=9.14.0",
     "pyinstrument>=5.0.1",
     "pyyaml>=6.0.2",
     "psutil>=7.0.0",
     "pandas>=2.2",
-    "numpy>=1.16.1,<2",
     "tables>=3.9",
     "orca>=1.8",
     "blosc2>=2.5.1",
@@ -44,6 +44,7 @@ testpaths = [
 
 [tool.setuptools.packages.find]
 where = ["."]
+include = ["populationsim*"]
 
 [tool.uv.sources]
 populationsim = { workspace = true }