Generic stochastic model

docs/source/index.rst

@@ -39,6 +39,7 @@ Contents
     optimisers/index
     noise_model_diagnostics
     toy/index
+    toy/stochastic/index
     transformations
     utilities
 

docs/source/toy/index.rst

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

docs/source/toy/stochastic/index.rst

@@ -0,0 +1,15 @@
+***********************
+Stochastic Toy Problems
+***********************
+
+The `stochastic` module provides toy :class:`models<pints.toy.stochastic.MarkovJumpModel>`,
+:class:`distributions<pints.LogPrior>` and
+:class:`error measures<pints.ErrorMeasure>` that can be used for tests and in
+examples.
+
+
+.. toctree::
+
+    markov_jump_model
+    stochastic_degradation_model
+    stochastic_michaelis_menten_model

docs/source/toy/stochastic/markov_jump_model.rst

@@ -0,0 +1,7 @@
+*****************
+Markov Jump Model
+*****************
+
+.. currentmodule:: pints.toy.stochastic
+
+.. autoclass:: MarkovJumpModel

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

@@ -2,6 +2,6 @@
 Stochastic degradation model
 ****************************
 
-.. currentmodule:: pints.toy
+.. currentmodule:: pints.toy.stochastic
 
-.. autoclass:: StochasticDegradationModel
+.. autoclass:: DegradationModel

docs/source/toy/stochastic/stochastic_michaelis_menten_model.rst

@@ -0,0 +1,7 @@
+*********************************
+Stochastic Michaelis Menten model
+*********************************
+
+.. currentmodule:: pints.toy.stochastic
+
+.. autoclass:: MichaelisMentenModel

examples/README.md

@@ -117,6 +117,7 @@ relevant code.
 - [SIR Epidemiology model](./toy/model-sir.ipynb)
 - [Stochastic Degradation model](./toy/model-stochastic-degradation.ipynb)
 - [Stochastic Logistic model](./toy/model-stochastic-logistic-growth.ipynb)
+- [Michaelis Menten model](./toy/model-stochastic-michaelis-menten.ipynb)
 
 ### Distributions
 - [Annulus](./toy/distribution-annulus.ipynb)

pints/tests/test_toy_markov_jump_model.py

@@ -0,0 +1,82 @@
+#!/usr/bin/env python3
+#
+# Tests if the markov jump 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 DegradationModel
+
+
+class TestMarkovJumpModel(unittest.TestCase):
+    """
+    Tests if the markov jump model works using
+    the degradation model.
+    """
+    def test_start_with_zero(self):
+        # Test the special case where the initial molecule count is zero
+        model = DegradationModel(0)
+        times = [0, 1, 2, 100, 1000]
+        parameters = [0.1]
+        values = model.simulate(parameters, times)
+        self.assertEqual(len(values), len(times))
+        self.assertTrue(np.all(values == np.zeros(5)))
+
+    def test_start_with_twenty(self):
+        # Run small simulation
+        model = DegradationModel(20)
+        times = [0, 1, 2, 100, 1000]
+        parameters = [0.1]
+        values = model.simulate(parameters, times)
+        self.assertEqual(len(values), len(times))
+        self.assertEqual(values[0], 20)
+        self.assertEqual(values[-1], 0)
+        self.assertTrue(np.all(values[1:] <= values[:-1]))
+
+    def test_simulate(self):
+        times = np.linspace(0, 100, 101)
+        model = DegradationModel(20)
+        time, mol_count = model.simulate_raw([0.1], 100)
+        values = model.interpolate_mol_counts(time, mol_count, times)
+        self.assertTrue(len(time), len(mol_count))
+        # Test output of Gillespie algorithm
+        expected = np.array([[x] for x in range(20, -1, -1)])
+        self.assertTrue(np.all(mol_count == expected))
+
+        # Check simulate function returns expected values
+        self.assertTrue(np.all(values[np.where(times < time[1])] == 20))
+
+        # Check interpolation function works as expected
+        temp_time = np.array([np.random.uniform(time[0], time[1])])
+        self.assertEqual(
+            model.interpolate_mol_counts(time, mol_count, temp_time)[0],
+            20)
+        temp_time = np.array([np.random.uniform(time[1], time[2])])
+        self.assertEqual(
+            model.interpolate_mol_counts(time, mol_count, temp_time)[0],
+            19)
+
+    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_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
+        parameters_3 = [0.1, 1]
+        self.assertRaises(ValueError, model.simulate, parameters_3, times)
+
+        # Initial value can't be negative
+        self.assertRaises(ValueError, DegradationModel, -1)
+
+
+if __name__ == '__main__':
+    unittest.main()

pints/tests/test_toy_stochastic_degradation_model.py

@@ -1,116 +1,88 @@
 #!/usr/bin/env python3
 #
-# Tests if the stochastic degradation (toy) model works.
+# Tests if the 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
-import pints
-import pints.toy
-from pints.toy import StochasticDegradationModel
+from pints.toy.stochastic import DegradationModel
 
 
-class TestStochasticDegradationModel(unittest.TestCase):
+class TestDegradationModel(unittest.TestCase):
     """
-    Tests if the stochastic degradation (toy) model works.
+    Tests if the degradation (toy) model works.
     """
-    def test_start_with_zero(self):
-        # Test the special case where the initial molecule count is zero
-        model = StochasticDegradationModel(0)
-        times = [0, 1, 2, 100, 1000]
-        parameters = [0.1]
-        values = model.simulate(parameters, times)
-        self.assertEqual(len(values), len(times))
-        self.assertTrue(np.all(values == np.zeros(5)))
-
-    def test_start_with_twenty(self):
-        # Run small simulation
-        model = pints.toy.StochasticDegradationModel(20)
-        times = [0, 1, 2, 100, 1000]
-        parameters = [0.1]
-        values = model.simulate(parameters, times)
-        self.assertEqual(len(values), len(times))
-        self.assertEqual(values[0], 20)
-        self.assertEqual(values[-1], 0)
-        self.assertTrue(np.all(values[1:] <= values[:-1]))
+    def test_n_parameters(self):
+        x_0 = 20
+        model = DegradationModel(x_0)
+        self.assertEqual(model.n_parameters(), 1)
+
+    def test_simulation_length(self):
+        x_0 = 20
+        model = DegradationModel(x_0)
+        times = np.linspace(0, 1, 100)
+        k = [0.1]
+        values = model.simulate(k, times)
+        self.assertEqual(len(values), 100)
+
+    def test_propensities(self):
+        x_0 = 20
+        k = [0.1]
+        model = DegradationModel(x_0)
+        self.assertTrue(
+            np.allclose(
+                model._propensities([x_0], k),
+                np.array([2.0])))
 
     def test_suggested(self):
-        model = pints.toy.StochasticDegradationModel(20)
+        model = DegradationModel(20)
         times = model.suggested_times()
         parameters = model.suggested_parameters()
         self.assertTrue(len(times) == 101)
         self.assertTrue(parameters > 0)
 
-    def test_simulate(self):
-        times = np.linspace(0, 100, 101)
-        model = StochasticDegradationModel(20)
-        time, mol_count = model.simulate_raw([0.1])
-        values = model.interpolate_mol_counts(time, mol_count, times)
-        self.assertTrue(len(time), len(mol_count))
-        # Test output of Gillespie algorithm
-        self.assertTrue(np.all(mol_count == np.array(range(20, -1, -1))))
-
-        # Check simulate function returns expected values
-        self.assertTrue(np.all(values[np.where(times < time[1])] == 20))
-
-        # Check interpolation function works as expected
-        temp_time = np.array([np.random.uniform(time[0], time[1])])
-        self.assertEqual(
-            model.interpolate_mol_counts(time, mol_count, temp_time)[0],
-            20)
-        temp_time = np.array([np.random.uniform(time[1], time[2])])
-        self.assertEqual(
-            model.interpolate_mol_counts(time, mol_count, temp_time)[0],
-            19)
-
     def test_mean_variance(self):
         # test mean
-        model = pints.toy.StochasticDegradationModel(10)
+        model = DegradationModel(10)
         v_mean = model.mean([1], [5, 10])
         self.assertEqual(v_mean[0], 10 * np.exp(-5))
         self.assertEqual(v_mean[1], 10 * np.exp(-10))
 
-        model = pints.toy.StochasticDegradationModel(20)
+        model = DegradationModel(20)
         v_mean = model.mean([5], [7.2])
         self.assertEqual(v_mean[0], 20 * np.exp(-7.2 * 5))
 
         # test variance
-        model = pints.toy.StochasticDegradationModel(10)
+        model = DegradationModel(10)
         v_var = model.variance([1], [5, 10])
         self.assertEqual(v_var[0], 10 * (np.exp(5) - 1.0) / np.exp(10))
         self.assertAlmostEqual(v_var[1], 10 * (np.exp(10) - 1.0) / np.exp(20))
 
-        model = pints.toy.StochasticDegradationModel(20)
+        model = DegradationModel(20)
         v_var = model.variance([2.0], [2.0])
         self.assertAlmostEqual(v_var[0], 20 * (np.exp(4) - 1.0) / np.exp(8))
 
     def test_errors(self):
-        model = pints.toy.StochasticDegradationModel(20)
+        model = DegradationModel(20)
         # parameters, times cannot be negative
         times = np.linspace(0, 100, 101)
         parameters = [-0.1]
-        self.assertRaises(ValueError, model.simulate, parameters, times)
         self.assertRaises(ValueError, model.mean, parameters, times)
         self.assertRaises(ValueError, model.variance, parameters, times)
 
         times_2 = np.linspace(-10, 10, 21)
         parameters_2 = [0.1]
-        self.assertRaises(ValueError, model.simulate, parameters_2, times_2)
         self.assertRaises(ValueError, model.mean, parameters_2, times_2)
         self.assertRaises(ValueError, model.variance, parameters_2, times_2)
 
         # this model should have 1 parameter
         parameters_3 = [0.1, 1]
-        self.assertRaises(ValueError, model.simulate, parameters_3, times)
         self.assertRaises(ValueError, model.mean, parameters_3, times)
         self.assertRaises(ValueError, model.variance, parameters_3, times)
 
-        # Initial value can't be negative
-        self.assertRaises(ValueError, pints.toy.StochasticDegradationModel, -1)
-
 
 if __name__ == '__main__':
     unittest.main()

pints/tests/test_toy_stochastic_michaelis_menten_model.py

@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+#
+# Tests if the Michaelis Menten (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 MichaelisMentenModel
+
+
+class TestMichaelisMentenModel(unittest.TestCase):
+    """
+    Tests if the Michaelis Menten (toy) model works.
+    """
+    def test_n_parameters(self):
+        x_0 = [1e4, 2e3, 2e4, 0]
+        model = MichaelisMentenModel(x_0)
+        self.assertEqual(model.n_parameters(), 3)
+
+    def test_simulation_length(self):
+        x_0 = [1e4, 2e3, 2e4, 0]
+        model = MichaelisMentenModel(x_0)
+        times = np.linspace(0, 1, 100)
+        k = [1e-5, 0.2, 0.2]
+        values = model.simulate(k, times)
+        self.assertEqual(len(values), 100)
+
+    def test_propensities(self):
+        x_0 = [1e4, 2e3, 2e4, 0]
+        k = [1e-5, 0.2, 0.2]
+        model = MichaelisMentenModel(x_0)
+        self.assertTrue(
+            np.allclose(
+                model._propensities(x_0, k),
+                np.array([200.0, 4000.0, 4000.0])))
+
+
+if __name__ == '__main__':
+    unittest.main()

pints/toy/__init__.py

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