<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 How to analyse your favourite dynamically reversible Markov model
This notebooks contains a simple example of how to analyse fluctuating currents in dynamically reversible Markov models.
NOTE
Step1:  Defining models and using the algorithm presented in Ref [2]
Step2:  A simple model
In this script, a model is a dictionary with
keys == edges (tuples of integers)
values == transition rates (symbolic expressions)
We start with the definition of simple model on four states with the following topology
Step3:  Parametrizing the transition rates
In our toy model, transitions along the outer edges occur with rate $r$ and $l$ in clockwise and counter-clockwise directions, respectively.
For convenience, we will use the substitutions 
$r = k\exp{\frac{f}{2}}$ and $l=k\exp{\frac{f}{2}}$,
which separates the anti-symmetric edge motance $f = \log{\frac{r}{l}}$ from the symmetric part $k=\sqrt{rl}$.
Transitions through the edge connecting states $0$ and $2$ have same symmetric part $b*k$, and a potentially different driving field $g$. The positive parameter $b$ can be used to turn off the center transition
Step4:  Calculating the first and second fundamental cumulants
The function getCumulants(model, chords, rates_params) from the cumulants library takes a dictionary model and a corresponding substitution list rates_params to calculate the fundamental cumulants specified by the fundamental chords given in chords.
While some consistency checks are given, you have to be sure that removing the edges contained in chords from the graph of the model yields a spanning tree, see Ref. [1].
Here, we chose the edges $(0,1)$ and $(0,2)$ as fundamental chords. The corresponding fundamental cycles are then the outer circuit and the lower half circuit (both in clockwise direction)
Step5:  Observables
The fundamental cumulants regard the counting statistics along the fundamental chords. As shown in Ref. [1], this is sufficient in order to determine the counting statistics of arbitrary observables.
In the following we consider three observables
Step6:  Exploring the parameter dependence
Our model depends on three symbolic parameters
Step7:  Plotting
Efficient (and numerically stable) plotting requires us to transform symbolic expressions into lambda functions (using the method lambdify()). The lambda functions are used to evaluate the expressions on a grid.
NOTE
<ASSISTANT_TASK:>
Python Code:
# inline plotting/interaction
%pylab inline
# replace the line above with the line below for command line scripts:
# from pylab import *
from sympy import * # symbolic python
init_printing() # pretty printing
import numpy as np # numeric python
import time # timing, for performance monitoring
# activate latex text rendering
from matplotlib import rc
rc('text', usetex=True)
Explanation: How to analyse your favourite dynamically reversible Markov model
This notebooks contains a simple example of how to analyse fluctuating currents in dynamically reversible Markov models.
NOTE: This script runs only on a python 2 kernel. 
References
[1] Fluctuating Currents in Stochastic Thermodynamics I. Gauge Invariance of Asymptotic Statistics. Wachtel, Altaner and Vollmer (2015)
[2] Fluctuating Currents in Stochastic Thermodynamics II. Energy Conversion and Nonequilibrium Response in Kinesin Models. Altaner, Wachtel and Vollmer (2015)
Set-up jupyter (ipython-notebook) and sympy libraries
End of explanation
import cumulants # cumulants.py implements the algorithm presented in Ref [2]
Explanation: Defining models and using the algorithm presented in Ref [2]
End of explanation
# We define the transition rates as symbolic expressions...
w01, w10, w02, w20, w03, w30, w12, w21, w23, w32 = \
    symbols("w_{01}, w_{10}, w_{02}, w_{20}, w_{03}, w_{30},"+\
            "w_{12}, w_{21}, w_{23}, w_{32}",\
            real=True, positive=True)
# ... and specify the model topology as a dictionary
model4State = {(0,1): w01, (1,0): w10,\
               (0,2): w02, (2,0): w20,\
               (0,3): w03, (3,0): w30,\
               (1,2): w12, (2,1): w21,\
               (2,3): w23, (3,2): w32 \
               }
Explanation: A simple model
In this script, a model is a dictionary with
keys == edges (tuples of integers)
values == transition rates (symbolic expressions)
We start with the definition of simple model on four states with the following topology:
0 – 1
| \ |
3 – 2
It features two independent cycles (cf. Ref [1]).
End of explanation
# Define the symbols for parameters, place-holders, etc.
b = symbols("b", positive=True)
f, g = symbols("f, g", real=True)
x, y      = symbols("x, y", real=True)
u, v, w   = symbols("u, v, w", positive=True)
# symmetric-antisymmetric substitutions
onehalf = Rational(1,2)
# Generate a substitution list to parametrize the symbolic transition rates
rate_cw =  exp(+onehalf*f)
rate_ccw = exp(-onehalf*f)
# Outer circle, clockwise
w01p = rate_cw
w12p = rate_cw
w23p = rate_cw
w30p = rate_cw
# Outer circle, counter-clockwise
w03p = rate_ccw
w32p = rate_ccw
w21p = rate_ccw
w10p = rate_ccw
# Center rates
w02p = b*exp(+onehalf*g)  
w20p = b*exp(-onehalf*g)
# Create the substitution list that gives the parametric dependence of the rates:
rates_parametrized = [(w01,w01p),(w10,w10p),(w02,w02p),(w20,w20p),\
           (w12,w12p),(w21,w21p),(w23,w23p),(w32,w32p),(w30,w30p),(w03,w03p)]
Explanation: Parametrizing the transition rates
In our toy model, transitions along the outer edges occur with rate $r$ and $l$ in clockwise and counter-clockwise directions, respectively.
For convenience, we will use the substitutions 
$r = k\exp{\frac{f}{2}}$ and $l=k\exp{\frac{f}{2}}$,
which separates the anti-symmetric edge motance $f = \log{\frac{r}{l}}$ from the symmetric part $k=\sqrt{rl}$.
Transitions through the edge connecting states $0$ and $2$ have same symmetric part $b*k$, and a potentially different driving field $g$. The positive parameter $b$ can be used to turn off the center transition:
$w^0_2 = kb\exp{\frac{g}{2}}$ and $w^2_0 = kb\exp{\frac{g}{2}}$
As a global factor of the transition matrix, $k$ can be absorbed into the time-scale of the model and we set $k\equiv 1$.
End of explanation
# calculate the cumulants (this takes only a few seconds)
c, C = cumulants.getCumulants( model4State, [(0,1),(0,2)], rates_parametrized)
## WARNING: READ BEFORE EXECUTING THIS CELL
##
## Additional simplifications may help if numerical problems are encountered (e.g. for very stiff rate matrices).
## They can be very time-consuming.
## Note: The time of simplification steps strongly increases with the number of free symbolic parameters
## 
## For the present example, they are NOT needed!
## 
t0 = time.time() # Time the simplification steps
c = simplify(factor(c))
C = factor(simplify(C))
display(time.time() - t0)
Explanation: Calculating the first and second fundamental cumulants
The function getCumulants(model, chords, rates_params) from the cumulants library takes a dictionary model and a corresponding substitution list rates_params to calculate the fundamental cumulants specified by the fundamental chords given in chords.
While some consistency checks are given, you have to be sure that removing the edges contained in chords from the graph of the model yields a spanning tree, see Ref. [1].
Here, we chose the edges $(0,1)$ and $(0,2)$ as fundamental chords. The corresponding fundamental cycles are then the outer circuit and the lower half circuit (both in clockwise direction):
```
0 –> 1
|    |   (0,1)
3 –– 2
0
| \
|  \     (0,2)
3 – 2
```
End of explanation
## Calculate cycle affinities:
# First fundamental cycle (for chord $(0,1)$): (0,1,2,3):
aff0 = simplify((log((w01*w12*w23*w30)/(w10*w21*w32*w03))).subs(rates_parametrized))
# Second fundamental cycle (for chord $(0,2)$ (0,2,3):
aff1 = simplify((log((w02*w23*w30)/(w03*w32*w20))).subs(rates_parametrized))
# affinities should be $4f$ and $2f+g$, respectively:
display((aff0,aff1))
## Define expressions for quantities of interest
topc = c[0] # average current through top edge
cenc = c[1] # average current through center edge
ep = simplify(c[0]*aff0 + c[1]*aff1) # entropy production
topd = onehalf*C[0,0] # two times variance yields diffusion constants
cend = onehalf*C[1,1]
cov = C[0,1] #co-variance of both currents
res = 2*c[0].diff(f)/C[0,0] # response of top current to the driving affinity $f$ divided by top diffusion constant
#res = simplify(res) 
Explanation: Observables
The fundamental cumulants regard the counting statistics along the fundamental chords. As shown in Ref. [1], this is sufficient in order to determine the counting statistics of arbitrary observables.
In the following we consider three observables:
The current flowing in the top edge (fundamental chord  $(0,1)$, i.e. c[0]).
The current through the center edge (fundamental chord  $(0,2)$, i.e. c[1]). 
The entropy production of the whole network
The entropy production is an observable that takes the value $\log{\frac{w_{\rightarrow}}{w_{\leftarrow}}}$ on each edge, where $w_{\rightarrow}$ and $w_{\leftarrow}$ are the corresponding forward and backward rates. The generalized Schnakenberg decomposition (cf. Ref [1]) ensures that for the calculation of the entropy statistics, we only need the affinities of the fundamental cycles.
End of explanation
# Show some analytical expressions for...
# ...top and center average currents
display('Average currents')
display(topc)
display(cenc)
# ...entropy production
display('Entropy production')
display(ep)
## ...diffusion constants and response (WARNING: these are longish expressions...)
#display('Diffusion constants')
#display(topd)
#display(cend)
#display('Normalized response')
#display(res)
Explanation: Exploring the parameter dependence
Our model depends on three symbolic parameters: $(f,g,b)$
End of explanation
## Lambdify all SymPy expressions into NumPy Expressions
topcL, cencL, epL, topdL, cendL, covL, resL\
    = [ lambdify( (f,g,b), N(thing), "numpy" )\
        for thing in ( topc, cenc, ep, topd, cend, cov, res\
        ) ]
## Prepare 2D plotting range
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show
# prepare the plotting grid
[xmin, xmax, ymin, ymax] = [-10,10,.01,20] # boundaries of the grid
resolution = 400 # plot resolution
plotarea = [xmin, xmax, ymin, ymax] # full plotarea
# prepare the plotting grid for kinesin 6 figures
xpts = linspace(xmin, xmax, resolution)
ypts = linspace(ymin, ymax, resolution)
X, Y  = meshgrid(xpts, ypts)
## General setup for figures
fig_size = (6,5) # in inch
fs = 22 # font size
colormap1 = cm.gist_earth # linear color gradient (blue) for densityplots
colormap2 = cm.coolwarm # color gradient (red-white-blue) for densityplots with highlighted center
# font setup
font = {'family' : 'serif',
        'color'  : 'black',
        'weight' : 'normal',
        'size'   : fs,
        }
ts = 16 # tick+contour label size
figdir = "toymodel/"
## This function takes a lambda function and creates a (logarithmic) 2D plot
##  g: lambda function with *exactly two* arguments
##  t: title string
##  x,y: x and y axis strings
##  logplot: flag whether to plot in logscale or not
##  crop: min/max values
def pplot(g, t, x='', y='', logplot=True, highlight=0, crop=[]):
    fig = figure(figsize=fig_size)
    
    if(logplot):
        G = np.log(np.abs(g(X,Y)))/np.log(10)
        GG = g(X,Y)
        ccmm = colormap1
    else:
        G = g(X,Y)
        GG = G
        ccmm = colormap2
    # the slicing parameter [::-1] reverses the y-axis before plotting
    
    im = imshow( G[::-1], cmap=ccmm, extent=plotarea ) # drawing the function
    if(len(crop)==2):
        im.set_clim(vmin=crop[0], vmax=crop[1])
        
    # adding the contour lines with labels
    cset1    = contour( X,Y, G, arange(-20,20,1),linewidths=1,linestyles="-",colors='black')
    stalling = contour( X,Y, GG,    [highlight], linewidths=3,linestyles="-",colors='white')
    
    # adding the colorbar on the right
    cb = colorbar(im)
    
    # latex fashion title
    title(t, fontdict=font)
    xlabel(x, fontdict=font)
    ylabel(y, fontdict=font)
    
    # Set tick label size
    tick_params(axis='both', which='major', labelsize=ts )
    #savefig(figdir+t+".png")
    return(fig)
## NOTE: Our lambda functions take three parameters (f,g,b).
## For plotting, we need to define a new anonymous lambda function, that takes only two parameters
# Currents and EP depending on two forces f,g with same symmetric contribution on the center edge
pplot(lambda f,g: topcL(f,g,1),"Average current through top edge (log scale)",'$f$','$g$')
pplot(lambda f,g: cencL(f,g,1),"Average current through center edge (log scale)",'$f$','$g$')
pplot(lambda f,g: epL(f,g,1),"Steady state entropy production (log scale)",'$f$','$g$')
# Currents and EP depending on one force f=g while increasing strength of center edge
pplot(lambda f,b: topcL(f,f,b),"Average current through top edge (log scale)",'$f$','$b$')
pplot(lambda f,b: cencL(f,f,b),"Average current through center edge (log scale)",'$f$','$b$')
pplot(lambda f,b: epL(f,f,b),"Steady state entropy production (log scale)",'$f$','$b$')
print("Done")
# Diffusion constant and normalized response depending on two forces f,g with same symmetric contribution on the center edge
pplot(lambda f,g: topdL(f,g,1),"Diffusion top edge (log scale)",'$f$','$g$')
pplot(lambda f,g: cendL(f,g,1),"Diffusion center edge (log scale)",'$f$','$g$')
pplot(lambda f,g: resL(f,g,1),"Normalized $f$-response in top edge",'$f$','$g$', False)
# Diffusion constant and normalized response  depending on one force f=g while increasing strength of center edge
pplot(lambda f,b: topdL(f,f,b),"Average current through top edge (log scale)",'$f$','$b$')
pplot(lambda f,b: cendL(f,f,b),"Average current through center edge (log scale)",'$f$','$b$')
pplot(lambda f,b: resL(f,f,b),"Normalized $f$-response in top edge",'$f$','$b$', False)
print("Done")
Explanation: Plotting
Efficient (and numerically stable) plotting requires us to transform symbolic expressions into lambda functions (using the method lambdify()). The lambda functions are used to evaluate the expressions on a grid.
NOTE: lambdify() has problems with symbolic expressions that have $\LaTeX$ commands as their representations.
If you use something like 
eps = symbols('\\varepsilon')
lambdify() may not work and you have to substitute the $\LaTeX$-style expressions by other (dummy) variables.
End of explanation
<END_TASK>