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

@@ -99,6 +99,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_stochastic_logistic_model.py

@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+#
+# Tests if the stochastic logistic growth (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 StochasticLogisticModel
+
+
+class TestStochasticLogistic(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
+        model = 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 small simulation
+        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):
+        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):
+        times = np.linspace(0, 100, 101)
+        model = StochasticLogisticModel(1)
+        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)
+
+    def test_mean_variance(self):
+        # test mean
+        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)))
+
+    def test_errors(self):
+        model = pints.toy.StochasticLogisticModel(1)
+        # parameters, times cannot be negative
+        times = np.linspace(0, 100, 101)
+        parameters = [-0.1, 50]
+        self.assertRaises(ValueError, model.simulate, parameters, times)
+        self.assertRaises(ValueError, model.mean, parameters, times)
+
+        parameters = [0.1, -50]
+        self.assertRaises(ValueError, model.simulate, parameters, times)
+        self.assertRaises(ValueError, model.mean, parameters, 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)
+
+        # this model should have 2 parameters
+        parameters_3 = [0.1]
+        self.assertRaises(ValueError, model.simulate, parameters_3, times)
+        self.assertRaises(ValueError, model.mean, parameters_3, times)
+
+        # model variance isn't implemented so we should throw a helpful error
+        parameters_4 = [0.1, 50]
+        self.assertRaises(NotImplementedError, model.variance,
+                          parameters_4, times)
+
+        # Initial value can't be negative
+        self.assertRaises(ValueError, pints.toy.StochasticLogisticModel, -1)
+
+
+if __name__ == '__main__':
+    unittest.main()

pints/toy/__init__.py

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

pints/toy/_stochastic_logistic_model.py

@@ -0,0 +1,132 @@
+#
+# 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 StochasticLogisticModel(pints.ForwardModel, ToyModel):
+    r"""
+
+    *Extends:* :class:`pints.ForwardModel`, :class:`pints.toy.ToyModel`.
+    """
+
+    def __init__(self, initial_molecule_count=2):
+        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 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"""
+        Returns the deterministic mean of infinitely many stochastic
+        simulations, which follows:
+        :math:`\frac{kC(0)}{C(0) + \exp{-kt}(k - C(0))}`.
+        """
+        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, 50])
+
+    def suggested_times(self):
+        """ See "meth:`pints.toy.ToyModel.suggested_times()`."""
+        return np.linspace(0, 100, 101)