diff --git a/CHANGELOG.md b/CHANGELOG.md index 870098a2..53683f1c 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/cpm/generators/__init__.py b/cpm/generators/__init__.py index b9271b50..9b8c3353 100644 --- a/cpm/generators/__init__.py +++ b/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 diff --git a/cpm/hierarchical/variational.py b/cpm/hierarchical/variational.py index a0453740..e0ecda65 100644 --- a/cpm/hierarchical/variational.py +++ b/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 diff --git a/cpm/optimisation/bads.py b/cpm/optimisation/bads.py index 5b972fcf..2ae8dc53 100644 --- a/cpm/optimisation/bads.py +++ b/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( diff --git a/cpm/optimisation/fmin.py b/cpm/optimisation/fmin.py index a916cf83..c506bdf4 100644 --- a/cpm/optimisation/fmin.py +++ b/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( diff --git a/cpm/optimisation/free.py b/cpm/optimisation/free.py index 304492ad..97393f3e 100644 --- a/cpm/optimisation/free.py +++ b/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( diff --git a/test/optimisation/test_parallel.py b/test/optimisation/test_parallel.py new file mode 100644 index 00000000..53c00977 --- /dev/null +++ b/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"] diff --git a/test/optimisation/test_starts.py b/test/optimisation/test_starts.py new file mode 100644 index 00000000..072c7488 --- /dev/null +++ b/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)