In [1]:
from pyabc import (ABCSMC, Distribution, RV,
History, Model,
ModelResult, MedianEpsilon)
from pyabc.populationstrategy import AdaptivePopulationSize
from pyabc.distance_functions import PNormDistance, WeightedPNormDistance
from pyabc.transition import MultivariateNormalTransition
from pyabc.visualization import plot_kde_matrix
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
%config Inline.Backend.figure_format = 'retina'
import os
import tempfile
import math
import pandas as pd
import scipy as sp
import numpy as np
import subprocess
from io import BytesIO
from pyabc_custom import simulate, voltage_dependence
In [2]:
myokit_python = ("/tmp/chouston/miniconda3/envs" +
"/ion_channel_ABC/bin/python")
args = [myokit_python, "get_measurements.py"]
args.append('ina_nrvm')
re = subprocess.run(args, stdout=subprocess.PIPE)
measurements = pd.read_table(BytesIO(re.stdout),
delim_whitespace=True,
header=0, index_col=False)
obs = measurements.to_dict()['y']
In [12]:
limits = dict(
p1=(0, 100),
p2=(-100, 0),
p3=(0, 1),
p4=(0, 100),
p5=(-100, 0), # Expanded from (-10, 0) after test run
p6=(0, 1),
p7=(0, 1000), # Expanded from (0, 500) after test run
q1=(0, 100),
q2=(0, 50))
# q3=(0, 10),
# q4=(-50, 50),
# q5=(-50, 50),
# q6=(0, 1),
# q7=(0, 100),
# q8=(-50, 50),
# q9=(0, 10),
# q10=(0, 1),
# q11=(0, 10),
# r1=(0, 100),
# r2=(-1, 0),
# r3=(-50, 50),
# r4=(-10, 0),
# r5=(0, 10),
# r6=(0, 100),
# r7=(0, 1),
# r8=(0, 100),
# r9=(-10, 0),
# r10=(0, 1),
# r11=(0, 10),
# r12=(-0.1, 0),
# r13=(0, 1),
# r14=(-1, 0),
# r15=(-1, 0),
# r16=(0, 100))
prior = Distribution(**{key: RV("uniform", a, b - a)
for key, (a,b) in limits.items()})
In [4]:
prior.rvs()
Out[4]:
In [5]:
distance = PNormDistance(p=2)
#distance = WeightedPNormDistance(p=2, adaptive=True)
In [6]:
class MyokitSimulation(Model):
def sample(self, pars):
try:
res = simulate('ina_nrvm', **pars).to_dict()['y']
except:
res = {key: float("inf") for key in range(len(obs))}
return res
def accept(self, pars, sum_stats_calculator, distance_calculator, eps) \
-> ModelResult:
distance_result = self.distance(pars, sum_stats_calculator,
distance_calculator)
if math.isnan(distance_result.distance):
distance_result.accepted = False
return distance_result
accepted = distance_result.distance <= eps
distance_result.accepted = accepted
return distance_result
In [7]:
MyokitSimulation().sample({})
Out[7]:
In [8]:
db_path = ("sqlite:///" +
os.path.join(tempfile.gettempdir(), "ina_korhonen.db"))
In [9]:
abc = ABCSMC(models=MyokitSimulation(),
parameter_priors=prior,
distance_function=distance,
population_size=AdaptivePopulationSize(
5000, 0.5,
min_population_size=500,
max_population_size=10000),
eps=MedianEpsilon(100, median_multiplier=1.0))
abc_id = abc.new(db_path, obs)
In [14]:
history = abc.run(max_nr_populations=20, minimum_epsilon=0.01)
In [15]:
df, w = history.get_distribution(m=0)
plot_kde_matrix(df, w, limits=limits)
Out[15]:
In [ ]:
history = abc.run(max_nr_populations=20, minimum_epsilon=0.01)
Note: from below assume the database file has been copied into the directory from /tmp/
In [7]:
history = History('sqlite:///pyabc-runs/ina_korhonen.db')
In [8]:
history.all_runs()
Out[8]:
In [ ]:
from IPython.display import HTML
from visualization_custom import animate_kde_matrix
ani = animate_kde_matrix(history, n_frames=50, limits=limits, colorbar=False)
HTML(ani.to_jshtml())
#ani.save('converging.gif', dpi=80)
In [11]:
from visualization_custom import plot_sim_results
g = plot_sim_results(history, 'ina_nrvm', n_samples=5, obs=measurements)
xlabels = ["Voltage (mV)", "Voltage (mV)", "Voltage (mV)", "Time (msec)", "Time (msec)", "Time (msec)"]
for ax, xl in zip(g.axes.flatten(), xlabels):
ax.set_xlabel(xl)
g.axes.flatten()[0].set_ylabel("Normalised output")
plt.savefig('korhonen_results.pdf', format='pdf', dpi=1000) # save to disk
In [19]:
from visualization_custom import plot_kde_2d_custom
g = plot_kde_2d_custom(history, "g_Na", "p1", times=[0, 1, 5, 10, 20], limits=limits)
#plt.savefig('convergence.pdf', format='pdf', dpi=1000)
In [20]:
df, w = history.get_distribution(m=0)
df['wt'] = w
print("Mean")
for key in limits.keys():
print(key, ":", sum(df[key] * df.wt))
print("Min")
for key in limits.keys():
print(key, ":", min(df[key]))
print("Max")
for key in limits.keys():
print(key, ":", max(df[key]))
In [27]:
out = pd.DataFrame({})
n_samples = 10
df, w = history.get_distribution(m=0)
par_samples = df.sample(n=n_samples, weights=w, replace=True).to_dict(orient='records')
for i, pars in enumerate(par_samples):
output = simulate('ina_nrvm', experiment=3, logvars=['environment.time', 'ina.i_Na'], **pars)
output['sample'] = i
out = out.append(output, ignore_index=True)
fig, ax = plt.subplots(figsize=(13, 4))
out.groupby("sample").plot(x="environment.time", y="ina.i_Na", ax=ax)
Out[27]:
In [28]:
out = simulate("ina_nrvm", experiment=0, logvars=['environment.time', 'ina.i_Na'])
fig,ax = plt.subplots(figsize=(13,13))
out.groupby("run").plot(x="environment.time", y="ina.i_Na",ax=ax)
min(out['ina.i_Na'])
Out[28]:
In [32]:
out = pd.DataFrame({})
n_samples=50
df, w = history.get_distribution(m=0)
variables = ['ina.m_ss', 'ina.tau_m', 'ina.j_ss', 'ina.tau_j', 'ina.tau_h']
par_samples = df.sample(n=n_samples, weights=w, replace=True).to_dict(orient='records')
for i, pars in enumerate(par_samples):
output = voltage_dependence('ina_nrvm', variables, **pars)
output['sample'] = i
out = out.append(output, ignore_index=True)
fig, ax = plt.subplots(figsize=(8, 4), ncols=2)
out.groupby("sample").plot(x="V", y="ina.m_ss", ax=ax[0])
out.groupby("sample").plot(x="V", y="ina.j_ss", ax=ax[0])
out.groupby("sample").plot(x="V", y="ina.tau_m", ax=ax[1])
Out[32]:
In [ ]:
abc_continued = ABCSMC(models=MyokitSimulation(),
parameter_priors=prior,
distance_function=distance,
population_size=AdaptivePopulationSize(
5000, 0.5,
min_population_size=500,
max_population_size=10000),
eps=MedianEpsilon(100, median_multiplier=1))
abc_continued.load(db_path, 2)
abc_continued.run(max_nr_populations=20, minimum_epsilon=0.01)
In [7]:
history = History('sqlite:///pyabc-runs/ina_korhonen.db')
In [8]:
history.all_runs()
Out[8]:
In [13]:
df, w = history.get_distribution(m=0)
df['wt'] = w
print("Mean")
for key in limits.keys():
print(key, ":", sum(df[key] * df.wt))
print("Min")
for key in limits.keys():
print(key, ":", min(df[key]))
print("Max")
for key in limits.keys():
print(key, ":", max(df[key]))
In [15]:
from visualization_custom import plot_sim_results
g = plot_sim_results(history, 'ina_nrvm', n_samples=100, obs=measurements)
xlabels = ["Voltage (mV)", "Voltage (mV)", "Voltage (mV)", "Time (msec)", "Time (msec)", "Time (msec)"]
for ax, xl in zip(g.axes.flatten(), xlabels):
ax.set_xlabel(xl)
g.axes.flatten()[0].set_ylabel("Normalised output")
plt.savefig('korhonen_results.pdf', format='pdf', dpi=1000) # save to disk
In [ ]: