In [1]:
import os, tempfile
import logging
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
In [2]:
from ionchannelABC import theoretical_population_size
from ionchannelABC import IonChannelDistance, EfficientMultivariateNormalTransition, IonChannelAcceptor
from ionchannelABC.experiment import setup
from ionchannelABC.visualization import plot_sim_results, plot_kde_matrix_custom
import myokit
In [3]:
from pyabc import Distribution, RV, History, ABCSMC
from pyabc.epsilon import MedianEpsilon
from pyabc.sampler import MulticoreEvalParallelSampler, SingleCoreSampler
from pyabc.populationstrategy import ConstantPopulationSize
In [4]:
from experiments.isus_wang import wang_act_and_kin
from experiments.isus_courtemanche import courtemanche_deact
from experiments.isus_firek import (firek_inact)
from experiments.isus_nygren import (nygren_inact_kin,
nygren_rec)
In [5]:
modelfile = 'models/standardised_isus.mmt'
Plot steady-state and tau functions
In [7]:
from ionchannelABC.visualization import plot_variables
In [8]:
sns.set_context('poster')
V = np.arange(-80, 40, 0.01)
sta_par_map = {'ri': 'isus.r_ss',
'si': 'isus.s_ss',
'rt': 'isus.tau_r',
'st': 'isus.tau_s'}
f, ax = plot_variables(V, sta_par_map, modelfile, figshape=(2,2))
Combine model and experiments to produce:
In [9]:
observations, model, summary_statistics = setup(modelfile,
wang_act_and_kin,
courtemanche_deact,
firek_inact,
nygren_inact_kin,
nygren_rec)
In [10]:
assert len(observations)==len(summary_statistics(model({})))
In [11]:
g = plot_sim_results(modelfile,
wang_act_and_kin,
courtemanche_deact,
firek_inact,
nygren_inact_kin,
nygren_rec)
In [12]:
limits = {'log_isus.p_1': (-7, 3),
'isus.p_2': (1e-7, 0.4),
'log_isus.p_3': (-7, 3),
'isus.p_4': (1e-7, 0.4),
'log_isus.p_5': (-7, 3),
'isus.p_6': (1e-7, 0.4),
'log_isus.p_7': (-7, 3),
'isus.p_8': (1e-7, 0.4)}
prior = Distribution(**{key: RV("uniform", a, b - a)
for key, (a,b) in limits.items()})
In [13]:
# Test this works correctly with set-up functions
assert len(observations) == len(summary_statistics(model(prior.rvs())))
In [16]:
db_path = ("sqlite:///" + os.path.join(tempfile.gettempdir(), "standardised_isus.db"))
In [17]:
logging.basicConfig()
abc_logger = logging.getLogger('ABC')
abc_logger.setLevel(logging.DEBUG)
eps_logger = logging.getLogger('Epsilon')
eps_logger.setLevel(logging.DEBUG)
In [18]:
pop_size = theoretical_population_size(2, len(limits))
print("Theoretical minimum population size is {} particles".format(pop_size))
In [19]:
abc = ABCSMC(models=model,
parameter_priors=prior,
distance_function=IonChannelDistance(
exp_id=list(observations.exp_id),
variance=list(observations.variance),
delta=0.05),
population_size=ConstantPopulationSize(1000),
summary_statistics=summary_statistics,
transitions=EfficientMultivariateNormalTransition(),
eps=MedianEpsilon(initial_epsilon=100),
sampler=MulticoreEvalParallelSampler(n_procs=16),
acceptor=IonChannelAcceptor())
In [20]:
obs = observations.to_dict()['y']
obs = {str(k): v for k, v in obs.items()}
In [21]:
abc_id = abc.new(db_path, obs)
In [ ]:
history = abc.run(minimum_epsilon=0., max_nr_populations=100, min_acceptance_rate=0.01)
In [ ]:
history = abc.run(minimum_epsilon=0., max_nr_populations=100, min_acceptance_rate=0.01)
In [14]:
history = History('sqlite:///results/standardised/isus/standardised_isus.db')
In [15]:
df, w = history.get_distribution(m=0)
In [16]:
df.describe()
Out[16]:
In [17]:
sns.set_context('poster')
mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14
g = plot_sim_results(modelfile,
wang_act_and_kin,
courtemanche_deact,
firek_inact,
nygren_inact_kin,
nygren_rec,
df=df, w=w)
plt.tight_layout()
In [18]:
import pandas as pd
N = 100
sta_par_samples = df.sample(n=N, weights=w, replace=True)
sta_par_samples = sta_par_samples.set_index([pd.Index(range(N))])
sta_par_samples = sta_par_samples.to_dict(orient='records')
In [19]:
sns.set_context('poster')
mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14
f, ax = plot_variables(V, sta_par_map,
modelfile,
[sta_par_samples],
figshape=(2,2))
plt.tight_layout()
In [20]:
m,_,_ = myokit.load(modelfile)
In [21]:
sns.set_context('paper')
g = plot_kde_matrix_custom(df, w, limits=limits)
plt.tight_layout()
In [ ]: