pypi sbmlsim 0.6.1

4 hours ago

Release notes for sbmlsim 0.6.1

sbmlsim

We are pleased to release the next version of sbmlsim. This release is a cleanup of the parameter fitting: sbmlsim.fit was gone through module by module, the bugs the reading surfaced are fixed, reporting a fit is separate from running it, and the fitting has a working reference problem and tests again. The library was pruned to its core at the same time: SED-ML and the COMBINE archives are gone, the generated model code of the simulator comparison is out of the package, and what only served an example lives with that example.

Breaking changes

  • library code does not show figures. sbmlsim.sensitivity.plots.heatmap and S1_ST_barplot return the matplotlib.figure.Figure they create and close it instead of calling plt.show(), the plots of the Morris and the sampling analysis close theirs; every figure is still saved to its path. On a machine with a display the calls opened a window and blocked, under a headless backend they were a warning per figure, and the figures accumulated in the pyplot registry for the length of an analysis. SimulationExperiment.run(show_figures=...) defaults to False, like the ExperimentRunner which calls it; show_figures=True still shows them
  • sbmlsim.utils.deprecated removed. It marked one function, add_data of the removed sbmlsim.plot.plotting_deprecated_matplotlib, and had no other user
  • SED-ML and COMBINE archives removed. sbmlsim.combine is gone: SEDMLReader, SEDMLParser, SEDMLSerializer, execute_sedml, the task tree, the data descriptions and the NuML parser. Nothing in the library used them, they were a partial implementation of L1V4 with a single smoke test, and a fit is exchanged as a PEtab problem instead, see sbmlsim.fit.petab_v2. combine/mathml.py is the one part the library needs, it evaluates the formula of a Data of type FUNCTION, and it is sbmlsim.mathml now. sbmlsim.simulation.kisaos, whose only consumer was the SED-ML parser, is gone as well; sbmlsim.simulation.algorithm still describes an integrator by its KISAO terms. python-libnuml is no longer a dependency, python-libsedml stays, sbmlsim.mathml parses L3 formulas with it. The duplicated combine.datagenerator and combine.sedml.report are sbmlsim.result.datagenerator and sbmlsim.result.report. The examples examples/sedml and examples/covid were SED-ML and COMBINE archive examples and are removed with it
  • sbmlsim.plot.plotting_deprecated_matplotlib removed. Its one function add_data draws a dataset onto a matplotlib axes; no library code called it and its only user was the dose response experiment of the glucose example, where it now lives as examples/glucose/experiments/plotting.py. A simulation experiment describes its figures with Figure, Plot and Curve (sbmlsim.plot.plotting)
  • sbmlsim.interpolation removed. It writes a dataset as an SBML model with interpolation rules, was never integrated with the rest of the package and had one example and one test. It is the example package examples/interpolation
  • the simulator comparison is examples/comparison, sbmlsim.comparison.diff stays. DataSetsComparison compares the results of two simulators and is library code; the scripts which produce such results with AMICI and COPASI (simulate.py, simulate_amici.py, simulate_copasi.py, simulate_roadrunner.py, example_comparison.py) import the optional dependencies and are examples. The 50 MB of generated AMICI model code and compiled extensions under sbmlsim/comparison/, which the wheel, ruff and ty all excluded, are deleted; they are regenerated by the scripts. With them out of the tree the comparison is linted and type checked like the rest of the package
  • reporting is separate from optimizing. A fit produces parameters, a report is made from parameters. OptimizationAnalysis is replaced by FitReport (sbmlsim.fit.report), which is created from the definition of an OptimizationProblem, the FitSettings and one or more ParameterSets, and report.create(output_dir, name) writes it. It does not need an optimization to have been run: ParameterSets.from_json reads the parameters a fit stored, so a report is created again later or for several parameter sets at once, which puts every set in the parameter table and every set as a curve in the plots. FitReport.from_optimization_result is the shortcut for the report of a fit, it reads the settings from the result and adds the plots which describe the runs rather than a parameter set. It reports the fitted parameters alone, with_model=True adds the values the model started from as the reference set. No parameter set is drawn in black any more, which is the color of the reference data of a fit mapping. sbmlsim.fit.analysis is gone
  • FitSettings (sbmlsim.fit.options) bundles residual, loss_function, weighting_curves, weighting_points, variable_step_size and the tolerances. run_optimization(problem, settings=...) and OptimizationProblem.initialize(settings) take one instead of the seven separate arguments, which callers had to hand carry as an untyped dictionary to get the same numbers into the fit and into its report. The settings are stored with the OptimizationResult and serialize with it
  • OptimizationProblem.initialize is a no-op when the problem already has the given settings, so a fit and its report resolve the data once instead of twice. problem.residual, loss_function, weighting_curves and weighting_points are read only properties of the settings
  • OptimizationAnalysis.show_plots is gone with the class; a report writes its figures, show_report=True opens the HTML in a browser. OptimizationAnalysis.plot_correlation was already removed
  • the loss functions of LossFunctionType are applied to the squared residuals, as in scipy.optimize.least_squares. They were applied to the signed residuals, which is not the definition of the loss and returns NaN for every residual r <= -1 (soft_l1 computed 2*(sqrt(1+r)-1), cauchy computed log(1+r)). The residuals are transformed to sign(r)*sqrt(rho(r**2)), so the cost is 0.5*sum(rho(r**2)); the new sbmlsim.fit.optimization.apply_loss_function does this. Fits with SOFT_L1, CAUCHY or ARCTAN give different results than in 0.6.0, fits with the default LINEAR are unchanged
  • the space the optimizer searches is FitSettings.parameter_scale, ParameterScaleType.LOG10 by default, with LOG and LINEAR for the problems which want them. It was hardcoded as log10, which a parameter that is not positive or that does not span orders of magnitude cannot use; the bounds, the start values and the fitted parameters stay on the linear scale, i.e. in the units of the model, and only the search happens in the scaled space. The scale is a property of the optimization and not of the problem, which is why it is part of the settings and why PEtab v2 removed the parameterScale of its parameter table. The scans of the profile likelihood follow the scale as well
  • the optimization runs in logarithmic parameter space, which requires finite positive bounds and start values. OptimizationProblem.initialize checks this and raises a ValueError naming the parameter; before, np.log10 of a non-positive or infinite bound silently produced NaN and the optimizer failed with an opaque error. The default bounds of a FitParameter (-inf, inf) are therefore not usable in a fit
  • OptimizationProblem raises for an empty list of fit parameters, which it only logged as an error, and for duplicate parameter ids. FitExperiment raises when weights are given together with use_mapping_weights=True, which it only logged; FitParameter raises for inverted bounds and a start value outside of the bounds
  • MappingMetaData is a keyword only dataclass. A subclass which added a field without a default, i.e., every application specific metadata, could not be created at all: the outlier field of the base class has a default and came first, so python raised TypeError: non-default argument follows default argument on the class definition. A subclass must not redeclare outlier
  • OptimizationAnalysis does not open a web browser any more, show_report=True does; show_plots defaults to False like the rest of the library, it was True and called plt.show() for every figure
  • OptimizationAnalysis.plot_correlation removed, run had it commented out behind a pass and it was the only user of seaborn in the fitting
  • sbmlsim.fit.sampling is library code again: plot_samples and example_sampling, which called plt.show(), are the new example examples/fit_sampling.py. sbmlsim.fit.rmse lost its __main__ block with the hardcoded numbers of another project, and aic(mse, N, k) is aic(mse, n, k)
  • examples/midazolam removed, examples/hctz_fitting replaces it as the reference problem of the fitting, see Features. The packaged model MIDAZOLAM_SBML of sbmlsim.resources stays, it is used in the documentation and the unit tests

Fixes

  • sbmlsim.fit.optimization: re-initializing an OptimizationProblem appended to the lists of the fit mappings instead of replacing them, so the standard workflow of run_optimization(serial=True) followed by an OptimizationAnalysis counted every curve twice, in the cost, in the report and in the plots
  • sbmlsim.fit.optimization: FitExperiment(mappings=None) is documented as "use all mappings of the experiment", but the constructor turned None into [] and the branch which resolved the mappings tested for None, so an experiment without explicit mappings contributed nothing to the fit. The mappings are resolved by the new FitExperiment.resolve_mappings
  • sbmlsim.fit.optimization: the check for NaN and INF in the reference data tested the x data instead of the errors, so non-finite errors reached the weights; the initial parameter values xmodel were read inside the loop over the fit mappings, where the loop variable k of the mappings was shadowed by the loop over the parameters; residuals(complete_data=True) raised an UnboundLocalError when the integration of the first mapping failed
  • sbmlsim.fit.sampling: an infinite upper bound was replaced by the hard lower bound 1e-10 instead of 1e10, which inverted the sampled interval; logarithmic sampling of a finite non-positive lower bound produced NaN start values. The seeding uses numpy.random.default_rng instead of the legacy global seed, and seed=0 is a seed
  • sbmlsim.fit.result and sbmlsim.fit.analysis: the warning for a parameter which ended up on its bound divided by the bound, which raised ZeroDivisionError for a bound of 0.0 and was meaningless for the infinite default. The shared bound_warnings measures the distance relative to the interval of the parameter, in logarithmic space, i.e., in the space the optimization runs in
  • sbmlsim.fit.result: OptimizationResult.from_json returned the parameter vectors as lists, xopt was not an array after a round trip; combine of an empty list raised an IndexError instead of a clear message
  • sbmlsim.fit.analysis: the HTML report linked svg images although image_format is configurable, and its tags were not closed; _datapoints_df stored the mapping id as a one element tuple; plot_traces tested len(array > 0), the length of a boolean array, instead of the length of the data; the plots directory was created without parents/exist_ok, so a second analysis in the same directory failed; the plots on logarithmic axes handled data points which are zero or negative, which produced divide by zero and "non-positive limits" warnings
  • sbmlsim.fit.petab_omex: every entry of the COMBINE archive was written with master=True, only the PEtab YAML is the master entry; the deprecated top level petab imports are petab.v1; missing files are reported with a FileNotFoundError instead of being silently skipped; the __main__ block pointed at a directory which does not exist
  • sbmlsim.fit.runner: the guards for the removed parameters fitting_type and weighting_local inspected locals(), where they could never appear, instead of kwargs; the seeds of the workers come from a numpy.random.SeedSequence and are distinct; the module level multiprocessing.Lock, which does nothing under the forkserver start method of python 3.14, is gone
  • the links a fit and a report print are file:// URIs, i.e., file:///C:/... with forward slashes on windows, where the path with the backslashes it printed is not a link a terminal or webbrowser.open opens
  • sbmlsim.fit.objects: FitMapping.weight raises a clear ValueError when neither a weight nor a count is available, it relied on an AttributeError of an attribute which was only set conditionally

Features

  • examples/comparison/diff_example.py is the example of sbmlsim.comparison.diff: it simulates the six timecourses of examples/comparison/diff/ with roadrunner, compares them against the results of the same simulations on JWS Online and writes the reports and figures into the working directory. It replaces tests/data/diff/simulate_examples.py, whose paths pointed at directories which do not exist, which is why both of its tests were skipped
  • parameter identifiability by profile likelihood. sbmlsim.fit.identifiability.profile_likelihood computes the profile of every fitted parameter around a parameter set: the parameter is fixed at values around the optimum, the other parameters are optimized again, and the profile is compared with the threshold of the likelihood ratio test, cost_min + chi2.ppf(alpha, df) / 2 in the cost convention of the fit (0.5 * Σ r², i.e., 1.92 above the minimal cost at 95% and one degree of freedom). The values at which the profile crosses the threshold are the confidence interval, and the shape of the profile classifies the parameter (Raue et al. 2009): identifiable, practically non-identifiable towards small and/or large values, or structurally non-identifiable when the profile is flat. The scans run in logarithmic space with adaptive steps, two per parameter in the worker pool of the fit runner, the other parameters start from the previous point and their paths along the profile are stored, so coupled parameters are seen. ProfileSettings holds the confidence level, the steps and reoptimize=False for a plain scan of the cost, IdentifiabilityResult the profiles with their intervals and classification, as table, text report and JSON. FitReport(identifiability=...) adds the identifiability section to the report, with the table of the parameters, the overview of the profiles and one figure per parameter with the paths of the other parameters; identifiability_cli is the command line tool for stored parameters and FitRun.identifiability the analysis of a finished fit. examples/hctz_fitting/fitting/identifiability.py runs a global optimization followed by the profiles of its best parameter set. The method and its references are documented in Parameter fitting
  • an initialized OptimizationProblem pickles as its definition: the experiment runner with the models and the unit registry, which cannot be pickled, is left out, so a problem which was initialized in the main process is handed to the workers of a pool without creating it again
  • the outliers of a fit are evaluated. A curve which a fit drops as unusable is resolved and simulated like the validation data, it stays out of the cost, and summary_df() has a row for it next to the training and the validation data, so a report says where the dropped data sits relative to the model and the decision to drop it is checkable: an outlier with an R² as good as the training data was dropped without reason. The row over all data points is gone, it pooled the data a fit was fitted on with the data it dropped; summary() without a kind still gives that number. A report shows what a fit evaluates and nothing else, i.e. the excluded data has no column in the overview of the data, no filter chip and no figures. EVALUATED_KINDS is training, validation and outlier now and UNUSED_KINDS is the excluded data alone, which is the only kind a problem does not resolve; outlier_indices is on the problem next to training_indices and validation_indices, and a PEtab problem carries the outliers with their kind, so the round trip keeps them
  • the report shows the goodness of fit and a Bland-Altman plot per subset of the data. goodness_of_fit is the prediction against the measurement on logarithmic axes with the identity line, bland_altman the agreement of the two as a ratio, log10(f(x)/y) over the geometric mean, with the bias and the limits of agreement bias ± 1.96 SD written as fold factors, e.g. an unbiased fit within a factor of two is bias 1.02x, LoA 0.51-2.04x. The limits are those of the training data and are the same in every panel, so the validation data and the outliers are read against what the fit agrees to. Both have one panel per kind of fit mapping, so the training, the validation and the outlier data are read next to each other, and their points carry the color of their study, i.e. of the simulation experiment a fit mapping belongs to, so a study is the same color in both figures and a study which the model describes badly is seen. Several parameter sets in one report are told apart by their marker. Both draw the same band in the same styles, i.e. a solid line for prediction = measurement, a dash-dotted line for the bias, dashed lines for the limits and the area between them filled, which on logarithmic axes is a horizontal band in the one figure and a band parallel to the diagonal in the other. goodness_of_fit replaces datapoint_scatter, which was the same plot over all data at once and carried a band of a fixed factor of ten instead, and residual_scatter, the relative residuals over the data, is gone: the Bland-Altman plot is the same information with a reference to read it against
  • the data of a fit is training data, validation data, an outlier or excluded. OUTLIER and EXCLUDED are both unused, and they say different things: an outlier is a decision about the data, i.e. the curve is not usable, and an exclusion is a decision about the model, i.e. the model does not describe what was measured, e.g. an arm of a study with a coadministration the model has no interaction for. EVALUATED_KINDS and UNUSED_KINDS are the two halves of MappingKind
  • the data of a fit is selected from the complete list of fit mappings in three steps. sbmlsim.fit.helpers.FitMappings instantiates the simulation experiments once and select(filters, outliers, validation) sets the kind of every fit mapping: the filters select the training data and a mapping which fails one is EXCLUDED, the outlier keys tag the training data which is not usable as OUTLIER, and the validation data is the part of what is left which is named by its keys or selected by a filter, so excluded beats outlier, outlier beats validation and validation beats training. The outliers are a decision about the data and are named once for all fits, the filters and the validation data are the definition of a fit problem. The result is a MappingSelection with the kind of every mapping, the overview table and one FitMappingCollection per experiment and kind, select_mapping_collections is the one call version. It replaces filtered_mapping_collections, f_collection and mapping_collections_by_kind, which took separate filters per kind and left it to the caller to keep them disjoint, and filter_keys, filter_not_keys and filter_empty, which are the keys, no filters and the complement the selection does itself
  • the data of a fit is training data, validation data or an outlier. The MappingKind of a FitExperiment decides what a fit does with the mappings it selects: TRAINING (the default) enters the cost of the optimization, VALIDATION is simulated and evaluated with the fit but not fitted, and OUTLIER is not used at all, an optimization problem skips it. The kind belongs to the selection of the data and not to the fit mappings of a simulation experiment, whose MappingMetaData describes the curve: the same curve is training data of one fit and validation data of another. The number of mappings per kind is part of the overview of the data and of the report of the problem. The metrics are calculated for the training data, for the validation data and over all data points, so a fit is evaluated on the data it was not fitted on. MappingMetaData.outlier is gone with the flag
  • the progress bar of a fit shows an estimate of the total runtime next to the elapsed time, e.g. ~ 0:12:30 total, from the runs which are done and the number of workers; the profile likelihood scans show the same. The PKIV fit of the HCTZ example is gone, PK is the reference problem; the tests build the iv subset of a single study themselves where a problem whose data does not determine the absorption parameters is needed
  • the metrics of a fit are scale free. The data of a fit spans orders of magnitude, so the absolute RMSE only reads the largest curves and a curve of small values which is missed by a factor of 20 still has a tiny error. NRES, the residual over the mean of its curve, and NRMSE, its root mean square, count every curve the same, and IWRES is the residual of the cost now, i.e. the residual of the residual type, weighted and with the loss function applied, so the cost of the training data is 0.5 * n * RMSE_w²; before it was the absolute residual times the weight. Both are columns of datapoints_df, mappings_df and summary, and the goodness of fit and the Bland-Altman plot cover a full row of the report with a box of key metrics per panel: , NRMSE and RMSE_w in the one, the bias and the SD of the panel as fold factors and the share of its points inside the limits of agreement in the other. The table of the fit mappings on the console carries the complete metadata of every mapping and greys out the excluded rows
  • the HTML report of the simulation experiments is interactive as well and shares its design with the report of a fit: report_base.html carries the style and the script and the reports fill in its blocks. The page of an experiment has the sections overview (its models and datasets), figures and code, the index lists the experiments with their figures, and both are searchable. They load nothing from the network any more, the templates pulled Vue, Vuetify, fonts and highlight.js from CDNs, so a report only rendered with a network connection; the source of an experiment is escaped now, code with a < in it broke the page
  • the HTML report of a fit is an interactive page with three sections: the overview of the fit, i.e., what the console reports, the results with the metrics and the plots, and one card per fit mapping. A search box filters the mappings and the tables, chips filter the training, validation and outlier data, the tables sort by any column and a figure opens full size when it is clicked. It is rendered from the jinja2 template resources/templates/fit_report.html, carries its own style and script and loads nothing from the network, so a report is archived or sent as it is
  • sbmlsim.fit.cli holds what every fit has in common, so that a model only defines its fits. A FitDefinition is the fit experiments which enter a fit, the parameters which are adjusted, the paths of the experiments and the settings; run_fit creates the optimization problems for an OptimizationStrategy (ALL fits the experiments together, SINGLE fits every experiment on its own) and runs them, and FitRun carries the problem with its result and creates the report. fit_cli and report_cli are the command line tools around this: a model passes its definitions by name and gets the tool. examples/hctz_fitting/fitting/fitting.py went from 301 lines to the definitions plus four lines of main
  • the report shows the Fisher information: FitReport(..., fisher=fim) adds its table of errors and intervals, its eigenvalues and the correlation of the parameters to the identifiability section and writes fisher.json and fisher.tsv, and a report of an information which is rank deficient says that its errors cannot be read. FitRun.fisher() computes it for a finished fit, like FitRun.identifiability()
  • every metric and every figure of a report explains itself: a ? next to a column of the metrics says what the value is, e.g. what the BIC charges a parameter and why R² of a non-linear model is not the square of a correlation, and one next to a figure says what to look for in it, e.g. that a flat plateau of the waterfall plot is the evidence that the multistart converged. The tooltips are css in the shared report_base.html, so the reports still load nothing from the network
  • sbmlsim.fit.fisher is the Fisher information of a fit, i.e. the local analysis next to the profile likelihood: FIM = J'J from one jacobian of the weighted residuals, and from it the standard errors, the correlations of the parameters, the confidence intervals θ ± t·SE in the space the optimizer searches, the eigenvalues, the rank and the condition number. A direction with a small eigenvalue is a combination of parameters the data does not determine, so the analysis finds a structural non-identifiability without scanning: the intravenous problem of the HCTZ example determines the renal excretion and carries no information about the absorption of an oral dose, and its information has rank 1 of 3. The Fisher information is the cheap local answer and the profile likelihood the exact one; they disagree where the cost is not a quadratic, which is the normal case for a non-linear model. FitRun.fisher() computes it for a finished fit and FitReport(fisher=...) writes it into the identifiability section next to the profiles, which identifiability_cli does unless --no-fisher is given
  • the metrics of a fit report the BIC next to the AIC. The two differ in what a parameter costs, k * ln(n) against 2 * k, so the BIC prefers the smaller of two models which describe the data equally well as soon as there are more than seven data points. sbmlsim.fit.metrics.bic and bic_from_mse are the functions
  • Parameter fitting explains what identifiability is before it explains the methods: structural non-identifiability as a property of the model and its observables, practical non-identifiability as a property of the model and the data at hand, sloppiness, what the profile likelihood and the Fisher information each cost and each find, and what to do with a parameter which comes back undetermined. The references of the two are in References
  • sbmlsim.fit.metrics (the former sbmlsim.fit.rmse) calculates the metrics of a parameter set on an optimization problem. FitMetrics.datapoints_df is the table of the data with the predictions in the convention of population pharmacokinetics (DV, PRED, IPRED, RES, IRES, IWRES), mappings_df the metrics per fit mapping and summary the metrics over all data points: the number of data points, the number of parameters, the cost, MSE, RMSE, RMSE_w, R2 and AIC. A report calculates them for every one of its parameter sets and writes metrics.tsv, metrics_mappings.tsv and datapoints.tsv, so parameter sets are compared by their AIC, RMSE and R². The functions sse, mse, rmse, aic and r_squared work on plain arrays; rmse and aic take the residuals now, rmse_from_mse and aic_from_mse take the mean squared error as rmse and aic did in rmse.py
  • sbmlsim.fit.parameters: ParameterSet is a named set of parameter values with their units, their cost and where they come from, ParameterSets holds one or more of them and stores them as JSON. OptimizationResult.parameter_set(k) and parameter_sets(size) hand out the best runs of a fit, OptimizationProblem.parameter_set_model() the values the models start from, which is the reference a fit is compared against
  • a fit keeps what it has. A single optimization which fails is a result with its message and the other repeats are unaffected, timeout gives every repeat a budget and a repeat which runs out of it keeps the best parameters it reached, and a worker of a parallel fit which dies loses only its own repeats. RuntimeErrorOptimizeResult is a scipy.optimize.OptimizeResult now, so a failed repeat is stored, serialized and reported like a successful one; it was a dataclass without the dictionary interface the serialization uses, which raised an AttributeError as soon as a repeat failed
  • FitExperiment is FitMappingCollection. A fit experiment was never an experiment: it selects the fit mappings of one SimulationExperiment, says what a fit does with them and how they are weighted, i.e. it is a collection of fit mappings, and the name collided with the simulation experiment it selects from and with the experiment of a PEtab problem. A collection carries an sid, Beermann1976_training by default, which names what it produces. OptimizationProblem(mapping_collections=...), FitDefinition.mapping_collections and FitDefinition.collections(study_ids), helpers.mapping_collections_by_kind and merge_mapping_collections follow the name, and examples/hctz_fitting/fitting/fit_experiments.py is mapping_collections.py
  • an observable of a PEtab problem which is a formula over the entities of a model, which is what the benchmark collection mostly has, is an entity of the model the fit simulates: sbmlsim.fit.petab_v2.observables writes a copy of the model in which every such observable is a parameter with an assignment rule, and the fit selects it. The math of PEtab is the math of the model, so the formula is the rule; the observables of Boehm_JProteomeRes2014 agree with the simulatedData of the collection to 2e-4. The parameters of a problem which are not estimated are applied to the model with the nominal values of the parameter table, which PEtab prescribes and which the reader ignored
  • examples/petab/benchmark.py fits problems sbmlsim did not write: Perelson_Science1996 and Boehm_JProteomeRes2014 of the PEtab benchmark collection, converted from PEtab 1.0 with petab.v2.petab1to2, read, fitted, analysed with the profile likelihood and reported. Reading a problem of another tool needed three things the HCTZ example does not have: a measurement which names no experiment is "use the model as is" and gets a simulation of the model over its measurements, the units of the data are the units of the model when there is no sbmlsim extension to say otherwise, and a parameter of the problem which is not an entity of a model, i.e. the standard deviation of an observable, is not fitted because the objective of sbmlsim weights the data instead of estimating its noise
  • sbmlsim.fit.petab_v2 is the PEtab v2 layer of the fitting. to_petab(problem, directory) writes an optimization problem as a PEtab v2 problem, built from the fit mapping collections: a collection whose mappings share a simulation is exactly one experiment of PEtab and carries its id, a collection over several simulations, e.g. the doses of a study, is one experiment per simulation numbered after it, and the reader gives one collection back per experiment, and from_petab(yaml_file) reads a PEtab v2 problem as an OptimizationProblem and the settings to run it with; PetabReader builds a SimulationExperiment from the tables the way the SED-ML parser builds one from a document. The problem which is written validates with petab itself
  • what PEtab does not express goes into the sbmlsim extension of the problem, i.e. the units of the parameters, the observables and the data, the FitSettings, whether a mapping is training, validation or outlier data, the weight it was given, the output grid of the timecourses and the settings of the integrator. The extension is required, which PEtab asks of an extension which changes the mathematical interpretation of a problem: the settings are the objective of the fit, so a tool which does not know sbmlsim rejects the problem instead of fitting the same data with another objective without saying so, and to_petab(..., required_extension=False) writes a problem for those tools. A round trip through sbmlsim keeps the fit: the HCTZ example comes back with its settings, its parameters and units, its mappings with their kinds and the reference data of every mapping, and its cost agrees to 7e-6
  • the layer follows the PEtab v2 specification where it differs from v1: the standard deviation of a measurement is the sd placeholder an observable declares in noisePlaceholders and not a noiseParameter${n}_${observableId} name, the parameter table has no parameterScale any more, and sbmlsim.fit.petab_v2.symbols converts the selections of roadrunner into the math of a model instead of dropping the brackets: S1 is the amount of a species and [S1] its concentration, while the identifier in the math of SBML is the amount for hasOnlySubstanceUnits=true and the concentration for false, so the concentration of an amount based species is written as S1 / compartment. A change which sets the concentration of an amount based species is not an identifier a condition can assign and the export raises
  • sbmlsim.fit.petab_v2.gaps is the catalogue of what PEtab v2 cannot express about an sbmlsim fit, and gaps_of_problem(problem) reports the gaps a problem runs into before it is written. A gap is carried by the extension, is lossy, i.e. the weights which are not a standard deviation and the pre-simulations which are not a pre-equilibration, or is unsupported, i.e. a structural model change, an observable which is a python function and a mapping over something else than time, for which the export raises rather than write a problem which means something else
  • a parallel fit hands every repeat to the worker which is free. The repeats were split over the workers before they ran, so a worker which drew the repeats which converge slowly kept the others waiting; every repeat is a task of the pool now, which also means a worker which dies loses the one repeat it was running instead of its whole share. Every worker initializes the problem once, as the initializer of the pool, and the runner stores the repeats, so the files of runs_dir no longer carry the index of the worker which ran them. The progress of the runs comes from the tasks, so the fit does not start a manager process for a queue any more
  • the start values of a fit do not depend on the number of workers. Every worker sampled its own start points from its own seed, so the same seed gave a different fit for a different number of cores, and latin hypercube sampling covered the repeats of a worker instead of the repeats of the fit. The runner samples all start points and hands one to every repeat: a fit with a seed gives the same start points for one worker or twelve, and a serial fit gives the same as a parallel one
  • every run of DIFFERENTIAL_EVOLUTION returned the same result when the fit had a seed: the seed was handed to the optimizer once and every repeat drew the same population from it, so size repeats were one repeat, size times. Every repeat gets its own seed, derived from the seed of the fit with a numpy.random.SeedSequence
  • a fit which is missing the if __name__ == "__main__": guard says so. The workers import the script, the script fits again and its workers fit again: python reports its own error in every worker, the pool replaces the workers which died, and the fit fills the terminal with tracebacks and does not end. The runner sees that the workers it started are gone and stops with the guard and how to write it; a worker which starts a parallel fit stops right there
  • the workers of a parallel fit are capped at the number of repeats, run_optimization raised the number of repeats to the number of cores instead, i.e., size=3 with 12 cores ran 12 optimizations
  • under the forkserver start method, the default on linux since python 3.14, the forkserver of a fit preloads sbmlsim and the modules of the experiments, so only the first fit of a process pays for importing them: three fits of the HCTZ PK problem in one process went from 7.4 s to 6.9 s, and the pools after the first start in 0.02 s instead of 1.4 s
  • an optimization stores what its report needs and nothing else. A repeat kept every field scipy.optimize.least_squares returns, i.e., the residuals fun, the jacobian jac, the gradient and the active mask, which grow with the data and are not read anywhere; it keeps x, x0, cost, success, status, message, duration, optimality and nfev now. A trajectory is the cost of every step, the parameter vector of every step was stored for the correlation plot which is gone. A repeat of the HCTZ PK fit went from 36 kB to 4 kB and its optimization_result.json from 72 kB to 8 kB, and a repeat which is interrupted keeps the best parameters it reached, which the trajectory carried before
  • the repeats of a fit are serialized on their own. run_optimization(runs_dir=...) writes every repeat as JSON the moment it finishes, so a fit which is interrupted or crashes leaves the repeats which are done, and OptimizationResult.from_directory reads them back into a result; OptimizationResult.run_result and write_run are the single repeat
  • the fit mappings which share a simulation are simulated together. Several mappings read different observables of one simulation, and the simulation ran once per mapping; it runs once per group now with the selections of all of them. The 30 fit mappings of the HCTZ example are 9 simulations, and an evaluation of the residuals went from 76 ms to 22 ms, i.e. the fit is about 3.4 times faster. The parameter quantities are created once per evaluation instead of once per mapping. Profiling says the remaining time is the integrator of libroadrunner (63% of a fit) and the unit handling of pint, not python arithmetic, so there is nothing for a jit to speed up
  • the console output of a fit is a sequence of sections, each with its own icon, which say what is being fitted: the fit with its strategy, algorithm and paths, a table of the parameters which are optimized with their bounds and units, a table of the settings, the data with the number of fit mappings per experiment and kind followed by the mappings themselves, the optimization and the report. sbmlsim.fit.display renders them, so the runner, the command line tools and an interactive session look the same
  • a fit gets its id when it starts, <problem>_<date>_<time>__<hash> from fit_id, and it is the id of the optimization problem, of its OptimizationResult, of its parameter sets and of the directory of its report, so everything a fit produces carries the same key and sorts by time
  • the runner shows the progress of the optimizations on a rich progress bar with a spinner, the elapsed time and the finished runs. The workers report every finished run through a queue, so the progress moves while they work; show_progress=False turns it off. The output of a fit is one block: the problem, the runs and the workers, the progress, how many runs converged with the best cost, and the link to the report
  • the messages of a fit are reported once. The mappings without errors in the reference data were reported one message per mapping and once per initialization, they are one message with the count now, and the worker processes, which resolve the same data on every core, only report their errors
  • examples/hctz_fitting/fitting/run_report.py creates the report of a finished fit from its parameters.json, and compares the parameters of several fits in one report
  • examples/hctz_fitting: a whole body pharmacokinetics model of hydrochlorothiazide with the simulation experiments of Beermann 1976 and Patel 1984, their data, the metadata of the fit mappings and the fit problems built on them. It is the reference problem of sbmlsim.fit: python -m examples.hctz_fitting.simulations runs the experiments, python -m examples.hctz_fitting.fitting.fitting the fit, and op_hctz creates the problem of a subset of the data
  • examples/fit_sampling.py: the four sampling types of the start values of a fit, compared on a figure
  • OptimizationAnalysis evaluates the residuals once per parameter vector and caches them. The analysis simulated every fit mapping eight times, once for each plot and table
  • sbmlsim.fit.rmse.mse completes rmse and aic; SamplingType.is_log and is_lhs; sbmlsim.fit.runner.resolve_n_cores, which uses os.process_cpu_count, i.e., the cores the process may use; MappingMetaData is exported from sbmlsim.fit
  • every module of sbmlsim.fit carries full type annotations and google style docstrings which name what is raised

Documentation

  • Parameter fitting rewritten against the HCTZ example, with the requirements on the parameters, the definition of the loss, the settings, the parameter sets and the separate report. Every code block was run
  • examples/README.md and CLAUDE.md describe the HCTZ example and how the fit tests use it
  • citation information in the README

Development

  • the pruning removes 606 files and about 166 000 lines from the repository, 50 MB of generated AMICI code and compiled extensions in the package tree and 18 MB of test data no test referenced (tests/data/combine, tests/data/sedml, tests/data/sedml_fit, tests/data/data/omex and the checked in outputs of tests/data/diff). The skipped tests tests/comparison/test_diff.py and tests/experiment/test_covid_examples.py are gone with the code they waited for
  • the branch protection is part of the repository as rulesets in .github/rulesets/, together with CODEOWNERS and a pull request template; main tracks the latest release and is fast-forwarded by the sync-main job of the release workflow, and the release is prepared on a branch and tagged on develop after the merge, so bump-my-version no longer tags
  • the tests of the fitting run again: tests/fit/ has 54 tests and no skips, they were all skipped with "no fit support" because they needed the midazolam example. tests/fit/conftest.py builds the problem from examples.hctz_fitting, tests/fit/test_petab_omex.py creates and reads an archive of the Boehm PEtab problem of examples/petab, and the analysis test asserts the files of the report
  • examples.hctz_fitting.simulations and examples.fit_sampling are part of the example tests
  • the tests run in parallel with pytest-xdist, pytest takes 40 seconds on a workstation instead of 3.5 minutes; pytest -n 0 runs them in one process
  • helpers.fit_experiments_by_kind initializes the simulation experiments once and filters them for every kind, instead of loading the models and the datasets once per kind, which makes the creation of a fit problem three times faster
  • examples/hctz_fitting/fitting/fitting.py and report.py run as scripts as well, i.e., through the "run file" of an IDE, which does not put the repository on sys.path; they failed with ModuleNotFoundError: No module named 'examples'

Your sbmlsim team

Don't miss a new sbmlsim release

NewReleases is sending notifications on new releases.