Added additional stochastic models[finished]

CHANGELOG.md

@@ -5,6 +5,7 @@ All notable changes to this project will be documented in this file.
 ## Unreleased
 
 ### Added
+- [#1432](https://github.com/pints-team/pints/pull/1432) Added 2 new stochastic models: production and degradation model, Schlogl's system of chemical reactions. Moved the stochastic logistic model into `pints.stochastic` to take advantage of the `MarkovJumpModel`.
 - [#1420](https://github.com/pints-team/pints/pull/1420) The `Optimiser` class now distinguishes between a best-visited point (`x_best`, with score `f_best`) and a best-guessed point (`x_guessed`, with approximate score `f_guessed`). For most optimisers, the two values are equivalent. The `OptimisationController` still tracks `x_best` and `f_best` by default, but this can be modified using the methods `set_f_guessed_tracking` and `f_guessed_tracking`.
 - [#1417](https://github.com/pints-team/pints/pull/1417) Added a module `toy.stochastic` for stochastic models. In particular, `toy.stochastic.MarkovJumpModel` implements Gillespie's algorithm for easier future implementation of stochastic models.
 

docs/source/toy/index.rst

@@ -37,5 +37,4 @@ Some toy classes provide extra functionality defined in the
     simple_egg_box_logpdf
     simple_harmonic_oscillator_model
     sir_model
-    stochastic_logistic_model
     twisted_gaussian_logpdf

docs/source/toy/stochastic/index.rst

@@ -12,4 +12,7 @@ examples.
 
     markov_jump_model
     stochastic_degradation_model
+    stochastic_logistic_model
     stochastic_michaelis_menten_model
+    stochastic_production_degradation_model
+    stochastic_schlogl_model

docs/source/toy/stochastic_logistic_model.rst → docs/source/toy/stochastic/stochastic_logistic_model.rst

@@ -2,6 +2,6 @@
 Stochastic Logistic Model
 *************************
 
-.. currentmodule:: pints.toy
+.. currentmodule:: pints.toy.stochastic
 
-.. autoclass:: StochasticLogisticModel
+.. autoclass:: LogisticModel

docs/source/toy/stochastic/stochastic_production_degradation_model.rst

@@ -0,0 +1,7 @@
+*******************************************
+Stochastic production and degradation model
+*******************************************
+
+.. currentmodule:: pints.toy.stochastic
+
+.. autoclass:: ProductionDegradationModel

docs/source/toy/stochastic/stochastic_schlogl_model.rst

@@ -0,0 +1,7 @@
+***************
+Schlogl's model
+***************
+
+.. currentmodule:: pints.toy.stochastic
+
+.. autoclass:: SchloglModel

examples/README.md

@@ -120,6 +120,8 @@ relevant code.
 - [Stochastic Degradation model](./toy/model-stochastic-degradation.ipynb)
 - [Stochastic Logistic model](./toy/model-stochastic-logistic-growth.ipynb)
 - [Stochastic Michaelis Menten model](./toy/model-stochastic-michaelis-menten.ipynb)
+- [Stochastic Production Degradation model](toy/model-stochastic-production-degradation.ipynb)
+- [Stochastic Schlogl model](toy/model-stochastic-schlogl.ipynb)
 
 ### Distributions
 - [Annulus](./toy/distribution-annulus.ipynb)

pints/tests/test_toy_markov_jump_model.py

@@ -61,16 +61,13 @@ def test_simulate(self):
 
     def test_errors(self):
         model = DegradationModel(20)
-        # parameters, times cannot be negative
-        times = np.linspace(0, 100, 101)
-        parameters = [-0.1]
-        self.assertRaises(ValueError, model.simulate, parameters, times)
-
+        # times cannot be negative
         times_2 = np.linspace(-10, 10, 21)
         parameters_2 = [0.1]
         self.assertRaises(ValueError, model.simulate, parameters_2, times_2)
 
         # this model should have 1 parameter
+        times = np.linspace(0, 100, 101)
         parameters_3 = [0.1, 1]
         self.assertRaises(ValueError, model.simulate, parameters_3, times)
 

pints/tests/test_toy_stochastic_logistic_model.py

@@ -10,6 +10,7 @@
 import numpy as np
 import pints
 import pints.toy
+import pints.toy.stochastic
 
 
 class TestStochasticLogisticModel(unittest.TestCase):
@@ -22,7 +23,7 @@ def test_start_with_zero(self):
         # Set seed for random generator
         np.random.seed(1)
 
-        model = pints.toy.StochasticLogisticModel(0)
+        model = pints.toy.stochastic.LogisticModel(0)
         times = [0, 1, 2, 100, 1000]
         parameters = [0.1, 50]
         values = model.simulate(parameters, times)
@@ -35,7 +36,7 @@ def test_start_with_one(self):
         # Set seed for random generator
         np.random.seed(1)
 
-        model = pints.toy.StochasticLogisticModel(1)
+        model = pints.toy.stochastic.LogisticModel(1)
         times = [0, 1, 2, 100, 1000]
         parameters = [0.1, 50]
         values = model.simulate(parameters, times)
@@ -46,7 +47,7 @@ def test_start_with_one(self):
 
     def test_suggested(self):
         # Check suggested values
-        model = pints.toy.StochasticLogisticModel(1)
+        model = pints.toy.stochastic.LogisticModel(1)
         times = model.suggested_times()
         parameters = model.suggested_parameters()
         self.assertTrue(len(times) == 101)
@@ -55,34 +56,32 @@ def test_suggested(self):
     def test_simulate(self):
         # Check each step in the simulation process
         np.random.seed(1)
-        model = pints.toy.StochasticLogisticModel(1)
+        model = pints.toy.stochastic.LogisticModel(1)
         times = np.linspace(0, 100, 101)
-        params = [0.1, 50]
-        time, raw_values = model._simulate_raw([0.1, 50])
-        values = model._interpolate_values(time, raw_values, times, params)
+        time, raw_values = model.simulate_raw([0.1, 50], 100)
+        values = model.interpolate_mol_counts(time, raw_values, times)
         self.assertTrue(len(time), len(raw_values))
 
         # Test output of Gillespie algorithm
-        self.assertTrue(np.all(raw_values == np.array(range(1, 51))))
+        raw_values = np.concatenate(raw_values)
+        self.assertTrue(np.all(raw_values == np.array(range(1, 28))))
 
         # Check simulate function returns expected values
         self.assertTrue(np.all(values[np.where(times < time[1])] == 1))
 
         # Check interpolation function works as expected
         temp_time = np.array([np.random.uniform(time[0], time[1])])
-        self.assertTrue(model._interpolate_values(time, raw_values, temp_time,
-                                                  params)[0] == 1)
+        self.assertTrue(model.interpolate_mol_counts(time, raw_values,
+                                                     temp_time)[0] == 1)
         temp_time = np.array([np.random.uniform(time[1], time[2])])
-        self.assertTrue(model._interpolate_values(time, raw_values, temp_time,
-                                                  params)[0] == 2)
+        self.assertTrue(model.interpolate_mol_counts(time, raw_values,
+                                                     temp_time)[0] == 2)
 
         # Check parameters, times cannot be negative
         parameters_0 = [-0.1, 50]
-        self.assertRaises(ValueError, model.simulate, parameters_0, times)
         self.assertRaises(ValueError, model.mean, parameters_0, times)
 
         parameters_1 = [0.1, -50]
-        self.assertRaises(ValueError, model.simulate, parameters_1, times)
         self.assertRaises(ValueError, model.mean, parameters_1, times)
 
         times_2 = np.linspace(-10, 10, 21)
@@ -96,21 +95,15 @@ def test_simulate(self):
         self.assertRaises(ValueError, model.mean, parameters_3, times)
 
         # Check initial value cannot be negative
-        self.assertRaises(ValueError, pints.toy.StochasticLogisticModel, -1)
+        self.assertRaises(ValueError, pints.toy.stochastic.LogisticModel, -1)
 
-    def test_mean_variance(self):
+    def test_mean(self):
         # Check the mean is what we expected
-        model = pints.toy.StochasticLogisticModel(1)
+        model = pints.toy.stochastic.LogisticModel(1)
         v_mean = model.mean([1, 10], [5, 10])
         self.assertEqual(v_mean[0], 10 / (1 + 9 * np.exp(-5)))
         self.assertEqual(v_mean[1], 10 / (1 + 9 * np.exp(-10)))
 
-        # Check model variance is not implemented
-        times = np.linspace(0, 100, 101)
-        parameters = [0.1, 50]
-        self.assertRaises(NotImplementedError, model.variance,
-                          parameters, times)
-
 
 if __name__ == '__main__':
     unittest.main()

pints/tests/test_toy_stochastic_michaelis_menten_model.py

@@ -37,6 +37,11 @@ def test_propensities(self):
                 model._propensities(x_0, k),
                 np.array([200.0, 4000.0, 4000.0])))
 
+    def test_n_outputs(self):
+        x_0 = [1e4, 2e3, 2e4, 0]
+        model = MichaelisMentenModel(x_0)
+        self.assertEqual(model.n_outputs(), 4)
+
 
 if __name__ == '__main__':
     unittest.main()

pints/tests/test_toy_stochastic_production_degradation_model.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env python3
+#
+# Tests if the production and degradation (toy) model works.
+#
+# This file is part of PINTS (https://github.com/pints-team/pints/) which is
+# released under the BSD 3-clause license. See accompanying LICENSE.md for
+# copyright notice and full license details.
+#
+import unittest
+import numpy as np
+from pints.toy.stochastic import ProductionDegradationModel
+
+
+class TestProductionDegradationModel(unittest.TestCase):
+    """
+    Tests if the degradation (toy) model works.
+    """
+    def test_n_parameters(self):
+        x_0 = 20
+        model = ProductionDegradationModel(x_0)
+        self.assertEqual(model.n_parameters(), 2)
+
+    def test_simulation_length(self):
+        x_0 = 20
+        model = ProductionDegradationModel(x_0)
+        times = np.linspace(0, 1, 100)
+        k = [0.1, 0.2]
+        values = model.simulate(k, times)
+        self.assertEqual(len(values), 100)
+
+    def test_propensities(self):
+        x_0 = 20
+        k = [0.1, 0.2]
+        model = ProductionDegradationModel(x_0)
+        self.assertTrue(
+            np.allclose(
+                model._propensities([x_0], k),
+                np.array([2.0, 0.2])))
+
+    def test_suggested(self):
+        model = ProductionDegradationModel(20)
+        times = model.suggested_times()
+        parameters = model.suggested_parameters()
+        self.assertTrue(len(times) == 101)
+        self.assertTrue(np.all(parameters > 0))
+
+
+if __name__ == '__main__':
+    unittest.main()

pints/tests/test_toy_stochastic_schlogl_model.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env python3
+#
+# Tests if the Schlogl (toy) model works.
+#
+# This file is part of PINTS (https://github.com/pints-team/pints/) which is
+# released under the BSD 3-clause license. See accompanying LICENSE.md for
+# copyright notice and full license details.
+#
+import unittest
+import numpy as np
+from pints.toy.stochastic import SchloglModel
+
+
+class TestSchloglModel(unittest.TestCase):
+    """
+    Tests if the degradation (toy) model works.
+    """
+    def test_n_parameters(self):
+        x_0 = 20
+        model = SchloglModel(x_0)
+        self.assertEqual(model.n_parameters(), 4)
+
+    def test_simulation_length(self):
+        x_0 = 20
+        model = SchloglModel(x_0)
+        times = np.linspace(0, 1, 100)
+        k = [0.1, 0.2, 0.3, 0.4]
+        values = model.simulate(k, times)
+        self.assertEqual(len(values), 100)
+
+    def test_propensities(self):
+        x_0 = 20
+        model = SchloglModel(x_0)
+        k = model.suggested_parameters()
+        self.assertTrue(
+            np.allclose(
+                model._propensities([x_0], k),
+                np.array([68.4, 1.71, 2200.0, 750.0])))
+
+    def test_suggested(self):
+        model = SchloglModel(20)
+        times = model.suggested_times()
+        parameters = model.suggested_parameters()
+        self.assertTrue(len(times) == 101)
+        self.assertTrue(np.all(parameters > 0))
+
+
+if __name__ == '__main__':
+    unittest.main()

pints/toy/__init__.py

@@ -32,4 +32,3 @@
 from ._simple_egg_box import SimpleEggBoxLogPDF
 from ._sir_model import SIRModel
 from ._twisted_gaussian_banana import TwistedGaussianLogPDF
-from ._stochastic_logistic_model import StochasticLogisticModel

pints/toy/stochastic/__init__.py

@@ -9,3 +9,6 @@
 from ._markov_jump_model import MarkovJumpModel  # noqa
 from ._michaelis_menten_model import MichaelisMentenModel # noqa
 from ._degradation_model import DegradationModel # noqa
+from ._logistic_model import LogisticModel # noqa
+from ._production_degradation_model import ProductionDegradationModel # noqa
+from ._schlogl_model import SchloglModel # noqa