Repository: DavAug/ErlotinibGefitinib

Issue:
SharedNoiseLogLikelihood
Implemement a pints.LogPDF that allows to share the noise parameters between problems.

Discussion:

Patch:
diff --git a/pkpd/__init__.py b/pkpd/__init__.py
index bc620e0..52948fa 100644
--- a/pkpd/__init__.py
+++ b/pkpd/__init__.py
@@ -9,7 +9,8 @@ import os
 
 from ._likelihoods import (  # noqa
     ConstantAndMultiplicativeGaussianLogLikelihood,
-    FixedEtaLogLikelihoodWrapper
+    FixedEtaLogLikelihoodWrapper,
+    PooledLogPDF
 )
 
 from ._models import (  # noqa
diff --git a/pkpd/_likelihoods.py b/pkpd/_likelihoods.py
index d5c1e93..043afa7 100644
--- a/pkpd/_likelihoods.py
+++ b/pkpd/_likelihoods.py
@@ -273,3 +273,173 @@ class FixedEtaLogLikelihoodWrapper(pints.LogPDF):
 
     def n_parameters(self):
         return self._log_pdf.n_parameters() - 1
+
+
+class PooledLogPDF(pints.LogPDF):
+    r"""
+    Calculates a sum of :class:`LogPDF` objects, while pooling selected
+    parameters across log-pdfs.
+
+    This :class:`LogPDF` requires that all :class:`LogPDF` objects have the
+    same number of parameters.
+
+    This is useful for e.g. modelling the time-series of multiple individuals
+    (each individual defines a separate :class:`LogPDF`), and some parameters
+    are expected to be the same across individuals. This may be a realistic
+    expectation for the noise parameter, if the experimental set up didn't vary
+    for the indiviudals.
+
+    For two :class:`LogPDF` objects :math:`L _1`  and :math:`L _2` with
+    parameters :math:`(\psi _1, \theta _1`) and :math:`(\psi _2, \theta _2`)
+    respectively, a pooling of the :math:`\theta _i` results in a pooled
+    log-pdf of the form
+
+    .. math::
+        \log L(\psi _1, \psi _2, \theta | D_1, D_2) =
+            \log L(\psi _1, \theta | D_1) + \log L(\psi _2, \theta | D_2),
+
+    where :math:`\theta := \theta _1 = \theta _2`, and $D_i$ is the measured
+    time-series of individual :math:`i`.
+
+    In general, the search of a :class:`PooledLogPDF` has dimensionality
+    :math:`nd_{\text{up}} + d_{\text{p}}`, where :math:`n` is the number of
+    individuals, :math:`d_{\text{up}}` is the dimension of the unpooled
+    :math:`psi` parameters, and :math:`d_{\text{up}}` is the dimension of the
+    pooled :math:`\theta` parameter.
+
+    The order of parameters for evaluation is largely kept, and the parameters
+    of individual log-pdfs are concatenated. However, pooled parameters are
+    always prepended. Consider for example two log-pdf with parameters
+    :math:`(\psi _{1,1}, \psi _{1, 2}, \psi _{1,3}, \psi _{1,4})` and
+    :math:`(\psi _{2,1}, \psi _{2, 2}, \psi _{2,3}, \psi _{2,4})`. Pooling the
+    first and the last parameter, i.e.
+    :math:`\theta _1 := \psi _{1,1} = \psi _{2, 1}` and
+    :math:`\theta _4 := \psi _{1,4} = \psi _{2, 4}` results in and expect
+    order of parameters
+    :math:`(\psi _{1, 2}, \psi _{1,3}, \psi _{2, 2}, \psi _{2,3}, \theta _1,
+    \theta _4)`.
+
+    Extends :class:`LogPDF`.
+
+    Parameters
+    ----------
+    log_pdfs
+        A sequence of :class:`LogPDF` objects.
+    pooled
+        An array-like object of dtype bool, indicating which parameters across
+        the likelihoods are pooled (`True`) or remain unpooled (`False`).
+
+    Example
+    -------
+    ::
+        pooled_log_likelihood = pints.PooledLogPDFs(
+            log_pdfs=[
+                pints.GaussianLogLikelihood(problem1),
+                pints.GaussianLogLikelihood(problem2)],
+            pooled=[False, True])
+    """
+    def __init__(self, log_pdfs, pooled):
+        super(PooledLogPDF, self).__init__()
+
+        # Check input arguments
+        if len(log_pdfs) < 2:
+            raise ValueError(
+                'PooledLogPDF requires at least two log-pdfs.')
+        for index, pdf in enumerate(log_pdfs):
+            if not isinstance(pdf, pints.LogPDF):
+                raise ValueError(
+                    'All log-pdfs passed to PooledLogPDFs must be instances of'
+                    ' pints.LogPDF (failed on argument '
+                    + str(index) + ').')
+
+        # Check parameter dimension across log-pdfs
+        self._log_pdfs = log_pdfs
+        n_parameters = self._log_pdfs[0].n_parameters()
+        for pdf in self._log_pdfs:
+            if pdf.n_parameters() != n_parameters:
+                raise ValueError(
+                    'All log-pdfs passed to PooledLogPDFs must have '
+                    'same dimension.')
+
+        # Check that pooled matches number of parameters
+        self._pooled = np.asarray(pooled)
+        if len(self._pooled) != n_parameters:
+            raise ValueError(
+                'The array-like input `pooled` needs to have the same length '
+                'as the number of parameters of the individual log-pdfs.')
+
+        # Check that pooled contains only booleans
+        for p in self._pooled:
+            if not isinstance(p, np.bool_):
+                raise ValueError(
+                    'The array-like input `pooled` passed to PooledLogPDFs '
+                    'has to contain booleans exclusively.')
+
+        # Get dimension of search space
+        self._n_pooled = np.sum(self._pooled)
+        n_individuals = len(self._log_pdfs)
+        self._n_unpooled = np.sum(~self._pooled)
+        self._n_parameters = \
+            self._n_pooled + n_individuals * self._n_unpooled
+
+    def __call__(self, parameters):
+        # Get parameters of pooled log-pdf
+        parameters = np.asarray(parameters)
+
+        # Create container for parameters of individuals log-pdf and fill with
+        # pooled parameters
+        params_ind = np.empty(shape=self._n_unpooled + self._n_pooled)
+        if self._n_pooled > 0:
+            params_ind[self._pooled] = parameters[-self._n_pooled:]
+
+        # Compute pdf score
+        total = 0
+        for idx, pdf in enumerate(self._log_pdfs):
+            # Get unpooled parameters for individual
+            params_ind[~self._pooled] = parameters[
+                idx * self._n_unpooled: (idx + 1) * self._n_unpooled]
+
+            # Compute pdf score contribution
+            total += pdf(params_ind)
+        return total
+
+    def evaluateS1(self, parameters):
+        """
+        See :meth:`LogPDF.evaluateS1()`.
+
+        *This method only works if all the underlying :class:`LogPDF` objects
+        implement the optional method :meth:`LogPDF.evaluateS1()`!*
+        """
+        # Get parameters of pooled log-pdf
+        parameters = np.asarray(parameters)
+
+        # Create container for parameters of individuals log-pdf and fill with
+        # pooled parameters
+        params_ind = np.empty(shape=self._n_unpooled + self._n_pooled)
+        if self._n_pooled > 0:
+            params_ind[self._pooled] = parameters[-self._n_pooled:]
+
+        # Compute pdf score and partials
+        total = 0
+        dtotal = np.zeros(shape=self._n_parameters)
+        for idx, pdf in enumerate(self._log_pdfs):
+            # Get unpooled parameters for individual
+            params_ind[~self._pooled] = parameters[
+                idx * self._n_unpooled: (idx + 1) * self._n_unpooled]
+
+            # Compute pdf score and partials for individual
+            score, partials = pdf.evaluateS1(params_ind)
+
+            # Add contributions to score and partials. Note that partials
+            # w.r.t. unpooled parameters receive only one contribution
+            total += score
+            dtotal[idx * self._n_unpooled: (idx + 1) * self._n_unpooled] = \
+                partials[~self._pooled]
+            if self._n_pooled > 0:
+                dtotal[-self._n_pooled:] += partials[self._pooled]
+
+        return total, dtotal
+
+    def n_parameters(self):
+        """ See :meth:`LogPDF.n_parameters()`. """
+        return self._n_parameters


Test patch:
diff --git a/pkpd/tests/test_likelihoods.py b/pkpd/tests/test_likelihoods.py
index 90fac88..9a5b564 100644
--- a/pkpd/tests/test_likelihoods.py
+++ b/pkpd/tests/test_likelihoods.py
@@ -617,5 +617,220 @@ class TestFixedEtaLogLikelihoodWrapper(unittest.TestCase):
         self.assertEqual(wrapped_log_likelihood.n_parameters(), 2)
 
 
+class TestPooledLogPDF(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        # Create test log-pdfs
+        model = pints.toy.ConstantModel(1)
+
+        problem = pints.SingleOutputProblem(
+            model=model, times=[1, 2, 3, 4], values=[1, 2, 3, 4])
+        cls.log_pdf_1 = pints.GaussianLogLikelihood(problem)
+
+        problem = pints.SingleOutputProblem(
+            model=model, times=[1, 2, 3, 4], values=[1, 1, 1, 1])
+        cls.log_pdf_2 = pints.GaussianLogLikelihood(problem)
+
+    def test_bad_number_log_pdfs(self):
+        log_pdfs = [self.log_pdf_1]
+        pooled = [True, True]
+        self.assertRaisesRegex(
+            ValueError, 'PooledLogPDF requires', pkpd.PooledLogPDF, log_pdfs,
+            pooled)
+
+    def test_bad_log_pdfs_objects(self):
+        log_pdfs = ['log_pdf_1', 'log_pdf_2']
+        pooled = [True, True]
+        self.assertRaisesRegex(
+            ValueError, 'All log-pdfs passed', pkpd.PooledLogPDF, log_pdfs,
+            pooled)
+
+    def test_bad_log_pdfs_parameters(self):
+        model = pints.toy.ConstantModel(1)
+        problem = pints.SingleOutputProblem(
+            model=model, times=[1, 2, 3, 4], values=[1, 2, 3, 4])
+        log_pdf = pkpd.ConstantAndMultiplicativeGaussianLogLikelihood(problem)
+
+        log_pdfs = [self.log_pdf_1, log_pdf]
+        pooled = [True, True]
+
+        self.assertRaisesRegex(
+            ValueError, 'All log-pdfs passed to PooledLogPDFs',
+            pkpd.PooledLogPDF, log_pdfs, pooled)
+
+    def test_bad_pooled_length(self):
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+        pooled = [True, True, True]
+        self.assertRaisesRegex(
+            ValueError, 'The array-like input `pooled` needs',
+            pkpd.PooledLogPDF, log_pdfs, pooled)
+
+    def test_bad_pooled_content(self):
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+        pooled = [True, 'Yes']
+        self.assertRaisesRegex(
+            ValueError, 'The array-like input `pooled` passed',
+            pkpd.PooledLogPDF, log_pdfs, pooled)
+
+    def test_n_parameters(self):
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+
+        # Pool nothing
+        pooled = [False, False]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        n_parameters = \
+            self.log_pdf_1.n_parameters() + self.log_pdf_2.n_parameters()
+        self.assertEqual(log_pdf.n_parameters(), n_parameters)
+
+        # Pool first parameter
+        pooled = [True, False]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        n_parameters = 2 * 1 + 1
+        self.assertEqual(log_pdf.n_parameters(), n_parameters)
+
+        # Pool second parameter
+        pooled = [False, True]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        n_parameters = 2 * 1 + 1
+        self.assertEqual(log_pdf.n_parameters(), n_parameters)
+
+        # Pool both parameters
+        pooled = [True, True]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        n_parameters = 2 * 1
+        self.assertEqual(log_pdf.n_parameters(), n_parameters)
+
+    def test_call_unpooled(self):
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+        pooled = [False, False]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        param_1 = [1, 0.2]
+        param_2 = [2, 0.1]
+
+        score = self.log_pdf_1(param_1) + self.log_pdf_2(param_2)
+        self.assertEqual(log_pdf(param_1 + param_2), score)
+
+    def test_call_partially_pooled(self):
+        # Pool first parameter
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+        pooled = [True, False]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        param_1 = [1, 0.2]
+        param_2 = [1, 0.1]
+        param = [0.2, 0.1, 1]
+
+        score = self.log_pdf_1(param_1) + self.log_pdf_2(param_2)
+        self.assertEqual(log_pdf(param), score)
+
+        # Pool second parameter
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+        pooled = [False, True]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        param_1 = [1, 0.1]
+        param_2 = [2, 0.1]
+        param = [1, 2, 0.1]
+
+        score = self.log_pdf_1(param_1) + self.log_pdf_2(param_2)
+        self.assertEqual(log_pdf(param), score)
+
+    def test_call_pooled(self):
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+        pooled = [True, True]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        param = [1, 0.1]
+
+        score = self.log_pdf_1(param) + self.log_pdf_2(param)
+        self.assertEqual(log_pdf(param), score)
+
+    def test_evaluateS1_unpooled(self):
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+        pooled = [False, False]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        param_1 = [1, 0.2]
+        param_2 = [2, 0.1]
+
+        score_1, dscore_1 = self.log_pdf_1.evaluateS1(param_1)
+        score_2, dscore_2 = self.log_pdf_2.evaluateS1(param_2)
+
+        score, dscore = log_pdf.evaluateS1(param_1 + param_2)
+        self.assertEqual(score, score_1 + score_2)
+
+        self.assertEqual(len(dscore), 4)
+        self.assertEqual(dscore[0], dscore_1[0])
+        self.assertEqual(dscore[1], dscore_1[1])
+        self.assertEqual(dscore[2], dscore_2[0])
+        self.assertEqual(dscore[3], dscore_2[1])
+
+    def test_evaluateS1_partially_pooled(self):
+        # Pool first parameter
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+        pooled = [True, False]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        param_1 = [1, 0.2]
+        param_2 = [1, 0.1]
+        param = [0.2, 0.1, 1]
+
+        score_1, dscore_1 = self.log_pdf_1.evaluateS1(param_1)
+        score_2, dscore_2 = self.log_pdf_2.evaluateS1(param_2)
+
+        score, dscore = log_pdf.evaluateS1(param)
+        self.assertEqual(score, score_1 + score_2)
+
+        self.assertEqual(len(dscore), 3)
+        self.assertEqual(dscore[0], dscore_1[1])
+        self.assertEqual(dscore[1], dscore_2[1])
+        self.assertEqual(dscore[2], dscore_1[0] + dscore_2[0])
+
+        # Pool second parameter
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+        pooled = [False, True]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        param_1 = [1, 0.1]
+        param_2 = [2, 0.1]
+        param = [1, 2, 0.1]
+
+        score_1, dscore_1 = self.log_pdf_1.evaluateS1(param_1)
+        score_2, dscore_2 = self.log_pdf_2.evaluateS1(param_2)
+
+        score, dscore = log_pdf.evaluateS1(param)
+        self.assertEqual(score, score_1 + score_2)
+
+        self.assertEqual(len(dscore), 3)
+        self.assertEqual(dscore[0], dscore_1[0])
+        self.assertEqual(dscore[1], dscore_2[0])
+        self.assertEqual(dscore[2], dscore_1[1] + dscore_2[1])
+
+    def test_evaluateS1_pooled(self):
+        log_pdfs = [self.log_pdf_1, self.log_pdf_2]
+        pooled = [True, True]
+        log_pdf = pkpd.PooledLogPDF(log_pdfs, pooled)
+
+        param_1 = [1, 0.2]
+        param_2 = [1, 0.2]
+        param = [1, 0.2]
+
+        score_1, dscore_1 = self.log_pdf_1.evaluateS1(param_1)
+        score_2, dscore_2 = self.log_pdf_2.evaluateS1(param_2)
+
+        score, dscore = log_pdf.evaluateS1(param)
+        self.assertEqual(score, score_1 + score_2)
+
+        self.assertEqual(len(dscore), 2)
+        self.assertEqual(dscore[0], dscore_1[0] + dscore_2[0])
+        self.assertEqual(dscore[1], dscore_1[1] + dscore_2[1])
+
+
 if __name__ == '__main__':
     unittest.main()