In [49]:
# 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 [50]:
# Custom imports
from ionchannelABC import (ion_channel_sum_stats_calculator,
IonChannelAcceptor,
IonChannelDistance,
EfficientMultivariateNormalTransition,
calculate_parameter_sensitivity,
plot_parameter_sensitivity,
plot_regression_fit,
plot_parameters_kde,
plot_sim_results,
plot_distance_weights)
In [51]:
# Other necessary imports
import numpy as np
import subprocess
import pandas as pd
import io
import os
import tempfile
In [52]:
# Plotting imports
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
%config Inline.Backend.figure_format = 'retina'
In [53]:
from channels.iha_generic import iha as model
#model.sample({})
In [54]:
measurements = model.get_experiment_data()
obs = measurements.to_dict()['y']
exp = measurements.to_dict()['exp']
errs = measurements.to_dict()['errs']
In [55]:
limits = dict(g_ha=(0, 0.2),
Vhalf_y=(-200, 200),
k_y=(-50, 0),
c_by=(0, 1000),
c_ay=(0, 1e6),
sigma_y=(0, 100),
Vmax_y=(-200, 200),
k_i_haNa=(0, 1.0))
In [8]:
parameters = ['iha.'+k for k in limits.keys()]
In [9]:
distance_fn=IonChannelDistance(
obs=obs,
exp_map=exp,
err_bars=errs,
err_th=0.1)
In [11]:
sns.set_context('talk')
g = plot_distance_weights(model, distance_fn)
In [12]:
g.savefig('results/iha-generic/dist-weights.pdf')
In [13]:
fitted, regression_fit, r2 = calculate_parameter_sensitivity(
model,
parameters,
distance_fn,
sigma=0.1,
n_samples=1000)
In [14]:
grid1 = plot_parameter_sensitivity(fitted, plot_cutoff=0.05)
In [15]:
grid2 = plot_regression_fit(regression_fit, r2)
In [16]:
grid1.savefig('results/iha-generic/sensitivity.pdf')
grid2.savefig('results/iha-generic/regression_fit.pdf')
In [ ]:
# Finding insensitive parameters
cutoff = 0.05
fitted_pivot = fitted.pivot(index='param',columns='exp')
insensitive_params = fitted_pivot[(abs(fitted_pivot['beta'][0])<cutoff) & (abs(fitted_pivot['beta'][1])<cutoff) &
(abs(fitted_pivot['beta'][2])<cutoff)].index.values
In [ ]:
insensitive_limits = dict((k, limits[k[4:]]) for k in insensitive_params)
insensitive_prior = Distribution(**{key: RV("uniform", a, b - a)
for key, (a,b) in insensitive_limits.items()})
In [ ]:
# Generate random samples for insensitive parameters
def generate_sample(insensitive_prior, n):
samples = [dict() for i in range(n)]
for i in range(n):
parameters = insensitive_prior.rvs()
sample = {key: value for key, value in parameters.items()}
samples[i].update(sample)
return samples
In [ ]:
samples = generate_sample(insensitive_prior, 1000)
In [ ]:
model.add_external_par_samples(samples)
In [ ]:
# parameters passed to ABC algorithm
limits = dict((k, limits[k]) for k in limits if k[4:] not in insensitive_params)
In [10]:
prior = Distribution(**{key: RV("uniform", a, b - a)
for key, (a,b) in limits.items()})
In [11]:
db_path = ('sqlite:///' +
os.path.join(tempfile.gettempdir(), "hl-1_iha-generic.db"))
print(db_path)
In [12]:
# 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 [13]:
from pyabc.populationstrategy import ConstantPopulationSize
In [14]:
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(2000),
#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=6),
acceptor=IonChannelAcceptor())
In [15]:
abc_id = abc.new(db_path, obs)
In [16]:
history = abc.run(minimum_epsilon=0.01, max_nr_populations=30, min_acceptance_rate=0.01)