PyBNF algorithms (pybnf.algorithms)

The pybnf.algorithms package holds PyBNF’s fitting and sampling algorithms. Every fit type is re-exported from the package facade (pybnf.algorithms.<Name>), but the implementations live in the submodules documented below: shared execution primitives in core and base, gradient-free and gradient-based optimizers under optimizers, Bayesian samplers under samplers, and the model-checking utility in model_check.

Core execution

Core execution primitives shared by every algorithm.

Holds the Result / FailedSimulation containers, the Job runner, the JobGroup smoothing/model-parallel folding family, run_job, result_from_completed, and the dask as_completed import.

This module is the patched seam (ADR-0001). Tests and the slow integration_harness patch run_job / as_completed / Job here, and the run loop resolves them as core.<name> so that “patch where it is defined” stays honest as the algorithm classes move out of the package facade in later commits.

Import direction is one-way: this module imports only leaf modules (data / pset / printing) and never anything from the algorithms package facade or any algorithm leaf.

class pybnf.algorithms.core.FailedSimulation(paramset, name, fail_type, einfo=(None, None, None))

normalize(settings)

Normalizes the Data object in this result, according to settings :param settings: Config value for ‘normalization’: a string representing the normalization type, a dict mapping exp files to normalization type, or None :return:

postprocess_data(settings)

Postprocess the Data objects in this result with a user-defined Python script :param settings: A dict that maps a tuple (model, suffix) to a Python filename to load. That file is expected to contain the definition for the function postprocess(data), which takes a Data object and returns a processed data object :return: None

class pybnf.algorithms.core.HybridJobGroup(job_id, replica_subjob_ids)

A JobGroup to handle combined smoothing and model-level parallelism.

average_results()

Merge model partitions within each smoothing replica, then average the complete replicas.

job_finished(res)

Called when one job in this group has finished. :param res: Result object for the completed job :return: Boolean, whether everything in this job group has finished

class pybnf.algorithms.core.Job(models, params, job_id, output_dir, timeout, calc_future, norm_settings, postproc_settings, delete_folder=False, replicate_index=0, stochastic_seed_policy='auto', model_slice=None)

Container for information necessary to perform a single evaluation in the fitting algorithm

run_simulation(debug=False, failed_logs_dir='')

Runs the simulation and reads in the result

class pybnf.algorithms.core.JobGroup(job_id, subjob_ids)

Represents a group of jobs that are identical replicates to be averaged together for smoothing

average_results()

To be called after all results are in for this group. Averages the results and returns a new Result object containing the averages

Returns:

New Result object with the job_id of this JobGroup and the averaged Data as the simdata

job_finished(res)

Called when one job in this group has finished :param res: Result object for the completed job :return: Boolean, whether everything in this job group has finished

class pybnf.algorithms.core.MultimodelJobGroup(job_id, subjob_ids)

A JobGroup to handle model-level parallelism

average_results()

To be called after all results are in for this group. Combines all results from the submodels into a single Result object :return:

class pybnf.algorithms.core.Result(paramset, simdata, name)

Container for the results of a single evaluation in the fitting algorithm

add_result(other)

Add simulation data of other models from another Result object into this Result object :param other: The other Result object :return:

normalize(settings)

Normalizes the Data object in this result, according to settings :param settings: Config value for ‘normalization’: a string representing the normalization type, a dict mapping exp files to normalization type, or None :return:

postprocess_data(settings)

Postprocess the Data objects in this result with a user-defined Python script :param settings: A dict that maps a tuple (model, suffix) to a Python filename to load. That file is expected to contain the definition for the function postprocess(data), which takes a Data object and returns a processed data object :return: None

pybnf.algorithms.core.result_from_completed(future, result, params, job_id)

Translate one completed item from as_completed(futures, with_results=True, raise_errors=False) into a Result.

With those flags, dask hands back per future:

• the job’s return value (a Result / FailedSimulation) for a finished job,

• the worker’s (type, exception, traceback) tuple for an errored job, and

• a CancelledError for a cancelled job.

A failed job becomes a FailedSimulation so the run continues, except a user-targeted PybnfError is re-raised (it would fail every job, so failing fast is better than burning the whole fit). A CancelledError is returned unchanged for the caller to treat as fatal, and anything unrecognized is treated as a failed simulation rather than crashing the run.

Future.status is part of dask’s public Future API, so this relies on no dask internals. See lanl/PyBNF#388 for the history (the previous approach subclassed as_completed and overrode a private method that dask later renamed).

pybnf.algorithms.core.run_job(j, debug=False, failed_logs_dir='', models=None, calc=None)

Runs the Job j. This function is passed to Dask instead of j.run_simulation because if you pass j.run_simulation, Dask leaks memory associated with j.

models / calc are the worker-resolved model_list and ObjectiveCalculator that Algorithm.run scattered once per fit and hands to client.submit as direct kwargs, so dask substitutes the Futures with their concrete values here on the worker (see _ResolvedFuture for why a Future stashed as a Job attribute cannot be used instead). When provided they are rebound onto the Job; when omitted – the in-process paths that call run_job directly (_rerun_best_fit_to_save_data, ModelCheck, the folder-free test harnesses) – the Job’s own models / calc_future attributes are used as-is.

The Algorithm base class and the module-level helpers it relies on.

Extracted from algorithms.py (M1 Step 2). The run loop resolves the three monkeypatched names as core.run_job / core.as_completed / core.Job (ADR-0001), so it imports the core module object rather than binding those names locally; the never-patched data classes are imported by name.

class pybnf.algorithms.base.Algorithm(config)

Base class for every PyBNF fit type (“method”); defines the run-loop contract.

The contract (ADR-0007). A method plugs into the framework by subclassing Algorithm and implementing exactly two abstract methods:

• start_run() () -> list[PSet] — the initial parameter sets to evaluate.

• got_result() (Result) -> list[PSet] | 'STOP' — given one completed result, the next parameter sets to evaluate, or the string 'STOP' to end the run.

plus registering itself with a config schema via @register_fit_type(..., schema=...) (ADR-0002, ADR-0005). The .name of any returned PSet, if set, MUST be unique across the run (it determines the simulation folder name; uniqueness is not checked elsewhere).

Everything else the framework needs has an overridable default (reset, add_iterations, cleanup, get_backup_every, should_pickle). The run loop itself — job submission, as_completed draining, backup cadence, best-fit save, sim-dir teardown — lives in run() and is shared by every method, not replaceable: a method customizes what to propose, never the outer loop.

add_iterations(n)

Adds n additional iterations to the algorithm. May be overridden in subclasses that don’t use self.max_iterations to track the iteration count

add_to_trajectory(res)

Adds the information from a Result to the Trajectory instance

backup(pending_psets=())

Create a backup of this algorithm object that can be reloaded later to resume the run

Parameters:

pending_psets – Iterable of PSets that are currently submitted as jobs, and will need to get re-submitted when resuming the algorithm

Returns:

cleanup()

Called before the program exits due to an exception. :return:

get_backup_every()

Returns a number telling after how many individual simulation returns should we back up the algorithm. Makes a good guess, but could be overridden in a subclass

abstractmethod got_result(res)

Called by the scheduler when a simulation is completed, with the pset that was run, and the resulting simulation data

Parameters:

res (Result) – result from the completed simulation

Returns:

List of PSet(s) to be run next or ‘STOP’ string.

make_job(params)

Creates a new Job using the specified params, and additional specifications that are already saved in the Algorithm object If smoothing or model-level parallelism is turned on, makes grouped subjobs.

Parameters:

params (PSet)

Returns:

list of Jobs

output_results(name='', no_move=False)

Tells the Trajectory to output a log file now with the current best fits.

This should be called periodically by each Algorithm subclass, and is called by the Algorithm class at the end of the simulation. :return: :param name: Custom string to add to the saved filename. If omitted, we just use a running counter of the number of times we’ve outputted. :param no_move: If True, overrides the config setting delete_old_files=2, and does not move the result to overwrite sorted_params.txt :type name: str

random_latin_hypercube_psets(n)

Generates n random PSets with a latin hypercube distribution. Variables with bounded initialization distributions follow the Latin hypercube; the others are randomized independently from their initializer.

Parameters:

n – Number of psets to generate

Returns:

random_pset()

Generates a random PSet based on the distributions and bounds for each parameter specified in the configuration

Returns:

reset(bootstrap)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

run(client, resume=None, debug=False)

Main loop for executing the algorithm

static should_pickle(k)

Checks to see if key ‘k’ should be included in pickling. Currently allows all entries in instance dictionary except for ‘trajectory’

Parameters:

k

Returns:

spawn_chain_rngs(n)

Return n independent Generators, one per parallel chain/replica.

Spawning is deterministic in (seed, n), so chain i always draws from the same stream regardless of the order in which dask returns its results. This is what makes the parallel samplers reproducible under nondeterministic scheduling – a single shared stream would interleave by completion order and so would not reproduce run to run.

abstractmethod start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

pybnf.algorithms.base.exp10(n)

Raise 10 to the power of a possibly user-defined value, and raise a helpful error if it overflows :param n: A float :return: 10.** n

pybnf.algorithms.base.latin_hypercube(nsamples, ndims, rng)

Latin hypercube sampling.

Returns a nsamples by ndims array, with entries in the range [0,1] You’ll have to rescale them to your actual param ranges.

rng is the caller’s np.random.Generator (the algorithm’s root rng).

Optimizers

Shared bases

Shared scaffolding for the start-point local optimizers (Powell + CMA-ES, #403).

Both Powell (conjugate-direction) and CMA-ES are derivative-free, black-box local optimizers that begin from a single point and search in sampling space u – log10 for log-scaled parameters, linear otherwise. That is the same space the prior and proposal arithmetic already operate in (FreeParameter._scale, ADR-0003/0010), so log parameters are optimized geometrically (a multiplicative step is an additive u step) exactly as Simplex does its log-space arithmetic.

StartPointOptimizer factors out the two pieces of plumbing they share:

• start-point resolution – the injected refiner start point (set by pybnf._refine_best_fit under START_POINT_KEY) when refining; the single-value var / logvar specs of a standalone point-start fit (the same start point Simplex parses, ADR-0015); or, for a start_from_box optimizer (CMA-ES) given bounded uniform_var / loguniform_var priors, the box center in sampling space u – its global-start mode (#404, ADR-0017);

• the u <-> PSet conversion, which maps each coordinate back to a stored value and reflects it into the box via FreeParameter.set_value() (a no-op for the unbounded var / logvar of a point-start fit; active when refining or globally searching a bounded fit’s parameters).

Simplex predates this and keeps its own byte-identical start-point parsing; the two new optimizers are the ≥2-member event (ADR-0009) that earns the shared base. These methods plug into the run loop through start_run / got_result only (ADR-0007); no method overrides run().

class pybnf.algorithms.optimizers.local_base.StartPointOptimizer(config)

Base for the start-point local optimizers. Subclasses implement start_run / got_result and set START_POINT_KEY.

START_POINT_KEY = None

The internal config key the refiner start point is injected under (mirrors SimplexAlgorithm’s 'simplex_start_point'). Set by each subclass; pybnf._refine_best_fit writes the best fit here so refinement starts from it instead of parsing the (refiner-irrelevant) variable specs.

Shared scaffolding for the gradient-based local optimizers (TRF/LM + L-BFGS-B, #386).

The metaheuristic fit types (de, pso, ss, cmaes, …) only ever ask each evaluated PSet for its scalar objective value. A gradient optimizer instead consumes the residual vector + residual-Jacobian (TRF / Levenberg–Marquardt) or the scalar gradient (L-BFGS-B) that #385 assembles from bngsim’s forward output-sensitivity tensor. GradientOptimizer factors out everything a new gradient method needs so a leaf (trf.py, lbfgs.py) implements only its step math – mirroring how StartPointOptimizer factors the start-point / u ↔ PSet plumbing out of Powell and CMA-ES.

What this base provides

• The edition + capability gate (_gate_gradient_supported()). Gradient fitting consumes the edition-2 surface (bind-by-id routing ADR-0034, the noise_model / measurement layer) and bngsim’s forward sensitivities, so the fit is refused – with a clear message pointing at a metaheuristic fit_type – on a legacy (edition < 2) config, a non-bngsim model, or a bngsim build without the output_sensitivities feature. Never a silent finite-difference fallback.

• The gradient path activation (_setup_gradient_path()). Builds each experiment’s ExperimentRouting once (it depends only on model structure, conditions, and free-parameter ids – never on the parameter values, #449) and apply_routings the union request onto every model, so each simulated Data carries its sensitivity tensor. Run before the model scatter (from start_run); the request rides the pickle to the workers (BngsimModel.__getstate__ keeps it, rebuilding only the engine).

• Master-side scoring (requires_master_scoring = True). The worker scoring path nulls res.simdata after scoring (#385/#388), which would discard the sensitivity tensors; the flag makes Algorithm.run keep scoring on the master so every Result returns with its full simdata for gradient_at().

• The per-evaluation assembly (gradient_at()). Aligns a Result’s simdata with exp_data and the prebuilt routings, and returns the assembled GradientResult (objective gradient + residual Jacobian, in sampling space u), folding in any constraint-penalty gradient.

• The ``u``-space box (_u_bounds()). The finite reflecting box for bounded (uniform_var / loguniform_var) priors, ±inf for an unbounded point start – the bounds the leaf’s step projects/reflects into.

The search runs in sampling space u (StartPointOptimizer), and #385 already delivers the gradient transformed into u once (ADR-0029), so a leaf never re-transforms. Leaves own their start_run / got_result state machine and must be picklable for backup/resume, exactly like Powell and CMA-ES (ADR-0007).

pybnf.algorithms.optimizers.gradient_base.DONE = <object object>

Sentinel a GradientRunner’s got returns when that start has converged or otherwise terminated (no further evaluation to propose). It is only ever returned synchronously into GradientOptimizer.got_result() and never stored in pickled state, so a plain module-level object (identity-checked) is enough.

class pybnf.algorithms.optimizers.gradient_base.GradientOptimizer(config, refine=False)

Base for the gradient-based local optimizers (#386).

A leaf subclass supplies only its per-start step math as a GradientRunner (Levenberg–Marquardt for trf, L-BFGS-B for lbfgs) via _make_runner(); this base owns the whole start_run / got_result orchestration – seeding the runners, the u``<->PSet plumbing, :meth:`gradient_at` assembly, name-tagging / routing, reporting, and the multi-start ``STOP coordination. The leaf must set START_POINT_KEY like any StartPointOptimizer, plus _method_label and _start_banner() for its messages.

Local multi-start (#386). A box-start gradient fit runs N independent starts concurrently (N reuses population_size) – start 0 from the box center (preserving the deterministic single-start behavior), the rest from Latin-hypercube samples across the prior box – and keeps the global best. The async run loop already supports this without change: start_run returns a list of N initial PSets, each got_result advances just the start that owns the returned Result (routed by PSet name through pending) and returns that start’s next PSet (or [] once it converges), and 'STOP' fires only when the last start finishes. Every evaluated PSet across all starts lands in the trajectory (add_to_trajectory runs before got_result), so trajectory.best_fit() is the global best for free – each runner only tracks its own best for its own convergence test.

add_iterations(n)

Extend every start’s per-start iteration budget by n (the -r resume path). Runners already exist when this is called on a resumed run (they ride the backup pickle), so bump each alongside the template budget.

got_result(res)

Route a completed Result to the start that owns it, advance just that start’s runner, and return its next PSet – [] once it converges (other starts keep going), or 'STOP' only when the last live start finishes.

gradient_at(res)

Assemble the objective gradient + residual-Jacobian at res’s point.

res is a master-scored Result, so res.simdata carries each experiment’s forward-sensitivity tensor. Aligns it with exp_data over the scored (model, suffix) pairs (the same intersection the objective scores), attaches each one’s prebuilt routing, and returns the assembled GradientResult – residual / Jacobian / scalar gradient in sampling space u (#385 transformed it once; the optimizer never re-transforms). Any constraint-penalty gradient is added to the scalar gradient (and clears least_squares_exact, since a penalty is not a sum of squares).

The free-parameter list (column order + current values) is read straight off the evaluated PSet, so the d theta/d u scale factors are taken at the point actually simulated. Converts a GradientNotSupported (an objective the assembly cannot differentiate) into a clear, fail-fast PybnfError pointing at a metaheuristic fit_type.

requires_master_scoring = True

Keep objective scoring on the master so every Result returns with its simdata (the sensitivity tensors the gradient assembly reads); see Algorithm.run. Without this the worker path nulls res.simdata.

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Seed the n_starts runners and return their initial PSets (one evaluation each), so all starts run concurrently from the first scheduler batch.

class pybnf.algorithms.optimizers.gradient_base.GradientRunner(u0, lower, upper, max_iterations)

Headless, picklable per-start step machine in sampling space u (#386).

A gradient leaf’s step math, factored out of the optimizer so a single fit can run N of them concurrently – local multi-start, the diversity a purely local gradient method otherwise lacks (it only ever descends into the one basin its start lands in). A runner owns one start’s entire mutable state – the iterate, the curvature / trust-region model, the reflecting box, the tunables – and is pure numpy: it knows nothing about PSets, the objective, routing, backup, or dask. GradientOptimizer drives it:

u0 = runner.start()                       # the first point to evaluate
nxt = runner.got(u_point, score, grad)    # consume one completed evaluation

where u_point is the realized (box-projected) u-vector of the completed evaluation, score its objective value, and grad the assembled GradientResult at it; got returns the next u to evaluate, or the DONE sentinel. Because a runner holds only plain float / ndarray / list, the optimizer that owns the list of runners pickles for backup/resume exactly like the single-start machine did (ADR-0007). Being backend-free, a runner is also unit-testable offline by feeding it scores + gradients from an analytic function (no bngsim) – how the step math is validated against a scipy oracle and how the multi-start win is demonstrated.

A leaf subclass (trf.py / lbfgs.py) sets up its own model state in __init__ and implements got(); the orchestrator only ever calls start(), got(), and reads iteration / fval / stop_reason / progress_detail() for reporting.

got(u_point, score, grad)

Consume one completed evaluation and return the next u or DONE.

progress_detail()

A short method-specific status suffix for the per-iteration report.

start()

The first u-point to evaluate (the start point).

Metaheuristic optimizers

ParticleSwarm optimizer (the pso fit type).

Extracted from the algorithms package (M1 Step 4). Subclasses Algorithm and inherits its run loop + execution seam; it makes no core.* call of its own.

class pybnf.algorithms.optimizers.particle_swarm.PSOConfig(*, particle_weight: float = 0.7, adaptive_n_max: float = 30, adaptive_n_stop: float = inf, adaptive_abs_tol: float = 0.0, adaptive_rel_tol: float = 0.0, cognitive: float = 1.5, social: float = 1.5, v_stop: float = 0.0)

Particle-swarm config fields, co-located with the method (ADR-0002, ADR-0006). These defaults previously lived in config_schema.GlobalConfig; moving them here makes pso’s config knowledge travel with its algorithm class. The values are byte-identical to the old global defaults – the int-literal adaptive_n_max = 30 stays an int (pydantic does not validate defaults), matching the golden-config net.

particle_weight_final is intentionally NOT here: it is not a default but a runtime fallback (set to particle_weight when absent, in __init__ below), so it stays a pass-through extra.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

class pybnf.algorithms.optimizers.particle_swarm.ParticleSwarm(config)

Implements particle swarm optimization.

The implementation roughly follows Moraes et al 2015, although is reorganized to better suit PyBNF’s format. Note the global convergence criterion discussed in that paper is not used (would require too long a computation), and instead uses ????

add_iterations(n)

Adds n additional iterations to the algorithm. May be overridden in subclasses that don’t use self.max_iterations to track the iteration count

got_result(res)

Updates particle velocity and position after a simulation completes.

Parameters:

res – Result object containing the run PSet and the resulting Data.

Returns:

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Start the run by initializing n particles at random positions and velocities :return:

The Differential Evolution optimizer family (de and ade fit types).

DifferentialEvolutionBase is the shared base; DifferentialEvolution (de) and AsynchronousDifferentialEvolution (ade) subclass it. Extracted byte-identical (M1 Step 4). The family makes no core.* call of its own — the run loop and the execution seam are inherited from Algorithm.

class pybnf.algorithms.optimizers.differential_evolution.AsynchronousDifferentialEvolution(config)

Implements a simple asynchronous differential evolution algorithm.

Contains no islands or migrations. Instead, each time a PSet finishes, proposes a new PSet at the same index using the standard DE formula and whatever the current population happens to be at the time.

got_result(res)

Called when a simulation run finishes

Parameters:

res – Result object

Returns:

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

class pybnf.algorithms.optimizers.differential_evolution.DEFamilyConfig(*, mutation_rate: float = 0.5, mutation_factor: float = 0.5, stop_tolerance: float = 0.002, de_strategy: str = 'rand1')

Config fields shared by the whole DE family, co-located with DifferentialEvolutionBase (ADR-0002, ADR-0006) – exactly the keys that base __init__ reads. ade registers against this base directly (it adds no config keys); de extends it (DifferentialEvolutionConfig). Values are byte-identical to the old GlobalConfig defaults.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

class pybnf.algorithms.optimizers.differential_evolution.DifferentialEvolution(config)

Implements the parallelized, island-based differential evolution algorithm described in Penas et al 2015.

In some cases, I had to make my own decisions for specifics I couldn’t find in the original paper. Namely: At each migration, a user-defined number of individuals are migrated from each island. For each individual, a random index is chosen; the same index for all islands. A random permutation is used to redistribute individuals with that index to different islands.

Each island performs its migration individually, on the first callback when all islands are ready for that migration. It receives individuals from the migration iteration, regardless of what the current iteration is. This can sometimes lead to wasted effort. For example, suppose migration is set to occur at iteration 40, but island 1 has reached iteration 42 by the time all islands reach 40. Individual j on island 1 after iteration 42 gets replaced with individual j on island X after iteration 40. Some other island Y receives individual j on island 1 after iteration 40.

got_result(res)

Called when a simulation run finishes

This is not thread safe - the Scheduler must ensure only one process at a time enters this function. (or, I should rewrite this function to make it thread safe)

Parameters:

res – Result object

Returns:

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

class pybnf.algorithms.optimizers.differential_evolution.DifferentialEvolutionBase(config)

got_result(res)

Called by the scheduler when a simulation is completed, with the pset that was run, and the resulting simulation data

Parameters:

res (Result) – result from the completed simulation

Returns:

List of PSet(s) to be run next or ‘STOP’ string.

new_individual(individuals, base_index=None)

Create a new individual for the specified island, according to the set strategy

Parameters:

base_index – The index to use for the new individual, or None for a random index.

Returns:

start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

class pybnf.algorithms.optimizers.differential_evolution.DifferentialEvolutionConfig(*, mutation_rate: float = 0.5, mutation_factor: float = 0.5, stop_tolerance: float = 0.002, de_strategy: str = 'rand1', islands: int = 1, migrate_every: int = 20, num_to_migrate: int = 3)

de-specific config: the island/migration fields only synchronous DE reads (ade is async and ignores them). Demonstrates the shared-base pattern the MCMC family reuses (ADR-0006).

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

ScatterSearch optimizer (the ss fit type).

Extracted byte-identical (M1 Step 4). Subclasses Algorithm and inherits its run loop + execution seam; it makes no core.* call of its own.

class pybnf.algorithms.optimizers.scatter_search.ScatterSearch(config)

Implements ScatterSearch as described in the introduction of Penas et al 2017 (but not the fancy parallelized version from that paper). Uses the individual combination method described in Egea et al 2009

get_backup_every()

Overrides base method because Scatter Search runs n*(n-1) PSets per iteration.

got_result(res)

Called when a simulation run finishes

Parameters:

res (Result)

Returns:

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

class pybnf.algorithms.optimizers.scatter_search.ScatterSearchConfig(*, local_min_limit: int = 5)

Scatter-search config fields, co-located with the method (ADR-0002, ADR-0006). local_min_limit is the only defaulted scatter key (it was the lone # --- scatter search --- entry in GlobalConfig); init_size and reserve_size are NOT here – like PSO’s particle_weight_final they are runtime-defaulted in __init__ (init_size → 10*len(variables), reserve_size → max_iterations when absent), so they stay pass-through extras. Value byte-identical to the old global default.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

SimplexAlgorithm — the parallelized Nelder-Mead simplex (the sim fit type).

Also used as the refinement step after another fit. Extracted byte-identical (M1 Step 4). Subclasses Algorithm and sets _is_simplex = True so the inherited run() teardown always clears the Simulations directory at end-of-fit. The start value is mapped out of sampling space by the parameter’s own scale (FreeParameter.from_sampling_space), so log10 and natural-log vars both work.

class pybnf.algorithms.optimizers.simplex.SimplexAlgorithm(config, refine=False)

Implements a parallelized version of the Simplex local search algorithm, as described in Lee and Wiswall 2007, Computational Economics

START_POINT_KEY = 'simplex_start_point'

Internal config key the refiner start point is injected under, read by __init__ below and written by pybnf._refine_best_fit (the registry-driven refiner dispatch, #403/ADR-0015). Simplex predates StartPointOptimizer and keeps its own byte-identical start-point parsing, so this is declared here.

a_plus_b_times_c_minus_d(a, b, c, d, v)

Performs the calculation a + b*(c-d), where a, c, and d are assumed to be in log space if v is in log space, and the final result respects the box constraints on v.

Parameters:

• a

• b

• c

• d

• v (FreeParameter)

Returns:

ab_plus_cd(a, b, c, d, v)

Performs the calculation ab + cd where b and d are assumed to be in log space if v is in log space, and the final result respects the box constraints on v :param a: :param b: :param c: :param d: :param v: :type v: FreeParameter :return:

get_sums()

Simplex helper function Returns a dict mapping parameter name p to the sum of the parameter value over the entire current simplex :return: dict

got_result(res)

Called by the scheduler when a simulation is completed, with the pset that was run, and the resulting simulation data

Parameters:

res (Result) – result from the completed simulation

Returns:

List of PSet(s) to be run next or ‘STOP’ string.

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

class pybnf.algorithms.optimizers.simplex.SimplexConfig(*, simplex_step: float = 1.0, simplex_reflection: float = 1.0, simplex_expansion: float = 1.0, simplex_contraction: float = 0.5, simplex_shrink: float = 0.5, simplex_stop_tol: float = 0.0)

Simplex (Nelder-Mead) config fields, co-located with the method (ADR-0002, ADR-0006) – the six defaulted simplex_* knobs __init__ reads unconditionally. NOT here (stay pass-through extras): simplex_max_iterations and simplex_log_step are runtime-guarded (if 'X' in config.config) so they default at runtime, and simplex_start_point is an internal key the algorithm sets itself. The Simplex var/logvar-only free-parameter rule stays in config.py (cross-config knowledge, ADR-0006 #5). Values are byte-identical to the old GlobalConfig defaults.

The refine->simplex cross-fit_type reach (ADR-0006/0013): _refine_best_fit runs Simplex on a non-simplex config and reads simplex_*. Under narrowing these keys are in a sim fit’s config via this schema, and in any other fit with refine == 1 via the whole-schema overlay _build_config applies (the _REFINER_SCHEMA seam in config.py).

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

Simulated Annealing optimizer (the sa fit type).

A clean rewrite-to-spec (M2.2, ADR-0008). sa was historically built on the Bayesian sampler base and evaluated the posterior (prior + likelihood) in its Metropolis accept – a silent MAP estimator for normal/lognormal priors, unlike every other PyBNF optimizer. It is now a true optimizer: it minimizes the raw objective, using the prior only for the initial random draw (like de/pso/ss/sim). Box constraints are enforced by proposal reflection (FreeParameter.add) plus the parameter bounds, exactly as before.

The algorithm runs population_size independent annealing chains in parallel. Each chain makes a fixed-magnitude random-walk proposal and a Metropolis accept at its own inverse temperature beta; an accepted uphill (objective-increasing) move cools the chain (beta += cooling). A chain finishes when its beta reaches beta_max or it runs max_iterations; the run stops once every chain has finished. sa is deprecated (it warns at dispatch) but still runs.

class pybnf.algorithms.optimizers.simulated_annealing.SimulatedAnnealing(config)

Simulated annealing as a true optimizer (ADR-0008): minimizes the raw objective across population_size independent annealing chains.

add_iterations(n)

Extend the budget by n iterations, restarting any chain that had already finished at max_iterations (mirrors the resume behavior of the other iterative methods); chains that finished by reaching beta_max stay finished.

choose_new_pset(oldpset)

Perturb oldpset by a fixed-magnitude (step_size) random-walk move. The direction is isotropic; any component that would leave the box is reflected back inside (FreeParameter.add defaults to reflect=True), so the symmetric proposal keeps the plain Metropolis ratio valid.

got_result(res)

Called by the scheduler when a simulation is completed, with the pset that was run, and the resulting simulation data

Parameters:

res (Result) – result from the completed simulation

Returns:

List of PSet(s) to be run next or ‘STOP’ string.

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

class pybnf.algorithms.optimizers.simulated_annealing.SimulatedAnnealingConfig(*, step_size: float = 0.2, beta: list = <factory>, cooling: float = 0.01, beta_max: float = inf)

Simulated-annealing config fields, co-located with the method (ADR-0002, ADR-0008). Standalone – sa is an optimizer, no longer part of the MCMC family – carrying only the four knobs the algorithm reads. The defaults are byte-identical to the values these keys held under the old MCMC-family schema; under narrowing (ADR-0013) they appear only in an sa fit’s effective config.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

Powell’s conjugate-direction optimizer (the powell fit type, #403).

A native, derivative-free local optimizer (Numerical Recipes §10.7), implemented inside PyBNF’s run-loop contract – start_run / got_result only, no run() override (ADR-0007). It is the second user-selectable refiner alongside Simplex (refine_method = powell) and also runs standalone (fit_type = powell), started from a single var / logvar point exactly like Simplex.

Why native (not scipy.optimize.minimize): scipy is a blocking driver, which would force either a run() override or a bridging thread; reimplementing the method as an explicit, picklable state machine keeps the one shared run loop (so backup/resume work like every other method) and adds no dependency.

The method, in sampling space u (StartPointOptimizer): hold a set of n search directions (initially the coordinate axes). One cycle line-minimizes along each direction in turn; then – per the Numerical Recipes criterion – the net move of the cycle may replace the direction of greatest decrease, building up mutually conjugate directions that converge quadratically on a quadratic bowl.

Line minimization is a serial bracketing + Brent 1-D search (ADR-0016, superseding the original fixed-step parabola): geometric (golden) expansion from ±powell_step brackets a minimum, then Brent’s method (parabolic interpolation with a golden-section fallback) refines it to powell_line_tol. The search is confined to the feasible t-interval – the range of step lengths that keep every parameter inside its box – so the bound reflection never folds the 1-D slice and Brent sees the true objective; a minimum that wants to lie past a bound lands on the boundary. Unlike the fixed-step parabola, this follows long curved (non-quadratic) valleys and adapts its step. It is serial – one objective evaluation per reactor batch – so Powell no longer evaluates probes concurrently (it was only ever two jobs on a single-point search); CMA-ES is the parallel derivative-free method.

All state is plain numpy / float / list – no generator, no thread, including the _BrentLineSearch sub-state-machine – so Algorithm.backup can pickle the optimizer mid-run and run(resume=...) continues it.

class pybnf.algorithms.optimizers.powell.PowellAlgorithm(config, refine=False)

Powell’s conjugate-direction method as a picklable reactor state machine.

START_POINT_KEY = 'powell_start_point'

Refiner start-point key (see StartPointOptimizer / pybnf._refine_best_fit).

add_iterations(n)

Adds n additional iterations to the algorithm. May be overridden in subclasses that don’t use self.max_iterations to track the iteration count

got_result(res)

Called by the scheduler when a simulation is completed, with the pset that was run, and the resulting simulation data

Parameters:

res (Result) – result from the completed simulation

Returns:

List of PSet(s) to be run next or ‘STOP’ string.

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

class pybnf.algorithms.optimizers.powell.PowellConfig(*, powell_step: float = 1.0, powell_line_tol: float = 0.0001, powell_stop_tol: float = 1e-05)

Powell config fields, co-located with the method (ADR-0002, ADR-0006).

powell_step is the initial bracketing step in sampling space u (a factor of 10**powell_step for a log parameter) – the first probe distance along each direction, from which the line search expands a bracket and refines. powell_line_tol is the fractional tolerance to which each 1-D (Brent) line minimum is resolved. powell_stop_tol ends the search when a whole cycle improves the objective by less than that fraction. Like Simplex’s simplex_max_iterations, powell_max_iterations (the cycle budget) is runtime-guarded – it defaults to the global max_iterations when unset – so it is a valid key but not a schema field. powell_start_point is internal (the refiner injects it), so it is not modeled here either.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

CMA-ES – Covariance Matrix Adaptation Evolution Strategy (cmaes fit type, #403).

A native, derivative-free, black-box optimizer (Hansen & Ostermeier 2001), the third PyBNF refiner alongside Simplex and Powell (refine_method = cmaes) and a first-class standalone optimizer (fit_type = cmaes). Implemented inside the run-loop contract – start_run / got_result only, no run() override (ADR-0007) – and natively (no cma dependency).

Unlike the serial Simplex/Powell line searches, CMA-ES is population-based: it draws lambda candidates per generation from a multivariate normal N(m, sigma**2 C), so the generation evaluates in parallel and CMA-ES actually exploits PyBNF’s cluster – the same generation-synchronized reactor pattern as Differential Evolution (accumulate the whole generation, then update). The search runs in sampling space u (StartPointOptimizer), so log parameters adapt geometrically.

Each generation ranks the candidates by objective and updates the distribution mean m (weighted recombination of the best mu), the step size sigma (cumulative step-length adaptation via the conjugate evolution path p_sigma), and the covariance C (rank-one p_c update + rank-mu update). The constants follow Hansen’s standard defaults (the positive-weights formulation). lambda is the configured population_size (floored at 4); cmaes_sigma0 is the initial step in u; the run ends at max_iterations generations or when the largest principal step sigma*sqrt(max eig C) falls below cmaes_stop_tol.

CMA-ES is also a strong global optimizer over a bounded box, so it gains a box / global-start mode (#404, ADR-0017; registry start_from_box=True): given bounded uniform_var / loguniform_var priors instead of a var / logvar start point, it begins at the box center (StartPointOptimizer) and seeds the covariance with the per-coordinate box widths in u – so the initial per-coordinate standard deviation is cmaes_sigma0 * (box width) and the first generation spans the whole box. There cmaes_sigma0 is read as a fraction of each box width (default 0.3); in point-start / refine mode it is the absolute initial step in u as before. Candidates are repaired into the box by FreeParameter.set_value and the reflected point enters the update (see got_result), so bounds are respected. As a refiner (refine_method = cmaes) the start point is the injected best fit and the covariance starts isotropic – box mode applies only to a standalone bounded-prior fit.

All state is plain numpy / float / list (mean, sigma, covariance, evolution paths, the pending generation) – picklable, so backup/resume work like every other method.

class pybnf.algorithms.optimizers.cmaes.CMAESAlgorithm(config, refine=False)

CMA-ES as a picklable, generation-synchronized reactor state machine.

START_POINT_KEY = 'cmaes_start_point'

Refiner start-point key (see StartPointOptimizer / pybnf._refine_best_fit).

add_iterations(n)

Adds n additional iterations to the algorithm. May be overridden in subclasses that don’t use self.max_iterations to track the iteration count

got_result(res)

Called by the scheduler when a simulation is completed, with the pset that was run, and the resulting simulation data

Parameters:

res (Result) – result from the completed simulation

Returns:

List of PSet(s) to be run next or ‘STOP’ string.

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

class pybnf.algorithms.optimizers.cmaes.CMAESConfig(*, cmaes_sigma0: float = 0.3, cmaes_stop_tol: float = 1e-11)

CMA-ES config fields, co-located with the method (ADR-0002, ADR-0006).

cmaes_sigma0 is the initial overall step size in sampling space u (a factor of 10**cmaes_sigma0 for a log parameter) in point-start / refine mode; in box / global-start mode (bounded priors, #404/ADR-0017) it is instead read as a fraction of each box width, so the initial per-coordinate standard deviation is cmaes_sigma0 * (box width) – one knob, one default (0.3), interpreted relative to the start the fit actually uses. cmaes_stop_tol ends the run when the largest principal standard deviation of the search distribution shrinks below it. The population size lambda is the global population_size (consistent with de/pso/sa), and the generation budget is the global max_iterations, so neither is duplicated here. cmaes_start_point is internal (the refiner injects it), so it is not a schema field.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

Gradient-based optimizers

Trust-Region-Reflective least-squares optimizer (trf fit type, #386 / #460).

The primary gradient-based local optimizer: the workhorse for Gaussian / sum-of- squares objectives (the common case in PyBNF). It consumes the residual vector + residual-Jacobian #385 assembles from bngsim’s forward output sensitivities, and approximates the objective Hessian as JᵀJ – far better conditioned and faster- converging than feeding a scalar gradient to a generic quasi-Newton method on a least-squares problem. This is the D2D (Data2Dynamics) workhorse step; see examples/becker_d2d_gradient/ for the standalone methodological reference.

Why native (not scipy.optimize.least_squares): scipy is a blocking driver that calls fun/jac synchronously, so it cannot farm its evaluations to PyBNF’s distributed propose/score loop (the same incompatibility powell.py documents). The method is reimplemented as an explicit, picklable step machine – here a headless GradientRunner (_TRFRunner) that GradientOptimizer drives inside the run-loop contract, no run() override (ADR-0007) – so backup/resume work like every other method and one evaluation is one scheduler job. Factoring the step machine into a per-start runner is also what lets a fit run N of them concurrently (local multi-start, the orchestration the base owns). scipy is the test oracle for the step math (tests/test_gradient_runner.py), never called in the production loop – the Branch–Coleman–Li math below is ported in-process (the _trf_* helpers).

The method (Trust-Region-Reflective, Branch–Coleman–Li 1999)

In sampling space u (StartPointOptimizer / GradientOptimizer), with the residual r and residual-Jacobian J from #385 at the current point x (gradient g = Jᵀr, Gauss–Newton Hessian A = JᵀJ), the bounds are handled by the Coleman–Li affine scaling, matching scipy.optimize.least_squares(method='trf'):

• the Coleman–Li scaling v(x): vᵢ is the distance from xᵢ to the bound the anti-gradient points at (ubᵢ-xᵢ when gᵢ<0, xᵢ-lbᵢ when gᵢ>0, else 1). With D = diag(√v) the first-order optimality condition is D²g = 0 – so the scaled gradient ‖v·g‖∞ reads as zero both in the interior (g=0) and on an active face (v=0 there). A bound-active optimum is therefore recognized as optimal, where a plain ‖g‖∞ never would be;

• the trust-region subproblem is solved in the scaled (“hat”) space x = D x̂, where the Newton model is ½ x̂ᵀ(DAD + C)x̂ + (Dg)ᵀx̂ with C = diag(g·dv) the curvature the scaling itself contributes (dv ∈ {-1,0,1}), via a single SVD of the augmented Jacobian [J D; C^{½}] and More’s root-find on the trust radius Δ;

• reflective step selection: a trust-region step that would leave the box is reflected off the first bound it crosses, and the best (by predicted model reduction) of {the constrained trust step, the reflected step, a scaled-gradient Cauchy step} is taken, kept strictly interior by a step-back factor θ (the Coleman–Li differentiability requirement). Reflection lets the iterate slide along an active face and converge cleanly onto a bound-active corner, where simply clipping the unconstrained step into the box can stall.

Δ is the trust-region knob: it shrinks when a step is rejected or over-predicted (Δ ← ¼‖step‖ when the gain ratio ρ < ¼) and grows when a good step pushes against the trust boundary (Δ ← 2Δ when ρ > ¾ and the step hit the boundary). The initial Δ₀ is derived from the (scaled) start point, not a tunable.

Each iteration costs one objective evaluation (the trial). The whole scaled-SVD subproblem solve + reflective step selection is in-process linear algebra on the cached r/J – no evaluations. On accept, the trial’s own residual/Jacobian – already assembled, since master-side scoring returns the simdata – become the next iterate’s, so an accepted step needs no re-evaluation; on reject, only Δ shrinks and the step is re-solved from the same cached SVD. The run stops when the scaled gradient is flat (‖v·g‖∞ ≤ trf_grad_tol – the first-order optimality test that respects active bounds), the step is negligible (‖δ‖ ≤ trf_step_tol·(‖x‖+trf_step_tol)), or the iteration budget is spent. In the box interior (no bound active at the solution) the scaling and reflection fall away and the method is an ordinary trust-region least-squares step, MINPACK-style.

Scope. TRF consumes the exact least-squares residual – the Gaussian (any scale/location) and the Student-t (#459). A fit whose objective is not an exact sum of squares (an estimated noise scale, a Laplace / count family, active constraints; GradientResult.least_squares_exact == False) has no faithful residual model, so this optimizer refuses it with a pointer to the L-BFGS-B path (fit_type = lbfgs, #386’s fallback). Local multi-start is provided by GradientOptimizer (the base runs N independent _TRFRunner starts concurrently and keeps the global best).

All runner state is plain numpy / float (the point, residual, Jacobian, the cached scaling + SVD, the trust radius) – picklable, so Algorithm.backup checkpoints the optimizer (and its list of runners) mid-run.

class pybnf.algorithms.optimizers.trf.TRFAlgorithm(config, refine=False)

Bounded Trust-Region-Reflective least-squares: a method-agnostic multi-start orchestrator (GradientOptimizer) over per-start _TRFRunner step machines.

START_POINT_KEY = 'trf_start_point'

The internal config key the refiner start point is injected under (mirrors SimplexAlgorithm’s 'simplex_start_point'). Set by each subclass; pybnf._refine_best_fit writes the best fit here so refinement starts from it instead of parsing the (refiner-irrelevant) variable specs.

fit_type = 'trf'

Message label + refiner start-point key (see StartPointOptimizer).

class pybnf.algorithms.optimizers.trf.TRFConfig(*, trf_grad_tol: float = 1e-08, trf_step_tol: float = 1e-08)

TRF (Trust-Region-Reflective) config fields, co-located with the method (ADR-0006).

trf_grad_tol ends the run when the largest component of the scaled least-squares gradient v·Jᵀr (in sampling space) falls below it – a first-order optimality test that reads as the ordinary ‖Jᵀr‖∞ in the interior and as zero on a bound the gradient pushes against (the Coleman–Li scaling v vanishes there). trf_step_tol ends it when an accepted step δ becomes negligible relative to the point (‖δ‖ ≤ trf_step_tol·(‖x‖+trf_step_tol)). The initial trust radius is derived from the (scaled) start point – there is no damping/radius tunable. Like Powell’s cycle budget, trf_max_iterations is runtime-guarded – it defaults to the global max_iterations when unset – so it is a valid key but not a schema field. trf_start_point is internal (the refiner injects it), so it is not modeled here either.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

Bounded limited-memory BFGS optimizer – L-BFGS-B (lbfgs fit type, #386).

The scalar-gradient sibling of the trust-region least-squares method (trf.py): where TRF consumes the residual vector + residual-Jacobian and is the workhorse for exact-least-squares (Gaussian / Student-t, fixed scale) objectives, this method consumes only the scalar objective value (res.score) and the scalar gradient (gradient_at(res).gradient) that #385 assembles from bngsim’s forward output sensitivities. That makes it the fallback for precisely the objectives TRF refuses – the ones whose GradientResult.least_squares_exact is False and so have no faithful sum-of-squares residual: an estimated noise scale (the retained +log σ normalizer is not a square), the Laplace / count families, and a fit with active constraint penalties. The assembly returns a finite scalar gradient for all of these, so the same gradient seam drives both methods; only the step math differs.

Why native (not scipy.optimize.minimize(method='L-BFGS-B')): scipy is a blocking driver that calls fun/jac synchronously, so it cannot farm its evaluations to PyBNF’s distributed propose/score loop (the same incompatibility powell.py and trf.py document). The method is reimplemented as an explicit, picklable step machine – here a headless GradientRunner (_LBFGSRunner) that GradientOptimizer drives inside the run-loop contract, no run() override (ADR-0007) – so backup/resume work like every other method and one evaluation is one scheduler job. Factoring the step machine into a per-start runner is also what lets a fit run N of them concurrently (local multi-start, the orchestration the base owns); the runner itself is pure u-space numpy, so its math is unit-testable offline against a scipy oracle with no backend.

The method (L-BFGS-B: Cauchy point + subspace min + line search)

The Byrd–Lu–Nocedal–Zhu (1995) algorithm. In sampling space u (StartPointOptimizer / GradientOptimizer), with the scalar gradient g = ∇F from #385 at the current point x:

• keep the m most recent curvature pairs (s, y) (s = xₖ₊₁ - xₖ, y = gₖ₊₁ - gₖ) and form the limited-memory Hessian approximation B in its direct (not inverse) form – here a dense θI seed updated by the cached pairs (_form_b()), since PyBNF fits have small n (the paper’s compact θI − W M Wᵀ products only pay off for large n);

• generalized Cauchy point (_generalized_cauchy_point(), their Algorithm CP): walk the piecewise-linear projected steepest-descent path P[x − t·g, l, u], sort the per-coordinate breakpoints, and stop at the first local minimizer of the quadratic model m along it – the Cauchy point xᶜ and the active set (the coordinates pinned to a bound on the way);

• subspace minimization (_subspace_min(), their §5 dense primal): minimize the quadratic model over the free variables from xᶜ, holding the active coordinates at their bounds, truncated to stay feasible – the subspace minimizer x̄. The search direction is d = x̄ − x. Both stages are pure in-process linear algebra on the cached pairs + gradient – no objective evaluations;

• backtracking (Armijo) line search along d: start at step length α = 1 (the full step lands at x̄) and project each trial P[x + α d] onto the box (_u_bounds()); accept the first α whose projected sufficient-decrease holds, F(P[x+αd]) ≤ F(x) + c₁ · gᵀ(P[x+αd] - x) (Bertsekas’ projected Armijo, testing the actual displacement so a clamped coordinate is handled correctly; it reduces to the usual α gᵀd when no bound is active), else shrink α ← backtrack·α and retry;

• on accept, fold the realized (s, y) into the history (skipping a non-positive- curvature pair, which would break the BFGS positive-definite invariant) and step.

Each line-search trial costs one objective evaluation; an accepted step’s gradient is assembled from that trial’s own simdata (master-side scoring returns it), so the step needs no re-evaluation. The run stops when the projected gradient is flat (‖P[x-g]-x‖∞ ≤ lbfgs_grad_tol – the first-order optimality test that respects active bounds), the accepted step is negligible (‖s‖ ≤ lbfgs_step_tol·(‖x‖+lbfgs_step_tol)), the line search cannot reduce the objective even along steepest descent (a stalled minimum), or the iteration budget is spent.

In the box interior (no bound active at the solution) the Cauchy point falls short of its first breakpoint and the subspace minimization runs over every variable, so the step reduces exactly to the unconstrained limited-memory quasi-Newton step −B⁻¹g; the active-set machinery only bites when bounds are active at the optimum.

All runner state is plain numpy / float / list (the point, gradient, the (s, y) history, the line-search scratch) – picklable, so Algorithm.backup checkpoints the optimizer (and its list of runners) mid-run, exactly like Powell / CMA-ES / TRF (ADR-0007). B and the Cauchy / subspace work are recomputed per iteration from the cached pairs, so no dense matrix joins the persisted state.

class pybnf.algorithms.optimizers.lbfgs.LBFGSAlgorithm(config, refine=False)

Bounded limited-memory BFGS (L-BFGS-B): a method-agnostic multi-start orchestrator (GradientOptimizer) over per-start _LBFGSRunner step machines.

START_POINT_KEY = 'lbfgs_start_point'

The internal config key the refiner start point is injected under (mirrors SimplexAlgorithm’s 'simplex_start_point'). Set by each subclass; pybnf._refine_best_fit writes the best fit here so refinement starts from it instead of parsing the (refiner-irrelevant) variable specs.

fit_type = 'lbfgs'

Message label + refiner start-point key (see StartPointOptimizer).

class pybnf.algorithms.optimizers.lbfgs.LBFGSConfig(*, lbfgs_grad_tol: float = 1e-06, lbfgs_step_tol: float = 1e-08, lbfgs_history: int = 10, lbfgs_c1: float = 0.0001, lbfgs_backtrack: float = 0.5)

L-BFGS-B config fields, co-located with the method (ADR-0006).

lbfgs_grad_tol ends the run when the largest component of the projected gradient P[x-g]-x (in sampling space) falls below it – a first-order optimality test that reads as the ordinary ‖g‖∞ in the interior and as zero on a bound the gradient pushes against. lbfgs_step_tol ends it when an accepted step s becomes negligible relative to the point (‖s‖ ≤ lbfgs_step_tol·(‖x‖+lbfgs_step_tol)). lbfgs_history is m, the number of recent curvature pairs the limited-memory Hessian retains (more pairs ⇒ a richer Hessian model at higher per-step cost). lbfgs_c1 is the Armijo sufficient-decrease constant (0 < c₁ < 1, conventionally 1e-4) and lbfgs_backtrack the step-length reduction factor on a rejected trial (0 < β < 1). Like Powell’s / TRF’s cycle budget, lbfgs_max_iterations is runtime-guarded – it defaults to the global max_iterations when unset – so it is a valid key but not a schema field. lbfgs_start_point is internal (the refiner injects it), so it is not modeled here either.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

Profile-likelihood identifiability + confidence intervals (profile_likelihood job type, #446 / #466).

Profile likelihood is the Data2Dynamics (D2D) method for parameter identifiability and confidence intervals (Raue et al., Bioinformatics 25(15):1923-1929, 2009). For each fitted parameter theta_k it fixes theta_k to a grid of values around the optimum theta* and re-optimizes all the other parameters at each grid point, tracing the profile chi2_PL(theta_k) = min_{j != k} chi2(theta). The confidence interval is the range where the profile stays below a Delta chi2 threshold (the chi-square quantile at the chosen confidence level, 1 dof); a profile that stays flat diagnoses a structurally non-identifiable parameter, one that rises on only one side (or reaches a bound without crossing) a practically non-identifiable one, and one that crosses the threshold on both sides an identifiable one with a finite CI.

A standalone new-era job, not a fit stage (ADR-0031, #446)

This is a self-contained job_type = profile_likelihood run selected on the modern (edition >= 2) surface, not a stage auto-triggered at the end of a fit. It subclasses GradientOptimizer, so it inherits the whole gradient path – the edition-2 / sensitivity-backend / differentiable-dynamics gates, the per-experiment routing setup, and the per-evaluation GradientOptimizer.gradient_at() assembly – and depends on it exactly as trf / lbfgs do (bngsim forward sensitivities, BNGSIM_HAS_OUTPUT_SENS).

The inner optimizer – the one that both polishes to theta* and re-optimizes the remaining free parameters at each grid point – mirrors the trf/lbfgs split so profile likelihood fits every objective the gradient methods do, not just the exact-least-squares ones. The kind is chosen automatically from the objective’s structure, read once from an assembled GradientResult (a preflight evaluation on the polish path, or the explicit-theta* evaluation): an exact sum of squares (a Gaussian / Student-t family with a fixed noise scale and no constraints, least_squares_exact) profiles with the reduced-dimension _TRFRunner (the D2D reference, consuming the residual / Jacobian); anything else – an estimated noise scale, a Laplace / count family, active constraint penalties – has no faithful residual and profiles with the scalar-gradient _LBFGSRunner instead (the same objectives job_type = lbfgs handles). least_squares_exact is a structural property of the objective, not the point, so the single read governs the whole run.

Obtaining theta* (the open design question of #466)

Two sources, resolved at start_run():

• Explicit override. If every free parameter declares an initial_value: in its parameter: record (so each FreeParameter carries a value), those values are theta*. The job evaluates that one point (for the reference objective) and skips straight to profiling – the fast, deterministic path for a fit you have already run.

• Integrated polish. Otherwise the job first runs a gradient polish – the base’s own multi-start Trust-Region-Reflective fit over the bounded-prior box (population_size starts, max_iterations budget) – to find theta*, then profiles around it. One .conf does the whole fit-and-profile end to end.

Either way the Delta chi2 reference is the global objective minimum found across the whole run (trajectory.best_score()), so a profile scan that improves slightly on the polish optimum re-references cleanly rather than reporting a spuriously negative rise.

The driver (this issue’s core deliverable)

Profiling walks outward from theta* in sampling space u (priors/scale.py, ADR-0029 – log10 for a log-scaled parameter, so a D2D log10 grid is the u grid; the Delta chi2 threshold is on the objective, never on the transformed parameter). Each profiled parameter runs two independent _ProfileTracks (one per direction). A track takes an adaptive step – shrinking where the profile steepens, growing where it is flat, targeting a fixed Delta chi2 rise per grid point – and at each grid point re-optimizes the remaining free parameters with a reduced-dimension inner runner of the selected kind (a _TRFRunner with the fixed column dropped from the Jacobian and the full residual kept, or an _LBFGSRunner on the scalar gradient with the fixed coordinate dropped), warm-started from the neighboring grid point’s optimum. A direction terminates on the threshold crossing, on a bound, or on a per-direction point cap.

Parallel across parameters (#446 sub-issue 2 / #467)

The per-parameter profiles – and the two directions of a single parameter – are independent, so they farm across the existing Dask scheduler rather than running serially, exactly as the base’s local multi-start does: up to profile_likelihood_max_parallel directional tracks (default: all of them) run concurrently, each an inherently sequential warm-started walk holding one evaluation in flight, routed back to its owner by PSet name. A cap only queues the excess tracks – they run as slots free – so coverage is never silently truncated; a direction that stops at its point budget is reported as such.

Serialized + resumable outputs (#467)

Every state object here is plain numpy / float / list (the tracks, the inner runners, the per-point cost / theta_opt / nfev / success accumulated in self._profiles), so the optimizer pickles through PyBNF’s ordinary backup/resume machinery exactly like the other gradient methods (ADR-0007) – a resumed run re-submits only the in-flight tracks and never recomputes a finished one. At the end the driver writes the finished artifacts to Results/: the per-parameter curve files, the CI + identifiability summary table, and – when matplotlib (the optional pybnf[plot] extra) is importable – the profile plots (Delta chi2 panels with the threshold + CI lines, the reference notebook’s Cell 9); a missing matplotlib skips only the plots, never the run. scipy stays out of the production loop: the chi-square threshold comes from a dependency-free probit approximation (_chi2_quantile_1dof()).

class pybnf.algorithms.optimizers.profile_likelihood.ProfileLikelihoodAlgorithm(config, refine=False)

Standalone profile-likelihood driver (job_type = profile_likelihood, #446/#466).

A two-phase job over the GradientOptimizer gradient path: an optional multi-start gradient polish to the optimum theta* (skipped when the config supplies theta* via initial_value: on every parameter), then the profile phase – one adaptive outward _ProfileTrack per parameter per direction, up to profile_likelihood_max_parallel of them farmed across the scheduler concurrently (#467), each re-optimizing the remaining free parameters with a reduced-dimension inner runner. Both the polish and the re-optimizations use the same inner optimizer, selected once from the objective’s structure (_select_runner_kind()): the trust-region-reflective _TRFRunner for an exact sum of squares, the scalar-gradient _LBFGSRunner otherwise – so profile likelihood covers the estimated-scale / Laplace / count / constrained objectives too, not only the exact-least-squares ones. At the end it extracts each parameter’s confidence interval and identifiability class and writes the finished artifacts to Results/ – the per-parameter curves, the CI/classification summary table, and (with matplotlib) the profile plots. All per-point state pickles through PyBNF’s backup/resume, so a resumed run continues without recomputing finished tracks.

got_result(res)

Route a completed Result to the start that owns it, advance just that start’s runner, and return its next PSet – [] once it converges (other starts keep going), or 'STOP' only when the last live start finishes.

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Seed the n_starts runners and return their initial PSets (one evaluation each), so all starts run concurrently from the first scheduler batch.

class pybnf.algorithms.optimizers.profile_likelihood.ProfileLikelihoodConfig(*, profile_likelihood_params: Any = None, profile_likelihood_confidence: float = 0.95, profile_likelihood_step: float = 0.1, profile_likelihood_min_step: float = 0.001, profile_likelihood_max_step: float = 1.0, profile_likelihood_dchi2_target: float = 0.0, profile_likelihood_max_points: int = 40, profile_likelihood_reopt_max_iterations: int = 50, profile_likelihood_grad_tol: float = 1e-08, profile_likelihood_step_tol: float = 1e-08, profile_likelihood_max_parallel: int = 0)

Profile-likelihood config fields, co-located with the method (ADR-0006).

profile_likelihood_params selects which free parameters to profile (a list of ids; absent -> profile every free parameter). profile_likelihood_confidence is the CI confidence level (the Delta chi2 threshold is its chi-square 1-dof quantile). profile_likelihood_step is the initial outward step in sampling space u (log10 decades for a log-scaled parameter), adapted between profile_likelihood_min_step and profile_likelihood_max_step to hold each grid point’s objective rise near profile_likelihood_dchi2_target (0 -> auto, one tenth of the threshold). profile_likelihood_max_points caps the grid points per direction. profile_likelihood_reopt_max_iterations caps the inner re-optimization’s iterations at each grid point; profile_likelihood_grad_tol / profile_likelihood_step_tol are its (and the polish’s) optimality / step tolerances, shared by whichever inner optimizer the objective selects (Trust-Region-Reflective or L-BFGS-B). profile_likelihood_max_parallel caps how many per-direction profile walks run concurrently on the scheduler (0 -> all of them, one per direction per profiled parameter – the fully-parallel default; the walks that do not fit the cap queue and run as slots free, they are never dropped). Like the other gradient methods’ cycle budgets, profile_likelihood_max_iterations (the polish budget) is runtime-guarded – it defaults to the global max_iterations when unset – so it is a valid key but not a schema field.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

Bayesian samplers

The BayesianAlgorithm base class — shared by all Bayesian sampler fit types.

Extracted from the algorithms package (M1 Step 3). Holds posterior/prior bookkeeping and the rank-normalized split-R-hat / ESS convergence diagnostics (Vehtari et al. 2021). The run loop and execution seam live in the Algorithm base (..base); BayesianAlgorithm inherits them and adds no core.* call of its own, so this module does not import the core seam.

class pybnf.algorithms.samplers.base.BayesianAlgorithm(config)

Superclass for Bayesian MCMC algorithms

check_convergence(iteration, max_rhat)

Check if R-hat has converged below threshold. Returns True if should stop.

cleanup()

Called when quitting due to error. Save the histograms in addition to the usual algorithm cleanup

compute_ess()

Bulk and tail effective sample size per parameter (Vehtari et al. 2021).

Thin glue over pybnf.diagnostics.ess() (ADR-0009): reads this instance’s chain history and delegates the pure math. Returns (bulk_ess, tail_ess) arrays of shape (n_dim,), or (None, None).

compute_rhat()

Rank-normalized split-R-hat per parameter (Vehtari et al. 2021).

Thin glue over pybnf.diagnostics.rhat() (ADR-0009): reads this instance’s chain history and delegates the pure math. Returns an (n_dim,) array, or None if there is insufficient history.

evaluate_constraints(simdata, chain_index)

Evaluate all constraints against simulation data and cache the pass/fail results for this chain.

Parameters:

• simdata – Simulation data dictionary

• chain_index – Index of the chain that was accepted

got_result(res)

Called by the scheduler when a simulation is completed, with the pset that was run, and the resulting simulation data

Parameters:

res (Result) – result from the completed simulation

Returns:

List of PSet(s) to be run next or ‘STOP’ string.

ln_prior(pset)

Returns the value of the prior distribution for the given parameter set

Parameters:

pset (PSet)

Returns:

float value of ln times the prior distribution

load_priors()

Builds the data structures for the priors, based on the variables specified in the config.

record_pointwise_loglik(res, chain_index)

Cache the just-accepted pset’s per-observation log-likelihoods for this chain (ADR-0056), mirroring evaluate_constraints(): computed here, where res.simdata for the accepted pset is in hand, and written later by sample_pset() for whichever pset is current at a sample iteration. A no-op unless output_inference_data is set and the objective is a likelihood; any failure is logged, never fatal (the run must not die for a diagnostics sidecar).

report_constraint_satisfaction(file_ext)

Read the constraint samples file and write a summary of constraint satisfaction percentages. :param file_ext: String to append to the output file name

report_convergence_diagnostics(iteration)

Compute and report R-hat, ESS, and ESS/evaluation. Called every 10 iterations. Returns max_rhat for convergence checking, or None.

sample_pset(pset, ln_prob, chain_index=None)

Adds this pset to the set of sampled psets for the final distribution. :param pset: :type pset: PSet :param ln_prob - The probability of this PSet to record in the samples file. :type ln_prob: float :param chain_index: Index of the chain, used to look up cached constraint results. :type chain_index: int or None

start_run(setup_samples=True)

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

update_histograms(file_ext)

Updates the files that contain histogram points for each variable :param file_ext: String to append to the save file names :type file_ext: str :return:

class pybnf.algorithms.samplers.base.MCMCFamilyConfig(*, step_size: float = 0.2, burn_in: int = 10000, sample_every: int = 100, output_hist_every: int = 100, hist_bins: int = 10, adaptive: int = 10000, credible_intervals: list = <factory>, beta: list = <factory>, continue_run: int = 0, rhat_threshold: float = 0.0, diagnostics_every: int = 0)

Config shared by the whole Bayesian/MCMC family (mh, pt, am, dream, p_dream), co-located with the family base (ADR-0002, ADR-0006). (sa was once in this family; M2.2/ADR-0008 made it a standalone optimizer.)

Holds the defaulted keys the shared BayesianAlgorithm base reads (so they are common to every MCMC method) plus the β-ladder postprocess hook. The per-method keys live on the leaf models that subclass this and co-locate with their algorithm class: BasicMCMCConfig (mh/pt/sa, in basic_mcmc.py), AdaptiveMCMCConfig (am, in adaptive_mcmc.py), DreamConfig (dream, in dream.py) and PDreamConfig (p_dream, in pdream.py). neg_bin_r stays in GlobalConfig – an objfunc param read regardless of fit_type. Values are byte-identical to the old global defaults; under narrowing (ADR-0013) these keys appear only in each MCMC fit_type’s own effective config. The per-fit_type validity of individual keys is reported by check_unused_keys in config.py.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

classmethod postprocess(conf_dict, fit_type)

The β-ladder (ported from Configuration.postprocess_mcmc_keys).

Algorithms ‘mh’, ‘pt’, ‘am’, ‘dream’, ‘p_dream’ have similar but non-identical valid config keys; this builds beta_list / reps_per_beta and reconciles exchange_every / population_size. Mutates conf_dict in place (operating on the RAW config dict, before defaults merge, so raw-presence checks like 'beta' not in conf_dict mean “user did not set it”). config._build_config dispatches this uniformly; non-MCMC models inherit the no-op postprocess from PyBNFConfigModel (ADR-0006 #3).

Only the config transformations live here now. The per-method unused-key warnings this hook used to also emit (exchange_every/reps_per_beta on non-pt, the sa-only cooling/beta_max, the DREAM-only crossover_number/zeta/lambda/gamma_prob on mh/pt/am) moved to the unified, schema-derived Configuration.check_unused_keys (#401, ADR-0014), which now warns about every foreign key each MCMC fit does not own – more precisely than this hook’s hand-listed branches did.

BasicBayesMCMCAlgorithm — Metropolis-Hastings and Parallel Tempering (the mh and pt fit types).

mh is deprecated but still runs; pt is a working method sharing this class. (Simulated annealing, formerly sa on this class, is now a standalone optimizer — see optimizers/simulated_annealing.py, M2.2/ADR-0008.) Subclasses the sampler base (BayesianAlgorithm) and inherits the run loop + execution seam from Algorithm.

class pybnf.algorithms.samplers.basic_mcmc.BasicBayesMCMCAlgorithm(config)

Implements a Bayesian Markov chain Monte Carlo simulation. This is essentially a non-parallel algorithm, but here, we run n instances in parallel, and pool all results. This will give you a best fit (which is maybe not great), but more importantly, generates an extra result file that gives the probability distribution of each variable. This distribution depends on the prior, which is specified according to the variable initialization rules.

add_iterations(n)

Adds n additional iterations to the algorithm. May be overridden in subclasses that don’t use self.max_iterations to track the iteration count

choose_new_pset(oldpset, index)

Helper function to perturb the old PSet, generating a new proposed PSet. The step is a fixed-magnitude (step_size) random-walk move; any component that would leave the box is reflected back inside (see FreeParameter._reflect). :param oldpset: The PSet to be changed :type oldpset: PSet :param index: The chain index, selecting that chain’s own Generator :type index: int :return: the new PSet

cleanup()

Called when quitting due to error. Save the histograms in addition to the usual algorithm cleanup

got_result(res)

Called by the scheduler when a simulation is completed, with the pset that was run, and the resulting simulation data :param res: PSet that was run in this simulation :type res: Result :return: List of PSet(s) to be run next.

replica_exchange()

Performs replica exchange for parallel tempering. Then proposes n new parameter sets to resume all chains after the exchange. :return: List of n PSets to run

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

should_sample(index)

Checks whether this replica index is one that gets sampled. For mcmc, always True. For pt, must be a replica at the max beta

start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run. :return: list of PSets

try_to_choose_new_pset(index)

Helper function Advances the iteration number, and tries to choose a new parameter set for chain index i If that fails (e.g. due to a box constraint), keeps advancing iteration number and trying again. If it hits an iteration where it has to stop and wait (a replica exchange iteration or the end), returns None Otherwise returns the new PSet. :param index: :return:

class pybnf.algorithms.samplers.basic_mcmc.BasicMCMCConfig(*, step_size: float = 0.2, burn_in: int = 10000, sample_every: int = 100, output_hist_every: int = 100, hist_bins: int = 10, adaptive: int = 10000, credible_intervals: list = <factory>, beta: list = <factory>, continue_run: int = 0, rhat_threshold: float = 0.0, diagnostics_every: int = 0, exchange_every: ~typing.Any = 20)

Config for the basic-MCMC codes mh/pt, co-located with the method (ADR-0006). Adds the one key those samplers read (exchange_every) on top of the shared family fields; the β-ladder postprocess hook is inherited. exchange_every is typed Any because the hook overwrites it with np.inf for the non-PT (mh) method (it holds both an int and an infinity).

(cooling/beta_max were sa-only; M2.2/ADR-0008 moved them to SimulatedAnnealingConfig when sa became a standalone optimizer.)

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

Adaptive_MCMC — the Adaptive Metropolis sampler (the am fit type).

PyBNF’s recommended Bayesian sampler. Extracted byte-identical (M1 Step 4). Subclasses the sampler base (BayesianAlgorithm) and inherits the run loop + execution seam from Algorithm.

class pybnf.algorithms.samplers.adaptive_mcmc.AdaptiveMCMCConfig(*, step_size: float = 0.2, burn_in: int = 10000, sample_every: int = 100, output_hist_every: int = 100, hist_bins: int = 10, adaptive: int = 10000, credible_intervals: list = <factory>, beta: list = <factory>, continue_run: int = 0, rhat_threshold: float = 0.0, diagnostics_every: int = 0, stablizingCov: float = 0.001, calculate_covari: ~typing.Any = None)

Config for adaptive MCMC (am), co-located with the method (ADR-0006). Adds the covariance-adaptation keys Adaptive_MCMC reads on top of the shared family fields; the β-ladder postprocess hook is inherited.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

class pybnf.algorithms.samplers.adaptive_mcmc.Adaptive_MCMC(config)

got_result(res)

Called by the scheduler when a simulation is completed, with the pset that was run, and the resulting simulation data :param res: PSet that was run in this simulation :type res: Result :return: List of PSet(s) to be run next.

pick_new_pset(idx)

Parameters:

idx – Index of PSet to update

Returns:

A mew

reset(bootstrap=None)

Resets the Algorithm, keeping loaded variables and models

Parameters:

bootstrap (int or None) – The bootstrap number (None if not bootstrapping)

Returns:

A rejected bootstrap replicate (objective over bootstrap_max_obj) is retried with bootstrap_attempt incremented by the caller; that advances the RNG sub-stream so the retry draws a fresh resample and fit rather than repeating the identical failing run. bootstrap_attempt == 0 (the first try, and every non-bootstrap reset) reproduces the historical replicate-only seeding byte for byte.

start_run()

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run. :return: list of PSets

update_histograms(file_ext)

Updates the files that contain histogram points for each variable :param file_ext: String to append to the save file names :type file_ext: str :return:

DreamAlgorithm — the DREAM(ZS) sampler (the dream fit type).

DiffeRential Evolution Adaptive Metropolis drawing proposal donors from a growing ZS archive (Vrugt 2016). Extracted byte-identical (M1 Step 4). Subclasses the sampler base (BayesianAlgorithm) and inherits the run loop + execution seam from Algorithm; it makes no core.* call of its own.

class pybnf.algorithms.samplers.dream.DreamAlgorithm(config)

Implements a variant of the DREAM algorithm as described in Vrugt (2016) Environmental Modelling and Software.

Adapts Bayesian MCMC to use methods from differential evolution for accelerated convergence and more efficient sampling of parameter space

calculate_new_pset(idx)

Uses differential evolution-like update to calculate new PSet. Returns (PSet, cr_idx) or (None, cr_idx) if the proposal is out of bounds.

Parameters:

idx – Index of PSet to update

Returns:

tuple of (PSet or None, int)

calculate_snooker_pset(idx)

Snooker update proposal (ter Braak & Vrugt, 2008). Projects archive points onto the line through the current state and a reference archive point, then jumps along that axis.

Returns (PSet or None, log_correction) where log_correction is the log of the Hastings correction factor (d-1)*log(||Xp - Zc|| / ||X - Zc||).

detect_and_reset_outliers()

Detect outlier chains using the configured method on mean log-posteriors (last 50% of history). Reset outlier chains to copies of randomly selected non-outlier chains.

Methods: ‘iqr’ (interquartile range), ‘grubbs’ (Grubbs test at alpha=0.01).

got_result(res)

Called by the scheduler when a simulation is completed, with the pset that was run, and the resulting simulation data

Parameters:

res (Result) – PSet that was run in this simulation

Returns:

List of PSet(s) to be run next.

start_run(setup_samples=True)

Called by the scheduler at the start of a fitting run. Must return a list of PSets that the scheduler should run.

Algorithm subclasses optionally may set the .name field of the PSet objects to give a meaningful unique identifier such as ‘gen0ind42’. If so, they MUST BE UNIQUE, as this determines the folder name. Uniqueness will not be checked elsewhere.

Returns:

list of PSets

class pybnf.algorithms.samplers.dream.DreamConfig(*, step_size: float = 0.2, burn_in: int = 10000, sample_every: int = 100, output_hist_every: int = 100, hist_bins: int = 10, adaptive: int = 10000, credible_intervals: list = <factory>, beta: list = <factory>, continue_run: int = 0, rhat_threshold: float = 0.0, diagnostics_every: int = 0, gamma_prob: float = 0.1, zeta: float = 1e-06, lambda_: float = 0.1, crossover_number: int = 3, adaptive_step_size: ~typing.Any = True, archive_size: int | None = None, archive_thin_rate: int = 10, snooker_prob: float = 0.1, delta: int = 1, outlier_method: str = 'iqr')

Config for DREAM(ZS) (dream), co-located with the method (ADR-0006). Adds the DREAM-specific keys on top of the shared family fields; the β-ladder postprocess hook is inherited. PDreamConfig (p_dream) subclasses this and adds precondition_adapt. adaptive_step_size is typed Any so a user int 0/1 stays an int (matching the old behavior). Values byte-identical to the old global defaults.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

classmethod postprocess(conf_dict, fit_type)

The MCMC β-ladder (inherited from MCMCFamilyConfig) plus the DREAM-family step_size coupling (ADR-0013): a user-set step_size pins adaptive_step_size off. Both keys are owned only by DREAM/P-DREAM, so the coupling is intra-family – it formerly lived as a global write in config.py that, post-narrowing, would have stamped an orphan adaptive_step_size onto any mh/am/sa fit that set step_size. PDreamConfig inherits this (it owns both keys too).

PDreamAlgorithm — Preconditioned DREAM (the p_dream fit type).

DREAM(ZS) with proposals computed in a covariance-whitened parameter space, for better sampling of correlated posteriors. Extracted byte-identical (M1 Step 4). Subclasses DreamAlgorithm (a sibling leaf in this sub-package).

class pybnf.algorithms.samplers.pdream.PDreamAlgorithm(config)

P-DREAM: Preconditioned DREAM.

Extends DREAM(ZS) by computing DE proposals in a covariance-whitened parameter space. An online covariance estimate C is learned from the chain history (as in Adaptive Metropolis). Donors are transformed to whitened coordinates z = L_inv @ x before computing DE differences, and crossover is applied in whitened space where dimensions are decorrelated.

The proposal remains symmetric (DE differences from an external archive), so standard Metropolis-Hastings acceptance is valid without additional Hastings correction.

After a configurable adaptation period, the covariance is updated every generation from the pooled chain history. Before adaptation begins, the sampler behaves identically to DREAM(ZS).

calculate_new_pset(idx)

DE proposal in whitened space.

When preconditioning is active: 1. Transform current state and archive donors to z = L_inv @ x 2. Compute DE difference in z-space 3. Apply crossover in z-space (dimensions are decorrelated) 4. Scale and add perturbation in z-space 5. Convert the total jump back: dx = L @ dz_total 6. Propose x_new = x_current + dx

Before preconditioning activates, falls back to standard DREAM proposals.

got_result(res)

Override to update covariance estimate after each generation sync.

class pybnf.algorithms.samplers.pdream.PDreamConfig(*, step_size: float = 0.2, burn_in: int = 10000, sample_every: int = 100, output_hist_every: int = 100, hist_bins: int = 10, adaptive: int = 10000, credible_intervals: list = <factory>, beta: list = <factory>, continue_run: int = 0, rhat_threshold: float = 0.0, diagnostics_every: int = 0, gamma_prob: float = 0.1, zeta: float = 1e-06, lambda_: float = 0.1, crossover_number: int = 3, adaptive_step_size: ~typing.Any = True, archive_size: int | None = None, archive_thin_rate: int = 10, snooker_prob: float = 0.1, delta: int = 1, outlier_method: str = 'iqr', precondition_adapt: ~typing.Any = None)

Config for preconditioned DREAM (p_dream), co-located with the method (ADR-0006). Extends DreamConfig with the one preconditioning key PDreamAlgorithm adds; everything else (incl. the β-ladder hook) is inherited.

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

HMCSampler – the hmc fit type: blackjax NUTS on the analytical model’s JAX log-density (issue #425, ADR-0059).

This is an opt-in reference sampler: a gradient-based NUTS run that exists to evaluate PyBNF’s gradient-free samplers (am / dream / p_dream) against a research-grade yardstick on the canonical stress geometries. It runs only on the analytical / bring-your-own-log-density model (pybnf/analytical_model.py), never on a simulator – a BNGL/SBML posterior provides no cheap gradient, and simulator-path HMC is rejected outright by the ADR.

Unlike every other sampler, hmc does not dispatch parameter sets through the dask -> execute() -> score-column -> DirectPassObjective loop: the gradient cannot survive the per-pset dask round-trip. Instead it builds a JAX logdensity_fn in process from (a) the model’s JAX NLL (AnalyticalModel.nll_jax) and (b) the priors’ JAX log-densities (Prior.logpdf_jax), hands it to blackjax NUTS with window adaptation, and writes the draws in the standard samples format so the ArviZ bridge (ADR-0055), the LOO/WAIC sidecar machinery, and the rank-normalized split-R-hat / bulk-tail ESS diagnostics all work unchanged.

It samples an unconstrained z in R^d (ADR-0059 item 5), mapped to sampling space u (ADR-0010) by a per-parameter support-aware bijection u = b(z) (pybnf.priors.bijector), with target

log pi(z) = sum_i [ prior_i.logpdf_jax(b_i(z_i)) + log|b_i’(z_i)| ] + ( -NLL( scale.inverse(b(z)) ) )

The prior is defined in u (ADR-0010), so there is no theta <-> u Jacobian – the only change-of-variables term is the unconstraining bijection’s log|b'(z)|. For an unbounded prior b is the identity (z == u) and this is exactly the density am samples, now differentiated; for a positive/bounded/truncated prior b is a log/logit reparameterization that lands strictly inside the open support for every finite z, so the -inf support wall (and the NUTS divergence it once caused) is unreachable. Diagnostics and the samples file report u (the recorded Ln_probability is un-Jacobianed back to log pi(u)), so HMC stays comparable to the gradient-free samplers on the same posterior coordinate.

It samples the closed-form-truth and stress-geometry menu targets (gaussian / rotated_gaussian / banana / multimodal), the full edition-2 prior catalog – every family supplies a hand-written, scipy-logpdf-oracle-checked logpdf_jax (ADR-0059 item 4), sampled divergence-free on any support via the bijection (item 5) – and log-scaled parameters (lognormal_var / loguniform_var), whose u -> theta inverse traces through JAX (Scale.inverse_jax). It samples every built-in menu target (rotated_quartic included now) and bring-your-own ``expression`` targets – a user’s own PEtab-math NLL lambdified to JAX (ExpressionModel.nll_jax, ADR-0059 item 2), so HMC is a reference sampler on arbitrary closed-form log-densities, not just the canned menu. The only HMC work still deferred is the Python callable = module:func BYO form (a general callable is not JAX-traceable; ADR-0050).

jax/blackjax are the optional pybnf[jax] extra (ADR-0019): only this module (and the lazily-imported nll_jax / logpdf_jax) touches them, and a missing install surfaces as a pointed PybnfError naming the extra, never a bare ImportError.

class pybnf.algorithms.samplers.hmc.HMCConfig(*, step_size: float = 0.2, burn_in: int = 10000, sample_every: int = 100, output_hist_every: int = 100, hist_bins: int = 10, adaptive: int = 10000, credible_intervals: list = <factory>, beta: list = <factory>, continue_run: int = 0, rhat_threshold: float = 0.0, diagnostics_every: int = 0, num_samples: int = 1000, num_warmup: int = 1000, target_accept: float = 0.8)

Config for the hmc sampler, co-located with the method (ADR-0006).

Reuses the global population_size as the number of independent NUTS chains (so the samples-file chain naming iter<draw>run<chain> and the per-chain diagnostics carry over unchanged) and adds the three NUTS knobs on top of the shared MCMC family fields. Window adaptation replaces the gradient-free burn_in/sample_every thinning, so those inherited keys are unused by hmc (NUTS draws are near-independent – every post-warmup draw is kept).

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'populate_by_name': True, 'validate_by_alias': True, 'validate_by_name': True}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

num_samples: int

Post-warmup draws kept per chain (each becomes one samples.txt row).

num_warmup: int

Window-adaptation (warmup) steps per chain – dual-averaging step size + mass matrix.

target_accept: float

NUTS dual-averaging target acceptance probability (Stan-like default).

class pybnf.algorithms.samplers.hmc.HMCSampler(config)

blackjax NUTS on the analytical model’s JAX log-density (ADR-0059).

Subclasses BayesianAlgorithm to reuse prior loading, the samples-file setup and writer (start_run / sample_pset), and the R-hat/ESS diagnostics; overrides run() to drive NUTS in process instead of the dask score-column loop.

divergences

Post-warmup divergent-transition count per chain (filled in by run()). A NUTS-specific reliability signal the gradient-free samplers have no analogue for: divergences flag regions the leapfrog integrator cannot traverse (a too-sharp curvature for the tuned step size), so a nonzero count – like a high split-R-hat – means HMC’s own draws are not yet trustworthy on this geometry (ADR-0059 gates the reference on its own diagnostics: split-R-hat / ESS / divergences).

run(client=None, resume=None, debug=False)

Run blackjax NUTS in process and write draws in the standard samples format.

client (the dask client the harness passes every sampler) is intentionally ignored: a single analytical NUTS chain is a tight numeric loop and the chains run as independent blackjax runs, so there is no per-pset dispatch (ADR-0059 item 1).

Model checking

ModelCheck – the check fit type, a first-class checking method.

A checking run (statistical model checking): evaluates the objective value and constraint satisfaction for the given parameters without searching parameter space. It registers in the checker family – a peer of optimizer and sampler, not a utility afterthought – but, being neither an optimizer nor a sampler, it lives at the algorithms package top level rather than under optimizers/ or samplers/. Extracted byte-identical (M1 Step 4). It does not subclass Algorithm, but it does use the execution machinery: the patched names are resolved as core.Job / core.run_job through the core module object (ADR-0001 seam), and ConstraintCounter is read from this module’s namespace (tests patch it here).

class pybnf.algorithms.model_check.ModelCheck(config)

An algorithm that just checks the fit quality for a job with no free parameters.

Does not subclass Algorithm. To run, instead call run_check() with no Cluster.

run_check(debug=False)

Main loop for executing the algorithm