Bayesian extend/resume

pyproject.toml

@@ -213,7 +213,7 @@ select = [
 # Ignore specific rules globally
 ignore = [
   'COM812', # https://docs.astral.sh/ruff/rules/missing-trailing-comma/
-  # The following is replaced by 'D'/[tool.ruff.lint.pydocstyle] and [tool.pydoclint]
+  # Replaced by 'D' plus pydocstyle and pydoclint
   'DOC', # https://docs.astral.sh/ruff/rules/#pydoclint-doc
   # Disable, as [tool.format_docstring] split one-line docstrings into the canonical multi-line layout
   'D200', # https://docs.astral.sh/ruff/rules/unnecessary-multiline-docstring/

src/easyscience/fitting/__init__.py

@@ -4,8 +4,10 @@
 from .available_minimizers import AvailableMinimizers
 from .fitter import Fitter
 from .minimizers.utils import FitResults
+from .sampler import Sampler
+from .sampler import SamplingResults
 
 # Causes circular import
 # from .multi_fitter import MultiFitter  # noqa: F401, E402
 
-all = [AvailableMinimizers, Fitter, FitResults]
+all = [AvailableMinimizers, Fitter, FitResults, Sampler, SamplingResults]

src/easyscience/fitting/minimizers/__init__.py

@@ -8,4 +8,11 @@
 from .utils import FitError
 from .utils import FitResults
 
-__all__ = [MinimizerBase, Bumps, DFO, LMFit, FitError, FitResults]
+__all__ = [
+    MinimizerBase,
+    Bumps,
+    DFO,
+    LMFit,
+    FitError,
+    FitResults,
+]

src/easyscience/fitting/minimizers/minimizer_bumps.py

@@ -4,6 +4,7 @@
 from __future__ import annotations
 
 import copy
+import math
 from typing import Any
 from typing import Callable
 
@@ -389,6 +390,7 @@
         burn: int = 2000,
         thin: int = 10,
         population: int | None = None,
+        resume_state: Any | None = None,
         sampler_kwargs: dict | None = None,
         progress_callback: Callable[[dict], bool | None] | None = None,
         abort_test: Callable[[], bool] | None = None,
@@ -411,13 +413,42 @@
         weights : np.ndarray
             Flattened weight array.
         samples : int, default=10000
-            Number of retained DREAM samples requested from BUMPS.
+            Number of raw samples to draw across all chains, before thinning.
+            A guaranteed minimum, not an exact count: DREAM advances in
+            blocks of 10 generations (one generation = one draw per chain)
+            and stops at the first block boundary at or past ``samples``.
         burn : int, default=2000
-            Burn-in steps.
+            Burn-in generations to discard. BUMPS counts ``burn`` in
+            generations while ``samples`` counts raw draws, so ``burn=500``
+            discards ``500 * n_chains`` raw samples.
         thin : int, default=10
-            Thinning interval.
+            Thinning interval — only every ``thin``-th generation is stored.
         population : int | None, default=None
-            BUMPS DREAM population count (number of parallel chains).
+            BUMPS DREAM population count per parameter (number of parallel
+            chains): BUMPS creates ``ceil(population * n_parameters)`` chains.
+        resume_state : Any | None, default=None
+            A BUMPS ``MCMCDraw`` state object from a previous
+            ``mcmc_sample()`` call (e.g. ``PosteriorResults.sampler_state``).
+            When provided, DREAM **continues** the saved chain instead of
+            starting cold.  The population, parameter count, and parameter
+            names must match the current model — a ``ValueError`` is raised
+            otherwise.
+
+            ``samples`` must be the **total** number of raw samples, not an
+            increment: to extend an existing chain of ``N`` raw samples by
+            ``M``, pass ``samples=N + M`` (DREAM keeps only the last
+            ``samples`` draws in its buffer). The `Sampler.extend` helper
+            computes this for you.
+
+            ``burn`` is forced to 0 on resume: a previously-converged chain is
+            never re-burned.
+
+            The ``population`` and ``initializer`` parameters
+            have **no effect** when ``resume_state`` is provided — they
+            are determined by the saved state.
+
+            Resuming against *different* data is undefined behaviour (the
+            chain's likelihood changes underneath it).
         sampler_kwargs : dict | None, default=None
             Additional keyword arguments forwarded to
             ``bumps.fitters.fit``.
@@ -439,8 +470,10 @@
         Raises
         ------
         ValueError
-            If the input shapes or weights are invalid, or if
-            ``progress_callback`` is not callable.
+            If the input shapes or weights are invalid, if
+            ``progress_callback`` is not callable, or if ``resume_state``
+            is incompatible with the current model (parameter count,
+            names/order, or population mismatch).
         FitError
             If DREAM sampling was aborted by the user (via
             ``abort_test``).
@@ -482,10 +515,100 @@
         curve = model_func(x, y, weights)
         problem = FitProblem(curve)
 
-        # Build DREAM kwargs
+        pop = population
+
+        # --- Resume-state compatibility validation ---
+        if resume_state is not None:
+            # Parameter count
+            n_fresh_params = len(problem._parameters)
+            state_nvar = resume_state.Nvar
+            if n_fresh_params != state_nvar:
+                raise ValueError(
+                    f'resume_state has {state_nvar} parameters but the current '
+                    f'model has {n_fresh_params}. The model must have the same '
+                    f'number of fitted parameters as when the saved chain was created.'
+                )
+
+            # Parameter names / order.
+            #
+            # BUMPS ``save_state``/``load_state`` does **not** preserve the
+            # original parameter labels: a reloaded ``MCMCDraw`` comes back
+            # with default labels like ``['P0', 'P1', ...]``. We can therefore
+            # only validate names when the state still carries our prefixed
+            # labels, i.e. an in-memory state from the current session. For a
+            # reloaded state we fall back to the parameter-count check (done
+            # above) plus an order-based warning, because BUMPS resumes
+            # positionally by column order, not by name.
+            fresh_names = [p.name[len(MINIMIZER_PARAMETER_PREFIX) :] for p in problem._parameters]
+            state_labels = list(resume_state.labels)
+            state_prefix = MINIMIZER_PARAMETER_PREFIX
+            labels_carry_our_names = bool(state_labels) and all(
+                lbl.startswith(state_prefix) for lbl in state_labels
+            )
+            if labels_carry_our_names:
+                state_stripped = [lbl[len(state_prefix) :] for lbl in state_labels]
+                if fresh_names != state_stripped:
+                    raise ValueError(
+                        f'Parameter names/order mismatch between the current model '
+                        f'and resume_state.\n'
+                        f'  Current model : {fresh_names}\n'
+                        f'  resume_state  : {state_stripped}'
+                    )
+            else:
+                from easyscience import global_object
+
+                global_object.log.getLogger('fitting.bumps').warning(
+                    'resume_state does not carry parameter names (it was most '
+                    'likely reloaded from disk, where BUMPS does not preserve '
+                    'labels). Parameter-name validation is skipped; the saved '
+                    'chain is matched to the current model by parameter order. '
+                    'Ensure this is the same model, with parameters in the same '
+                    'order, used to create the chain.'
+                )
+
+            # Population: BUMPS interprets ``pop`` as a **scale factor**.
+            # ``initpop.generate(problem, pop=P)`` creates
+            # ``ceil(P * n_params)`` chains, which becomes ``state.Npop``.
+            # On resume we must recover the scale factor from the state.
+            n_params = len(problem._parameters)
+            if pop is not None:
+                # Caller explicitly set population — validate it
+                # produces the same Npop as the saved state.
+                expected_npop = int(math.ceil(pop * n_params))
+                if expected_npop != resume_state.Npop:
+                    raise ValueError(
+                        f'Requested population ({pop}) would produce '
+                        f'{expected_npop} chains but the saved state has '
+                        f'{resume_state.Npop} chains. The population cannot '
+                        f'be changed on resume.'
+                    )
+            else:
+                # Recover the original scale factor from the state's Npop.
+                # Try ceil-division first, then fall back to floor.
+                recovered = int(math.ceil(resume_state.Npop / n_params))
+                if int(math.ceil(recovered * n_params)) != resume_state.Npop:
+                    recovered = int(resume_state.Npop / n_params)
+                pop = recovered
+
+            # Force burn=0 on resume: re-burning a previously converged
+            # chain discards good draws and is almost always a mistake.
+            if burn > 0:
+                from easyscience import global_object
+
+                global_object.log.getLogger('fitting.bumps').warning(
+                    f'burn={burn} ignored on resume: a previously converged '
+                    f'chain is not re-burned. Forcing burn=0.'
+                )
+            burn = 0
+
+        # Build DREAM kwargs. Use the resolved ``pop`` (the alias result, or
+        # the value recovered from the saved state on resume), not the raw
+        # ``population`` argument — on resume ``population`` is ``None`` while
+        # ``pop`` holds the scale factor that reproduces the saved state's
+        # population, which BUMPS requires to match.
         dream_kwargs: dict = {'samples': samples, 'burn': burn, 'thin': thin}
-        if population is not None:
-            dream_kwargs['pop'] = population
+        if pop is not None:
+            dream_kwargs['pop'] = pop
         if sampler_kwargs:
             dream_kwargs.update(sampler_kwargs)
 
@@ -497,7 +620,7 @@
             # Compute total DREAM steps for progress display (burn + sampling generations).
             # BUMPS DREAM default population count is 10 when not specified by the user.
             _dream_default_pop = 10
-            pop_val = population if population is not None else _dream_default_pop
+            pop_val = pop if pop is not None else _dream_default_pop
             _total_steps = burn + (samples + pop_val - 1) // pop_val
             monitors.append(
                 BumpsProgressMonitor(
@@ -525,7 +648,16 @@
         global_object.stack.enabled = False
 
         try:
-            x_opt, fx = driver.fit()
+            fit_kwargs = {}
+            if resume_state is not None:
+                # Defensive copy: BUMPS mutates the state object in-place
+                # (via MCMCDraw.resize() — see bumps/dream/core.py allocate_state)
+                # during resume.  Without a copy, the caller's original state
+                # object is silently altered, making it impossible to compare
+                # pre- and post-resume state (shape mismatch).  See
+                # https://github.com/easyscience/core/pull/257
+                fit_kwargs['fit_state'] = copy.deepcopy(resume_state)
+            x_opt, fx = driver.fit(**fit_kwargs)
             result_state = getattr(driver.fitter, 'state', None)
             if result_state is None:
                 raise FitError('Sampling aborted by user')
@@ -535,9 +667,10 @@
         finally:
             global_object.stack.enabled = stack_status
 
-        draws = result_state.draw().points
+        _draw = result_state.draw()
+        draws = _draw.points
         param_names = [p.name[len(MINIMIZER_PARAMETER_PREFIX) :] for p in problem._parameters]
-        logp = getattr(result_state, 'logp', None)
+        logp = _draw.logp
 
         return {
             'draws': draws,

src/easyscience/fitting/multi_fitter.py

@@ -1,8 +1,6 @@
 # SPDX-FileCopyrightText: 2026 EasyScience contributors <https://github.com/easyscience>
 # SPDX-License-Identifier: BSD-3-Clause
 
-from __future__ import annotations
-
 from typing import Callable
 
 import numpy as np