#!/usr/bin/env python
# -*- coding: utf-8 -*-
# This file is part of curveball.
# https://github.com/yoavram/curveball
# Licensed under the MIT license:
# http://www.opensource.org/licenses/MIT-license
# Copyright (c) 2015, Yoav Ram <yoav@yoavram.com>
from __future__ import division
from builtins import range
import warnings
import copy
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import curve_fit
import pandas as pd
import lmfit
import curveball
import seaborn as sns
sns.set_style("ticks")
def _alfa(t, q0, v):
if np.isinf(q0) or np.isinf(v):
return 1.0
return q0 / (q0 + np.exp(-v * t))
[docs]def double_baranyi_roberts_ode0(y, t, r, K, nu, q0, v):
r"""A two species Baranyi-Roberts model [1]_. The function calculates the population growth rate at given time points.
.. math::
\frac{dN_i}{dt} = r_i \alpha_i(t) N_i \Big(1 - \Big(\frac{\sum_{j}{N_j}}{K_i}\Big)^{\nu_i}\Big)
- :math:`N_i`: population size of strain *i*.
- :math:`r_i`: initial per capita growth rate of strain *i*
- :math:`K_i`: maximum population size of strain *i*
- :math:`\nu_i`: curvature of the logsitic term of strain *i*
- :math:`\alpha_i(t)= \frac{q_{0,i}}{q_{0,i} + e^{-v_i t}}`
- :math:`q_{0,i}`: initial adjustment of strain *i* to current environment
- :math:`v_i`: adjustment rate of strain *i*
Parameters
----------
y : float, float
population size
t : float
time, usually in hours
r : float, float
initial per capita growth rate
K : float, float
maximum population size (:math:`K>0`)
nu : float, float
curvature of the logsitic term (:math:`\nu>0`)
q0 : float, float
initial adjustment to current environment (:math:`0<q_0<1`)
v : float, float
adjustment rate (:math:`v>0`)
Returns
-------
float, float
population growth rate.
References
----------
.. [1] Baranyi, J., Roberts, T. A., 1994. `A dynamic approach to predicting bacterial growth in food <www.ncbi.nlm.nih.gov/pubmed/7873331>`_. Int. J. Food Microbiol.
"""
alfa = _alfa(t, q0[0], v[0]), _alfa(t, q0[1], v[1])
dydt = alfa[0] * r[0] * y[0] * (1 - ((y[0] + y[1]) / K[0])**nu[0]), alfa[1] * r[1] * y[1] * (1 - ((y[0] + y[1]) / K[1])**nu[1])
return dydt
[docs]def double_baranyi_roberts_ode1(y, t, r, K, nu, q0, v):
r"""A two species Baranyi-Roberts model. The function calculates the population growth rate at given time points.
.. math::
\frac{dN_i}{dt} = r_i \alpha_i(t) N_i \Big(1 - \Big(\sum_{j}{\frac{N_j}{K_j}}\Big)^{\nu_i}\Big)
See also
--------
curveball.competitions.double_baranyi_roberts_ode0
"""
alfa = _alfa(t, q0[0], v[0]), _alfa(t, q0[1], v[1])
dydt = alfa[0] * r[0] * y[0] * (1 - (y[0] / K[0] + y[1] / K[1])**nu[0]), alfa[1] * r[1] * y[1] * (1 - (y[0] / K[0] + y[1]/K[1])**nu[1])
return dydt
[docs]def double_baranyi_roberts_ode2(y, t, r, K, nu, q0, v):
r"""A two species Baranyi-Roberts model. The function calculates the population growth rate at given time points.
.. math::
\frac{dN_i}{dt} = r_i \alpha_i(t) N_i \Big(1 - \sum_{j}{\Big(\frac{N_j}{K_j}\Big)^{\nu_j}}\Big)
See also
--------
curveball.competitions.double_baranyi_roberts_ode0
"""
alfa = _alfa(t, q0[0], v[0]), _alfa(t, q0[1], v[1])
dydt = alfa[0] * r[0] * y[0] * (1 - (y[0] / K[0])**nu[0] - (y[1] / K[1])**nu[1]), alfa[1] * r[1] * y[1] * (1 - (y[0] / K[0])**nu[0] - (y[1] / K[1])**nu[1])
return dydt
[docs]def double_baranyi_roberts_ode3(y, t, r, K, nu, q0, v):
r"""A two species Baranyi-Roberts model. The function calculates the population growth rate at given time points.
.. math::
\frac{dN_i}{dt} = r_i \alpha_i(t) N_i \Big( 1 - \Big(\frac{\sum_j{N_j}}{\bar{K}}\Big)^{\nu_i} \Big)
See also
--------
curveball.competitions.double_baranyi_roberts_ode0
"""
alfa = _alfa(t, q0[0], v[0]), _alfa(t, q0[1], v[1])
Kmean = (K[0] * y[0] + K[1] * y[1]) / (y[0] + y[1])
dydt = alfa[0] * r[0] * y[0] * (1 - ((y[0] + y[1]) / Kmean)**nu[0]), alfa[1] * r[1] * y[1] * (1 - ((y[0] + y[1]) / Kmean)**nu[1])
return dydt
[docs]def double_baranyi_roberts_ode4(y, t, r, K, nu, q0, v):
r"""A two species Baranyi-Roberts model. The function calculates the population growth rate at given time points.
.. math::
\frac{dN_i}{dt} = r_i \alpha_i(t) N_i
See also
--------
curveball.competitions.double_baranyi_roberts_ode0
"""
alfa = _alfa(t, q0[0], v[0]), _alfa(t, q0[1], v[1])
dydt = alfa[0] * r[0] * y[0], alfa[1] * r[1] * y[1]
return dydt
[docs]def double_baranyi_roberts_ode5(y, t, r, K, nu, q0, v):
r"""A two species Baranyi-Roberts model. The function calculates the population growth rate at given time points.
.. math::
\frac{dN_i}{dt} = r_i \alpha_i(t) N_i \Big( 1 - \Big(\frac{N_i}{K_i}\Big)^{\nu_i}\Big)
See also
--------
curveball.competitions.double_baranyi_roberts_ode0
"""
alfa = _alfa(t, q0[0], v[0]), _alfa(t, q0[1], v[1])
dydt = alfa[0] * r[0] * y[0] * (1 - (y[0]/K[0])**nu[0]), alfa[1] * r[1] * y[1] * (1 - (y[1]/K[1])**nu[1])
return dydt
[docs]def double_baranyi_roberts_ode6(y, t, r, K, nu, q0, v):
r"""A two species Baranyi-Roberts model. The function calculates the population growth rate at given time points.
.. math::
\frac{dN_i}{dt} = r_i \alpha_i(t) N_i \Big(1 - \Big(\frac{\sum_j{N_j}}{K_i}\Big)^{\bar{\nu}}\Big)
See also
--------
curveball.competitions.double_baranyi_roberts_ode0
"""
alfa = _alfa(t, q0[0], v[0]), _alfa(t, q0[1], v[1])
numean = (nu[0]*y[0] + nu[1]*y[1]) / (y[0] + y[1])
dydt = alfa[0] * r[0] * y[0] * (1 - ((y[0]+y[1])/K[0])**numean), alfa[1] * r[1] * y[1] * (1 - ((y[0]+y[1])/K[1])**numean)
return dydt
[docs]def double_baranyi_roberts_ode7(y, t, r, K, nu, q0, v):
r"""A two species Baranyi-Roberts model. The function calculates the population growth rate at given time points.
.. math::
\frac{dN_1}{dt} = 0 \;\;\; \frac{dN_2}{dt} = r_2 \alpha_2(t) N_2 \Big(1 - \Big(\frac{N_2}{K_2}\Big)^{\nu_1}\Big)
See also
--------
curveball.competitions.double_baranyi_roberts_ode0
"""
alfa = _alfa(t, q0[0], v[0]), _alfa(t, q0[1], v[1])
numean = (nu[0]*y[0] + nu[1]*y[1]) / (y[0] + y[1])
dydt = 0, alfa[1] * r[1] * y[1] * (1 - (y[1]/K[1])**nu[1])
return dydt
[docs]def double_baranyi_roberts_ode8(y, t, r, K, nu, q0, v):
r"""A two species Baranyi-Roberts model. The function calculates the population growth rate at given time points.
.. math::
\frac{dN_i}{dt} = r_i N_i \Big(1 - \Big(\frac{\sum_{j}{N_j}}{K_i}\Big)^{\nu_i}\Big)
See also
--------
curveball.competitions.double_baranyi_roberts_ode0
"""
dydt = r[0] * y[0] * (1 - ((y[0] + y[1]) / K[0])**nu[0]), r[1] * y[1] * (1 - ((y[0] + y[1]) / K[1])**nu[1])
return dydt
[docs]def double_baranyi_roberts_gimenez_delgado_ode(y, t, r, K, nu, q0, v):
r"""A two species Baranyi-Roberts model with competition model inspired by Gimenez and Delgado (2004).
The function calculates the population growth rate at given time points.
.. math::
\frac{dN_i}{dt} = r_i \alpha_i(t) N_i \prod_j{\Big(1 - \Big(\frac{N_j}{K_j}\Big)^{\nu_j}\Big)}
See also
--------
curveball.competitions.double_baranyi_roberts_ode0
"""
alfa = _alfa(t, q0[0], v[0]), _alfa(t, q0[1], v[1])
prod = (1 - (y[0]/K[0])**(nu[0])) * (1 - (y[1]/K[1])**(nu[1]))
dydt = alfa[0] * r[0] * y[0] * prod, alfa[1] * r[1] * y[1] * prod
return dydt
[docs]def compete(m1, m2, p0=(0.5, 0.5), y0=None, t=None, hours=24, num_of_points=100,
nsamples=1, lag_phase=True, ode=double_baranyi_roberts_ode1,
params1=None, params2=None, ci=95, colors=None, ax=None, PLOT=False,
sampler='covar', df1=None, df2=None):
"""Simulate competitions between two strains using growth parameters estimated
by fitting growth models to growth curves data.
Integrates a 2-dimensional ODE system given by `ode`
with parameters extracted from :py:class:`lmfit.model.ModelResult` instances `m1` and `m2`.
This implementation includes plotting (if required by `PLOT`);
resampling from the distribution of model parameters (when `nsamples` > 1);
changing the ODE system (by providing a different function in `ode`);
and more.
The function competes two strains/species;
therefore it expects two :py:class:`lmfit.model.ModelResult` objects,
two initial values in `y0`, etc.
Parameters
----------
m1, m2 : lmfit.model.ModelResult
model fitting results of growth models defined in :py:mod:`curveball.models`.
p0, y0 : tuple, optional
`p0` is the initial relative frequencies; `y0` is the initial population sizes.
if `y0` is given than ``y0[0]`` and ``y0[0]`` are the initial population size for `m1` and `m2`, respectively.
if `y0` is not given, then it will be set to the average estimated `y0` in `m1` and `m2` multiplied by `p0`.
t : numpy.ndarray or None, optional
array of time points in which to calculate the population sizes.
hours : int, optional
if `t` is not given, determines how many hours should the competition proceed, defaults to 24.
num_of_points : int, optional
if `t` is not given, determines the number of time points to use, defaults to 100.
nsamples : int, optional
how many replicates of the competition should be simulated;
if `nsamples` = 1, only one competition is simulated with the estimated parameters;
otherwise `nsamples` competitions are simulated with parameters drawn from a parameter distribution
determined by `sampler`. Defaults to 1.
lag_phase : bool, optional
if `True`, use lag phase as given by `m1` and `m2`. Otherwise, override the lag phase parameters to prevent a lag phase. Defaults to :const:`True`.
ode : func, optional
an ordinary differential systems system defined by a function that accepts ``y``, ``t``, and additional arguments, and returns the derivate of ``y`` with respect to ``t``. Defaults to :py:func:`.double_baranyi_roberts_ode0`.
params1, params2 : dict, optional
dictionaries of model parameter values; if given, overrides values from `m1` and `m2`.
ci : float, optional
confidence interval size, in (0, 100), only applicable when `PLOT` is :const:`True`, defaults to 95%.
colors : sequence of str, optional
if `PLOT` is :const:`True`, this sets the colors of the drawn lines. `colors[0]` will be used for `m1`; `colors[1]` for `m2`. If not provided, defaults to the current pallete.
ax : matplotlib.axes.Axes, optional
if `PLOT` is :const:`True`, an axes to plot into; if not provided, a new one is created.
PLOT : bool, optional
if :const:`True`, the function will plot the curves of *y* as a function of *t*. Defaults to :const:`False`.
sampler : str, optional
if ``covar``, the parameters will be sampled using the covariance matrix (:py:func:`curveball.models.sample_params`);
if ``bootstrap`` the parameters will be sampled by resampling the growth curves (:py:func:`curveball.models.bootstrap_params`).
df1, df2 : pandas.DataFrame, optional
the data frames used to fit `m1` and `m2`, only used when `sampler` is ``bootstrap``.
Returns
-------
t : numpy.ndarray
1d (or 2d, if `nsamples`>1) array of time points, in hours.
y: numpy.ndarray
2d (or 3d, if `nsamples`>1) array of strain frequencies.
First axis is time, second axis is strain, third axis (if applicable) is sample.
fig : matplotlib.figure.Figure
figure object.
ax : numpy.ndarray
array of :py:class:`matplotlib.axes.Axes` objects.
Raises
------
TypeError
if `m1` or `m2` are not :py:class:`lmfit.model.ModelResult`.
AssertionError
if an intermediate calculation produced an invalid result (not guaranteed).
Examples
--------
>>> import pandas as pd
>>> import curveball
>>> plate = pd.read_csv('plate_templates/G-RG-R.csv')
>>> df = curveball.ioutils.read_tecan_xlsx('data/Tecan_280715.xlsx', label='OD', plate=plate)
>>> green = curveball.models.fit_model(df[df.Strain == 'G'], PLOT=False, PRINT=False)[0]
>>> red = curveball.models.fit_model(df[df.Strain == 'R'], PLOT=False, PRINT=False)[0]
>>> t, y = curveball.competitions.compete(green, red, PLOT=False)
Notes
-----
To debug, uncomment lines to return the ``infodict`` that is returned from :py:func:`scipy.integrate.odeint` is run with ``full_output=1``.
"""
if not isinstance(m1, lmfit.model.ModelResult):
raise TypeError("m1 must be %s, instead it is %s", lmfit.model.ModelResult, type(m1))
if not isinstance(m2, lmfit.model.ModelResult):
raise TypeError("m2 must be %s, instead it is %s", lmfit.model.ModelResult, type(m2))
if t is None:
t = np.linspace(0, hours, num_of_points)
if nsamples > 1:
# draw random params
sampler = sampler.lower()
if sampler == 'covar':
m1_samples = curveball.models.sample_params(m1, nsamples, params=params1)
m2_samples = curveball.models.sample_params(m2, nsamples, params=params2)
elif sampler == 'bootstrap':
if params1 or params2:
warnings.warn("Bootstrap sampling doesn't support params1 and params2 arguments")
if df1 is None or df2 is None:
raise ValueError("Bootstrap sampling requires kwargs df1 and df2")
m1_fixed = {pname for pname,p in m1.params.items() if not p.vary}
m2_fixed = {pname for pname,p in m2.params.items() if not p.vary}
# FIXME bootstrap_params has changed
m1_samples = curveball.models.bootstrap_params(df1, m1, nsamples,
fit_kws={'param_fix': m1_fixed})
m2_samples = curveball.models.bootstrap_params(df2, m2, nsamples,
fit_kws={'param_fix': m2_fixed})
else:
raise ValueError("Unknow sampler method: {0}".format(sampler))
min_nsamples = min(len(m1_samples), len(m2_samples))
if nsamples > min_nsamples:
warnings.warn("{0} resamples lost".format(nsamples - min_nsamples))
nsamples = min_nsamples
else:
nsamples = 1
# override model result params with arguments params1 and params2
if params1:
_params = copy.copy(m1.best_values)
_params.update(params1)
params1 = _params
else:
params1 = m1.best_values
if params2:
_params = copy.copy(m2.best_values)
_params.update(params2)
params2 = _params
else:
params2 = m2.best_values
# param samples contain the model fit estimated params
m1_samples = pd.DataFrame([params1])
m2_samples = pd.DataFrame([params2])
assert len(m1_samples) == len(m2_samples)
y = np.empty((len(t), 2, nsamples))
#infodict = [None]*nsamples # DEBUG
# simulate the ode for each param sample
for i in range(nsamples):
if y0 is None:
p0 = p0[0] / (p0[0] + p0[1]), p0[1] / (p0[0] + p0[1])
y0 = m1_samples.iloc[i]['y0'] * p0[0], m2_samples.iloc[i]['y0'] * p0[1]
r = m1_samples.iloc[i]['r'], m2_samples.iloc[i]['r']
K = m1_samples.iloc[i]['K'], m2_samples.iloc[i]['K']
nu = m1_samples.iloc[i].get('nu', 1.0), m2_samples.iloc[i].get('nu', 1.0)
if lag_phase:
q0 = m1_samples.iloc[i].get('q0', np.inf), m2_samples.iloc[i].get('q0', np.inf)
v = m1_samples.iloc[i].get('v', r[0]), m2_samples.iloc[i].get('v', r[1])
else:
q0 = np.inf, np.inf
v = np.inf, np.inf
args = (r, K, nu, q0, v)
y[:,:,i] = odeint(ode, y0=y0, t=t, args=args)
# DEBUG
#_y_,info = odeint(double_baranyi_roberts_ode0, y0, t, args=args, full_output=1)
#info['args'] = (y0,) + args
#infodict[i] = info
#if info['message'] == 'Integration successful.':
# y[:,:,i] = _y_
if PLOT:
if ax is None:
fig,ax = plt.subplots(1, 1)
else:
fig = ax.get_figure()
df = pd.DataFrame()
for i in range(y.shape[1]):
_df = pd.DataFrame(y[:,i,:])
_df['Time'] = t
_df = pd.melt(_df, id_vars='Time', var_name='Replicate', value_name='y')
_df['Strain'] = i
df = pd.concat((df, _df))
if colors is not None:
colors = {i:c for i,c in enumerate(colors)}
sns.lineplot(data=df, x='Time', hue='Strain', y='y',
ci=ci, palette=colors, ax=ax)
ax.set_xlabel('Time (hour)')
ax.set_ylabel('OD')
sns.despine()
return t,y,fig,ax
return t,y#,infodict
[docs]def selection_coefs_ts(t, y, ax=None, PLOT=False):
r"""Calculate selection coefficient according to the following formula[2]_,
where :math:`A(t), B(t)` are population densities of assay strain *A* and reference strain *B* at time *t*:
.. math::
s = \frac{d}{dt} \log{\frac{A(t)}{B(t)}}
Parameters
----------
t : numpy.ndarray
array of time points, as produced by :py:func:`compete`
y : numpy.ndarray
array of population densities, as produced by :py:func:`compete`, where the first axis is time and the second axis is strain.
ax : matplotlib.axes.Axes, optional
if `PLOT` is :const:`True`, an axes to plot into; if not provided, a new one is created.
PLOT : bool, optional
if :const:`True`, the function will plot the curve of *s* as a function of *t*.
Returns
-------
svals : numpy.ndarray
the selection coefficients of the assay strain relative to the reference strain over time.
fig : matplotlib.figure.Figure
figure object.
ax : numpy.ndarray of matplotlib.axes.Axes
array of axes objects.
Notes
-----
This formula assumes that the frequencies of the strains follow a logistic curve.
Lag phases, interactions, etc. may cause this formula to become irrelevant.
References
----------
.. [12] Chevin, L-M. 2011. `On Measuring Selection in Experimental Evolution <http://dx.doi.org/10.1098/rsbl.2010.0580>`_. Biology Letters.
"""
dt = np.ediff1d(t, to_end=1)
dt[-1] = dt[-2]
svals = np.gradient(np.log(y[:,0] / y[:,1]), dt)
svals[np.isinf(svals)] = svals[np.isfinite(svals)].max()
if PLOT:
if ax is None:
fig,ax = plt.subplots(1,1)
else:
fig = ax.get_figure()
ax.plot(t, svals)
ax.set_ylabel('Selection coefficient')
ax.set_xlabel('Time (hour)')
sns.despine()
return svals, fig, ax
return svals
[docs]def fitness_LTEE(y, ref_strain=0, assay_strain=1, t0=0, t1=-1, ci=0):
r"""Calculate relative fitness according to the definition used in the *Long Term Evolutionary Experiment* (LTEE) [3]_,
where :math:`A(t), B(t)` are population densities of assay strain *A* and reference strain *B* at time *t*:
.. math::
\omega = \frac{\log{(A(t) / A(0))}}{\log{(B(t) / B(0))}}
If `ci` > 0, treats the third axis of `y` as replicates (from a parameteric boostrap) to calculate the confidence interval on the fitness.
Parameters
----------
y : numpy.ndarray
array of population densities, as produced by :py:func:`compete`, where the first axis is time and the second axis is strain. A third axis is also applicable if `ci`>0.
ref_strain : int, optional
the index of the reference strain within `y`. This strain's fitness is set to 1 by definition. Defaults to 0 (first).
assay_strain : int, optional
the index of the assay strain within `y`. The result will be the fitness of this strain relative to the fitness of the reference strain. Defaults to 1 (second).
t0 : int, optional
the index of the time point from which to start the calculation of the relative fitness, defaults to 0 (first).
t1 : int, optional
the index of the time point at which to end the calculation of the relative fitness, defaults to -1 (last).
ci : float between 0 and 1, optional
if not zero, a confidence interval will be calculated using the third axis of `y` as replicates.
Returns
-------
w : float
the fitness of the assay strain relative to the reference strain.
low : float
if `ci` > 0 and `y.ndim` = 3, this is the low limit of the confidence interval of the fitness
high : float
if `ci` > 0 and `y.ndim` = 3, this is the higher limit of the confidence interval of the fitness
Raises
------
ValueError
if confidence interval is requested (`ci` > 0) but there are no replicates (`y.ndim` != 3).
Notes
-----
The result may depend on the choice of `t0` and `t1` as well as the strain designations (`ref_strain` and `assay_strain`).
References
----------
.. [3] Wiser, M. J. & Lenski, R. E. 2015 `A Comparison of Methods to Measure Fitness in Escherichia coli <http://dx.plos.org/10.1371/journal.pone.0126210>`_. PLoS One.
"""
if ci == 0:
y = y.reshape(y.shape + (1,))
elif y.ndim != 3:
raise ValueError()
w = np.zeros(y.shape[2])
for i in range(y.shape[2]):
At0, Bt0 = y[t0, assay_strain, i], y[t0, ref_strain, i]
At1, Bt1 = y[t1, assay_strain, i], y[t1, ref_strain, i]
w[i] = np.log(At1 / At0) / np.log(Bt1 / Bt0)
if ci == 0:
return w.mean()
else:
margin = (1 - ci) * 50
return w.mean(), np.percentile(w, margin), np.percentile(w, ci * 100 + margin)
[docs]def baranyi_roberts_gd(y, t, *args):
r"""
Fujikawa, Hiroshi, and Mohammad Z Sakha. 2014. Prediction of Microbial Growth in Mixed Culture with a Competition Model. Biocontrol Science 19 (2): 89-92.
"""
y1, y2 = y
K, r, nu, q0, v, a = args
r1, r2 = r
K1, K2 = K
nu1, nu2 = nu
q01, q02 = q0
v1, v2 = v
a1, a2 = a
alfa1 = _alfa(t, q01, v1)
alfa2 = _alfa(t, q02, v2)
dy1dt = r1 * alfa1 * y1 * (1 - (y1/K1)**nu1) * (1 - (y2/K2)**(nu2*a2))
dy2dt = r2 * alfa2 * y2 * (1 - (y1/K1)**(nu1*a1)) * (1 - (y2/K2)**nu2)
return [dy1dt, dy2dt]
[docs]def baranyi_roberts_lv(y, t, *args):
r"""
Fujikawa, Hiroshi, and Mohammad Z Sakha. 2014. Prediction of Microbial Growth in Mixed Culture with a Competition Model. Biocontrol Science 19 (2): 89-92.
"""
y1, y2 = y
K, r, nu, q0, v, a = args
r1, r2 = r
K1, K2 = K
nu1, nu2 = nu
q01, q02 = q0
v1, v2 = v
a1, a2 = a
alfa1 = _alfa(t, q01, v1)
alfa2 = _alfa(t, q02, v2)
Kmax = max(K1**nu1, K2**nu2)
dy1dt = r1 * alfa1 * y1 * (1 - (y1/K1)**nu1) * (1 - (y1**a1 + y2**a2)/Kmax)
dy2dt = r2 * alfa2 * y2 * (1 - (y2/K2)**nu2) * (1 - (y1**a1 + y2**a2)/Kmax)
return [dy1dt, dy2dt]
def baranyi_roberts_yr2(y, t, *args):
y1, y2 = y
K, r, nu, q0, v, a = args
r1, r2 = r
K1, K2 = K
nu1, nu2 = nu
q01, q02 = q0
v1, v2 = v
a1, a2 = a
alfa1 = _alfa(t, q01, v1)
alfa2 = _alfa(t, q02, v2)
dy1dt = r1 * alfa1 * y1 * (1 - (y1/K1)**nu1 - (a2 * y2/K2)**nu2)
dy2dt = r2 * alfa2 * y2 * (1 - (a1 * y1/K1)**nu1 - (y2/K2)**nu2)
return [dy1dt, dy2dt]
def baranyi_roberts_yr(y, t, *args):
y1, y2 = y
K, r, nu, q0, v, a = args
r1, r2 = r
K1, K2 = K
nu1, nu2 = nu
q01, q02 = q0
v1, v2 = v
a1, a2 = a
alfa1 = _alfa(t, q01, v1)
alfa2 = _alfa(t, q02, v2)
dy1dt = r1 * alfa1 * y1 * (1 - (y1**nu1) / (K1**nu1) - a2 * (y2**nu2) / (K1**nu1))
dy2dt = r2 * alfa2 * y2 * (1 - a1 * (y1**nu1) / (K2**nu2) - (y2**nu2) / (K2**nu2))
return [dy1dt, dy2dt]
# obsolete - just call baranyi_roberts_yr(y, t, *args, (1, 1))
# def baranyi_roberts_yr_a1(y, t, *args):
# return baranyi_roberts_yr(y, t, *args, (1, 1))
def fit_and_compete(m1, m2, df_mixed, y0=None, aguess=(1, 1), fixed=False,
ode=baranyi_roberts_yr, num_of_points=100, method='nelder',
value_var = 'OD', time_var = 'Time',
PLOT=False, colors=sns.color_palette('Set1', 3)):
best_values1 = m1.best_values if hasattr(m1, 'best_values') else m1
best_values2 = m2.best_values if hasattr(m2, 'best_values') else m2
K = best_values1['K'], best_values2['K']
r = best_values1['r'], best_values2['r']
nu = best_values1.get('nu',1), best_values2.get('nu',1)
q0 = best_values1.get('q0', np.inf), best_values2.get('q0', np.inf)
v = best_values1.get('v', r[0]), best_values2.get('v', r[1])
if y0 is None:
y0 = best_values1['y0']/2, best_values2['y0']/2
y_mixed = df_mixed.groupby(time_var)[value_var].mean().values
y_mixed_std = df_mixed.groupby(time_var)[value_var].std().values
t_mixed = np.unique(df_mixed[time_var])
t = np.linspace(0, t_mixed.max(), num_of_points)
if fixed:
a = aguess
a1, a2 = a
else:
def mixed_model(t, a1, a2):
y = odeint(ode, y0, t, args=(K, r, nu, q0, v, (a1, a2)))
return y.sum(axis=1)
model = lmfit.Model(mixed_model)
params = model.make_params(a1=aguess[0], a2=aguess[1])
params['a1'].set(min=1e-1, max=10)
params['a2'].set(min=1e-1, max=10)
result = model.fit(data=y_mixed, t=t_mixed, params=params, weights=1.0 / y_mixed_std,
method=method, nan_policy='propagate')
a = result.best_values['a1'], result.best_values['a2']
a1, a2 = a
y = odeint(ode, y0, t, args=(K, r, nu, q0, v, a))
if PLOT:
ysum = y.sum(axis=1)
p1 = y[:, 0] / ysum
p2 = y[:, 1] / ysum
MRSE = ((odeint(ode, y0, t_mixed, args=(K, r, nu, q0, v, a)).sum(axis=1) - y_mixed)**2).mean()
w, h = plt.rcParams['figure.figsize']
fig, ax = plt.subplots(1, 2, figsize=(w * 2, h), sharex=True)
ax[0].plot(t, y[:, 0], color=colors[0])
ax[0].plot(t, y[:, 1], color=colors[1])
ax[0].plot(t, ysum, color=colors[2])
if not fixed:
ax[0].plot(t_mixed, result.init_fit, ls='-.', color=colors[2])
ax[0].errorbar(t_mixed, y_mixed, y_mixed_std, fmt='o', color=colors[2])
ax[0].set(xlabel='Time (hour)', ylabel='OD',
xlim=(-0.5, t.max() + 0.5),
title="MRSE: {:.2g}".format(MRSE))
ax[1].plot(t, p1, color=colors[0])
ax[1].plot(t, p2, color=colors[1])
ax[1].set(xlabel='Time (hour)', ylabel='Frequency',
xlim=(-0.5, t.max() + 0.5), ylim=(0, 1),
title="a1={:.2g}, a2={:.2g}".format(*a))
sns.despine()
fig.tight_layout()
return t, y, a, fig, ax
return t, y, a
def fit_and_compete_ci(param_samples1, param_samples2, df_mixed, y0=None, ci=0.95,
colors=('g', 'r'), line_kws=None, ci_kws=None, ax=None, PLOT=False):
if line_kws is None:
line_kws = dict()
if ci_kws is None:
ci_kws = dict(color='gray', alpha=0.5)
nsamples = param_samples1.shape[0]
assert param_samples1.shape[0] == param_samples2.shape[0], "Parameters samples should have the same length"
t = [None] * nsamples
y = [None] * nsamples
a = [None] * nsamples
for i in range(nsamples):
sample1 = param_samples1.iloc[i, :]
sample2 = param_samples2.iloc[i, :]
t[i], y[i], a[i] = curveball.competitions.fit_and_compete(sample1, sample2, df_mixed, y0=y0)
t = np.array(t)
y = np.array(y)
a = np.array(a)
ysum = y.sum(axis=2)
f1 = y[:, :, 0] / ysum
f2 = y[:, :, 1] / ysum
f = np.array((f1, f2)).T
margin = (1.0 - ci) * 50.0
low_f = np.percentile(f, margin, axis=1)
high_f = np.percentile(f, ci * 100.0 + margin, axis=1)
avg_f = f.mean(axis=1)
low_a = np.percentile(a, margin, axis=0)
avg_a = a.mean(axis=0)
high_a = np.percentile(a, ci * 100.0 + margin, axis=0)
if PLOT:
if ax is None:
fig, ax = plt.subplots()
else:
fig = ax.figure
ax.fill_between(t[0, :], low_f[:, 0], high_f[:, 0], **ci_kws)
ax.plot(t[0,:], avg_f[:, 0], color=colors[0], **line_kws)
ax.fill_between(t[0, :], low_f[:, 1], high_f[:, 1], **ci_kws)
ax.plot(t[0,:], avg_f[:, 1], color=colors[1], **line_kws)
ax.set(xlabel='Time', ylabel='Frequency')
fig.tight_layout()
sns.despine()
return low_a, avg_a, high_a, low_f, avg_f, high_f, fig, ax
return low_a, avg_a, high_a, low_f, avg_f, high_f