In [1]:
# PyABC imports
from pyabc import (ABCSMC, Distribution, RV,
History, MedianEpsilon)
from pyabc.populationstrategy import AdaptivePopulationSize
from pyabc.epsilon import MedianEpsilon
from pyabc.sampler import MulticoreEvalParallelSampler
In [2]:
# Custom imports
from ionchannelABC import (ion_channel_sum_stats_calculator,
IonChannelAcceptor,
IonChannelDistance,
EfficientMultivariateNormalTransition,
calculate_parameter_sensitivity,
plot_parameter_sensitivity,
plot_regression_fit)
In [3]:
# Other necessary imports
import numpy as np
import subprocess
import pandas as pd
import io
import os
import tempfile
In [4]:
# Plotting imports
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
%config Inline.Backend.figure_format = 'retina'
In [5]:
from channels.icat_generic import icat as model
model.sample({})
Out[5]:
In [6]:
measurements = model.get_experiment_data()
obs = measurements.to_dict()['y']
exp = measurements.to_dict()['exp']
errs = measurements.to_dict()['errs']
In [7]:
limits = dict(g_CaT=(0, 2),
v_offset=(0, 500),
Vhalf_b=(-100, 100),
k_b=(0, 10),
c_bb=(0, 10),
c_ab=(0, 100),
sigma_b=(0, 20),
Vmax_b=(-100, 100),
Vhalf_g=(-100, 100),
k_g=(-10, 0),
c_bg=(0, 50),
c_ag=(0, 200),
sigma_g=(0, 100),
Vmax_g=(-100, 100))
prior = Distribution(**{key: RV("uniform", a, b - a)
for key, (a,b) in limits.items()})
In [8]:
distance_fn=IonChannelDistance(
obs=obs,
exp_map=exp,
err_bars=errs,
err_th=0.1)
In [11]:
parameters = ['icat.'+k for k in limits.keys()]
In [12]:
fitted, regression_fit, r2 = calculate_parameter_sensitivity(
model,
parameters,
distance_fn,
sigma=0.1,
n_samples=1000)
In [13]:
sns.set_context('talk')
grid = plot_parameter_sensitivity(fitted, plot_cutoff=0.05)
In [14]:
grid.savefig('results/icat/sensitivity.pdf')
In [15]:
grid2 = plot_regression_fit(regression_fit, r2)
In [16]:
grid2.savefig('results/icat/sensitivity_fit.pdf')
In [17]:
from ionchannelABC import plot_distance_weights
grid = plot_distance_weights(model, distance_fn)
In [18]:
grid.savefig('results/icat/dist_weights.pdf')
In [9]:
db_path = ("sqlite:///" +
os.path.join(tempfile.gettempdir(), "hl-1_icat5000.db"))
print(db_path)
In [10]:
# Let's log all the sh!t
import logging
logging.basicConfig()
abc_logger = logging.getLogger('ABC')
abc_logger.setLevel(logging.DEBUG)
eps_logger = logging.getLogger('Epsilon')
eps_logger.setLevel(logging.DEBUG)
cv_logger = logging.getLogger('CV Estimation')
cv_logger.setLevel(logging.DEBUG)
In [11]:
from pyabc.populationstrategy import ConstantPopulationSize
In [13]:
abc = ABCSMC(models=model,
parameter_priors=prior,
distance_function=IonChannelDistance(
obs=obs,
exp_map=exp,
err_bars=errs,
err_th=0.1),
population_size=ConstantPopulationSize(5000),
#population_size=AdaptivePopulationSize(
# start_nr_particles=5000,
# mean_cv=0.2,
# max_population_size=5000,
# min_population_size=100),
summary_statistics=ion_channel_sum_stats_calculator,
transitions=EfficientMultivariateNormalTransition(),
eps=MedianEpsilon(),
sampler=MulticoreEvalParallelSampler(n_procs=6),
acceptor=IonChannelAcceptor())
In [14]:
abc_id = abc.new(db_path, obs)
In [ ]:
history = abc.run(minimum_epsilon=0.1, max_nr_populations=100, min_acceptance_rate=0.01)
In [32]:
history = [History('sqlite:////storage/hhecm/cellrotor/chouston/results/convergence/icat-generic/hl-1_icat50.db'),
History('sqlite:////storage/hhecm/cellrotor/chouston/results/convergence/icat-generic/hl-1_icat100.db'),
History('sqlite:////storage/hhecm/cellrotor/chouston/results/convergence/icat-generic/hl-1_icat500.db'),
History('sqlite:////storage/hhecm/cellrotor/chouston/results/convergence/icat-generic/hl-1_icat750.db'),
History('sqlite:////storage/hhecm/cellrotor/chouston/results/convergence/icat-generic/hl-1_icat1000.db'),
History('sqlite:////storage/hhecm/cellrotor/chouston/results/convergence/icat-generic/hl-1_icat1250.db'),
History('sqlite:////storage/hhecm/cellrotor/chouston/results/convergence/icat-generic/hl-1_icat2000.db'),
History('sqlite:////scratch/cph211/tmp/hl-1_icat5000.db')]
particle_num = [50,100,500,750,1000,1250,2000,5000]
In [33]:
n_samples=5000
th_samples = pd.DataFrame({})
for h, p in zip(history, particle_num):
h.id = len(h.all_runs())
df, w = h.get_distribution(m=0)
th = pd.DataFrame(df.sample(n=n_samples, weights=w, replace=True).to_dict(orient='records'))
th['particle_num'] = p
th_samples = th_samples.append(th)
In [34]:
th_samples = th_samples.melt(value_vars=th_samples.columns[:-1], id_vars=['particle_num'])
In [35]:
sns.set_context('talk')
#grid = sns.lmplot(x='particle_num', y='value', hue='variable', data=th_normalised, fit_reg=False,
# x_estimator=np.mean, x_ci='sd', sharex=True, sharey=True, aspect=2)
grid = sns.relplot(x='particle_num',y='value',row='variable',#hue='variable',#style='variable',
data=th_samples,
kind='line',n_boot=5000,ci='sd',aspect=3,height=5,
facet_kws={'sharex': 'row', 'sharey': False})
In [36]:
grid.savefig('results/convergence/icat-generic/icat_convergence.pdf')
In [ ]: