fix: parallel fits ignoring updated priors, and multi-start comparison

CHANGELOG.md

@@ -20,6 +20,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `cpm.optimisation.FminBound` with `approx_grad=True` now computes the finite-difference gradient itself, together with the objective, exactly as SciPy does, with identical results; with `LogLikelihood.bernoulli`, this makes the hierarchical fits of the speed comparison 19% faster
 - The defaults of the new arguments, `hessian="numdifftools"` and `start="random"`, give the same results as before for the same seed; the faster options have to be chosen
 
+### Fixed
+
+- Fixed fits with `parallel=True` and `prior=True` going back to the priors the model was created with after the first evaluation in every worker, so that priors updated with `update_prior`, such as the population priors of `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`, were ignored; parallel fits now equal serial ones, and parallel results with updated priors change. `cpm.generators` no longer deletes its submodules, which made dill copy its classes into the workers instead of referencing them
+- Fixed `Fmin`, `FminBound`, `Minimize` and `Bads` comparing every later start with the first instead of the best so far, which with three or more starts could replace a participant's best fit with a worse one; results with three or more starts can change
+- Fixed the `cpm.hierarchical.VariationalBayes` docstring giving 1000 as the default of `iteration`, which is 50
+
 ## [0.26.0.dev2] - 2026-10-06
 
 ### Added

cpm/generators/__init__.py

@@ -8,5 +8,3 @@
 from .wrapper import Wrapper
 from .session import SessionWrapper
 from .parameters import Parameters, Value, LogParameters
-
-del simulator, wrapper, session, parameters

cpm/hierarchical/variational.py

@@ -22,7 +22,7 @@ class VariationalBayes:
     objective : str
         The objective of the optimisation, either 'maximise' or 'minimise'. Default is 'minimise'. Only affects how we arrive at the participant-level *a posteriori* parameter estimates.
     iteration : int, optional
-        The maximum number of iterations. Default is 1000.
+        The maximum number of iterations. Default is 50.
     tolerance_lme : float, optional
         The tolerance for convergence with respect to the log model evidence. Default is 1e-3.
     tolerance_param : float, optional

cpm/optimisation/bads.py

@@ -348,8 +348,9 @@ def __record(results):
                 __record(results)
             else:
                 nll = __extract_nll(results)
-                # check if ppt fit is better than the previous fit
+                # check if ppt fit is better than the best fit so far
                 indices = np.where(nll < old_nll)[0]
+                old_nll[indices] = nll[indices]
                 for ppt in indices:
                     self.details[ppt] = copy.deepcopy(results[ppt])
                     for i in range(len(self.parameter_names)):
@@ -362,12 +363,14 @@ def __record(results):
 
         if self.hessian == "finite_differences":
             ## choose each participant's start as above: a later start replaces
-            ## the first where it is better
+            ## the best so far where it is better
             best = list(starts[0])
-            first_nll = __extract_nll(starts[0])
+            best_nll = __extract_nll(starts[0])
             for results in starts[1:]:
-                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                nll = __extract_nll(results)
+                for ppt in np.where(nll < best_nll)[0]:
                     best[ppt] = results[ppt]
+                    best_nll[ppt] = nll[ppt]
             items = list(zip(self.data, best))
             if self.__parallel__:
                 results = execute_parallel(

cpm/optimisation/fmin.py

@@ -379,8 +379,9 @@ def __record(results):
                 __record(results)
             else:
                 nll = __extract_nll(results)
-                # check if ppt fit is better than the previous fit
+                # check if ppt fit is better than the best fit so far
                 indices = np.where(nll < old_nll)[0]
+                old_nll[indices] = nll[indices]
                 for ppt in indices:
                     self.details[ppt] = copy.deepcopy(results[ppt])
                     for i in range(len(self.parameter_names)):
@@ -393,12 +394,14 @@ def __record(results):
 
         if self.hessian == "finite_differences":
             ## choose each participant's start as above: a later start replaces
-            ## the first where it is better
+            ## the best so far where it is better
             best = list(starts[0])
-            first_nll = __extract_nll(starts[0])
+            best_nll = __extract_nll(starts[0])
             for results in starts[1:]:
-                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                nll = __extract_nll(results)
+                for ppt in np.where(nll < best_nll)[0]:
                     best[ppt] = results[ppt]
+                    best_nll[ppt] = nll[ppt]
             items = list(zip(self.data, best))
             if self.__parallel__:
                 results = execute_parallel(
@@ -739,8 +742,9 @@ def __record(results):
                 __record(results)
             else:
                 nll = __extract_nll(results)
-                # check if ppt fit is better than the previous fit
+                # check if ppt fit is better than the best fit so far
                 indices = np.where(nll < old_nll)[0]
+                old_nll[indices] = nll[indices]
                 for ppt in indices:
                     self.details[ppt] = copy.deepcopy(results[ppt])
                     for i in range(len(self.parameter_names)):
@@ -753,12 +757,14 @@ def __record(results):
 
         if self.hessian == "finite_differences":
             ## choose each participant's start as above: a later start replaces
-            ## the first where it is better
+            ## the best so far where it is better
             best = list(starts[0])
-            first_nll = __extract_nll(starts[0])
+            best_nll = __extract_nll(starts[0])
             for results in starts[1:]:
-                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                nll = __extract_nll(results)
+                for ppt in np.where(nll < best_nll)[0]:
                     best[ppt] = results[ppt]
+                    best_nll[ppt] = nll[ppt]
             items = list(zip(self.data, best))
             if self.__parallel__:
                 results = execute_parallel(

cpm/optimisation/free.py

@@ -282,8 +282,9 @@ def __record(results):
                 __record(results)
             else:
                 nll = __extract_nll(results)
-                # check if ppt fit is better than the previous fit
+                # check if ppt fit is better than the best fit so far
                 indices = np.where(nll < old_nll)[0]
+                old_nll[indices] = nll[indices]
                 for ppt in indices:
                     self.details[ppt] = copy.deepcopy(results[ppt])
                     for i in range(len(self.parameter_names)):
@@ -296,12 +297,14 @@ def __record(results):
 
         if self.hessian == "finite_differences":
             ## choose each participant's start as above: a later start replaces
-            ## the first where it is better
+            ## the best so far where it is better
             best = list(starts[0])
-            first_nll = __extract_nll(starts[0])
+            best_nll = __extract_nll(starts[0])
             for results in starts[1:]:
-                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                nll = __extract_nll(results)
+                for ppt in np.where(nll < best_nll)[0]:
                     best[ppt] = results[ppt]
+                    best_nll[ppt] = nll[ppt]
             items = list(zip(self.data, best))
             if self.__parallel__:
                 results = execute_parallel(

test/optimisation/test_parallel.py

@@ -0,0 +1,55 @@
+"""
+`parallel=True` sends the model to the workers with dill. dill copied the classes
+of `cpm.generators` by value instead of by reference, because the package
+deleted its submodules, so `isinstance` checks failed in the workers: after the
+first run, `Wrapper.reset` replaced the priors updated with `update_prior` with
+the priors the model was created with, and fits with `prior=True` differed from
+serial ones.
+"""
+
+import dill
+import numpy as np
+import pytest
+
+from cpm.generators import LogParameters, Parameters, SessionWrapper, Simulator, Value, Wrapper
+from cpm.optimisation import FminBound
+
+from test.optimisation.test_fminbound_gradient import LOSS, bandit, rlrw
+
+
+@pytest.mark.parametrize("cls", [Value, Parameters, LogParameters, Wrapper, SessionWrapper, Simulator])
+def test_dill_pickles_the_classes_by_reference(cls):
+    assert dill.loads(dill.dumps(cls)) is cls
+
+
+def test_a_copy_keeps_updated_priors_across_runs():
+    data = bandit(participants=1)
+    model = dill.loads(dill.dumps(rlrw(data)))
+    for x in ([0.3, 3.0], [0.05, 1.0]):
+        model.reset(parameters=np.array(x))
+        model.run()
+        assert model.parameters.alpha.prior.kwds["scale"] == 0.5
+
+
+def fit(parallel):
+    data = bandit(participants=4, trials=60, seed=3)
+    optimiser = FminBound(
+        model=rlrw(data),
+        data=data.groupby("ppt"),
+        initial_guess=[[0.3, 3.0]],
+        minimisation=LOSS,
+        prior=True,
+        approx_grad=True,
+        parallel=parallel,
+        cl=2,
+        ppt_identifier="ppt",
+    )
+    optimiser.optimise()
+    return optimiser
+
+
+def test_parallel_fits_with_updated_priors_equal_serial_ones():
+    serial, parallel = fit(False), fit(True)
+    for a, b in zip(serial.fit, parallel.fit):
+        np.testing.assert_array_equal(a["x"], b["x"])
+        assert a["fun"] == b["fun"]

test/optimisation/test_starts.py

@@ -0,0 +1,75 @@
+"""
+With several starts, every optimiser keeps each participant's best start. Later
+starts used to be compared only with the first, so with three or more starts, a
+start that was worse than the second but better than the first replaced the
+second.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from cpm.generators import Parameters, Value, Wrapper
+from cpm.optimisation import Fmin, FminBound, Minimize
+
+try:
+    from cpm.optimisation import Bads
+
+    HAS_BADS = True
+except ImportError:  # pragma: no cover - pybads is an optional dependency
+    HAS_BADS = False
+
+## three wells, each start at the bottom of one: -1 at 0.15, -3 at 0.5, -2 at 0.85
+CENTRES = np.array([0.15, 0.5, 0.85])
+DEPTHS = np.array([1.0, 3.0, 2.0])
+
+
+def wells(predicted, observed, **kwargs):
+    alpha = np.ravel(predicted)[0]
+    return -np.sum(DEPTHS * np.exp(-(((alpha - CENTRES) / 0.05) ** 2)))
+
+
+def identity(parameters, trial):
+    return {"dependent": np.array([parameters.alpha.value])}
+
+
+def build(optimiser_class, **kwargs):
+    data = pd.DataFrame({"ppt": [1, 2], "observed": [0.0, 0.0]})
+    model = Wrapper(
+        model=identity,
+        data=data[data.ppt == 1],
+        parameters=Parameters(alpha=Value(value=0.5, lower=0.0, upper=1.0, prior="norm", args={"mean": 0.5, "sd": 1.0})),
+    )
+    return optimiser_class(
+        model=model,
+        data=data,
+        minimisation=wells,
+        ppt_identifier="ppt",
+        initial_guess=[[0.15], [0.5], [0.85]],
+        **kwargs,
+    )
+
+
+CASES = [
+    pytest.param(Fmin, {}, id="Fmin"),
+    pytest.param(FminBound, {"approx_grad": True}, id="FminBound"),
+    pytest.param(Minimize, {}, id="Minimize"),
+]
+if HAS_BADS:
+    CASES.append(
+        pytest.param(
+            Bads,
+            {"options": {"uncertainty_handling": False, "display": "off", "max_fun_evals": 20}},
+            id="Bads",
+        )
+    )
+
+
+@pytest.mark.parametrize("hessian", ["numdifftools", "finite_differences"])
+@pytest.mark.parametrize("optimiser_class, kwargs", CASES)
+def test_the_best_of_three_starts_is_kept(optimiser_class, kwargs, hessian):
+    optimiser = build(optimiser_class, hessian=hessian, **kwargs)
+    optimiser.optimise()
+    for fit, parameters in zip(optimiser.fit, optimiser.parameters):
+        assert fit["fun"] == pytest.approx(-3.0, abs=1e-3)
+        assert parameters["alpha"] == pytest.approx(0.5, abs=1e-2)