Stabilization: modernize deps and fix the instabilities that surfaced

.github/workflows/python-package.yml

@@ -11,12 +11,13 @@ on:
 
 jobs:
   build:
-
-    runs-on: ubuntu-latest
+    name: Test (${{ matrix.os }}, Python ${{ matrix.python-version }})
+    runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        os: [ubuntu-latest, windows-latest, macos-latest]
+        python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"]
 
     steps:
     - uses: actions/checkout@v4
@@ -25,7 +26,7 @@ jobs:
       uses: astral-sh/setup-uv@v5
       with:
         python-version: ${{ matrix.python-version }}
-        version: "0.6.14"
+        version: "0.11.8"
 
     - name: Install the project
       run: uv sync --all-extras --dev

.gitignore

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

README.md

@@ -1,7 +1,7 @@
 PopulationSim
 =============
 
-[![Build Status](https://travis-ci.org/activitysim/populationsim.svg?branch=master)](https://travis-ci.org/ActivitySim/populationsim) [![Coverage Status](https://coveralls.io/repos/ActivitySim/populationsim/badge.png?branch=master)](https://coveralls.io/r/ActivitySim/populationsim?branch=master)<a href="https://medium.com/zephyrfoundation/populationsim-the-synthetic-commons-670e17383048"><img src="https://github.com/ZephyrTransport/zephyr-website/blob/gh-pages/img/badging/project_pages/populationsim/PopulationSim.png" width="72.6" height="19.8"></a>
+[![Python package](https://github.com/ActivitySim/populationsim/actions/workflows/python-package.yml/badge.svg?branch=master)](https://github.com/ActivitySim/populationsim/actions/workflows/python-package.yml)<a href="https://medium.com/zephyrfoundation/populationsim-the-synthetic-commons-670e17383048"><img src="https://github.com/ZephyrTransport/zephyr-website/blob/gh-pages/img/badging/project_pages/populationsim/PopulationSim.png" width="72.6" height="19.8"></a>
 
 
 PopulationSim is an open platform for population synthesis.  It emerged

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

docs/software.rst

@@ -57,7 +57,7 @@ being implemented in the ActivitySim framework means:
 
   * A protected master branch that can only be written to after tests have passed
   * `pytest <https://docs.pytest.org/en/latest/>`__ for tests
-  * `TravisCI <https://travis-ci.org>`__ for building and testing with each commit
+  * `GitHub Actions <https://github.com/ActivitySim/populationsim/actions>`__ for building and testing with each commit
 
 PopulationSim also requires an optimization library for balancing and integerizing.  The software makes
 use of the open source and easy to install `ortools <https://github.com/google/or-tools>`__ package.  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/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.15"
 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 }