Release notes for sbmlsim 0.6.1
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.heatmapandS1_ST_barplotreturn thematplotlib.figure.Figurethey create and close it instead of callingplt.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 toFalse, like theExperimentRunnerwhich calls it;show_figures=Truestill shows them sbmlsim.utils.deprecatedremoved. It marked one function,add_dataof the removedsbmlsim.plot.plotting_deprecated_matplotlib, and had no other user- SED-ML and COMBINE archives removed.
sbmlsim.combineis 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, seesbmlsim.fit.petab_v2.combine/mathml.pyis the one part the library needs, it evaluates the formula of aDataof typeFUNCTION, and it issbmlsim.mathmlnow.sbmlsim.simulation.kisaos, whose only consumer was the SED-ML parser, is gone as well;sbmlsim.simulation.algorithmstill describes an integrator by its KISAO terms.python-libnumlis no longer a dependency,python-libsedmlstays,sbmlsim.mathmlparses L3 formulas with it. The duplicatedcombine.datageneratorandcombine.sedml.reportaresbmlsim.result.datageneratorandsbmlsim.result.report. The examplesexamples/sedmlandexamples/covidwere SED-ML and COMBINE archive examples and are removed with it sbmlsim.plot.plotting_deprecated_matplotlibremoved. Its one functionadd_datadraws 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 asexamples/glucose/experiments/plotting.py. A simulation experiment describes its figures withFigure,PlotandCurve(sbmlsim.plot.plotting)sbmlsim.interpolationremoved. 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 packageexamples/interpolation- the simulator comparison is
examples/comparison,sbmlsim.comparison.diffstays.DataSetsComparisoncompares 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 undersbmlsim/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.
OptimizationAnalysisis replaced byFitReport(sbmlsim.fit.report), which is created from the definition of anOptimizationProblem, theFitSettingsand one or moreParameterSets, andreport.create(output_dir, name)writes it. It does not need an optimization to have been run:ParameterSets.from_jsonreads 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_resultis 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=Trueadds 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.analysisis gone FitSettings(sbmlsim.fit.options) bundlesresidual,loss_function,weighting_curves,weighting_points,variable_step_sizeand the tolerances.run_optimization(problem, settings=...)andOptimizationProblem.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 theOptimizationResultand serialize with itOptimizationProblem.initializeis 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_curvesandweighting_pointsare read only properties of the settingsOptimizationAnalysis.show_plotsis gone with the class; a report writes its figures,show_report=Trueopens the HTML in a browser.OptimizationAnalysis.plot_correlationwas already removed- the loss functions of
LossFunctionTypeare applied to the squared residuals, as inscipy.optimize.least_squares. They were applied to the signed residuals, which is not the definition of the loss and returns NaN for every residualr <= -1(soft_l1computed2*(sqrt(1+r)-1),cauchycomputedlog(1+r)). The residuals are transformed tosign(r)*sqrt(rho(r**2)), so the cost is0.5*sum(rho(r**2)); the newsbmlsim.fit.optimization.apply_loss_functiondoes this. Fits withSOFT_L1,CAUCHYorARCTANgive different results than in 0.6.0, fits with the defaultLINEARare unchanged - the space the optimizer searches is
FitSettings.parameter_scale,ParameterScaleType.LOG10by default, withLOGandLINEARfor 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 theparameterScaleof 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.initializechecks this and raises aValueErrornaming the parameter; before,np.log10of a non-positive or infinite bound silently produced NaN and the optimizer failed with an opaque error. The default bounds of aFitParameter(-inf,inf) are therefore not usable in a fit OptimizationProblemraises for an empty list of fit parameters, which it only logged as an error, and for duplicate parameter ids.FitExperimentraises whenweightsare given together withuse_mapping_weights=True, which it only logged;FitParameterraises for inverted bounds and a start value outside of the boundsMappingMetaDatais 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: theoutlierfield of the base class has a default and came first, so python raisedTypeError: non-default argument follows default argumenton the class definition. A subclass must not redeclareoutlierOptimizationAnalysisdoes not open a web browser any more,show_report=Truedoes;show_plotsdefaults toFalselike the rest of the library, it wasTrueand calledplt.show()for every figureOptimizationAnalysis.plot_correlationremoved,runhad it commented out behind apassand it was the only user of seaborn in the fittingsbmlsim.fit.samplingis library code again:plot_samplesandexample_sampling, which calledplt.show(), are the new exampleexamples/fit_sampling.py.sbmlsim.fit.rmselost its__main__block with the hardcoded numbers of another project, andaic(mse, N, k)isaic(mse, n, k)examples/midazolamremoved,examples/hctz_fittingreplaces it as the reference problem of the fitting, see Features. The packaged modelMIDAZOLAM_SBMLofsbmlsim.resourcesstays, it is used in the documentation and the unit tests
Fixes
sbmlsim.fit.optimization: re-initializing anOptimizationProblemappended to the lists of the fit mappings instead of replacing them, so the standard workflow ofrun_optimization(serial=True)followed by anOptimizationAnalysiscounted every curve twice, in the cost, in the report and in the plotssbmlsim.fit.optimization:FitExperiment(mappings=None)is documented as "use all mappings of the experiment", but the constructor turnedNoneinto[]and the branch which resolved the mappings tested forNone, so an experiment without explicit mappings contributed nothing to the fit. The mappings are resolved by the newFitExperiment.resolve_mappingssbmlsim.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 valuesxmodelwere read inside the loop over the fit mappings, where the loop variablekof the mappings was shadowed by the loop over the parameters;residuals(complete_data=True)raised anUnboundLocalErrorwhen the integration of the first mapping failedsbmlsim.fit.sampling: an infinite upper bound was replaced by the hard lower bound1e-10instead of1e10, which inverted the sampled interval; logarithmic sampling of a finite non-positive lower bound produced NaN start values. The seeding usesnumpy.random.default_rnginstead of the legacy global seed, andseed=0is a seedsbmlsim.fit.resultandsbmlsim.fit.analysis: the warning for a parameter which ended up on its bound divided by the bound, which raisedZeroDivisionErrorfor a bound of0.0and was meaningless for the infinite default. The sharedbound_warningsmeasures the distance relative to the interval of the parameter, in logarithmic space, i.e., in the space the optimization runs insbmlsim.fit.result:OptimizationResult.from_jsonreturned the parameter vectors as lists,xoptwas not an array after a round trip;combineof an empty list raised anIndexErrorinstead of a clear messagesbmlsim.fit.analysis: the HTML report linkedsvgimages althoughimage_formatis configurable, and its tags were not closed;_datapoints_dfstored the mapping id as a one element tuple;plot_tracestestedlen(array > 0), the length of a boolean array, instead of the length of the data; the plots directory was created withoutparents/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 produceddivide by zeroand "non-positive limits" warningssbmlsim.fit.petab_omex: every entry of the COMBINE archive was written withmaster=True, only the PEtab YAML is the master entry; the deprecated top levelpetabimports arepetab.v1; missing files are reported with aFileNotFoundErrorinstead of being silently skipped; the__main__block pointed at a directory which does not existsbmlsim.fit.runner: the guards for the removed parametersfitting_typeandweighting_localinspectedlocals(), where they could never appear, instead ofkwargs; the seeds of the workers come from anumpy.random.SeedSequenceand are distinct; the module levelmultiprocessing.Lock, which does nothing under theforkserverstart 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 orwebbrowser.openopens sbmlsim.fit.objects:FitMapping.weightraises a clearValueErrorwhen neither a weight nor a count is available, it relied on anAttributeErrorof an attribute which was only set conditionally
Features
examples/comparison/diff_example.pyis the example ofsbmlsim.comparison.diff: it simulates the six timecourses ofexamples/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 replacestests/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_likelihoodcomputes 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) / 2in the cost convention of the fit (0.5 * Σ r², i.e.,1.92above 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.ProfileSettingsholds the confidence level, the steps andreoptimize=Falsefor a plain scan of the cost,IdentifiabilityResultthe 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_cliis the command line tool for stored parameters andFitRun.identifiabilitythe analysis of a finished fit.examples/hctz_fitting/fitting/identifiability.pyruns a global optimization followed by the profiles of its best parameter set. The method and its references are documented in Parameter fitting - an initialized
OptimizationProblempickles 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_KINDSis training, validation and outlier now andUNUSED_KINDSis the excluded data alone, which is the only kind a problem does not resolve;outlier_indicesis on the problem next totraining_indicesandvalidation_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_fitis the prediction against the measurement on logarithmic axes with the identity line,bland_altmanthe agreement of the two as a ratio,log10(f(x)/y)over the geometric mean, with the bias and the limits of agreementbias ± 1.96 SDwritten as fold factors, e.g. an unbiased fit within a factor of two isbias 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 forprediction = 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_fitreplacesdatapoint_scatter, which was the same plot over all data at once and carried a band of a fixed factor of ten instead, andresidual_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.
OUTLIERandEXCLUDEDare 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_KINDSandUNUSED_KINDSare the two halves ofMappingKind - the data of a fit is selected from the complete list of fit mappings in three steps.
sbmlsim.fit.helpers.FitMappingsinstantiates the simulation experiments once andselect(filters, outliers, validation)sets the kind of every fit mapping: the filters select the training data and a mapping which fails one isEXCLUDED, the outlier keys tag the training data which is not usable asOUTLIER, 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 aMappingSelectionwith the kind of every mapping, the overview table and oneFitMappingCollectionper experiment and kind,select_mapping_collectionsis the one call version. It replacesfiltered_mapping_collections,f_collectionandmapping_collections_by_kind, which took separate filters per kind and left it to the caller to keep them disjoint, andfilter_keys,filter_not_keysandfilter_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
MappingKindof aFitExperimentdecides what a fit does with the mappings it selects:TRAINING(the default) enters the cost of the optimization,VALIDATIONis simulated and evaluated with the fit but not fitted, andOUTLIERis 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, whoseMappingMetaDatadescribes 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.outlieris 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. ThePKIVfit of the HCTZ example is gone,PKis 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, andNRMSE, its root mean square, count every curve the same, andIWRESis 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 is0.5 * n * RMSE_w²; before it was the absolute residual times the weight. Both are columns ofdatapoints_df,mappings_dfandsummary, 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:R²,NRMSEandRMSE_win 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.htmlcarries 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.cliholds what every fit has in common, so that a model only defines its fits. AFitDefinitionis the fit experiments which enter a fit, the parameters which are adjusted, the paths of the experiments and the settings;run_fitcreates the optimization problems for anOptimizationStrategy(ALLfits the experiments together,SINGLEfits every experiment on its own) and runs them, andFitRuncarries the problem with its result and creates the report.fit_cliandreport_cliare the command line tools around this: a model passes its definitions by name and gets the tool.examples/hctz_fitting/fitting/fitting.pywent from 301 lines to the definitions plus four lines ofmain- 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 writesfisher.jsonandfisher.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, likeFitRun.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 sharedreport_base.html, so the reports still load nothing from the network sbmlsim.fit.fisheris the Fisher information of a fit, i.e. the local analysis next to the profile likelihood:FIM = J'Jfrom one jacobian of the weighted residuals, and from it the standard errors, the correlations of the parameters, the confidence intervalsθ ± t·SEin 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 andFitReport(fisher=...)writes it into the identifiability section next to the profiles, whichidentifiability_clidoes unless--no-fisheris given- the metrics of a fit report the
BICnext to theAIC. The two differ in what a parameter costs,k * ln(n)against2 * 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.bicandbic_from_mseare 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 formersbmlsim.fit.rmse) calculates the metrics of a parameter set on an optimization problem.FitMetrics.datapoints_dfis the table of the data with the predictions in the convention of population pharmacokinetics (DV,PRED,IPRED,RES,IRES,IWRES),mappings_dfthe metrics per fit mapping andsummarythe metrics over all data points: the number of data points, the number of parameters, the cost,MSE,RMSE,RMSE_w,R2andAIC. A report calculates them for every one of its parameter sets and writesmetrics.tsv,metrics_mappings.tsvanddatapoints.tsv, so parameter sets are compared by their AIC, RMSE and R². The functionssse,mse,rmse,aicandr_squaredwork on plain arrays;rmseandaictake the residuals now,rmse_from_mseandaic_from_msetake the mean squared error asrmseandaicdid inrmse.pysbmlsim.fit.parameters:ParameterSetis a named set of parameter values with their units, their cost and where they come from,ParameterSetsholds one or more of them and stores them as JSON.OptimizationResult.parameter_set(k)andparameter_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,
timeoutgives 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.RuntimeErrorOptimizeResultis ascipy.optimize.OptimizeResultnow, 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 anAttributeErroras soon as a repeat failed FitExperimentisFitMappingCollection. A fit experiment was never an experiment: it selects the fit mappings of oneSimulationExperiment, 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 ansid,Beermann1976_trainingby default, which names what it produces.OptimizationProblem(mapping_collections=...),FitDefinition.mapping_collectionsandFitDefinition.collections(study_ids),helpers.mapping_collections_by_kindandmerge_mapping_collectionsfollow the name, andexamples/hctz_fitting/fitting/fit_experiments.pyismapping_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.observableswrites 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 ofBoehm_JProteomeRes2014agree with thesimulatedDataof the collection to2e-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.pyfits problemssbmlsimdid not write:Perelson_Science1996andBoehm_JProteomeRes2014of the PEtab benchmark collection, converted from PEtab 1.0 withpetab.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 nosbmlsimextension 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 ofsbmlsimweights the data instead of estimating its noisesbmlsim.fit.petab_v2is 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, andfrom_petab(yaml_file)reads a PEtab v2 problem as anOptimizationProblemand the settings to run it with;PetabReaderbuilds aSimulationExperimentfrom the tables the way the SED-ML parser builds one from a document. The problem which is written validates withpetabitself- what PEtab does not express goes into the
sbmlsimextension of the problem, i.e. the units of the parameters, the observables and the data, theFitSettings, 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 isrequired, 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 knowsbmlsimrejects the problem instead of fitting the same data with another objective without saying so, andto_petab(..., required_extension=False)writes a problem for those tools. A round trip throughsbmlsimkeeps 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 to7e-6 - the layer follows the PEtab v2 specification where it differs from v1: the standard deviation of a measurement is the
sdplaceholder an observable declares innoisePlaceholdersand not anoiseParameter${n}_${observableId}name, the parameter table has noparameterScaleany more, andsbmlsim.fit.petab_v2.symbolsconverts the selections of roadrunner into the math of a model instead of dropping the brackets:S1is the amount of a species and[S1]its concentration, while the identifier in the math of SBML is the amount forhasOnlySubstanceUnits=trueand the concentration forfalse, so the concentration of an amount based species is written asS1 / 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.gapsis the catalogue of what PEtab v2 cannot express about ansbmlsimfit, andgaps_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_dirno 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_EVOLUTIONreturned 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, sosizerepeats were one repeat,sizetimes. Every repeat gets its own seed, derived from the seed of the fit with anumpy.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_optimizationraised the number of repeats to the number of cores instead, i.e.,size=3with 12 cores ran 12 optimizations - under the
forkserverstart 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_squaresreturns, i.e., the residualsfun, the jacobianjac, the gradient and the active mask, which grow with the data and are not read anywhere; it keepsx,x0,cost,success,status,message,duration,optimalityandnfevnow. 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 itsoptimization_result.jsonfrom 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, andOptimizationResult.from_directoryreads them back into a result;OptimizationResult.run_resultandwrite_runare 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.displayrenders 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>fromfit_id, and it is the id of the optimization problem, of itsOptimizationResult, 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=Falseturns 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.pycreates the report of a finished fit from itsparameters.json, and compares the parameters of several fits in one reportexamples/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 ofsbmlsim.fit:python -m examples.hctz_fitting.simulationsruns the experiments,python -m examples.hctz_fitting.fitting.fittingthe fit, andop_hctzcreates the problem of a subset of the dataexamples/fit_sampling.py: the four sampling types of the start values of a fit, compared on a figureOptimizationAnalysisevaluates the residuals once per parameter vector and caches them. The analysis simulated every fit mapping eight times, once for each plot and tablesbmlsim.fit.rmse.msecompletesrmseandaic;SamplingType.is_logandis_lhs;sbmlsim.fit.runner.resolve_n_cores, which usesos.process_cpu_count, i.e., the cores the process may use;MappingMetaDatais exported fromsbmlsim.fit- every module of
sbmlsim.fitcarries 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.mdandCLAUDE.mddescribe 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/omexand the checked in outputs oftests/data/diff). The skipped teststests/comparison/test_diff.pyandtests/experiment/test_covid_examples.pyare gone with the code they waited for - the branch protection is part of the repository as rulesets in
.github/rulesets/, together withCODEOWNERSand a pull request template;maintracks the latest release and is fast-forwarded by thesync-mainjob of the release workflow, and the release is prepared on a branch and tagged ondevelopafter the merge, sobump-my-versionno 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.pybuilds the problem fromexamples.hctz_fitting,tests/fit/test_petab_omex.pycreates and reads an archive of the Boehm PEtab problem ofexamples/petab, and the analysis test asserts the files of the report examples.hctz_fitting.simulationsandexamples.fit_samplingare part of the example tests- the tests run in parallel with pytest-xdist,
pytesttakes 40 seconds on a workstation instead of 3.5 minutes;pytest -n 0runs them in one process helpers.fit_experiments_by_kindinitializes 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 fasterexamples/hctz_fitting/fitting/fitting.pyandreport.pyrun as scripts as well, i.e., through the "run file" of an IDE, which does not put the repository onsys.path; they failed withModuleNotFoundError: No module named 'examples'
Your sbmlsim team