In [1]:

    

# PyABC imports
from pyabc import (ABCSMC, Distribution, RV,
                   History, MedianEpsilon)
from pyabc.populationstrategy import ConstantPopulationSize, 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,
                           plot_parameter_sensitivity,
                           plot_parameters_kde,
                           plot_distance_weights)


    

    




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'


    

In [5]:

    

from channels.ik1 import ik1 as model
#model.sample({})


    

In [6]:

    

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

# Normalise as conductance is not fit in this notebook
obs = {k: v / abs(max(obs.values(), key=abs)) for k, v in obs.items()}


    

In [29]:

    

measurements['y'] = obs.values()


    

In [7]:

    

limits = dict(k_1=(-1000, 1000),
              k_2=(0, 100),
              k_3=(0, 1),
              k_4=(0, 0.1))
prior = Distribution(**{key: RV("uniform", a, b - a)
                        for key, (a,b) in limits.items()})


    

In [8]:

    

parameters = ['ik1.'+k for k in limits.keys()]


    

In [9]:

    

distance_fn=IonChannelDistance(
    obs=obs,
    exp_map=exp,
    err_bars=errs,
    err_th=0.1)


    

In [10]:

    

sns.set_context('talk')
g = plot_distance_weights(model, distance_fn)
g.savefig('results/ik1/dist_weights.pdf')


    

In [12]:

    

grid1, grid2 = plot_parameter_sensitivity(
    model,
    parameters,
    distance_fn,
    sigma=0.1,
    n_samples=1000,
    plot_cutoff=0.05)


    

In [13]:

    

grid1.savefig('results/ik1/sensitivity.pdf')
grid2.savefig('results/ik1/sensitivity_fit.pdf')


    

In [14]:

    

# Update limits from study
limits = dict(k_2=(0, 100),
              k_3=(0, 1),
              k_4=(0, 0.1))
prior = Distribution(**{key: RV("uniform", a, b - a)
                        for key, (a,b) in limits.items()})


    

In [15]:

    

db_path = ('sqlite:///' + 
           os.path.join(tempfile.gettempdir(), "hl-1_ik1.db"))
print(db_path)


    

    




sqlite:////var/folders/18/lbp63zt51g90fcjs73tbh1x80000gn/T/hl-1_ik1.db

In [16]:

    

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


    

In [17]:

    

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(
             #    nr_particles=10000),
             population_size=AdaptivePopulationSize(
                 start_nr_particles=5000,
                 mean_cv=0.2,
                 max_population_size=5000,
                 min_population_size=2500),
             summary_statistics=ion_channel_sum_stats_calculator,
             transitions=EfficientMultivariateNormalTransition(),
             eps=MedianEpsilon(),
             sampler=MulticoreEvalParallelSampler(n_procs=1),
             acceptor=IonChannelAcceptor())


    

    




DEBUG:ABC:ion channel weights: {0: 1.0000000000000002, 1: 1.0000000000000002, 2: 1.0000000000000002, 3: 1.0000000000000002, 4: 1.0000000000000002, 5: 1.0000000000000002, 6: 1.0000000000000002, 7: 1.0000000000000002, 8: 1.0000000000000002, 9: 1.0000000000000002, 10: 1.0000000000000002, 11: 1.0000000000000002, 12: 1.0000000000000002, 13: 1.0000000000000002, 14: 1.0000000000000002, 15: 1.0000000000000002, 16: 1.0000000000000002, 17: 1.0000000000000002, 18: 1.0000000000000002}
DEBUG:Epsilon:init quantile_epsilon initial_epsilon=from_sample, quantile_multiplier=1

In [18]:

    

abc_id = abc.new(db_path, obs)


    

    




INFO:History:Start <ABCSMC(id=1, start_time=2018-09-11 11:30:59.780491, end_time=None)>
/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/pyabc/epsilon.py:321: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  distances = weighted_distances.distance.as_matrix()
/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/pyabc/epsilon.py:325: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  weights = weighted_distances.w.as_matrix()
INFO:Epsilon:initial epsilon is 0.5793157463461251

In [19]:

    

history = abc.run(minimum_epsilon=0.01, max_nr_populations=30, min_acceptance_rate=0.01)


    

    




INFO:ABC:t:0 eps:0.5793157463461251
DEBUG:ABC:now submitting population 0
DEBUG:ABC:population 0 done
DEBUG:ABC:
total nr simulations up to t =0 is 10136
/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/pyabc/epsilon.py:321: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  distances = weighted_distances.distance.as_matrix()
/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/pyabc/epsilon.py:325: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  weights = weighted_distances.w.as_matrix()
DEBUG:Epsilon:new eps, t=1, eps=0.269629969974653
INFO:ABC:t:1 eps:0.269629969974653
/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/pyabc/storage/history.py:200: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  w_arr = w.w.as_matrix()
/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/pyabc/transition/multivariatenormal.py:64: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  self._X_arr = X.as_matrix()
/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/pyabc/smc.py:729: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  self.history.max_t)["p"].as_matrix()
INFO:Adaptation:Change nr particles 5000 -> 2500
DEBUG:ABC:now submitting population 1
DEBUG:ABC:population 1 done
DEBUG:ABC:
total nr simulations up to t =1 is 14999
DEBUG:Epsilon:new eps, t=2, eps=0.16882501475510972
INFO:ABC:t:2 eps:0.16882501475510972
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 2
DEBUG:ABC:population 2 done
DEBUG:ABC:
total nr simulations up to t =2 is 19855
DEBUG:Epsilon:new eps, t=3, eps=0.13107281302042798
INFO:ABC:t:3 eps:0.13107281302042798
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 3
DEBUG:ABC:population 3 done
DEBUG:ABC:
total nr simulations up to t =3 is 25016
DEBUG:Epsilon:new eps, t=4, eps=0.11292366403290414
INFO:ABC:t:4 eps:0.11292366403290414
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 4
DEBUG:ABC:population 4 done
DEBUG:ABC:
total nr simulations up to t =4 is 30388
DEBUG:Epsilon:new eps, t=5, eps=0.09616173077857348
INFO:ABC:t:5 eps:0.09616173077857348
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 5
DEBUG:ABC:population 5 done
DEBUG:ABC:
total nr simulations up to t =5 is 35572
DEBUG:Epsilon:new eps, t=6, eps=0.08187938374641765
INFO:ABC:t:6 eps:0.08187938374641765
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 6
DEBUG:ABC:population 6 done
DEBUG:ABC:
total nr simulations up to t =6 is 40729
DEBUG:Epsilon:new eps, t=7, eps=0.07213287408895694
INFO:ABC:t:7 eps:0.07213287408895694
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 7
DEBUG:ABC:population 7 done
DEBUG:ABC:
total nr simulations up to t =7 is 45771
DEBUG:Epsilon:new eps, t=8, eps=0.06495080865338006
INFO:ABC:t:8 eps:0.06495080865338006
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 8
DEBUG:ABC:population 8 done
DEBUG:ABC:
total nr simulations up to t =8 is 50791
DEBUG:Epsilon:new eps, t=9, eps=0.06007023776234479
INFO:ABC:t:9 eps:0.06007023776234479
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 9
DEBUG:ABC:population 9 done
DEBUG:ABC:
total nr simulations up to t =9 is 55900
DEBUG:Epsilon:new eps, t=10, eps=0.05709946777333504
INFO:ABC:t:10 eps:0.05709946777333504
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 10
DEBUG:ABC:population 10 done
DEBUG:ABC:
total nr simulations up to t =10 is 60909
DEBUG:Epsilon:new eps, t=11, eps=0.05526965735883426
INFO:ABC:t:11 eps:0.05526965735883426
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 11
DEBUG:ABC:population 11 done
DEBUG:ABC:
total nr simulations up to t =11 is 66043
DEBUG:Epsilon:new eps, t=12, eps=0.0542213206619916
INFO:ABC:t:12 eps:0.0542213206619916
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 12
DEBUG:ABC:population 12 done
DEBUG:ABC:
total nr simulations up to t =12 is 71066
DEBUG:Epsilon:new eps, t=13, eps=0.05350879494690451
INFO:ABC:t:13 eps:0.05350879494690451
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 13
DEBUG:ABC:population 13 done
DEBUG:ABC:
total nr simulations up to t =13 is 76214
DEBUG:Epsilon:new eps, t=14, eps=0.05304893131857799
INFO:ABC:t:14 eps:0.05304893131857799
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 14
DEBUG:ABC:population 14 done
DEBUG:ABC:
total nr simulations up to t =14 is 81478
DEBUG:Epsilon:new eps, t=15, eps=0.052788110179169366
INFO:ABC:t:15 eps:0.052788110179169366
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 15
DEBUG:ABC:population 15 done
DEBUG:ABC:
total nr simulations up to t =15 is 86622
DEBUG:Epsilon:new eps, t=16, eps=0.05260564668181428
INFO:ABC:t:16 eps:0.05260564668181428
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 16
DEBUG:ABC:population 16 done
DEBUG:ABC:
total nr simulations up to t =16 is 91757
DEBUG:Epsilon:new eps, t=17, eps=0.05248802336143141
INFO:ABC:t:17 eps:0.05248802336143141
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 17
DEBUG:ABC:population 17 done
DEBUG:ABC:
total nr simulations up to t =17 is 96976
DEBUG:Epsilon:new eps, t=18, eps=0.05241994878331866
INFO:ABC:t:18 eps:0.05241994878331866
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 18
DEBUG:ABC:population 18 done
DEBUG:ABC:
total nr simulations up to t =18 is 102084
DEBUG:Epsilon:new eps, t=19, eps=0.05237417858691864
INFO:ABC:t:19 eps:0.05237417858691864
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 19
DEBUG:ABC:population 19 done
DEBUG:ABC:
total nr simulations up to t =19 is 107207
DEBUG:Epsilon:new eps, t=20, eps=0.0523454454425288
INFO:ABC:t:20 eps:0.0523454454425288
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 20
DEBUG:ABC:population 20 done
DEBUG:ABC:
total nr simulations up to t =20 is 112373
DEBUG:Epsilon:new eps, t=21, eps=0.05232708717581424
INFO:ABC:t:21 eps:0.05232708717581424
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 21
DEBUG:ABC:population 21 done
DEBUG:ABC:
total nr simulations up to t =21 is 117600
DEBUG:Epsilon:new eps, t=22, eps=0.052315796226907636
INFO:ABC:t:22 eps:0.052315796226907636
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 22
DEBUG:ABC:population 22 done
DEBUG:ABC:
total nr simulations up to t =22 is 122696
DEBUG:Epsilon:new eps, t=23, eps=0.05230858880982126
INFO:ABC:t:23 eps:0.05230858880982126
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 23
DEBUG:ABC:population 23 done
DEBUG:ABC:
total nr simulations up to t =23 is 127711
DEBUG:Epsilon:new eps, t=24, eps=0.0523043485199133
INFO:ABC:t:24 eps:0.0523043485199133
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 24
DEBUG:ABC:population 24 done
DEBUG:ABC:
total nr simulations up to t =24 is 132857
DEBUG:Epsilon:new eps, t=25, eps=0.05230154874720897
INFO:ABC:t:25 eps:0.05230154874720897
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 25
DEBUG:ABC:population 25 done
DEBUG:ABC:
total nr simulations up to t =25 is 138104
DEBUG:Epsilon:new eps, t=26, eps=0.05229987121641676
INFO:ABC:t:26 eps:0.05229987121641676
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 26
DEBUG:ABC:population 26 done
DEBUG:ABC:
total nr simulations up to t =26 is 143246
DEBUG:Epsilon:new eps, t=27, eps=0.05229879612654123
INFO:ABC:t:27 eps:0.05229879612654123
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 27
DEBUG:ABC:population 27 done
DEBUG:ABC:
total nr simulations up to t =27 is 148449
DEBUG:Epsilon:new eps, t=28, eps=0.052298111669235016
INFO:ABC:t:28 eps:0.052298111669235016
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 28
DEBUG:ABC:population 28 done
DEBUG:ABC:
total nr simulations up to t =28 is 153686
DEBUG:Epsilon:new eps, t=29, eps=0.05229771059485619
INFO:ABC:t:29 eps:0.05229771059485619
INFO:Adaptation:Change nr particles 2500 -> 2500
DEBUG:ABC:now submitting population 29
DEBUG:ABC:population 29 done
DEBUG:ABC:
total nr simulations up to t =29 is 158742
DEBUG:Epsilon:new eps, t=30, eps=0.05229745336337088
INFO:History:Done <ABCSMC(id=1, start_time=2018-09-11 11:30:59.780491, end_time=2018-09-11 12:51:32.807986)>

In [21]:

    

db_path = 'sqlite:///results/ik1/hl-1_ik1.db'
history = History(db_path)
history.all_runs()


    

    
Out[21]:





[<ABCSMC(id=1, start_time=2018-09-11 11:30:59.780491, end_time=2018-09-11 12:51:32.807986)>]

In [22]:

    

history.id = 1


    

In [23]:

    

sns.set_context('talk')
evolution = history.get_all_populations()
grid = sns.relplot(x='t', y='epsilon', size='samples', data=evolution[evolution.t>=0])
grid.savefig('results/ik1/eps_evolution.pdf')


    

In [24]:

    

df, w = history.get_distribution(m=0)


    

    




/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/pyabc/storage/history.py:200: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  w_arr = w.w.as_matrix()

In [43]:

    

g = plot_parameters_kde(df, w, limits, aspect=5, height=0.8)
g.savefig('results/ik1/parameters_kde.pdf')


    

    




/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/pyabc/transition/multivariatenormal.py:64: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  self._X_arr = X.as_matrix()






    

In [26]:

    

# Generate parameter samples
n_samples = 100
df, w = history.get_distribution(m=0)
th_samples = df.sample(n=n_samples, weights=w, replace=True).to_dict(orient='records')


    

    




/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/pyabc/storage/history.py:200: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  w_arr = w.w.as_matrix()

In [27]:

    

# Generate sim results samples
samples = pd.DataFrame({})
for i, th in enumerate(th_samples):
    output = model.sample(pars=th, n_x=50)
    output['sample'] = i
    output['distribution'] = 'post'
    samples = samples.append(output, ignore_index=True)


    

In [35]:

    

from ionchannelABC import plot_sim_results
sns.set_context('talk')
g = plot_sim_results(samples, obs=measurements)

# Set axis labels
xlabels = ["voltage, mV"]
ylabels = ["normalised current density"]
for ax, xl in zip(g.axes.flatten(), xlabels):
    ax.set_xlabel(xl)
for ax, yl in zip(g.axes.flatten(), ylabels):
    ax.set_ylabel(yl)


    

    




/Users/charles/miniconda3/envs/ionchannelABC/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval






    

In [36]:

    

g.savefig('results/ik1/ik1_sim_results.pdf')


    

In [ ]: