In [1]:
%pylab inline
import scipy
import sys
# make nice plots
import plt_fmt
if 'marktools' in sys.modules:
del sys.modules['marktools']
# from marktools import *
import marktools
from marktools.markpy import *
from marktools.marksim import *
from marktools.marknet import *
from marktools.renorm_neq import *
from marktools.lamben import *
# from msmbuilder.example_datasets import fetch_alanine_dipeptide
In [2]:
# # Reload misbehaving modules
# import sys
# if 'markpy' in sys.modules:
# del sys.modules['markpy']
# from markpy import *
# if 'marksim' in sys.modules:
# del sys.modules['marksim']
# from renorm_neq import *
# if 'marknet' in sys.modules:
# del sys.modules['marknet']
# from marknet import *
# if 'renorm_neq' in sys.modules:
# del sys.modules['renorm_neq']
# from renorm_neq import *
# if 'lamben' in sys.modules:
# del sys.modules['lamben']
# from renorm_neq import *
In [2]:
# load a sample transition matrix pickled from MSMBuilder double well tutorial
try:
import cPickle as pickle
except:
import pickle
f = open('msm_data/dbwell_10stat.pkl', 'rb')
tmat = pickle.load(f)
matprint(tmat)
figure()
imshow(tmat)
figure()
p0 = rand(10)
p0 = p0/sum(p0)
# p0 = getss(tmat)
hold(True)
for ii in range(10):
if mod(ii,1)==0:
plot( p0.dot(matrix_power(tmat, ii)) )
else:
pass
figure()
plot(getss(tmat))
Out[2]:
In [5]:
# # Enforce diagonal elements generally being larger than the rest (metastable)
# num_states = 10
# tmat = zeros([num_states, num_states])
# tmat += .1*rand(num_states,num_states)
# for ii in range(tmat.shape[0]):
# for jj in range(tmat.shape[1]):
# if ii == jj+1:
# tmat[ii, jj] += 0.0
# tmat = tmat/sum(tmat, axis=1)
# for ii in range(tmat.shape[0]):
# for jj in range(tmat.shape[1]):
# if ii == jj:
# tmat[ii, jj] = 0.0
# tmat[ii, jj] = 1.0 - sum(tmat[ii,:])
# from numpy.linalg import matrix_power
# p0 = rand(num_states)
# p0 = p0/sum(p0)
# plot(p0)
# plot( matrix_power(tmat, 30).dot(p0) )
In [6]:
# import networkx as nx
# # Make a scale-free graph (MultiDiGraph allows parallel pathways)
# G=nx.scale_free_graph(10)
# pos=nx.spring_layout(G)
# pos=nx.spectral_layout(G)
# nx.draw(G,pos)
# adj = nx.adjacency_matrix(G)
# adj = adj.toarray()
# # adj = adj/tile(sum(adj,axis=1),adj.shape[0])
# figure()
# imshow(adj)
In [4]:
# from renorm_neq import *
# qq=copy(tmat)
# figure()
# imshow(renorm_col(qq, its=1))
# figure()
# imshow( renorm_neq(qq, its=1) )
In [3]:
figure()
imshow(tmat)
figure()
imshow(renorm(tmat))
print (sum(renorm(tmat,k=2,its=2),axis=1))
In [4]:
import networkx as nx
mat2 = array([[.1,2,1],[1e-3,.1,1],[1,.1,1]])
mat2 = array([[.5,.5,0],[.9,.1,0],[.6,.3,.1]])
# matTg = mat2DiGraph(mat2)
# pos=nx.spring_layout(matTg)
# nx.draw(matTg,pos)
markov_plot(tmat)
figure()
markov_plot(renorm(tmat))
# print (getss(mat2))
In [7]:
from scipy.linalg import eig
from numpy.linalg import matrix_power
# from lamben import *
dbwl = LambEn(tmat,.3)
dbwl.fen(200)
hold(True)
for tt in range(10):
gg = list()
for ii in range(10):
kk = dbwl.px(200, 93, ii)
gg.append(kk)
plot(array(gg),'-')
In [8]:
test_traj = simtraj(tmat,300,stt=0)
plot(test_traj,'-')
figure()
plot(simtraj(renorm(tmat, k=3),300,stt=0),'-')
Out[8]:
In [61]:
figure()
imshow(make_bethe(3,5,stay_rate=0.0))
sum(make_bethe(3,2,stay_rate=0.0),axis=1)
tmat3 = make_bethe(3,5,stay_rate=0.0)
figure()
hist(log(1+simtraj(tmat3,5000)))
figure()
hist(log(1+2*simtraj(renorm(tmat3),5000)))
Out[61]:
In [190]:
tmat4 = make_bethe(3,4,stay_rate=1.0)
# tmat5 = LambEn(tmat4,.82).tiltmat
tmat5 = copy(tmat4)
tmat5 = tmat5 + tril(tmat5)*1.5
tpert = zeros(tmat5.shape)
tpert[-1,-2] = 1.0
tpert[-2,-1] = 5.0
# tpert = rand(*tmat5.shape)
# tpert[12,13] = 1e-1
# tpert[13,12] = 4e-1
tmat5 = tmat5 + tpert
tmat5 = rownorm(tmat5)
print (real(ss_ent(tmat4)))
print (real(ss_ent(tmat5)))
figure()
plot(getss(tmat4),'.')
figure()
plot(getss(tmat5),'.')
figure()
markov_plot(tmat4,scale=10000)
figure()
markov_plot(tmat5,scale=10000)
In [11]:
tmat5 = make_bethe(3,5,stay_rate=1.0)
tmat5 = tmat5 + tril(tmat5)*10.5
tmat5 = rownorm(tmat5)
nsteps = 1000
ttraj = simtraj(tmat5,nsteps,stt=8)
figure()
plot(ttraj)
figure()
cum_ent = path_ent(tmat5, ttraj)
# plot(cum_ent)
print (sum(cumsum(cum_ent)/nsteps))
print(sum(cum_ent))
print (std(cum_ent))
In [34]:
# from numpy.random import random_sample
# def weighted_values(values, probabilities, size):
# bins = cumsum(probabilities)
# return values[digitize(random_sample(size), bins)]
# values = array([1,2,3])
# probabilities = array([0.2, 0.5, 0.3])
# kk = (weighted_values(values, probabilities, 10000))
# #Sample output:
# # [ 2.2 2.2 1.1 2.2 2.2 3.3 3.3 2.2 3.3 3.3]
In [12]:
Z = 3
k= 3
latt = make_bethe(Z, k)
# ind_t = ceil((1+Z*((Z-1.)**k - 1.)/(Z-2.))/(Z-1))-2-1
# latt[ind_t, ind_t] = 0.0
imshow(latt[:20,:20])
hold(True)
# matprint( make_bethe(3,3) )
In [23]:
import networkx as nx
G=nx.scale_free_graph(30)
pos=nx.spring_layout(G)
nx.draw(G,pos)
adj = nx.adjacency_matrix(G)
adj = adj.toarray()
# adj = adj/tile(sum(adj,axis=1),adj.shape[0])
# matprint(adj)
# figure()
# imshow(tmat)
# figure()
# nx.draw(mat2DiGraph(tmat))
# figure()
# qq = renorm2(renorm1(tmat,its=1),its=1)
# imshow(qq)
# figure()
# nx.draw(mat2DiGraph(qq))
In [22]:
tmat_lambda = LambEn(tmat_unb,.01).tiltmat
tmat_lambda = rownorm2(tmat_lambda)
imshow(tmat_lambda)
Out[22]:
In [3]:
bias = 8.5
Z = 3
k = 5
tmat_unb = make_bethe(Z,k,stay_rate=0.0)
tmat_b = tmat_unb + tril(tmat_unb)*bias
tmat_b = rownorm(tmat_b)
rand_conn_frac = 1.0
noise_amp = .1
noise = (floor(rand(*tmat_b.shape)+ rand_conn_frac)) * (noise_amp*rand(*tmat_b.shape))
tmat_noise = rownorm(tmat_b + noise)
imshow(tmat_noise)
Out[3]:
Non-equilibrium steady state characterization
In [2]:
biases = linspace(1.0,30,100)
noise_amps = 10.**linspace(-4,1,11)
rand_conn_frac = 1.0 # singular entropy production if this parameter is less than one
Z = 3
k = 5
all_ent_prod = list()
all_frac_fold = list()
for noise_amp in noise_amps:
frac_fold = list()
ent_prod = list()
for bias in biases:
tmat_unb = make_bethe(Z,k,stay_rate=0.0)
tmat_b = tmat_unb + tril(tmat_unb)*bias
tmat_b = rownorm(tmat_b)
noise = (floor(rand(*tmat_b.shape)+ rand_conn_frac)) * (noise_amp*rand(*tmat_b.shape))
tmat_noise = rownorm(tmat_b + noise)
kk = real(getss(tmat_noise))
ent_prod.append(real(ss_ent(tmat_noise)))
frac_fold.append(kk[0])
all_ent_prod.append(ent_prod)
all_frac_fold.append(frac_fold)
In [148]:
dosave = True
fig1 = figure()
hold(True)
for item in all_ent_prod:
loglog(biases, array(item))
xlim([biases[0],biases[-1]])
xlabel("Bias of rates towards native state")
ylabel("Entropy production")
if dosave:
savefig('ent_v.bias.pdf')
fig2 = figure()
hold(True)
for item in all_frac_fold:
plot(biases, array(item))
xlabel("Bias of rates towards native state")
ylabel("relative occupation of native state in equilibrium")
if dosave:
savefig('occ_v_bias.pdf')
## Make legend
fig3 = figure()
tableau20 = array([(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),
(44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),
(148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),
(227, 119, 194), (247, 182, 210), (127, 127, 127), (199, 199, 199),
(188, 189, 34), (219, 219, 141), (23, 190, 207), (158, 218, 229)]);
tableau20 =tableau20/255.
for (ind, kk) in enumerate(noise_amps):
eventplot( array([log10(noise_amps[ind])]), colors = array([tableau20[ind]]),orientation='vertical',linewidths=18 )
xlim([.5,1.5])
gca().set_aspect(3)
gca().xaxis.set_visible(False)
ylabel('log(A)')
if dosave:
savefig('bias_legend.pdf')
Lambda ensemble
In [422]:
# This normalization portion needs to be checked since
# it doesn't quite seem right for the discrete case
lmda = .12
tilted = copy(( tmat**lmda )*( tmat.T**(1.0-lmda) ))
qq2 = copy(tilted)
eigset = lefteig(qq2)
max_w = eigset[0][0]
ss = eigset[1][0]
plot(ss)
umat = zeros(qq2.shape)
for (ind, row) in enumerate(umat):
umat[ind,ind] = ss[ind]
compo = real(max_w)*(umat.dot(qq2.dot(inv(umat))))
compo = real(compo).T
# compo = unhollow(compo)
figure()
imshow(tmat)
figure()
imshow(compo)
sum(compo, axis=1)
Out[422]:
In [5]:
import networkx as nx
figure()
G=nx.scale_free_graph(15)
pos=nx.spring_layout(G)
nx.draw(G,pos)
adj = nx.adjacency_matrix(G)
adj = adj.toarray()
figure()
new_adj = adj
for ii in range(5):
to_merge = unravel_index(nanmax(new_adj), new_adj.shape)
new_adj = merge_states(new_adj, to_merge)
G3=nx.from_numpy_matrix(new_adj,create_using=nx.MultiDiGraph())
pos=nx.spring_layout(G3)
nx.draw(G3,pos)
First-hitting time distribution for different values of biasing:
In [7]:
make_bethe(6,4).shape
Out[7]:
In [226]:
all_lens = list()
# uniformly sample from an exponential distribution
biases = exp( linspace(0.0,log(10), 8) )
for bias in biases:
Z = 3
k = 8
tmat_unb = make_bethe(Z,k,stay_rate=1.0)
tmat_b = tmat_unb + tril(tmat_unb)*bias
tmat_b = rownorm(tmat_b)
leaves = len(tmat_b)
ht = sim_hits(tmat_b,list(range(383,765)),[0], ntraj=20000, cutoff=1500)
traj_lens = array([len(item) for item in ht])
print ('dead ends:')
print (sum(any(isnan(traj_lens))))
print ('\n')
all_lens.append(traj_lens)
fig = figure()
hold(True)
for traj in all_lens:
hh = histogram(traj,100,density=True)
loglog(hh[1][:-1], hh[0])
xlabel('Number of timesteps')
ylabel('Fraction of trajectories')
xlim([6, 800])
# savefig('xxsemilog_fht_v_bias.pdf')
First hitting distribution for different values of connectivity:
In [22]:
all_lens = list()
# uniformly sample from an exponential distribution
Zs = range(3,7)
bias = 4.0
for Z in Zs:
k = 4
tmat_unb = make_bethe(Z,k,stay_rate=1.0)
tmat_b = tmat_unb + tril(tmat_unb)*bias
tmat_b = rownorm(tmat_b)
leaves = len(tmat_b)
sz = tmat_b.shape[0]
ht = sim_hits(tmat_b,[sz-1],[0], ntraj=20000, cutoff=1500)
traj_lens = array([len(item) for item in ht])
print ('dead ends:')
print (sum(any(isnan(traj_lens))))
print ('\n')
all_lens.append(traj_lens)
fig = figure()
hold(True)
for traj in all_lens:
hh = histogram(traj,100,density=True)
loglog(hh[1][:-1], hh[0])
xlabel('Number of timesteps')
ylabel('Fraction of trajectories')
xlim([4, 700])
# savefig('xloglog_fht_v_Z.pdf')
Plots of occupation of native state as a function of driving (which is like temperature)
In [ ]:
In [54]:
# bias = 1.25
bias = 1.0
Z = 3
k = 10
tmat_unb = make_bethe(Z,k,stay_rate=1.0)
tmat_b = tmat_unb + tril(tmat_unb)*bias
tmat_b = rownorm(tmat_b)
# rand_conn_frac = 1.0
# noise_amp = .8
# noise = (floor(rand(*tmat_b.shape)+ rand_conn_frac)) * (noise_amp*rand(*tmat_b.shape))
# tmat_noise = rownorm(tmat_b + noise)
In [25]:
ht = sim_hits(tmat_b,list(range(46,93)),[0], ntraj=1000, cutoff=2500)
In [61]:
figure()
hist(array([len(item) for item in ht]), sqrt(len(ht)))
figure()
pent = list()
for item in ht:
pent.append(path_ent(tmat_b, item))
hold(True)
for item in pent:
plot(item)
In [60]:
hold(True)
for item in ht:
plot(item,'-')
In [55]:
print (tmat_b.shape)
# imshow(tmat_b[:,:])
In [57]:
# imshow(tmat_b)
tmat2 = copy(tmat_b)
# tmat2 = copy(tmat_unb)
# tmat2[9,5] = .2
# tmat2[5,9] = .28
# (tmat2[10,23], tmat2[23,10]) = (tmat2[23,10], tmat2[10,23])
# (tmat2[11,24], tmat2[24,11]) = (tmat2[24,11], tmat2[11,24])
# (tmat2[13,29], tmat2[29,13]) = (tmat2[29,13], tmat2[13,29])
tmat2=rownorm2(tmat2)
print (real(ss_ent(tmat2)))
# ht = sim_hits(tmat2,list(range(46,93)),[0], ntraj=10000, cutoff=1500)
# ht = sim_hits(tmat2,list(range(400,700)),[0], ntraj=10000, cutoff=1500)
# ht = sim_hits(tmat2,list(range(10,21)),[0], ntraj=10000, cutoff=1500)
ht = sim_hits(tmat2,list(range(1534,3069)),[0], ntraj=10000, cutoff=1500)
traj_lens = array([len(item) for item in ht])
print ('dead ends:')
print (sum(any(isnan(traj_lens))))
print ('\n')
print(mean(traj_lens))
print(std(traj_lens))
print ('\n')
pent = list()
for item in ht:
pent.append(mean(path_ent(tmat_b, item)))
print(mean(array(pent)))
print(std(array(pent)))
In [129]:
def levy(x,c=.5,mu=0.0):
dd = x - mu
return sqrt(c/(2*pi))*exp(-c/(2*dd))/(dd**1.5)
def invgauss(x,c=2.5,mu=1.0):
dd = x - mu
return sqrt(c/(2*pi*x**3))*exp( -c*(dd**2)/(2*x*mu**2) )
semilogy( linspace(.01,10,1000), invgauss(linspace(.01,10,1000)) )
xlim([2,4])
ylim([.01,.1])
Out[129]:
In [2]:
plot(array([3,4,5,6,7,8,9,10]), array([20,30,42,57,71,88,105,114]),'.',markersize=30)
xlabel('k')
ylabel('mean number of steps')
title('mean first hitting time versus shells for driven network')
Out[2]:
In [64]:
# yerr=array([17,27,41,56,75,95,119]
plot(log(array([3,4,5,6,7,8,9])), log(array([22,34,50,69,92,117,145])),'.',markersize=30)
xlabel('k')
ylabel('mean number of steps')
title('mean first hitting time versus shells for driven network')
Out[64]:
In [4]:
def cayley_shell_inds(Z,k):
"""
For a transition matrix corresponding to the topology of a
regular Cayley tree, returns the indices of the start of each new shell
shell (numbering convention is 0 for the root, clockwise spiral
out so that the largest index node is in the furthest out shell)
"""
inds = [0,1]
tot = 0
for kk in range(1,k):
tot += Z*(Z-1)**(kk-1)
inds.append(tot+1)
inds = array(inds)
return inds
In [7]:
hist(to_merge)
Out[7]:
In [9]:
from random import choice
Z = 3
k = 6
bias = 2.0
tmat = make_bethe(Z,k,stay_rate=1.0)
tmat = copy(tmat) + tril(tmat)*bias
tmat = rownorm(tmat)
# Perturb the lattice by merging states at the same level
n_perts = 50
shell_inds = cayley_shell_inds(Z,k+1)
shells = list()
for (loc, shell_ind) in enumerate(shell_inds):
if shell_ind == 0:
pass
elif shell_ind == shell_inds[-1]:
pass
else:
shells.append(range(shell_ind, shell_inds[loc+1]-1))
tmat_def = copy(tmat)
# PROBLEM: THE INDICES CHANGE AROUND AS YOU MERGE DERP DERP DERP
to_merge = list()
for ii in range(n_perts):
inds_in_shell = array(choice(shells))
print (inds_in_shell)
groups_in_shell = inds_in_shell[ mod(inds_in_shell, Z-1)==0 ]
group_to_merge = choice(groups_in_shell)
gset = list( range(group_to_merge,group_to_merge+Z-1) )
to_merge.append(gset)
to_merge = array(to_merge)
for loc in range(len(to_merge)):
tmat_def = merge_states(tmat_def, list(to_merge[loc]) )
to_merge[to_merge>to_merge[loc]] -=1
imshow(tmat_def)
Out[9]:
In [11]:
leaf_inds = list(range(shell_inds[-2],shell_inds[-1]-1))
ht = sim_hits(tmat,leaf_inds,[0], ntraj=20000, cutoff=2500)
# ht2 = sim_hits(tmat_def,list(range(46,93)),[0], ntraj=1000, cutoff=2500)
traj_lens = array([len(item) for item in ht])
hh = histogram(traj_lens,100,density=True)
loglog(hh[1][:-1], hh[0])
hold(True)
leaf_inds = list(range(shell_inds[-2],shell_inds[-1]-1))
ht = sim_hits(tmat_def,list(range(55,139)),[0], ntraj=20000, cutoff=2500)
# ht2 = sim_hits(tmat_def,list(range(46,93)),[0], ntraj=1000, cutoff=2500)
traj_lens = array([len(item) for item in ht])
hh = histogram(traj_lens,100,density=True)
loglog(hh[1][:-1], hh[0])
Out[11]:
In [ ]: