diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index def5998..594c32f 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -31,9 +31,25 @@ This is a constraint on the math, not the encoding: `v2/` files may (and often s structurally different from a mechanical v1-to-v2 port to make idiomatic use of v2 features (experiments table, v2 prior/distribution syntax, mapping tables, ...) where they fit better. They just have to evaluate to the same log-likelihood / log-posterior / gradient at nominal -parameters as `v1/`. +parameters as `v1/`, up to the reparametrization described below. -For problems with priors (maximum a posteriori estimation), the policy is not decided yet: -v1 priors on log-scaled parameters cannot, in general, be expressed as v2 priors without -changing either the optimum or the prior distribution. See -[#333](https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/issues/333). +PEtab v2 has no `parameterScale`: priors and the objective are always expressed in terms of the +parameter itself. To preserve the v1 optimization problem (and, for maximum a posteriori problems, +the v1 prior distribution), estimated `log`/`log10`-scaled v1 parameters are reparametrized in +`v2/`: an estimated `log10_` (or `log_`) replaces ``, and `10^log10_` +(or `exp(log_)`) is substituted everywhere `` appeared (model, observable and noise +formulas, condition and measurement table overrides). A v1 `parameterScale*` prior becomes the +corresponding prior on the new parameter; the v1 default prior (uniform on parameter scale) +becomes the v2 default prior (uniform within the transformed bounds). +Parameters with a prior on linear scale (e.g. a v1 `objectivePriorType` of `normal`, `laplace`, +...) are not reparametrized; their prior is already the same in both encodings. + +Measurement table overrides cannot be expressions. Therefore, for a reparametrized parameter +used as an override, the placeholder is transformed instead (`` becomes +`10^` in the observable or noise formula), and all overrides of that placeholder are +log-transformed. If the same placeholder also receives non-positive numbers or parameters that +stay on linear scale, the observable is split, so that each copy's placeholder only receives +overrides that can be log-transformed, or none at all. + +See [#333](https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/issues/333) for +the rationale (option C there). diff --git a/problems/Boehm_JProteomeRes2014/v2/measurements.tsv b/problems/Boehm_JProteomeRes2014/v2/measurements.tsv new file mode 100644 index 0000000..376c33b --- /dev/null +++ b/problems/Boehm_JProteomeRes2014/v2/measurements.tsv @@ -0,0 +1,49 @@ +observableId experimentId time measurement noiseParameters datasetId +pSTAT5A_rel 0.0 7.90107299873911 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 2.5 66.3634939682528 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 5.0 81.1713239165337 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 10.0 94.7303080627876 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 15.0 95.1164830493244 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 20.0 91.44171655392 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 30.0 91.2570992263545 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 40.0 93.672297840992 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 50.0 88.754232823019 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 60.0 85.2697032178279 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 80.0 81.1323953369608 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 100.0 76.135928475796 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 120.0 65.2480591250229 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 160.0 42.5996587097166 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 200.0 25.1577975445976 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5A_rel 240.0 15.4301823994994 log10_sd_pSTAT5A_rel model1_data1_pSTAT5A_rel +pSTAT5B_rel 0.0 4.59653334288927 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 2.5 29.6345459874189 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 5.0 46.0438064695701 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 10.0 81.9747336191685 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 15.0 80.5716092952489 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 20.0 79.0357196371442 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 30.0 75.6723803658298 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 40.0 71.6247198618303 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 50.0 69.0628632755742 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 60.0 67.1473843168934 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 80.0 60.8994762938824 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 100.0 54.8092577721492 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 120.0 43.9812899778025 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 160.0 29.7714581636675 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 200.0 20.089016562771 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +pSTAT5B_rel 240.0 10.961845166698 log10_sd_pSTAT5B_rel model1_data1_pSTAT5B_rel +rSTAT5A_rel 0.0 14.7231682180584 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 2.5 33.762342293107 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 5.0 36.7998512942493 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 10.0 49.7176022859654 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 15.0 46.9281200952955 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 20.0 47.8365745550135 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 30.0 46.9287272497385 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 40.0 40.5977529351181 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 50.0 43.783663885071 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 60.0 44.4573876523855 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 80.0 41.3271592618297 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 100.0 41.0627332148332 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 120.0 39.2358300289377 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 160.0 36.6194605442418 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 200.0 34.8937144010561 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel +rSTAT5A_rel 240.0 32.2110771608676 log10_sd_rSTAT5A_rel model1_data1_rSTAT5A_rel diff --git a/problems/Boehm_JProteomeRes2014/v2/model.xml b/problems/Boehm_JProteomeRes2014/v2/model.xml new file mode 100755 index 0000000..f886f4f --- /dev/null +++ b/problems/Boehm_JProteomeRes2014/v2/model.xml @@ -0,0 +1,350 @@ + + + + + + 
PEtab implementation of the model from Boehm et al. (2014), J. Proteome Res., 13, 12, 5685-5694
+ + + + + + + + + + Boehm + Martin + + + + + + 2019-11-29T09:41:48Z + + + 2019-11-29T09:41:48Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 207.6 + ratio + + + + + + + + 207.6 + + + 207.6 + ratio + + + + + + + + + + + 1.25 -7 + + + + + -1 + Epo_degradation_BaF3 + t + + + + + + + + + + 10 + log10_Epo_degradation_BaF3 + + + + + + + + 10 + log10_k_exp_hetero + + + + + + + + 10 + log10_k_exp_homo + + + + + + + + 10 + log10_k_imp_hetero + + + + + + + + 10 + log10_k_imp_homo + + + + + + + + 10 + log10_k_phos + + + + + + + + + + + + + + + + + cyt + BaF3_Epo + + + STAT5A + 2 + + k_phos + + + + + + + + + + + + + + + + + cyt + BaF3_Epo + STAT5A + STAT5B + k_phos + + + + + + + + + + + + + + + + cyt + BaF3_Epo + + + STAT5B + 2 + + k_phos + + + + + + + + + + + + + + + + cyt + k_imp_homo + pApA + + + + + + + + + + + + + + + + cyt + k_imp_hetero + pApB + + + + + + + + + + + + + + + + cyt + k_imp_homo + pBpB + + + + + + + + + + + + + + + + nuc + k_exp_homo + nucpApA + + + + + + + + + + + + + + + + + nuc + k_exp_hetero + nucpApB + + + + + + + + + + + + + + + + nuc + k_exp_homo + nucpBpB + + + + + + + diff --git a/problems/Boehm_JProteomeRes2014/v2/observables.tsv b/problems/Boehm_JProteomeRes2014/v2/observables.tsv new file mode 100644 index 0000000..3e74b60 --- /dev/null +++ b/problems/Boehm_JProteomeRes2014/v2/observables.tsv @@ -0,0 +1,4 @@ +observableId observableFormula noiseFormula noiseDistribution noisePlaceholders +pSTAT5A_rel (200*pApA*specC17 + 100*pApB)/(STAT5A*specC17 + 2*pApA*specC17 + pApB) 10 ^ noiseParameter1_pSTAT5A_rel normal noiseParameter1_pSTAT5A_rel +pSTAT5B_rel (-100*pApB + 200*pBpB*(specC17 - 1))/(STAT5B*(specC17 - 1) - pApB + 2*pBpB*(specC17 - 1)) 10 ^ noiseParameter1_pSTAT5B_rel normal noiseParameter1_pSTAT5B_rel +rSTAT5A_rel (100*STAT5A*specC17 + 200*pApA*specC17 + 100*pApB)/(STAT5A*specC17 - STAT5B*(specC17 - 1) + 2*pApA*specC17 + 2*pApB - 2*pBpB*(specC17 - 1)) 10 ^ noiseParameter1_rSTAT5A_rel normal noiseParameter1_rSTAT5A_rel diff --git a/problems/Boehm_JProteomeRes2014/v2/parameters.tsv b/problems/Boehm_JProteomeRes2014/v2/parameters.tsv new file mode 100644 index 0000000..2c151a6 --- /dev/null +++ b/problems/Boehm_JProteomeRes2014/v2/parameters.tsv @@ -0,0 +1,12 @@ +parameterId parameterName lowerBound upperBound nominalValue estimate +log10_Epo_degradation_BaF3 $\log_{10}(\mathrm{EPO}_\mathrm{degradation,BaF3})$ -5.0 5.0 -1.5689175883999986 true +log10_k_exp_hetero $\log_{10}(k_\mathrm{exp,hetero})$ -5.0 5.0 -4.999704893599997 true +log10_k_exp_homo $\log_{10}(k_\mathrm{exp,homo})$ -5.0 5.0 -2.209698781700016 true +log10_k_imp_hetero $\log_{10}(k_\mathrm{imp,hetero})$ -5.0 5.0 -1.78600654750001 true +log10_k_imp_homo $\log_{10}(k_\mathrm{imp,homo})$ -5.0 5.0 4.990114008799999 true +log10_k_phos $\log_{10}(k_\mathrm{phos})$ -5.0 5.0 4.1977354885 true +ratio ratio 0.0 5.0 0.693 false +log10_sd_pSTAT5A_rel $\log_{10}(\sigma_\mathrm{pSTAT5A,rel})$ -5.0 5.0 0.5857552705999998 true +log10_sd_pSTAT5B_rel $\log_{10}(\sigma_\mathrm{pSTAT5B,rel})$ -5.0 5.0 0.8189828191999998 true +log10_sd_rSTAT5A_rel $\log_{10}(\sigma_\mathrm{rSTAT5A,rel})$ -5.0 5.0 0.4986844040000004 true +specC17 specC17 -5.0 5.0 0.107 false diff --git a/problems/Boehm_JProteomeRes2014/v2/problem.yaml b/problems/Boehm_JProteomeRes2014/v2/problem.yaml new file mode 100644 index 0000000..81144d8 --- /dev/null +++ b/problems/Boehm_JProteomeRes2014/v2/problem.yaml @@ -0,0 +1,11 @@ +format_version: 2.0.0 +parameter_files: +- parameters.tsv +model_files: + model: + location: model.xml + language: sbml +measurement_files: +- measurements.tsv +observable_files: +- observables.tsv diff --git a/src/python/benchmark_models_petab/check_petablint.py b/src/python/benchmark_models_petab/check_petablint.py index 655a478..d081238 100644 --- a/src/python/benchmark_models_petab/check_petablint.py +++ b/src/python/benchmark_models_petab/check_petablint.py @@ -2,25 +2,44 @@ import sys from pathlib import Path -from . import MODELS, get_problem_yaml_path +from . import MODELS, v1, v2 +from .C import MODELS_DIR def main(): - """Check all PEtab problems with petablint.""" - num_failures = 0 + """Check all PEtab problems with petablint. + Checks all encodings (`v1/`, `v2/`) of each problem. A problem without + any encoding counts as a failure. + """ + num_checked = 0 + num_failures = 0 for petab_problem_id in MODELS: - print(petab_problem_id, flush=True) - petab_yaml = get_problem_yaml_path(petab_problem_id) - ret = os.system(f"petablint -v {petab_yaml}") - print("=" * 100, flush=True) # just for output readability - - if ret: + petab_yamls = [] + for version in (v1, v2): + try: + petab_yamls.append( + version.get_problem_yaml_path(petab_problem_id) + ) + except ValueError: + pass + if not petab_yamls: + print(f"{petab_problem_id}: no PEtab problem found", flush=True) + num_checked += 1 num_failures += 1 - num_passed = len(MODELS) - num_failures + for petab_yaml in petab_yamls: + print(petab_yaml.relative_to(MODELS_DIR), flush=True) + ret = os.system(f"petablint -v {petab_yaml}") + print("=" * 100, flush=True) # just for output readability + + num_checked += 1 + if ret: + num_failures += 1 + + num_passed = num_checked - num_failures print(f"Result: {Path(__file__).stem}") - print(f"{num_passed} out of {len(MODELS)} passed.") - print(f"{num_failures} out of {len(MODELS)} failed.") - # Fail unless all models passed + print(f"{num_passed} out of {num_checked} passed.") + print(f"{num_failures} out of {num_checked} failed.") + # Fail unless all problems passed sys.exit(num_failures) diff --git a/src/python/test/test_v2.py b/src/python/test/test_v2.py index ced6147..1a2d6e8 100644 --- a/src/python/test/test_v2.py +++ b/src/python/test/test_v2.py @@ -1,25 +1,33 @@ -"""Tests for the v2 accessors. - -No problem in the collection has a `v2/` encoding yet, so only the -"not found" error path can be tested here. Once real `problems//v2/` -directories exist, add a happy-path test (`get_problem`/`get_simulation_df` -on an actual v2 problem) alongside these. -""" +"""Tests for the v2 accessors.""" import pytest import benchmark_models_petab as models +NO_V2_PROBLEM_ID = "NotAProblem_2000" + def test_get_problem_yaml_path_raises_for_missing_v2(): with pytest.raises(ValueError, match="v2"): - models.v2.get_problem_yaml_path(models.MODELS[0]) + models.v2.get_problem_yaml_path(NO_V2_PROBLEM_ID) def test_get_problem_raises_for_missing_v2(): with pytest.raises(ValueError, match="v2"): - models.v2.get_problem(models.MODELS[0]) + models.v2.get_problem(NO_V2_PROBLEM_ID) def test_get_simulation_df_returns_none_for_missing_v2(): - assert models.v2.get_simulation_df(models.MODELS[0]) is None + assert models.v2.get_simulation_df(NO_V2_PROBLEM_ID) is None + + +def test_get_problem_yaml_path(): + path = models.v2.get_problem_yaml_path("Boehm_JProteomeRes2014") + assert path.is_file() + + +def test_get_problem(): + problem = models.v2.get_problem("Boehm_JProteomeRes2014") + assert "log10_k_phos" in problem.parameter_df.index + assert "k_phos" not in problem.parameter_df.index + assert len(problem.measurement_df) > 0