The goal of this script is to map between the results obtained using the Fourstate class and those obtained with the BestMSM package. In particular we care about the estimation of fluxes and paths for the four state model. First we create the instances of the Fourstate and BestMSM.MSM classes.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import fourstate
import itertools
import networkx as nx
import numpy as np
bhs = fourstate.FourState()
This much we knew already. Now we just have to produce an instance of the BestMSM.MSM class that has the same count and transition matrices.
In [2]:
import bestmsm.msm as msm
bhsmsm = msm.MSM(keys = range(4))
bhsmsm.count = bhs.count
bhsmsm.keep_states = range(4)
bhsmsm.keep_keys = range(4)
bhsmsm.trans = bhs.trans
bhsmsm.rate = bhs.K
bhsmsm.tauK, bhsmsm.peqK = bhsmsm.calc_eigsK()
So now we already have a transition matrix, eigenvalues and eigenvectors which are the same for the Fourstate and BestMSM.MSM class. The important bits start next, with the flux and pfold estimations.
In [3]:
bhs.run_commit()
bhsmsm.do_pfold(UU=[0], FF=[3])
print " These are the flux matrices"
print bhs.J
print bhsmsm.J
print " And these are the total flux values"
print " from Fourstate: %g"%bhs.sum_flux
print " from BestMSM.MSM: %g"%bhs.sum_flux
It really looks like we have got the same exact thing, as intended. Next comes the generation of paths. This implies defining a function in the case of the Fourstate object.
In [4]:
def gen_path_lengths(keys, J, pfold, flux, FF, UU):
nkeys = len(keys)
I = [x for x in range(nkeys) if x not in FF+UU]
Jnode = []
# calculate flux going through nodes
for i in range(nkeys):
Jnode.append(np.sum([J[i,x] for x in range(nkeys) \
if pfold[x] < pfold[i]]))
# define matrix with edge lengths
Jpath = np.zeros((nkeys, nkeys), float)
for i in UU:
for j in I + FF:
if J[j,i] > 0:
Jpath[j,i] = np.log(flux/J[j,i]) + 1 # I add 1
for i in I:
for j in [x for x in FF+I if pfold[x] > pfold[i]]:
if J[j,i] > 0:
Jpath[j,i] = np.log(Jnode[j]/J[j,i]) + 1 # I add 1
return Jnode, Jpath
Jnode, Jpath = gen_path_lengths(range(4), bhs.J, bhs.pfold, \
bhs.sum_flux, [3], [0])
JpathG = nx.DiGraph(Jpath.transpose())
tot_flux = 0
for path in nx.all_simple_paths(JpathG, 0, 3):
print path
f = bhs.J[path[1],path[0]]
print "%2i -> %2i: %10.4e "%(path[0], path[1], \
bhs.J[path[1],path[0]])
for i in range(2, len(path)):
print "%2i -> %2i: %10.4e %10.4e"%(path[i-1], path[i], \
bhs.J[path[i],path[i-1]], Jnode[path[i-1]])
f *= bhs.J[path[i],path[i-1]]/Jnode[path[i-1]]
tot_flux +=f
print " J(path) = %10.4e"%f
print
print " Commulative flux: %10.4e"%tot_flux
For BestMSM.MSM everything should be built in. First we obtain the 4 highest flux paths.
In [5]:
bhsmsm.do_dijkstra(UU=[0], FF=[3], npath=3)
Out[5]:
Then we use the alternative mechanism of giving a cutoff for the flux left. Here we want to account for 80% of the flux.
In [11]:
bhsmsm.do_pfold(UU=[0], FF=[3])
bhsmsm.do_dijkstra(UU=[0], FF=[3], cut=0.2)
Out[11]:
So at least for this simple example, things are working perfectly fine.