From 8a6f6afe874217ace3ed1182a817fd80d3157760 Mon Sep 17 00:00:00 2001 From: DavAug Date: Sun, 20 Mar 2022 16:04:13 +0000 Subject: [PATCH] #255 first draft of docs quick overview --- chi/__init__.py | 4 +- chi/_inference.py | 6 +- chi/_mechanistic_models.py | 375 +++++++++--------- chi/library/_model_library_api.py | 30 +- docs/requirements.txt | 1 + docs/source/api/index.rst | 8 +- docs/source/api/mechanistic_models.rst | 33 +- .../getting_started/code/1_simulation_1.py | 341 ++++++++++++++++ docs/source/getting_started/error_model.rst | 5 + .../images/1_simulation_1.html | 71 ++++ .../images/1_simulation_2.html | 71 ++++ .../images/1_simulation_3.html | 71 ++++ .../images/1_simulation_4.html | 71 ++++ .../images/1_simulation_5.html | 71 ++++ docs/source/getting_started/index.rst | 20 +- .../source/getting_started/log_likelihood.rst | 5 + docs/source/getting_started/log_posterior.rst | 5 + .../getting_started/mechanistic_model.rst | 5 + .../source/getting_started/quick_overview.rst | 276 +++++++++++++ docs/source/index.rst | 2 +- 20 files changed, 1231 insertions(+), 240 deletions(-) create mode 100644 docs/source/getting_started/code/1_simulation_1.py create mode 100644 docs/source/getting_started/error_model.rst create mode 100644 docs/source/getting_started/images/1_simulation_1.html create mode 100644 docs/source/getting_started/images/1_simulation_2.html create mode 100644 docs/source/getting_started/images/1_simulation_3.html create mode 100644 docs/source/getting_started/images/1_simulation_4.html create mode 100644 docs/source/getting_started/images/1_simulation_5.html create mode 100644 docs/source/getting_started/log_likelihood.rst create mode 100644 docs/source/getting_started/log_posterior.rst create mode 100644 docs/source/getting_started/mechanistic_model.rst create mode 100644 docs/source/getting_started/quick_overview.rst diff --git a/chi/__init__.py b/chi/__init__.py index dc8620e0..07883bbc 100644 --- a/chi/__init__.py +++ b/chi/__init__.py @@ -29,8 +29,8 @@ from ._mechanistic_models import ( # noqa MechanisticModel, - PharmacodynamicModel, - PharmacokineticModel, + SBMLModel, + PKPDModel, ReducedMechanisticModel ) diff --git a/chi/_inference.py b/chi/_inference.py index 52e8f409..5198fbeb 100644 --- a/chi/_inference.py +++ b/chi/_inference.py @@ -748,7 +748,7 @@ def __init__(self, log_posterior): super(SamplingController, self).__init__(log_posterior) # Set default sampler - self._sampler = pints.HaarioACMC + self._sampler = pints.HaarioBardenetACMC def _format_chains(self, chains, names, ids, divergent_iters): """ @@ -776,8 +776,8 @@ def _format_chains(self, chains, names, ids, divergent_iters): trajectories occured, or None """ # Broadcast IDs to length of names, if posterior has only one ID - if isinstance(ids, str): - ids = np.broadcast_to(ids, shape=len(names)) + if (not ids) or isinstance(ids, str): + ids = [ids] * len(names) # Convert names and ids to numpy arrays ids = np.asarray(ids) diff --git a/chi/_mechanistic_models.py b/chi/_mechanistic_models.py index 85ff8dbe..123cccb8 100644 --- a/chi/_mechanistic_models.py +++ b/chi/_mechanistic_models.py @@ -13,18 +13,114 @@ class MechanisticModel(object): + r""" + A base class for mechanistic time series models of the form + + .. math: + \bar{y} = g(\bar{y}, \psi, t), + + where :math:`\bar{y}` are the model outputs, :math`\psi` are the model + parameters and :math:`t` is the time. :math:`g` can be any (deterministic) + function. + """ + def __init__(self): + super(MechanisticModel, self).__init__() + + def copy(self): + """ + Returns a deep copy of the mechanistic model. + + .. note: + The default implementation may have to be adapted based on the + implementation of :meth:`simulate`. + """ + return copy.deepcopy(self) + + def enable_sensitivities(self, enabled): + r""" + Enables the computation of the model output sensitivities to the model + parameters if set to ``True``. + + The sensitivities of the model outputs are defined as the partial + derviatives of the ouputs :math:`\bar{y}` with respect to the model + parameters :math:`\psi` + + .. math: + \frac{\del \bar{y}}{\del \psi}. + + :param enabled: A boolean flag which enables (``True``) / disables + (``False``) the computation of sensitivities. + :type enabled: bool + """ + raise NotImplementedError + + def has_sensitivities(self): + """ + Returns a boolean indicating whether sensitivities have been enabled. + """ + raise NotImplementedError + + def n_outputs(self): + """ + Returns the number of output dimensions. + + By default this is the number of states. + """ + raise NotImplementedError + + def n_parameters(self): + """ + Returns the number of parameters in the model. + + Parameters of the model are initial state values and structural + parameter values. + """ + raise NotImplementedError + + def outputs(self): + """ + Returns the output names of the model. + """ + raise NotImplementedError + + def parameters(self): + """ + Returns the parameter names of the model. + """ + raise NotImplementedError + + def simulate(self, parameters, times): + """ + Returns the numerical solution of the model outputs (and optionally + the sensitivites) for the specified parameters and times. + + The model outputs are returned as a 2 dimensional NumPy array of shape + (n_outputs, n_times). If sensitivities are enabled, a tuple is returned + with the NumPy array of the model outputs and a NumPy array of the + sensitivities of shape (n_times, n_outputs, n_parameters). + + :param parameters: An array-like object with values for the model + parameters. + :type parameters: list, numpy.ndarray + :param times: An array-like object with time points at which the output + values are returned. + :type times: list, numpy.ndarray + """ + raise NotImplementedError + + +class SBMLModel(MechanisticModel): """ - A base class for models that are specified by sbml files. + Instantiates a mechanistic model from a SBML specification. - Parameters - ---------- - sbml_file - A path to the SBML model file that specifies the model. + Extends :class:`MechanisticModel`. + :param sbml_file: A path to the SBML model file that specifies the model. + :type sbml_file: str """ def __init__(self, sbml_file): - super(MechanisticModel, self).__init__() + super(SBMLModel, self).__init__() # Import model self._model = sbml.SBMLImporter().model(sbml_file) @@ -37,7 +133,7 @@ def __init__(self, sbml_file): # Create simulator without sensitivities # (intentionally public property) - self.simulator = myokit.Simulation(self._model) + self._simulator = myokit.Simulation(self._model) self._has_sensitivities = False def _get_time_unit(self): @@ -58,7 +154,7 @@ def _set_const(self, parameters): Sets values of constant model parameters. """ for id_var, var in enumerate(self._const_names): - self.simulator.set_constant(var, float(parameters[id_var])) + self._simulator.set_constant(var, float(parameters[id_var])) def _set_state(self, parameters): """ @@ -66,7 +162,7 @@ def _set_state(self, parameters): """ parameters = np.array(parameters) parameters = parameters[self._original_order] - self.simulator.set_state(parameters) + self._simulator.set_state(parameters) def _set_number_and_names(self): """ @@ -119,7 +215,7 @@ def copy(self): """ # Copy model manually and get protocol myokit_model = self._model.clone() - protocol = self.simulator._protocol + protocol = self._simulator._protocol # Copy the mechanistic model model = copy.deepcopy(self) @@ -154,13 +250,13 @@ def enable_sensitivities(self, enabled, parameter_names=None): enabled = bool(enabled) # Get dosing regimen from existing simulator - protocol = self.simulator._protocol + protocol = self._simulator._protocol if not enabled: if self._has_sensitivities: # Disable sensitivities sim = myokit.Simulation(self._model, protocol) - self.simulator = sim + self._simulator = sim self._has_sensitivities = False return None @@ -199,7 +295,7 @@ def enable_sensitivities(self, enabled, parameter_names=None): sim = myokit.Simulation(self._model, protocol, sensitivities) # Update simulator and sensitivity state - self.simulator = sim + self._simulator = sim self._has_sensitivities = True def has_sensitivities(self): @@ -271,7 +367,7 @@ def set_outputs(self, outputs): # Check that outputs are valid for output in outputs: try: - var = self.simulator._model.get(output) + var = self._simulator._model.get(output) if not (var.is_state() or var.is_intermediary()): raise ValueError( 'Outputs have to be state or intermediary variables.') @@ -396,7 +492,7 @@ def simulate(self, parameters, times): :type times: list, numpy.ndarray """ # Reset simulation - self.simulator.reset() + self._simulator.reset() # Set initial conditions self._set_state(parameters[:self._n_states]) @@ -406,13 +502,13 @@ def simulate(self, parameters, times): # Simulate if not self._has_sensitivities: - output = self.simulator.run( + output = self._simulator.run( times[-1] + 1, log=self._output_names, log_times=times) output = np.array([output[name] for name in self._output_names]) return output - output, sensitivities = self.simulator.run( + output, sensitivities = self._simulator.run( times[-1] + 1, log=self._output_names, log_times=times) output = np.array([output[name] for name in self._output_names]) sensitivities = np.array(sensitivities) @@ -426,87 +522,27 @@ def time_unit(self): return self._time_unit -class PharmacodynamicModel(MechanisticModel): +class PKPDModel(SBMLModel): """ - Converts a pharmacodynamic model specified by an SBML file into a forward - model that can be solved numerically. - - Extends :class:`MechanisticModel`. + Instantiates a PKPD model from a SBML specification. - Parameters - ---------- - sbml_file - A path to the SBML model file that specifies the pharmacodynamic model. + Extends :class:`SBMLModel`. + :param sbml_file: A path to the SBML model file that specifies the + PKPD model. + :type sbml_file: str """ def __init__(self, sbml_file): - super(PharmacodynamicModel, self).__init__(sbml_file) - - # Set default pharmacokinetic input variable - # (Typically drug concentration) - self._pk_input = None - if self._model.has_variable('myokit.drug_concentration'): - self._pk_input = 'myokit.drug_concentration' - - def pk_input(self): - """ - Returns the pharmacokinetic input variable. In most models this will be - the concentration of the drug. - - Defaults to ``None`` or ``myokit.drug_concentration`` if the latter is - among the model parameters. - """ - return self._pk_input - - def set_pk_input(self, name): - """ - Sets the pharmacokinetic input variable. In most models this will be - the concentration of the drug. - - The name has to match a parameter of the model. - """ - if name not in self._parameter_names: - raise ValueError( - 'The name does not match a model parameter.') - - self._pk_input = name - - -class PharmacokineticModel(MechanisticModel): - """ - Converts a pharmacokinetic model specified by an SBML file into a forward - model that can be solved numerically. - - Extends :class:`MechanisticModel`. - - Parameters - ---------- - sbml_file - A path to the SBML model file that specifies the pharmacokinetic model. - - """ - - def __init__(self, sbml_file): - super(PharmacokineticModel, self).__init__(sbml_file) + super(PKPDModel, self).__init__(sbml_file) # Set default dose administration self._administration = None + self._dosing_regimen = None # Safe vanilla model self._vanilla_model = self._model.clone() - # Set default output variable that interacts with the pharmacodynamic - # model - # (Typically drug concentration in central compartment) - self._pd_output = None - if self._model.has_variable('central.drug_concentration'): - self._pd_output = 'central.drug_concentration' - - # Set default output to pd output if not None - if self._pd_output is not None: - self.set_outputs([self._pd_output]) - def _add_dose_compartment(self, model, drug_amount): """ Adds a dose compartment to the model with a linear absorption rate to @@ -599,7 +635,7 @@ def dosing_regimen(self): :class:`myokit.Protocol`. If the protocol has not been set, ``None`` is returned. """ - return self.simulator._protocol + return self._dosing_regimen def set_administration( self, compartment, amount_var='drug_amount', direct=True): @@ -684,7 +720,7 @@ def set_administration( # Update model and simulator # (otherwise simulator won't know about pace bound variable) self._model = model - self.simulator = myokit.Simulation(model) + self._simulator = myokit.Simulation(model) self._has_sensitivities = False # Remember type of administration @@ -692,7 +728,7 @@ def set_administration( {'compartment': compartment, 'direct': direct}) def set_dosing_regimen( - self, dose, start, duration=0.01, period=None, num=None): + self, dose, start=0, duration=0.01, period=None, num=None): """ Sets the dosing regimen with which the compound is administered. @@ -708,23 +744,23 @@ def set_dosing_regimen( By default the dose is administered once. To apply multiple doses provide a dose administration period. - Parameters - ---------- - dose - The amount of the compound that is injected at each administration. - start - Start time of the treatment. - duration - Duration of dose administration. For a bolus injection, a dose - duration of 1% of the time unit should suffice. By default the - duration is set to 0.01 (bolus). - period - Periodicity at which doses are administered. If ``None`` the dose - is administered only once. - num - Number of administered doses. If ``None`` and the periodicity of - the administration is not ``None``, doses are administered - indefinitely. + :param dose: The amount of the compound that is injected at each + administration, or a myokit.Protocol instance that defines the + dosing regimen. + :type dose: float or myokit.Protocol + :param start: Start time of the treatment. By default the + administration starts at t=0. + :type start: float, optional + :param duration: Duration of dose administration. By default the + duration is set to 0.01 of the time unit (bolus). + :type duration: float, optional + :param period: Periodicity at which doses are administered. If ``None`` + the dose is administered only once. + :type period: float, optional + :param num: Number of administered doses. If ``None`` and the + periodicity of the administration is not ``None``, doses are + administered indefinitely. + :type num: int, optional """ if self._administration is None: raise ValueError( @@ -740,6 +776,11 @@ def set_dosing_regimen( period = 0 num = 0 + if isinstance(dose, myokit.Protocol): + self._simulator.set_protocol(dose) + self._dosing_regimen = dose + return None + # Translate dose to dose rate dose_rate = dose / duration @@ -747,61 +788,20 @@ def set_dosing_regimen( dosing_regimen = myokit.pacing.blocktrain( period=period, duration=duration, offset=start, level=dose_rate, limit=num) - self.simulator.set_protocol(dosing_regimen) - - def pd_output(self): - """ - Returns the variable which interacts with the pharmacodynamic model. - In most models this will be the concentration of the drug in the - central compartment. + self._simulator.set_protocol(dosing_regimen) + self._dosing_regimen = dosing_regimen - This variable is mapped to the - :meth:`chi.PharmacodynamicModel.pk_input` variable when a - :class:`PKPDModel` is instantiated. - - Defaults to ``None`` or ``central.drug_concentration`` if the latter is - among the model parameters. - """ - return self._pd_output - def set_pd_output(self, name): - """ - Sets the variable which interacts with the pharmacodynamic model. - In most models this will be the concentration of the drug in the - central compartment. - - The name has to match a parameter of the model. - - This variable is mapped to the - :meth:`chi.PharmacodynamicModel.pk_input` variable when a - :class:`PKPDModel` is instantiated. - """ - # Get intermediate variable names - inter_names = [ - var.qname() for var in self._model.variables(inter=True)] - - names = inter_names + self._parameter_names - if name not in names: - raise ValueError( - 'The name does not match a model variable.') - - self._pd_output = name - - -class ReducedMechanisticModel(object): +class ReducedMechanisticModel(MechanisticModel): """ - A class that can be used to permanently fix model parameters of a - :class:`MechanisticModel` instance. + A wrapper class for a :class:`MechanisticModel` instance that can be used + to fix model parameters to fixed values. - This may be useful to explore simplified versions of a model before - defining a new SBML file. + Extends :class:`MechanisticModel`. - Parameters - ---------- - mechanistic_model - An instance of a :class:`MechanisticModel`. + :param mechanistic_model: A mechanistic model. + :type mechanistic_model: chi.MechanisticModel """ - def __init__(self, mechanistic_model): super(ReducedMechanisticModel, self).__init__() @@ -812,7 +812,6 @@ def __init__(self, mechanistic_model): 'chi.MechanisticModel') self._mechanistic_model = mechanistic_model - self.simulator = mechanistic_model.simulator # Set defaults self._fixed_params_mask = None @@ -837,7 +836,6 @@ def copy(self): # Replace mechanistic model and simulator model._mechanistic_model = mechanistic_model - model.simulator = mechanistic_model.simulator return model @@ -990,41 +988,6 @@ def parameters(self): return copy.copy(names) - def pd_output(self): - """ - Returns the variable which interacts with the pharmacodynamic model. - In most models this will be the concentration of the drug in the - central compartment. - - This variable is mapped to the - :meth:`chi.PharmacodynamicModel.pk_input` variable when a - :class:`PKPDModel` is instantiated. - - Defaults to ``None`` or ``central.drug_concentration`` if the latter is - among the model parameters. - - If the model does not support a pd output, ``None`` is returned. - """ - try: - return self._mechanistic_model.pd_output() - except AttributeError: - return None - - def pk_input(self): - """ - Returns the pharmacokinetic input variable. In most models this will be - the concentration of the drug. - - Defaults to ``None`` or ``myokit.drug_concentration`` if the latter is - among the model parameters. - - If the model does not support a pk input, ``None`` is returned. - """ - try: - return self._mechanistic_model.pk_input() - except AttributeError: - return None - def set_dosing_regimen( self, dose, start, duration=0.01, period=None, num=None): """ @@ -1077,7 +1040,11 @@ def set_outputs(self, outputs): A list of quantifiable variable names of the :class:`myokit.Model`, e.g. `compartment.variable`. """ - self._mechanistic_model.set_outputs(outputs) + try: + self._mechanistic_model.set_outputs(outputs) + except AttributeError: + raise AttributeError( + 'The mechanistic model does not support setting outputs.') def set_output_names(self, names): """ @@ -1089,7 +1056,11 @@ def set_output_names(self, names): names A dictionary that maps the current output names to new names. """ - self._mechanistic_model.set_output_names(names) + try: + self._mechanistic_model.set_output_names(names) + except AttributeError: + raise AttributeError( + 'The mechanistic model does not support setting output names.') def set_parameter_names(self, names): """ @@ -1101,9 +1072,13 @@ def set_parameter_names(self, names): names A dictionary that maps the current parameter names to new names. """ - # Set parameter names - self._mechanistic_model.set_parameter_names(names) - self._parameter_names = self._mechanistic_model.parameters() + try: + self._mechanistic_model.set_parameter_names(names) + self._parameter_names = self._mechanistic_model.parameters() + except AttributeError: + raise AttributeError( + 'The mechanistic model does not support setting parameter ' + 'names.') def simulate(self, parameters, times): """ @@ -1134,4 +1109,10 @@ def time_unit(self): """ Returns the model's unit of time. """ - return self._mechanistic_model.time_unit() + try: + time_unit = self._mechanistic_model.time_unit() + except AttributeError: + raise AttributeError( + 'The mechanistic model does not define a time unit.') + + return time_unit diff --git a/chi/library/_model_library_api.py b/chi/library/_model_library_api.py index 8d545bf1..2f135b9b 100644 --- a/chi/library/_model_library_api.py +++ b/chi/library/_model_library_api.py @@ -7,6 +7,8 @@ import os +import chi + class ModelLibrary(object): """ @@ -30,22 +32,17 @@ def __init__(self): def erlotinib_tumour_growth_inhibition_model(self): """ - .. warning:: - This model is going to be deprecated soon in favour of a dynamic - PK model-PD model composition. - This model is a combination of a :meth:`ModelLibrary.one_compartment_pk_model` and a :meth:`tumour_growth_inhibition_model_koch_reparametrised`. """ file_name = 'temporary_full_pkpd_model.xml' - return self._path + file_name + return chi.PKPDModel(self._path + file_name) def one_compartment_pk_model(self): r""" - Returns the absolute path to a SBML file, specifying a one compartment - pharmacokinetic model. + Returns an instantiation of a 1-compartment PK model. In this model the distribution of the drug is modelled by one compartment with a linear elimination rate :math:`k_e` @@ -62,13 +59,16 @@ def one_compartment_pk_model(self): the drug through the liver may be approximated by an exponential decay with the rate :math:`k_e`. - With a :class:`erlotinib.PharmacokineticModel` the drug may be either - directly administered to :math:`A` or indirectly through a dosing - compartment. + The drug may be either directly administered to :math:`A` or indirectly + through a dosing compartment. + + :rtype: chi.PKPDModel """ file_name = 'pk_one_comp.xml' + model = chi.PKPDModel(self._path + file_name) + model.set_outputs(['central.drug_concentration']) - return self._path + file_name + return model def tumour_growth_inhibition_model_koch(self): r""" @@ -100,10 +100,12 @@ def tumour_growth_inhibition_model_koch(self): .. math:: V_{\text{crit}} = \frac{\lambda _1}{2\lambda _0}. + + :rtype: chi.SBMLModel """ file_name = 'tgi_Koch_2009.xml' - return self._path + file_name + return chi.SBMLModel(self._path + file_name) def tumour_growth_inhibition_model_koch_reparametrised(self): r""" @@ -136,7 +138,9 @@ def tumour_growth_inhibition_model_koch_reparametrised(self): .. math:: V_{\text{crit}} = \frac{\lambda _1}{2\lambda _0} \quad \text{and} \quad \lambda = 2\lambda _1 . + + :rtype: chi.SBMLModel """ file_name = 'tgi_Koch_2009_reparametrised.xml' - return self._path + file_name + return chi.SBMLModel(self._path + file_name) diff --git a/docs/requirements.txt b/docs/requirements.txt index eef37541..e9a0202c 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,3 +1,4 @@ . +kaleido sphinx>=1.5, !=1.7.3 sphinx-rtd-theme \ No newline at end of file diff --git a/docs/source/api/index.rst b/docs/source/api/index.rst index 718d2b36..7dab94a0 100644 --- a/docs/source/api/index.rst +++ b/docs/source/api/index.rst @@ -9,7 +9,7 @@ ordered list of all classes and functions. The menu bar can be used to explore chi's functionality in a more structured way. .. toctree:: - :hidden: + :caption: Table of Contents :maxdepth: 1 covariate_models @@ -23,6 +23,8 @@ used to explore chi's functionality in a more structured way. predictive_models problems +Summary of all functions and classes in chi + .. autosummary:: @@ -49,8 +51,7 @@ used to explore chi's functionality in a more structured way. chi.MultiplicativeGaussianErrorModel chi.OptimisationController chi.PAMPredictiveModel - chi.PharmacodynamicModel - chi.PharmacokineticModel + chi.PKPDModel chi.plots.MarginalPosteriorPlot chi.plots.ParameterEstimatePlot chi.plots.PDTimeSeriesPlot @@ -70,4 +71,5 @@ used to explore chi's functionality in a more structured way. chi.ReducedMechanisticModel chi.ReducedPopulationModel chi.SamplingController + chi.SBMLModel chi.TruncatedGaussianModel diff --git a/docs/source/api/mechanistic_models.rst b/docs/source/api/mechanistic_models.rst index 0b18eb49..25223ec8 100644 --- a/docs/source/api/mechanistic_models.rst +++ b/docs/source/api/mechanistic_models.rst @@ -1,4 +1,3 @@ -.. _chi: https://github.com/DavAug/chi .. _SBML: http://sbml.org ****************** @@ -7,41 +6,35 @@ Mechanistic Models .. currentmodule:: chi -Mechanistic models in chi_ model the pharmacokinetics and +Mechanistic models in chi model the pharmacokinetics and pharmacodynamics of patients based on models specified by SBML files (System Biology Markup Language (SBML_)). Some SBML files relevant to the modelling of chi are provided in the :class:`ModelLibrary`. -Functional classes ------------------- +Classes +------- -- :class:`PharmacodynamicModel` -- :class:`PharmacokineticModel` +- :class:`MechanisticModel` +- :class:`SBMLModel` +- :class:`PKPDModel` +- :class:`ReducedMechanisticModel` Detailed API ^^^^^^^^^^^^ -.. autoclass:: PharmacodynamicModel +.. autoclass:: MechanisticModel :members: - :inherited-members: -.. autoclass:: PharmacokineticModel +.. autoclass:: SBMLModel :members: :inherited-members: -Base classes ------------- - -- :class:`MechanisticModel` -- :class:`ReducedMechanisticModel` - -Detailed API -^^^^^^^^^^^^ - -.. autoclass:: MechanisticModel +.. autoclass:: PKPDModel :members: + :inherited-members: .. autoclass:: ReducedMechanisticModel - :members: \ No newline at end of file + :members: + :inherited-members: \ No newline at end of file diff --git a/docs/source/getting_started/code/1_simulation_1.py b/docs/source/getting_started/code/1_simulation_1.py new file mode 100644 index 00000000..0156a732 --- /dev/null +++ b/docs/source/getting_started/code/1_simulation_1.py @@ -0,0 +1,341 @@ +## 1. Start example one +import chi.library + +# Define 1-compartment pharmacokinetic model +mechanistic_model = chi.library.ModelLibrary().one_compartment_pk_model() + +# Define parameters +parameters = [ + 10, # Initial drug amount + 2, # Volume of the compartment + 1, # Elimination rate +] + +# Define evaluation times +times = [0, 0.2, 0.5, 0.6, 1] + +# Simulate the model +result = mechanistic_model.simulate(parameters, times) +## end 1. + +## Start 2. +import numpy as np +import plotly.graph_objects as go + +# Set administration and dosing regimen +mechanistic_model.set_administration( + compartment='central', amount_var='drug_amount') +mechanistic_model.set_dosing_regimen(dose=1, period=0.5) + +# Simulate the model +times = np.linspace(start=0, stop=10, num=1000) +result = mechanistic_model.simulate(parameters, times) + +# Visualise results +fig = go.Figure() +fig.add_trace(go.Scatter( + x=times, + y=result[0], + mode='lines' +)) +fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' +) +fig.show() +## End 2. + +# Export for docs +import os + +directory = os.path.dirname(os.path.dirname(__file__)) +fig.write_html(directory + '/images/1_simulation_1.html') + +## Start 3. +import chi + +# Define error model and error model parameters +error_model = chi.GaussianErrorModel() +parameters = [0.2] # Sigma + +# Down sample times and mechanistic model evaluations +measurement_times = times[::25] +corresponding_evaluations = result[0, ::25] + +# Simulate measurements +measurements = error_model.sample( + parameters, model_output=corresponding_evaluations, seed=1)[:, 0] + +# Visualise results +fig = go.Figure() +fig.add_trace(go.Scatter( + x=times, + y=result[0], + mode='lines', + line_color='black', + name='Mechanistic model' +)) +fig.add_trace(go.Scatter( + x=measurement_times, + y=measurements, + mode='markers', + name='Measurements' +)) +fig.update_layout( + xaxis_title='Time', + yaxis_title='Drug concentration', + template='plotly_white' +) +fig.show() +# End 3. + +fig.write_html(directory + '/images/1_simulation_2.html') + +# Start 4. +# Define log-likelihood +log_likelihood = chi.LogLikelihood( + mechanistic_model, + error_model, + observations=measurements, + times=measurement_times +) + +# Evaluate log-likelihood for made-up parameters +madeup_parameters = [2, 7, 1.4, 3] +score_1 = log_likelihood(madeup_parameters) + +# Evaluate log-likelihood for data-generating parameters +true_parameters = [10, 2, 1, 0.2] +score_2 = log_likelihood(true_parameters) +# End 4. + +# Start 5. +import pints + +# Run optimisation +initial_parameters = [9, 3, 5, 1] +parameters_mle, score = pints.optimise( + log_likelihood, initial_parameters, method=pints.CMAES) +# End 5. + +parameters_mle = [10.26564936, 2.01524534, 1.00148417, 0.18456719] +# Start 6. +mle_output = mechanistic_model.simulate(parameters_mle, times) +fig.add_trace(go.Scatter( + x=times, + y=mle_output[0], + mode='lines', + line_dash='dash', + name='Inferred model' +)) +fig.show() +# End 6. +fig.write_html(directory + '/images/1_simulation_3.html') + +# Start 7. +# Define log-posterior +log_prior = pints.ComposedLogPrior( + pints.UniformLogPrior(0, 20), # Initial drug amount + pints.UniformLogPrior(0, 20), # Compartment volume + pints.UniformLogPrior(0, 20), # Elimination rate + pints.UniformLogPrior(0, 20) # Sigma +) +log_posterior = chi.LogPosterior(log_likelihood, log_prior) + +# Evaluate log-posterior +score_1 = log_posterior(madeup_parameters) +score_2 = log_posterior(true_parameters) +# End 7. + +# Start 8. +from plotly.subplots import make_subplots + +# Run inference +controller = chi.SamplingController(log_posterior) +n_iterations = 5000 +posterior_samples = controller.run(n_iterations) + +# Discard warmup iterations +warmup = 3000 +posterior_samples = posterior_samples.sel(draw=slice(warmup, n_iterations)) + +# Visualise Markov chains +fig = make_subplots(rows=4, cols=1, shared_xaxes=True) + +# Add MCMC chains for initial drug amount +for idc, chain in enumerate(posterior_samples['central.drug_amount'].values): + fig.add_trace( + go.Scatter( + legendgroup='Initial drug amount', + legendgrouptitle_text='Initial drug amount', + name='Chain %d' % idc, + x=np.arange(warmup, n_iterations) + 1, + y=chain, + showlegend=False + ), + row=1, + col=1 + ) +# Add MCMC chains for compartment volume +for idc, chain in enumerate(posterior_samples['central.size'].values): + fig.add_trace( + go.Scatter( + legendgroup='Compartment volume', + legendgrouptitle_text='Compartment volume', + name='Chain %d' % idc, + x=np.arange(warmup, n_iterations) + 1, + y=chain, + showlegend=False + ), + row=2, + col=1 + ) +# Add MCMC chains of elimination rate +for idc, chain in enumerate( + posterior_samples['myokit.elimination_rate'].values): + fig.add_trace( + go.Scatter( + legendgroup='Elimination rate', + legendgrouptitle_text='Elimination rate', + name='Chain %d' % idc, + x=np.arange(warmup, n_iterations) + 1, + y=chain, + showlegend=False + ), + row=3, + col=1 + ) +# Add MCMC chains of sigma +for idc, chain in enumerate(posterior_samples['Sigma'].values): + fig.add_trace( + go.Scatter( + legendgroup='Sigma', + legendgrouptitle_text='Sigma', + name='Chain %d' % idc, + x=np.arange(warmup, n_iterations) + 1, + y=chain, + showlegend=False + ), + row=4, + col=1 + ) + +fig.update_layout( + xaxis4_title='MCMC iteration', + yaxis_title='a_0', + yaxis2_title='v', + yaxis3_title='k_e', + yaxis4_title='sigma', + template='plotly_white' +) +fig.show() +# End 8. +fig.write_html(directory + '/images/1_simulation_4.html') + +# Start 9. +fig = make_subplots(rows=2, cols=2, shared_yaxes=True) +fig.add_trace( + go.Histogram( + name='Posterior samples', + x=posterior_samples['central.drug_amount'].values.flatten(), + histnorm='probability', + showlegend=False + ), + row=1, + col=1 +) +fig.add_trace( + go.Scatter( + name='Data-generating parameters', + x=[true_parameters[0]]*2, + y=[0, 0.04], + mode='lines', + line_color='black', + showlegend=False + ), + row=1, + col=1 +) + +fig.add_trace( + go.Histogram( + name='Posterior samples', + x=posterior_samples['central.size'].values.flatten(), + histnorm='probability', + showlegend=False + ), + row=1, + col=2 +) +fig.add_trace( + go.Scatter( + name='Data-generating parameters', + x=[true_parameters[1]]*2, + y=[0, 0.04], + mode='lines', + line_color='black', + showlegend=False + ), + row=1, + col=2 +) + +fig.add_trace( + go.Histogram( + name='Posterior samples', + x=posterior_samples['myokit.elimination_rate'].values.flatten(), + histnorm='probability', + showlegend=False + ), + row=2, + col=1 +) +fig.add_trace( + go.Scatter( + name='Data-generating parameters', + x=[true_parameters[2]]*2, + y=[0, 0.04], + mode='lines', + line_color='black', + showlegend=False + ), + row=2, + col=1 +) + +fig.add_trace( + go.Histogram( + name='Posterior samples', + x=posterior_samples['Sigma'].values.flatten(), + histnorm='probability', + showlegend=False + ), + row=2, + col=2 +) +fig.add_trace( + go.Scatter( + name='Data-generating parameters', + x=[true_parameters[3]]*2, + y=[0, 0.04], + mode='lines', + line_color='black', + showlegend=False + ), + row=2, + col=2 +) + +fig.update_layout( + xaxis_title='Initial drug amount', + xaxis2_title='Compartment volume', + xaxis3_title='Elimination rate', + xaxis4_title='Sigma', + yaxis_title='Probability', + yaxis3_title='Probability', + template='plotly_white' +) +fig.show() +# End 9. +fig.write_html(directory + '/images/1_simulation_5.html') diff --git a/docs/source/getting_started/error_model.rst b/docs/source/getting_started/error_model.rst new file mode 100644 index 00000000..4ff178e2 --- /dev/null +++ b/docs/source/getting_started/error_model.rst @@ -0,0 +1,5 @@ +*********************** +Defining an error model +*********************** + +To do. \ No newline at end of file diff --git a/docs/source/getting_started/images/1_simulation_1.html b/docs/source/getting_started/images/1_simulation_1.html new file mode 100644 index 00000000..c804a5ef --- /dev/null +++ b/docs/source/getting_started/images/1_simulation_1.html @@ -0,0 +1,71 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/1_simulation_2.html b/docs/source/getting_started/images/1_simulation_2.html new file mode 100644 index 00000000..806779ba --- /dev/null +++ b/docs/source/getting_started/images/1_simulation_2.html @@ -0,0 +1,71 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/1_simulation_3.html b/docs/source/getting_started/images/1_simulation_3.html new file mode 100644 index 00000000..19555853 --- /dev/null +++ b/docs/source/getting_started/images/1_simulation_3.html @@ -0,0 +1,71 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/1_simulation_4.html b/docs/source/getting_started/images/1_simulation_4.html new file mode 100644 index 00000000..f9241176 --- /dev/null +++ b/docs/source/getting_started/images/1_simulation_4.html @@ -0,0 +1,71 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/images/1_simulation_5.html b/docs/source/getting_started/images/1_simulation_5.html new file mode 100644 index 00000000..082a32eb --- /dev/null +++ b/docs/source/getting_started/images/1_simulation_5.html @@ -0,0 +1,71 @@ + + + +
+
+ + \ No newline at end of file diff --git a/docs/source/getting_started/index.rst b/docs/source/getting_started/index.rst index 669c8097..8ccd582f 100644 --- a/docs/source/getting_started/index.rst +++ b/docs/source/getting_started/index.rst @@ -1,7 +1,25 @@ +.. _README: https://github.com/DavAug/chi/README.md + *************** Getting started *************** .. currentmodule:: chi -To do. \ No newline at end of file +This part of the documentation gets you started using chi. Each section will +give a brief introduction into the modelling framework and show examples of +how to implement models in chi. The covered topics include **simulation** and +**inference** of e.g. ODE models, PKPD models with drug administration and +non-linear mixed effects models / hierarchical models. For installation +instructions please refer to the README_. + +.. toctree:: + :numbered: + :caption: Table of Contents + :maxdepth: 1 + + quick_overview + mechanistic_model + error_model + log_likelihood + log_posterior \ No newline at end of file diff --git a/docs/source/getting_started/log_likelihood.rst b/docs/source/getting_started/log_likelihood.rst new file mode 100644 index 00000000..eb856338 --- /dev/null +++ b/docs/source/getting_started/log_likelihood.rst @@ -0,0 +1,5 @@ +************************* +Defining a log-likelihood +************************* + +To do. \ No newline at end of file diff --git a/docs/source/getting_started/log_posterior.rst b/docs/source/getting_started/log_posterior.rst new file mode 100644 index 00000000..bbe5d599 --- /dev/null +++ b/docs/source/getting_started/log_posterior.rst @@ -0,0 +1,5 @@ +************************ +Defining a log-posterior +************************ + +To do. \ No newline at end of file diff --git a/docs/source/getting_started/mechanistic_model.rst b/docs/source/getting_started/mechanistic_model.rst new file mode 100644 index 00000000..ba6bd21c --- /dev/null +++ b/docs/source/getting_started/mechanistic_model.rst @@ -0,0 +1,5 @@ +**************************** +Defining a mechanistic model +**************************** + +To do. \ No newline at end of file diff --git a/docs/source/getting_started/quick_overview.rst b/docs/source/getting_started/quick_overview.rst new file mode 100644 index 00000000..50e09245 --- /dev/null +++ b/docs/source/getting_started/quick_overview.rst @@ -0,0 +1,276 @@ +.. _Pints: https://github.com/pints-team/pints + +.. currentmodule:: chi + +************** +Quick overview +************** + +The quick overview displays some of chi's main features to help you decide +whether chi is the software that you want to use. The example model we will be +using is a 1-compartment pharmacokinetic model (any other deterministic time +series model would have been an equally good choice, but the 1-compartment +pharmacokinetic model happens to be predefined in chi's model library; +:class:`chi.library.ModelLibrary`). The model is defined by an ordinary +differential equation for the drug amount and an algebraic equation for the +drug concentration + +.. math:: + \frac{\mathrm{d}a}{\mathrm{d}t} = -k_e a, + \quad c = \frac{a}{v}, + \quad a(t=0) = a_0, + +where :math:`a` is the drug amount in the compartment, :math:`t` is the time, +:math:`k_e` is the elimination rate of the drug, :math:`c` is the drug +concentration, :math:`v` is the volume of the compartment and :math:`a_0` is +the initial drug amount (note that during this quick overview we will omit +units for simplicity). + +Simulating a mechanistic model +****************************** + +Using the 1-compartment pharmacokinetic model from chi's model library, +simulation of the model simply involves specifying the set of model parameters +:math:`\psi = (a_0, v, k_e)` and the time points of interest + +.. literalinclude:: code/1_simulation_1.py + :lines: 2-18 + +The simulation returns a :class:`numpy.ndarray` with the simulated drug +concentrations at the specified times +:math:`[c(t=0), c(t=0.2), c(t=0.5), c(t=0.6), c(t=1)]` + +.. code-block:: console + + >>> result + array([[5. , 4.09359579, 3.03157658, 2.74324885, 1.83903834]]) + +The simulation results are of shape ``(n_outputs, n_times)`` which in this case +is ``(1, 5)`` because we simuated the drug concentration for 5 different times. +For details on how to implement your own :class:`chi.MechanisticModel` and +many other details concerning the mechanistic model in chi, we refer to the +next section :doc:`mechanistic_model`. + +Visualisation of the simulation +******************************* + +The results of the simulation can be visualised with any of the standard Python +visualisation libraries (e.g. matplotlib, seaborn, etc.). +We will use plotly and make the simulation results slightly more interesting +by administering bolus doses to the compartment in intervals of 0.5 + +.. literalinclude:: code/1_simulation_1.py + :lines: 22-46 + +.. raw:: html + :file: images/1_simulation_1.html + + +Simulation of measurements +************************** + +Measurements are in chi modelled using a :class:`chi.ErrorModel` of the form + +.. math:: + y(\psi, t) = \bar{y}(\psi, t) + \epsilon (\bar{y}, \sigma, t), + +where :math:`y` models the measurements and :math:`\epsilon` models the +residual error of the mechanistic model output with respect to the +measurements. :math:`\sigma` are the error model parameters. + +The simplest error model is a Gaussian error model +:math:`\epsilon \sim \mathcal{N}(\cdot | 0, \sigma ^2)`. So let us use a +Gaussian error model to simulate measurements of the drug concentration using +the simulation results from the previous section (we will down sample the +measurement times relative to the evaluation times of the mechanistic model for +aestetic reasons) + +.. literalinclude:: code/1_simulation_1.py + :lines: 56-90 + +.. raw:: html + :file: images/1_simulation_2.html + +For details on how to implement a :class:`chi.ErrorModel` and +many other details concerning error models in chi, we refer to the section +:doc:`error_model`. + + +Inference of model parameters +***************************** + +While the simulation of mechanistic model outputs and measurements is an +interesting feature of chi in its own right, the inference of model parameters +from real-world measurements is arguably an even more interesting feature of +chi. Here, we will use the simulated measurements to infer the model +parameters, but the simulated measurement can be straightforwardly replaced +by real-world measurements. + +Inference in chi leverages the fact that the modelled measurements :math:`y` +follow a distribution that is defined by the mechanistic model and error model. +In the case of a Gaussian error model, as in the previous example, the +distribution of the measurements is a Gaussian distribution centered at the +mechanistic model output + +.. math:: + p(y | \psi , \sigma , t) = \mathcal{N}\left( + y \, |\, \bar{y}(\psi, t), \sigma ^2 + \right) . + +This distribution can be used to define the log-likelihood of different +parameter values for the measurements +:math:`\mathcal{D}=\{ (y_1, t_1), (y_2, t_2), \ldots (y_n, t_n)\}` + +.. math:: + \log p(\mathcal{D} | \psi , \sigma) = + \sum _{j=1}^n \log p(y_j | \psi , \sigma , t_j). + +In chi the log-likelihood of model parameters can be defined using +:class:`chi.LogLikelihood`. A :class:`chi.LogLikelihood` defines the +distribution of the measurements using a :class:`MechanisticModel` and a +:class:`ErrorModel` for each mechanistic model output, and couples the +distribution to the measurements as defined above. The log-likelihood of +different parameter values can now be evaluated using the +:meth:`chi.LogLikelihood.__call__` method + +.. literalinclude:: code/1_simulation_1.py + :lines: 96-110 + +.. code-block:: console + + >>> score_1 + -86.14214936785024 + >>> score_2 + 10.324673731287309 + + +We can see that the data-generating parameter values have a larger +log-likelihood than the made-up parameter values (which should intuitively make +sense). For details on how to define a :class:`chi.LogLikelihood` and +many other details concerning log-likelihoods in chi, we refer to the section +:doc:`log_likelihood`. + +Maximum likelihood estimation +----------------------------- + +A popular approach to infer the model parameters that best describe the +measurements is to find the parameter values that maximise the log-likelihood, +a.k.a. maximum likelihood estimation + +.. math:: + \hat{\psi}, \hat{\sigma} = + \mathop{\text{arg}\,\text{max}}_{\psi, \sigma} + \log p(\mathcal{D} | \psi , \sigma). + +In chi the maximum likelihood estimates of the model parameters can be found +using any of the standard Python optimisation libraries such as scipy.optimize. +We will use Pints_ and its implementation of the Covariance Matrix +Adaption-Evolution Strategy (CMA-ES) optimisation algorithm to optimise the +log-likelihood + +.. literalinclude:: code/1_simulation_1.py + :lines: 114-119 + +.. code-block:: console + + >>> parameters_mle + array([10.26564936, 2.01524534, 1.00148417, 0.18456719]) + >>> score + 10.832126448890321 + +The optimisation algorithm finds a set of parameter values that is different, +but close to the data-generating parameters. Note, however, that the inferred +parameters have a larger log-likelihood than the data-generating parameters, +which may seem surprising at first. How can +a different set of parameters be more likely to describe the measurements than +the set of parameters that generated the measurments? +A thorough discussion of the shortcomings of maximum likelihood estimation +is beyond the scope of this documentation, but let us plot the +mechanstic model output for the inferred parameters to confirm that +inferred model is indeed marginally closer to the measurements + +.. literalinclude:: code/1_simulation_1.py + :lines: 124-132 + +.. raw:: html + :file: images/1_simulation_3.html + + +Bayesian inference +------------------ + +The reason why maximum likelihood estimation can differ from the +data-generating parameters is that a finite number of +measurements does not uniquely define the data-generating parameters of a +probabilistic model. The +maximum likelihood estimates can therefore be far away from the +data-generating parameters when the measurements are limited (in this example +the measurements estimate the data-generating parameters clearly quite well). +In other words, a finite number of measurements +leaves uncertainty about the model parameters, a.k.a. *parametric uncertainty*. +While for real-world measurements the notion of +data-generating parameters may seem alien since models only +approximate the real data-generating processes, we can generalise the notion of +data-generating parameters to the set of parameters that capture the most about +the data-generating process for a given model. Here, the maximum +likelihood estimates can analogously differ significantly from the sought after +data-generating parameters. + +In chi the uncertainty of parameter estimates can be estimated using Bayesian +inference. In Bayesian inference Bayes' rule is used to define a distribution +of likely parameter values conditional on the measurements, a.k.a. the +posterior parameter distribution + +.. math:: + \log p(\psi, \sigma | \mathcal{D}) = + \log p(\mathcal{D} | \psi , \sigma) + \log p(\psi , \sigma ) + + \text{constant}. + +Here, :math:`\log p(\psi, \sigma | \mathcal{D})` is the log-posterior, +:math:`\log p(\mathcal{D} | \psi , \sigma)` is the log-likelihood as defined +above and :math:`\log p(\psi , \sigma )` is the log-prior distribution of the +model parameters. The prior distribution of the model parameters is a +modelling choice. + +In chi the log-posterior can be defined using :class:`chi.LogPosterior` which +is instiated with a :class:`chi.LogLikelihood` and a :class:`pints.LogPrior`. +For simplicity, we will use uniform priors that constrain the parameters to +values between 0 and 100. The log-posterior can be evaluated up to the constant +term in the equation above similar to the :class:`chi.LogLikelihood` using +:meth:`chi.LogPosterior.__call__` + +.. literalinclude:: code/1_simulation_1.py + :lines: 138-149 + +.. code-block:: console + + >>> score_1 + -104.56283011180261 + >>> score_2 + -8.096007012665059 + +For details on how to define a :class:`chi.LogPosterior` and +many other details concerning log-posetriors in chi, we refer to the section +:doc:`log_posterior`. + +While the :class:`chi.LogPosterior` allows us to evaluate the log-posterior +up to constant term for different parameter values it does not yet tell us +how the posterior distribution of likely parameter values looks like. +This distribution can be inferred from :class:`chi.LogPosterior` using MCMC +sampling algorithms. In particular, we will use the implementation of Haario +and Bardenet's Adaptive Covariance Matrix Marcov Chain Monte Carlo, +:class:`pints.HaarioBardenetACMC` to infer the posterior distribution + +.. literalinclude:: code/1_simulation_1.py + :lines: 152-232 + +.. raw:: html + :file: images/1_simulation_4.html + +Some more text + +.. literalinclude:: code/1_simulation_1.py + :lines: 237-339 + +.. raw:: html + :file: images/1_simulation_5.html diff --git a/docs/source/index.rst b/docs/source/index.rst index fb65701f..2b0ef0bf 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -1,5 +1,5 @@ -.. Root of all pints docs +.. Root of all chi docs .. _GitHub: https://github.com/DavAug/chi