(* ::Package:: *)

(* ::Section::Closed:: *)
(*Functions to build a 2d grid*)


(* GenerateGrid generates a grid of xbins by ybins grid points spanning [xmin, xmax] x [ymin, ymax]. 
The function returns a list with seven elements:
1. A list of the grid spacing in the x and y directions
2. A list of lists (a matrix that is) of the (x,y) points on the grid. {{{x1, y1}, {x1, y2}, ...}, {{x2, y1}, {x2, y2}, ...}, ...}
3. A flattened and rounded version of those same grid points.  {{x1, y1}, {x1, y2}, ..., {x2, y1}, {x2, y2}, ...}
4. A list of the x coordinates of grid points. {x1, x2, ...}
5. A list of the y coordinates of grid points. {y1, y2, ...}
6. A list of the vertex numbers included in each plaquette.
7. A list of the location of the center of each plaquette.
*)
GenerateGrid[xmax_, xmin_,xbins_, ymax_, ymin_, ybins_]:=Module[{h, points, fpoints, xpoints, ypoints, pvertices, pcenters, 
lenxpts, lenypts, lenxplaquettes, lenyplaquettes, plaquettexyindices, llvertices, ulvertices, lrvertices, urvertices},
h = {(xmax-xmin)/(xbins-1), (ymax-ymin)/(ybins-1)};
(* points is the location of the grid points *)
points = N[Table[{x, y}, {x, xmin, xmax, h[[1]]}, {y, ymin, ymax, h[[2]]}]];
(* Rounded and flattened version of points *)
fpoints = Round[Flatten[points, 1],h[[2]]/400000];
(* List of the x coordinates *)
xpoints = N[Table[x, {x, xmin, xmax, h[[1]]}]];
(* List of the y coordinates *)
ypoints = N[Table[y, {y, ymin, ymax, h[[2]]}]];

lenxpts = Length[xpoints];
lenypts = Length[ypoints];
lenxplaquettes = lenxpts - 1;
lenyplaquettes = lenypts - 1;
plaquettexyindices = Table[{Floor[(i-1)/ lenyplaquettes]+1, Mod[i-1, lenyplaquettes]+1}, {i, 1, lenxplaquettes*lenyplaquettes}];
llvertices = ((#[[1]]-1)*lenypts + #[[2]])&/@plaquettexyindices;
ulvertices = #+1&/@llvertices;
lrvertices = # + lenypts &/@llvertices;
urvertices = #+1&/@lrvertices;
(* List of the vertices associated with each (ordered) plaquette number *)
pvertices = Table[{llvertices[[i]], ulvertices[[i]], lrvertices[[i]], urvertices[[i]]}, {i, 1, lenxplaquettes * lenyplaquettes}];
(* List of the (x,y) coordinate to associate with each plaquette *)
pcenters = Table[{Mean[fpoints[[#, 1]]&/@pvertex], Mean[fpoints[[#,2]]&/@pvertex]}, {pvertex, pvertices}];
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters}
]

(* GenerateGridPeriodicx generates a grid of xbins by ybins grid points spanning [xmin, xmax] x [ymin, ymax] with 
periodic boundary conditions in the x direction. 
The function returns a list with seven elements:
1. A list of the grid spacing in the x and y directions
2. A list of lists (a matrix that is) of the (x,y) points on the grid. {{{x1, y1}, {x1, y2}, ...}, {{x2, y1}, {x2, y2}, ...}, ...}
3. A flattened and rounded version of those same grid points.  {{x1, y1}, {x1, y2}, ..., {x2, y1}, {x2, y2}, ...}
4. A list of the x coordinates of grid points. {x1, x2, ...}
5. A list of the y coordinates of grid points. {y1, y2, ...}
6. A list of the vertex numbers included in each plaquette.
7. A list of the location of the center of each plaquette.
*)
GenerateGridPeriodicx[xmax_, xmin_,xbins_, ymax_, ymin_, ybins_]:=Module[{h, points, fpoints, xpoints, ypoints, pvertices, pcenters, lenxpts, lenypts, lenxplaquettes, lenyplaquettes, plaquettexyindices, llvertices, ulvertices, lrvertices, urvertices},
h = {(xmax-xmin)/(xbins), (ymax-ymin)/(ybins-1)};
(* points is the location of the grid points *)
points = N[Table[{x, y}, {x, xmin, xmax-h[[1]], h[[1]]}, {y, ymin, ymax, h[[2]]}]];
(* Rounded and flattened version of points *)
fpoints = Round[Flatten[points, 1],h[[2]]/400000];
xpoints = N[Table[x, {x, xmin, xmax-h[[1]], h[[1]]}]];
ypoints = N[Table[y, {y, ymin, ymax, h[[2]]}]];

lenxpts = Length[xpoints];
lenypts = Length[ypoints];
lenxplaquettes = lenxpts;
lenyplaquettes = lenypts - 1;
plaquettexyindices = Table[{Floor[(i-1)/ lenyplaquettes]+1, Mod[i-1, lenyplaquettes]+1}, {i, 1, lenxplaquettes*lenyplaquettes}];
llvertices = ((#[[1]]-1)*lenypts + #[[2]])&/@plaquettexyindices;
ulvertices = #+1&/@llvertices;
lrvertices = Mod[# + lenypts,Length[fpoints]] &/@llvertices;
urvertices = #+1&/@lrvertices;
(* List of the vertices associated with each (ordered) plaquette number *)
pvertices = Table[{llvertices[[i]], ulvertices[[i]], lrvertices[[i]], urvertices[[i]]}, {i, 1, lenxplaquettes * lenyplaquettes}];
(* List of the (x,y) coordinate to associate with each plaquette *)
pcenters = Table[{If[pvertex[[3]]<=lenyplaquettes, 0.5, 0]+Mean[fpoints[[#, 1]]&/@pvertex], Mean[fpoints[[#,2]]&/@pvertex]}, {pvertex, pvertices}];
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters}
]


(* This funciton is nasty, but it was a bottleneck for big matrices *)
(* It returns a list of all of the indices for all of the nearest neighbors of i with the second and third arguments giving the 
number of sites in the x and y directions *)
AllNeighbors = Compile[{{i, _Integer}, {numx, _Integer}, {numy, _Integer}},
(* First handle the four corners *)
If[i==1, {2, numy+1},
If[i==numy, {numy - 1, 2 numy},
If[i==numx*numy-numy+1, {numx*numy-2numy+1, numx*numy-numy+2},
If[i==numx*numy, {numx*numy-1, numx*numy-numy},
(* Then handle the left and right edge *)
If[1<i<numy, {i-1, i+1, i+numy},
If[numx*numy-numy+1<i<numx*numy,{i-1, i+1, i-numy},
(* Then handle the top and bottom boundary *)
If[Mod[i,numy]==1, {i-numy, i+1, i+numy},
If[Mod[i, numy]==0, {i-numy, i-1, i+numy},
(* Otherwise you have four neighbors *)
{i-numy, i-1, i+1, i+numy}
]
]
]
]
]
]
]
]
];


(* This funciton is nasty, but it was a bottleneck for big matrices *)
(* It returns a list of all of the indices for all of the nearest neighbors of i with the second and third arguments giving 
the number of sites in the x and y directions *)
AllNeighborsPeriodicx = Compile[{{i, _Integer}, {numx, _Integer}, {numy, _Integer}},
(* First handle the four corners *)
If[i==1, {2, numy+1, numx*numy-numy+1},
If[i==numy, {numy - 1, 2 numy, numx*numy},
If[i==numx*numy-numy+1, {1, numx*numy-2numy+1, numx*numy-numy+2},
If[i==numx*numy, {numy, numx*numy-1, numx*numy-numy},
(* Then handle the left and right edge *)
If[1<i<numy, {i-1, i+1, i+numy, i-numy+numx*numy},
If[numx*numy-numy+1<i<numx*numy,{i-1, i+1, i-numy, i+numy-numx*numy},
(* Then handle the top and bottom boundary *)
If[Mod[i,numy]==1, {i-numy, i+1, i+numy},
If[Mod[i, numy]==0, {i-numy, i-1, i+numy},
(* Otherwise you have four neighbors *)
{i-numy, i-1, i+1, i+numy}
]
]
]
]
]
]
]
]
];


(* 
Convert x and y coordinate to the index number of the list. 
Suppose you want to know the index for the grid point at position {x1, y1}.
You have already constructed the grid, so you have access to fpoints (a flattened, rounded list of all the grid points).
xytoi[{x1, y1}, fpoints] will return the index.
*)
xytoi = Compile[{{point, _Real, 1}, {fpoints, _Real, 2}},Position[fpoints, point][[1,1]]];


(* ::Section::Closed:: *)
(*Functions to build the two species model rate matrix*)


(* These functions generate a rate matrix for two species birth-death problem. I'll write the commments later since they're
not the current focus. *)
GenerateWTwoSpecies[grid_, Wplus_Symbol, Wminus_Symbol, \[Tau]_]:=Module[
{W, i, j, point1, point2, h, tot,  points, fpoints, xpoints, ypoints, edges, pvertices, pcenters, numx, numy},

{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;
numx = Length[xpoints];
numy = Length[ypoints];

(* edges is a list of pairs of vertex numbers that are connected by an edge *)
edges =Flatten[
Table[{i, j}, 
{i,1,Length[fpoints]}, 
{j, AllNeighbors[i, numx, numy]}],
1
];

W = SparseArray[
Table[
(* i and j are indices of the vertices -- they are integers *)
{i,j}=edge;
(* point 1 and point2 are the coordinates of those indices -- they aren't integers *)
point1 = fpoints[[j]];
point2 = fpoints[[i]];
(* This sets the W[[i, j]] matrix element,
which is the rate from index j to index i *)
{i,j}->If[point1 - point2 == {0, h[[2]]},Wminus[point1[[1]], point1[[2]], \[Tau][[1]], \[Tau][[2]]][[2]], 
If[point1 - point2 == {0, -h[[2]]},Wplus[point1[[1]], point1[[2]]][[2]], 
If[point1 - point2 == {h[[1]], 0}, Wminus[point1[[1]], point1[[2]], \[Tau][[1]], \[Tau][[2]]][[1]],
Wplus[point1[[1]], point1[[2]]][[1]]]]],
{edge, edges}
]
];
tot = Total[W, {1}];
Do[
W[[i,i]]=-tot[[i]],
{i, 1, Length[fpoints]}
];
{edges, W}
]


(* ::Section::Closed:: *)
(*Functions to build the two bead model rate matrix*)


(* These functions generate a rate matrix for a model with two coupled beads at different temperatures. 
I'll write the commments later since they're not the current focus. *)

(*Force Function*)
f = Compile[{{x, _Real}, {y, _Real},{k, _Real} },{-2.*k*x+k*y,k*x-2*k*y}];

(* Transition rate to the right *)
Wright = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {k, _Real}},(f[x,y,k][[1]] / (2h)) + (Diffusion / h^2)];
(* Transition rate to the left *)
Wleft = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {k, _Real}},-(f[x,y,k] [[1]]/ (2h)) + (Diffusion / h^2)];
(* Transition rate up *)
Wup = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {k, _Real}},(f[x,y,k] [[2]]/ (2h)) + (Diffusion / h^2)];
(* Transition rate down *)
Wdown = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {k, _Real}},-(f[x,y,k] [[2]]/ (2h)) + (Diffusion / h^2)];


GenerateWTwoBeads[grid_, Diffusion_, k_]:=Module[{W, i, j, point1, point2, h, tot,  points, fpoints, xpoints, ypoints, edges, pvertices, pcenters},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;

(* edges is a list of pairs of vertex numbers that are connected by an edge *)
edges =Flatten[Table[{i, j}, {i,1,Length[fpoints]}, {j, AllNeighbors[i, Length[xpoints], Length[ypoints]]}],1];

W = SparseArray[
Table[
(* i and j are indices of the vertices -- they are integers *)
{i,j}=edge;
(* point 1 and point2 are the coordinates of those indices -- they aren't integers *)
point1 = fpoints[[j]];
point2 = fpoints[[i]];
(* This sets the W[[i, j]] matrix element, which is the rate from index j to index i *)
{i,j}->If[point1 - point2 == {0, h[[2]]},Wdown[point1[[1]], point1[[2]], h[[2]], Diffusion[[2]], k], 
If[point1 - point2 == {0, -h[[2]]},Wup[point1[[1]], point1[[2]], h[[2]], Diffusion[[2]], k],If[point1 - point2 == {h[[1]], 0}, 
Wleft[point1[[1]], point1[[2]], h[[1]], Diffusion[[1]], k],
Wright[point1[[1]], point1[[2]], h[[1]], Diffusion[[1]],k]]]],
{edge, edges}
]
];
tot = Total[W, {1}];
Do[
W[[i,i]]=-tot[[i]],
{i, 1, Length[fpoints]}
];
{edges, W}
]


(* ::Section:: *)
(*Functions to build the four well rate matrix*)


(* These functions generate a rate matrix for a model with a single particle moving on a 4-well potential. 
I'll write the commments later since they're not the current focus. *)

(* The output of this block was copy-pasted into the function for f *)
(* First hardcode the forces *)
(* 
wellcenters = {{1, 1}, {1, -1}, {-1, -1}, {-1, 1}};
\[Sigma] = 1;
U[x_, y_]:=-B Sum[Exp[-((x-xv[[1]])^2+(y-xv[[2]])^2)/\[Sigma]^2], {xv, wellcenters}];
gradU=Grad[U[x,y], {x,y}];
f1 = -gradU;
f2=A (x^2+y^2)Exp[-3Sqrt[x^2+y^2]]{y, -x};
f = f1 + f2 
*)

(* Hardcoded forces *)
f4well = Compile[{{x, _Real}, {y, _Real},{A, _Real}, {B, _Real}},{B (-2 E^(-(-1+x)^2-(-1+y)^2) (-1+x)-2 E^(-(-1+x)^2-(1+y)^2) (-1+x)-2 E^(-(1+x)^2-(-1+y)^2) (1+x)-2 E^(-(1+x)^2-(1+y)^2) (1+x))+A E^(-3 Sqrt[x^2+y^2]) y (x^2+y^2),-A E^(-3 Sqrt[x^2+y^2]) x (x^2+y^2)+B (-2 E^(-(-1+x)^2-(-1+y)^2) (-1+y)-2 E^(-(1+x)^2-(-1+y)^2) (-1+y)-2 E^(-(-1+x)^2-(1+y)^2) (1+y)-2 E^(-(1+x)^2-(1+y)^2) (1+y))}];
(* Transition rate to the right *)
Wright4well = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {B, _Real}},(f4well[x,y,A,B][[1]] / (2h)) + (Diffusion / h^2)];
(* Transition rate to the left *)
Wleft4well = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {B, _Real}},-(f4well[x,y,A,B] [[1]]/ (2h)) + (Diffusion / h^2)];
(* Transition rate up *)
Wup4well = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {B, _Real}},(f4well[x,y,A,B] [[2]]/ (2h)) + (Diffusion / h^2)];
(* Transition rate down *)
Wdown4well = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {B, _Real}},-(f4well[x,y,A,B] [[2]]/ (2h)) + (Diffusion / h^2)];


GenerateW4well[grid_, Diffusion_, A_, B_]:=Module[{W, i, j, point1, point2, h, tot,  points, fpoints, xpoints, ypoints, edges, pvertices, pcenters},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;

(* edges is a list of pairs of vertex numbers that are connected by an edge *)
edges =Flatten[Table[{i, j}, {i,1,Length[fpoints]}, {j, AllNeighbors[i, Length[xpoints], Length[ypoints]]}],1];

W = SparseArray[
Table[
(* i and j are indices of the vertices -- they are integers *)
{i,j}=edge;
(* point 1 and point2 are the coordinates of those indices -- they aren't integers *)
point1 = fpoints[[j]];
point2 = fpoints[[i]];
(* This sets the W[[i, j]] matrix element, which is the rate from index j to index i *)
{i,j}->If[point1 - point2 == {0, h[[2]]},Wdown4well[point1[[1]], point1[[2]], h[[2]], Diffusion[[2]], A, B], 
If[point1 - point2 == {0, -h[[2]]},Wup4well[point1[[1]], point1[[2]], h[[2]], Diffusion[[2]], A, B],If[point1 - point2 == {h[[1]], 0}, 
Wleft4well[point1[[1]], point1[[2]], h[[1]], Diffusion[[1]], A, B],
Wright4well[point1[[1]], point1[[2]], h[[1]], Diffusion[[1]],A, B]]]],
{edge, edges}
]
];
tot = Total[W, {1}];
Do[
W[[i,i]]=-tot[[i]],
{i, 1, Length[fpoints]}
];
{edges, W}
]

GenerateW4wellDiscreteTime[grid_, Diffusion_, A_, B_]:=Module[{W, i, j, point1, point2, h, tot,  points, fpoints, xpoints, ypoints, edges, pvertices, pcenters},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;

(* edges is a list of pairs of vertex numbers that are connected by an edge *)
edges =Flatten[Table[{i, j}, {i,1,Length[fpoints]}, {j, AllNeighbors[i, Length[xpoints], Length[ypoints]]}],1];

W = SparseArray[
Table[
(* i and j are indices of the vertices -- they are integers *)
{i,j}=edge;
(* point 1 and point2 are the coordinates of those indices -- they aren't integers *)
point1 = fpoints[[j]];
point2 = fpoints[[i]];
(* This sets the W[[i, j]] matrix element, which is the rate from index j to index i *)
{i,j}->If[point1 - point2 == {0, h[[2]]},Wdown4well[point1[[1]], point1[[2]], h[[2]], Diffusion[[2]], A, B], 
If[point1 - point2 == {0, -h[[2]]},Wup4well[point1[[1]], point1[[2]], h[[2]], Diffusion[[2]], A, B],If[point1 - point2 == {h[[1]], 0}, 
Wleft4well[point1[[1]], point1[[2]], h[[1]], Diffusion[[1]], A, B],
Wright4well[point1[[1]], point1[[2]], h[[1]], Diffusion[[1]],A, B]]]],
{edge, edges}
]
];
tot = Total[W, {1}];
Do[
W[[i,i]]=-tot[[i]],
{i, 1, Length[fpoints]}
];
{edges, W}
]


(* ::Section::Closed:: *)
(*Functions to build the time-dependent rate matrices*)


(* eratchet is a funciton that computes the energy at point {x, y} and time t for a given A, freq, a1, a2, L. *)
(* You will need to write a function that will first solve for the energy of each site when provided the boundary conditions.
I expect that will look somewhat like this: http://www.public.asu.edu/~hhuang38/pde_slides_numerical_laplace.pdf.
You should then have a function like eratchet which accesses the energy at position {x,y}. *)

(* rachet energy *)
eratchet = Compile[
{{x, _Real}, {y, _Real},{A, _Real}, {freq, _Real}, {t, _Real}, {a1, _Real}, {a2, _Real}, {L, _Real}},
A Sin[2\[Pi] freq t](a1 Sin[2 \[Pi] x / L] + a2 Sin[4 \[Pi] x / L])/(0.5 + y)^2];

Wrightratchet = Compile[
{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real}, {a1, _Real}, {a2, _Real}, {L, _Real}},
(Diffusion / h^2)Exp[-(1/(2 Diffusion))(eratchet[x+h,y,A,freq,t,a1,a2,L] - eratchet[x,y,A,freq,t,a1,a2,L])]];

(* Transition rate to the left *)
Wleftratchet = Compile[
{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real}, {a1, _Real}, {a2, _Real}, {L, _Real}},
(Diffusion / h^2)Exp[-(1/(2 Diffusion))(eratchet[x-h,y,A,freq,t,a1,a2,L] - eratchet[x,y,A,freq,t,a1,a2,L])]];

(* Transition rate up *)
Wupratchet = Compile[
{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real},  {a1, _Real}, {a2, _Real}, {L, _Real}},
(Diffusion / h^2)Exp[-(1/(2 Diffusion))(eratchet[x,y+h,A,freq,t,a1,a2,L] - eratchet[x,y,A,freq,t,a1,a2,L])]];

(* Transition rate down *)
Wdownratchet = Compile[
{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real},  {a1, _Real}, {a2, _Real}, {L, _Real}},
(Diffusion / h^2)Exp[-(1/(2 Diffusion))(eratchet[x,y-h,A,freq,t,a1,a2,L] - eratchet[x,y,A,freq,t,a1,a2,L])]];


(* This block is a different strategy that is deprecated and commented out. *)

(*Clear[x,y,a1,a2,L,A,freq,t]
U[x_, y_, a1_, a2_, L_,A_,freq_, t_]:=A Sin[2\[Pi] freq t](a1 Sin[2 \[Pi] x / L] + a2 Sin[4 \[Pi] x / L])/(0.5 + y)^2;
gradU=Grad[U[x,y, a1, a2, L,A,freq,t], {x,y}];
force = -gradU *)

(* ratchet force came from the above time-dependent potential *)
(*fratchet = Compile[{{x, _Real}, {y, _Real},{A, _Real}, {freq, _Real}, {t, _Real}, {a1, _Real}, {a2, _Real}, {L, _Real}},{-((A ((2 a1 \[Pi] Cos[(2 \[Pi] x)/L])/L+(4 a2 \[Pi] Cos[(4 \[Pi] x)/L])/L) Sin[2 freq \[Pi] t])/(0.5` +y)^2),(2 A Sin[2 freq \[Pi] t] (a1 Sin[(2 \[Pi] x)/L]+a2 Sin[(4 \[Pi] x)/L]))/(0.5` +y)^3}];*)

(* Transition rate to the right *)
(*Wrightratchet = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real}, {a1, _Real}, {a2, _Real}, {L, _Real}},(fratchet[x,y,A,freq,t,a1,a2,L][[1]] / (2h)) + (Diffusion / h^2)];
(* Transition rate to the left *)
Wleftratchet = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real}, {a1, _Real}, {a2, _Real}, {L, _Real}},-(fratchet[x,y,A,freq,t,a1,a2,L][[1]] / (2h)) + (Diffusion / h^2)];
(* Transition rate up *)
Wupratchet = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real},  {a1, _Real}, {a2, _Real}, {L, _Real}},(fratchet[x,y,A,freq,t,a1,a2,L][[2]] / (2h)) + (Diffusion / h^2)];
(* Transition rate down *)
Wdownratchet = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real},  {a1, _Real}, {a2, _Real}, {L, _Real}},-(fratchet[x,y,A,freq,t,a1,a2,L][[2]] / (2h)) + (Diffusion / h^2)];*)

(*Wrightratchet = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real}, {a1, _Real}, {a2, _Real}, {L, _Real}},(eratchet[x,y,A,freq,t,a1,a2,L] - eratchet[x+h,y,A,freq,t,a1,a2,L] / (2h^2)) + (Diffusion / h^2)];
(* Transition rate to the left *)
Wleftratchet = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real}, {a1, _Real}, {a2, _Real}, {L, _Real}},(eratchet[x,y,A,freq,t,a1,a2,L] - eratchet[x-h,y,A,freq,t,a1,a2,L] / (2h^2)) + (Diffusion / h^2)];
(* Transition rate up *)
Wupratchet = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real},  {a1, _Real}, {a2, _Real}, {L, _Real}},(eratchet[x,y,A,freq,t,a1,a2,L] - eratchet[x,y+h,A,freq,t,a1,a2,L] / (2h^2)) + (Diffusion / h^2)];
(* Transition rate down *)
Wdownratchet = Compile[{{x, _Real}, {y, _Real}, {h, _Real}, {Diffusion, _Real}, {A, _Real}, {freq, _Real}, {t, _Real},  {a1, _Real}, {a2, _Real}, {L, _Real}},(eratchet[x,y,A,freq,t,a1,a2,L] - eratchet[x,y-h,A,freq,t,a1,a2,L] / (2h^2)) + (Diffusion / h^2)];*)


(* GenerateWratchet generates the continuous time rate matrices {W_1, W_2, ... W_NStages} for the time-dependent driving.
The function first calculates edges, a list of pairs of vertices that neighbor each other.
It then loops over those edges, determines if the edge corresponds to a step upward, downward, leftward, or rightward, and
computes the rate for the transition accordingly.
This is repeated for each stage of the time-dependent driving, and a list of all the W's is constructed.
The function outputs the list of pairs of vertices defining the edges as well as the list of W's
*)

GenerateWratchet[grid_, Diffusion_, A_, freq_, NStages_, a1_, a2_, L_]:=Module[
{W, i, j, point1, point2, h, tot,  points, fpoints, xpoints, ypoints, edges, period, deltat, Ws, pvertices, pcenters},

{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;
(* edges is a list of pairs of vertex numbers that are connected by an edge *)
edges =Flatten[Table[{i, j}, {i,1,Length[fpoints]}, {j, AllNeighborsPeriodicx[i, Length[xpoints], Length[ypoints]]}],1];

period = 1 / freq;
deltat = period / NStages;

Ws = Table[
W = SparseArray[
Table[
(* i and j are indices of the vertices -- they are integers *)
{i,j}=edge;
(* point 1 and point2 are the coordinates of those indices -- they aren't integers *)
point1 = fpoints[[j]];
point2 = fpoints[[i]];
(* This sets the W[[i, j]] matrix element, which is the rate from index j to index i *)
{i,j}->If[Round[(point1 - point2)/h[[2]]] == {0, 1},
Wdownratchet[point1[[1]], point1[[2]], h[[2]], Diffusion[[2]], A, freq, t, a1, a2,L], 
If[Round[(point1 - point2)/h[[2]]] == {0, -1},
Wupratchet[point1[[1]], point1[[2]], h[[2]], Diffusion[[2]], A, freq, t, a1, a2,L],
If[Round[(Mod[point1 - point2,L])/h[[1]]] == {1, 0}||Round[(Mod[point1+L-point2,L])/h[[1]]]=={1, 0}, 
Wleftratchet[point1[[1]], point1[[2]], h[[1]], Diffusion[[1]], A, freq, t, a1, a2,L],
Wrightratchet[point1[[1]], point1[[2]], h[[1]], Diffusion[[1]],A, freq, t, a1, a2,L]]]],
{edge, edges}
]
];
tot = Total[W, {1}];
Do[
W[[i,i]]=-tot[[i]],
{i, 1, Length[fpoints]}
];
W,
{t, 0, period-deltat, deltat}
];
{edges, Ws}
]

(* This function is exactly like GenerateWratchet except it is for the case that x is not periodically replicated. *)
GenerateWratchetNotPeriodic[grid_, Diffusion_, A_, freq_, NStages_, a1_, a2_, L_]:=Module[
{W, i, j, point1, point2, h, tot,  points, fpoints, xpoints, ypoints, edges, period, deltat, Ws, pvertices, pcenters},

{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;
(* edges is a list of pairs of vertex numbers that are connected by an edge *)
edges =Flatten[Table[{i, j}, {i,1,Length[fpoints]}, {j, AllNeighbors[i, Length[xpoints], Length[ypoints]]}],1];

period = 1 / freq;
deltat = period / NStages;

Ws = Table[
W = SparseArray[
Table[
(* i and j are indices of the vertices -- they are integers *)
{i,j}=edge;
(* point 1 and point2 are the coordinates of those indices -- they aren't integers *)
point1 = fpoints[[j]];
point2 = fpoints[[i]];
(* This sets the W[[i, j]] matrix element, which is the rate from index j to index i *)
{i,j}->If[point1 - point2 == {0, h[[2]]},
Wdownratchet[point1[[1]], point1[[2]], h[[2]], Diffusion[[2]], A, freq, t, a1, a2,L], 
If[point1 - point2 == {0, -h[[2]]},
Wupratchet[point1[[1]], point1[[2]], h[[2]], Diffusion[[2]], A, freq, t, a1, a2,L],
If[point1 - point2 == {h[[1]], 0}, 
Wleftratchet[point1[[1]], point1[[2]], h[[1]], Diffusion[[1]], A, freq, t, a1, a2,L],
Wrightratchet[point1[[1]], point1[[2]], h[[1]], Diffusion[[1]],A, freq, t, a1, a2,L]]]],
{edge, edges}
]
];
tot = Total[W, {1}];
Do[
W[[i,i]]=-tot[[i]],
{i, 1, Length[fpoints]}
];
W,
{t, 0, period-deltat, deltat}
];
{edges, Ws}
]


(* ::Section::Closed:: *)
(*Functions for generating Gillespie trajectories*)


(* Get accumulated transition matrix accW, each row of which has the transition rates up to that point, not counting the diagonal elements of W. This way of storing the rate matrix is useful for the binary searches. *)
accW[Wmatrix_] := Module[
	{WmatTr},
	WmatTr = Transpose[Wmatrix];
	Accumulate /@ (WmatTr - DiagonalMatrix[Diagonal[WmatTr]])
];


(* Copied from the Combinatoria package and then compiled *)
(* Compute
\ the index of the first element in l to exceed k *)
CompiledBinarySearch = Compile[{{l, _Real, 1}, {k, _Real}},
	Module[
		{lo=1,mid,hi=Length[l],el},
			While[lo<=hi,If[(el=Identity[l[[mid=Floor[(lo+hi)/2]]]])===k,Return[mid]];
			If[el>k,hi=mid-1,lo=mid+1]];
			Return[Ceiling[lo-1/2]]
	],
	CompilationTarget->"C",
	RuntimeOptions->"Speed"
]


(* Carry out a step of the kinetic Monte Carlo Dynamics, assuming the rate matrix is dense *)
 (* 
Inputs: 
-accumulatedWmatrix is an accumulated form of the rate matrix, with the accumulation making it easy to do a binary search
-oldsite is the index for the site that was occupied
-\[Xi]step is a random number from the unit interval that determines which site to hop to next
-\[Xi]times is a random number from the unit interval that determines the wait time
-weights is a matrix that determines the contribution from a jump from oldsite to newsite (a d vector)

Outputs:
-newsite is the new site that is entered after the hop
-timeelapsed is how much time elapsed between hops
-op is what the hop contributes to the order parameter
*)
StepOfDynamicsDense = Compile[
{{accumulatedWmatrix, _Real, 2},{oldsite, _Integer}, {\[Xi]step, _Real}, {\[Xi]times, _Real}, {weights, _Real, 2}},

	Module[
		{rates,newsite,op, forwardrate, reversedrate, timeelapsed, outgoingrates},

		outgoingrates = accumulatedWmatrix[[oldsite]][[-1]];

		newsite = CompiledBinarySearch[accumulatedWmatrix[[oldsite]],\[Xi]step * outgoingrates];

(* We passed the accumulated transition matrix rather than the transition matrix, so we compute the transition rates by finding the difference between neighboring matrix elements.  This accumulated transition matrix is like a cumulative distribution function, we're adding up all the transition rates up to that point *)
forwardrate = If[newsite==1, 
accumulatedWmatrix[[oldsite]][[newsite]],accumulatedWmatrix[[oldsite]][[newsite]]-accumulatedWmatrix[[oldsite]][[newsite-1]]
];
op = weights[[newsite]][[oldsite]];

timeelapsed = -Log[\[Xi]times] / outgoingrates;
{newsite, timeelapsed, op}
	],
	RuntimeAttributes->{Listable}, 
	Parallelization->True, 
	CompilationTarget->"C"
];


SampleInitialState[ssdist_]:=CompiledBinarySearch[Accumulate[ssdist], RandomReal[]];


(* Generate a kinetic Monte Carlo trajectory, assuming the rate matrix is dense*)
(*
Inputs:
-accumulatedWmatrix is an accumulated form of the rate matrix, with the accumulation making it easy to do a binary search
-Tobs is the duration of time before the observation ends (or the constant rate matrix is changed to a new rate matrix)
-initialsite is the site at the intitial time
-weights is a matrix that determines the contribution from a jump from oldsite to newsite
-numx is the number of sites in the x direction
-numy is the number of sites in the y direction

Outputs:
-List of the site, time, and accumulated current for a trajectory
*)
KMCTrajectoryDense[accumulatedWmatrix_, Tobs_, initialsite_, weights_, numx_, numy_] := Module[
	{ site,time, op, counter, timeelapsed, opincrement, traj},
	time = 0;
	op = 0;
site = initialsite;
	counter = 1;

(* This generates a trajectory that takes one step too many (the while loop stops only after time has passed Tobs *)
	traj=Reap[Sow[{site, time, op}];
				While[time < Tobs, 
						{site, timeelapsed, opincrement} = 
							StepOfDynamicsDense[accumulatedWmatrix, site,RandomReal[], RandomReal[],weights, numx, numy]; 
						time += timeelapsed; 
						op += opincrement; 
						Sow[{site, time, op}];
						counter = counter +1;
				];
		][[2]][[1]];

(* By only taking the trajectory up to the -2 element, we don't take the extra step *)
	traj[[;;-2]]
]


(* Carry out a step of the kinetic Monte Carlo Dynamics, assuming the rate matrix is sparse *)
 (* 
Inputs: 
-Wmatrix is the sparse rate matrix
-oldsite is the index for the site that was occupied
-\[Xi]step is a random number from the unit interval that determines which site to hop to next
-\[Xi]times is a random number from the unit interval that determines the wait time
-fweights is the contribution to the dynamic order parameter for occupying a state
-gweights is the contribution to the dynamic order parameter for jumping from state to state
-numx is the number of states in the x direction
-numy is the number of states in the y direction

Outputs:
-newsite is the new site that is entered after the hop
-timeelapsed is how much time elapsed between hops
-op is what the hop contributes to the order parameter (the entropy production)
*)
SetAttributes[StepOfDynamicsSparse, HoldFirst];
StepOfDynamicsSparse[Wmatrix_, oldsite_, \[Xi]step_, \[Xi]times_, fweights_Symbol, gweights_Symbol, numx_, numy_]:=Module[
		{neighbors,newsite,op,timeelapsed, outgoingrates, totaloutgoingrates},

neighbors = AllNeighbors[oldsite, numx, numy];
outgoingrates =Wmatrix[[#, oldsite]]&/@neighbors;
totaloutgoingrates = Total[outgoingrates];
 

newsite = neighbors[[CompiledBinarySearch[Accumulate[outgoingrates],\[Xi]step * totaloutgoingrates]]];

op = timeelapsed*fweights[oldsite] + gweights[newsite, oldsite];

timeelapsed = -Log[\[Xi]times] / totaloutgoingrates;
{newsite, timeelapsed, op}
];


(* Generate a kinetic Monte Carlo trajectory, assuming the rate matrix is sparse*)
(*
Inputs:
-Wmatrix is the sparse rate matrix
-Tobs is the duration of time before the observation ends (or the constant rate matrix is changed to a new rate matrix)
-initialsite is the site at the intitial time
-fweights is the contribution to the dynamic order parameter for occupying a state
-gweights is the contribution to the dynamic order parameter for jumping from state to state
-numx is the number of states in the x direction
-numy is the number of states in the y direction

Outputs:
-List of the site, time, and accumulated current for a trajectory
*)
KMCTrajectorySparse[Wmatrix_, Tobs_, initialsite_, fweights_Symbol, gweights_Symbol, numx_, numy_] := Module[
	{ site,time, op, counter, timeelapsed, opincrement, traj},
	time = 0;
	op = 0;
site = initialsite;
	counter = 1;

(* This generates a trajectory that takes one step too many (the while loop stops only after time has passed Tobs *)
	traj=Reap[Sow[{site, time, op}];
				While[time < Tobs, 
						{site, timeelapsed, opincrement} = 
							StepOfDynamicsSparse[Wmatrix, site,RandomReal[], RandomReal[], fweights, gweights, numx, numy]; 
						time += timeelapsed; 
						op += opincrement; 
						Sow[{site, time, op}];
						counter = counter +1;
				];
		][[2]][[1]];

(* By only taking the trajectory up to the -2 element, we don't take the extra step *)
	traj[[;;-2]]
]


(* Carry out a step of the kinetic Monte Carlo Dynamics, assuming the rate matrix is sparse *)
 (* 
Inputs: 
-Wmatrix is the sparse rate matrix
-oldsite is the index for the site that was occupied
-\[Xi]step is a random number from the unit interval that determines which site to hop to next
-\[Xi]times is a random number from the unit interval that determines the wait time
-fweights is the contribution to the dynamic order parameter for occupying a state
-gweights is the contribution to the dynamic order parameter for jumping from state to state
-numx is the number of states in the x direction
-numy is the number of states in the y direction

Outputs:
-newsite is the new site that is entered after the hop
-timeelapsed is how much time elapsed between hops
-op is what the hop contributes to the order parameter (the entropy production)
*)
SetAttributes[StepOfDynamicsSparsePeriodicx, HoldFirst];
StepOfDynamicsSparsePeriodicx[Wmatrix_, oldsite_, \[Xi]step_, \[Xi]times_, fweights_Symbol, gweights_Symbol, numx_, numy_]:=Module[
		{neighbors,newsite,op,timeelapsed, outgoingrates, totaloutgoingrates},

neighbors = AllNeighborsPeriodicx[oldsite, numx, numy];
outgoingrates =Wmatrix[[#, oldsite]]&/@neighbors;
totaloutgoingrates = Total[outgoingrates];
 

newsite = neighbors[[CompiledBinarySearch[Accumulate[outgoingrates],\[Xi]step * totaloutgoingrates]]];

op = timeelapsed*fweights[oldsite] + gweights[newsite, oldsite];

timeelapsed = -Log[\[Xi]times] / totaloutgoingrates;
{newsite, timeelapsed, op}
];


(* Generate a kinetic Monte Carlo trajectory, assuming the rate matrix is sparse*)
(*
Inputs:
-Wmatrix is the sparse rate matrix
-Tobs is the duration of time before the observation ends (or the constant rate matrix is changed to a new rate matrix)
-initialsite is the site at the intitial time
-fweights is the contribution to the dynamic order parameter for occupying a state
-gweights is the contribution to the dynamic order parameter for jumping from state to state
-numx is the number of states in the x direction
-numy is the number of states in the y direction

Outputs:
-List of the site, time, and accumulated current for a trajectory
*)
KMCTrajectorySparsePeriodicx[Wmatrix_, Tobs_, initialsite_, fweights_Symbol, gweights_Symbol, numx_, numy_] := Module[
	{ site,time, op, counter, timeelapsed, opincrement, traj},
	time = 0;
	op = 0;
site = initialsite;
	counter = 1;

(* This generates a trajectory that takes one step too many (the while loop stops only after time has passed Tobs *)
	traj=Reap[Sow[{site, time, op}];
				While[time < Tobs, 
						{site, timeelapsed, opincrement} = 
							StepOfDynamicsSparsePeriodicx[Wmatrix, site,RandomReal[], RandomReal[], fweights, gweights, numx, numy]; 
						time += timeelapsed; 
						op += opincrement; 
						Sow[{site, time, op}];
						counter = counter +1;
				];
		][[2]][[1]];

(* By only taking the trajectory up to the -2 element, we don't take the extra step *)
	traj[[;;-2]]
]


(* ::Section::Closed:: *)
(*Functions for mixed raster/vector graphics and combining graphics objects*)


(* These functions can basically be ignored. It's just dealing with the graphics for when I want to format my publication-quality plots. *)
Clear[undistortedGraphicsColumn];
Options[undistortedGraphicsColumn]=Join[{Spacings->0},Options[Graphics]];
undistortedGraphicsColumn[list_, opts:OptionsPattern[]]:=Module[{sizes=ImageDimensions/@list,width, xsizes, ysizes},width=Max@sizes[[All,1]];
ysizes=sizes[[All,2]]+Join[{0},ConstantArray[OptionValue[Spacings],Length@list-1]];
xsizes = sizes[[All, 1]];
Graphics[Table[Inset[list[[i]],{(width - xsizes[[i]]) / 2,-Plus@@ysizes[[;;i]]},ImageScaled[{0,0}]],{i,Length[list]}],ImageSize->{width,Plus@@ysizes},ImagePadding->None,PlotRange->{{0,width},{-Plus@@ysizes,0}},AspectRatio->Plus@@ysizes/width,PlotRangePadding->None]]

Clear[undistortedGraphicsRow];
Options[undistortedGraphicsRow]=Join[{Spacings->0},Options[Graphics]];
undistortedGraphicsRow[list_, opts:OptionsPattern[]]:=Module[{sizes=ImageDimensions/@list,height, xsizes, ysizes},height=Max@sizes[[All,2]];
xsizes=sizes[[All,1]]+Join[ConstantArray[OptionValue[Spacings],Length@list-1],{0}];
ysizes = sizes[[All, 2]];
Graphics[Table[Inset[list[[i]],{If[i>1, Plus@@xsizes[[;;i-1]],0], (height-ysizes[[i]])/2},ImageScaled[{0,0}]],{i,Length[list]}],ImageSize->{Plus@@xsizes, height},ImagePadding->None,PlotRange->{{0, Plus@@xsizes}, {0, height}},AspectRatio->height/(Plus@@xsizes),PlotRangePadding->None]]


rasterListContourPlot[pList_,opts:OptionsPattern[]]:=Module[{img,cont,contL,plotRangeRule,contourOptions,frameOptions,rangeCoords},contourOptions=Join[FilterRules[{opts},FilterRules[Options[ListContourPlot],Except[{Background,Frame,Axes}]]],{Frame->None,Axes->None}];
contL=ListContourPlot[pList,Evaluate@Apply[Sequence,contourOptions]];
cont=First[Cases[{contL},Graphics[__],Infinity]];
img=Rasterize[Graphics[GraphicsComplex[cont[[1,1]],cont[[1,2,1]]],PlotRangePadding->None,ImagePadding->None,Options[cont,PlotRange]],"Image",ImageSize->With[{size=Total[{2,0} (ImageSize/.{opts})/.{ImageSize->CurrentValue[ImageSize]}]},If[NumericQ[size],size,First[WindowSize/.Options[EvaluationNotebook[]]]]]];
plotRangeRule=FilterRules[AbsoluteOptions[cont],PlotRange];
rangeCoords=Transpose[PlotRange/.plotRangeRule];
frameOptions=Join[FilterRules[{opts},FilterRules[Options[Graphics],Except[{PlotRangeClipping,PlotRange}]]],{plotRangeRule,Frame->True,PlotRangeClipping->True}];
If[Head[contL]===Legended,Legended[#,contL[[2]]],#]&@Show[Graphics[{Inset[Show[SetAlphaChannel[img,"ShadingOpacity"/.{opts}/.{"ShadingOpacity"->1}],AspectRatio->Full],rangeCoords[[1]],{0,0},rangeCoords[[2]]-rangeCoords[[1]]]},PlotRangePadding->None],Graphics[GraphicsComplex[cont[[1,1]],cont[[1,2,2]]]],Evaluate@Apply[Sequence,frameOptions]]]


rasterContourPlot[f_, rx_, ry_, opts : OptionsPattern[]] := Module[
  {
   img,
   cont,
   contL,
   plotRangeRule,
   contourOptions,
   frameOptions,
   rangeCoords
   },
  contourOptions =
   Join[
    FilterRules[{opts},
     FilterRules[
      Options[ContourPlot],
      Except[{Background, Frame, Axes}]
      ]
     ],
    {Frame -> None, Axes -> None}
    ];
  contL = ContourPlot[f, rx, ry,
    Evaluate@Apply[Sequence, contourOptions]
    ];
  cont = First[Cases[{contL}, Graphics[__], Infinity]];
  img = Rasterize[
    Graphics[
     GraphicsComplex[cont[[1, 1]], cont[[1, 2, 1]]],
     PlotRangePadding -> None, ImagePadding -> None,
     Options[cont, PlotRange]
     ], "Image",
    ImageSize -> With[
      {size =
        Total[{2, 0} (ImageSize /. {opts}) /.
          {ImageSize -> 
            CurrentValue[ImageSize]}]},
      If[NumericQ[size],
       size,
       First[WindowSize /. Options[EvaluationNotebook[]]]
       ]
      ]
    ];
  plotRangeRule = FilterRules[Quiet@AbsoluteOptions[cont], PlotRange];
  rangeCoords = Transpose[PlotRange /. plotRangeRule];
  frameOptions = Join[
    FilterRules[{opts},
     FilterRules[Options[Graphics],
      Except[{PlotRangeClipping, PlotRange}]
      ]
     ],
    {plotRangeRule, Frame -> True, PlotRangeClipping -> True}
    ];
  If[Head[contL] === Legended, Legended[#, contL[[2]]], #] &@
  Show[
   Graphics[
    { 
     Inset[
      Show[
       SetAlphaChannel[img,
        "ShadingOpacity" /. {opts} /. {"ShadingOpacity" -> 1}
        ], AspectRatio -> Full
       ],
      rangeCoords[[1]], {0, 0}, rangeCoords[[2]] - rangeCoords[[1]]
      ]
     },
    PlotRangePadding -> None
    ],
   Graphics[
    GraphicsComplex[cont[[1, 1]], cont[[1, 2, 2]]],
    PlotRangePadding -> None,
    ImagePadding -> None
    ],
   Evaluate@Apply[Sequence, frameOptions]
   ]
  ]


frameticks[xmin_, xmax_, xmajor_, xminor_, ymin_, ymax_, ymajor_, yminor_, majorlen_, minorlen_]:=
{Flatten[{Table[{i,MaTeX[i, FontSize->16],{majorlen,0}},{i,xmin,xmax,xmajor}],Table[{i,"",{minorlen,0}},{i,xmin,xmax,xminor}]},1],Flatten[{Table[{j,MaTeX[j, FontSize->16],{majorlen,0}},{j,ymin,ymax,ymajor}],Table[{j,"",{minorlen,0}},{j,ymin,ymax,yminor}]},1],Flatten[{Table[{i,"",{majorlen,0}},{i,xmin,xmax,xmajor}],Table[{i,"",{minorlen,0}},{i,xmin,xmax,xminor}]},1],Flatten[{Table[{j,"",{majorlen,0}},{j,ymin,ymax,ymajor}],Table[{j,"",{minorlen,0}},{j,ymin,ymax,yminor}]},1]};


(* ::Section::Closed:: *)
(*Functions for computing steady-state properties from a rate matrix*)


(* SteadyState accepts a rate matrix as an argument and returns the eigenvector associated to the top eigenvalue (which should be 0).
The function also normalizes that eigenvector so that the sum of the elements is one, in which case the vector is the 
steady-state probability of occupying each site.
The function returns a list {ev, evseed}, where ev is the normalized eigenvector and evseed is the unnormalized eigenvector which can
be retained to use as a seed for future eigenvalue calculations.
*)

SetAttributes[SteadyState, HoldFirst];
SteadyState[W_]:=Module[{evseed, vertexdegrees,keepindices, ev, lev},
evseed = Eigenvectors[W, 1, Method->{"Arnoldi", "Criteria"->"RealPart", "MaxIterations"->250000}][[1]];
(* This part is not really necessary because the graph is already fully connected *)
vertexdegrees = VertexDegree[AdjacencyGraph[HeavisideTheta[Sign[W]-0.5]]];
keepindices =Reap[Do[If[vertexdegrees[[i]]!= 0, Sow[i]], {i, 1, Length[vertexdegrees]}]][[2]][[1]];W= W[[keepindices, keepindices]];

ev = evseed;
lev = Table[1, {i, 1, Length[W]}];
lev = lev *Length[lev]/ Total[lev];
(* ev is the vector of steady state probabilities *)
ev = Re[evseed / (lev.evseed)];
{ev, evseed}
]


(* thermodynamicforces returns a sparse array whose {i,j} element is the thermodynamic force associated to the edge 
from vertex j to vertex i. *)

SetAttributes[thermodynamicforces, HoldFirst];
thermodynamicforces[W_, pss_, edges_] := Module[{i,j},
SparseArray[
Reap[
Do[
{i, j} = edge;
Sow[{j,i}->Re[Log[(W[[j,i]]pss[[i]])/(pss[[j]]W[[i,j]])]]];,{edge, edges}
]
][[2,1]]
]
]


(* currents returns a sparse array whose {i,j} element is the current through the edge from j to i. *)
SetAttributes[currents, HoldFirst];
currents[W_, pss_, edges_]:=Module[{i,j},
SparseArray[
Reap[
Do[
{i,j} = edge;
Sow[{j,i}->Chop[(W[[j,i]]pss[[i]])-(pss[[j]]W[[i,j]])]];,
{edge,edges}
]
][[2,1]]
]
]


(* EPdvec returns a sparse array whose {i,j} element is the entropy production (EP) "dvector", which is essentially the thermodynamic
force without the p terms *)
SetAttributes[EPdvec, HoldFirst];
EPdvec[W_, edges_]:= Module[{i,j},
SparseArray[
Table[
{i, j} = edge;
{j, i}-> Log[W[[j,i]] / W[[i,j]]]
,{edge, edges}
]
]
];


(* Consider a dynamic order parameter (as in Hugo's paper) which uses a function f to weight the occupancy of each state and 
a function g to weight the jumps between states *)
gvec[g_Symbol, edges_]:= Module[{i,j},
SparseArray[
Table[
{i, j} = edge;
{j, i}-> g[j,i]
,{edge, edges}
]
]
];

fvec[f_Symbol, NSites_]:= Module[{i},
SparseArray[
Table[
{i,i}->f[i]
,{i, 1, NSites}
]
]
];


(* Compute the top eigenvalue and top eigenvector with the Arnoldi method of a tilted rate matrix.
The tilt is performed with a dynamic order parameter that weights both occupancies (through f) and jumps (through g).*)
SetAttributes[TiltedTopEigenSystem, HoldFirst];
TiltedTopEigenSystem[W_, fv_, gv_, \[Lambda]_, evseed_]:=Module[{Wtilt, eigenval, eigenvec},
Wtilt = N[W Exp[-\[Lambda] gv] - \[Lambda] fv];
{eigenval, eigenvec} = Eigensystem[Wtilt, 1, Method->{"Arnoldi", "Criteria"->"RealPart", "StartingVector"->ev,"MaxIterations"->250000}];
{First[eigenval], First[eigenvec]}
]


(* Compute the mean and variance of the dynamical order parameter from the maximum eigenvalue of a tilted matrix *)
SetAttributes[MeanVar, HoldFirst];
Options[MeanVar] = {\[Lambda]->0.0001};
MeanVar[W_, fv_, gv_, evseed_, OptionsPattern[]]:=Module[{evl, evr},
evl = TiltedTopEigenSystem[W, fv, gv, -OptionValue[\[Lambda]], evseed][[1]];
evr = TiltedTopEigenSystem[W, fv, gv, OptionValue[\[Lambda]], evseed][[1]];
{-(evr - evl) / (2 OptionValue[\[Lambda]]), (evr + evl) / (OptionValue[\[Lambda]]^2)}
]


(* Compute the scaled cumulant generating function for the dynamic order parameter *)
(* I assume the eigenvector seed is the eigenvector with \[Lambda] = 0 *)

SetAttributes[scgf, HoldFirst];
Options[scgf] = {\[Lambda]min->-0.0001, \[Lambda]max->0.0001, \[Lambda]step->0.0001};
scgf[W_, fv_, gv_, evseed_, OptionsPattern[]]:=Module[{evec, \[Lambda], \[Lambda]values, ascendingtable, descendingtable},
evec = evseed;

\[Lambda]values = Table[\[Lambda], {\[Lambda], OptionValue[\[Lambda]min], OptionValue[\[Lambda]max], OptionValue[\[Lambda]step]}];
negative\[Lambda]values = Reverse[\[Lambda]values[[;;CompiledBinarySearch[\[Lambda]values, 0]-1]]];
positive\[Lambda]values = \[Lambda]values[[CompiledBinarySearch[\[Lambda]values, 0]-1;;]];

ascending = Table[{eigenval, evec} = TiltedTopEigenSystem[W,fv,gv,\[Lambda],evec]; {\[Lambda], eigenval}, {\[Lambda], positive\[Lambda]values}];

evec = evseed;
descending = Table[{eigenval, evec} = TiltedTopEigenSystem[W,fv,gv,\[Lambda],evec]; {\[Lambda], eigenval}, {\[Lambda], negative\[Lambda]values}];

Join[Reverse[descending], ascending]
]


(* Compute the mean and the variance of the entropy production rate using the maximum eigenvalue of a tilted operator *)
SetAttributes[EPRMeanVar, HoldFirst];
EPRMeanVar[W_, evseed_]:=Module[{Wtiltr, Wtiltl, z, eigenval, eigenvec, epr, varepr, left, right, \[Delta]z},
Wtiltr = W^(1-\[Delta]z) * Transpose[W]^\[Delta]z;
Wtiltl = W^(1+\[Delta]z) * Transpose[W]^-\[Delta]z;
\[Delta]z = 0.001;
{eigenval, eigenvec} = 
Chop[
Eigensystem[Wtiltr,
1,
Method->{"Arnoldi", "Criteria"->"RealPart", "StartingVector"->evseed,"MaxIterations"->250000}
]
];
right = First[eigenval];

{eigenval, eigenvec} = 
Chop[
Eigensystem[Wtiltl,
1,
Method->{"Arnoldi", "Criteria"->"RealPart", "StartingVector"->evseed,"MaxIterations"->250000}
]
];

left = First[eigenval];

epr = -(right - left) / (2 \[Delta]z);
varepr = (right + left) / (\[Delta]z^2);
{epr, varepr}
]


(* ::Section::Closed:: *)
(*Functions for computing steady-state properties from a transition (discrete-time) matrix*)


(* Use the same steady state function as for the continuous time rate matrix *)

(* The discrete time matrix is assumed to be on a fully connected graph *)
SetAttributes[thermodynamicforces, HoldFirst];
thermodynamicforcesT[T_, ev_] := Table[Re[Log[(T[[j,i]]ev[[i]])/(ev[[j]]T[[i,j]])]],
{j, 1, Length[T]},
{i, 1, Length[T]}
]


(* The discrete time matrix is assumed to be on a fully connected graph *)
currentsT[T_, ev_]:=Table[(T[[j,i]]ev[[i]])-(ev[[j]]T[[i,j]]),
{j, 1, Length[T]},
{i, 1, Length[T]}
]


EPdvecT[T_]:= Table[Log[T[[j,i]]/T[[i,j]]],
{j, 1, Length[T]},
{i, 1, Length[T]}
]


(* Consider a dynamic order parameter (as in Hugo's paper) which uses a function f to weight the occupancy of each state and a function g to weight the jumps between states *)
gvecT[g_Symbol, NSites_]:= Table[g[j,i],
{j, 1, NSites},
{i, 1, NSites}
]

fvecT[f_Symbol, NSites_]:= Table[If[i==j,f[i],0],
{j,1,NSites},
{i,1,NSites}
]


(* ::Section::Closed:: *)
(*Functions for computing density, local dissipation, current, force in a diffusive limit assuming a sufficiently dense grid*)


(* Compute a list {x, y, steadystate(x,y)/product of grid spacings} for every point {x,y} on the grid.
The steady-state density of the site is scaled by the product of the grid spacings since this scaling will return the continuous
steady-state density in the limit that grid spacing is getting infinitesimally small. *)
ComputeDensity[grid_,pss_]:=Module[
{h, points, fpoints, xpoints, ypoints, sqrtpts, verticies, \[Rho], densities, t, val, xcoord, ycoord, pvertices, pcenters},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;
sqrtpts = Sqrt[Length[fpoints]];
verticies=Table[{i, i}, {i,1,sqrtpts*sqrtpts}];

\[Rho] = ArrayReshape[pss, {sqrtpts, sqrtpts}];

densities= Table[
t={Floor[(#-1)/ sqrtpts]+1, Mod[#-1, sqrtpts]+1}&/@vertex;

val=\[Rho][[t[[1,1]],t[[1,2]]]];

xcoord = points[[1,1,1]] + h[[1]] * t[[1,1]]-h[[1]];
ycoord = points[[1,1,2]] + h [[2]]* t[[1,2]]-h[[2]];
{xcoord, ycoord,val / (h[[1]] * h[[2]])},
{vertex, verticies}
];
densities
]


(* This function is like ComputeDensity but for the case that x is a periodically replicated direction. *)
ComputeDensityPeriodicx[grid_,pss_]:=Module[
{h, points, fpoints, xpoints, ypoints, numx, numy, verticies, \[Rho], densities, t, val, xcoord, ycoord, pvertices, pcenters},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;
numx = Length[xpoints];
numy = Length[ypoints];
verticies=Table[{i, i}, {i,1,numx*numy}];

\[Rho] = ArrayReshape[pss, {numx, numy}];

densities= Table[
t={Floor[(#-1)/ numy]+1, Mod[#-1, numy]+1}&/@vertex;

val=\[Rho][[t[[1,1]],t[[1,2]]]];

xcoord = points[[1,1,1]] + h[[1]] * t[[1,1]]-h[[1]];
ycoord = points[[1,1,2]] + h [[2]]* t[[1,2]]-h[[2]];
{xcoord, ycoord,val / (h[[1]] * h[[2]])},
{vertex, verticies}
];
densities = Join[densities,  {numx * h[[1]], 0, 0} + # &/@densities[[1;;numy]]];
densities
]


(* This function returns the x and y coordinates for the plaquette (taken to be the center of the plaquette) along with
the average of the dissipation on all four edges of the plaquette. *)
ComputeLocalDissipation[grid_, currents_, F_]:=Module[
{h, points, fpoints, xpoints, ypoints, plaquettedissipations, val, xcoord, ycoord, pvertices, pcenters, ll, ul, lr, ur},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;
plaquettedissipations = Table[
xcoord = pcenters[[i,1]];
ycoord = pcenters[[i,2]];
{ll, ul, lr, ur} = pvertices[[i]];
val = currents[[ll, lr]]*F[[ll, lr]] + currents[[ll, ul]]*F[[ll, ul]] + currents[[ul, ur]]*F[[ul, ur]] + currents[[lr, ur]]*F[[lr, ur]];
{xcoord, ycoord, val/(2*h[[1]]*h[[2]])},
{i, 1, Length[pvertices]}
]
]


(* This function returns the x and y coordinates for the plaquette (taken to be the center of the plaquette) along with
the average current vector computed by assigning the x component of the current to the average of top and bottom edge currents
and the y component of the current to the average of the left and right edge currents.*)

ComputePlaquetteCurrents[grid_, currents_]:=Module[
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters, plaquettecurrents, xcur, ycur, xcoord, ycoord, ll, ul, lr, ur},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;

plaquettecurrents = Table[
xcoord = pcenters[[i,1]];
ycoord = pcenters[[i,2]];
{ll, ul, lr, ur} = pvertices[[i]];
xcur = -(currents[[ll, lr]] + currents[[ul, ur]])/2;
ycur = -(currents[[ll, ul]] + currents[[lr, ur]])/2 ;
{{xcoord, ycoord},{xcur/h[[1]], ycur/h[[2]]}},
{i, 1, Length[pvertices]}
]
]


(* This function returns the x and y coordinates for the plaquette (taken to be the center of the plaquette) along with
the average force vector computed by assigning the x component of the current to the average of top and bottom edge thermodynamic forces
and the y component of the current to the average of the left and right edge thermodynamic forces.*)

ComputePlaquetteForces[grid_, F_]:=Module[
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters, plaquetteforces, xfor, yfor, xcoord, ycoord, ll, ul, lr, ur},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;

plaquetteforces = Table[
xcoord = pcenters[[i,1]];
ycoord = pcenters[[i,2]];
{ll, ul, lr, ur} = pvertices[[i]];
xfor = -(F[[ll, lr]] + F[[ul, ur]])/2;
yfor = -(F[[ll, ul]] + F[[lr, ur]])/2 ;
{{xcoord, ycoord},{xfor/h[[1]], yfor/h[[2]]}},
{i, 1, Length[pvertices]}
]
]


(* ::Section::Closed:: *)
(*Functions for plotting density, local dissipation, current, force in a diffusive limit assuming a sufficiently dense grid*)


(* Don't worry about this function yet. *)
Clear[format]
Options[format] = Options[PaddedForm];
format[x_, sd_]:=Module[{m,e},
{m,e} = MantissaExponent[x];
If[x == 0, "0",
If[-2<e<2,ToString[PaddedForm[x, sd]],ToString[PaddedForm[ToExpression[ToString[m]]*10, sd]]<>"\\times10^{"<>ToString[e-1]<>"}"]]
]


(* Make a fully-formatted plot of the steady-state density *)
Options[PlotDensity] = {ColorBarNumberFormat->{2,1}, ColorBarNumberLabelShift->0, ColorBarLabel->"\\rho^\\pi(\\boldsymbol{x})", 
xlabel->"x_1", ylabel->"x_2"};
PlotDensity[densities_, OptionsPattern[]]:=Module[
{cmapdensity, mind, maxd, ticklocationsd, stepd, cmapdensity2, ymin, ymax, xmin, xmax, xmajor, xminor, ymajor, yminor, 
majorlen, minorlen, ft, localdensities, ticklabelsd, densitycolorbar},
cmapdensity = (Blend[{RGBColor[{255, 247, 243}/255],RGBColor[{145, 0, 63}/255]}, #]&);
mind= Round[Min[Transpose[densities][[-1]]]];
maxd = Max[Transpose[densities][[-1]]];

ticklocationsd = N[FindDivisions[{mind, maxd}, 5]];

mind = ticklocationsd[[1]];
maxd = ticklocationsd[[-1]];
stepd = (maxd - mind) /10.;
cmapdensity2 = cmapdensity[(#-mind) / (10.*stepd)]&;
ymin=densities[[1,2]];
ymax=densities[[-1,2]];
xmin=densities[[1,1]];
xmax=densities[[-1,1]];
xmajor=Round[(xmax-xmin)/4];
If[xmajor == 0, xmajor = (xmax-xmin)/4];
xminor=xmajor / 5;
ymajor=Round[(ymax-ymin)/4];
If[ymajor == 0, ymajor = (ymax-ymin)/4];
yminor=ymajor / 5;
majorlen=0.0;
minorlen=0.0;
localdensities = Quiet[
Show[
rasterListContourPlot[densities,
ColorFunction->cmapdensity2,
PlotRange->Full,
ColorFunctionScaling->False
],
Frame->True,
FrameTicks->frameticks[xmin, xmax, xmajor, xminor, ymin, ymax, ymajor, yminor, majorlen, minorlen],
FrameLabel->{MaTeX[OptionValue[xlabel],FontSize->16],MaTeX[OptionValue[ylabel],FontSize->16]},
ImageSize->300,
AspectRatio->1,
ImagePadding->{{60, 10}, {50, 10}}
]
];

ticklabelsd=MaTeX[format[#,OptionValue[ColorBarNumberFormat]], FontSize->12]&/@ticklocationsd;
densitycolorbar = Graphics[
Table[Inset[ticklabelsd[[i]],
Scaled[{0.67+OptionValue[ColorBarNumberLabelShift],  ((0.873-0.1)/(Length[ticklocationsd]-1))(i-1) + 0.1}]], 
{i, 1, Length[ticklocationsd]}],
ImageSize->{110, 230},
Epilog->Inset[
BarLegend[{cmapdensity2, {mind, maxd}},
LegendMarkerSize->{30, 200},
"Ticks"->ticklocationsd,
"TickLabels"->(""&/@ticklocationsd),
Method->{FrameStyle->Black,AxesStyle->None,TicksStyle->Black}]/.AbsoluteThickness[_]->AbsoluteThickness[2],
Scaled[{0, 0}], Scaled[{0, 0}]
],
Prolog->Text[MaTeX[OptionValue[ColorBarLabel], FontSize->16],Scaled[{0.3, 0.95}]],
AspectRatio->Full
];

undistortedGraphicsRow[{localdensities,undistortedGraphicsColumn[{ densitycolorbar, Graphics[{}, ImageSize->{110, 30}]}]}]
]


(* Make a fully-formatted plot of the local dissipation rate *)
Options[PlotLocalDissipation] = {ColorBarNumberFormat->{2,1}, ColorBarNumberLabelShift->0, 
ColorBarLabel->"\\sigma^\\pi(\\boldsymbol{x})", xlabel->"x_1", ylabel->"x_2"};
PlotLocalDissipation[localdissipation_, grid_, OptionsPattern[]]:=Module[
{cmap, minv, maxv, h, points, fpoints, xpoints, ypoints, ticklocations, stepv, cmap2, xmin, xmax, ymin, ymax, xmajor, xminor, 
ymajor, yminor, majorlen, minorlen, ft, ticklabels, localdiss, dissipationcolorbar, pvertices, pcenters},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;
(* Set up color map *)
cmap = ColorData["M10DefaultDensityGradient"];
minv = Min[Chop[Transpose[localdissipation][[-1]]]];
maxv = Max[Chop[Transpose[localdissipation][[-1]]]];

(* Decide how to break up the scale into 5 nice chunks. Ex: if the color map was going to go from 0 to 0.8 it would be nice to have a tick mark at 0, 0.2, 0.4, 0.6, 0.8 *)
ticklocations = N[FindDivisions[{minv, maxv}, 5]];

minv = ticklocations[[1]];
maxv = ticklocations[[-1]];
stepv = (maxv - minv) /10.;
cmap2 = cmap[(#-minv) / (10.*stepv)]&;

ymin=Min[ypoints];
ymax=Max[ypoints];
xmin=Min[xpoints];
xmax=Max[xpoints];
xmajor=Round[(xmax-xmin)/4];
If[xmajor == 0, xmajor = (xmax-xmin)/4];
xminor=xmajor / 5;
ymajor=Round[(ymax-ymin)/4];
If[ymajor == 0, ymajor = (ymax-ymin)/4];
yminor=ymajor / 5;
majorlen=0.0;
minorlen=0.0;
localdiss = Quiet[
Show[
rasterListContourPlot[localdissipation,
(*PlotLegends\[Rule]True,*)
PlotRange->{{xmin, xmax}, {ymin, ymax}},
ColorFunctionScaling->False,
ColorFunction->cmap2,
ClippingStyle->Automatic
],
Frame->True,
FrameTicks->frameticks[xmin, xmax, xmajor, xminor, ymin, ymax, ymajor, yminor, majorlen, minorlen],
FrameLabel->{MaTeX[OptionValue[xlabel],FontSize->16],MaTeX[OptionValue[ylabel],FontSize->16]},
ImageSize->300,
AspectRatio->1,
ImagePadding->{{60, 10}, {50, 10}}
]
];

ticklabels=MaTeX[format[#,OptionValue[ColorBarNumberFormat]], FontSize->12]&/@ticklocations;
dissipationcolorbar = Graphics[
Table[Inset[ticklabels[[i]],
Scaled[{0.65+OptionValue[ColorBarNumberLabelShift], ((0.873-0.10)/(Length[ticklocations]-1))(i-1) + 0.10}]], 
{i, 1, Length[ticklocations]}],
ImageSize->{110, 230},
Epilog->Inset[
BarLegend[{cmap2, {minv, maxv}},
LegendMarkerSize->{30, 200},
"Ticks"->ticklocations,
"TickLabels"->(""&/@ticklocations),
Method->{FrameStyle->Black,AxesStyle->None,TicksStyle->Black}]/.AbsoluteThickness[_]->AbsoluteThickness[2],
Scaled[{0, 0}], Scaled[{0, 0}]
],
Prolog->Text[MaTeX[OptionValue[ColorBarLabel], FontSize->16],Scaled[{0.3, 0.96}]],
AspectRatio->Full
];

undistortedGraphicsRow[{localdiss,undistortedGraphicsColumn[{ dissipationcolorbar, Graphics[{}, ImageSize->{110, 30}]}]}]
]


(* Make a fully-formatted plot of the local currents *)
Options[PlotLocalCurrents] = {ColorBarNumberFormat->{2,1}, ColorBarNumberLabelShift->0,
ColorBarLabel->"\\left|\\boldsymbol{\\mathcal{J}}^\\pi(\\boldsymbol{x})\\right|", xlabel->"x_1", ylabel->"x_2"};
PlotLocalCurrents[plaquettecurrents_, grid_, OptionsPattern[]]:=Module[
{cmapcurrent, minc, maxc, h, points, fpoints, xpoints, ypoints, ticklocationsc, stepc, cmapcurrent2, xmin, xmax, ymin, ymax, 
xmajor, xminor, ymajor, yminor, majorlen, minorlen, ft, ticklabelsc, localcurrents, currentcolorbar, pvertices, pcenters},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;
cmapcurrent = (Blend[{RGBColor[{255, 245, 235}/255], RGBColor[{204,76,2}/255]},#]&);
minc= Min[Norm[#]&/@(Transpose[plaquettecurrents][[-1]])];
maxc =Max[Norm[#]&/@(Transpose[plaquettecurrents][[-1]])];

ticklocationsc = N[FindDivisions[{minc, maxc}, 5]];

minc = ticklocationsc[[1]];
maxc = ticklocationsc[[-1]];
stepc = (maxc - minc) /10.;
cmapcurrent2 = cmapcurrent[(#-minc) / (10.*stepc)]&;

ymin=Min[ypoints];
ymax=Max[ypoints];
xmin=Min[xpoints];
xmax=Max[xpoints];
If[Max[Transpose[pcenters][[1]]]>xmax, xmax+=h[[1]]];
xmajor=Round[(xmax-xmin)/4];
If[xmajor == 0, xmajor = (xmax-xmin)/4];
xminor=xmajor / 5;
ymajor=Round[(ymax-ymin)/4];
If[ymajor == 0, ymajor = (ymax-ymin)/4];
yminor=ymajor / 5;
majorlen=0.0;
minorlen=0.0;
localcurrents = Show[
ListVectorPlot[plaquettecurrents,
PlotRange->{{xmin, xmax}, {ymin, ymax}},
VectorColorFunction->(cmapcurrent2[#5]&),
VectorColorFunctionScaling->False,
VectorScale->{0.04, Automatic, None}
],
Frame->True,
FrameTicks->frameticks[xmin, xmax, xmajor, xminor, ymin, ymax, ymajor, yminor, majorlen, minorlen],
FrameLabel->{MaTeX[OptionValue[xlabel],FontSize->16],MaTeX[OptionValue[ylabel],FontSize->16]},
ImageSize->300,
AspectRatio->1,
ImagePadding->{{60, 10}, {50, 10}}
];

ticklabelsc=MaTeX[format[#,OptionValue[ColorBarNumberFormat]], FontSize->12]&/@ticklocationsc;
currentcolorbar = Graphics[
Table[Inset[ticklabelsc[[i]],
Scaled[{0.58+OptionValue[ColorBarNumberLabelShift], ((0.873-0.1)/(Length[ticklocationsc]-1))(i-1) + 0.1}]], 
{i, 1, Length[ticklocationsc]}],
ImageSize->{110, 230},
Epilog->Inset[
BarLegend[{cmapcurrent2, {minc, maxc}},
LegendMarkerSize->{30, 200},
"Ticks"->ticklocationsc,
"TickLabels"->(""&/@ticklocationsc),
Method->{FrameStyle->Black,AxesStyle->None,TicksStyle->Black}]/.AbsoluteThickness[_]->AbsoluteThickness[2],
Scaled[{0, 0}], Scaled[{0, 0}]
],
Prolog->Text[MaTeX[OptionValue[ColorBarLabel], FontSize->16],Scaled[{0.3, 0.95}]],
AspectRatio->Full
];

undistortedGraphicsRow[{localcurrents,undistortedGraphicsColumn[{ currentcolorbar, Graphics[{}, ImageSize->{110, 30}]}]}]
]


(* Make a fully-formatted plot of the thermodynamic forces *)
Options[PlotLocalForces] = {ColorBarNumberFormat->{2,1}, ColorBarNumberLabelShift->0,
ColorBarLabel->"\\left|\\boldsymbol{\\mathcal{F}}^\\pi(\\boldsymbol{x})\\right|", xlabel->"x_1", ylabel->"x_2"};
PlotLocalForces[plaquetteforces_, grid_, OptionsPattern[]]:=Module[{cmapforce, minf, maxf, h, points, fpoints, xpoints, ypoints, 
ticklocationsf, stepf, cmapforce2, xmin, xmax, ymin, ymax, xmajor, xminor, ymajor, yminor, majorlen, minorlen, ft, ticklabelsf, 
localforces, forcecolorbar, pvertices, pcenters},
{h, points, fpoints, xpoints, ypoints, pvertices, pcenters} = grid;
cmapforce = (Blend[{RGBColor[{255, 245, 235}/255], RGBColor[{204,76,2}/255]},#]&);
minf= Min[Norm[#]&/@(Transpose[plaquetteforces][[-1]])];
maxf =Max[Norm[#]&/@(Transpose[plaquetteforces][[-1]])];

ticklocationsf = N[FindDivisions[{minf, maxf}, 5]];

minf = ticklocationsf[[1]];
maxf = ticklocationsf[[-1]];
stepf = (maxf - minf) /10.;
cmapforce2 = cmapforce[(#-minf) / (10.*stepf)]&;

ymin=Min[ypoints];
ymax=Max[ypoints];
xmin=Min[xpoints];
xmax=Max[xpoints];
xmajor=Round[(xmax-xmin)/4];
If[xmajor == 0, xmajor = (xmax-xmin)/4];
xminor=xmajor / 5;
ymajor=Round[(ymax-ymin)/4];
If[ymajor == 0, ymajor = (ymax-ymin)/4];
yminor=ymajor / 5;
majorlen=0.0;
minorlen=0.0;
localforces = Show[
ListVectorPlot[plaquetteforces,
PlotRange->{{xmin, xmax}, {ymin, ymax}},
VectorColorFunction->(cmapforce2[#5]&),
VectorColorFunctionScaling->False,
VectorScale->{0.04, Automatic, None}
],
Frame->True,
FrameTicks->frameticks[xmin, xmax, xmajor, xminor, ymin, ymax, ymajor, yminor, majorlen, minorlen],
FrameLabel->{MaTeX[OptionValue[xlabel],FontSize->16],MaTeX[OptionValue[ylabel],FontSize->16]},
ImageSize->300,
AspectRatio->1,
ImagePadding->{{60, 10}, {50, 10}}
];

ticklabelsf=MaTeX[format[#,OptionValue[ColorBarNumberFormat]], FontSize->12]&/@ticklocationsf;
forcecolorbar = Graphics[
Table[Inset[ticklabelsf[[i]],
Scaled[{0.58+OptionValue[ColorBarNumberLabelShift], ((0.873-0.1)/(Length[ticklocationsf]-1))(i-1) + 0.1}]], 
{i, 1, Length[ticklocationsf]}],
ImageSize->{110, 230},
Epilog->Inset[
BarLegend[{cmapforce2, {minf, maxf}},
LegendMarkerSize->{30, 200},
"Ticks"->ticklocationsf,
"TickLabels"->(""&/@ticklocationsf),
Method->{FrameStyle->Black,AxesStyle->None,TicksStyle->Black}]/.AbsoluteThickness[_]->AbsoluteThickness[2],
Scaled[{0, 0}], Scaled[{0, 0}]
],
Prolog->Text[MaTeX[OptionValue[ColorBarLabel], FontSize->16],Scaled[{0.3, 0.95}]],
AspectRatio->Full
];

undistortedGraphicsRow[{localforces,undistortedGraphicsColumn[{ forcecolorbar, Graphics[{}, ImageSize->{110, 30}]}]}]
]


(* ::Section::Closed:: *)
(*Functions to estimate a rate function from sampled data*)


Options[RateFunctionList]={NumberBins->30};
RateFunctionList[data_, OptionsPattern[]]:=Module[{histdata, bincenters, binheights, ratefunctionheights},
histdata=HistogramList[data, OptionValue[NumberBins], "PDF"];
bincenters = (histdata[[1]][[;;-2]] + histdata[[1]][[2;;]])/2;
binheights = histdata[[2]];
ratefunctionheights = -Log[binheights]/Tobs;
ratefunctionheights -= Min[ratefunctionheights];
Transpose[{bincenters, ratefunctionheights}]
]


(* Numerically Legendre transform a cumulant generating function into a 
 * rate function *)
RateFunction[cgf_]:=Module[
{},
	lt = Table[
			slope=(cgf[[i+1]][[2]]-cgf[[i-1]][[2]]) / 
					(cgf[[i+1]][[1]]-cgf[[i-1]][[1]]);
			{-slope, cgf[[i]][[1]]slope - cgf[[i]][[2]]},
			{i, 2, Length[cgf]-1}
		];
lt
]


(* ::Section::Closed:: *)
(*Temporal Coarse Graining (in progress)*)


SetAttributes[CGEPRMeanVar, HoldFirst];
CGEPRMeanVar[W_, \[CapitalDelta]t_, evseed_]:=Module[{Wtiltr, Wtiltl, W\[CapitalDelta]t, z, eigenval, eigenvec, epr, varepr, left, right, \[Delta]z},
W\[CapitalDelta]t = MatrixExp[W \[CapitalDelta]t];
Wtiltr = W\[CapitalDelta]t^(1-\[Delta]z)*Transpose[W\[CapitalDelta]t]^\[Delta]z;(*MatrixExp[W^(1-\[Delta]z) * Transpose[W]^\[Delta]z \[CapitalDelta]t];*)
Wtiltl = W\[CapitalDelta]t^(1+\[Delta]z)*Transpose[W\[CapitalDelta]t]^-\[Delta]z;(*MatrixExp[W^(1+\[Delta]z) * Transpose[W]^-\[Delta]z \[CapitalDelta]t];*)
\[Delta]z = 0.001;
{eigenval, eigenvec} = 
Chop[
Eigensystem[Wtiltr,
1,
Method->{"Arnoldi", "Criteria"->"RealPart", "StartingVector"->evseed,"MaxIterations"->250000}
]
];
right = Log[First[eigenval]]/\[CapitalDelta]t;

{eigenval, eigenvec} = 
Chop[
Eigensystem[Wtiltl,
1,
Method->{"Arnoldi", "Criteria"->"RealPart", "StartingVector"->evseed,"MaxIterations"->250000}
]
];

left = Log[First[eigenval]]/\[CapitalDelta]t;

epr = -(right - left) / (2 \[Delta]z);
varepr = (right + left) / (\[Delta]z^2);
{epr, varepr}
]


(* Numerically Legendre transform a cumulant generating function into a 
 * rate function *)
RateFunction[cgf_]:=Module[
{},
	lt = Table[
			slope=(cgf[[i+1]][[2]]-cgf[[i-1]][[2]]) / 
					(cgf[[i+1]][[1]]-cgf[[i-1]][[1]]);
			{-slope, cgf[[i]][[1]]slope - cgf[[i]][[2]]},
			{i, 2, Length[cgf]-1}
		];
lt
]
