Merge pull request #884 from pints-team/toy_models-daniel

Add stochastic degradation toy model. See #890

docs/source/toy/index.rst

@@ -30,8 +30,8 @@ Some toy classes provide extra functionality defined in the
     parabolic_error
     repressilator_model
     rosenbrock
-    sir_model
     simple_egg_box_logpdf
+    sir_model
+    stochastic_degradation_model
     toy_classes
     twisted_gaussian_logpdf
-

docs/source/toy/stochastic_degradation_model.rst

@@ -0,0 +1,7 @@
+****************************
+Stochastic degradation model
+****************************
+
+.. currentmodule:: pints.toy
+
+.. autoclass:: StochasticDegradationModel

examples/README.md

@@ -84,6 +84,7 @@ relevant code.
 - [Lotka-Volterra predator-prey model](./toy-model-lotka-volterra.ipynb)
 - [Repressilator model](./toy-model-repressilator.ipynb)
 - [SIR Epidemiology model](./toy-model-sir.ipynb)
+- [Stochastic Degradation model](./toy-model-stochastic-degradation.ipynb)
 
 ### Distributions
 - [Annulus](./toy-distribution-annulus.ipynb)

pints/tests/test_toy_stochastic_degradation_model.py

@@ -0,0 +1,116 @@
+#!/usr/bin/env python3
+#
+# Tests if the stochastic degradation (toy) model works.
+#
+# This file is part of PINTS.
+#  Copyright (c) 2017-2019, University of Oxford.
+#  For licensing information, see the LICENSE file distributed with the PINTS
+#  software package.
+#
+import unittest
+import numpy as np
+import pints
+import pints.toy
+from pints.toy import StochasticDegradationModel
+
+
+class TestStochasticDegradation(unittest.TestCase):
+    """
+    Tests if the stochastic 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_suggested(self):
+        model = pints.toy.StochasticDegradationModel(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.assertTrue(model.interpolate_mol_counts(time, mol_count,
+                                                     temp_time)[0] == 20)
+        temp_time = np.array([np.random.uniform(time[1], time[2])])
+        self.assertTrue(model.interpolate_mol_counts(time, mol_count,
+                                                     temp_time)[0] == 19)
+
+    def test_mean_variance(self):
+        # test mean
+        model = pints.toy.StochasticDegradationModel(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)
+        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)
+        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)
+        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)
+        # 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/toy/__init__.py

@@ -32,3 +32,4 @@
 from ._simple_egg_box import SimpleEggBoxLogPDF                     # noqa
 from ._sir_model import SIRModel                                    # noqa
 from ._twisted_gaussian_banana import TwistedGaussianLogPDF         # noqa
+from ._stochastic_degradation_model import StochasticDegradationModel  # noqa

pints/toy/_stochastic_degradation_model.py

@@ -0,0 +1,179 @@
+#
+# Stochastic degradation toy model.
+#
+# This file is part of PINTS.
+#  Copyright (c) 2017-2019, University of Oxford.
+#  For licensing information, see the LICENSE file distributed with the PINTS
+#  software package.
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+import numpy as np
+from scipy.interpolate import interp1d
+import pints
+
+from . import ToyModel
+
+
+class StochasticDegradationModel(pints.ForwardModel, ToyModel):
+    r"""
+    Stochastic degradation model of a single chemical reaction starting from
+    an initial molecule count :math:`A(0)` and degrading to 0 with a fixed rate
+    :math:`k`:
+
+    .. math::
+        A \xrightarrow{k} 0
+
+    Simulations are performed using the Gillespie algorithm [1, 2]:
+
+    1. Sample a random value :math:`r` from a uniform distribution
+
+    .. math::
+        r \sim U(0,1)
+
+    2. Calculate the time :math:`\tau` until the next single reaction as
+
+    .. math::
+        \tau = \frac{-\ln(r)}{A(t) k}
+
+    3. Update the molecule count :math:`A` at time :math:`t + \tau` as:
+
+    .. math::
+        A(t + \tau) = A(t) - 1
+
+    4. Return to step (1) until the molecule count reaches 0
+
+    The model has one parameter, the rate constant :math:`k`.
+
+    The initial molecule count :math:`A(0)` can be set using the optional
+    constructor argument ``initial_molecule_count``
+
+    [1] A Practical Guide to Stochastic Simulations of Reaction Diffusion
+    Processes. Erban, Chapman, Maini (2007). arXiv:0704.1908v2 [q-bio.SC]
+    https://arxiv.org/abs/0704.1908
+
+    [2] A general method for numerically simulating the stochastic time
+    evolution of coupled chemical reactions. Gillespie (1976).
+    Journal of Computational Physics
+    https://doi.org/10.1016/0021-9991(76)90041-3
+
+    *Extends:* :class:`pints.ForwardModel`, :class:`pints.toy.ToyModel`.
+    """
+    def __init__(self, initial_molecule_count=20):
+        super(StochasticDegradationModel, self).__init__()
+        self._n0 = float(initial_molecule_count)
+        if self._n0 < 0:
+            raise ValueError('Initial molecule count cannot be negative.')
+
+    def n_parameters(self):
+        """ See :meth:`pints.ForwardModel.n_parameters()`. """
+        return 1
+
+    def simulate_raw(self, parameters):
+        """
+        Returns raw times, mol counts when reactions occur
+        """
+        parameters = np.asarray(parameters)
+        if len(parameters) != self.n_parameters():
+            raise ValueError('This model should have only 1 parameter.')
+        k = parameters[0]
+
+        if k <= 0:
+            raise ValueError('Rate constant must be positive.')
+
+        # Initial time and count
+        t = 0
+        a = self._n0
+
+        # Run stochastic degradation algorithm, calculating time until next
+        # reaction and decreasing molecule count by 1 at that time
+        mol_count = [a]
+        time = [t]
+        while a > 0:
+            r = np.random.uniform(0, 1)
+            t += -np.log(r) / (a * k)
+            a = a - 1
+            time.append(t)
+            mol_count.append(a)
+        return time, mol_count
+
+    def interpolate_mol_counts(self, time, mol_count, output_times):
+        """
+        Takes raw times and inputs and mol counts and outputs interpolated
+        values at output_times
+        """
+        # Interpolate as step function, decreasing mol_count by 1 at each
+        # reaction time point
+        interp_func = interp1d(time, mol_count, kind='previous')
+
+        # Compute molecule count values at given time points using f1
+        # at any time beyond the last reaction, molecule count = 0
+        values = interp_func(output_times[np.where(output_times <= time[-1])])
+        zero_vector = np.zeros(
+            len(output_times[np.where(output_times > time[-1])])
+        )
+        values = np.concatenate((values, zero_vector))
+        return values
+
+    def simulate(self, parameters, times):
+        """ See :meth:`pints.ForwardModel.simulate()`. """
+        times = np.asarray(times)
+        if np.any(times < 0):
+            raise ValueError('Negative times are not allowed.')
+        if self._n0 == 0:
+            return np.zeros(times.shape)
+
+        # run Gillespie
+        time, mol_count = self.simulate_raw(parameters)
+
+        # interpolate
+        values = self.interpolate_mol_counts(time, mol_count, times)
+        return values
+
+    def mean(self, parameters, times):
+        r"""
+        Returns the deterministic mean of infinitely many stochastic
+        simulations, which follows :math:`A(0) \exp(-kt)`.
+        """
+        parameters = np.asarray(parameters)
+        if len(parameters) != self.n_parameters():
+            raise ValueError('This model should have only 1 parameter.')
+        k = parameters[0]
+
+        if k <= 0:
+            raise ValueError('Rate constant must be positive.')
+
+        times = np.asarray(times)
+        if np.any(times < 0):
+            raise ValueError('Negative times are not allowed.')
+
+        mean = self._n0 * np.exp(-k * times)
+        return mean
+
+    def variance(self, parameters, times):
+        r"""
+        Returns the deterministic variance of infinitely many stochastic
+        simulations, which follows :math:`\exp(-2kt)(-1 + \exp(kt))A(0)`.
+        """
+        parameters = np.asarray(parameters)
+        if len(parameters) != self.n_parameters():
+            raise ValueError('This model should have only 1 parameter.')
+        k = parameters[0]
+
+        if k <= 0:
+            raise ValueError('Rate constant must be positive.')
+
+        times = np.asarray(times)
+        if np.any(times < 0):
+            raise ValueError('Negative times are not allowed.')
+
+        variance = np.exp(-2 * k * times) * (-1 + np.exp(k * times)) * self._n0
+        return variance
+
+    def suggested_parameters(self):
+        """ See :meth:`pints.toy.ToyModel.suggested_parameters()`. """
+        return np.array([0.1])
+
+    def suggested_times(self):
+        """ See "meth:`pints.toy.ToyModel.suggested_times()`."""
+        return np.linspace(0, 100, 101)