PEtab-dev · LeonardSchmiester · Nov 26, 2020 · Nov 25, 2020 · Nov 25, 2020 · Nov 25, 2020
diff --git a/petabtests/cases/0010/0010.py b/petabtests/cases/0010/0010.py
@@ -14,7 +14,7 @@
 condition_df = pd.DataFrame(data={
     CONDITION_ID: ['preeq_c0', 'c0'],
     'k1': [0.3, 0.8],
-    'B': [0, 0],
+    'B': [0, 1],
 }).set_index([CONDITION_ID])
 
 measurement_df = pd.DataFrame(data={
@@ -49,7 +49,7 @@
 steady_state_a = analytical_a(1000, 1, 0, 0.3, 0.6)
 # use steady state as initial state
 simulation_df[SIMULATION] = [
-    analytical_a(t, steady_state_a, 0, 0.8, 0.6)
+    analytical_a(t, steady_state_a, 1, 0.8, 0.6)
     for t in simulation_df[TIME]]
 
 chi2 = petab.calculate_chi2(

diff --git a/petabtests/cases/0010/_0010_solution.yaml b/petabtests/cases/0010/_0010_solution.yaml
@@ -1,5 +1,5 @@
-chi2: 0.548440802855135
-llh: -0.7258031067170223
+chi2: 1.5094134279439457
+llh: -1.2062894192614277
 simulation_files:
 - _simulations.tsv
 tol_chi2: 0.001

diff --git a/petabtests/cases/0010/_conditions.tsv b/petabtests/cases/0010/_conditions.tsv
@@ -1,3 +1,3 @@
 conditionId	k1	B
 preeq_c0	0.3	0
-c0	0.8	0
+c0	0.8	1
diff --git a/petabtests/cases/0010/_simulations.tsv b/petabtests/cases/0010/_simulations.tsv
@@ -1,3 +1,3 @@
 observableId	preequilibrationConditionId	simulationConditionId	time	simulation
-obs_a	preeq_c0	c0	1	0.3796559862634692
-obs_a	preeq_c0	c0	10	0.28571460248713115
+obs_a	preeq_c0	c0	1	0.7025430017170664
+obs_a	preeq_c0	c0	10	0.7142856746891086
diff --git a/petabtests/cases/0017/0017.md b/petabtests/cases/0017/0017.md
@@ -0,0 +1,19 @@
+# PEtab test case 0017
+
+## Objective 
+
+This case tests support for partial preequilibration with `NaN`'s in the 
+condition file.
+
+The model is to be simulated for a preequilibration condition and a
+simulation condition.
+For preequilibration, species `B` is initialized with `0`. For simulation,
+`B` is set to `NaN`, meaning that it is initialized with the result from
+preequilibration.
+`A` is reinitialized to the value in the condition table after 
+preequilibration.
+
+## Model
+
+A simple conversion reaction `A <=> B` in a single compartment, following
+mass action kinetics.
diff --git a/petabtests/cases/0017/0017.py b/petabtests/cases/0017/0017.py
@@ -0,0 +1,61 @@
+from petabtests import *
+from petab.C import *
+import petab
+
+import pandas as pd
+
+
+test_id = 17
+
+# problem --------------------------------------------------------------------
+
+model = DEFAULT_SBML_FILE
+
+condition_df = pd.DataFrame(data={
+    CONDITION_ID: ['preeq_c0', 'c0'],
+    'k1': [0.3, 0.8],
+    'B': [0, 'NaN'],
+    'A': [0, 1],
+}).set_index([CONDITION_ID])
+
+measurement_df = pd.DataFrame(data={
+    OBSERVABLE_ID: ['obs_a', 'obs_a'],
+    PREEQUILIBRATION_CONDITION_ID: ['preeq_c0', 'preeq_c0'],
+    SIMULATION_CONDITION_ID: ['c0', 'c0'],
+    TIME: [1, 10],
+    MEASUREMENT: [0.7, 0.1]
+})
+
+observable_df = pd.DataFrame(data={
+    OBSERVABLE_ID: ['obs_a'],
+    OBSERVABLE_FORMULA: ['A'],
+    NOISE_FORMULA: [0.5]
+}).set_index([OBSERVABLE_ID])
+
+parameter_df = pd.DataFrame(data={
+    PARAMETER_ID: ['k2'],
+    PARAMETER_SCALE: [LIN],
+    LOWER_BOUND: [0],
+    UPPER_BOUND: [10],
+    NOMINAL_VALUE: [0.6],
+    ESTIMATE: [1],
+}).set_index(PARAMETER_ID)
+
+
+# solutions ------------------------------------------------------------------
+
+simulation_df = measurement_df.copy(deep=True).rename(
+    columns={MEASUREMENT: SIMULATION})
+# simulate for far time point as steady state
+steady_state_b = analytical_b(1000, 0, 0, 0.3, 0.6)
+# use steady state as initial state
+simulation_df[SIMULATION] = [
+    analytical_a(t, 1, steady_state_b, 0.8, 0.6)
+    for t in simulation_df[TIME]]
+
+chi2 = petab.calculate_chi2(
+    measurement_df, simulation_df, observable_df, parameter_df)
+
+llh = petab.calculate_llh(
+    measurement_df, simulation_df, observable_df, parameter_df)
+print(llh)
diff --git a/petabtests/cases/0017/_0017.yaml b/petabtests/cases/0017/_0017.yaml
@@ -0,0 +1,11 @@
+format_version: 1
+parameter_file: _parameters.tsv
+problems:
+- condition_files:
+  - _conditions.tsv
+  measurement_files:
+  - _measurements.tsv
+  observable_files:
+  - _observables.tsv
+  sbml_files:
+  - _model.xml
diff --git a/petabtests/cases/0017/_0017_solution.yaml b/petabtests/cases/0017/_0017_solution.yaml
@@ -0,0 +1,7 @@
+chi2: 0.4999757102229031
+llh: -0.7015705604009064
+simulation_files:
+- _simulations.tsv
+tol_chi2: 0.001
+tol_llh: 0.001
+tol_simulations: 0.001
diff --git a/petabtests/cases/0017/_conditions.tsv b/petabtests/cases/0017/_conditions.tsv
@@ -0,0 +1,3 @@
+conditionId	k1	B	A
+preeq_c0	0.3	0	0
+c0	0.8	NaN	1
diff --git a/petabtests/cases/0017/_measurements.tsv b/petabtests/cases/0017/_measurements.tsv
@@ -0,0 +1,3 @@
+observableId	preequilibrationConditionId	simulationConditionId	time	measurement
+obs_a	preeq_c0	c0	1	0.7
+obs_a	preeq_c0	c0	10	0.1
diff --git a/petabtests/cases/0017/_model.xml b/petabtests/cases/0017/_model.xml
@@ -0,0 +1,94 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<sbml xmlns="http://www.sbml.org/sbml/level2/version4" level="2" version="4">
+  <model id="conversion_reaction_0" name="Conversion Reaction 0">
+
+    <listOfUnitDefinitions>
+      <unitDefinition id="volume" name="volume">
+        <listOfUnits>
+          <unit kind="litre" exponent="1" scale="-3" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+      <unitDefinition id="substance" name="substance">
+        <listOfUnits>
+          <unit kind="mole" exponent="1" scale="-3" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+    </listOfUnitDefinitions>
+
+    <listOfCompartments>
+      <compartment id="compartment" name="compartment" spatialDimensions="3" size="1" constant="true">
+      </compartment>
+    </listOfCompartments>
+
+    <listOfSpecies>
+      <species id="A" name="A" compartment="compartment" initialConcentration="1" boundaryCondition="false" constant="false">
+      </species>
+      <species id="B" name="B" compartment="compartment" initialConcentration="1" boundaryCondition="false" constant="false">
+      </species>
+    </listOfSpecies>
+
+    <listOfParameters>
+      <parameter id="a0" name="a0" value="1" constant="true">
+      </parameter>
+      <parameter id="b0" name="b0" value="1" constant="true">
+      </parameter>
+      <parameter id="k1" name="k1" value="0" constant="true">
+      </parameter>
+      <parameter id="k2" name="k2" value="0" constant="true">
+      </parameter>
+    </listOfParameters>
+
+    <listOfInitialAssignments>
+      <initialAssignment symbol="A">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <ci> a0 </ci>
+        </math>
+      </initialAssignment>
+      <initialAssignment symbol="B">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <ci> b0 </ci>
+        </math>
+      </initialAssignment>
+    </listOfInitialAssignments>
+
+    <listOfReactions>
+      <reaction id="fwd" name="fwd" reversible="false">
+        <listOfReactants>
+          <speciesReference species="A" stoichiometry="1"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference species="B" stoichiometry="1"/>
+        </listOfProducts>
+        <kineticLaw>
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <times/>
+              <ci> compartment </ci>
+              <ci> k1 </ci>
+              <ci> A </ci>
+            </apply>
+          </math>
+        </kineticLaw>
+      </reaction>
+      <reaction id="rev" name="rev" reversible="false">
+        <listOfReactants>
+          <speciesReference species="B" stoichiometry="1"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference species="A" stoichiometry="1"/>
+        </listOfProducts>
+        <kineticLaw>
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <times/>
+              <ci> compartment </ci>
+              <ci> k2 </ci>
+              <ci> B </ci>
+            </apply>
+          </math>
+        </kineticLaw>
+      </reaction>
+    </listOfReactions>
+
+  </model>
+</sbml>
diff --git a/petabtests/cases/0017/_observables.tsv b/petabtests/cases/0017/_observables.tsv
@@ -0,0 +1,2 @@
+observableId	observableFormula	noiseFormula
+obs_a	A	0.5
diff --git a/petabtests/cases/0017/_parameters.tsv b/petabtests/cases/0017/_parameters.tsv
@@ -0,0 +1,2 @@
+parameterId	parameterScale	lowerBound	upperBound	nominalValue	estimate
+k2	lin	0	10	0.6	1
diff --git a/petabtests/cases/0017/_simulations.tsv b/petabtests/cases/0017/_simulations.tsv
@@ -0,0 +1,3 @@
+observableId	preequilibrationConditionId	simulationConditionId	time	simulation
+obs_a	preeq_c0	c0	1	0.5694839793952038
+obs_a	preeq_c0	c0	10	0.42857190373069665
diff --git a/petabtests/cases/pysb/0010/_0010_solution.yaml b/petabtests/cases/pysb/0010/_0010_solution.yaml
@@ -1,5 +1,5 @@
-chi2: 0.548440802855135
-llh: -0.7258031067170223
+chi2: 1.5094134279439457
+llh: -1.2062894192614277
 simulation_files:
 - _simulations.tsv
 tol_chi2: 0.001

diff --git a/petabtests/cases/pysb/0010/_conditions.tsv b/petabtests/cases/pysb/0010/_conditions.tsv
@@ -1,3 +1,3 @@
 conditionId	k1	B
 preeq_c0	0.3	0
-c0	0.8	0
+c0	0.8	1
diff --git a/petabtests/cases/pysb/0010/_simulations.tsv b/petabtests/cases/pysb/0010/_simulations.tsv
@@ -1,3 +1,3 @@
 observableId	preequilibrationConditionId	simulationConditionId	time	simulation
-obs_a	preeq_c0	c0	1	0.3796559862634692
-obs_a	preeq_c0	c0	10	0.28571460248713115
+obs_a	preeq_c0	c0	1	0.7025430017170664
+obs_a	preeq_c0	c0	10	0.7142856746891086
diff --git a/petabtests/cases/pysb/0017/_0017.yaml b/petabtests/cases/pysb/0017/_0017.yaml
@@ -0,0 +1,11 @@
+format_version: 1
+parameter_file: _parameters.tsv
+problems:
+- condition_files:
+  - _conditions.tsv
+  measurement_files:
+  - _measurements.tsv
+  observable_files:
+  - _observables.tsv
+  sbml_files:
+  - _model.py
diff --git a/petabtests/cases/pysb/0017/_0017_solution.yaml b/petabtests/cases/pysb/0017/_0017_solution.yaml
@@ -0,0 +1,7 @@
+chi2: 0.4999757102229031
+llh: -0.7015705604009064
+simulation_files:
+- _simulations.tsv
+tol_chi2: 0.001
+tol_llh: 0.001
+tol_simulations: 0.001
diff --git a/petabtests/cases/pysb/0017/_conditions.tsv b/petabtests/cases/pysb/0017/_conditions.tsv
@@ -0,0 +1,3 @@
+conditionId	k1	B	A
+preeq_c0	0.3	0	0
+c0	0.8	NaN	1
diff --git a/petabtests/cases/pysb/0017/_measurements.tsv b/petabtests/cases/pysb/0017/_measurements.tsv
@@ -0,0 +1,3 @@
+observableId	preequilibrationConditionId	simulationConditionId	time	measurement
+obs_a	preeq_c0	c0	1	0.7
+obs_a	preeq_c0	c0	10	0.1
diff --git a/petabtests/cases/pysb/0017/_model.py b/petabtests/cases/pysb/0017/_model.py
@@ -0,0 +1,21 @@
+from pysb import Model, Monomer, Parameter, Compartment, Rule, Initial, Observable
+
+Model()
+
+Compartment('compartment')
+
+Monomer('A_')
+Monomer('B_')
+
+Parameter('a0', 1.0)
+Parameter('b0', 1.0)
+Parameter('k1', 0.0)
+Parameter('k2', 0.0)
+
+Rule('conversion', A_() ** compartment | B_() ** compartment, k1, k2)
+
+Initial(A_() ** compartment, a0)
+Initial(B_() ** compartment, b0)
+
+Observable("A", A_())
+Observable("B", B_())
diff --git a/petabtests/cases/pysb/0017/_observables.tsv b/petabtests/cases/pysb/0017/_observables.tsv
@@ -0,0 +1,2 @@
+observableId	observableFormula	noiseFormula
+obs_a	A	0.5
diff --git a/petabtests/cases/pysb/0017/_parameters.tsv b/petabtests/cases/pysb/0017/_parameters.tsv
@@ -0,0 +1,2 @@
+parameterId	parameterScale	lowerBound	upperBound	nominalValue	estimate
+k2	lin	0	10	0.6	1
diff --git a/petabtests/cases/pysb/0017/_simulations.tsv b/petabtests/cases/pysb/0017/_simulations.tsv
@@ -0,0 +1,3 @@
+observableId	preequilibrationConditionId	simulationConditionId	time	simulation
+obs_a	preeq_c0	c0	1	0.5694839793952038
+obs_a	preeq_c0	c0	10	0.42857190373069665