Bayesian post-analysis using BUMPS DREAM

src/easyscience/fitting/minimizers/minimizer_bumps.py

@@ -334,6 +334,90 @@ def _make_func(x, y, weights):
 
         return _outer(self)
 
+    def sample(
+        self,
+        x: np.ndarray,
+        y: np.ndarray,
+        weights: np.ndarray,
+        samples: int = 10000,
+        burn: int = 2000,
+        thin: int = 10,
+        chains: int | None = None,
+        population: int | None = None,
+        seed: int | None = None,
+        sampler_kwargs: dict | None = None,
+    ) -> dict:
+        """Run Bayesian MCMC sampling using the BUMPS DREAM sampler.
+
+        Builds a BUMPS :class:`~bumps.names.FitProblem` from the current
+        model and runs the DREAM sampler.  This is the public minimizer-level
+        entry point for Bayesian sampling; the higher-level
+        :meth:`easyscience.fitting.multi_fitter.MultiFitter.sample` delegates
+        to this method after flattening multi-dataset arrays.
+
+        :param x: Flattened independent variable array.
+        :type x: np.ndarray
+        :param y: Flattened dependent variable array.
+        :type y: np.ndarray
+        :param weights: Flattened weight array.
+        :type weights: np.ndarray
+        :param samples: Number of retained DREAM samples requested from BUMPS.
+        :type samples: int
+        :param burn: Burn-in steps.
+        :type burn: int
+        :param thin: Thinning interval.
+        :type thin: int
+        :param chains: User-friendly alias for BUMPS DREAM population count.
+        :type chains: int | None
+        :param population: BUMPS DREAM population count (``pop``) for advanced users.
+        :type population: int | None
+        :param seed: Best-effort random seed passed to ``numpy.random.seed``.
+            BUMPS DREAM may use additional internal RNG state that is not
+            controlled by this seed, so exact reproducibility is not guaranteed.
+        :type seed: int | None
+        :param sampler_kwargs: Additional keyword arguments forwarded to
+            :func:`bumps.fitters.fit`.
+        :type sampler_kwargs: dict | None
+        :return: Dictionary with keys ``'draws'``, ``'param_names'``, ``'state'``,
+            and ``'logp'``.
+        :rtype: dict
+        """
+        from bumps.fitters import fit as bumps_fit
+        from bumps.names import FitProblem
+
+        # Build the BUMPS Curve model using the minimizer's existing machinery
+        model_func = self._make_model()
+        x_flat = np.linspace(0, y.size - 1, y.size)
+        curve = model_func(x_flat, y, weights)
+        problem = FitProblem(curve)
+
+        # Best-effort seed: sets numpy's global RNG state just before DREAM starts.
+        # BUMPS DREAM may have its own internal RNG paths that are not fully
+        # controlled by this, so exact reproducibility is not guaranteed.
+        if seed is not None:
+            np.random.seed(seed)
+
+        # Run DREAM sampler
+        dream_kwargs: dict = {'samples': samples, 'burn': burn, 'thin': thin}
+        if chains is not None or population is not None:
+            pop = chains if chains is not None else population
+            dream_kwargs['pop'] = pop
+        if sampler_kwargs:
+            dream_kwargs.update(sampler_kwargs)
+        result = bumps_fit(problem, method='dream', **dream_kwargs)
+
+        # Extract posterior
+        draws = result.state.draw().points
+        param_names = [p.name[len(MINIMIZER_PARAMETER_PREFIX) :] for p in problem._parameters]
+        logp = getattr(result.state, 'logp', None)
+
+        return {
+            'draws': draws,
+            'param_names': param_names,
+            'state': result.state,
+            'logp': logp,
+        }
+
     def _set_parameter_fit_result(
         self,
         fit_result,

src/easyscience/fitting/multi_fitter.py

@@ -2,6 +2,7 @@
 # SPDX-License-Identifier: BSD-3-Clause
 
 from typing import Callable
+from typing import Dict
 from typing import List
 from typing import Optional
 
@@ -149,3 +150,109 @@
             fit_results_list.append(current_results)
             sp = ep
         return fit_results_list
+
+    def sample(
+        self,
+        x: List[np.ndarray],
+        y: List[np.ndarray],
+        weights: List[np.ndarray],
+        samples: int = 10000,
+        burn: int = 1000,
+        thin: int = 10,
+        chains: int | None = None,
+        population: int | None = None,
+        seed: int | None = None,
+        vectorized: bool = False,
+        sampler_kwargs: dict | None = None,
+    ) -> Dict:
+        """Run Bayesian MCMC sampling using the BUMPS DREAM sampler.
+
+        Requires that the current minimizer is a BUMPS instance (i.e. the
+        minimizer was switched to ``AvailableMinimizers.Bumps`` or equivalent).
+
+        :param x: List of independent variable arrays (one per dataset).
+        :type x: List[np.ndarray]
+        :param y: List of dependent variable arrays (one per dataset).
+        :type y: List[np.ndarray]
+        :param weights: List of weight arrays (one per dataset).
+        :type weights: List[np.ndarray]
+        :param samples: Number of retained DREAM samples requested from BUMPS.
+        :type samples: int
+        :param burn: Burn-in steps.
+        :type burn: int
+        :param thin: Thinning interval.
+        :type thin: int
+        :param chains: User-friendly alias for BUMPS DREAM population count.
+        :type chains: int | None
+        :param population: BUMPS DREAM population count (``pop``) for advanced users.
+        :type population: int | None
+        :param seed: Best-effort random seed. BUMPS DREAM may use additional
+            internal RNG state that is not controlled by this seed, so exact
+            reproducibility is not guaranteed.
+        :type seed: int | None
+        :param vectorized: Whether the fit function expects vectorized
+            (multidimensional) input. Defaults to ``False``.
+        :type vectorized: bool
+        :param sampler_kwargs: Additional keyword arguments forwarded to the
+            BUMPS DREAM sampler via :func:`bumps.fitters.fit`.
+        :type sampler_kwargs: dict | None
+        :return: Dictionary with keys ``'draws'``, ``'param_names'``, ``'state'``,
+            and ``'logp'``.
+        :rtype: dict
+        :raises RuntimeError: If the current minimizer is not a BUMPS instance.
+        """
+        # --- Alias resolution ---
+        if chains is not None and population is not None:
+            if chains != population:
+                raise ValueError(
+                    f'Conflicting population arguments: chains={chains}, population={population}. '
+                    'Only provide one.'
+                )
+            pop = chains
+        elif chains is not None:
+            pop = chains
+        elif population is not None:
+            pop = population
+        else:
+            pop = None
+
+        # Flatten multi-dataset arrays
+        _, x_new, y_new, w_new, _dims = self._precompute_reshaping(
+            x, y, weights, vectorized=vectorized
+        )
+        self._dependent_dims = _dims
+
+        # Wrap fit functions for multi-dataset flattening, mirroring the
+        # ``Fitter.fit`` lifecycle: use the property setter so the minimizer
+        # is re-created with the wrapped fit function.
+        original_fit_func = self.fit_function
+        fit_fun_wrap = self._fit_function_wrapper(x_new, flatten=True)
+        self.fit_function = fit_fun_wrap
+
+        try:
+            minimizer = self.minimizer
+
+            # Verify it's a BUMPS minimizer (sampling only works with BUMPS/DREAM)
+            if not (hasattr(minimizer, 'package') and minimizer.package == 'bumps'):
+                raise RuntimeError(
+                    'Bayesian sampling requires a BUMPS minimizer. '
+                    'Use ``fitter.switch_minimizer(AvailableMinimizers.Bumps)`` first.'
+                )
+
+            # Delegate to the BUMPS minimizer's public sample method
+            result = minimizer.sample(
+                x=x_new,
+                y=y_new,
+                weights=w_new,
+                samples=samples,
+                burn=burn,
+                thin=thin,
+                chains=None,  # alias already resolved into `pop`
+                population=pop,
+                seed=seed,
+                sampler_kwargs=sampler_kwargs,
+            )
+        finally:
+            self.fit_function = original_fit_func
+
+        return result