WIP - Stochastic Logistic Growth (ABC Toy problem)

docs/source/toy/index.rst

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

docs/source/toy/stochastic_logistic_model.rst

@@ -0,0 +1,7 @@
+*************************
+Stochastic Logistic Model
+*************************
+
+.. currentmodule:: pints.toy
+
+.. autoclass:: StochasticLogisticModel

examples/README.md

@@ -117,6 +117,7 @@ relevant code.
 - [Simple Harmonic Oscillator model](./toy/model-simple-harmonic-oscillator.ipynb)
 - [SIR Epidemiology model](./toy/model-sir.ipynb)
 - [Stochastic Degradation model](./toy/model-stochastic-degradation.ipynb)
+- [Stochastic Logistic model](./toy/model-stochastic-logistic-growth.ipynb)
 
 ### Distributions
 - [Annulus](./toy/distribution-annulus.ipynb)

pints/tests/test_toy_logistic_model.py

@@ -12,7 +12,7 @@
 import pints.toy
 
 
-class TestLogistic(unittest.TestCase):
+class TestLogisticModel(unittest.TestCase):
     """
     Tests if the logistic (toy) model works.
     """

pints/tests/test_toy_stochastic_degradation_model.py

@@ -13,7 +13,7 @@
 from pints.toy import StochasticDegradationModel
 
 
-class TestStochasticDegradation(unittest.TestCase):
+class TestStochasticDegradationModel(unittest.TestCase):
     """
     Tests if the stochastic degradation (toy) model works.
     """

pints/tests/test_toy_stochastic_logistic_model.py

@@ -0,0 +1,116 @@
+#!/usr/bin/env python3
+#
+# Tests if the stochastic logistic growth (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
+
+
+class TestStochasticLogisticModel(unittest.TestCase):
+    """
+    Tests if the stochastic logistic growth (toy) model works.
+    """
+    def test_start_with_zero(self):
+        # Test the special case where the initial population count is zero
+
+        # Set seed for random generator
+        np.random.seed(1)
+
+        model = pints.toy.StochasticLogisticModel(0)
+        times = [0, 1, 2, 100, 1000]
+        parameters = [0.1, 50]
+        values = model.simulate(parameters, times)
+        self.assertEqual(len(values), len(times))
+        self.assertTrue(np.all(values == np.zeros(5)))
+
+    def test_start_with_one(self):
+        # Run a small simulation and check it runs properly
+
+        # Set seed for random generator
+        np.random.seed(1)
+
+        model = pints.toy.StochasticLogisticModel(1)
+        times = [0, 1, 2, 100, 1000]
+        parameters = [0.1, 50]
+        values = model.simulate(parameters, times)
+        self.assertEqual(len(values), len(times))
+        self.assertEqual(values[0], 1)
+        self.assertEqual(values[-1], 50)
+        self.assertTrue(np.all(values[1:] >= values[:-1]))
+
+    def test_suggested(self):
+        # Check suggested values
+        model = pints.toy.StochasticLogisticModel(1)
+        times = model.suggested_times()
+        parameters = model.suggested_parameters()
+        self.assertTrue(len(times) == 101)
+        self.assertTrue(np.all(parameters > 0))
+
+    def test_simulate(self):
+        # Check each step in the simulation process
+        np.random.seed(1)
+        model = pints.toy.StochasticLogisticModel(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)
+        self.assertTrue(len(time), len(raw_values))
+
+        # Test output of Gillespie algorithm
+        self.assertTrue(np.all(raw_values == np.array(range(1, 51))))
+
+        # 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)
+        temp_time = np.array([np.random.uniform(time[1], time[2])])
+        self.assertTrue(model._interpolate_values(time, raw_values, temp_time,
+                                                  params)[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)
+        parameters_2 = [0.1, 50]
+        self.assertRaises(ValueError, model.simulate, parameters_2, times_2)
+        self.assertRaises(ValueError, model.mean, parameters_2, times_2)
+
+        # Check this model takes 2 parameters
+        parameters_3 = [0.1]
+        self.assertRaises(ValueError, model.simulate, parameters_3, times)
+        self.assertRaises(ValueError, model.mean, parameters_3, times)
+
+        # Check initial value cannot be negative
+        self.assertRaises(ValueError, pints.toy.StochasticLogisticModel, -1)
+
+    def test_mean_variance(self):
+        # Check the mean is what we expected
+        model = pints.toy.StochasticLogisticModel(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/toy/__init__.py

@@ -36,3 +36,4 @@
 from ._sir_model import SIRModel
 from ._twisted_gaussian_banana import TwistedGaussianLogPDF
 from ._stochastic_degradation_model import StochasticDegradationModel
+from ._stochastic_logistic_model import StochasticLogisticModel

pints/toy/_hes1_michaelis_menten.py

@@ -34,20 +34,20 @@ class Hes1Model(pints.ForwardModel, ToyModel):
 
     Extends :class:`pints.ForwardModel`, :class:`pints.toy.ToyModel`.
 
-    References
-    ----------
-    .. [1] Silk, D., el al. 2011. Designing attractive models via automated
-           identification of chaotic and oscillatory dynamical regimes. Nature
-           communications, 2, p.489.
-           https://doi.org/10.1038/ncomms1496
-
     Parameters
     ----------
     y0 : float
         The initial condition of the observable. Requires ``y0 >= 0``.
     implicit_parameters
         The implicit parameter of the model that is not inferred, given as a
         vector ``[p1_0, p2_0, k_deg]`` with ``p1_0, p2_0, k_deg >= 0``.
+
+    References
+    ----------
+    .. [1] Silk, D., el al. 2011. Designing attractive models via automated
+           identification of chaotic and oscillatory dynamical regimes. Nature
+           communications, 2, p.489.
+           https://doi.org/10.1038/ncomms1496
     """
     def __init__(self, y0=None, implicit_parameters=None):
         if y0 is None:

pints/toy/_stochastic_logistic_model.py

@@ -0,0 +1,163 @@
+#
+# Stochastic logistic model.
+#
+# 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.
+#
+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 StochasticLogisticModel(pints.ForwardModel, ToyModel):
+    r"""
+    This model describes the growth of a population of individuals, where the
+    birth rate per capita, initially :math:`b_0`, decreases to :math:`0` as the
+    population size, :math:`\mathcal{C}(t)`, starting from an initial
+    population size, :math:`n_0`, approaches a carrying capacity, :math:`k`.
+    This process follows a rate according to [1]_
+
+    .. math::
+       A \xrightarrow{b_0(1-\frac{\mathcal{C}(t)}{k})} 2A.
+
+    The model is simulated using the Gillespie stochastic simulation algorithm
+    [2]_, [3]_.
+
+    *Extends:* :class:`pints.ForwardModel`, :class:`pints.toy.ToyModel`.
+
+    Parameters
+    ----------
+    initial_molecule_count : float
+        Sets the initial population size :math:`n_0`.
+
+    References
+    ----------
+    .. [1] Simpson, M. et al. 2019. Process noise distinguishes between
+           indistinguishable population dynamics. bioRxiv.
+           https://doi.org/10.1101/533182
+    .. [2] Gillespie, D. 1976. A General Method for Numerically Simulating the
+           Stochastic Time Evolution of Coupled Chemical Reactions.
+           Journal of Computational Physics. 22 (4): 403-434.
+           https://doi.org/10.1016/0021-9991(76)90041-3
+    .. [3] Erban R. et al. 2007. A practical guide to stochastic simulations
+           of reaction-diffusion processes. arXiv.
+           https://arxiv.org/abs/0704.1908v2
+    """
+
+    def __init__(self, initial_molecule_count=50):
+        super(StochasticLogisticModel, 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 2
+
+    def _simulate_raw(self, parameters):
+        """
+        Returns tuple (raw times, population sizes) when reactions occur.
+        """
+        parameters = np.asarray(parameters)
+        if len(parameters) != self.n_parameters():
+            raise ValueError('This model should have only 2 parameters.')
+        b = parameters[0]
+        k = parameters[1]
+        if b <= 0:
+            raise ValueError('Rate constant must be positive.')
+
+        # Initial time and count
+        t = 0
+        a = self._n0
+
+        # Run stochastic logistic birth-only algorithm, calculating time until
+        # next reaction and increasing population count by 1 at that time
+        mol_count = [a]
+        time = [t]
+        while a < k:
+            r = np.random.uniform(0, 1)
+            t += np.log(1 / r) / (a * b * (1 - a / k))
+            a = a + 1
+            time.append(t)
+            mol_count.append(a)
+        return time, mol_count
+
+    def _interpolate_values(self, time, pop_size, output_times, parameters):
+        """
+        Takes raw times and population size values as inputs and outputs
+        interpolated values at output_times.
+        """
+        # Interpolate as step function, increasing pop_size by 1 at each
+        # event time point
+        interp_func = interp1d(time, pop_size, kind='previous')
+
+        # Compute population size values at given time points using f1
+        # at any time beyond the last event, pop_size = k
+        values = interp_func(output_times[np.where(output_times <= time[-1])])
+        zero_vector = np.full(
+            len(output_times[np.where(output_times > time[-1])]),
+            parameters[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, pop_size = self._simulate_raw(parameters)
+
+        # interpolate
+        values = self._interpolate_values(time, pop_size, times, parameters)
+        return values
+
+    def mean(self, parameters, times):
+        r"""
+        Computes the deterministic mean of infinitely many stochastic
+        simulations with times :math:`t` and parameters (:math:`b`, :math:`k`),
+        which follows:
+        :math:`\frac{kC(0)}{C(0) + (k - C(0)) \exp(-bt)}`.
+
+        Returns an array with the same length as `times`.
+        """
+        parameters = np.asarray(parameters)
+        if len(parameters) != self.n_parameters():
+            raise ValueError('This model should have only 2 parameters.')
+
+        b = parameters[0]
+        if b <= 0:
+            raise ValueError('Rate constant must be positive.')
+
+        k = parameters[1]
+        if k <= 0:
+            raise ValueError("Carrying capacity must be positive")
+
+        times = np.asarray(times)
+        if np.any(times < 0):
+            raise ValueError('Negative times are not allowed.')
+        c0 = self._n0
+        return (c0 * k) / (c0 + np.exp(-b * times) * (k - c0))
+
+    def variance(self, parameters, times):
+        r"""
+        Returns the deterministic variance of infinitely many stochastic
+        simulations.
+        """
+        raise NotImplementedError
+
+    def suggested_parameters(self):
+        """ See :meth:`pints.toy.ToyModel.suggested_parameters()`. """
+        return np.array([0.1, 500])
+
+    def suggested_times(self):
+        """ See :meth:`pints.toy.ToyModel.suggested_times()`."""
+        return np.linspace(0, 100, 101)