Multiple shooting (pybnf.shooting)

pybnf.transcription (ADR-0109) is the reusable half of issue #563 — an augmented variable layout, an equality residual/Jacobian interface, and an optimizer-agnostic augmented-Lagrangian outer loop with its homotopy and its certification. It makes no simulator call and defines no job_type. This package is the other half: what it takes to state this fit as that kind of problem.

Split each scored experiment’s time course at m - 1 knots. Segment j is integrated from its own start state — segment 0’s is the model’s own initial conditions, which are often fitted parameters, and each interior knot carries an auxiliary state z_j that is searched, bounded and differentiated but is never a fit result. Continuity is imposed as c_j = Phi_j(z_j, theta) - z_{j+1} = 0 and enforced by the augmented Lagrangian, whose subproblem is solved by gntr’s own Gauss-Newton trust-region step machine. Every reported score comes from discarding z, re-simulating theta with ordinary single shooting, and scoring that, so a run that leaves continuity unconverged scores as what it actually is.

The fit type that runs all of it is job_type = ms (pybnf.algorithms.optimizers.multiple_shooting). The design, its measurements, and what this cut deliberately leaves out are recorded in ADR-0110.

Knot placement

Knot placement: where a time course is cut, and which data point lands in which piece (#563).

Multiple shooting’s first decision is where to break the horizon. SegmentGrid owns that decision and nothing else – no simulator, no objective, no optimizer – so the two things that are easy to get subtly wrong are testable in isolation: which segment a data point belongs to, and whether a coarser grid’s knots are recognisably the same knots as a finer grid’s.

Three placements, and one of them is the default for a measured reason

The issue asks for “a segment count or explicit knots; default to generic equal-time or equal-observation segments”, and all three are here:

equal_time (the default)

[start, horizon] cut into m equal spans. This is what the #563 prototype solved the motivating problem with, and it is the default because it places the knots using only the experiment’s own time axis – a fact rather than a guess. The obvious refinements that read a trajectory (knots at a burst, at a peak) are start-point dependent: they place the transcription’s structure using dynamics the fit has not established, and on the motivating problem those dynamics are exactly what is in question (an oscillator whose period is wrong everywhere except a 3 % window).

equal_observations

Cut so each segment owns the same number of measurements. It reads the data’s own time axis and nothing else, so it is not start-point dependent either – what it uses is the sampling, not the dynamics. It is the right placement when a time course is sampled unevenly, where equal spans leave some segments with nothing to fit and the auxiliary states of those segments are determined by continuity alone (ADR-0109 finding 5.2’s under-determination, arrived at through the sampling rather than through m). It needs at least two measurements per segment, so its own segment ceiling is half equal_time’s.

explicit knots

The caller supplies the times. Nothing here second-guesses them; they are validated for order and for lying strictly inside the horizon, and the run reports them.

Names, because the homotopy transfers by name

carry_over() moves an auxiliary block from one stage to the next iff the next stage declares a block of the same name and size – the layer never learns what a knot is, which is what keeps the rule generic. So the naming here is load-bearing: a knot must get the same name at every segment count that has it, or the 4 -> 2 -> 1 ladder reseeds instead of continuing and the coarsening (the mechanism, ADR-0109 finding 5.2) buys nothing.

Naming a knot by its exact fraction of the segment count does that: knot i of m is Fraction(i, m). At m = 4 the knots are 1/4, 1/2, 3/4; at m = 2, 1/2 – and Fraction(2, 4) == Fraction(1, 2), so the surviving knot carries its solved state down the ladder while 1/4 and 3/4 are discarded, which is what coarsening is. Exact rational arithmetic rather than a rounded float, so 1/3 and 0.333333 can never be two names for one knot (or one name for two).

The fraction is an ordinal, not a claim about where the knot sits. Under equal_time the two coincide, which is why the naming reads so naturally there; under the other two placements exp1@1/2 is “the knot halfway along this grid’s knot list”, which is exactly what carry-over needs to match on and is independent of what time it landed at. Every placement therefore maps the same fraction to the same knot, at every rung – which is the whole property the ladder rests on.

Which segment owns a data point

Half-open [start_j, start_{j+1}), with the final segment owning the horizon endpoint – the same convention tests/test_transcription.py::ShootingProblem pins offline. A point lying exactly on a knot is therefore read at dt = 0 from that knot’s own auxiliary state, so its prediction is the auxiliary variable itself: the one row where the data sees a segment-start state directly rather than through an integration. The alternative (a knot point belonging to the segment that ends there) reads it through the previous segment’s whole span, which is the same number only when continuity has already converged.

pybnf.shooting.grid.EQUAL_OBSERVATIONS = 'equal_observations'

Knots placed so every segment owns the same number of measurements.

pybnf.shooting.grid.EQUAL_TIME = 'equal_time'

Knots at equal spans of the experiment’s own time axis. The default, and what the #563 prototype solved the motivating problem with.

pybnf.shooting.grid.EXPLICIT = 'explicit'

Knots supplied by the caller (ms_knots), rather than derived from a rule.

pybnf.shooting.grid.KNOT = '@'

Separator between an experiment’s label and a knot’s fraction in a block name ('exp1@1/2'). A block name may not contain the layout’s '::' qualifier, and this one does not.

pybnf.shooting.grid.OBSERVATIONS_PER_SEGMENT = 2

Measurements a segment needs under EQUAL_OBSERVATIONS. Two rather than one, because a knot is placed at a measurement (that point then belongs to the later segment, read at dt = 0), so a one-observation-per-segment grid would put the last knot on the horizon itself and leave a zero-length final segment.

class pybnf.shooting.grid.SegmentGrid(times, n_segments, label='exp', start=None, horizon=None, placement='equal_time', knots=None)

The knots of one experiment’s time course at one segment count.

Parameters:

• times – The experiment’s measurement times, in the independent variable’s units. Need not be sorted; duplicates are allowed (repeat measurements at one time).

• n_segments – m, the number of segments. m = 1 is the ordinary unsegmented problem: no knots, no auxiliary variables, no continuity constraints.

• label – The experiment’s label, used to build block names. Two experiments in one fit get different labels so their knots never collide in one layout.

• start – The horizon’s start. Defaults to min(times), but a time course whose first measurement is after the model’s t = 0 still integrates from 0, so a caller that knows the simulation’s own start passes it.

• horizon – The horizon’s end. Defaults to max(times).

• placement – One of PLACEMENTS. See the module docstring.

• knots – Explicit knot times, which force placement = 'explicit'. These are the finest rung’s knots: a coarser rung of the same ladder keeps the sublist its own fractions select, so a 4 -> 2 -> 1 ladder over three explicit knots keeps the middle one at m = 2 – by the same fraction identity every other placement uses.

property block_names

One auxiliary-block name per interior knot, in segment order.

'<label>@<fraction>' – ('exp1@1/4', 'exp1@1/2', 'exp1@3/4') at m = 4. The m = 2 grid of the same experiment declares ('exp1@1/2',), which is a subset by name, so the ladder carries that knot’s solved state over.

describe()

One line for the run log.

ends

Each segment’s end time.

fractions

Each interior knot’s exact fraction of the segment count – the identity carry_over() matches on. An ordinal, not a claim about where the knot sits (see the module docstring); under EQUAL_TIME the two coincide.

knot_times

Interior knot times (m - 1 of them; empty at m = 1).

row_positions(segment)

Where each of segment segment’s data points lands in its output grid.

Returned alongside the row indices as (output rows, data row indices) so a caller reads the simulated trajectory at exactly the measured times without re-searching the grid per point.

rows_in(segment)

Indices into times of the data points segment segment owns.

sample_times(segment)

(output times, data row indices) for one segment.

The output grid is this segment’s start knot, its own data times, and its end knot, de-duplicated and sorted. The start knot is present because integration begins there; the end knot because the continuity defect is read off the same run that produced the data rows, rather than from a second simulation of the same span.

seed_times()

One output grid covering the whole horizon: every knot and every data time.

What a stage’s auxiliary variables are seeded from – one ordinary single-shoot simulation whose state is read off at each knot. Seeding this way makes the transcription feasible at iteration zero, so every discontinuity the run subsequently holds is the optimizer’s own choice rather than an artifact of the start (the prototype’s seed_aux).

segment_of

Which segment each entry of times belongs to.

starts

Each segment’s start time: the horizon’s start, then every interior knot.

pybnf.shooting.grid.max_segments(times, placement='equal_time', knots=None)

The finest segment count times can support under placement.

Above it a rung is not the method: knots fall between observations everywhere and the auxiliary states of segments with no data are determined by continuity alone. feasible_ladder() drops rungs above this and reports it.

The segment-simulation seam

The segment-simulation seam: the one place multiple shooting touches a simulator (#563).

Everything else in pybnf.shooting is arithmetic over what this seam returns – knot placement, the IC-routed objective assembly, the continuity block, the inner solver, the homotopy. Narrowing the simulator coupling to a single method is what lets all of that be verified offline against a closed-form backend, exactly as pybnf.transcription is verified against a closed-form transcription: an offline implementation of SegmentBackend exercises the same code path bngsim does.

What a segment simulation is

One span [t_j, t_{j+1}] of one experiment, integrated from a state the transcription supplies rather than from the model’s own initial conditions, with output at the segment’s own data times and at its end knot. Two properties make it different from an ordinary PyBNF simulation, and both are why this cannot ride the propose/score loop:

• the initial state is an argument, not a property of the parameter set – segment j’s start is the auxiliary variable z_j, which is never a fit result and never enters a PSet; and

• the run must come back carrying its forward sensitivities on both axes. The parameter axis gives dPhi/dtheta; the initial-condition axis gives dPhi/dz_j, which is the block the continuity Jacobian is built from and the reason #563 needed sensitivity_ic to exist at all.

The prototype’s structural finding is what keeps the objective half cheap: a segment-start state enters the data fit as an IC route with chain-rule factor 1, so the same tensor this seam already returns feeds assemble_gradient_and_fisher_hessian() with no new residual math. This module therefore returns the tensor as-is, in the native units the assembly expects, and does no differentiation of its own.

Segment 0 is not special-cased

The first segment starts from the model’s own initial conditions – which are frequently fitted (init_Z_state and friends), so they are reported free parameters rather than auxiliary ones. A caller passes initial_state = None for it, and the backend simulates the model as configured. Everything downstream reads that back as “this segment has no auxiliary block”, which is the same shape as the m = 1 stage having none at all.

class pybnf.shooting.backend.SegmentBackend

What multiple shooting needs from a simulator, and nothing else.

One instance per scored experiment (one (model, condition) pair), built once per fit and reused across every evaluation – the #563 prototype measured that constructing a sensitivity-bearing simulator costs ~17 ms warm against ~50 ms for the integration itself, so at m segments per evaluation construction would dominate.

n_simulations = 0

Segment integrations this backend has performed. The cost the #563 acceptance benchmark reports and the prototype’s paired sweeps measured multiple shooting’s 2-7x overhead in – and the number a run reports as its completed simulations, which is deliberately not the count of augmented-model evaluations: one evaluation is m integrations, so reporting evaluations would understate the cost by the very factor the method is being judged on. An implementation increments it per simulate().

abstract property nominal_state

A representative magnitude for each state, in state_names order.

Used for two things a fit cannot do without: the strictly positive EqualityModel scales – a continuity defect is a difference of states, so a model whose species span six orders of magnitude would otherwise hand the penalty term a condition number for free – and the floor under an auxiliary variable’s box.

open_lanes(pset, n_lanes)

Prepare up to n_lanes independent simulation contexts at pset; return how many are actually available (at least 1).

A lane is whatever state simulate() restarts to run one span – for the bngsim backends, a warm engine model and Simulator. Two segments cannot share one, because a segment is run by resetting that state to its own start knot, so a second segment on the same object would integrate from the first’s start. Lanes are what makes SegmentPool able to run L segments of one experiment at once.

Called on the calling thread, before any worker starts, because preparing a lane is where a backend touches the model it owns. The default implementation offers one lane, which is the serial behaviour and needs no preparation of its own.

abstractmethod simulate(pset, sample_times, initial_state=None, lane=0)

Integrate one span and return its Data.

Parameters:

• pset – The reported parameter set, exactly as an ordinary evaluation would apply it (conditions, bind-by-id, initial assignments and all).

• sample_times – Strictly increasing output times. sample_times[0] is where integration starts and sample_times[-1] is the segment’s end knot.

• initial_state – {state name: value} overriding the model’s own initial conditions, or None for the first segment, which starts from the model as configured.

• lane – Which of the contexts open_lanes() prepared to run in. 0 is the only one a serial pass ever uses, and the only one a backend that offers no lanes has to honour. A caller must never run two segments in one lane at once.

Raises SegmentSimulationFailed for a point that does not integrate.

abstract property state_names

The model state carried across a knot, in a fixed order.

For an ODE model this is its species. It is the vector an auxiliary block holds, the vector a continuity row is written on, and the axis of d_end_ic.

exception pybnf.shooting.backend.SegmentSimulationFailed

One segment did not integrate.

Raised by a backend rather than returned, because a non-integrable segment is a property of the point, not of the run: the caller converts it into a non-finite local model, the inner solver’s trust region shrinks, and the search backs off – the same way a failed simulation is handled everywhere else on the gradient path (#492). It is not a fit-ending error, which is exactly the robustness the #563 prototype measured (a segment that fails does not kill the whole trajectory: median -166.95 against single shooting’s -105.50 from an uninformed box draw).

class pybnf.shooting.backend.SegmentTrace(data, end_state, d_end_param=None, d_end_ic=None, param_axis=(), ic_axis=())

One simulated segment: its scored rows, and its end-knot state with derivatives.

Parameters:

• data – The segment’s Data over its own output grid, carrying the OutputSensitivities payload. This is the object the objective assembly consumes, unchanged.

• end_state – The model state at the segment’s end knot, in the order of state_names. The left half of the continuity defect c_j = Phi_j(z_j, theta) - z_{j+1}.

• d_end_param – d(end state)/d(native parameter), shape (n_state, len(param_axis)), or None when no parameter axis was requested.

• d_end_ic – d(end state)/d(initial state), shape (n_state, len(ic_axis)), or None. For an interior segment this is the dPhi_j/dz_j block; for segment 0 it is the derivative with respect to the model’s own initials, which is how a fitted initial condition reaches the continuity rows.

• param_axis – Native parameter ids labelling d_end_param’s columns.

• ic_axis – Species labelling d_end_ic’s columns.

pybnf.shooting.backend.trace_from_data(data, state_names, selector_prefix='species')

Read a SegmentTrace off a simulated Data.

The end-knot state is the last row of the state columns, and its derivatives are the last row of the sensitivity tensor – which is the point of putting the end knot in the segment’s own output grid rather than simulating the span twice.

A Data with no output_sensitivities yields a value-only trace, which is what the certification path (an ordinary single-shoot score, no derivatives) and a scalar evaluation need.

The bngsim SBML/Antimony segment backend (#563).

The one place pybnf.shooting touches a simulator. It implements SegmentBackend against the same private seams an ordinary SBML evaluation uses – _engine_model_for_action to build the point’s engine model, _run_simulation to integrate one span at explicit output times, _result_to_data to convert the result – so a segment is simulated by exactly the machinery a whole experiment is, with two differences: it starts at a knot instead of at t = 0, and its initial state is supplied rather than read off the model.

Why the SBML/Antimony path, and what the .net path would need

Multiple shooting is written on the model’s state: a knot carries the ODE state vector, a continuity row is a difference of states, and the continuity Jacobian is d(state)/d(state). The bngsim SBML/Antimony path hands all of that over on one run – its trajectory’s columns are the species, and its forward-sensitivity selectors are species:<name> on both axes.

A .net model is a reaction network with exactly the same kind of state, and bngsim returns both its species trajectory and its d(species)/d(species_0) when asked. What is missing is on PyBNF’s side: the net backend’s Data carries time + observables + expressions (_build_data()) and its sensitivity request names observable: / expression: selectors, so neither the state at a knot nor its derivative is in what a segment simulation would return. Closing that is an adapter change rather than a modelling obstacle – and it is not free: an experiment scores observables, so a net segment needs both selector families on one run, and on a combinatorially expanded network the auxiliary block is (m - 1) x n_species wide and the initial-condition sensitivity system is n_species wide, so the transcription’s cost scales with the expanded species count rather than with the number of fitted parameters.

This cut does not close it. job_type = ms refuses the net backend up front, naming the gap, rather than failing later at a missing selector.

Simulator reuse

The #563 prototype measured that constructing a sensitivity-bearing bngsim.Simulator costs ~230 ms cold and ~17 ms warm against ~50 ms for the integration itself, so at m segments per evaluation construction would dominate. It also verified the fix: mutating the model behind an existing Simulator and save_concentrations() + reset()-ing gives bit-identical states and sensitivities to a freshly built one. This backend therefore keeps one engine model and one Simulator per parameter point and restarts them at each knot. The restart is what distinguishes this from the pre-equilibration protocol (ADR-0052), which runs a second phase on the same simulator without a reset precisely so the state carries over.

It is also why running two segments at once needs two of them: the restart is a mutation of the simulator’s own state, so a second segment sharing it would integrate from the first’s start knot. open_lanes builds one engine+simulator pair per lane at a point, on the calling thread, and each lane is then claimed by one segment for the whole of its integration (pybnf.shooting.parallel, which also carries the cost model that decides whether extra lanes are worth their preparation).

class pybnf.shooting.bngsim_backend.BngsimSegmentBackend(model, action, mutant, suffix, timeout=None, method='ode')

One scored (model, condition) pair’s segment simulator.

Parameters:

• model – The PyBNF model object, already built and with its sensitivity request applied (enable_output_sensitivities). This backend owns it: it assigns the parameter set in place rather than deep-copying per evaluation, which is sound because the multiple-shooting driver runs on the master (a segment is not a PSet evaluation and never reaches a worker) and because the two places that write to the model – the parameter set and the action suffix – are done once per point under this backend’s lock, before any segment runs. The integration path itself only reads the model, which is what lets several segments of one point run at once on separate lanes.

• action – The TimeCourse action this experiment is measured by.

• mutant – The MutationSet (condition) it is measured under.

• suffix – The full output suffix, action.suffix + mutant.suffix.

• timeout – Per-segment wall-clock bound, from wall_time_sim. A pathological parameter point can otherwise spend minutes inside CVODE’s analytical-Jacobian failure and finite-difference retry, and a multi-start sweep is then dominated by points that were never going to score.

property nominal_state

A representative magnitude for each state, in state_names order.

Used for two things a fit cannot do without: the strictly positive EqualityModel scales – a continuity defect is a difference of states, so a model whose species span six orders of magnitude would otherwise hand the penalty term a condition number for free – and the floor under an auxiliary variable’s box.

open_lanes(pset, n_lanes)

Build up to n_lanes engine+simulator pairs at pset; return how many.

Called on the master before a parallel pass, so every lane a worker can claim is already built and no thread ever enters the construction path. A lane that cannot be built is not an error here – the pass simply runs at the width that succeeded, and the failure resurfaces as a failed segment if it was going to.

simulate(pset, sample_times, initial_state=None, lane=0)

Integrate one span and return its Data.

Parameters:

• pset – The reported parameter set, exactly as an ordinary evaluation would apply it (conditions, bind-by-id, initial assignments and all).

• sample_times – Strictly increasing output times. sample_times[0] is where integration starts and sample_times[-1] is the segment’s end knot.

• initial_state – {state name: value} overriding the model’s own initial conditions, or None for the first segment, which starts from the model as configured.

• lane – Which of the contexts open_lanes() prepared to run in. 0 is the only one a serial pass ever uses, and the only one a backend that offers no lanes has to honour. A caller must never run two segments in one lane at once.

Raises SegmentSimulationFailed for a point that does not integrate.

property state_names

The model state carried across a knot, in a fixed order.

For an ODE model this is its species. It is the vector an auxiliary block holds, the vector a continuity row is written on, and the axis of d_end_ic.

The bngsim .net segment backend – multiple shooting on a reaction network (#577).

The SBML/Antimony backend (pybnf.shooting.bngsim_backend) got the state for free: that path reports its species as the trajectory’s columns and labels its forward-sensitivity selectors species:<name> on both axes, so one run hands over everything a knot needs. The .net path does not, and for a while job_type = ms refused it on that basis – wrongly framed as a property of the model. It is not one. A .net file is a fully expanded reaction network with exactly the same kind of ODE state, and bngsim returns both its species trajectory and its d(species)/d(species_0) when asked. What was missing was that nobody asked: _build_data() assembles time + observables + expressions, and the net backend’s sensitivity request names observable: / expression: selectors.

This module asks. It is the one place in pybnf.shooting that needs both selector families on one run, and that is the whole of what makes it different from its SBML peer:

• an experiment scores observables, a continuity row is a difference of species. On the SBML path those are the same columns; here they are not. OutputSensitivities.selectors is a plain list, so one tensor carries the observable/expression rows the data terms read and the species rows the continuity block reads – requested together, off one integration, because asking twice would mean integrating twice.

• the segment’s ``Data`` carries both too, the ordinary observable/expression columns the objective scores plus the species columns trace_from_data() reads the end-knot state from. Species names carry parentheses (A(b!1).B(a!1)) and observables do not, so a collision is vanishingly unlikely – and is refused rather than silently overwritten, because a species column quietly shadowing an observable would change what the fit scores.

Nothing in the net backend is modified: this composes its existing pieces (_build_data, the engine clone, the mutant copy, the species-initializer sync) from outside, so every ordinary .net fit is untouched.

The cost, which is real and is the model’s, not the method’s

The auxiliary block is (m - 1) x n_species wide and the initial-condition sensitivity system is n_species wide, so the transcription scales with the expanded species count rather than with the number of fitted parameters. On a small network that is nothing; on a combinatorially expanded one it is the dominant term – egfr_ground.net (356 species) at m = 4 adds ~1068 auxiliary variables. That is a property of writing multiple shooting on the state of a rule-based model, not something this backend can arrange away, so job_type = ms reports the added width when it starts rather than letting a user discover it from the run time.

class pybnf.shooting.net_backend.NetSegmentBackend(model, sim_params, mutant, suffix, timeout=None)

One scored (model, condition) pair’s segment simulator, on the .net path.

Parameters:

• model – The BngsimModel. This backend owns it, assigning the parameter set in place rather than deep-copying per evaluation – sound because the multiple-shooting driver runs on the master (a segment is not a PSet evaluation and never reaches a worker), and because the writes to it happen once per point under this backend’s lock, before any segment of that point runs. A lane is one such per-point clone plus its engine and Simulator, and two segments cannot share one (see pybnf.shooting.parallel).

• sim_params – The parsed simulate() action this experiment is measured by, as _parse_simulate_action() returns it.

• mutant – The MutationSet (condition) it is measured under.

• suffix – The full output suffix, action suffix + mutant suffix.

• timeout – Per-segment wall-clock bound, from wall_time_sim.

property nominal_state

A representative magnitude for each state, in state_names order.

Used for two things a fit cannot do without: the strictly positive EqualityModel scales – a continuity defect is a difference of states, so a model whose species span six orders of magnitude would otherwise hand the penalty term a condition number for free – and the floor under an auxiliary variable’s box.

open_lanes(pset, n_lanes)

Build up to n_lanes model-copy + engine + simulator triples at pset.

The .net peer of open_lanes(), and the more expensive of the two: a lane here is a whole cloned PyBNF model, so the cost model in pybnf.shooting.parallel tilts further toward serial on a small network and further toward parallel on a large one, where the n_species-wide initial-condition sensitivity system dominates a segment.

simulate(pset, sample_times, initial_state=None, lane=0)

Integrate one span and return its Data.

Parameters:

• pset – The reported parameter set, exactly as an ordinary evaluation would apply it (conditions, bind-by-id, initial assignments and all).

• sample_times – Strictly increasing output times. sample_times[0] is where integration starts and sample_times[-1] is the segment’s end knot.

• initial_state – {state name: value} overriding the model’s own initial conditions, or None for the first segment, which starts from the model as configured.

• lane – Which of the contexts open_lanes() prepared to run in. 0 is the only one a serial pass ever uses, and the only one a backend that offers no lanes has to honour. A caller must never run two segments in one lane at once.

Raises SegmentSimulationFailed for a point that does not integrate.

property state_names

The model state carried across a knot, in a fixed order.

For an ODE model this is its species. It is the vector an auxiliary block holds, the vector a continuity row is written on, and the axis of d_end_ic.

The segment pass

Running one point’s segments: serially, or across a pool of lanes (#563).

Issue #563’s implementation proposal ends “Segment simulations can run in parallel”, and they are the one embarrassingly parallel thing in the method: one augmented-model evaluation is m spans of one trajectory integrated from m states the transcription already knows, with no data flowing between them. This module is that scheduler.

Threads, because the integration releases the GIL

The choice between threads and processes is not a matter of taste here, and it was measured rather than assumed. On the motivating model (Borghans_BiophysChem1997: 3 species, 21 sensitivity parameters, both axes requested), four warm engine+simulator replicas driven from four threads integrated 160 segments in 0.057 s against 0.153 s serially – a 2.7x speedup on 4 workers – and every trajectory column and every entry of d(y)/d(theta) came back bit-identical to the serial run. So bngsim drops the GIL inside CVODE, and the arithmetic above this seam does not have to change to accommodate a scheduler.

Processes were never a real option: a segment costs one integration, and the #563 prototype measured that constructing a sensitivity-bearing simulator costs ~230 ms cold against ~50 ms for an integration. A process pool would pay that per worker per point (or pay to pickle an engine model, which does not pickle), which is why SegmentBackend keeps warm state per parameter point in the first place.

A lane is a warm engine + simulator, and it is not free

Two segments cannot share one Simulator: the backend restarts it from the knot’s state with save_concentrations() + reset(), so a second segment on the same object would be integrating from the first’s start state. A lane is therefore an independent engine+simulator pair, and running L segments at once needs L of them at that parameter point – so the point pays L preparations instead of one.

That is the whole cost model, and it decides the default:

parallel wins when (m - 1) * t_integrate > (L - 1) * t_prepare

On Borghans, measured, t_prepare is ~4.1 ms and t_integrate is ~1-2.3 ms, so at m = L = 4 parallel segments would lose – the model is small enough that preparing the extra lanes costs more than the integrations they save. On a model where a segment is the expensive thing, it wins nearly linearly, and both cost terms move with the model rather than with the fit: the initial-condition sensitivity system is n_species wide, so a segment’s integration grows with the state while a lane’s preparation does not grow as fast.

So ms_parallel_segments defaults to 1 (serial), which is the measured right answer for the model this feature was built for, and a run that sets it reports what it measured (SegmentPool.describe()) rather than leaving a user to guess whether it helped.

What parallel gives up

The serial pass stops at the first segment that fails to integrate – the rest of that point’s segments are work whose answer is already decided, which is what keeps a search that has wandered into a non-integrable corner from paying m simulations per rejected trial. A submitted future cannot be un-run, so the parallel pass pays all m. The answer is identical; the simulation count a run reports is not, and on a multi-start sweep over an uninformed box (where most points do not integrate) that difference is the dominant cost. Stated here rather than discovered from a benchmark table.

pybnf.shooting.parallel.SERIAL = SegmentPool(lanes=1)

The pool a problem built without one uses – serial, so nothing changes for a caller that does not ask for lanes.

class pybnf.shooting.parallel.SegmentPool(n_lanes=1)

The segment pass, at one degree of parallelism.

Parameters:

n_lanes – How many segments to integrate at once. 1 (the default) runs them on the calling thread, in order, short-circuiting at the first failure – byte-for- byte the behaviour that shipped with ADR-0110.

One pool is built per fit and shared by every rung of the ladder and every start, so a thread pool is created at most once per run.

close()

Shut the thread pool down. Idempotent; a serial pool never opened one.

run(pset, tasks)

Integrate every task at pset.

Returns (traces, ok): the per-task SegmentTrace list and whether every one integrated to a finite trajectory. A pass with ok = False returns None for the traces, because the caller turns any failure into a non-finite local model and never reads the partial result – and returning a half-filled list would invite someone to.

class pybnf.shooting.parallel.SegmentTask(backend, times, initial_state, state_names)

One span to integrate: everything SegmentPool.run() needs and nothing else.

The transcription

The multiple-shooting transcription: #563’s first consumer of ADR-0109’s layer.

MultipleShootingProblem implements the two abstract methods TranscriptionProblem declares, plus certify, and gets the augmented-Lagrangian outer loop, the penalty schedule, the homotopy, the best-iterate certification, and the reporting for free. Everything specific to multiple shooting is here:

• which data point belongs to which segment (pybnf.shooting.grid);

• how a segment is simulated from a state the transcription supplies (pybnf.shooting.backend);

• how the fit’s own objective is assembled over the pieces; and

• the continuity constraints c_j = Phi_j(z_j, theta) - z_{j+1} and their Jacobian.

The objective half needs no new residual math

This is the structural finding the #563 prototype established, and it is what keeps this module small. A segment-start state enters the data fit as an IC route with chain-rule factor 1 – sensitivity_ic is dy(t)/dy0, and the transcription sets y0 directly – so an auxiliary variable is, to assemble_gradient_and_fisher_hessian(), an ordinary free parameter with an ordinary route. Each segment is presented to the assembly as its own experiment: its own simulated Data, its own slice of the observations, and its own ExperimentRouting. The assembly already sums across experiments, so the segmented data fit is the unsegmented one rearranged, and its gradient and Fisher block come out over the augmented column list in one pass.

Two things that rearrangement has to get right, and both are structural rather than cosmetic:

Segment ``j > 0`` does not read the model’s own initial conditions. They were overridden by z_j. A reported free parameter that is a fitted initial condition therefore has no effect on that segment, and its IC contribution is dropped from that segment’s routing. Keeping it would credit init_Z_state with a derivative it does not have on m - 1 of the m segments – a wrong column the fit cannot detect from its own objective. Segment 0 keeps it, which is exactly how a fitted initial condition reaches the first continuity row.

A quantity profiled or normalised across a whole series cannot be cut. An analytic per-series scale (ADR-0066) and a Data-level normalization (ADR-0053/0102) are functions of the series they are computed over, so splitting the series changes them; a cumulative-to-incident transform (ADR-0051) is a difference between neighbouring rows, and the row before a knot is in another piece. Those are refused by the ms fit type’s gates rather than silently rearranged. An analytically profiled noise scale (ADR-0108) is the opposite case and needs no gate: it is profiled over the pooled residuals of every supplied experiment, so cutting one series into m pieces pools exactly the same residuals and gives exactly the same sigma_hat.

The constraint terms never enter f

objective_value is the fit’s own objective, scored by the fit’s own objective function on the segmented trajectory. The penalty lives strictly outside it (ADR-0109), which is what keeps an estimated noise scale from absorbing continuity violation as measurement noise – and 13 of the 23 slugs in the motivating benchmark corpus estimate at least one.

pybnf.shooting.problem.AUX_DECADES = 6.0

Half-width, in decades, of an auxiliary state’s box around its own magnitude. Wide enough that the box is not a constraint on the search (the #563 prototype used a fixed 1e-6 .. 1e3 window and never reported a knot pinned at a bound), finite because the inner optimizers this layer feeds are bound-constrained and a segment-start concentration that is allowed to go negative is not a state the simulator can restart from.

class pybnf.shooting.problem.MultipleShootingProblem(experiments, objective, variables, pset_from_u, blocks, scales, name=None, pool=None)

One rung of the segment ladder: the fit, transcribed at a fixed segment count.

Parameters:

• experiments – The SegmentedExperiments, all at the same segment count.

• objective – The fit’s own objective function.

• variables – The fit’s reported free parameters, in Configuration.variables order – the same order every existing PyBNF seam uses.

• pset_from_u – (u_reported) -> PSet: the algorithm’s own sampling-space to parameter-set bridge (_pset_from_u()), so this module never re-derives the theta <-> u transform.

• blocks – The auxiliary VariableBlocks, one per knot of every experiment, already seeded (see seed_stage()).

• scales – {experiment key: per-state constraint scale} – the magnitude each continuity defect is measured against.

• name – The stage label for the trace ('m=4').

• pool – The SegmentPool that runs this stage’s segment passes. Defaults to the serial one, so a caller that does not ask for lanes gets exactly the behaviour that shipped with ADR-0110.

certify(reported)

Reconstruct reported through the fit’s ordinary unsegmented path.

Every auxiliary state is discarded and each experiment is simulated once over its whole horizon from the model’s own initial conditions – which is the single-shoot problem – and scored by the fit’s own objective. That score is the only one comparable with an ordinary PyBNF fit’s, and under ADR-0109 it is the only one this run may report: the augmented objective at an infeasible point is computed on trajectories that do not join up.

property constraint_names

One label per constraint, in residual order.

describe()

The stage in one line: what it added, and where.

equality_at(u)

The EqualityModel at u.

property layout

This stage’s AugmentedLayout.

property name

A short label for the stage trace (multiple shooting uses 'm=4').

objective_at(u)

The fit’s objective and its derivatives at u, memoised on the point.

The memo is not an optimisation first — it is what keeps this method scoring each simulated trajectory exactly once (#578). Scoring goes through objective.evaluate_multiple, which asks the measurement layer to materialise every observable: <id>, formula: ... column into the trajectory in place (ADR-0036), and that materialisation deliberately refuses a column that already exists. Every ordinary fit satisfies it for free: the propose/score loop scores a freshly simulated Data each time. Multiple shooting does not, because _traces() caches the segment trajectories per point so that one augmented-model evaluation costs one pass of segment simulations rather than two – and the outer loop then re-evaluates at the point the inner solver finished at, which is a cache hit on the very same objects.

Memoising the assembled model rather than copying the trajectories fixes that at the cause: the second call returns the first call’s answer instead of re-scoring anything. It also removes a redundant gradient/Fisher assembly per outer iteration, which is the larger of the two costs. Sound because this method does not depend on the multipliers – only AugmentedModel combines them – so the objective at a point is a property of the point alone.

pybnf.shooting.problem.STATE_FLOOR = 1e-12

Floor under a state magnitude, so a species that is identically zero over the whole horizon still gets a strictly positive constraint scale and a finite auxiliary box. The layer requires strictly positive scales precisely so this decision is made once, here, rather than being discovered as a divide-by-zero inside the penalty term.

class pybnf.shooting.problem.SegmentedExperiment(spec, grid)

One ShootingExperiment cut at the knots of one stage.

certification_grid

The m = 1 grid of the same experiment – what certify() simulates. Not a special case: the unsegmented problem is the one-segment transcription.

segment_exp_data(segment)

This experiment’s observations restricted to one segment’s rows.

Cached: the slice depends only on the grid, never on the fit point, so it is built once per stage rather than once per evaluation.

segment_key(segment)

The dictionary key one segment is scored under.

A synthetic suffix, distinct per segment, because the objective scores a {model: {suffix: Data}} mapping and the m pieces of one experiment must be m entries rather than one that overwrites the others.

class pybnf.shooting.problem.ShootingExperiment(key, backend, exp_data, routing, times=None, label=None, start=None, horizon=None, placement='equal_time', knots=None)

One scored (model, condition) pair, at the fit level: everything about it that does not depend on the segment count.

Built once per fit and shared by every rung of the ladder, which is what makes a stage cheap to construct: a homotopy builds one SegmentedExperiment per rung around the same backend, observations, and routing.

Parameters:

• key – (model name, suffix) – the identity this experiment has in the simulated/experimental data dictionaries the objective consumes.

• backend – Its SegmentBackend.

• exp_data – Its observations, as an ordinary PyBNF Data.

• routing – The experiment’s prebuilt ExperimentRouting over the reported free parameters, exactly as the gradient optimizers build it once per fit.

• times – Its measurement times. Defaults to exp_data’s independent-variable column.

• label – The label its knots are named under. Defaults to the suffix.

• start – Where the simulation starts. Defaults to the first measurement time; a time course whose first measurement is after t = 0 should pass 0.0.

• placement – How this experiment’s knots are placed – one of PLACEMENTS. Fixed for the whole fit, so every rung of the ladder places its knots by the same rule and their fractions mean one thing.

• knots – Explicit knot times for the finest rung, which force placement = 'explicit'.

grid(n_segments)

This experiment’s knots at one segment count.

property max_segments

The finest rung this experiment’s data supports under its own placement.

pybnf.shooting.problem.seed_stage(specs, n_segments, objective, variables, pset_from_u, reported, aux_decades=6.0, name=None, pool=None)

Build one rung of the ladder, with its auxiliary states seeded from reported.

Each experiment is simulated once, unsegmented, at the incoming parameters, and the state is read off at each knot. Seeding this way makes the transcription feasible at iteration zero: every continuity defect is exactly zero at the start point, so any discontinuity the run subsequently holds is the optimizer’s own choice rather than an artifact of how the stage was built. It is also why a homotopy stage is a callable in run_homotopy() – the seeds are not knowable until the previous stage has finished moving theta.

A start point that does not simulate is not a failure: the knots fall back to the model’s own declared state magnitudes, the first inner solve sees a large defect, and the run proceeds (or stops with a stated reason) rather than dying at construction.

aux_decades sets each auxiliary variable’s box: +/- aux_decades around the state’s own magnitude, in log10. Wide enough not to constrain the search, finite because the inner optimizers are bound-constrained and a segment-start concentration that is allowed to go negative is not a state a simulator can restart from.

The inner solver

The inner solver: gntr’s step machine, driven against an augmented subproblem (#563).

ADR-0109’s layer ships no inner solver, deliberately – shipping one would state a preference the optimizer-agnostic contract exists to avoid. A consumer has to choose one, and for multiple shooting the choice is measured rather than assumed. Over 30 data seeds x 2 starts on the layer’s offline shooting problem (120 runs), a trust-region least-squares solver converged 60/60 while a quasi-Newton one converged 36/60 and stalled out on the rest. The reason is structural: the KKT stop needs the scaled defect and the first-order optimality below tolerance in one iterate, and a method built from gradient differences handles an augmented Lagrangian whose penalty term carries a large rho far less well than one that sees rho J_c^T J_c explicitly. So the MVP steps from the Gauss-Newton form, and a consumer should not treat the inner optimizer as free.

Reusing the fit type rather than reimplementing it

The Gauss-Newton form is (gradient, PSD hessian) – exactly what job_type = gntr already consumes, and _GNTRRunner is already a headless, backend-free step machine over that pair: ridge-regularise H, eigen-factor it into the pseudo-Jacobian that reproduces the Coleman-Li-scaled Newton step, and run trf’s bound-constrained trust-region-reflective accept/reject state machine unchanged. So this module is a driver, not a method: it feeds that runner the augmented Lagrangian instead of the fit’s own objective, and nothing about the step math is new or separately tuned.

Driving it synchronously is the one difference from gntr’s own use, and it is forced by the interface above rather than chosen: ADR-0109’s contract is solve(subproblem, u0, tolerance) -> InnerOutcome, a blocking call, because an inner solver “never calls back into the outer loop”. Since a segment simulation is not a PSet evaluation either (pybnf.shooting.backend), nothing is lost: the propose/score loop was never available to this path.

Two tolerances, and why the outer loop’s is only a floor away from being obeyed

tolerance is omega_k, the outer loop’s inner-optimality target: loose at first and tightening as the penalty rises. It becomes the runner’s grad_tol, so an early subproblem is solved roughly and a late one tightly – which is the entire economic argument for the augmented-Lagrangian frame. It is floored (GaussNewtonSolver.grad_tol_floor) because omega decays geometrically and will eventually pass below what any solver can demonstrate on a finite-precision Hessian; past that point the runner would simply spend its whole iteration budget every outer iteration, which is the same waste ADR-0109 finding 5.1 measured on a too-loose penalty, arrived at from the other side.

class pybnf.shooting.solver.GaussNewtonSolver(max_iterations=50, ridge=1e-10, step_tol=1e-10, grad_tol_floor=1e-10, stop_check=None)

An inner solver on ADR-0109’s contract, stepping from the Gauss-Newton form.

Parameters:

• max_iterations – Trust-region iterations per inner solve. Bounded per outer iteration rather than per run: an approximate inner minimisation is what the augmented-Lagrangian method is designed around, and the outer loop’s stall detector is what notices a solver that stops achieving anything.

• ridge – The relative Levenberg ridge added to the curvature before the pseudo-Jacobian factorisation, as in gntr_ridge. It matters more here than in an ordinary fit: with one observed state of three, the auxiliary states of the unobserved two carry no data term at all and their data-fit curvature is exactly zero, so the constraint block rho J_c^T J_c is the only curvature they have.

• step_tol – Negligible-step tolerance, as in gntr_step_tol.

• grad_tol_floor – Floor under the outer loop’s omega_k (see the module docstring).

• stop_check – Zero-argument callable; True truncates the inner solve at whichever iterate it has reached. The wall-clock-budget seam (ADR-0093/0107) – a truncated inner solve is a normal outcome here, not a failure, and the outer loop keeps the iterate.

n_evaluations

Model evaluations spent across every inner solve this object has driven – the cost accounting the #563 acceptance benchmark reports, and the quantity the prototype’s paired sweep measured multiple shooting’s 2-7x overhead in.

The ladder

Driving the ladder: one fit, transcribed at several segment counts in turn (#563).

The whole of multiple shooting’s outer control flow is run_homotopy(); this module supplies the three things it needs that are specific to segments – the ladder of segment counts, a stage factory per rung, and the defaults ADR-0109’s findings fixed.

Why the stages are callables

run_homotopy accepts a stage as a TranscriptionProblem or as a (reported) -> TranscriptionProblem callable, and multiple shooting needs the second form: a stage’s auxiliary variables are seeded by reading a nominal trajectory at the incoming parameters at each knot, and the incoming parameters are whatever the previous rung finished at. Seeding that way makes each stage feasible at its own iteration zero, so every discontinuity a run holds is the optimizer’s choice rather than an artifact of how the stage was built.

The ladder starts in the middle

(4, 2, 1), from coarsening_ladder(), and the middle is the measured part. Starting at many short segments – the easiest landscape, and what the original formulation proposed – is the wrong end on the motivating problem: with one observed state of three and ~14 points per segment, m = 8 is under-determined, the data term is satisfiable without correct dynamics, and the stage routinely certifies worse than its own start. Over eight paired starts, 4-2-1 had the best tail at moderate cost and 8-4-2-1 did not. The ladder always ends at m = 1, which is the ordinary unsegmented fit – so a multiple-shooting run finishes by solving the problem the user actually asked about, and its last rung is the one whose objective needs no certification because it is the certificate.

pybnf.shooting.driver.coarsening_stages(specs, rungs, objective, variables, pset_from_u, aux_decades=6.0, pool=None)

One stage factory per rung, in the order run_homotopy will run them.

pybnf.shooting.driver.feasible_ladder(specs, ladder=None)

The requested ladder, clamped to what the data can support.

A segment count above what an experiment’s data supports places knots between observations everywhere and leaves segments with nothing to fit – a transcription whose auxiliary states are determined by continuity alone in every segment, which is not the method, it is an ODE solver with extra variables. Rungs above that are dropped and the caller is told which (a silent cap would read as “we ran the ladder you asked for”).

The ceiling is the experiment’s own, because it depends on the knot placement: equal time needs one measurement per segment, equal observations needs two, and an explicit knot list is the finest rung (max_segments()).

Returns (rungs, dropped).

pybnf.shooting.driver.run_multiple_shooting(specs, objective, variables, pset_from_u, reported_start, ladder=None, schedule=None, inner_solver=None, max_outer=25, aux_decades=6.0, stop_check=None, on_iterate=None, on_stage=None, pool=None)

Solve one fit by multiple shooting, down the coarsening ladder.

Parameters:

• specs – The ShootingExperiments – one per scored (model, condition) pair, built once for the whole fit.

• objective – The fit’s own objective function.

• variables – Its reported free parameters, in Configuration.variables order.

• pset_from_u – The algorithm’s sampling-space to PSet bridge.

• reported_start – The fit’s start point in sampling space.

• ladder – Segment counts, finest first. Defaults to (4, 2, 1).

• schedule – The PenaltySchedule. Its defaults carry ADR-0109 finding 5.1 – the penalty starts tight (rho_0 = 10, gamma = 5), because on the motivating problem a loose start was both worse and twice as expensive: the inner solve on a nearly-unconstrained subproblem never converges and burns its whole budget every outer iteration.

• inner_solver – Defaults to a GaussNewtonSolver.

• stop_check – The wall-clock-budget seam, passed to the outer loop and to the default inner solver, so a deadline lands inside a long inner solve rather than after it.

• pool – The SegmentPool every rung’s segment passes run through. Defaults to the serial one.

Returns the HomotopyResult, whose best is the best certified iterate over the whole ladder – not its last, which on one prototype start held -147.0 while an earlier iterate certified at -196.3.