Twinner

The MAP optimiser behind ReplayBG.twin() (multi-start Powell, polish of the n_polish best distinct starts, optional parallelism). Used directly, Twinner defaults to parallelize=True (unlike ReplayBG.twin(), which defaults to False), exposes polish=False to skip the polish, and twin() returns {'fun', 'x'} (the best negative log-posterior and parameter vector).

A class that estimates model parameters by maximizing a posterior objective.

The twinner performs MAP (maximum a posteriori) estimation: it minimizes the negative log-posterior (-(log_prior + log_likelihood)) over the unknown parameters using a multi-start Powell optimisation, optionally parallelised across CPUs with multiprocessing.

...

Attributes:

Name	Type	Description
parallelize	bool	A boolean that specifies whether to run the multi-start optimisation in parallel.
n_jobs	int	The number of parallel jobs to use (-1 means all available CPUs).
n_starts	int	The number of multi-start optimisations.
verbose	bool	A boolean that specifies the verbosity of the twinner.
polish	bool	A boolean that specifies whether the best multi-start results are refined with a tight-tolerance polish.
n_polish	int	The number of best distinct multi-start results that are polished.
log_history	bool	A boolean that specifies whether the optimisation history is recorded.
history	dict	The recorded optimisation history (only present when log_history is True), with keys theta, log_prior, log_likelihood and log_posterior.

Methods:

Name	Description
twin	Runs the multi-start MAP estimation and returns the best parameters.

Source code in src/py_replay_bg/twinner/twinner.py

class Twinner():
    """A class that estimates model parameters by maximizing a posterior objective.

    The twinner performs MAP (maximum a posteriori) estimation: it minimizes the
    negative log-posterior (``-(log_prior + log_likelihood)``) over the unknown
    parameters using a multi-start Powell optimisation, optionally parallelised
    across CPUs with multiprocessing.

    ...
    Attributes
    ----------
    parallelize : bool
        A boolean that specifies whether to run the multi-start optimisation in
        parallel.
    n_jobs : int
        The number of parallel jobs to use (``-1`` means all available CPUs).
    n_starts : int
        The number of multi-start optimisations.
    verbose : bool
        A boolean that specifies the verbosity of the twinner.
    polish : bool
        A boolean that specifies whether the best multi-start results are
        refined with a tight-tolerance polish.
    n_polish : int
        The number of best distinct multi-start results that are polished.
    log_history : bool
        A boolean that specifies whether the optimisation history is recorded.
    history : dict
        The recorded optimisation history (only present when ``log_history`` is
        ``True``), with keys ``theta``, ``log_prior``, ``log_likelihood`` and
        ``log_posterior``.

    Methods
    -------
    twin(model, rbg_data, unknown_parameters_prior, correlations=None):
        Runs the multi-start MAP estimation and returns the best parameters.
    """

    def __init__(self,
                 parallelize: bool = True,
                 n_jobs: int | None = None,
                 n_starts: int = 64,
                 log_history: bool = False,
                 verbose: bool = True,
                 polish: bool = True,
                 n_polish: int = 8,
    ):
        """Constructs all the necessary attributes for the Twinner object.

        Parameters
        ----------
        parallelize : bool, optional, default : True
            A boolean that specifies whether to run the multi-start optimisation
            in parallel using multiprocessing.
        n_jobs : int or None, optional, default : None
            The number of parallel jobs to use. If ``None``, all available CPUs
            are used (stored internally as ``-1``).
        n_starts : int, optional, default : 64
            An integer that specifies the number of multi-start optimisations.
        log_history : bool, optional, default : False
            A boolean that specifies whether to record the optimisation history.
        verbose : bool, optional, default : True
            A boolean that specifies the verbosity of the twinner.
        polish : bool, optional, default : True
            A boolean that specifies whether to run a deterministic tight-tolerance
            refinement of the ``n_polish`` best distinct multi-start results.
            Powell's default line-search tolerance stops early on the
            ill-conditioned (sloppy) physiological likelihood; the polish restarts
            Powell from each of them with tight tolerances until the objective
            stops improving, and the best polished result is kept. Off preserves
            the raw multi-start behaviour.
        n_polish : int, optional, default : 8
            The number of best multi-start results (with distinct objectives)
            that are polished; the best polished result is kept. The raw
            default-tolerance ranking is a poor guide to which start polishes
            into the best optimum, so polishing several starts is far more
            reliable than polishing only the best one. The polishes run in
            parallel when ``parallelize`` is ``True``. ``1`` polishes only the
            best start. Ignored when ``polish`` is ``False``.
        """
        if n_polish < 1:
            raise ValueError("n_polish must be at least 1")
        self.parallelize = parallelize
        self.n_jobs = -1 if n_jobs is None else n_jobs
        self.n_starts = n_starts
        self.verbose = verbose
        self.polish = polish
        self.n_polish = n_polish

        self.log_history = log_history
        if self.log_history:
            self.history = dict()
            self.history['theta'] = []
            self.history['log_prior'] = []
            self.history['log_likelihood'] = []
            self.history['log_posterior'] = []

    def twin(self, model : Any, rbg_data, unknown_parameters_prior, correlations=None) -> dict:
        """Runs the multi-start MAP estimation for a model.

        Builds ``n_starts`` initial guesses by sampling the parameter priors,
        runs a Powell optimisation of the negative log-posterior (bounds enforced
        through the prior) from each start (sequentially or in parallel),
        polishes the ``n_polish`` best starts with distinct objectives (when
        ``polish`` is on, sequentially or in parallel), and selects the result
        with the lowest objective. Integer
        parameters are rounded and all parameters are clipped to their bounds
        before being returned.

        Parameters
        ----------
        model : object
            A model instance being fit, implementing ``reset(theta0, g0, u0)``,
            ``step(u, t)`` and ``output(t)``.
        rbg_data : object
            A data object used to compute the likelihood.
        unknown_parameters_prior : dict
            A dictionary describing the priors and bounds for each parameter to
            estimate (keys ``prior``, ``min``, ``max`` and optional ``integer``).
        correlations : dict or None, optional, default : None
            Optional pairwise prior correlations specified as
            ``{(name_a, name_b): rho}``. When given, the joint prior is a
            Gaussian copula over the named parameters (marginals unchanged); when
            ``None`` the parameters are treated as independent.

        Returns
        -------
        dict
            A dictionary with keys ``fun`` (the best negative log-posterior
            value) and ``x`` (the best-fit parameter vector, rounded and clipped
            to the bounds).
        """

        # Build the (validated) Gaussian-copula correlation structure once, up
        # front, so any error surfaces before the optimisation starts.
        correlation_structure = _build_correlation_structure(unknown_parameters_prior, correlations)

        # Set the worker arguments
        global _worker_args
        _worker_args = (self._neg_log_posterior, self._log_posterior, unknown_parameters_prior, model, rbg_data,
                        correlation_structure)

        # Build initial guesses for the parameters using their priors
        start_guesses = []
        for i in range(self.n_starts):
            start_guess = np.array([v['prior'].sample(v['min'], v['max'], i)
                                    for v in unknown_parameters_prior.values()])
            start_guesses.append((i, start_guess))

        # Resolve the effective number of jobs (used both for running and for the header)
        n_jobs = multiprocessing.cpu_count() if self.n_jobs == -1 else self.n_jobs

        # Print a run header so the user can sanity-check the twinning setup
        if self.verbose:
            param_names = list(unknown_parameters_prior.keys())
            mode = f'parallel ({n_jobs} jobs)' if self.parallelize else 'sequential'
            print('--- Twinning (MAP estimation) ---')
            print(f'  starts          : {self.n_starts}')
            print(f'  mode            : {mode}')
            print(f'  polish          : ' + (f'top {self.n_polish} distinct starts' if self.polish else 'off'))
            print(f'  parameters ({len(param_names)})  : {", ".join(param_names)}')
            if correlation_structure is not None:
                pairs = ", ".join(f'{a}~{b}: {rho:g}'
                                  for (a, b), rho in correlations.items())
                print(f'  correlations    : {pairs}')

        t0 = time.perf_counter()

        if self.parallelize:
            # Set the context
            ctx = multiprocessing.get_context('fork' if os.name != 'nt' else 'spawn')
            # Run the optimization in parallel
            with ctx.Pool(processes=n_jobs) as pool:
                raw = list(tqdm(
                    pool.imap(_run_optimization, start_guesses),
                    total=self.n_starts,
                    desc='Twinning using MAP',
                    disable=not self.verbose
                ))
            results = [r for r, _ in raw]
            if self.log_history:
                for _, h in raw:
                    if h is not None:
                        for k in self.history:
                            self.history[k].extend(h[k])
        else:
            # Run the optimization sequentially; history accumulates in-process via _neg_log_posterior
            raw = [_run_optimization(a)
                   for a in tqdm(start_guesses, desc='Twinning', disable=not self.verbose)]
            results = [r for r, _ in raw]

        # Refine the best starts with a deterministic tight-tolerance polish. The
        # multi-start Powell runs at scipy's default tolerance, which halts early
        # on the sloppy (ill-conditioned) physiological likelihood ridge; a
        # tight-tolerance restart loop walks the rest of the way down to the
        # optimum. The raw ranking is a poor guide to which start polishes into
        # the best optimum, so the n_polish best distinct starts are polished and
        # the best polished result is kept. No RNG is involved, so twinning stays
        # reproducible.
        candidates = _top_distinct(results, self.n_polish) if self.polish else []
        if candidates:
            # Copies, so the raw results (used in the summary) stay untouched
            polish_args = [copy.copy(results[i]) for i in candidates]
            if self.parallelize and len(candidates) > 1:
                ctx = multiprocessing.get_context('fork' if os.name != 'nt' else 'spawn')
                with ctx.Pool(processes=min(n_jobs, len(candidates))) as pool:
                    raw = list(tqdm(
                        pool.imap(_run_polish, polish_args),
                        total=len(candidates),
                        desc='Polishing',
                        disable=not self.verbose
                    ))
                polished = [r for r, _ in raw]
                if self.log_history:
                    for _, h in raw:
                        if h is not None:
                            for k in self.history:
                                self.history[k].extend(h[k])
            else:
                polished = [self._polish(r, model, rbg_data, unknown_parameters_prior, correlation_structure)
                            for r in tqdm(polish_args, desc='Polishing', disable=not self.verbose)]
            # min keeps the first on ties, i.e. the better raw rank
            best_rank = min(range(len(polished)), key=lambda j: polished[j].fun)
            best = polished[best_rank]
        else:
            best_rank = None
            best = min(results, key=lambda r: r.fun)

        elapsed = time.perf_counter() - t0

        # Print a convergence summary across all starts
        if self.verbose:
            funs = np.array([r.fun for r in results], dtype=float)
            finite = funs[np.isfinite(funs)]
            n_failed = len(funs) - len(finite)
            print('--- Twinning summary ---')
            print(f'  elapsed         : {elapsed:.1f} s')
            print(f'  failed starts   : {n_failed} / {len(funs)} (non-finite objective)')
            if len(finite):
                print(f'  objective (neg log-posterior, lower is better)')
                print(f'    raw starts    : min {finite.min():.3f} / med {np.median(finite):.3f} / '
                      f'max {finite.max():.3f}')
                if candidates:
                    pol = np.array([r.fun for r in polished], dtype=float)
                    winner = candidates[best_rank]
                    print(f'    polished ({len(candidates)})  : best {pol.min():.3f} / worst {pol.max():.3f}')
                    print(f'    final         : {best.fun:.3f}  (start {winner}, raw {results[winner].fun:.3f}, '
                          f'{_ordinal(best_rank + 1)} best distinct raw start)')
                else:
                    print(f'    final         : {best.fun:.3f}  (not polished)')
            else:
                print('  WARNING: no start produced a finite objective '
                      '(check priors/bounds and the model).')

        # Round integer parameters and clip to bounds (Powell has no native bound support)
        clipped_x = []
        for x, (k, v) in zip(best.x, unknown_parameters_prior.items()):
            if v.get('integer', False):
                x = float(round(x))
            clipped_x.append(np.clip(x, v['min'], v['max']))
        clipped_x = np.array(clipped_x)

        # Print the best-fit parameters, flagging any that sit at their bounds
        if self.verbose:
            print('--- Best-fit parameters ---')
            for x, (k, v) in zip(clipped_x, unknown_parameters_prior.items()):
                lo, hi = v['min'], v['max']
                tol = 1e-6 * max(1.0, abs(hi - lo))
                flags = []
                if x <= lo + tol:
                    flags.append('at min bound')
                elif x >= hi - tol:
                    flags.append('at max bound')
                if v.get('integer', False):
                    flags.append('integer')
                suffix = f'   [{", ".join(flags)}]' if flags else ''
                print(f'  {k:<14} = {x:.4g}   (bounds: {lo:g}..{hi:g}){suffix}')

        # Return the best result
        ret = dict()
        ret['fun'] = best.fun
        ret['x'] = clipped_x
        return ret

    def _polish(self, best, model, rbg_data, unknown_parameters_prior,
                correlation_structure=None, max_rounds: int = 5, tol: float = 1e-8,
                rel_improve: float = 1e-7, evals_per_dim: int = 1000,
                max_budget: int = 8000):
        """Refines the best multi-start result with a tight-tolerance Powell loop.

        Powell resets its search-direction set on each fresh call, so restarting
        it from the incumbent with tight tolerances lets it keep descending a
        narrow, ill-conditioned valley that a single default-tolerance run stops
        short of. The loop repeats until a round fails to improve the objective by
        a meaningful *relative* amount, so it stops promptly once the optimum is
        reached. The starting point is fixed and no randomness is involved, so the
        result is deterministic.

        Cost is bounded on purpose: each restart's function-evaluation budget is
        ``min(evals_per_dim * n_params, max_budget)`` and the number of restarts is
        capped, with an early stop once a round stops improving. The absolute
        ``max_budget`` cap keeps a single restart affordable on high-dimensional
        models (e.g. the 25-parameter extended model, where an uncapped tight
        Powell runs for many minutes), while the ``max_rounds`` restarts let the
        cheaper low-dimensional models still converge fully. The bounds depend only
        on fixed problem sizes, so the result stays deterministic.

        Parameters
        ----------
        best : scipy.optimize.OptimizeResult
            The best multi-start result; its ``x``/``fun`` seed the refinement.
        model : object
            The model instance being fit (mutated in place by each evaluation).
        rbg_data : object
            The data object used to compute the likelihood.
        unknown_parameters_prior : dict
            The priors/bounds spec (fixes parameter order and integer handling).
        correlation_structure : dict or None, optional, default : None
            The precomputed Gaussian-copula structure, or ``None``.
        max_rounds : int, optional, default : 5
            Maximum number of tight-tolerance Powell restarts (the early stop
            often triggers before this).
        tol : float, optional, default : 1e-8
            The ``xtol``/``ftol`` passed to each Powell restart. Far tighter than
            scipy's 1e-4 default (which caused the under-convergence).
        rel_improve : float, optional, default : 1e-7
            Minimum relative objective improvement ``(fun - r.fun)/max(1, |fun|)``
            a round must yield to keep going; below it the loop stops.
        evals_per_dim : int, optional, default : 1000
            Per-restart function-evaluation budget per parameter (before the cap).
        max_budget : int, optional, default : 8000
            Absolute cap on a single restart's function-evaluation budget; bounds
            the per-restart cost on high-dimensional models.

        Returns
        -------
        scipy.optimize.OptimizeResult
            ``best`` with ``x``/``fun`` updated to the refined optimum (unchanged
            when no restart improves on the incumbent).
        """
        x = np.asarray(best.x, dtype=float)
        fun = best.fun
        budget = min(evals_per_dim * len(x), max_budget)
        for _ in range(max_rounds):
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                r = minimize(self._neg_log_posterior, x, method='Powell',
                             args=(model, rbg_data, unknown_parameters_prior,
                                   correlation_structure),
                             options={'maxiter': budget, 'maxfev': budget,
                                      'xtol': tol, 'ftol': tol, 'disp': False})
            if not np.isfinite(r.fun) or r.fun >= fun:
                break
            gain = (fun - r.fun) / max(1.0, abs(fun))
            x, fun = r.x, r.fun
            if gain < rel_improve:
                break
        best.x = x
        best.fun = fun
        return best

    def _log_prior(self, model, unknown_parameters_prior, correlation_structure=None):
        """Computes the log prior probability of the model parameters.

        The base term is the sum of marginal log-densities. When a Gaussian-copula
        ``correlation_structure`` is supplied, a copula correction term is added so
        that the named parameters are jointly correlated while keeping their
        marginals unchanged.

        Parameters
        ----------
        model : object
            A model instance containing the current parameter values.
        unknown_parameters_prior : dict
            A dictionary describing the priors for each parameter. Each prior
            must provide a ``log_pdf(value)`` method.
        correlation_structure : dict or None, optional, default : None
            Precomputed Gaussian-copula structure (see
            :func:`_build_correlation_structure`). When ``None`` the parameters
            are treated as independent.

        Returns
        -------
        float
            The log prior probability, or ``-inf`` when any parameter is outside
            its valid range.
        """
        # Iterate over the parameters and compute the marginal log prior
        lp = 0.0
        for up, v in unknown_parameters_prior.items():
            parameter_value = getattr(model, up)
            a = v["min"]
            b = v["max"]
            # If the parameter is outside the valid range, return -inf
            if parameter_value <= a or parameter_value >= b:
                return -np.inf
            # Otherwise, add the marginal log prior contribution (corrected for the truncation)
            log_density = v["prior"].log_pdf(parameter_value)
            z = v["prior"].cdf(b) - v["prior"].cdf(a)
            if log_density == -np.inf or z <= 0:
                return -np.inf
            lp += log_density - np.log(z)

        # Add the Gaussian-copula correction for correlated parameters (if any)
        if correlation_structure is not None:
            names = correlation_structure['names']
            # Map each correlated parameter to a standard-normal score z = Phi^-1(F(x))
            z_scores = np.empty(len(names))
            for i, name in enumerate(names):
                v = unknown_parameters_prior[name]
                x = getattr(model, name)
                fa = v["prior"].cdf(v["min"])
                fb = v["prior"].cdf(v["max"])
                fx = v["prior"].cdf(x)
                u = (fx - fa) / (fb - fa)
                u = np.clip(u, 1e-12, 1.0 - 1e-12)
                z_scores[i] = ndtri(u)
            # log copula density: -0.5*logdet(R) - 0.5*z^T (R^-1 - I) z
            lp += -0.5 * correlation_structure['logdet_R'] \
                  - 0.5 * float(z_scores @ correlation_structure['R_inv_minus_I'] @ z_scores)

        # Return the (possibly correlation-corrected) log prior
        return lp

    def _log_likelihood(self, model, rbg_data):
        """Computes the log likelihood of the observed data under the model.

        The model is simulated forward using the input sequence in ``rbg_data``.
        The predicted output is then compared against the observed glucose
        measurements using a Gaussian error model with a constant coefficient of
        variation.

        Parameters
        ----------
        model : object
            A model instance to simulate.
        rbg_data : object
            A data object containing inputs, output timestamps, and observed
            glucose values.

        Returns
        -------
        float
            The log likelihood value for the simulated trajectory.
        """
        # Simulate the model forward and get the output
        out = np.zeros(rbg_data.tsteps, )
        out[0] = model.output(0)
        for k in np.arange(1,out.shape[0]):
            model.step(rbg_data.u[k], k)
            out[k] = model.output(k)
        # Subsample the output to match the sampling rate of the data
        out = out[0::rbg_data.yts]

        # TODO: enable the choice of multiple ll shapes
        # Two-term noise model: σ² = σ_add² + (cv·|ŷ|)²
        #sigma_add = 5.0          # mg/dL, literature default
        #cv        = 0.05         # 5%
        #sdn = np.sqrt(sigma_add**2 + (cv * np.abs(out[rbg_data.y_idxs]))**2)
        #residuals = out[rbg_data.y_idxs] - rbg_data.y[rbg_data.y_idxs]

        # Calculate the log-likelihood with a Gaussian error model (constant coefficient of variation)
        cv = 0.05  # constant coefficient of variation (5%) #TODO: make this a parameter
        residuals = out[rbg_data.y_idxs] - rbg_data.y[rbg_data.y_idxs]
        sdn = cv * np.abs(out[rbg_data.y_idxs])
        subsampling = 1 #TODO: decide where to evaluate ll (now every 30 minutes but maybe we can do something smarter)
        sdn = sdn[0::subsampling]
        residuals = residuals[0::subsampling]
        n = len(residuals)
        return -0.5 * np.sum((residuals / sdn) ** 2) - np.sum(np.log(sdn)) - 0.5 * n * np.log(2 * np.pi)

    def _neg_log_posterior(self, theta, model, rbg_data, unknown_parameters_prior, correlation_structure=None):
        """Returns the negative log posterior for optimization.

        Parameters
        ----------
        theta : numpy.ndarray
            A parameter vector in the natural (constrained) parameter space,
            ordered to match ``unknown_parameters_prior``.
        model : object
            A model instance being fit.
        rbg_data : object
            A data object used to compute the likelihood.
        unknown_parameters_prior : dict
            A dictionary describing priors and bounds for the parameters.

        Returns
        -------
        float
            The negative log posterior value.
        """
        # Just return the negative log-posterior
        log_prior, log_likelihood, log_post = self._log_posterior(theta, model, rbg_data, unknown_parameters_prior,
                                                                  correlation_structure)

        # log the history
        if self.log_history:
            self.history['theta'].append(theta)
            self.history['log_prior'].append(log_prior)
            self.history['log_likelihood'].append(log_likelihood)
            self.history['log_posterior'].append(log_post)

        return -log_post

    def _log_posterior(self, theta, model, rbg_data, unknown_parameters_prior, correlation_structure=None):
        """Compute log-prior, log-likelihood, and log-posterior in one model pass.

        The parameter vector is used directly in the natural (constrained)
        parameter space: each entry is written into the model's typed parameter
        dict (integer parameters are rounded) before ``model.reset``. Every reset
        also re-seeds the cold-start glucose and the t=0 input row from
        ``rbg_data`` (see ``Model.reset``), since ``reset`` clears them each call.
        Bounds are enforced by the prior, which returns ``-inf`` outside
        ``[min, max]``.

        Parameters
        ----------
        theta : numpy.ndarray
            A parameter vector in the natural (constrained) parameter space,
            ordered to match ``unknown_parameters_prior``.
        model : object
            A model instance being updated with candidate parameters.
        rbg_data : object
            A data object used to compute the likelihood.
        unknown_parameters_prior : dict
            A dictionary describing priors and bounds for the parameters.

        Returns
        -------
        tuple
            The triple ``(log_prior, log_likelihood, log_posterior)``. Any
            component that is invalid is returned as ``-np.inf``.
        """

        # Create the theta input dictionary (note: order is maintained by construction)
        theta_dict = Dict.empty(key_type=types.unicode_type, value_type=float64)
        for i, k in enumerate(unknown_parameters_prior.keys()):
            val = theta[i]
            #val = np.clip(theta[i], unknown_parameters_prior[k]['min'], unknown_parameters_prior[k]['max'])
            is_int =  unknown_parameters_prior[k].get('integer', False)
            if is_int:
                val = int(round(val))
            theta_dict[k] = val

        # Reset the model with the new parameters, re-seeding the cold-start
        # glucose and the t=0 input row from data every time (reset() clears
        # both, so they can't just be seeded once outside this loop).
        g0 = float(rbg_data.y[rbg_data.y_idxs[0]]) if len(rbg_data.y_idxs) > 0 else np.nan
        model.reset(theta_dict, g0, rbg_data.u[0])

        # Calculate log-prior
        lp = self._log_prior(model, unknown_parameters_prior, correlation_structure)
        if lp == -np.inf or np.isnan(lp):
            return -np.inf, -np.inf, -np.inf
        # Calculate log-likelihood
        ll = self._log_likelihood(model, rbg_data)
        if ll == -np.inf or np.isnan(ll):
            return lp, -np.inf, -np.inf

        # Return log-posterior
        return lp, ll, lp + ll

parallelize = parallelize instance-attribute

n_jobs = -1 if n_jobs is None else n_jobs instance-attribute

n_starts = n_starts instance-attribute

verbose = verbose instance-attribute

polish = polish instance-attribute

n_polish = n_polish instance-attribute

log_history = log_history instance-attribute

history = dict() instance-attribute

__init__(parallelize=True, n_jobs=None, n_starts=64, log_history=False, verbose=True, polish=True, n_polish=8)

Constructs all the necessary attributes for the Twinner object.

Parameters:

Name	Type	Description	Default
parallelize	bool	A boolean that specifies whether to run the multi-start optimisation in parallel using multiprocessing.	: True
n_jobs	int or None	The number of parallel jobs to use. If None, all available CPUs are used (stored internally as -1).	: None
n_starts	int	An integer that specifies the number of multi-start optimisations.	: 64
log_history	bool	A boolean that specifies whether to record the optimisation history.	: False
verbose	bool	A boolean that specifies the verbosity of the twinner.	: True
polish	bool	A boolean that specifies whether to run a deterministic tight-tolerance refinement of the n_polish best distinct multi-start results. Powell's default line-search tolerance stops early on the ill-conditioned (sloppy) physiological likelihood; the polish restarts Powell from each of them with tight tolerances until the objective stops improving, and the best polished result is kept. Off preserves the raw multi-start behaviour.	: True
n_polish	int	The number of best multi-start results (with distinct objectives) that are polished; the best polished result is kept. The raw default-tolerance ranking is a poor guide to which start polishes into the best optimum, so polishing several starts is far more reliable than polishing only the best one. The polishes run in parallel when parallelize is True. 1 polishes only the best start. Ignored when polish is False.	: 8

Source code in src/py_replay_bg/twinner/twinner.py

def __init__(self,
             parallelize: bool = True,
             n_jobs: int | None = None,
             n_starts: int = 64,
             log_history: bool = False,
             verbose: bool = True,
             polish: bool = True,
             n_polish: int = 8,
):
    """Constructs all the necessary attributes for the Twinner object.

    Parameters
    ----------
    parallelize : bool, optional, default : True
        A boolean that specifies whether to run the multi-start optimisation
        in parallel using multiprocessing.
    n_jobs : int or None, optional, default : None
        The number of parallel jobs to use. If ``None``, all available CPUs
        are used (stored internally as ``-1``).
    n_starts : int, optional, default : 64
        An integer that specifies the number of multi-start optimisations.
    log_history : bool, optional, default : False
        A boolean that specifies whether to record the optimisation history.
    verbose : bool, optional, default : True
        A boolean that specifies the verbosity of the twinner.
    polish : bool, optional, default : True
        A boolean that specifies whether to run a deterministic tight-tolerance
        refinement of the ``n_polish`` best distinct multi-start results.
        Powell's default line-search tolerance stops early on the
        ill-conditioned (sloppy) physiological likelihood; the polish restarts
        Powell from each of them with tight tolerances until the objective
        stops improving, and the best polished result is kept. Off preserves
        the raw multi-start behaviour.
    n_polish : int, optional, default : 8
        The number of best multi-start results (with distinct objectives)
        that are polished; the best polished result is kept. The raw
        default-tolerance ranking is a poor guide to which start polishes
        into the best optimum, so polishing several starts is far more
        reliable than polishing only the best one. The polishes run in
        parallel when ``parallelize`` is ``True``. ``1`` polishes only the
        best start. Ignored when ``polish`` is ``False``.
    """
    if n_polish < 1:
        raise ValueError("n_polish must be at least 1")
    self.parallelize = parallelize
    self.n_jobs = -1 if n_jobs is None else n_jobs
    self.n_starts = n_starts
    self.verbose = verbose
    self.polish = polish
    self.n_polish = n_polish

    self.log_history = log_history
    if self.log_history:
        self.history = dict()
        self.history['theta'] = []
        self.history['log_prior'] = []
        self.history['log_likelihood'] = []
        self.history['log_posterior'] = []

twin(model, rbg_data, unknown_parameters_prior, correlations=None)

Runs the multi-start MAP estimation for a model.

Builds n_starts initial guesses by sampling the parameter priors, runs a Powell optimisation of the negative log-posterior (bounds enforced through the prior) from each start (sequentially or in parallel), polishes the n_polish best starts with distinct objectives (when polish is on, sequentially or in parallel), and selects the result with the lowest objective. Integer parameters are rounded and all parameters are clipped to their bounds before being returned.

Parameters:

Name	Type	Description	Default
model	object	A model instance being fit, implementing reset(theta0, g0, u0), step(u, t) and output(t).	required
rbg_data	object	A data object used to compute the likelihood.	required
unknown_parameters_prior	dict	A dictionary describing the priors and bounds for each parameter to estimate (keys prior, min, max and optional integer).	required
correlations	dict or None	Optional pairwise prior correlations specified as {(name_a, name_b): rho}. When given, the joint prior is a Gaussian copula over the named parameters (marginals unchanged); when None the parameters are treated as independent.	: None

Returns:

Type	Description
dict	A dictionary with keys fun (the best negative log-posterior value) and x (the best-fit parameter vector, rounded and clipped to the bounds).

Source code in src/py_replay_bg/twinner/twinner.py

def twin(self, model : Any, rbg_data, unknown_parameters_prior, correlations=None) -> dict:
    """Runs the multi-start MAP estimation for a model.

    Builds ``n_starts`` initial guesses by sampling the parameter priors,
    runs a Powell optimisation of the negative log-posterior (bounds enforced
    through the prior) from each start (sequentially or in parallel),
    polishes the ``n_polish`` best starts with distinct objectives (when
    ``polish`` is on, sequentially or in parallel), and selects the result
    with the lowest objective. Integer
    parameters are rounded and all parameters are clipped to their bounds
    before being returned.

    Parameters
    ----------
    model : object
        A model instance being fit, implementing ``reset(theta0, g0, u0)``,
        ``step(u, t)`` and ``output(t)``.
    rbg_data : object
        A data object used to compute the likelihood.
    unknown_parameters_prior : dict
        A dictionary describing the priors and bounds for each parameter to
        estimate (keys ``prior``, ``min``, ``max`` and optional ``integer``).
    correlations : dict or None, optional, default : None
        Optional pairwise prior correlations specified as
        ``{(name_a, name_b): rho}``. When given, the joint prior is a
        Gaussian copula over the named parameters (marginals unchanged); when
        ``None`` the parameters are treated as independent.

    Returns
    -------
    dict
        A dictionary with keys ``fun`` (the best negative log-posterior
        value) and ``x`` (the best-fit parameter vector, rounded and clipped
        to the bounds).
    """

    # Build the (validated) Gaussian-copula correlation structure once, up
    # front, so any error surfaces before the optimisation starts.
    correlation_structure = _build_correlation_structure(unknown_parameters_prior, correlations)

    # Set the worker arguments
    global _worker_args
    _worker_args = (self._neg_log_posterior, self._log_posterior, unknown_parameters_prior, model, rbg_data,
                    correlation_structure)

    # Build initial guesses for the parameters using their priors
    start_guesses = []
    for i in range(self.n_starts):
        start_guess = np.array([v['prior'].sample(v['min'], v['max'], i)
                                for v in unknown_parameters_prior.values()])
        start_guesses.append((i, start_guess))

    # Resolve the effective number of jobs (used both for running and for the header)
    n_jobs = multiprocessing.cpu_count() if self.n_jobs == -1 else self.n_jobs

    # Print a run header so the user can sanity-check the twinning setup
    if self.verbose:
        param_names = list(unknown_parameters_prior.keys())
        mode = f'parallel ({n_jobs} jobs)' if self.parallelize else 'sequential'
        print('--- Twinning (MAP estimation) ---')
        print(f'  starts          : {self.n_starts}')
        print(f'  mode            : {mode}')
        print(f'  polish          : ' + (f'top {self.n_polish} distinct starts' if self.polish else 'off'))
        print(f'  parameters ({len(param_names)})  : {", ".join(param_names)}')
        if correlation_structure is not None:
            pairs = ", ".join(f'{a}~{b}: {rho:g}'
                              for (a, b), rho in correlations.items())
            print(f'  correlations    : {pairs}')

    t0 = time.perf_counter()

    if self.parallelize:
        # Set the context
        ctx = multiprocessing.get_context('fork' if os.name != 'nt' else 'spawn')
        # Run the optimization in parallel
        with ctx.Pool(processes=n_jobs) as pool:
            raw = list(tqdm(
                pool.imap(_run_optimization, start_guesses),
                total=self.n_starts,
                desc='Twinning using MAP',
                disable=not self.verbose
            ))
        results = [r for r, _ in raw]
        if self.log_history:
            for _, h in raw:
                if h is not None:
                    for k in self.history:
                        self.history[k].extend(h[k])
    else:
        # Run the optimization sequentially; history accumulates in-process via _neg_log_posterior
        raw = [_run_optimization(a)
               for a in tqdm(start_guesses, desc='Twinning', disable=not self.verbose)]
        results = [r for r, _ in raw]

    # Refine the best starts with a deterministic tight-tolerance polish. The
    # multi-start Powell runs at scipy's default tolerance, which halts early
    # on the sloppy (ill-conditioned) physiological likelihood ridge; a
    # tight-tolerance restart loop walks the rest of the way down to the
    # optimum. The raw ranking is a poor guide to which start polishes into
    # the best optimum, so the n_polish best distinct starts are polished and
    # the best polished result is kept. No RNG is involved, so twinning stays
    # reproducible.
    candidates = _top_distinct(results, self.n_polish) if self.polish else []
    if candidates:
        # Copies, so the raw results (used in the summary) stay untouched
        polish_args = [copy.copy(results[i]) for i in candidates]
        if self.parallelize and len(candidates) > 1:
            ctx = multiprocessing.get_context('fork' if os.name != 'nt' else 'spawn')
            with ctx.Pool(processes=min(n_jobs, len(candidates))) as pool:
                raw = list(tqdm(
                    pool.imap(_run_polish, polish_args),
                    total=len(candidates),
                    desc='Polishing',
                    disable=not self.verbose
                ))
            polished = [r for r, _ in raw]
            if self.log_history:
                for _, h in raw:
                    if h is not None:
                        for k in self.history:
                            self.history[k].extend(h[k])
        else:
            polished = [self._polish(r, model, rbg_data, unknown_parameters_prior, correlation_structure)
                        for r in tqdm(polish_args, desc='Polishing', disable=not self.verbose)]
        # min keeps the first on ties, i.e. the better raw rank
        best_rank = min(range(len(polished)), key=lambda j: polished[j].fun)
        best = polished[best_rank]
    else:
        best_rank = None
        best = min(results, key=lambda r: r.fun)

    elapsed = time.perf_counter() - t0

    # Print a convergence summary across all starts
    if self.verbose:
        funs = np.array([r.fun for r in results], dtype=float)
        finite = funs[np.isfinite(funs)]
        n_failed = len(funs) - len(finite)
        print('--- Twinning summary ---')
        print(f'  elapsed         : {elapsed:.1f} s')
        print(f'  failed starts   : {n_failed} / {len(funs)} (non-finite objective)')
        if len(finite):
            print(f'  objective (neg log-posterior, lower is better)')
            print(f'    raw starts    : min {finite.min():.3f} / med {np.median(finite):.3f} / '
                  f'max {finite.max():.3f}')
            if candidates:
                pol = np.array([r.fun for r in polished], dtype=float)
                winner = candidates[best_rank]
                print(f'    polished ({len(candidates)})  : best {pol.min():.3f} / worst {pol.max():.3f}')
                print(f'    final         : {best.fun:.3f}  (start {winner}, raw {results[winner].fun:.3f}, '
                      f'{_ordinal(best_rank + 1)} best distinct raw start)')
            else:
                print(f'    final         : {best.fun:.3f}  (not polished)')
        else:
            print('  WARNING: no start produced a finite objective '
                  '(check priors/bounds and the model).')

    # Round integer parameters and clip to bounds (Powell has no native bound support)
    clipped_x = []
    for x, (k, v) in zip(best.x, unknown_parameters_prior.items()):
        if v.get('integer', False):
            x = float(round(x))
        clipped_x.append(np.clip(x, v['min'], v['max']))
    clipped_x = np.array(clipped_x)

    # Print the best-fit parameters, flagging any that sit at their bounds
    if self.verbose:
        print('--- Best-fit parameters ---')
        for x, (k, v) in zip(clipped_x, unknown_parameters_prior.items()):
            lo, hi = v['min'], v['max']
            tol = 1e-6 * max(1.0, abs(hi - lo))
            flags = []
            if x <= lo + tol:
                flags.append('at min bound')
            elif x >= hi - tol:
                flags.append('at max bound')
            if v.get('integer', False):
                flags.append('integer')
            suffix = f'   [{", ".join(flags)}]' if flags else ''
            print(f'  {k:<14} = {x:.4g}   (bounds: {lo:g}..{hi:g}){suffix}')

    # Return the best result
    ret = dict()
    ret['fun'] = best.fun
    ret['x'] = clipped_x
    return ret

_polish(best, model, rbg_data, unknown_parameters_prior, correlation_structure=None, max_rounds=5, tol=1e-08, rel_improve=1e-07, evals_per_dim=1000, max_budget=8000)

Refines the best multi-start result with a tight-tolerance Powell loop.

Powell resets its search-direction set on each fresh call, so restarting it from the incumbent with tight tolerances lets it keep descending a narrow, ill-conditioned valley that a single default-tolerance run stops short of. The loop repeats until a round fails to improve the objective by a meaningful relative amount, so it stops promptly once the optimum is reached. The starting point is fixed and no randomness is involved, so the result is deterministic.

Cost is bounded on purpose: each restart's function-evaluation budget is min(evals_per_dim * n_params, max_budget) and the number of restarts is capped, with an early stop once a round stops improving. The absolute max_budget cap keeps a single restart affordable on high-dimensional models (e.g. the 25-parameter extended model, where an uncapped tight Powell runs for many minutes), while the max_rounds restarts let the cheaper low-dimensional models still converge fully. The bounds depend only on fixed problem sizes, so the result stays deterministic.

Parameters:

Name	Type	Description	Default
best	OptimizeResult	The best multi-start result; its x/fun seed the refinement.	required
model	object	The model instance being fit (mutated in place by each evaluation).	required
rbg_data	object	The data object used to compute the likelihood.	required
unknown_parameters_prior	dict	The priors/bounds spec (fixes parameter order and integer handling).	required
correlation_structure	dict or None	The precomputed Gaussian-copula structure, or None.	: None
max_rounds	int	Maximum number of tight-tolerance Powell restarts (the early stop often triggers before this).	: 5
tol	float	The xtol/ftol passed to each Powell restart. Far tighter than scipy's 1e-4 default (which caused the under-convergence).	: 1e-8
rel_improve	float	Minimum relative objective improvement (fun - r.fun)/max(1, |fun|) a round must yield to keep going; below it the loop stops.	: 1e-7
evals_per_dim	int	Per-restart function-evaluation budget per parameter (before the cap).	: 1000
max_budget	int	Absolute cap on a single restart's function-evaluation budget; bounds the per-restart cost on high-dimensional models.	: 8000

Returns:

Type	Description
OptimizeResult	best with x/fun updated to the refined optimum (unchanged when no restart improves on the incumbent).

Source code in src/py_replay_bg/twinner/twinner.py

def _polish(self, best, model, rbg_data, unknown_parameters_prior,
            correlation_structure=None, max_rounds: int = 5, tol: float = 1e-8,
            rel_improve: float = 1e-7, evals_per_dim: int = 1000,
            max_budget: int = 8000):
    """Refines the best multi-start result with a tight-tolerance Powell loop.

    Powell resets its search-direction set on each fresh call, so restarting
    it from the incumbent with tight tolerances lets it keep descending a
    narrow, ill-conditioned valley that a single default-tolerance run stops
    short of. The loop repeats until a round fails to improve the objective by
    a meaningful *relative* amount, so it stops promptly once the optimum is
    reached. The starting point is fixed and no randomness is involved, so the
    result is deterministic.

    Cost is bounded on purpose: each restart's function-evaluation budget is
    ``min(evals_per_dim * n_params, max_budget)`` and the number of restarts is
    capped, with an early stop once a round stops improving. The absolute
    ``max_budget`` cap keeps a single restart affordable on high-dimensional
    models (e.g. the 25-parameter extended model, where an uncapped tight
    Powell runs for many minutes), while the ``max_rounds`` restarts let the
    cheaper low-dimensional models still converge fully. The bounds depend only
    on fixed problem sizes, so the result stays deterministic.

    Parameters
    ----------
    best : scipy.optimize.OptimizeResult
        The best multi-start result; its ``x``/``fun`` seed the refinement.
    model : object
        The model instance being fit (mutated in place by each evaluation).
    rbg_data : object
        The data object used to compute the likelihood.
    unknown_parameters_prior : dict
        The priors/bounds spec (fixes parameter order and integer handling).
    correlation_structure : dict or None, optional, default : None
        The precomputed Gaussian-copula structure, or ``None``.
    max_rounds : int, optional, default : 5
        Maximum number of tight-tolerance Powell restarts (the early stop
        often triggers before this).
    tol : float, optional, default : 1e-8
        The ``xtol``/``ftol`` passed to each Powell restart. Far tighter than
        scipy's 1e-4 default (which caused the under-convergence).
    rel_improve : float, optional, default : 1e-7
        Minimum relative objective improvement ``(fun - r.fun)/max(1, |fun|)``
        a round must yield to keep going; below it the loop stops.
    evals_per_dim : int, optional, default : 1000
        Per-restart function-evaluation budget per parameter (before the cap).
    max_budget : int, optional, default : 8000
        Absolute cap on a single restart's function-evaluation budget; bounds
        the per-restart cost on high-dimensional models.

    Returns
    -------
    scipy.optimize.OptimizeResult
        ``best`` with ``x``/``fun`` updated to the refined optimum (unchanged
        when no restart improves on the incumbent).
    """
    x = np.asarray(best.x, dtype=float)
    fun = best.fun
    budget = min(evals_per_dim * len(x), max_budget)
    for _ in range(max_rounds):
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            r = minimize(self._neg_log_posterior, x, method='Powell',
                         args=(model, rbg_data, unknown_parameters_prior,
                               correlation_structure),
                         options={'maxiter': budget, 'maxfev': budget,
                                  'xtol': tol, 'ftol': tol, 'disp': False})
        if not np.isfinite(r.fun) or r.fun >= fun:
            break
        gain = (fun - r.fun) / max(1.0, abs(fun))
        x, fun = r.x, r.fun
        if gain < rel_improve:
            break
    best.x = x
    best.fun = fun
    return best

_log_prior(model, unknown_parameters_prior, correlation_structure=None)

Computes the log prior probability of the model parameters.

The base term is the sum of marginal log-densities. When a Gaussian-copula correlation_structure is supplied, a copula correction term is added so that the named parameters are jointly correlated while keeping their marginals unchanged.

Parameters:

Name	Type	Description	Default
model	object	A model instance containing the current parameter values.	required
unknown_parameters_prior	dict	A dictionary describing the priors for each parameter. Each prior must provide a log_pdf(value) method.	required
correlation_structure	dict or None	Precomputed Gaussian-copula structure (see :func:_build_correlation_structure). When None the parameters are treated as independent.	: None

Returns:

Type	Description
float	The log prior probability, or -inf when any parameter is outside its valid range.

Source code in src/py_replay_bg/twinner/twinner.py

def _log_prior(self, model, unknown_parameters_prior, correlation_structure=None):
    """Computes the log prior probability of the model parameters.

    The base term is the sum of marginal log-densities. When a Gaussian-copula
    ``correlation_structure`` is supplied, a copula correction term is added so
    that the named parameters are jointly correlated while keeping their
    marginals unchanged.

    Parameters
    ----------
    model : object
        A model instance containing the current parameter values.
    unknown_parameters_prior : dict
        A dictionary describing the priors for each parameter. Each prior
        must provide a ``log_pdf(value)`` method.
    correlation_structure : dict or None, optional, default : None
        Precomputed Gaussian-copula structure (see
        :func:`_build_correlation_structure`). When ``None`` the parameters
        are treated as independent.

    Returns
    -------
    float
        The log prior probability, or ``-inf`` when any parameter is outside
        its valid range.
    """
    # Iterate over the parameters and compute the marginal log prior
    lp = 0.0
    for up, v in unknown_parameters_prior.items():
        parameter_value = getattr(model, up)
        a = v["min"]
        b = v["max"]
        # If the parameter is outside the valid range, return -inf
        if parameter_value <= a or parameter_value >= b:
            return -np.inf
        # Otherwise, add the marginal log prior contribution (corrected for the truncation)
        log_density = v["prior"].log_pdf(parameter_value)
        z = v["prior"].cdf(b) - v["prior"].cdf(a)
        if log_density == -np.inf or z <= 0:
            return -np.inf
        lp += log_density - np.log(z)

    # Add the Gaussian-copula correction for correlated parameters (if any)
    if correlation_structure is not None:
        names = correlation_structure['names']
        # Map each correlated parameter to a standard-normal score z = Phi^-1(F(x))
        z_scores = np.empty(len(names))
        for i, name in enumerate(names):
            v = unknown_parameters_prior[name]
            x = getattr(model, name)
            fa = v["prior"].cdf(v["min"])
            fb = v["prior"].cdf(v["max"])
            fx = v["prior"].cdf(x)
            u = (fx - fa) / (fb - fa)
            u = np.clip(u, 1e-12, 1.0 - 1e-12)
            z_scores[i] = ndtri(u)
        # log copula density: -0.5*logdet(R) - 0.5*z^T (R^-1 - I) z
        lp += -0.5 * correlation_structure['logdet_R'] \
              - 0.5 * float(z_scores @ correlation_structure['R_inv_minus_I'] @ z_scores)

    # Return the (possibly correlation-corrected) log prior
    return lp

_log_likelihood(model, rbg_data)

Computes the log likelihood of the observed data under the model.

The model is simulated forward using the input sequence in rbg_data. The predicted output is then compared against the observed glucose measurements using a Gaussian error model with a constant coefficient of variation.

Parameters:

Name	Type	Description	Default
model	object	A model instance to simulate.	required
rbg_data	object	A data object containing inputs, output timestamps, and observed glucose values.	required

Returns:

Type	Description
float	The log likelihood value for the simulated trajectory.

Source code in src/py_replay_bg/twinner/twinner.py

def _log_likelihood(self, model, rbg_data):
    """Computes the log likelihood of the observed data under the model.

    The model is simulated forward using the input sequence in ``rbg_data``.
    The predicted output is then compared against the observed glucose
    measurements using a Gaussian error model with a constant coefficient of
    variation.

    Parameters
    ----------
    model : object
        A model instance to simulate.
    rbg_data : object
        A data object containing inputs, output timestamps, and observed
        glucose values.

    Returns
    -------
    float
        The log likelihood value for the simulated trajectory.
    """
    # Simulate the model forward and get the output
    out = np.zeros(rbg_data.tsteps, )
    out[0] = model.output(0)
    for k in np.arange(1,out.shape[0]):
        model.step(rbg_data.u[k], k)
        out[k] = model.output(k)
    # Subsample the output to match the sampling rate of the data
    out = out[0::rbg_data.yts]

    # TODO: enable the choice of multiple ll shapes
    # Two-term noise model: σ² = σ_add² + (cv·|ŷ|)²
    #sigma_add = 5.0          # mg/dL, literature default
    #cv        = 0.05         # 5%
    #sdn = np.sqrt(sigma_add**2 + (cv * np.abs(out[rbg_data.y_idxs]))**2)
    #residuals = out[rbg_data.y_idxs] - rbg_data.y[rbg_data.y_idxs]

    # Calculate the log-likelihood with a Gaussian error model (constant coefficient of variation)
    cv = 0.05  # constant coefficient of variation (5%) #TODO: make this a parameter
    residuals = out[rbg_data.y_idxs] - rbg_data.y[rbg_data.y_idxs]
    sdn = cv * np.abs(out[rbg_data.y_idxs])
    subsampling = 1 #TODO: decide where to evaluate ll (now every 30 minutes but maybe we can do something smarter)
    sdn = sdn[0::subsampling]
    residuals = residuals[0::subsampling]
    n = len(residuals)
    return -0.5 * np.sum((residuals / sdn) ** 2) - np.sum(np.log(sdn)) - 0.5 * n * np.log(2 * np.pi)

_neg_log_posterior(theta, model, rbg_data, unknown_parameters_prior, correlation_structure=None)

Returns the negative log posterior for optimization.

Parameters:

Name	Type	Description	Default
theta	ndarray	A parameter vector in the natural (constrained) parameter space, ordered to match unknown_parameters_prior.	required
model	object	A model instance being fit.	required
rbg_data	object	A data object used to compute the likelihood.	required
unknown_parameters_prior	dict	A dictionary describing priors and bounds for the parameters.	required

Returns:

Type	Description
float	The negative log posterior value.

Source code in src/py_replay_bg/twinner/twinner.py

def _neg_log_posterior(self, theta, model, rbg_data, unknown_parameters_prior, correlation_structure=None):
    """Returns the negative log posterior for optimization.

    Parameters
    ----------
    theta : numpy.ndarray
        A parameter vector in the natural (constrained) parameter space,
        ordered to match ``unknown_parameters_prior``.
    model : object
        A model instance being fit.
    rbg_data : object
        A data object used to compute the likelihood.
    unknown_parameters_prior : dict
        A dictionary describing priors and bounds for the parameters.

    Returns
    -------
    float
        The negative log posterior value.
    """
    # Just return the negative log-posterior
    log_prior, log_likelihood, log_post = self._log_posterior(theta, model, rbg_data, unknown_parameters_prior,
                                                              correlation_structure)

    # log the history
    if self.log_history:
        self.history['theta'].append(theta)
        self.history['log_prior'].append(log_prior)
        self.history['log_likelihood'].append(log_likelihood)
        self.history['log_posterior'].append(log_post)

    return -log_post

_log_posterior(theta, model, rbg_data, unknown_parameters_prior, correlation_structure=None)

Compute log-prior, log-likelihood, and log-posterior in one model pass.

The parameter vector is used directly in the natural (constrained) parameter space: each entry is written into the model's typed parameter dict (integer parameters are rounded) before model.reset. Every reset also re-seeds the cold-start glucose and the t=0 input row from rbg_data (see Model.reset), since reset clears them each call. Bounds are enforced by the prior, which returns -inf outside [min, max].

Parameters:

Name	Type	Description	Default
theta	ndarray	A parameter vector in the natural (constrained) parameter space, ordered to match unknown_parameters_prior.	required
model	object	A model instance being updated with candidate parameters.	required
rbg_data	object	A data object used to compute the likelihood.	required
unknown_parameters_prior	dict	A dictionary describing priors and bounds for the parameters.	required

Returns:

Type	Description
tuple	The triple (log_prior, log_likelihood, log_posterior). Any component that is invalid is returned as -np.inf.

Source code in src/py_replay_bg/twinner/twinner.py

def _log_posterior(self, theta, model, rbg_data, unknown_parameters_prior, correlation_structure=None):
    """Compute log-prior, log-likelihood, and log-posterior in one model pass.

    The parameter vector is used directly in the natural (constrained)
    parameter space: each entry is written into the model's typed parameter
    dict (integer parameters are rounded) before ``model.reset``. Every reset
    also re-seeds the cold-start glucose and the t=0 input row from
    ``rbg_data`` (see ``Model.reset``), since ``reset`` clears them each call.
    Bounds are enforced by the prior, which returns ``-inf`` outside
    ``[min, max]``.

    Parameters
    ----------
    theta : numpy.ndarray
        A parameter vector in the natural (constrained) parameter space,
        ordered to match ``unknown_parameters_prior``.
    model : object
        A model instance being updated with candidate parameters.
    rbg_data : object
        A data object used to compute the likelihood.
    unknown_parameters_prior : dict
        A dictionary describing priors and bounds for the parameters.

    Returns
    -------
    tuple
        The triple ``(log_prior, log_likelihood, log_posterior)``. Any
        component that is invalid is returned as ``-np.inf``.
    """

    # Create the theta input dictionary (note: order is maintained by construction)
    theta_dict = Dict.empty(key_type=types.unicode_type, value_type=float64)
    for i, k in enumerate(unknown_parameters_prior.keys()):
        val = theta[i]
        #val = np.clip(theta[i], unknown_parameters_prior[k]['min'], unknown_parameters_prior[k]['max'])
        is_int =  unknown_parameters_prior[k].get('integer', False)
        if is_int:
            val = int(round(val))
        theta_dict[k] = val

    # Reset the model with the new parameters, re-seeding the cold-start
    # glucose and the t=0 input row from data every time (reset() clears
    # both, so they can't just be seeded once outside this loop).
    g0 = float(rbg_data.y[rbg_data.y_idxs[0]]) if len(rbg_data.y_idxs) > 0 else np.nan
    model.reset(theta_dict, g0, rbg_data.u[0])

    # Calculate log-prior
    lp = self._log_prior(model, unknown_parameters_prior, correlation_structure)
    if lp == -np.inf or np.isnan(lp):
        return -np.inf, -np.inf, -np.inf
    # Calculate log-likelihood
    ll = self._log_likelihood(model, rbg_data)
    if ll == -np.inf or np.isnan(ll):
        return lp, -np.inf, -np.inf

    # Return log-posterior
    return lp, ll, lp + ll