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).

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

The gradient-based leg of the concurrent multi-start base (#386/#500).

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(); the shared ConcurrentMultiStartOptimizer owns the start_run / got_result orchestration – seeding the runners, the name routing, reporting, and the multi-start STOP coordination – and this class fills in what the gradient path does differently: the pre-flight gates + master scoring, the u <-> PSet plumbing, the sensitivity-path activation (_setup_gradient_path(), hung on _pre_seed()), the gradient_at() assembly, and the gradient-consuming _advance(). 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 – the gradient path predates the n_starts field, hence _n_starts_key) – 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. 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.

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.

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.

The orchestrator passes grad = None with a non-finite score for a failed simulation (a non-integrable candidate point). A leaf’s got must tolerate that: mid-search its isfinite(score) guard already rejects the trial without dereferencing the gradient (the fit backs off); at the start point (phase == 'init', grad is None) it must terminate via _failed_start(). See GradientOptimizer._advance() (#492).

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(*, n_starts: int = 1, 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), plus the shared n_starts multi-start field (MultiStartConfig, #498). 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

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.

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.AsyncDEConfig(*, mutation_rate: float = 0.5, mutation_factor: float = 0.5, stop_tolerance: float = 0.002, de_strategy: str = 'rand1', n_starts: int = 1)

ade-specific config: the shared DE family fields plus the n_starts multi-start field (MultiStartConfig, #498/ADR-0071), and nothing else. ade opts into multi-start through its own subclass – mirroring how de extends the shared DEFamilyConfig base (DifferentialEvolutionConfig) – so the ADR-0006 “ade adds no keys to the family base” seam stays intact: n_starts rides this subclass, not the base de and ade share. Unlike de, ade has no islands or migrations, so this adds no fields of its own – it is exactly DEFamilyConfig + n_starts (#501).

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.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.

Opts into n_starts sequential-restart multi-start (#498/#501): the MultiStartOptimizer mixin (before the family base in the MRO) runs n_starts independent searches and keeps the global best. ade is the async one-in-one-out case the mixin’s draining path was built for – a full population stays in flight at each inner STOP, so the mixin drains those stragglers (already scored into the trajectory) before seeding the next start. n_starts == 1 (the default) is byte-identical to the single-run behavior.

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.

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, and nothing more. Neither de nor ade registers against this base directly: each extends it with its own subclass carrying the shared n_starts multi-start field (MultiStartConfig) – de also adds the island/migration fields (DifferentialEvolutionConfig), ade adds nothing else (AsyncDEConfig). Keeping the shared base itself key-minimal preserves the ADR-0006 “ade adds no keys to the family base” seam. 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.

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.

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', n_starts: int = 1, 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), plus the shared n_starts multi-start field (MultiStartConfig, #498). The n_starts key rides each method’s own subclass, not the shared DEFamilyConfig base, so the ADR-0006 “ade adds no keys to the family base” seam stays intact – ade gets n_starts through its own AsyncDEConfig (#501). 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.

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.

class pybnf.algorithms.optimizers.scatter_search.ScatterSearchConfig(*, n_starts: int = 1, 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. Inherits the shared n_starts multi-start field (MultiStartConfig, #498).

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. Subclasses StartPointOptimizer (the shared start-point / u <-> PSet plumbing it now shares with Powell and CMA-ES) and sets _is_simplex = True so the inherited run() teardown always clears the Simulations directory at end-of-fit. It runs from a single var / logvar start point or, over a bounded-prior box, from the box center – its global-start mode (start_from_box, #404/ADR-0017), which also unlocks concurrent multi-start (n_starts, #498). The start value is mapped out of sampling space by the parameter’s own scale, so log10 and natural-log vars both work.

Concurrency (#498). One simplex fans out only parallel_count (roughly n_variables - 1) probes per generation. SimplexRunner factors the whole Nelder-Mead state machine into a headless, picklable per-start object so LocalMultiStartOptimizer can run n_starts of them concurrently (n_starts x parallel_count probes in flight), keeping the global best. A point-start / refine fit is a single runner, byte-identical to the pre-multi-start behavior; each runner keeps Simplex’s exact box-clamped reflection / expansion / contraction arithmetic.

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

The parallelized Nelder-Mead simplex (Lee and Wiswall 2007), as n_starts concurrent SimplexRunners.

A standalone box fit (bounded priors) runs n_starts independent starts concurrently – start 0 from the box center, the rest from Latin-hypercube samples – and keeps the global best (LocalMultiStartOptimizer). A point-start (var / logvar) or refiner-injected fit is a single runner, byte-identical to the pre-multi-start behavior. Sets _is_simplex so the inherited run() teardown clears the Simulations directory.

START_POINT_KEY = 'simplex_start_point'

Internal config key the refiner start point is injected under, resolved by the shared StartPointOptimizer._resolve_start_pset (box center / point / injected) and written by pybnf._refine_best_fit (the registry-driven refiner dispatch, #403/ADR-0015).

class pybnf.algorithms.optimizers.simplex.SimplexConfig(*, n_starts: int = 1, 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-vs-prior free-parameter rule stays in config.py (cross-config knowledge, ADR-0006 #5). Inherits the shared n_starts multi-start field (MultiStartConfig, #498): sim in box / global-start mode runs n_starts concurrent starts, keeping the global best. 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].

class pybnf.algorithms.optimizers.simplex.SimplexRunner(variables, rng, start_point, start_steps, max_iterations, parallel_count, alpha, gamma, beta, tau, stop_tol)

Headless, picklable per-start Nelder-Mead simplex state machine (#498).

Simplex’s parallelized state machine (Lee and Wiswall 2007), factored out of the optimizer so a single fit can run N of them concurrently – local multi-start, the diversity a purely local method otherwise lacks. A runner owns one start’s entire mutable state – the simplex vertices, the per-iteration reflection / expansion / contraction bookkeeping, the pending map – and works in PSet space with Simplex’s exact box-clamped arithmetic; it knows nothing about the trajectory, backup, or dask. The orchestrator (LocalMultiStartOptimizer) drives it:

psets = runner.start()                 # the initial n+1 simplex vertices
nxt = runner.got(pset, score)          # -> list[PSet] | DONE

where got consumes one completed evaluation (routed to its vertex by PSet name) and returns the next batch of PSets to evaluate, or the DONE sentinel when this start has spent its iteration budget or its moves fell below simplex_stop_tol. All state is plain list / dict / int plus PSets, so the optimizer that owns the list of runners pickles for backup/resume exactly like the single-start machine did (ADR-0007).

Index convention (unchanged from the original single-simplex code): in the per-iteration lists, index 0 corresponds to the process from the worst point on the simplex, simplex[-1], index 1 to simplex[-2], etc.

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

start()

The initial n+1 simplex vertices: the start point, plus one per variable moved by its step size.

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 or, over a bounded-prior box, from the box center – its global-start mode (start_from_box, #404/ADR-0017), which also unlocks concurrent multi-start (n_starts, #498).

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.

Concurrency (#498). One Powell search is strictly serial – one line-search probe in flight at a time – so a single run uses one worker. PowellRunner factors that serial search into a headless, picklable per-start step machine so LocalMultiStartOptimizer can run n_starts of them concurrently (one worker each), the diversity a purely local method otherwise lacks. A point-start / refine fit is a single runner, byte-identical to the pre-multi-start behavior.

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: n_starts concurrent PowellRunners.

A standalone box fit (bounded priors) runs n_starts independent starts concurrently – start 0 from the box center, the rest from Latin-hypercube samples – and keeps the global best (LocalMultiStartOptimizer). A point-start (var / logvar) or refiner-injected fit is a single runner, byte-identical to the pre-multi-start behavior.

START_POINT_KEY = 'powell_start_point'

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

class pybnf.algorithms.optimizers.powell.PowellConfig(*, n_starts: int = 1, 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. Inherits the shared n_starts multi-start field (MultiStartConfig, #498): powell in box / global-start mode runs n_starts concurrent starts, keeping the global best.

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.powell.PowellRunner(u0, lower, upper, max_iterations, step, line_tol, stop_tol)

Headless, picklable per-start Powell state machine in sampling space u (#498).

Powell’s conjugate-direction search, factored out of the optimizer so a single fit can run N of them concurrently – local multi-start, the diversity a purely local 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 search directions, the _BrentLineSearch scratch, the reflecting box, the tunables – and is pure numpy / float / list: it knows nothing about PSets, the objective, backup, or dask. The orchestrator (LocalMultiStartOptimizer) drives it one evaluation at a time (Powell is serial):

(u, label), = runner.start()            # the first point to evaluate
nxt = runner.got(u_point, score)        # -> [(u, label), ...] or DONE

where each (u, label) is a next u-point to evaluate (label distinguishes the init / line / extrap phase in the submitted pset name), u_point is the realized (box-projected) u of the completed evaluation, and score its objective. got returns the next point(s) to evaluate, or the DONE sentinel when this start has converged or spent its cycle budget. 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).

got(u, score)

Consume one completed evaluation and return the next point(s) or DONE.

start()

The first point to evaluate (the start point, init phase).

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.

Restart for multimodal search (cmaes_restarts, IPOP / BIPOP; #498, ADR-0070). A single CMA-ES run descends into the one basin its start lands in, so on a multimodal objective it reaches only a local minimum. Opting into cmaes_restarts > 0 turns on the canonical multimodal-CMA-ES restart: whenever a run finishes – it converges (Hansen’s stopping battery: the search distribution shrinks below cmaes_stop_tol, the best objective stagnates, every coordinate step falls below tolerance, the covariance goes ill-conditioned, or the step size degenerates) or reaches cmaes_run_maxgen – as distinct from exhausting the global generation budget max_iterations – CMA-ES reinitializes from a fresh random point in the box with a rescaled population and keeps searching, until cmaes_restarts restarts have run or the global budget is spent. (The stagnation / ill-conditioning triggers are the #506 fix: a single principal-step test never fires on an ill-conditioned landscape, where max(d) grows as sigma shrinks, so IPOP/BIPOP would otherwise degenerate to one trapped run.) The global best is kept for free: every evaluated PSet across every restart lands in the trajectory, so trajectory.best_fit() is the best over all runs. Two schedules (cmaes_restart_strategy): IPOP grows the population geometrically each restart (population_size * cmaes_ipop_factor**k, Auger & Hansen 2005); BIPOP interleaves that increasing-population regime with a small-population regime, launching whichever has spent fewer evaluations so far (Hansen 2009), which balances broad and fine-grained search across the budget. cmaes_restarts == 0 (the default) is a single run – byte-identical to the pre-restart behavior. Restarts only draw fresh box points in box / global-start mode; a point-start / refine fit has no box to resample.

An optional cmaes_run_maxgen bounds the generations spent in each run (initial and restarted) before yielding to the next restart (#507, ADR-0085). This prevents one steadily improving local basin from monopolizing the global max_iterations budget. It defaults to unset (unbounded), preserving the existing schedule. BIPOP small runs retain their automatic evaluation-balancing cap; when the user cap is set, the smaller of the two applies.

All state is plain numpy / float / list (mean, sigma, covariance, evolution paths, the pending generation, the restart bookkeeping) – 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, cmaes_restarts: int = 0, cmaes_restart_strategy: str = 'ipop', cmaes_ipop_factor: float = 2.0, cmaes_run_maxgen: Annotated[int | None, Ge(ge=1)] = None)

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; in restart mode (cmaes_restarts > 0) it is also the tolerance for the restart battery (#506) – the relative TolFun (best-objective) stagnation threshold and the absolute TolX coordinate-step threshold. The population size lambda is the global population_size (consistent with de/pso/sa), and max_iterations remains the global generation budget across all runs. cmaes_run_maxgen optionally adds a smaller per-run generation cap; unset means no user cap. cmaes_start_point is internal (the refiner injects it), so it is not a schema field.

cmaes_restarts is the maximum number of IPOP / BIPOP restarts (#498, ADR-0070); 0 (the default) is a single run. cmaes_restart_strategy selects the restart schedule – ipop (increasing population, the default) or bipop (interleaved small / large populations). cmaes_ipop_factor is the geometric population growth per restart (2.0 = the standard population doubling), used by IPOP and by BIPOP’s large regime. cmaes_run_maxgen (#507, ADR-0085) is an optional positive per-run generation cap applied to the initial run and every restart.

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 (by PSet name), advance just that start’s runner, and return its next job(s) – [] once it terminates (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 every start’s initial job(s), 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 the accepted pset’s simulation data 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).

Reads the simdata through _result_simdata() so it works whether the objective was scored worker-side (data on res.out) or master-side (data on res.simdata); see lanl/PyBNF#480.

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, base=None)

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

Dispatches on the configured proposal operator (ADR-0067 axis 1):

• ‘whitened’ runs once its preconditioner has warmed up (else falls through to ‘de’), and

• ‘kalman’ runs during the burn-in window (_kalman_active(), else falls through to ‘de’), returning cr_idx = None since it uses no crossover;

otherwise (and always for ‘de’) the classic parallel-direction DE proposal.

base overrides the point the jump is applied to (default: the current chain state). Multi-try (ADR-0067 Stage 2) passes the selected candidate Y as base to draw the reference set from T(Y, .); with the default base=None this is byte-identical to the classic single-try proposal. (‘kalman’ is single-try only, so it never receives a base.)

Parameters:

idx – Index of PSet to update

Returns:

tuple of (PSet or None, int or None)

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.

record_accepted_output(res, index)

Cache the just-accepted state’s model output vector f(x) for chain index (ADR-0067 Stage 3), mirroring record_pointwise_loglik(): computed here, where the accepted pset’s simulation data is in hand, and consumed later at archive growth (_advance_generation()) and by the Kalman proposal.

A no-op unless the proposal declared it needs archived outputs (_archive_stores_outputs; the ‘kalman’ proposal). Reads the simdata through _result_simdata() so it works whether the objective was scored worker-side or master-side (lanl/PyBNF#480), and extracts the aligned output vector via the objective’s aligned_prediction_data() seam. Any failure is logged, never fatal – the proposal falls back to de when a chain lacks a cached output.

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', proposal: str = 'de', kalman_burnin_frac: float = 0.3, n_try: int = 1)

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).

As of ADR-0067 (Stage 1) the preconditioned proposal is no longer a separate algorithm: it is the 'whitened' proposal operator on DreamAlgorithm itself. p_dream is therefore just DREAM(ZS) with the proposal pinned to 'whitened' — expressed here as a one-line config default plus a thin subclass kept so the public job_type = p_dream code, the PDreamAlgorithm class name (re-exported by pybnf.algorithms), and the existing oracle tests all continue to resolve unchanged.

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

P-DREAM: Preconditioned DREAM — DREAM(ZS) with the 'whitened' proposal.

All behaviour lives in DreamAlgorithm, selected by proposal = 'whitened' (pinned by PDreamConfig): DE proposals are computed in a covariance-whitened space z = L_inv @ x learned online from the chain history (Adaptive-Metropolis style), so correlated posteriors are sampled in decorrelated coordinates. The proposal stays symmetric (DE differences from an external archive), so standard Metropolis-Hastings acceptance is valid.

This subclass adds no code of its own; it exists only to keep p_dream a distinct, named public job_type over the unified engine (ADR-0067).

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', proposal: str = 'whitened', kalman_burnin_frac: float = 0.3, n_try: int = 1, precondition_adapt: ~typing.Any = None)

Config for preconditioned DREAM (p_dream), co-located with the method (ADR-0006). Pins the DREAM proposal operator to 'whitened' and keeps the one preconditioning key precondition_adapt that p_dream owns (see [[test_config_unused_keys]]); everything else is inherited from DreamConfig.

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