In [8]:
%matplotlib inline
import matplotlib.pyplot as plt
import parallel
import itertools
import networkx as nx
import numpy as np
bhs = parallel.Parallel()
Again, the system exhibits a time scale separation between the first and second modes as we would expect.
In [9]:
fig, ax = plt.subplots()
ax.bar([x+0.9 for x in range(5)], -1./bhs.evals[1:], width=0.2)
ax.set_xlabel(r'Eigenvalue', fontsize=16)
ax.set_ylabel(r'$\tau_i$', fontsize=18)
ax.set_xlim([0.5,5.5])
plt.show()
In [10]:
bhs.run_commit()
We check for the property of flux conservation, i.e. $J_j=\sum_{p_{fold}(i)<p_{fold}(j)}J_{i\rightarrow j}=\sum_{p_{fold}(i)>p_{fold}(j)}J_{j\rightarrow i}$.
In [11]:
print " j J_j(<-) J_j(->)"
print " - -------- --------"
for i in range(1,5):
print "%2i %10.4e %10.4e"%(i, np.sum([bhs.J[i,x] for x in range(6) if bhs.pfold[x] < bhs.pfold[i]]),\
np.sum([bhs.J[x,i] for x in range(6) if bhs.pfold[x] > bhs.pfold[i]]))
Once that we have checked for this, we can now generate paths.
In [12]:
import tpt_functions
Jnode, Jpath = tpt_functions.gen_path_lengths(range(6), bhs.J, bhs.pfold, \
bhs.sum_flux, [5], [0])
JpathG = nx.DiGraph(Jpath.transpose())
In [13]:
tot_flux = 0
for path in nx.all_simple_paths(JpathG, 0, 5):
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
According to the enumeration of paths, the highest flux paths would be:
In [14]:
while True:
Jnode, Jpath = tpt_functions.gen_path_lengths(range(6), bhs.J, bhs.pfold, \
bhs.sum_flux, [5], [0])
# generate nx graph from matrix
JpathG = nx.DiGraph(Jpath.transpose())
# find shortest path
try:
path = nx.dijkstra_path(JpathG, 0, 5)
pathlength = nx.dijkstra_path_length(JpathG, 0, 5)
print " shortest path:", path, pathlength
except nx.NetworkXNoPath:
print " No path for %g -> %g\n Stopping here"%(0, 5)
break
# calculate contribution to flux
f = bhs.J[path[1],path[0]]
print "%2i -> %2i: %10.4e "%(path[0], path[1], bhs.J[path[1],path[0]])
path_fluxes = [f]
for j in range(2, len(path)):
i = j - 1
print "%2i -> %2i: %10.4e %10.4e"%(path[i], path[j], \
bhs.J[path[j],path[i]], \
bhs.J[path[j],path[i]]/Jnode[path[i]])
f *= bhs.J[path[j],path[i]]/Jnode[path[i]]
path_fluxes.append(bhs.J[path[j],path[i]])
# find bottleneck
ib = np.argmin(path_fluxes)
print "bottleneck: %2i -> %2i"%(path[ib],path[ib+1])
# remove flux from edges
for j in range(1,len(path)):
i = j - 1
bhs.J[path[j],path[i]] -= f
#bhs.J[path[ib],path[ib+1]] -= f
# numerically there may be some leftover flux in bottleneck
#bhs.J[path[ib],path[ib+1]] = 0.
bhs.sum_flux -= f
print ' flux from path ', path, ': %10.4e'%f
print ' fluxes', path_fluxes
print ' leftover flux: %10.4e; %10.2f \n'%(bhs.sum_flux, bhs.sum_flux/tot_flux)
if bhs.sum_flux/tot_flux < 0.01:
break
Here we find that the elimination of flux from edges fails to reproduce the correct order of paths as obtained from the exhaustive enumeration.