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.ical_mewes import mewes_act
from experiments.ical_li import (li_inact_1000,
li_inact_kin_80,
li_recov)
Load the myokit modelfile for this channel.
In [5]:
modelfile = 'models/nygren_ical.mmt'
In [6]:
from ionchannelABC.visualization import plot_variables
In [7]:
sns.set_context('poster')
V = np.arange(-140, 50, 0.01)
nyg_par_map = {'di': 'ical.d_inf',
'f1i': 'ical.f_inf',
'f2i': 'ical.f_inf',
'dt': 'ical.tau_d',
'f1t': 'ical.tau_f_1',
'f2t': 'ical.tau_f_2'}
f, ax = plot_variables(V, nyg_par_map, 'models/nygren_ical.mmt', figshape=(3,2))
Combine model and experiments to produce:
In [8]:
observations, model, summary_statistics = setup(modelfile,
mewes_act)
assert len(observations)==len(summary_statistics(model({})))
In [9]:
g = plot_sim_results(modelfile,
mewes_act)
In [10]:
limits = {'ical.p1': (-50, 50),
'ical.p2': (0, 20),
'log_ical.p3': (-7, 3),
'ical.p4': (-50, 50),
'ical.p5': (0, 20),
'log_ical.p6': (-7, 3)}
prior = Distribution(**{key: RV("uniform", a, b - a)
for key, (a,b) in limits.items()})
In [11]:
# Test this works correctly with set-up functions
assert len(observations) == len(summary_statistics(model(prior.rvs())))
In [12]:
db_path = ("sqlite:///" + os.path.join(tempfile.gettempdir(), "nygren_ical_dgate_original.db"))
In [13]:
logging.basicConfig()
abc_logger = logging.getLogger('ABC')
abc_logger.setLevel(logging.DEBUG)
eps_logger = logging.getLogger('Epsilon')
eps_logger.setLevel(logging.DEBUG)
Test theoretical number of particles for approximately 2 particles per dimension in the initial sampling of the parameter hyperspace.
In [14]:
pop_size = theoretical_population_size(2, len(limits))
print("Theoretical minimum population size is {} particles".format(pop_size))
Initialise ABCSMC (see pyABC documentation for further details).
IonChannelDistance calculates the weighting applied to each datapoint based on the experimental variance.
In [16]:
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(500),
summary_statistics=summary_statistics,
transitions=EfficientMultivariateNormalTransition(),
eps=MedianEpsilon(),#initial_epsilon=100),
sampler=MulticoreEvalParallelSampler(n_procs=16),
acceptor=IonChannelAcceptor())
In [17]:
# Convert observations to dictionary format for calibration
obs = observations.to_dict()['y']
obs = {str(k): v for k, v in obs.items()}
In [18]:
abc_id = abc.new(db_path, obs)
Run calibration with stopping criterion of particle 1\% acceptance rate.
In [19]:
history = abc.run(minimum_epsilon=0., max_nr_populations=100, min_acceptance_rate=0.01)
In [17]:
history = History(db_path)
In [11]:
history.all_runs()
Out[11]:
In [8]:
df, w = history.get_distribution()
In [9]:
df.describe()
Out[9]:
In [23]:
sns.set_context('poster')
mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14
g = plot_sim_results(modelfile,
mewes_act,
df=df, w=w)
plt.tight_layout()
In [24]:
import pandas as pd
N = 100
nyg_par_samples = df.sample(n=N, weights=w, replace=True)
nyg_par_samples = nyg_par_samples.set_index([pd.Index(range(N))])
nyg_par_samples = nyg_par_samples.to_dict(orient='records')
In [25]:
sns.set_context('talk')
mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14
f, ax = plot_variables(V, nyg_par_map,
'models/nygren_ical.mmt',
[nyg_par_samples],
figshape=(3,2))
In [26]:
from ionchannelABC.visualization import plot_kde_matrix_custom
import myokit
import numpy as np
In [27]:
m,_,_ = myokit.load(modelfile)
In [28]:
originals = {}
for name in limits.keys():
if name.startswith("log"):
name_ = name[4:]
else:
name_ = name
val = m.value(name_)
if name.startswith("log"):
val_ = np.log10(val)
else:
val_ = val
originals[name] = val_
In [29]:
sns.set_context('notebook')
g = plot_kde_matrix_custom(df, w, limits=limits, refval=originals)
In [19]:
observations, model, summary_statistics = setup(modelfile,
li_inact_1000,
li_inact_kin_80,
li_recov)
In [20]:
assert len(observations)==len(summary_statistics(model({})))
In [21]:
g = plot_sim_results(modelfile,
li_inact_1000,
li_inact_kin_80,
li_recov)
In [22]:
limits = {'ical.q1': (0, 100),
'ical.q2': (0, 50),
'log_ical.q3': (-7, 3),
'ical.q4': (0, 100),
'ical.q5': (0, 50),
'log_ical.q6': (-7, 3),
'log_ical.r1': (-7, 3),
'ical.r2': (0, 100),
'ical.r3': (0, 100),
'log_ical.r4': (-7, 3)}
prior = Distribution(**{key: RV("uniform", a, b - a)
for key, (a,b) in limits.items()})
In [23]:
# Test this works correctly with set-up functions
assert len(observations) == len(summary_statistics(model(prior.rvs())))
In [24]:
db_path = ("sqlite:///" + os.path.join(tempfile.gettempdir(), "nygren_ical_fgates_original.db"))
In [36]:
logging.basicConfig()
abc_logger = logging.getLogger('ABC')
abc_logger.setLevel(logging.DEBUG)
eps_logger = logging.getLogger('Epsilon')
eps_logger.setLevel(logging.DEBUG)
In [37]:
pop_size = theoretical_population_size(2, len(limits))
print("Theoretical minimum population size is {} particles".format(pop_size))
In [38]:
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(2000),
summary_statistics=summary_statistics,
transitions=EfficientMultivariateNormalTransition(),
eps=MedianEpsilon(initial_epsilon=100),
sampler=MulticoreEvalParallelSampler(n_procs=8),
acceptor=IonChannelAcceptor())
In [39]:
obs = observations.to_dict()['y']
obs = {str(k): v for k, v in obs.items()}
In [40]:
abc_id = abc.new(db_path, obs)
In [41]:
history = abc.run(minimum_epsilon=0., max_nr_populations=10, min_acceptance_rate=0.01)
In [42]:
history = abc.run(minimum_epsilon=0., max_nr_populations=100, min_acceptance_rate=0.01)
In [25]:
history = History(db_path)
In [26]:
df, w = history.get_distribution(m=0)
In [27]:
df.describe()
Out[27]:
In [28]:
sns.set_context('poster')
mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14
g = plot_sim_results(modelfile,
li_inact_1000,
li_inact_kin_80,
li_recov,
df=df, w=w)
plt.tight_layout()
In [47]:
from ionchannelABC.visualization import plot_kde_matrix_custom
import myokit
import numpy as np
In [48]:
m,_,_ = myokit.load(modelfile)
In [49]:
originals = {}
for name in limits.keys():
if name.startswith("log"):
name_ = name[4:]
else:
name_ = name
val = m.value(name_)
if name.startswith("log"):
val_ = np.log10(val)
else:
val_ = val
originals[name] = val_
In [50]:
sns.set_context('paper')
g = plot_kde_matrix_custom(df, w, limits=limits, refval=originals)
plt.tight_layout()
In [31]:
import pandas as pd
N = 100
nyg_par_samples = df.sample(n=N, weights=w, replace=True)
nyg_par_samples = nyg_par_samples.set_index([pd.Index(range(N))])
nyg_par_samples = nyg_par_samples.to_dict(orient='records')
In [32]:
sns.set_context('poster')
mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14
f, ax = plot_variables(V, nyg_par_map,
'models/nygren_ical.mmt',
[nyg_par_samples],
figshape=(3,2))
plt.tight_layout()
In [ ]: