Repository navigation
feat: job_type = lwf, the Liu–West filter - #985
Open
elyfmiller wants to merge 8 commits into
Open
elyfmiller wants to merge 8 commits into
elyfmiller wants to merge 8 commits into
Conversation
…e math constant _build_safe_eval_namespace laid the math names over the model's parameter values, on the premise that BNG2.pl reserves them all. It does not: BNG2.pl 2.9.2 refuses 15 of the namespace's 22 names as parameter names (exp, log10, log2, sqrt, abs, sin, cos, tan, asin, acos, atan, min, max, if, rint) and accepts the other 7 (e, pi, log, atan2, ceil, floor, pow); its own constants are _e and _pi. A model that declares e 4 and seeds B() e starts from B = 4 under BNG2.pl and under bngsim alone, but the lanl#450 sync in BngsimModel.execute() reseeded B with math.e. The network-free path's _evaluate_bngl_params kept e out of its namespace the same way, so [e = 4, k = 2*e] gave k = 2 math.e instead of 8. One set, _BNG_RESERVED_EVAL_NAMES, now holds the names BNG2.pl refuses. A model parameter takes precedence over every other math name, both in _build_safe_eval_namespace (the lanl#450 sync, _eval_model_expression and _eval_numeric) and in _evaluate_bngl_params. A reserved name still means the built-in whatever the seed holds. test_evaluate_bngl_params_param_named_like_builtin_does_not_shadow pinned x = math.e for [e = 5, x = e]; it now pins the same rule with a name BNG2.pl does reserve (sqrt). New tests read each of the seven names as the parameter, evaluate a parameter block as BNG2.pl does, and, on bngsim, check execute() for e and for pi against the seed by hand and against bngsim alone, with and without a setConcentration expression. Signed-off-by: Ely F. Miller <efm46@nau.edu>
… the expected count 0 Under location = mean, NegBinomial.data_fit took the prediction as the mean in prob = r/(r + mean) and clipped prob to [1e-10, 1 - 1e-10]. A negative prediction takes prob out of (0, 1]. In (-r, 0) the clip happened to give the score of a zero prediction; at exactly -r it divided by zero; below -r prob turned negative and was clipped to the other end, so an observed 0 scored -r ln(1e-10), 46 at r = 2. A cumulative output that falls by more than r over an interval gives such a prediction, as does any model output below zero. _mean now clamps a MEAN prediction at 0, as _mean_for_median already clamps the median's target. No score of a prediction above -r changes, and neither does any derivative: _mean_for_slope floored the mean above 0 before and still does. The tests compare data_fit with scipy.stats.nbinom at the clipped prob, for counts of 0 and 3: predictions of -1e6, -r and -1 score as the mean 0, and predictions of -1, 0 and 4 score as they did. Signed-off-by: Ely F. Miller <efm46@nau.edu>
… cluster and -r Three hooks in main() for a method that runs in the PyBNF process and keeps its own state between runs. They have no user until the Liu–West filter commit later in this series: no existing job type sets them, and nothing changes for an existing job. - pybnf/run_directories.py names the folders and files main() deletes before every run that is not a --resume. _prepare_run_directories reads them from there, so a method can refuse to keep a file in a place main() will delete. - Algorithm.needs_cluster, True by default. main() starts the dask cluster, and hands run() its client, only when it is set. - Algorithm.resume_refusal, None by default: a (message, hint) pair with which _resolve_continue_file refuses -r, and with which main() no longer offers to continue an in-progress pickled run. The -t and -s flags are now written into the configuration before a fresh run's algorithm is built (_apply_cluster_flags), and right after a resumed one is loaded, so a method that cannot honour them can refuse them before any cluster starts. Every existing job type reads them when it starts its cluster, as before. Tests: _apply_cluster_flags writes each given flag over its key and nothing else, and _prepare_run_directories deletes exactly the names run_directories lists, each spelled out in the test too, so a name dropped from the shared lists is caught. Signed-off-by: Ely F. Miller <efm46@nau.edu>
…es state (lanl#973) A sequential method that carries each particle's species state from one observation to the next needs one primitive: set a state, integrate one interval under one parameter set, and read back the new state and the outputs. The Liu–West filter later in this series and the IBIS sampler of lanl#973 both need it, and lanl#973 names its property: a chain of segments and one full simulation over the same times agree to solver tolerance. pybnf/bngsim_model/segment.py adds SegmentIntegrator(model, parameter_names, *, suffix=None, timeout=None), with initial_state(values) and integrate(state, values, t0, t1, times=None) -> Segment(state, data), and SegmentFailed. It is built from BngsimModel's own pieces, as NetSegmentBackend is (lanl#577): the pickling hooks for an engine with the codegen artifact attached, _codegen_kwargs, execute()'s initial-state preamble with the lanl#450 sync, the simulate action's tolerances, wall_time_sim and _build_data. net_model.py is unchanged. One engine is reused across calls, and each call writes every parameter and then the whole state, in that order, since writing a parameter re-derives the species it seeds. An initial state is computed on a fresh clone, since the seed sync pins a reused engine's reset(); that costs 0.05 to 0.15 ms on 2- to 4-species networks and 2.7 ms on the 356-species EGFR network. A call made while another runs is refused: eight threads sharing one integrator returned wrong chains and, when reproduced, crashed the interpreter. A copy or an unpickled integrator has engines of its own. At construction it refuses, by name, what it cannot integrate the way the ordinary path simulates: a model that is not a bngsim .net model; a simulate action that is not ODE or sets continue, stop_if or steady_state; an experiment that does not start from the seed state; a model with conditions; and a parameter that is unknown, a function's slot, or defined by an expression (a reused engine would carry the bngsim lanl#188 override between parameter sets). Edition 2's later experiments on one model, which reset parameters and species first, are accepted. The oracles do not go through the integrator: the closed form of A -> B, also with its rate changed between segments and against the ordinary path's setParameter and continue=>1 run; a chain against one full simulation, the integrator's own and execute()'s, within 100 units of rtol |y| + atol at rtol = 1e-6 and 1e-10; and the first segment against execute(), byte for byte. Solver tolerance means the two runs' global error, not a fixed multiple of rtol, since every segment restarts the solver: on the Brusselator over [0, 100] at bngsim's defaults, a chain of 10 segments differs from one run by 141 units and a chain of 1,000 by 4,471, while at 10 segments both are equally far from a reference integrated at rtol = 1e-13. Signed-off-by: Ely F. Miller <efm46@nau.edu>
…ly through its initial state A method that carries each particle's species state gains nothing by moving an id that only sets the initial state: the carried state already holds its effect, and a moved value would only part the reported value from the state it produced. The Liu–West filter, next in this series, leaves such ids unmoved. expressions.initial_state_only_ids(net_lines, ids, exclude=()) returns, in order, each id that appears in a species initializer, itself or through a parameter derived from it at any depth, and on no line outside the parameters and species blocks. A block it does not know counts as a later use, so an unknown use can only keep an id moving. Numbers are read before names, so the e of 1e-9 is not a name. exclude names ids never to report, such as a noise parameter the likelihood reads at every row. It reads the network with the existing _net_blocks and _net_param_definitions. The unit tests use hand-written .net text: seeds through derived parameters, an id that seeds and also reaches a rate law, uses in a function, a rate and an unknown block, exclude, and a number with an exponent. The oracle is the engine: from a carried state, one segment under two values of each id is byte-identical exactly for the ids reported. Signed-off-by: Ely F. Miller <efm46@nau.edu>
job_type = lwf forecasts count data that arrive one interval at a time from a mechanistic model. Each particle carries its free parameters, its own model state and a weight. At each data row the Liu and West (2001) kernel moves the free parameters, x <- a x + (1 - a) m + h L z with a = sqrt(1 - h^2), in the logit of each box (lwf_bounds = logit) or in sampling space folded back into the box (reflect); SegmentIntegrator integrates each particle's state over the interval; the job's own cumulative neg_bin line weighs each particle on the interval's increment; and the population is resampled systematically when the Weight ESS, 1 / sum(w^2), falls below lwf_resample_threshold times the number of particles. After the last row it forecasts lwf_forecast_intervals intervals. A free parameter that acts only through the initial state is never moved. The move is a transition of the parameters, so the output under Results/LWF/ is a forecasting sample of a model whose free parameters drift by the kernel, not the posterior of fixed parameters, and every output file and the run summary say so. The filter integrates in the PyBNF process, as hmc and ms do, so it sets needs_cluster to False; lwf_independent_runs = N runs N Independent Runs, each from its own prior draw, as tasks on PyBNF's cluster, and stacks their files. Every random stream is keyed on the run, its purpose and the row's time, so with lwf_continue = 1 each run assimilates only the rows after those its state file holds and writes what one uninterrupted run writes. A state file carries a fingerprint of everything the arithmetic depends on, and a continuation under a configuration that differs, or over rows that were revised, is refused naming what changed. main() refuses -r for lwf through resume_refusal, naming lwf_continue. wall_time_fit stops a run between rows and leaves a state file to continue from. The filter needs edition 2, one bngsim .net model and one experiment with one count column, its rows one step apart and the first one step after t = 0. The column is declared by a neg_bin line with location = mean, a fitted dispersion and the cumulative flag, and a whole-fit line must say the same. Every key, declaration or value it would not honour is refused by name when the configuration loads. The oracles do not go through the filter's own arithmetic. With no move and no resample, the weights are the likelihood the fit path's objective gives an ordinary execute() of each particle; each weight is scipy.stats.nbinom of the particle's own increment; on A -> B each increment, forecast intervals included, is the closed form A0 (exp(-k t0) - exp(-k t1)); the initial population passes a Kolmogorov–Smirnov test against each prior; and the forecast counts have the negative-binomial mean and variance. The kernel keeps each weighted mean and variance and scales each covariance by a_i a_j + h_i h_j, exactly, on draws built for it; systematic resampling gives each particle floor(N w) or ceil(N w) copies; the fold into the box is FreeParameter.set_value's. A continuation, and a wall-time stop followed by a continuation, write byte for byte what the uninterrupted run writes, and several runs on a dask LocalCluster write byte for byte what the same runs write in sequence in this process. The refusal tests check that each message names its cause. Signed-off-by: Ely F. Miller <efm46@nau.edu>
docs/algorithms.rst gains a Liu–West filter section: the four steps of an update with their formulas; why the output is a forecasting sample of a model whose free parameters drift and not a posterior; the output files; Independent Runs and the cluster; the random streams, the state file and continuation; the wall-time budget; what the filter refuses; and its references (Gordon et al. 1993, Kong et al. 1994, Kitagawa 1996, Carpenter et al. 1999, Liu and West 2001, Chopin 2002). docs/config_keys.rst documents job_type = lwf and each lwf_* key with its default and an example, and says how wall_time_fit stops the filter. The docs index, the module reference and the CONTEXT.md glossary name the filter and its family, filter. Signed-off-by: Ely F. Miller <efm46@nau.edu>
examples/tutorial/50_liu_west_filter forecasts daily counts from an SIR outbreak with job_type = lwf. make_data.py draws 42 negative-binomial counts around the model's daily increases, integrating the model by hand so it needs only numpy; the README runs the filter, reads its files and continues the run with lwf_continue = 1 when a new row arrives. Lesson 50 joins the tutorial table and docs/tutorial.rst under Forecasting. tests/test_tutorial_liu_west.py checks make_data.py's means against an ordinary bngsim simulation of the lesson's model, that the committed data file is the one the script writes, that a short run of the lesson through pybnf writes the files the README lists, and that the README's continuation writes byte for byte what an uninterrupted run writes. It joins the pre-push bngsim test list. Signed-off-by: Ely F. Miller <efm46@nau.edu>
elyfmiller
added a commit
to elyfmiller/flubnf
that referenced
this pull request
Oct 13, 2026
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This adds
job_type = lwf, a Liu–West filter that forecasts count data arriving one interval at a time from a mechanistic.bnglmodel. It keeps a weighted population of particles, each with its own free-parameter values and model state, and assimilates each new row instead of refitting. Its output is a forecasting sample of a model whose free parameters drift by the kernel, not a posterior. It runs on the segment integrator #973 asks for, which the exact sequential sampler planned there can reuse.No ADR is proposed: the design decisions are recorded in the commit messages and in the three documentation subsections (the algorithm, the forecasting-sample stance, what it refuses).
What changed
Commits 1–3 change PyBNF outside the filter. Commits 1 and 2 are separable; commit 3 is what the filter needs from
main().fix(bngsim). A model parameter namede,pi,log,atan2,ceil,floororpowis read as the parameter, as BNG2.pl and bngsim read it, wherever PyBNF evaluates a model expression. This can change an existing fit of such a model.fix(noise).neg_binatlocation = meanscores a negative prediction as the expected count 0. Before, a prediction below −r was scored wrongly, and one of exactly −r divided by zero. Scores of predictions above −r are unchanged.feat(core).pybnf/run_directories.pynames whatmain()deletes, and two hooks,Algorithm.needs_clusterandAlgorithm.resume_refusal, default to today's behaviour.-t/-sare written into the configuration before the algorithm is built.feat(bngsim):SegmentIntegrator(Add an IBIS sampler (job_type = ibis) that updates the exact posterior as observations arrive, absorbs revised observations, and finishes within wall_time_fit #973). Integrates one interval from a carried species state. A chain of segments gives one full simulation's outputs within integration error. It holds bngsim's logger at warning level while a segment integrates.feat(bngsim):initial_state_only_ids. Finds the free parameters that act on a.netmodel only through its initial state, conservatively. The filter never moves them.feat(filter):job_type = lwf, itslwf_*keys, a new registry familyfilter, and its tests.docs(filter). The algorithm section, the key reference and the glossary entry; the algorithm-development guide gains theanalysisandfilterfamilies and the two hooks, and the troubleshooting page a note on a budget-stopped run.docs(filter). Tutorial lesson 50, which runs in CI.What the filter is and is not
neg_binline on the increment, the same term an ordinary fit scores; the population is resampled systematically when the Weight ESS falls below a threshold. After the last row it forecastslwf_forecast_intervalssteps.lwf_continue = 1assimilates only new rows and writes what an uninterrupted run writes, byte for byte.am; an exact sequential sampler is Add an IBIS sampler (job_type = ibis) that updates the exact posterior as observations arrive, absorbs revised observations, and finishes within wall_time_fit #973's subject.lwf_independent_runsis neithern_starts(ADR-0071) nor the gradient methods'population_size: those keep the best of competing searches, which is why ADR-0071 scopesn_startsto the methods that do; independent runs are combined, not selected.Errors, warnings and shared paths
am: same objective, no distributional comparison, since the targets differ. Several Independent Runs in sequence and on a local cluster: byte-identical. A continued run: byte-identical to the uninterrupted one, in process, on a cluster and throughmain(). Multiple shooting and gradients: not used; a segment refuses an experiment whose start state an earlier action moved, where multiple shooting replays none of the list (ADR-0151, job_type = ms on a .net model integrates only the scored simulate() line and skips the rest of the action protocol (pre-equilibration, setParameter/setConcentration, carried state), so it fits and certifies the wrong experiment #910). The shared run loop and its_finalize_runtail: not used (design choice 6). BNG2.pl, network-free, SBML and RoadRunner: refused.Tests
The PR adds 265 tests to the CI selection (6,444 on
mainat b5fa894, 6,709 here).Oracles, none going through the code under test: weights equal an ordinary fit's likelihood (no move, no resample) and scipy's negative binomial; increments match a first-order conversion's closed form; the kernel keeps each mean and variance and scales each covariance by a_i a_j + h_i h_j, checked by hand; the prior draw passes a KS test; systematic resampling gives ⌊Nw⌋ or ⌈Nw⌉ copies; forecast counts have the negative binomial's mean and variance; a segment chain matches a closed form and one full simulation, and the first segment equals
execute()byte for byte; continuation and cluster runs are byte-identical to uninterrupted sequential runs.Run on macOS arm64 (Accelerate), numpy 2.5.3 (2.5.4 on the 3.13 and 3.14 legs), scipy 1.18.1, BioNetGen 2.9.2, with the dependency set CI's setup action installs for each leg:
pytest -m "not slow and not recovery": 6,691 passed, 12 skipped, 6 xfailed, 0 failed.PYBNF_NO_BNGSIM=1as CI's third leg runs: 5,903 passed, 804 skipped, 2 xfailed, 0 failed.ruff0.15.14,tools/changelog.py check, andsphinx-build -Wwithout bngsim as the docs job runs it: passed.Not run: Linux (OpenBLAS) and Windows.
Left out
v1 refuses by name: anything but one
.bnglmodel on the bngsim.netpath withmethod: ode; more than one experiment, data file or column; conditions, mutants, pre-equilibration, constraints, measurement models,normalization,time_error; any observation model but the cumulativeneg_binatlocation = meanwith a fitted or fixed dispersion; rows not one step apart (a missing row is ananrow), and a count below 0 or not finite; on a continuation, a revised assimilated row;samples.txt,-r,-t/-sand the cluster keys. Best-fit outputs andoutput_inference_dataare not written, and their keys are reported as unused. Choosing h from the data is not offered. Hand-edited values in a state file are not checked; its format, fingerprint, shapes, finiteness and rows are. A count below 0 is refused rather than excluded with a warning, which is what ADR-0090 has an ordinary fit do with the same value and what the filter does with ananrow: in a series of increments a negative row means an earlier cumulative total was revised, which an assimilated state file cannot take back, so the exclusion is the user's explicit choice, a row written asnan, rather than the filter's. NoResults/<model>_bestfit.bnglis written: ADR-0048's artifact is the recorded minimum-objective point, and a weighted population has none to label, sosave_best_dataandembed_best_fit_dataare reported as unused, as ADR-0014 has every job type do for a key it does not read.Noticed, not changed
All three predate this PR; the first two are not reachable from an
lwfconfiguration.objective.pymatches a simulated row to a data row withnp.iscloseat its default relative tolerance, so fromt / step >= 1e5on it picks the wrong row and scores the wrong interval, for any fit. The filter would reach it at row 100,000._sync_species_initial_concentrationsreads^in a hand-written.netspecies initializer as Python's XOR (2^3gives 1), silently. BNG2.pl hoists such expressions into parameters, so a.bnglmodel never meets it; only the Python API with a.netdoes.main()deletesResults/before the algorithm is built, so a run refused at load has already lost the previous run's outputs (state files are untouched). This applies to every job type.Design choices
Each follows what PyBNF already does.
lwfregisters in a family of its own,filter, asdesigndid withanalysiswhen it fit neither optimizers nor samplers: the family is what keeps a method out of the benchmark harness and PEtabjob_type = all, and a forecasting sample is not a posterior.Results/stop_reason.txtis the only difference, andmethod_chain.jsonrecords the phase aswall_time_expired. The budget is checked between rows, so the overrun is at most one row oflwf_particlessegments, each bounded bywall_time_sim.lwf_continue = 1finishes the job from the end of the data.main()records the fit phase afterrun()returns, such a run writesResults/LWF/but nomethod_chain.json, as any run that ends in an error does (ADR-0107).dispersion = fit <r>anddispersion = fix_at <value>are both accepted, as the edition-2 noise grammar offers them; the whole-fit line and the column's must agree. Underfix_atthe kernel moves the model's free parameters only, and the weights and the forecast draws use the fixed value.samples.txt; no run is weighted by its evidence, which nothing in PyBNF does.lwfoverridesAlgorithm.run, ashmcdoes, and likehmcdoes not call_finalize_run, which would report a missing best fit for a run that finished; it uses the shared stop-reason and Simulations-folder tails (ADR-0093, ADR-0110).refineandbootstrapother than 0, whichmain()would act on, and the cluster keys, since Independent Runs use the local cluster, are refused. A key no schema owns is the shared check's._bngsim_capsscopes its own handler.Interactions under
job_type = lwf.bnglmodel on the bngsim.netpath.target;bngl_backend = bionetgen; network-free.bnglwhose actions are not only its network definitionexperiment:lines;condition:,mutant,preequilibrate:,equil_t_end:type:other thantime_course; first data column nottime; a time ofinfor not after 0; a count below 0 or not finitemethod:other thanode; severaldata:files;.con/.propconstraints;t_start:,t_end:,n_steps:observable: <entity>, column: <header>measurement_params:, formula observables,gml:,complex:noise_model <obs> = neg_bin, dispersion = fit <r>orfix_at <value>,location = mean, cumulative, with the matching whole-fit line edition 2 requireslocation = median, nocumulative, a dispersion that differs between the two lines, a second columnnormalization,time_errorrefine,bootstrapmain()would act on them.smoothing,noise_profiling,linear_profiling,output_inference_data,save_best_data,embed_best_fit_data,smooth_plot_points,ind_var_rounding,refine_methodvar/logvar, or aparameter:record with no priorinitialization = lh/rand; other values;initialization_distribution = boundsstart_point,initial_value:,starting_paramsmax_iterationspopulation_sizelwf_independent_runs.lwf_independent_runs,lwf_state_file,lwf_continuerandom_seedwall_time_fit;wall_time_sim;wall_time_genparallel_count-t/-s;-r-rnameslwf_continue).lwf_*values; unknown, repeated or start-onlylwf_parameter_jitternames; anlwf_jitterno moved parameter usesmain()would delete, or an existing file that is not a state filelwfjob