Generic stochastic model

CHANGELOG.md

@@ -5,6 +5,7 @@ All notable changes to this project will be documented in this file.
 ## Unreleased
 
 ### Added
+- [#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.
 - [#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`.
 
 ### Changed

LICENSE.md

@@ -1,6 +1,6 @@
 BSD 3-Clause License
 
-Copyright (c) 2017-2021, University of Oxford (University of Oxford means the
+Copyright (c) 2017-2022, University of Oxford (University of Oxford means the
 Chancellor, Masters and Scholars of the University of Oxford, having an
 administrative office at Wellington Square, Oxford OX1 2JD, UK).
 All rights reserved.

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

@@ -103,7 +103,7 @@ relevant code.
 
 ## Toy problems
 
-### Models
+### Deterministic Models
 - [Beeler-Reuter action potential model](./toy/model-beeler-reuter-ap.ipynb)
 - [Constant model](./toy/model-constant.ipynb)
 - [Fitzhugh-Nagumo model](./toy/model-fitzhugh-nagumo.ipynb)
@@ -115,8 +115,11 @@ relevant code.
 - [Repressilator model](./toy/model-repressilator.ipynb)
 - [Simple Harmonic Oscillator model](./toy/model-simple-harmonic-oscillator.ipynb)
 - [SIR Epidemiology model](./toy/model-sir.ipynb)
+
+### Stochastic Models
 - [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)
 
 ### 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