Imports


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)


INFO:myokit:Loading Myokit version 1.27.4

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'

Create ion channel model


In [5]:
from channels.icat_generic import icat as model
model.sample({})


Out[5]:
x y exp
0 -75.000000 -3.193597e-02 0
1 -70.000000 -7.633718e-02 0
2 -65.000000 -1.776760e-01 0
3 -60.000000 -3.521354e-01 0
4 -55.000000 -5.540940e-01 0
5 -50.000000 -9.073693e-01 0
6 -45.000000 -1.631584e+00 0
7 -40.000000 -2.647995e+00 0
8 -35.000000 -3.511219e+00 0
9 -30.000000 -3.895625e+00 0
10 -25.000000 -3.888756e+00 0
11 -20.000000 -3.654232e+00 0
12 -15.000000 -3.326613e+00 0
13 -10.000000 -2.962956e+00 0
14 -5.000000 -2.585967e+00 0
15 0.000000 -2.204157e+00 0
16 5.000000 -1.820552e+00 0
17 10.000000 -1.435933e+00 0
18 15.000000 -1.049758e+00 0
19 20.000000 -6.611336e-01 0
20 25.000000 -2.738717e-01 0
21 30.000000 9.173555e-02 0
22 35.000000 3.996382e-01 0
0 -75.000000 4.007594e-03 1
1 -70.000000 1.006477e-02 1
2 -65.000000 2.467611e-02 1
3 -60.000000 5.166265e-02 1
4 -55.000000 8.614934e-02 1
5 -50.000000 1.500400e-01 1
6 -45.000000 2.881004e-01 1
... ... ... ...
2 2.997984 -8.689363e-01 4
3 4.843272 -9.920756e-01 4
4 6.382910 -9.853281e-01 4
5 8.468817 -9.093501e-01 4
6 10.844118 -7.787269e-01 4
7 13.939650 -6.059574e-01 4
8 20.527411 -3.289630e-01 4
9 29.407557 -1.376171e-01 4
10 41.002796 -4.332088e-02 4
11 54.664434 -1.115123e-02 4
12 71.377707 -2.098443e-03 4
13 84.473564 -5.762544e-04 4
14 96.158223 -1.765617e-04 4
15 113.487676 -3.251330e-05 4
16 131.028484 -8.226175e-06 4
17 145.343695 -6.051083e-06 4
18 157.436431 -1.153113e-06 4
0 -80.000000 -4.115745e-01 5
1 -70.000000 -1.642188e-01 5
2 -60.000000 -2.169239e-01 5
3 -50.000000 -1.455369e-02 5
4 -40.000000 -3.177721e-04 5
5 -30.000000 -2.792228e-04 5
6 -20.000000 -2.441804e-04 5
7 -10.000000 -1.973035e-04 5
8 0.000000 -1.482416e-04 5
9 10.000000 -9.883162e-05 5
10 20.000000 -4.921528e-05 5
11 30.000000 1.345231e-08 5
12 40.000000 4.569364e-05 5

94 rows × 3 columns

Get experimental measurements


In [6]:
measurements = model.get_experiment_data()
obs = measurements.to_dict()['y']
exp = measurements.to_dict()['exp']
errs = measurements.to_dict()['errs']

Set limits and generate uniform initial priors


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()})

Parameter sensitivity analysis


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')

Weights of distance function


In [17]:
from ionchannelABC import plot_distance_weights
grid = plot_distance_weights(model, distance_fn)



In [18]:
grid.savefig('results/icat/dist_weights.pdf')

Initialize pyabc database


In [9]:
db_path = ("sqlite:///" +
           os.path.join(tempfile.gettempdir(), "hl-1_icat5000.db"))
print(db_path)


sqlite:////scratch/cph211/tmp/hl-1_icat5000.db

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())


DEBUG:ABC:ion channel weights: {0: 0.3437158292831627, 1: 0.3437158292831627, 2: 0.3437158292831627, 3: 0.3437158292831627, 4: 0.3437158292831627, 5: 0.3437158292831627, 6: 0.23359946175355675, 7: 0.19709954585456368, 8: 0.19709954585456368, 9: 0.19709954585456368, 10: 0.19406724514910909, 11: 0.20345759572083993, 12: 0.22525662383378656, 13: 0.21023951557820128, 14: 0.23800699876777479, 15: 0.27422545510200186, 16: 0.30034216511171646, 17: 0.32344540858184795, 18: 0.3437158292831627, 19: 0.3437158292831627, 20: 0.3437158292831627, 21: 0.3437158292831627, 22: 0.3437158292831627, 23: 0.9694434159168854, 24: 0.9694434159168854, 25: 0.9694434159168854, 26: 0.9694434159168854, 27: 0.9694434159168854, 28: 0.9694434159168854, 29: 0.9694434159168854, 30: 0.9694434159168854, 31: 0.9694434159168854, 32: 0.9694434159168854, 33: 0.9694434159168854, 34: 0.9694434159168854, 35: 0.9694434159168854, 36: 0.9694434159168854, 37: 0.9694434159168854, 38: 0.9694434159168854, 39: 0.9694434159168854, 40: 1.175363791989361, 41: 1.175363791989361, 42: 1.175363791989361, 43: 1.175363791989361, 44: 1.175363791989361, 45: 1.175363791989361, 46: 1.175363791989361, 47: 1.175363791989361, 48: 1.175363791989361, 49: 1.175363791989361, 50: 1.175363791989361, 51: 1.175363791989361, 52: 2.5647022547867855, 53: 2.168846037200134, 54: 1.8139439583855623, 55: 1.6355232411673097, 56: 1.8475355131704811, 57: 1.6355232411673097, 58: 2.558126095159128, 59: 2.5647022547867855, 60: 2.5647022547867855, 61: 2.5647022547867855, 62: 1.113390163577478, 63: 1.113390163577478, 64: 1.113390163577478, 65: 1.113390163577478, 66: 1.113390163577478, 67: 1.113390163577478, 68: 1.113390163577478, 69: 1.113390163577478, 70: 1.113390163577478, 71: 1.113390163577478, 72: 1.113390163577478, 73: 1.113390163577478, 74: 1.113390163577478, 75: 1.113390163577478, 76: 1.113390163577478, 77: 1.113390163577478, 78: 1.113390163577478, 79: 1.113390163577478, 80: 1.113390163577478, 81: 1.0590433063764761, 82: 1.0590433063764761, 83: 1.0590433063764761, 84: 1.0590433063764761, 85: 1.0590433063764761, 86: 1.0590433063764761, 87: 1.0590433063764761, 88: 1.0590433063764761, 89: 1.0590433063764761, 90: 1.0590433063764761, 91: 1.0590433063764761, 92: 1.0590433063764761, 93: 1.0590433063764761}
DEBUG:Epsilon:init quantile_epsilon initial_epsilon=from_sample, quantile_multiplier=1

In [14]:
abc_id = abc.new(db_path, obs)


INFO:History:Start <ABCSMC(id=2, start_time=2018-11-27 10:36:14.070994, end_time=None)>
INFO:Epsilon:initial epsilon is 60.43595344679339

In [ ]:
history = abc.run(minimum_epsilon=0.1, max_nr_populations=100, min_acceptance_rate=0.01)


INFO:ABC:t:70 eps:0.894952595980924

Results analysis


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 [ ]: