Particle weights are now stored using the global normalization

CHANGELOG.rst

@@ -7,7 +7,7 @@ Release Notes
 0.12 Series
 ...........
 
-0.12.19 (2026-07-27)
+0.12.19 (2026-07-28)
 --------------------
 
 General:
@@ -24,12 +24,16 @@ Visualization:
   per-generation walltimes and no longer include the idle time that passed
   between a stored analysis and a later resumed run. Resolves #667.
 
-Storage:
+Storage (breaking):
 
 * Store the per-generation walltime in the database (new ``wall_time`` column,
   database version 2). Databases created with older pyABC versions must be
   migrated via ``abc-migrate`` before they can be resumed; for such databases
   the walltime plots fall back to the previous, end-time-based behavior.
+* Particle weights are now stored using the global normalization (weights sum
+  to 1 across all particles of all models, matching the in-memory
+  ``Population`` representation). This is bundled into database version 2 and
+  handled by ``abc-migrate`` for existing databases. Resolves #47.
 
 0.12.18 (2026-04-14)
 --------------------

README.md

@@ -26,4 +26,5 @@
 - 📄 **Cite**: [https://pyabc.rtfd.io/en/latest/cite.html](https://pyabc.rtfd.io/en/latest/cite.html)
 
 #### Related Projects
-- 🧠 **Neural Posterior Estimation**: [BayesFlow](https://bayesflow.org/main/_examples/From_ABC_to_BayesFlow.html)
+- **Parameter Estimation with Likelihoods**: [pyPESTO](https://github.com/ICB-DCM/pyPESTO)
+- **Neural Posterior Estimation**: [BayesFlow](https://bayesflow.org/main/_examples/From_ABC_to_BayesFlow.html)

pyabc/acceptor/pdf_norm.py

@@ -32,7 +32,10 @@ def pdf_norm_max_found(
         prev_pdf_norm = -np.inf
 
     # take maximum over all normalizations
-    pdf_norm = max(prev_pdf_norm, *pdfs)
+    if len(pdfs) == 0:
+        pdf_norm = prev_pdf_norm
+    else:
+        pdf_norm = max(prev_pdf_norm, float(pdfs.max()))
 
     return pdf_norm
 

pyabc/distance/distance.py

@@ -105,7 +105,9 @@ def _dict_to_vect(self, x):
 
     def _calculate_whitening_transformation_matrix(self, sum_stats):
         # create data matrix, shape (n_sample, n_y)
-        x = np.asarray([self._dict_to_vect(x) for x in sum_stats])
+        # force float dtype so in-place centering below works also for
+        # integer-valued summary statistics (e.g. counts)
+        x = np.asarray([self._dict_to_vect(x) for x in sum_stats], dtype=float)
         # center
         mean = np.mean(x, axis=0)
         x -= mean

pyabc/distance/pnorm.py

@@ -148,7 +148,7 @@ def get_weights(self, t: int) -> np.ndarray:
 
     @staticmethod
     def format_dict(
-        vals: dict[str, float] | dict[int, dict[str, float]],
+        vals: dict[str, float] | dict[int, dict[str, float]] | None,
         t: int,
         s_ids: list[str],
     ) -> dict[int, float | np.ndarray]:
@@ -172,8 +172,7 @@ def format_dict(
             vals = {t: vals}
 
         # convert dicts to arrays
-        for _t, dct in vals.items():
-            vals[_t] = dict2arr(dct, keys=s_ids)
+        vals = {_t: dict2arr(dct, keys=s_ids) for _t, dct in vals.items()}
 
         return vals
 
@@ -322,7 +321,9 @@ def __init__(
         # call p-norm constructor
         super().__init__(p=p, fixed_weights=fixed_weights, sumstat=sumstat)
 
-        self.initial_scale_weights: dict[str, float] = initial_scale_weights
+        self.initial_scale_weights: dict[str, float] | None = (
+            initial_scale_weights
+        )
 
         self.scale_weights: dict[int, np.ndarray] = {}
 

pyabc/distance/scale.py

@@ -141,7 +141,7 @@ def root_mean_square_deviation(
     rmse = np.sqrt(mse)
 
     # debugging
-    warn_obs_off(off_ixs=np.flatnonzero(bs > 2 * std), s_ids=s_ids)
+    warn_obs_off(off_ixs=np.flatnonzero(np.abs(bs) > 2 * std), s_ids=s_ids)
 
     return rmse
 
@@ -161,15 +161,15 @@ def std_or_rmsd(
     bs = bias(samples=samples, s0=s0)
     std = standard_deviation(samples=samples)
 
-    if sum(bs > 2 * std) > 1 / 3 * len(std):
+    if sum(np.abs(bs) > 2 * std) > 1 / 3 * len(std):
         logger.info('Too many high-bias values, correcting only for scale.')
         return std
 
     mse = bs**2 + std**2
     rmse = np.sqrt(mse)
 
     # debugging
-    warn_obs_off(off_ixs=np.flatnonzero(bs > 2 * std), s_ids=s_ids)
+    warn_obs_off(off_ixs=np.flatnonzero(np.abs(bs) > 2 * std), s_ids=s_ids)
 
     return rmse
 

pyabc/epsilon/temperature.py

@@ -86,6 +86,9 @@ def __init__(
         enforce_exact_final_temperature: bool = True,
         log_file: str | None = None,
     ):
+        # normalize a single callable to a list, as all consumers iterate
+        if schemes is not None and callable(schemes):
+            schemes = [schemes]
         self.schemes = schemes
 
         if aggregate_fun is None:

pyabc/inference_util/inference_util.py

@@ -690,10 +690,10 @@ def create_analysis_id():
     return str(uuid.uuid4())
 
 
-def eps_from_hist(history: History, t: int | None = None) -> float:
+def eps_from_hist(history: History, t: int | None = None) -> float | None:
     """Read epsilon value for time `t` from `history`. Defaults to latest."""
     pops = history.get_all_populations()
-    if len(pops) == 0 or (t is not None and t not in pops.t):
+    if len(pops) == 0 or (t is not None and t not in pops.t.values):
         return None
     if t is None:
         return pops.epsilon.to_numpy()[-1]

pyabc/predictor/predictor.py

@@ -732,6 +732,8 @@ def __init__(
 
         if f_score is None:
             self.f_score = root_mean_square_error
+        else:
+            self.f_score = f_score
 
         # holds the chosen predictor model
         self.chosen_one: Predictor | None = None

pyabc/sampler/multicore.py

@@ -36,7 +36,10 @@ def work(feed_q, result_q, simulate_one, max_eval, single_core_sampler):
             break
 
         res = single_core_sampler.sample_until_n_accepted(
-            1, simulate_one, max_eval
+            1,
+            simulate_one,
+            t=0,  # t is not used in this sampler
+            max_eval=max_eval,
         )
         result_q.put((res, single_core_sampler.nr_evaluations_))
 

pyabc/sampler/redis_eps/cli.py

@@ -199,7 +199,8 @@ def _work(
             logger.info('Received stop signal. Shutdown redis worker.')
             return
 
-        # TODO other messages (some integers?) are ignored
+        # Any other messages on the channel (e.g. redis subscription-count
+        # notifications) are ignored.
 
         # check total time condition
         elapsed_time = time() - start_time

pyabc/sge/db.py

@@ -49,11 +49,9 @@ def wait_for_job(self, ID, max_run_time_h):
         Return true if we should still wait for the job.
         Return false otherwise
         """
-        # TODO Possible SQL injection error should be fixed, e.g. via
-        #  pre-calculated expressions
         with self.connection:
             results = self.connection.execute(
-                'SELECT status, time from status WHERE ID=' + str(ID)
+                'SELECT status, time from status WHERE ID=?', (ID,)
             ).fetchall()
             nr_rows = len(results)
 

pyabc/storage/bytes_storage.py

@@ -16,7 +16,7 @@ def r_to_py(object_):
                 py_object_ = conv.rpy2py(object_)
             # Ensure factor columns are converted to strings
             for col in py_object_.columns:
-                if isinstance(py_object_[col], pd.CategoricalDtype):
+                if isinstance(py_object_[col].dtype, pd.CategoricalDtype):
                     py_object_[col] = py_object_[col].astype(str)
             return py_object_
         return object_

pyabc/storage/db_model.py

@@ -169,6 +169,8 @@ class SummaryStatistic(Base):
     value = Column(BytesStorage)
 
 
-def datetime2str(datetime: datetime.datetime) -> str:
+def datetime2str(datetime: datetime.datetime | None) -> str:
     """Format print datetime."""
+    if datetime is None:
+        return 'None'
     return datetime.strftime('%Y-%m-%d %H:%M:%S')

pyabc/storage/history.py

@@ -344,10 +344,10 @@ def get_distribution(
         ).sort_index()
         w = df[['id', 'w']].drop_duplicates().set_index('id').sort_index()
         w_arr = w.w.values
-        if w_arr.size > 0 and not np.isclose(w_arr.sum(), 1):
-            raise AssertionError(
-                f'Weight not close to 1, w.sum()={w_arr.sum()}'
-            )
+        # Stored weights are global (summing to 1 across all models); within a
+        # single model they sum to the model probability.
+        if w_arr.size > 0:
+            w_arr = w_arr / w_arr.sum()
         return pars, w_arr
 
     @with_session
@@ -735,18 +735,16 @@ def _save_to_population_db(
             # append model
             population.models.append(model)
 
-            # TODO This normalization is different than in the in-memory
-            #  population. It would be cleaner to update the db too.
-            total_model_weight = sum(p.weight for p in model_population)
-
             # iterate over model population of particles
             for py_particle in model_population:
                 # a store_item is a Particle
                 py_parameter = py_particle.parameter
 
-                # create new particle
+                # create new particle. The stored weight is the global
+                # particle weight (normalized across all particles of all
+                # models so that they sum to 1)
                 particle = Particle(
-                    w=py_particle.weight / total_model_weight,
+                    w=py_particle.weight,
                     proposal_id=py_particle.proposal_id,
                 )
                 # append particle to model
@@ -858,14 +856,13 @@ def get_model_probabilities(self, t: int | None = None) -> pd.DataFrame:
             .order_by(Model.m)
             .all()
         )
-        # TODO this is a mess
+        # Two return shapes: for a single t, a frame indexed by model id `m`
+        # with a `p` column; for t=None, a t-by-model pivot table (below).
         if t is not None:
             p_models_df = pd.DataFrame(
                 [p[:2] for p in p_models], columns=['p', 'm']
             ).set_index('m')
-            # TODO the following line is redundant
-            # only models with no-zero weight are stored for each population
-            p_models_df = p_models_df[p_models_df.p >= 0]
+            # Only models with non-zero weight are stored per population
             return p_models_df
         else:
             p_models_df = (
@@ -933,7 +930,7 @@ def get_weighted_distances(self, t: int | None = None) -> pd.DataFrame:
         distances = []
         for model in models:
             for particle in model.particles:
-                weight = particle.w * model.p_model
+                weight = particle.w
                 for sample in particle.samples:
                     weights.append(weight)
                     distances.append(sample.distance)
@@ -1018,14 +1015,31 @@ def get_weighted_sum_stats_for_model(
             .filter(ABCSMC.id == self.id)
             .filter(Population.t == t)
             .filter(Model.m == m)
+            .options(
+                subqueryload(Particle.samples).subqueryload(
+                    Sample.summary_statistics
+                )
+            )
             .all()
         )
 
+        # model probability, used to renormalize the stored global weights
+        # back to the within-model posterior (weights summing to 1)
+        p_model = (
+            self._session.query(Model.p_model)
+            .join(Population)
+            .join(ABCSMC)
+            .filter(ABCSMC.id == self.id)
+            .filter(Population.t == t)
+            .filter(Model.m == m)
+            .scalar()
+        )
+
         results = []
         weights = []
         for particle in particles:
             for sample in particle.samples:
-                weights.append(particle.w)
+                weights.append(particle.w / p_model if p_model else particle.w)
                 sum_stats = {}
                 for ss in sample.summary_statistics:
                     sum_stats[ss.name] = ss.value
@@ -1075,7 +1089,7 @@ def get_weighted_sum_stats(
 
         for model in models:
             for particle in model.particles:
-                weight = particle.w * model.p_model
+                weight = particle.w
                 for sample in particle.samples:
                     # extract sum stats
                     sum_stats = {}
@@ -1125,7 +1139,7 @@ def get_population(self, t: int | None = None):
             py_m = model.m
             for particle in model.particles:
                 # weight
-                py_weight = particle.w * model.p_model
+                py_weight = particle.w
 
                 # parameter
                 py_parameter = {}
@@ -1134,7 +1148,9 @@ def get_population(self, t: int | None = None):
                 py_parameter = PyParameter(**py_parameter)
 
                 # simulations
-                # TODO this is legacy from when there were multiple
+                # NOTE: samples is a one-to-many relationship for legacy
+                # reasons (a particle could once store multiple samples);
+                # today exactly one is expected.
                 if len(particle.samples) != 1:
                     raise AssertionError('There should be exactly one sample.')
                 sample = particle.samples[0]

pyabc/storage/json.py

@@ -19,7 +19,7 @@ def save_dict_to_json(dct: dict, file_: str):
     for key, val in dct.items():
         # cannot handle ndarrays
         if isinstance(val, np.ndarray):
-            dct[key] = list(val)
+            dct[key] = val.tolist()
     with open(file_, 'w') as f:
         json.dump(dct, f)