<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Parallel pathways model
We give another go at the estimation of highest flux paths using a slightly more complicated model involving six states. I adapt the Fourstate class so that we can do the same exact thing as before. The model system looks like the following.
<img src="files/parallel.png">
Initialization
As before we initialize the system producing the transition count matrix ($\mathbf{N}$), transition probability matrix ($\mathbf{T}$) and rate matrix ($\mathbf{K}$). I slightly modify the values in the count matrix preserving its topology, so that not all numbers are equal.
Step1:  Again, the system exhibits a time scale separation between the first and second modes as we would expect.
Step2:  Committors and reactive flux
Step3:  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}$.
Step4:  Once that we have checked for this, we can now generate paths.
Step5:  According to the enumeration of paths, the highest flux paths would be
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
import parallel
import itertools
import networkx as nx
import numpy as np
bhs = parallel.Parallel()
Explanation: Parallel pathways model
We give another go at the estimation of highest flux paths using a slightly more complicated model involving six states. I adapt the Fourstate class so that we can do the same exact thing as before. The model system looks like the following.
<img src="files/parallel.png">
Initialization
As before we initialize the system producing the transition count matrix ($\mathbf{N}$), transition probability matrix ($\mathbf{T}$) and rate matrix ($\mathbf{K}$). I slightly modify the values in the count matrix preserving its topology, so that not all numbers are equal.
End of explanation
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()
Explanation: Again, the system exhibits a time scale separation between the first and second modes as we would expect.
End of explanation
bhs.run_commit()
Explanation: Committors and reactive flux
End of explanation
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]]))
Explanation: 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}$.
End of explanation
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())
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
Explanation: Once that we have checked for this, we can now generate paths.
End of explanation
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
Explanation: According to the enumeration of paths, the highest flux paths would be:
0$\rightarrow$1$\rightarrow$4$\rightarrow$5; $J=$9.87e-06
0$\rightarrow$2$\rightarrow$4$\rightarrow$5; $J=$6.93e-06
0$\rightarrow$1$\rightarrow$3$\rightarrow$5; $J=$6.19e-06
0$\rightarrow$2$\rightarrow$3$\rightarrow$5; $J=$2.18e-06
Highest flux paths
Again, we check whether we can get the hierarchy of paths by removing the flux from the edges of the graph involved in the highest flux paths.
End of explanation
<END_TASK>