diff --git a/CLAUDE.md b/CLAUDE.md index 3bb76eff..99dee8c0 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -56,7 +56,7 @@ Documentation is [Zensical](https://zensical.org/): markdown sources in `docs/`, **`units.py` - units.** One registry for the package, `sbmlsim.units.ureg`, and its quantity `sbmlsim.Q`; the registry adds `none`, `item` and `percent` to pint and no international unit, which is specific to a substance (a model or a user defines its own, e.g. `IU_insulin`). `UnitsInformation` maps identifiers to unit expressions; `from_sbml` reads the units of a model as expressions pint parses (`Units.udef_to_str`, prefixes and names of pint, e.g. `mmol/min`), the derived unit of a rule or a kinetic law is the first unit definition of the same unit (`unit_signature`: the kinds with their summed exponents and the factor of the multipliers and scales, instead of `libsbml.UnitDefinition.areIdentical`, which keeps about 11 to 31 kB per call for a unit with a multiplier, 13 MB per model, sbmlteam/libsbml#503), the ids of unit definitions are not defined in the registry except for a unit with a factor which is no prefix, under a collision free name. `sbmlsim.Q` takes the units of pint only (its `UndefinedUnitError` carries a note naming `model.Q` and `experiment.Q`); the unit ids of a model are kept in its `UnitsInformation` (`unit_ids`, unit id to expression, from `model_uid_dict`), and `UnitsInformation.quantity` replaces every whole identifier of a unit string which is a unit id by its expression in brackets (`replace_unit_ids`, a number like `1e-3` is no identifier, a unit id wins over a pint unit of the same name), so `mg_per_min/kg` works; `RoadrunnerSBMLModel.Q` is that function, `SimulationExperiment.Q(value, unit, model=None)` resolves the ids of the models the runner loaded before `initialize` (without `model=` an id must have one expression in every model which has it), and a `FitParameter.unit` may be a unit id of the model of its target (`ParameterMapping.parameter_units` resolves it in the models of the simulations it reaches and compares the units of versions as parsed units); the parameter keeps its unit for the display and the PEtab export, while the values of a fit carry the resolved units of pint (`OptimizationProblem.parameter_units`, the units of a `ParameterSet`, `OptimizationResult.units`, `FisherInformation.units`, the parameter set of an `IdentifiabilityResult`), which the designs of `simulation/sampling` read without a model. `normalize_changes` converts quantities to model units, and `Quantity` (an alias of pint's `PlainQuantity`, the base of every quantity a registry creates) is the type used in all annotations. -**`experiment/`, `data.py`, `task/`, `plot/`, `report/` - simulation experiments.** `SimulationExperiment` (`experiment/experiment.py`) is subclassed with `models()`, `datasets()`, `simulations()`, `observables()` (the `Formula`, `PK` and `Custom` of the scan core by their id), `tasks()`, `data()`, `figures()` and `fit_mappings()`, each returning a dict keyed by id; `initialize()` collects them and `run(simulator)` runs the tasks, evaluates the data and renders the figures. `ExperimentRunner` (`experiment/runner.py`) resolves the abstract models into roadrunner models, loaded once per `model_key` (the resolved path of a file, the content of inline SBML and, for a hierarchical one, its resolved base directory, the changes, the added parameters and the selections), runs several experiments and writes results, figures and the JSON serialization (`to_json`, encoded before the file is written, raises for a value `ObjectJSONEncoder` has no codec for; the run of an experiment reports it, names the type and writes no JSON instead of failing). A run keeps what the report needs and not more: an experiment keeps only the keys of its matplotlib figures once they are written (a figure holds the pixel buffer of its last rendering), and the runner releases the results of an experiment once its outputs are written (`keep_results=False`, `SimulationExperiment.results` raises then; `SimulationExperiment.run` keeps them by default), so a run of many experiments holds only the one it runs. `Data` (`data.py`) is a promise for data of a task (`Data("[X]", task=...)`, an observable `Data("pk.cmax", task=...)`), a dataset column or a function of other data; `DataSet` is a `pandas.DataFrame` with `uinfo`; `DataSet.from_df` brings the rows of a column to the unit of its first row and reads a unit with a numeric factor, which PK-DB writes (`ml/min/(1.73*m^2)`, `micromol/(0.5*hr)`), as a number of the values (`_parse_unit`, #294): the numbers of its rows are multiplied by the factor and the column gets the unit without it; a `*_unit` column applies to its column and its `_sd`, `_se`, `_min` and `_max`, the `unit` column to all numbers of a row (`mean`, `value`, `median`, `sd`, `se`, `min`, `max`, not `cv`), converted in one pass (`_unify_units`), so the statistics get one unit and shared numbers are converted once. `Task` pairs a model and a `Simulation` or a `Scan`; the results are `ScanResult`s, written as netCDF; the datasets of a run go as TSV into the folder `datasets` of its output path (`DATASETS_DIR`, `dataset_file`). The index of every task data of `data()`, the fit mappings and the figures is classified at `initialize()` as `time`, an observable, a coordinate of the scan or a selection of the model (an unknown one raises), and a task whose data read observables runs with the observables they need and keeps them next to its selections in one run. `Data.get_data` returns an `xarray.DataArray` with the dimensions of its source, their coordinates and `attrs["units"]` (`to_quantity` gives the pint quantity): a task variable over `(*dims, time)`, `(*dims, _point)` for a ragged result padded with `NaN`, or `(*dims)`, a dataset column over `row`, a function broadcast by dimension name; the plain name of a symbol is its timecourse, every array of a task names the values a dimension sets `.`, the plain name of an unkept changed target raises; `sel=` selects labels (a dimension of the scan the data has not is skipped) or the rows of a dataset by column values, and `sel` of task data is checked against the scan at `initialize()`. A curve names the dimensions of a scan it draws one line per point of (`over=`) and a band (`Plot.band`, `Band`) reduces a dimension of draws to two quantiles and the median; `plot/points.py` splits the broadcast arrays into the lines (`curve_lines`, `band_lines`, without the padding) and gives their colours (shades of the colour of the curve or viridis, a second dimension sets the line style), legend labels (the value and unit of the changed target) and from eleven points a colour bar, which both serializers draw; `initialize()` checks from the definitions that every dimension of a scan a curve or band reads is its axis, in `over`/`across` or selected, and a band of a ragged scan is interpolated onto the time points of a figure (`plot/points.py` `onto_grid`: `result/timecourse.py` `plot_times`, `PLOT_POINTS` 500 points, three quarters at the quantiles of all steps, a quarter equidistant, the times of the changes; the state after a change at its time, `NaN` outside of a simulation); its quantiles are `nan_quantile`, vectorized, since `numpy.nanquantile` reduces an array with a NaN in python per position. A log axis shows at most `LOG_DECADES` (6) below its largest value (`plotting.log_limits`, matplotlib and plotly; `limit_log_axis(ax)` for figures drawn with matplotlib directly), so the tiny first steps of the integrator do not stretch it; data within the decades keep their limits. `plot/plotting.py` is the backend independent figure model (`Figure`, `Plot`, `Axis`, `Curve`, `Band`, `ShadedArea`, `Style` with `Line`, `Marker`, `Fill`); `plot/serialization_matplotlib.py` renders it as the static images and `plot/serialization_plotly.py` as the interactive pages; the format chooses the backend, i.e. `figure_formats=["svg", "html"]`, `STATIC_FORMAT`/`INTERACTIVE_FORMAT` name the two and a run which asks only for `html` never builds a matplotlib figure. `FigureSettings` (`plot/plotting.py`) holds the settings of the rendering (dpi, panel size, fonts, legend), which cascade, the more specific wins: its defaults, the run (`figure_settings=` of `SimulationExperiment.run`, `ExperimentRunner.run_experiments` and `run_experiments`, checked when it starts), the experiment (the read-only class attribute `SimulationExperiment.figure_settings`, inherited by subclasses, checked at `initialize`) and the figure (`Figure(settings=...)`); both serializers read only the resolved settings (`Figure.resolved_settings`), a figure without an explicit `width`/`height` gets its size from the panels when it is rendered (`Figure.size`), the margins of a figure are the ones of a figure of one panel whatever its number of panels, the fractions of matplotlib of the size of a panel (`MatplotlibFigureSerializer._margins`, #323), grown by what the ticks and labels of the drawn panels need beyond them, at least `EDGE_PAD`, so a figure saved without a tight box keeps its labels, and a label longer than its axis is wrapped at its words (`_fit_margins`, `_wrap_label`, #362), every text is measured with one renderer of a pixel per resolution (`_measuring_renderer`; a text measured without a renderer drew the figure into a renderer of its size and kept it, 10 GB), and a figure the serializer built is saved with its tight box measured with it (`saved_box`, `bbox_inches` of `SimulationExperiment.save_mpl_figures`) instead of `"tight"`, which draws the figure once more to find it, and the plots of a figure with a title begin right below the band of the titles whatever its height (#312), a legend outside (`legend_position="outside"`) is placed right of its plot by `_place_right` together with the colour bars, its column gets the room of its widest legend and the figure is that much wider (no other gap grows; plotly has one legend of the figure), a legend taller than its panel gives the row the room of its overflow below it and the figure is that much taller, and the title of a widened figure is centered over the panels; `legend_scope` (`LEGEND_SCOPES`: `plot`, the default, and `row`, #327) collects the entries of the plots of a row into one legend right of the grid (`_row_legends`, a figure legend placed as a `_Slot` without a panel by `_place_right`, at the top of the cell of its row, right of the right y axes and colour bars of every panel of the last column, so the legends of all rows align, entries merged by `_entry_key`, label and style, as plotly merges them; `Plot.legend` decides whether a plot's entries enter; plotly gives each row a legend of its own, `legend`, `legend2`, ...), and the names of the settings on `Figure`, on a figure or in the body of a subclass raise an `AttributeError` with the migration (the settings were class attributes which leaked into later experiments, #260). plotly is in the `dev` extra and not a dependency; matplotlib keeps the static images because plotly writes them through a headless browser at 301 ms per figure against 117 ms. `report/experiment_report.py` renders HTML/markdown reports from jinja2 templates in `resources/templates/`; `report/templates.py` provides the shared jinja2 environment and `report_base.html`, the base the HTML reports of the experiments and of a fit extend, which carries the style, the search, the filter chips, the sortable tables and the figure overlay. +**`experiment/`, `data.py`, `task/`, `plot/`, `report/` - simulation experiments.** `SimulationExperiment` (`experiment/experiment.py`) is subclassed with `models()`, `datasets()`, `simulations()`, `observables()` (the `Formula`, `PK` and `Custom` of the scan core by their id), `tasks()`, `data()`, `figures()` and `fit_mappings()`, each returning a dict keyed by id; `initialize()` collects them and `run(simulator)` runs the tasks, evaluates the data and renders the figures. `ExperimentRunner` (`experiment/runner.py`) resolves the abstract models into roadrunner models, loaded once per `model_key` (the resolved path of a file, the content of inline SBML and, for a hierarchical one, its resolved base directory, the changes, the added parameters and the selections), runs several experiments and writes results, figures and the JSON serialization (`to_json`, encoded before the file is written, raises for a value `ObjectJSONEncoder` has no codec for; the run of an experiment reports it, names the type and writes no JSON instead of failing). A run keeps what the report needs and not more: an experiment keeps only the keys of its matplotlib figures once they are written (a figure holds the pixel buffer of its last rendering), and the runner releases the results of an experiment once its outputs are written (`keep_results=False`, `SimulationExperiment.results` raises then; `SimulationExperiment.run` keeps them by default), so a run of many experiments holds only the one it runs. A run of several experiments runs them in a pool (`run_experiments(n_workers=None)`: a worker per `EXPERIMENTS_PER_WORKER` (4) experiments, at most one per CPU, serial below and in a worker): a worker (`_run_job`) creates an experiment again from its class and the settings of the runner (`_Job`, pickled; a class defined in a function, at runtime or in the main script, which a worker imports without its guard, a shown figure, kept results or a unit registry of its own make the run serial, `_serial_reason`, and so does by default a simulator which runs the scans in a pool, `_scan_pool`), keeps the loaded models of its experiments in `parallel.worker_cache`, runs the experiment with the matplotlib parameters the run changed (`_changed_rc_params`) and a simulator of one worker, and sends back the fields of its `ExperimentResult` and its log records (`_Outcome`, collected by the only handler of the root logger of the worker while the experiment runs, `_Collect`), which the parent logs in the order of the experiments and releases the results of its experiment as a serial run does; the parent keeps its own initialized experiments for the report, submits the most expensive experiments first (`_cost`: the points of the scans and `FIGURE_COST` per figure), the forkserver preloads the modules of the experiments, and an experiment whose worker died runs in the parent again. `Data` (`data.py`) is a promise for data of a task (`Data("[X]", task=...)`, an observable `Data("pk.cmax", task=...)`), a dataset column or a function of other data; `DataSet` is a `pandas.DataFrame` with `uinfo`; `DataSet.from_df` brings the rows of a column to the unit of its first row and reads a unit with a numeric factor, which PK-DB writes (`ml/min/(1.73*m^2)`, `micromol/(0.5*hr)`), as a number of the values (`_parse_unit`, #294): the numbers of its rows are multiplied by the factor and the column gets the unit without it; a `*_unit` column applies to its column and its `_sd`, `_se`, `_min` and `_max`, the `unit` column to all numbers of a row (`mean`, `value`, `median`, `sd`, `se`, `min`, `max`, not `cv`), converted in one pass (`_unify_units`), so the statistics get one unit and shared numbers are converted once. `Task` pairs a model and a `Simulation` or a `Scan`; the results are `ScanResult`s, written as netCDF; the datasets of a run go as TSV into the folder `datasets` of its output path (`DATASETS_DIR`, `dataset_file`). The index of every task data of `data()`, the fit mappings and the figures is classified at `initialize()` as `time`, an observable, a coordinate of the scan or a selection of the model (an unknown one raises), and a task whose data read observables runs with the observables they need and keeps them next to its selections in one run. `Data.get_data` returns an `xarray.DataArray` with the dimensions of its source, their coordinates and `attrs["units"]` (`to_quantity` gives the pint quantity): a task variable over `(*dims, time)`, `(*dims, _point)` for a ragged result padded with `NaN`, or `(*dims)`, a dataset column over `row`, a function broadcast by dimension name; the plain name of a symbol is its timecourse, every array of a task names the values a dimension sets `.`, the plain name of an unkept changed target raises; `sel=` selects labels (a dimension of the scan the data has not is skipped) or the rows of a dataset by column values, and `sel` of task data is checked against the scan at `initialize()`. A curve names the dimensions of a scan it draws one line per point of (`over=`) and a band (`Plot.band`, `Band`) reduces a dimension of draws to two quantiles and the median; `plot/points.py` splits the broadcast arrays into the lines (`curve_lines`, `band_lines`, without the padding) and gives their colours (shades of the colour of the curve or viridis, a second dimension sets the line style), legend labels (the value and unit of the changed target) and from eleven points a colour bar, which both serializers draw; `initialize()` checks from the definitions that every dimension of a scan a curve or band reads is its axis, in `over`/`across` or selected, and a band of a ragged scan is interpolated onto the time points of a figure (`plot/points.py` `onto_grid`: `result/timecourse.py` `plot_times`, `PLOT_POINTS` 500 points, three quarters at the quantiles of all steps, a quarter equidistant, the times of the changes; the state after a change at its time, `NaN` outside of a simulation); its quantiles are `nan_quantile`, vectorized, since `numpy.nanquantile` reduces an array with a NaN in python per position. A log axis shows at most `LOG_DECADES` (6) below its largest value (`plotting.log_limits`, matplotlib and plotly; `limit_log_axis(ax)` for figures drawn with matplotlib directly), so the tiny first steps of the integrator do not stretch it; data within the decades keep their limits. `plot/plotting.py` is the backend independent figure model (`Figure`, `Plot`, `Axis`, `Curve`, `Band`, `ShadedArea`, `Style` with `Line`, `Marker`, `Fill`); `plot/serialization_matplotlib.py` renders it as the static images and `plot/serialization_plotly.py` as the interactive pages; the format chooses the backend, i.e. `figure_formats=["svg", "html"]`, `STATIC_FORMAT`/`INTERACTIVE_FORMAT` name the two and a run which asks only for `html` never builds a matplotlib figure. `FigureSettings` (`plot/plotting.py`) holds the settings of the rendering (dpi, panel size, fonts, legend), which cascade, the more specific wins: its defaults, the run (`figure_settings=` of `SimulationExperiment.run`, `ExperimentRunner.run_experiments` and `run_experiments`, checked when it starts), the experiment (the read-only class attribute `SimulationExperiment.figure_settings`, inherited by subclasses, checked at `initialize`) and the figure (`Figure(settings=...)`); both serializers read only the resolved settings (`Figure.resolved_settings`), a figure without an explicit `width`/`height` gets its size from the panels when it is rendered (`Figure.size`), the margins of a figure are the ones of a figure of one panel whatever its number of panels, the fractions of matplotlib of the size of a panel (`MatplotlibFigureSerializer._margins`, #323), grown by what the ticks and labels of the drawn panels need beyond them, at least `EDGE_PAD`, so a figure saved without a tight box keeps its labels, and a label longer than its axis is wrapped at its words (`_fit_margins`, `_wrap_label`, #362), every text is measured with one renderer of a pixel per resolution (`_measuring_renderer`; a text measured without a renderer drew the figure into a renderer of its size and kept it, 10 GB), and a figure the serializer built is saved with its tight box measured with it (`saved_box`, `bbox_inches` of `SimulationExperiment.save_mpl_figures`) instead of `"tight"`, which draws the figure once more to find it, and the plots of a figure with a title begin right below the band of the titles whatever its height (#312), a legend outside (`legend_position="outside"`) is placed right of its plot by `_place_right` together with the colour bars, its column gets the room of its widest legend and the figure is that much wider (no other gap grows; plotly has one legend of the figure), a legend taller than its panel gives the row the room of its overflow below it and the figure is that much taller, and the title of a widened figure is centered over the panels; `legend_scope` (`LEGEND_SCOPES`: `plot`, the default, and `row`, #327) collects the entries of the plots of a row into one legend right of the grid (`_row_legends`, a figure legend placed as a `_Slot` without a panel by `_place_right`, at the top of the cell of its row, right of the right y axes and colour bars of every panel of the last column, so the legends of all rows align, entries merged by `_entry_key`, label and style, as plotly merges them; `Plot.legend` decides whether a plot's entries enter; plotly gives each row a legend of its own, `legend`, `legend2`, ...), and the names of the settings on `Figure`, on a figure or in the body of a subclass raise an `AttributeError` with the migration (the settings were class attributes which leaked into later experiments, #260). plotly is in the `dev` extra and not a dependency; matplotlib keeps the static images because plotly writes them through a headless browser at 301 ms per figure against 117 ms. `report/experiment_report.py` renders HTML/markdown reports from jinja2 templates in `resources/templates/`; `report/templates.py` provides the shared jinja2 environment and `report_base.html`, the base the HTML reports of the experiments and of a fit extend, which carries the style, the search, the filter chips, the sortable tables and the figure overlay. **`fit/` - parameter fitting.** `FitMappingCollection` (the former `FitExperiment`) names a simulation experiment class, the mappings of it a fit uses together, their `MappingKind` and an `sid`; it is the unit a fit is defined in and the unit a PEtab problem is built from, `FitMapping` pairs a reference `FitData` (dataset) with an observable `FitData` (task), optionally with an `ObservableModel` when the observable is a formula of selections, `FitParameter` carries bounds, unit and an optional `Prior` (the priors of PEtab v2, truncated at the bounds, which `petab_v2.likelihood.log_prior` evaluates and the optimizer does not use), plus a `target` (the model entity it writes, `pid` by default) and a `mappings` selector so several parameters can estimate one target separately for different parts of the data, e.g. a dissolution rate once for oral and once for intravenous dosing; a selector must be a module level function, not a lambda or a closure, because `OptimizationProblem.__getstate__` pickles the parameters for the workers of a parallel fit. A fit mapping is a timecourse, a value per simulation (`FitData(xid=None)`, e.g. `pk.cmax`, against one or several rows of a reference) or values over one dimension of a scan (`xid="."`, interpolated along the values of the dimension), decided at `initialize` from the dimensions of the observable after `sel` (`OptimizationProblem.observation_kinds`); a simulation group of a scan computes the values and times of its points once (`simulator.scan_point_values`, which `scan_point_plans` shares) and simulates every point on the plan of an evaluation with the fitted values applied first (`worker.point_plan`), so a dimension wins over a fitted parameter of its target (`initialize` warns), and evaluates its observables with `simulator.worker.observe`, as the workers of a run do; a scan without dimensions is its simulation, while a group of timecourses of a plain simulation runs as before; the PEtab export raises for such mappings (the gaps `scalar-observable` and `scan-task`). `fit/parameter_mapping.py` holds `ParameterMapping`, which resolves the selectors to the simulation groups of an initialized problem, validates the binding (a selector must not split a simulation, two parameters must not write one target in one simulation, the versions of a target must agree on their unit) and reports the coverage of every parameter, which the console and the HTML report show as a table once a parameter is versioned. `OptimizationProblem` (`fit/optimization.py`) is picklable, it reads the start values of the parameters from the model as loaded (`RoadrunnerSBMLModel.loaded_value`, #315); `initialize(settings)` creates the runner, resolves the data and the weights and is a no-op when the `FitSettings` did not change, `residuals`/`cost_least_square` are the objective, `optimize` runs scipy least squares or differential evolution. `FitSettings` (`fit/options.py`, which also enumerates residual, loss, weighting and parameter scale types) is the bundle both a fit and its report need; `parameter_scale` (`LOG10` by default) is the space the optimizer searches, `problem.to_scale`/`from_scale` convert and the bounds and fitted parameters stay on the linear scale. `data_at_change` (`DataAtChangeType`, `before` by default) compares a data point at the time of a change after the start, within the resolution of the time of `merge_times` (`time_tolerance`, `Plan.at_change`), with the state before the change, except the first point of a curve (the point or points with the smallest time of a mapping which has data at a later time, `_first_after`; a curve of a single time takes the state before), which takes the state after it: the output of a fit is the one of a run, with two rows at a change after the start (an output time of the data at a change is the time of the change), and a data point takes one of them (`_rows` of a group which outputs the times of its data, whose output has the times of the changes within them as well; the mask `before` of `_interpolate_after_jumps`, the rule of `grid_weights`, for a group which keeps its output, which interpolates the data up to a change toward the state before it and from it from the state after it); a group which keeps its output of explicit times and has a timecourse outputs the times of the changes within them too (`_group_added`), computes its values per simulation on the output without them, as a run has it, and its timecourses on the output with them (`_timecourse_graph`), so its observables, e.g. the `PK` parameters, are the ones of the run, and the curve of a mapping in the report jumps at a change; a PEtab problem read without the extension and the PEtab suites use `after`, and the export reports the gap `measurement-at-change`. `fit/runner.py` samples the start points of all repeats (`fit/sampling.py`), so they do not depend on the number of workers, and runs the repeats serially in the calling process for one worker (the default, `n_cores=1`, or `serial=True`) and otherwise as tasks of a pool of `sbmlsim.parallel.start_pool` which the fit stops after it, a worker initializes the problem once per fit and keeps it in `parallel.worker_cache`; `FitPool` keeps at most `n_cores` repeats in flight and collects them as they finish, ordered by their index, a pool which breaks because a worker dies is started again (at most `MAX_POOL_RESTARTS` times) and the repeats which were running, the suspects, run again one at a time, so only a repeat which breaks the pool alone fails; it shows the progress of the runs with rich. A parallel fit needs the `if __name__ == "__main__":` guard, the runner reports it when the workers die while they start; `fit/result.py` collects the fits into an `OptimizationResult`, which stores the settings and hands out the fitted parameters as the `ParameterSets` of `fit/parameters.py`. Reporting is separate from optimizing: `FitReport` (`fit/report.py`) is built from a problem, its settings and one or more parameter sets, so a report is created from stored parameters without running a fit and several sets are compared in one report. The figure of every fit mapping (`FitReport.plot_mapping`: its curve, linear and logarithmic, and its goodness of fit and Bland-Altman plot, drawn by the panel functions of the summary figures, `_draw_goodness_of_fit` and `_draw_bland_altman`, with the metrics of the mapping) is most of the cost of a report, `mapping_figures=False` (and `--no-mapping-figures` on the command line tools) leaves them out, which is what the tests do. Its `index.html` is rendered from the jinja2 template `resources/templates/fit_report.html` into the sections overview, results and fit mappings; the page carries its own css and javascript (search, filter chips, sortable tables, click to enlarge) and loads nothing from the network. The `MappingKind` of a `FitExperiment` decides what a fit does with the mappings it selects: `TRAINING` enters the cost, `VALIDATION` (data a fit is checked against) and `OUTLIER` (unusable data) are resolved, simulated and evaluated but not fitted, `EXCLUDED` (data the model does not describe) is not resolved at all; `EVALUATED_KINDS` and `UNUSED_KINDS` are the two halves, and `training_indices`/`validation_indices`/`outlier_indices` are the positions in the resolved data. The kind belongs to the selection of the data, not to the fit mappings of a simulation experiment, whose `MappingMetaData` describes the curve; `fit/helpers.py` is the selection: `FitMappings` is the complete list of the fit mappings of the experiments, loaded once, and `FitMappings.select(filters, outliers, validation)` sets the kind of every mapping in three ordered steps, the filters select the training data and exclude the rest, the outlier keys tag the training data which is not usable, and the validation keys or filter take part of what is left, so excluded beats outlier beats validation beats training; the `MappingSelection` holds the kinds, the overview table and one `FitMappingCollection` per experiment and kind, `select_mapping_collections` is the one call version a `FitDefinition` uses. `fit/metrics.py` (the former `rmse.py`) calculates the metrics of a parameter set on a problem, i.e., the `DV`/`PRED`/`IPRED` table with the residuals and MSE, RMSE, R² and AIC per mapping (keyed by experiment and mapping), per kind (training, validation and outlier) and overall, the AIC and BIC of a perfect fit are `-inf` and R² of fewer than two points `NaN`; a report writes them as TSV. `FitReport.points`/`point_kinds` are the data points with their kind and the kinds the problem has, which `plot_goodness_of_fit` (prediction against measurement) and `plot_bland_altman` (the ratio over the geometric mean with the bias and the limits of agreement of `agreement`, i.e. of the training data, as fold factors and the same lines in every panel) draw one panel each for; both draw that band in the same styles (`IDENTITY_STYLE`, `BIAS_STYLE`, `LIMITS_STYLE`, `BAND_ALPHA`), horizontal in the one and parallel to the diagonal in the other; there is no `all` panel and no `all` row in `summary_df`, which would pool fitted with dropped data. The points of both carry the color of their study (`studies`/`study_color`, `STUDY_COLORS`) and the marker of their parameter set (`set_marker`, `SET_MARKERS`). A report only shows what a fit evaluates, i.e. `EVALUATED_KINDS`; the excluded data has no rows, no chip and no figures. `fit/display.py` renders the console output of a fit as sections (the problem, the parameters, the settings and the data with its kinds), so that the runner, the command line tools and an interactive session look the same. `fit/cli.py` holds what every fit shares: a `FitDefinition` is the fit experiments, the parameters, the paths and the settings of a problem, `run_fit` builds the problems for an `OptimizationStrategy` and runs them, and `fit_cli`/`report_cli` are the command line tools a model reuses, so `examples/hctz_fitting/fitting/fitting.py` is only the definitions. `fit/identifiability.py` is the profile likelihood: `profile_likelihood(problem, settings, parameter_set)` scans every parameter around the optimum in the space of the `parameter_scale` with adaptive steps, re-optimizes the other parameters at every point in the worker pool of the runner, and classifies the parameters from the shape of the profile against the threshold `cost_min + chi2.ppf(alpha, df) / 2` (the cost is `0.5 * Σ r²`); `FitReport(identifiability=...)` adds the section with the profiles and `cli.identifiability_cli` is the tool for stored parameters. `fit/fisher.py` is the local analysis next to it: `fisher_information` builds `FIM = J'J` from one jacobian of the weighted residuals in the space of the `parameter_scale` and gives the standard errors, the correlations, the confidence intervals `θ ± t·SE`, the eigenvalues, the rank and the condition number, so a direction the data does not constrain is found without scanning; `FitRun.fisher()` computes it for a finished fit and `FitReport(fisher=...)` writes it into the same identifiability section, which `identifiability_cli` does unless `--no-fisher` is given. Pickling an initialized `OptimizationProblem` gives its uninitialized definition (`__getstate__`), which is what the workers need. `fit/petab_omex.py` packages PEtab problems as archives. `fit/petab_v2/` is the PEtab v2 layer: `export.py` writes an `OptimizationProblem` as a PEtab v2 problem (models, conditions/experiments from the timecourses, observables and measurements from the fit mappings, parameters, a versioned parameter as a condition which assigns its target the value of the estimated parameter), `reader.py` reads one back and builds a `SimulationExperiment` subclass at runtime, resolving such a condition into a `FitParameter` selected by `helpers.filter_keys`, i.e. by the mapping ids the condition covers rather than by the selector which chose them, `extension.py` is the `sbmlsim` block of the problem YAML which carries what PEtab drops (units, `FitSettings`, mapping kinds, mapping weights, the output grid) so a round trip is exact, an observable of the problem which is a formula rather than an entity is read as an `ObservableModel` (`fit/objects.py`, the formula in the selections of roadrunner with placeholders and their values per measurement), which the fit evaluates on the simulation at the data, so the model is not rewritten, and `gaps.py` is the catalogue of the differences with `gaps_of_problem` reporting the ones a given problem hits. `testsuite.py` runs the PEtab v2 test suite (a pinned commit, cached like the SBML Test Suite, `SBMLSIM_PETAB_SUITE_PATH`): a case of a model is read, simulated at the nominal values of its parameter table and compared on llh, chi2, simulations, log prior and log posterior, a math case compiles its expression with `compile_formula`; `scripts/petab_testsuite.py` writes and checks `tests/data/petab_baseline.json` (with `sbmlsim.testsuite.baseline`, which the SciML suite shares) and the cases are tests with the marker `petab_testsuite`. `benchmark.py` is the same for the benchmark collection (a pinned commit converted with `petab1to2`, `SBMLSIM_BENCHMARK_PATH`): a problem is simulated once at its nominal values and compared with the simulations of the collection, matched by row or by their keys, and with the log-likelihood of AMICI; `scripts/petab_benchmark.py` runs the problems in parallel with a time limit and keeps `tests/data/benchmark_baseline.json`, marker `petab_benchmark`. A problem which is read builds its `SimulationExperiment` at runtime as a subclass of `reader.PetabExperiment` bound to its `_PetabExperimentDefinition` (the tables, the paths and the compiled models, without the reader), which pickles as that definition and builds the reader again in a worker, so such a fit runs in parallel as well. `examples/hctz_fitting` is the reference problem, `tests/fit/` runs against it. diff --git a/docs/experiments.md b/docs/experiments.md index 9e9ae7ce..3729a419 100644 --- a/docs/experiments.md +++ b/docs/experiments.md @@ -93,6 +93,8 @@ res = experiment.results["task_tc"] print(res["[X]"].values[-3:]) ``` +The experiments of a run are independent, so a run of several experiments runs them in a pool of worker processes: `run_experiments(n_workers=None)` starts a worker per four experiments, at most one per CPU, and runs a few experiments in the calling process; `n_workers=1` runs every experiment in the calling process and a number in a pool of that many workers. A worker creates an experiment again from its class with the settings of the runner, loads its models once for all its experiments, runs it with the parameters of matplotlib of the run and writes its outputs; the messages of the workers are logged in the order of the experiments, the most expensive experiments (the points of their scans and their figures) start first. An experiment runs its scans serially in its worker, so by default a runner whose simulator runs its scans in a pool (`Simulator(n_workers=4)`, or `n_workers=None` and a scan of `POOL_THRESHOLD` points) runs its experiments in the calling process, and its scans get the pool of the simulator; a number of workers runs the experiments in a pool all the same. The 28 experiments of the HCTZ model of pkdb_models take 12 s instead of 21 s. A run whose figures are shown, whose results are kept (`keep_results=True`), whose runner has a unit registry of its own or whose experiments cannot be pickled or imported by a worker, e.g. a class defined in a function or in the main script, runs in the calling process: a worker imports the main script without the block of its guard, so an experiment is defined in a module to run in a pool. The workers import the script again, so a script which runs a pool needs the guard `if __name__ == "__main__":`; a worker which dies, e.g. by a crash of roadrunner, has its experiment run in the calling process again. + `run_experiments(reduced_selections=True)` records for every task only the variables its data refer to, which speeds up large experiments; `reduced_selections=False` records every variable of the model and the ones the data of the task refer to. One failing experiment does not stop a run. An error in the definition of an experiment (its models, datasets, simulations, tasks, data, figures or fit mappings, or a key which is no valid SId) is logged with its traceback when the runner is created, and the experiment is not run; an error while an experiment runs ends that experiment; a failing figure is skipped and the other figures are written. The run ends with a summary of the failures in the log, every `ExperimentResult` records them (`error`, `failed_figures` and `failed`) and the report lists them. `run_experiments(raise_on_failure=True)` raises one `ExperimentRunError` after every experiment ran, `ExperimentRunner(on_error="raise")` raises the error of a definition at once, which is what a fit does, and `SimulationExperiment.run` called directly raises its error. A figure which failed, or which a failed experiment did not write, leaves no file of an earlier run in the output directory, so no figure there looks like one of this run. diff --git a/src/sbmlsim/experiment/runner.py b/src/sbmlsim/experiment/runner.py index 26415eb3..02bb261e 100644 --- a/src/sbmlsim/experiment/runner.py +++ b/src/sbmlsim/experiment/runner.py @@ -6,16 +6,35 @@ - loading of models - running tasks (simulation on models) - creating outputs + +The experiments of a run are independent of each other, a run of several +experiments runs them in a pool of worker processes (`n_workers` of +`ExperimentRunner.run_experiments`): every worker initializes an experiment +again from its class, loads its models once and keeps them for the next +experiment, runs it serially (its scans included, a worker starts no pool) and +writes its outputs; the parent keeps its own initialized experiments, which +the report is created from, and logs the messages of the workers in the order +of the experiments. A run is serial when an experiment cannot be sent to a +worker, see `ExperimentRunner._serial_reason`. """ +import copy import dataclasses import inspect import logging +import math +import os +import pickle import time from collections.abc import Iterable, Mapping +from concurrent.futures import Future +from concurrent.futures.process import BrokenProcessPool from pathlib import Path -from typing import Any, Literal, cast +from typing import Any, Literal, cast, override + +import matplotlib +from sbmlsim import parallel from sbmlsim.console import console from sbmlsim.experiment.experiment import ( DATASETS_DIR, @@ -32,6 +51,7 @@ ) from sbmlsim.plot import FigureSettings from sbmlsim.report.experiment_report import ExperimentReport, ReportResults +from sbmlsim.simulation import Scan from sbmlsim.simulator import Simulator from sbmlsim.units import UnitRegistry from sbmlsim.units import ureg as package_ureg @@ -39,6 +59,20 @@ logger = logging.getLogger(__name__) +#: the experiments per worker of a pool of `run_experiments` for +#: `n_workers=None`: a worker imports sbmlsim and the modules of the +#: experiments and loads their models before its first experiment, about 2 s, +#: which it wins back over a few experiments. The 28 experiments of the HCTZ +#: model of pkdb_models take 20.8 s serially, 12.9 s on 4 workers, 11.9 s on +#: 8 and 12.8 s on 20, where the most expensive experiment is most of the time +#: and every worker keeps 330 MB. +EXPERIMENTS_PER_WORKER: int = 4 + +#: the cost of a figure in simulations, which orders the experiments of a pool +#: by their cost, see `_cost`: drawing and writing a figure takes about as +#: long as 30 simulations of a model of a whole body +FIGURE_COST: int = 30 + #: the model of an experiment, i.e. what decides whether two experiments can #: share one loaded roadrunner instance @@ -148,6 +182,9 @@ def __init__( self.simulator: Simulator | None = None # the sids of the experiments in the order of their classes self._order: list[str] = [] + # the settings every experiment was created with, which a worker + # creates it with again + self._settings: dict[str, dict[str, Any]] = {} classes: list[type[SimulationExperiment]] = ( list(experiment_classes) @@ -203,6 +240,7 @@ def initialize( sid = experiment.sid if sid not in self._order: self._order.append(sid) + self._settings[sid] = kwargs # an experiment of the same sid replaces an earlier one self.experiments.pop(sid, None) self.failed.pop(sid, None) @@ -260,6 +298,7 @@ def run_experiments( keep_results: bool = False, raise_on_failure: bool = False, figure_settings: FigureSettings | Mapping[str, Any] | None = None, + n_workers: int | None = None, ) -> list[ExperimentResult]: """Run the experiments and write their outputs. @@ -281,15 +320,27 @@ def run_experiments( a `FigureSettings` or overrides of its defaults, which the `figure_settings` of an experiment and the settings of a figure override. + n_workers: the processes which run the experiments: `1` runs them + in the calling process, a number in a pool of that many + workers (at most one per experiment), `None` in a pool of a + worker per `EXPERIMENTS_PER_WORKER` experiments, at most one + per CPU, i.e. serially for a few experiments, and serially in a + worker. In a pool every experiment + runs serially in its worker, its scans included; a run with + `show_figures`, `keep_results`, a unit registry of its own or + an experiment which cannot be pickled is serial, see the + module. The pool starts worker processes, a script which runs + it needs the guard `if __name__ == "__main__":`. Returns: The results of the experiments, which the report is created from. Raises: ExperimentRunError: with `raise_on_failure`, if something failed. - ValueError: if a key of `figure_settings` is no setting, before - any experiment runs. + ValueError: if a key of `figure_settings` is no setting or + `n_workers` is less than 1, before any experiment runs. """ + parallel.check_workers(n_workers) settings = ( FigureSettings.resolve(run=figure_settings) if figure_settings is not None @@ -303,36 +354,202 @@ def run_experiments( output_path.mkdir(parents=True) formats = figure_formats if figure_formats is not None else [STATIC_FORMAT] - exp_results = [] - for sid in self._order: - console.rule(style="white") - if sid in self.failed: - # the definition raised, the experiment is reported, not run + options: dict[str, Any] = { + "show_figures": show_figures, + "save_results": save_results, + "reduced_selections": reduced_selections, + "keep_results": keep_results, + "figure_settings": settings, + } + jobs = self._jobs(output_path, formats, figure_formats, options) + workers = self._workers(n_workers, jobs, options) + executor = ( + parallel.start_pool(workers, preload=_modules(jobs.values())) + if workers > 1 + else None + ) + try: + # the most expensive experiments first, so that none of them + # starts last and keeps the others waiting + futures: dict[str, Future[_Outcome]] = ( + { + sid: executor.submit(_run_job, jobs[sid]) + for sid in sorted( + jobs, key=lambda sid: -_cost(self.experiments[sid]) + ) + } + if executor is not None + else {} + ) + exp_results = [] + for sid in self._order: + console.rule(style="white") + if sid in self.failed: + # the definition raised, the experiment is reported, not run + exp_results.append( + self._failed_result( + self.failed[sid], output_path / sid, formats + ) + ) + continue + logger.info("Running SimulationExperiment: '%s'", sid) + # ExperimentResult used to create report; an error of the + # experiment or of a figure is logged and recorded in it + outcome = self._outcome(futures.get(sid)) exp_results.append( - self._failed_result(self.failed[sid], output_path / sid, formats) + self._run_experiment( + sid, + output_path / sid, + formats, + figure_formats=figure_formats, + **options, + ) + if outcome is None + else outcome.result(self.experiments[sid], output_path / sid) ) - continue - logger.info("Running SimulationExperiment: '%s'", sid) - # ExperimentResult used to create report; an error of the experiment - # or of a figure is logged and recorded in it - exp_results.append( - self._run_experiment( - sid, - output_path / sid, - formats, - show_figures=show_figures, - save_results=save_results, - figure_formats=figure_formats, - reduced_selections=reduced_selections, - keep_results=keep_results, - figure_settings=settings, - ) - ) + finally: + if executor is not None: + parallel.stop(executor) self._log_summary(exp_results) if raise_on_failure and any(r.failed for r in exp_results): raise ExperimentRunError(exp_results) return exp_results + def _jobs( + self, + output_path: Path, + formats: list[str], + figure_formats: list[str] | None, + options: dict[str, Any], + ) -> dict[str, "_Job"]: + """Get what a worker needs to run every experiment, see `_run_job`.""" + rc, levels = _changed_rc_params(), _logger_levels() + return { + sid: _Job( + experiment_class=type(self.experiments[sid]), + settings=self._settings[sid], + base_path=self.base_path, + data_path=self.data_path, + jit_optimization=self.jit_optimization, + on_error=self.on_error, + simulator=self.simulator, + output_path=output_path / sid, + formats=formats, + figure_formats=figure_formats, + options=options, + rc=rc, + levels=levels, + ) + for sid in self._order + if sid in self.experiments + } + + def _workers( + self, n_workers: int | None, jobs: dict[str, "_Job"], options: dict[str, Any] + ) -> int: + """Get the number of processes which run the experiments, see `run_experiments`.""" + if n_workers is None: + if parallel.in_worker(): + return 1 + n_workers = min( + os.process_cpu_count() or 1, + math.ceil(len(jobs) / EXPERIMENTS_PER_WORKER), + ) + if n_workers > 1 and (reason := self._scan_pool()) is not None: + logger.info("The experiments run in this process: %s", reason) + return 1 + n_workers = min(n_workers, len(jobs)) + if n_workers <= 1: + return 1 + reason = self._serial_reason(jobs, options) + if reason is not None: + logger.info("The experiments run in this process: %s", reason) + return 1 + return n_workers + + def _scan_pool(self) -> str | None: + """Get why the scans need the pool of the simulator, `None` if they do not. + + An experiment runs its scans serially in a worker: a simulator which + runs every scan in a pool of its own (`n_workers` above one), or a + scan which the simulator runs in a pool of every CPU (`n_workers=None` + and `parallel.POOL_THRESHOLD` points), is faster in this process. + """ + if self.simulator is None or self.simulator.n_workers == 1: + return None + if self.simulator.n_workers is not None: + return ( + f"the simulator runs the scans in a pool of " + f"{self.simulator.n_workers} workers" + ) + largest = max( + ( + len(simulation) + for experiment in self.experiments.values() + for simulation in experiment._simulations.values() + if isinstance(simulation, Scan) + ), + default=0, + ) + if largest >= parallel.POOL_THRESHOLD: + return f"the simulator runs a scan of {largest} points in a pool" + return None + + def _serial_reason( + self, jobs: dict[str, "_Job"], options: dict[str, Any] + ) -> str | None: + """Get why the experiments cannot run in workers, `None` if they can. + + A worker creates an experiment from its class, which pickle sends by + its module and name: a class defined in a function or at runtime, or + in the main module of an interactive session, which a worker cannot + import, is not found; nor is one of the main script which is defined + behind its guard `if __name__ == "__main__":`, which the worker does + not run, so an experiment of the main script runs here. The figures + shown and the results kept are in the worker, and a worker has the + unit registry of the package. + """ + if parallel.in_worker(): + return "a worker starts no pool" + if options["show_figures"]: + return "the figures are shown" + if options["keep_results"]: + return "the results of the experiments are kept" + if self.ureg is not package_ureg: + return "the runner has a unit registry of its own" + for job in jobs.values(): + cls = job.experiment_class + if cls.__module__ == "__main__": + return ( + f"'{cls.__name__}' is defined in the main script, define it " + f"in a module to run it in a worker" + ) + try: + pickle.dumps(job) + except Exception as err: + return f"'{cls.__name__}' cannot be sent to a worker: {err}" + return None + + @staticmethod + def _outcome(future: "Future[_Outcome] | None") -> "_Outcome | None": + """Get the outcome of an experiment in a worker and log its messages. + + `None` without a worker or when the worker died, e.g. by a crash of + roadrunner; the experiment then runs in this process. + """ + if future is None: + return None + try: + outcome = future.result() + except BrokenProcessPool: + logger.error( + "The worker of the experiment died, it runs in this process again" + ) + return None + for record in outcome.records: + logging.getLogger(record.name).handle(record) + return outcome + def _run_experiment( self, sid: str, @@ -531,6 +748,204 @@ def _log_summary(results: list[ExperimentResult]) -> None: ) +@dataclasses.dataclass(frozen=True) +class _Job: + """What a worker needs to run an experiment, see `_run_job`.""" + + experiment_class: type[SimulationExperiment] + #: the settings the experiment is created with + settings: dict[str, Any] + base_path: Path | None + data_path: Path | Iterable[Path] | None + jit_optimization: JitOptimization + on_error: Literal["raise", "log"] + simulator: Simulator | None + output_path: Path + formats: list[str] + figure_formats: list[str] | None + #: the other arguments of `SimulationExperiment.run` + options: dict[str, Any] + #: the parameters of matplotlib which the run changed, see `_changed_rc_params` + rc: dict[Any, Any] + #: the levels of the loggers of the run, see `_logger_levels` + levels: dict[str, int] + + +@dataclasses.dataclass(frozen=True) +class _Outcome: + """What a worker sends back of an experiment: its result and messages.""" + + error: str | None + failed_figures: dict[str, str] + figures: dict[str, list[str]] + datasets: list[str] + records: list[logging.LogRecord] + + def result( + self, experiment: SimulationExperiment, output_path: Path + ) -> ExperimentResult: + """Get the result of the experiment of the run. + + Its results stayed in the worker, so they are released, as the ones of + an experiment of a serial run are. + """ + experiment._results_released = True + return ExperimentResult( + experiment=experiment, + output_path=output_path, + error=self.error, + failed_figures=self.failed_figures, + figures=self.figures, + datasets=self.datasets, + ) + + +class _Collect(logging.Handler): + """Collect the messages of a worker, which the parent logs, see `_Outcome`.""" + + def __init__(self) -> None: + """Collect every message the loggers pass on.""" + super().__init__(logging.NOTSET) + self.records: list[logging.LogRecord] = [] + + @override + def emit(self, record: logging.LogRecord) -> None: + """Keep the message with its arguments and traceback as text, which pickle sends. + + A message with attributes which cannot be pickled, e.g. an `extra`, + is kept with the attributes of a plain message. + """ + record.msg = record.getMessage() + record.args = None + if record.exc_info: + record.exc_text = logging.Formatter().formatException(record.exc_info) + record.exc_info = None + try: + pickle.dumps(record) + except Exception: + plain = logging.LogRecord( + record.name, + record.levelno, + record.pathname, + record.lineno, + record.msg, + None, + None, + record.funcName, + ) + plain.exc_text = record.exc_text + plain.created, plain.msecs = record.created, record.msecs + record = plain + self.records.append(record) + + +def _cost(experiment: SimulationExperiment) -> int: + """Estimate the cost of an experiment: its simulations and figures. + + Every point of a scan is a simulation, a figure costs `FIGURE_COST`. + """ + points = 0 + for task in experiment._tasks.values(): + simulation = experiment._simulations.get(task.simulation_id) + points += len(simulation) if isinstance(simulation, Scan) else 1 + return points + FIGURE_COST * len(experiment._figures) + + +def _modules(jobs: Iterable[_Job]) -> list[str]: + """Get the modules which the forkserver imports once for all workers. + + The runner and the modules of the experiments, with what they import + (sbmlsim, matplotlib, the packages of a model); the main module is not + preloaded, see `sbmlsim.parallel`. + """ + modules = {job.experiment_class.__module__ for job in jobs} - {"__main__"} + return sorted({__name__, *modules}) + + +def _changed_rc_params() -> dict[Any, Any]: + """Get the parameters of matplotlib which differ from the ones as loaded. + + A caller sets them for the figures of a run, e.g. with `rc_context`; the + worker loads the same `matplotlibrc` and sets them again. The backend is + the one of the worker. + """ + # `copy` reads the values without the warnings of deprecated parameters + current, loaded = matplotlib.rcParams.copy(), matplotlib.rcParamsOrig.copy() + return { + key: value + for key, value in current.items() + if key != "backend" and value != loaded.get(key) + } + + +def _logger_levels() -> dict[str, int]: + """Get the levels of the loggers which have one, the root logger included.""" + loggers = logging.root.manager.loggerDict.values() + levels = { + logger.name: logger.level + for logger in loggers + if isinstance(logger, logging.Logger) and logger.level != logging.NOTSET + } + levels[""] = logging.root.level + return levels + + +def _run_job(job: _Job) -> _Outcome: + """Run an experiment in a worker, see the module. + + The models of the experiments of the worker are kept for the next + experiment, as the runner of the parent keeps them. + """ + models: dict[ModelKey, RoadrunnerSBMLModel] = parallel.worker_cache( + ("experiment models", job.jit_optimization), dict + ) + for name, level in job.levels.items(): + logging.getLogger(name or None).setLevel(level) + # the messages go to the run only: a handler of the worker, e.g. of a + # script which configures logging when it is imported, would print them + # a second time and out of order + handler = _Collect() + handlers = logging.root.handlers[:] + logging.root.handlers[:] = [handler] + simulator = copy.copy(job.simulator) + if simulator is not None: + # a worker runs the scans of its experiment itself + simulator.n_workers = 1 + try: + with matplotlib.rc_context(job.rc): + runner = ExperimentRunner( + [], + base_path=job.base_path, + data_path=job.data_path, + simulator=simulator, + on_error=job.on_error, + jit_optimization=job.jit_optimization, + ) + runner.models = models + runner.initialize([job.experiment_class], **job.settings) + (sid,) = runner._order + result = ( + runner._failed_result(runner.failed[sid], job.output_path, job.formats) + if sid in runner.failed + else runner._run_experiment( + sid, + job.output_path, + job.formats, + figure_formats=job.figure_formats, + **job.options, + ) + ) + finally: + logging.root.handlers[:] = handlers + return _Outcome( + error=result.error, + failed_figures=result.failed_figures, + figures=result.figures, + datasets=result.datasets, + records=handler.records, + ) + + def run_experiments( experiments: type[SimulationExperiment] | list[type[SimulationExperiment]], output_path: Path, @@ -539,6 +954,7 @@ def run_experiments( show_report: bool = False, raise_on_failure: bool = False, figure_settings: FigureSettings | Mapping[str, Any] | None = None, + n_workers: int | None = None, ) -> list[ExperimentResult]: """Run simulation experiments and write their report to the output path. @@ -553,6 +969,8 @@ def run_experiments( written. figure_settings: the settings of the figures of every experiment, see `ExperimentRunner.run_experiments`. + n_workers: the processes which run the experiments, see + `ExperimentRunner.run_experiments`. Returns: The results of the experiments. @@ -574,6 +992,7 @@ def run_experiments( output_path=output_path, show_figures=False, figure_settings=figure_settings, + n_workers=n_workers, ) report_results = ReportResults() for exp_result in results: diff --git a/tests/experiment/test_runner_parallel.py b/tests/experiment/test_runner_parallel.py new file mode 100644 index 00000000..7b536cd2 --- /dev/null +++ b/tests/experiment/test_runner_parallel.py @@ -0,0 +1,273 @@ +"""Tests of the experiments of a run in a pool of worker processes.""" + +import json +import logging +import os +import pickle +import sys +from pathlib import Path +from typing import Any, override + +import matplotlib +import pytest + +from sbmlsim import parallel +from sbmlsim.experiment import ExperimentResult, ExperimentRunner, SimulationExperiment +from sbmlsim.experiment.runner import EXPERIMENTS_PER_WORKER, _Collect, _run_job +from sbmlsim.simulator import Simulator +from tests.experiment.test_runner_errors import ( + BrokenTask, + GoodExperiment, + _bad_figure, + _Base, + _figure, +) + +#: the file in which an experiment records where it ran +PROCESS_FILE = "process.json" + + +class _Recorded(_Base): + """An experiment which records its process and a setting of matplotlib.""" + + @override + def run(self, *args: Any, **kwargs: Any) -> ExperimentResult: + result = super().run(*args, **kwargs) + output_path = kwargs["output_path"] + (output_path / PROCESS_FILE).write_text( + json.dumps( + { + "pid": os.getpid(), + "linewidth": matplotlib.rcParams["lines.linewidth"], + } + ) + ) + return result + + +class First(_Recorded): + """An experiment which works.""" + + +class Second(_Recorded): + """An experiment which works.""" + + +class Third(_Recorded): + """An experiment which works.""" + + +class BadFigure(_Recorded): + """An experiment with a figure which fails.""" + + @override + def figures(self) -> dict: + ok, bad = _figure("fig_ok"), _bad_figure() + ok.experiment = bad.experiment = self + return {"fig_ok": ok, "fig_bad": bad} + + +class DiesInAWorker(_Recorded): + """An experiment whose worker dies, it runs again in the parent.""" + + @override + def run(self, *args: Any, **kwargs: Any) -> ExperimentResult: + if parallel.in_worker(): + os._exit(1) + return super().run(*args, **kwargs) + + +CLASSES: list[type[SimulationExperiment]] = [ + First, + BrokenTask, + Second, + BadFigure, + GoodExperiment, + Third, +] + + +def _runner( + *classes: type[SimulationExperiment], simulator: Simulator | None = None +) -> ExperimentRunner: + return ExperimentRunner( + experiment_classes=list(classes), + simulator=simulator or Simulator(n_workers=1), + base_path=Path("."), + data_path=Path("."), + ) + + +def _process(path: Path, sid: str) -> dict[str, Any]: + return json.loads((path / sid / PROCESS_FILE).read_text()) + + +def _outcome(result: ExperimentResult) -> tuple[Any, ...]: + return ( + result.experiment.sid, + result.error, + sorted(result.failed_figures), + result.figures, + result.datasets, + ) + + +def _run( + path: Path, + n_workers: int | None, + classes: list[type[SimulationExperiment]], + simulator: Simulator | None = None, +) -> list[ExperimentResult]: + with matplotlib.rc_context({"lines.linewidth": 4.25}): + return _runner(*classes, simulator=simulator).run_experiments( + output_path=path, n_workers=n_workers + ) + + +def test_the_experiments_of_a_pool_give_the_results_of_a_serial_run( + tmp_path: Path, caplog: pytest.LogCaptureFixture +) -> None: + """The experiments run in workers, with the settings of matplotlib of the run.""" + with caplog.at_level(logging.INFO): + pooled = _run(tmp_path / "pooled", 3, CLASSES) + serial = _run(tmp_path / "serial", 1, CLASSES) + + assert [_outcome(r) for r in pooled] == [_outcome(r) for r in serial] + assert [r.experiment.sid for r in pooled] == [c.__name__ for c in CLASSES] + # the experiments of the run, not copies of a worker, are in the results + assert all( + isinstance(r.experiment, c) for r, c in zip(pooled, CLASSES, strict=True) + ) + for sid in ("First", "Second", "Third", "BadFigure"): + assert _process(tmp_path / "pooled", sid)["pid"] != os.getpid() + assert _process(tmp_path / "serial", sid)["pid"] == os.getpid() + assert _process(tmp_path / "pooled", sid)["linewidth"] == 4.25 + for fig_format in ("svg",): + assert (tmp_path / "pooled" / sid / f"{sid}_fig_ok.{fig_format}").exists() + # the messages of the workers reach the logging of the run, in its order + messages = [r.getMessage() for r in caplog.records] + running = [m for m in messages if m.startswith("Running SimulationExperiment")] + assert running == [f"Running SimulationExperiment: '{c.__name__}'" for c in CLASSES] + assert any( + r.levelno == logging.ERROR and "fig_bad" in r.getMessage() + for r in caplog.records + ) + assert "experiment broken" in caplog.text + assert any("2 of 6 experiments failed" in m for m in messages) + + +def test_an_experiment_whose_worker_dies_runs_in_the_parent(tmp_path: Path) -> None: + results = _run(tmp_path, 2, [First, DiesInAWorker, Second]) + + assert [r.error for r in results] == [None, None, None] + assert _process(tmp_path, "DiesInAWorker")["pid"] == os.getpid() + + +def test_a_run_of_few_experiments_is_serial_by_default(tmp_path: Path) -> None: + _run(tmp_path / "few", None, CLASSES[:EXPERIMENTS_PER_WORKER]) + _run(tmp_path / "more", None, CLASSES) + + assert _process(tmp_path / "few", "First")["pid"] == os.getpid() + assert _process(tmp_path / "more", "First")["pid"] != os.getpid() + + +def test_experiments_which_cannot_be_pickled_run_serially(tmp_path: Path) -> None: + """A class of a function cannot reach a worker, the run is serial.""" + + class Local(First): + pass + + results = _run(tmp_path, 2, [Local, Second]) + + assert [r.error for r in results] == [None, None] + assert _process(tmp_path, "Local")["pid"] == os.getpid() + assert _process(tmp_path, "Second")["pid"] == os.getpid() + + +def test_a_number_of_workers_below_one_raises(tmp_path: Path) -> None: + with pytest.raises(ValueError, match="at least 1"): + _run(tmp_path, 0, [First]) + + +def test_a_simulator_with_a_pool_runs_the_experiments_serially_by_default( + tmp_path: Path, +) -> None: + """The scans of an experiment run serially in a worker. + + A simulator with a pool of its own runs every scan on all its workers, a + pool of experiments would run them serially; an explicit number of + workers runs the experiments in a pool all the same. + """ + _run(tmp_path / "default", None, CLASSES, Simulator(n_workers=4)) + _run(tmp_path / "explicit", 2, CLASSES, Simulator(n_workers=4)) + + assert _process(tmp_path / "default", "First")["pid"] == os.getpid() + assert _process(tmp_path / "explicit", "First")["pid"] != os.getpid() + + +def test_an_experiment_of_the_main_script_runs_serially( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch, caplog: pytest.LogCaptureFixture +) -> None: + """A worker imports the script without its guard, a class in it may be missing. + + A class defined in the guard pickles in the script and broke the pool. + """ + main = type("InTheScript", (First,), {"__module__": "__main__"}) + monkeypatch.setattr(sys.modules["__main__"], "InTheScript", main, raising=False) + + with caplog.at_level(logging.INFO): + results = _run(tmp_path, 2, [main, Second]) + + assert [r.error for r in results] == [None, None] + assert _process(tmp_path, "InTheScript")["pid"] == os.getpid() + assert "died" not in caplog.text + + +def test_the_results_of_a_pooled_experiment_are_released(tmp_path: Path) -> None: + (result, _) = _run(tmp_path, 2, [First, Second]) + + with pytest.raises(RuntimeError, match="released"): + _ = result.experiment.results + + +def test_a_worker_logs_its_messages_only_through_the_run(tmp_path: Path) -> None: + """A handler of the worker, e.g. of a script which configures logging, is + silent while an experiment runs: the run logs its messages once.""" + seen: list[logging.LogRecord] = [] + + class Seen(logging.Handler): + def emit(self, record: logging.LogRecord) -> None: + seen.append(record) + + runner = _runner(BadFigure) + job = runner._jobs( + tmp_path, + ["svg"], + None, + { + "show_figures": False, + "save_results": False, + "reduced_selections": True, + "keep_results": False, + "figure_settings": None, + }, + )["BadFigure"] + handler = Seen() + logging.root.addHandler(handler) + try: + outcome = _run_job(job) + finally: + logging.root.removeHandler(handler) + + assert seen == [] + assert any("fig_bad" in r.getMessage() for r in outcome.records) + + +def test_a_message_whose_arguments_cannot_be_pickled_is_sent_as_text() -> None: + collect = _Collect() + record = logging.LogRecord("x", logging.INFO, __file__, 1, "%s", (object(),), None) + record.unpicklable = lambda: None + collect.emit(record) + + (sent,) = collect.records + assert pickle.loads(pickle.dumps(sent)).getMessage().startswith("