In [ ]:
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
import os
import glob
import pickle
import numpy as np
import pandas as pd
from collections import OrderedDict
In [3]:
from dcpyps import mechanism
from dcpyps.samples import samples
In [4]:
def constrain(mec):
for i in range(len(mec.Rates)):
mec.Rates[i].fixed = False
# Constrained rates.
mec.Rates[21].is_constrained = True
mec.Rates[21].constrain_func = mechanism.constrain_rate_multiple
mec.Rates[21].constrain_args = [20, 1.5]
mec.Rates[18].is_constrained = True
mec.Rates[18].constrain_func = mechanism.constrain_rate_multiple
mec.Rates[18].constrain_args = [19, 2]
mec.Rates[14].is_constrained = True
mec.Rates[14].constrain_func = mechanism.constrain_rate_multiple
mec.Rates[14].constrain_args = [12, 3]
mec.Rates[13].is_constrained = True
mec.Rates[13].constrain_func = mechanism.constrain_rate_multiple
mec.Rates[13].constrain_args = [12, 2]
mec.Rates[15].is_constrained = True
mec.Rates[15].constrain_func = mechanism.constrain_rate_multiple
mec.Rates[15].constrain_args = [17, 3]
mec.Rates[16].is_constrained = True
mec.Rates[16].constrain_func = mechanism.constrain_rate_multiple
mec.Rates[16].constrain_args = [17, 2]
mec.update_constrains()
mec.set_mr(True, 9, 0)
mec.set_mr(True, 11, 1)
mec.update_constrains()
return mec
In [5]:
mec_true = samples.GlyR_flip()
ig = [4200, 28000, 130000, 3400, 2100, 6700, 180, 6800, 22000,
29266, 18000, 948, 302, 604, 906, 1.77e6, 1.18e6, 0.59e6, 300e6, 150e6,
2500, 3750]
mec_true.set_rateconstants(ig)
mec_true = constrain(mec_true)
mec = samples.GlyR_flip()
mec = constrain(mec)
In [6]:
os.chdir("results")
In [7]:
r = []
x = []
lik = []
cputime = []
for file in glob.glob("*.result"):
with open(file, 'rb') as f:
content = pickle.load(f)
#print (np.exp(content[0].x))
mec.theta_unsqueeze(np.exp(content[0].x))
mec.update_constrains()
x.append(np.exp(content[0].x))
r.append(mec.unit_rates())
lik.append(-content[0].fun)
cputime.append(content[1])
rates = np.array(r)
In [8]:
rate_names = []
true_rates = []
for i in range(len(mec.unit_rates())):
rate_names.append(mec.Rates[i].name)
true_rates.append(mec_true.Rates[i].unit_rate())
In [69]:
trueratesdf = pd.DataFrame({rate_names[i]: true_rates[i] for i in range(len(rate_names))}, index=np.arange(1))
datadict = {rate_names[i]: rates[:,i] for i in range(len(rate_names))}
datadict['likelihood'] = lik
datadict['cputime'] = cputime
datadf = pd.DataFrame(datadict)
trueratesser = pd.Series(true_rates, index=rate_names)
In [72]:
trueratesser
Out[72]:
In [73]:
datadf
datadf.head()
Out[73]:
In [74]:
meandf = pd.DataFrame({'Mean' : datadf[rate_names].mean(),
'CV%' : 100*datadf[rate_names].std(ddof=0)/datadf[rate_names].mean(),
'Bias%' : (100*(datadf[rate_names].mean()/trueratesser - 1)),
'True' : trueratesser})
meandf.reindex(columns=['True', 'Mean', 'CV%', 'Bias%'])
Out[74]:
In [78]:
fig,axes = plt.subplots(11, 2, figsize=(12,35))
for i,axis in enumerate(axes.ravel()):
rate = rate_names[i]
axis.hist(datadf[rate])
axis.set_title(rate)
axis.axvline(trueratesdf[rate][0], color='r')
axis.set_title(rate)
axis.set_xlim(trueratesdf[rate][0] * 0.5, 1.5 * trueratesdf[rate][0])
fig.tight_layout()
print('Red line - true value')
In [51]:
def plot_intersection(name1, name2):
fig, ax = plt.subplots(1,1, figsize=(10,3))
ax.plot(name1, name2, 'b.', data=datadf)
ax.axvline(trueratesdf[name1][0], color='r')
ax.axhline(trueratesdf[name2][0], color='r')
ax.set_xlabel(name1)
ax.set_ylabel(name2)
fig.tight_layout()
In [52]:
name1 = 'beta3'
name2 = 'alpha3'
plot_intersection(name1, name2)
In [53]:
plot_intersection(rate_names[12], rate_names[17])
In [77]:
fig, axes = plt.subplots(1,2, figsize=(12,3))
axes[0].hist(datadf['likelihood'])
axes[0].set_xlabel('LogLikelihood')
axes[1].hist(datadf['cputime'])
axes[1].set_xlabel('CPU time, s')
print('Mean CPU time = ', np.mean(cputime))