(***  
Package to carry out probabilistic cellular automata 
simulations on arbitrary lattices with Mathematica.

Started sometime in '98 
versions 1,2, the mist of time...
version 3, 23-12-1999, changed tracking function
version 4, ..-01-2000, changed event determining system
version 5, 14-01-2000, re-changed tracking function
                       introduced posibility to refer to
                       p[{state}], p[{state1,state2}] and
                       q[state1,state2] in rate specs.
version 6, 17-01-2000, simplified lattice defs
version 7, 20-01-2000, changed ref to lattice stats
version 8, 21-01-2000, changed treatment of history variables
version 9, 21-01-2000, again changed treatment of history variables
version 9.1 06-06-2000, minor change in calc of event rates
version 10   2-7-2000,   changed history variables again  
version 10.1 06-12-2000 disabled lattice checking
version 10.2 05-03-2001 changed logging of the state of the lattice
                        (no longer BIG-G-G-G files!)
version 10.3 08-03-2001 slight change in location of saved files
version 11 31-03-2001   many changes to integrate with the SimPCA package
version 11.2 30-04-2001 added export of network to file readable by SimPCA
version 12 2-05-2001    hopefully removed ambiguity in model def &*^%!
version 13 2-05-2001    hopefully removed ambiguity in model def
version 14 28-06-2001   added tools to create movies
version 15 1-9-2001     added history[x] function
version 16 25-10-2001   changed ---again--- history lookup scheme
version 17 02-11-2001   small adjustments to make it work under Math.3
version 18 02-11-2001   changed interface with SimPCA
version 19 02-11-2001   experimental
version 20 28-02-2002   added Moore neighbourhoods, minor additions
version 21 3-3-2003     minor changes to file names to support unix mathematica
version 22 4-3-2003     include global density events
version 23 1-4-2003     MAJOR change in how models are defined:
                           all events in eventList that fit are added
                           (before it was just the first that fits)
version 24 29-5-2003    attempt to revive the Mathematica PCA code
version 25 2-3-2004		PCA no longer tests existence of previous run output
version 26 7-3-2004		simpca now requires an input file with the initial
                        state of the lattics (INCOMPATIBLE with previous versions)
version 27 29-4-2004    - addition of distribution-generating routines 
                          (experimental)
                        - gives warning if rates depend nonlinearly on environment
                          (capacity to deal with this not currently incorporated)
version 27.2 23-6-2004  small bug fixes
version 28 8-2-2005     allow transfer of arbitrary (non-linear) rate expressions
                        to simPCA.v29
***)


PCAversion = "PCA.v28 February 2, 2005"



If[AtomQ[omega],omega={}]

stateQ[x_] := MemberQ[omega,x]


displayCols = 12
displayRows = 12 

logEvents = False
logLattice = True

plotTimeSwitch=True



(*** Setting up the lattice ***)

createLattice::usage =
"createLattice[type,parameters] sets up the memory structure that represents the lattice"

pos::usage = "pos[site] is the real spatial location of site"
loc::usage = "loc[site] is the relative location of the site on a square grid"




layoutSquares[rows_,cols_] :=
   (noOfSites = rows cols;
    site[{xpos_,ypos_}] := site[{xpos,ypos}] =
       site[Mod[xpos,cols] + cols Mod[ypos,rows]];
    loc[site[x_]] := loc[site[x]] =
       {Mod[x, cols], Quotient[x, cols]};
    loc[x,s_site] := loc[x,s] = loc[s][[1]];
    loc[y,s_site] := loc[y,s] = loc[s][[2]];
    findSite[s_, direction_] := site[loc[s] + direction];)

layoutSquares::usage = "Auxiliary function that defines relative positions of sites"


latticeDefs[square,rows_,cols_] :=
   (layoutSquares[rows,cols];
 	$nbh[s_] := $nbh[s] =
       {edge[s, findSite[s, {0, -1}]],
        edge[s, findSite[s, {-1, 0}]],
        edge[s, findSite[s, {+1, 0}]],
        edge[s, findSite[s, {0, +1}]]};
    findSite[s_, direction_] := site[loc[s] + direction];
    pos[s_site] := pos[s] = loc[s];)
 
 latticeDefs[vonNeumann,rows_,cols_] := latticeDefs[square,rows,cols]
latticeDefs[moore,rows_,cols_] :=
   (layoutSquares[rows,cols];
 	$nbh[s_] := $nbh[s] =
       {edge[s, findSite[s, { 0, -1}]],
        edge[s, findSite[s, { 1, -1}]],
        edge[s, findSite[s, { 1,  0}]],
        edge[s, findSite[s, { 1,  1}]],
        edge[s, findSite[s, { 0,  1}]],
        edge[s, findSite[s, {-1,  1}]],
        edge[s, findSite[s, {-1,  0}]],
        edge[s, findSite[s, {-1, -1}]]};
    findSite[s_, direction_] := site[loc[s] + direction];
    pos[s_site] := pos[s] = loc[s];)
    
latticeDefs[triangular,rows_,cols_] :=
   (layoutSquares[rows,cols];
    $nbh[s_] := $nbh[s] =
       If[EvenQ[loc[y,s]],
                  {edge[s, findSite[s, {0,-1}]],
                   edge[s, findSite[s, {1,-1}]],
                   edge[s, findSite[s, {-1,0}]],
                   edge[s, findSite[s, {1,0}]],
                   edge[s, findSite[s, {0,1}]],
                   edge[s, findSite[s, {1,1}]]},
                  {edge[s, findSite[s, {0,-1}]],
                   edge[s, findSite[s, {-1,-1}]],
                   edge[s, findSite[s, {-1,0}]],
                   edge[s, findSite[s, {1,0}]],
                   edge[s, findSite[s, {0,1}]],
                   edge[s, findSite[s, {-1,1}]]}];
     pos[s_site] := pos[s] =
       If[EvenQ[loc[y,s]],
            {loc[x,s]+0.5,
             0.5 Sqrt[3] loc[y,s]},
            {loc[x,s],
             0.5 Sqrt[3] loc[y,s]}];
     )

trimap[{x_, y_}] :=
  If[EvenQ[y],
    {3 Quotient[x, 2] + Mod[x, 2] + 0.5 ,
      Sqrt[3]/2 y},
    {3 Quotient[x, 2] + 2 Mod[x, 2] ,
      Sqrt[3]/2 y}]
latticeDefs[hexagonal,rows_,cols_] :=
   (layoutSquares[rows,cols];
    $nbh[s_] := $nbh[s] =
       {edge[s,findSite[s,{0,-1}]],
        edge[s,findSite[s,{0,+1}]],
        If[EvenQ[loc[x,s]+loc[y,s]],
           edge[s,findSite[s,{+1,0}]],
           edge[s,findSite[s,{-1,0}]]]};
     pos[s_site] := pos[s] = trimap[loc[s]];)


hasDoubleLink[s_site] :=
   (Length[neighbourhood[s]] != Length[Union[neighbourhood[s]]])


edgeAllowed[x_site,y_site] :=
   (x =!= y) && 
   Not[MemberQ[neighbourhood[x],edge[x,y]]] 
 
latticeDefs[randomLattice[n_],rows_,cols_] := latticeDefs[{random,n},rows,cols]

latticeDefs[{random,nNeighbours_}, rows_, cols_] :=
  Module[{available, nSites,act,test},
    noOfSites = rows cols;
    edgesToDo = noOfSites nNeighbours;
    layoutSquares[rows,cols];
    siteList = Table[site[i],{i,0,noOfSites-1}];
    neighbourhood[s_site] := $nbh[s];
    Map[($nbh[#]={})&,siteList];
    available = 
      Apply[Join,Map[Table[#, {i,nNeighbours}]&,siteList]];
    While[available =!= {},
      act = Length[available];
      i1 = Random[Integer, {1, act}];
      i2 = i1;
      test = 1;
      While[(Not[edgeAllowed[available[[i1]],available[[i2]]]])&&
            (test<act),
            (i2 = Random[Integer, {1, Length[available]}];
             ++act)];
      If[test==act, Print[Help]];
      AppendTo[
        $nbh[available[[i1]]],
        edge[available[[i1]], available[[i2]]]];
      AppendTo[
        $nbh[available[[i2]]],
        edge[available[[i2]], available[[i1]]]];
      --edgesToDo;
      available =
        Delete[available, {{i1}, {i2}}]
    ];
   (* removeMultipleEdges[];*)
    pos[s_site] := pos[s] = 
      {Random[Real,Sqrt[noOfSites]],Random[Real,Sqrt[noOfSites]]};
    ]

neighbourhood[s_site] := $nbh[s]
neighbourhood[e_edge] := $enbh[e] =
   {Complement[neighbourhood[e[[1]]],{e}],
    Complement[neighbourhood[e[[2]]],{Reverse[e]}]};

makeLatticeLists[] :=
   (siteList = Table[site[x], {x, 0, noOfSites - 1}];
    neighbourhoodList :=
      TableOverSet[neighbourhood[x], {x, siteList}];
    edgeList = Apply[Join, neighbourhoodList];
    edgeNeighbourhoodList :=
      TableOverSet[neighbourhood[e], {e, edgeList}];
    noOfEdges = Length[edgeList];
    )

averageNoOfNeighbours := noOfEdges/noOfSites

   
createLattice::usage = 
"createLattice[type,rows,columns] creates a lattice of the desired type"
createLattice[latticeType_, r_, k_] :=
  (Clear[$nbh,$enbh,site,pos,loc,edgeList];
   rows = r;
   cols = k;
   latticeDefs[latticeType,rows,cols];
   makeLatticeLists[])

edge[i_Integer] := edgeList[[i]]


(* to randomize regular lattice *)

multipleEdges[] := 
  Length[Union[edgeList]] =!= Length[edgeList]
findMultipleEdge[] :=
  findEdge[First[edgeList], Rest[edgeList]]
findEdge[ed_, edLst_] :=
  If[MemberQ[edLst, ed], ed, findEdge[First[edLst], Rest[edLst]]]
removeMultipleEdges[] :=
  While[multipleEdges[],
        breakEdge[findMultipleEdge[]]]
breakEdge[e_edge] :=
  crossEdge[e,pickRandomEdge[]]


removeOne[l_, e_] :=
  If[MemberQ[l, e],
    Drop[l, Position[l, e][[1]]]]

removeEdge[edge[v1_,v2_]] :=
   ($nbh[v1] =
      removeOne[$nbh[v1],edge[v1,v2]];
    $nbh[v2] =
      removeOne[$nbh[v2],edge[v2,v1]];
    edgeList = 
      removeOne[removeOne[edgeList,edge[v1,v2]],edge[v2,v1]])
    
addEdge[edge[v1_,v2_]] :=
   (AppendTo[$nbh[v1],edge[v1,v2]];
    AppendTo[$nbh[v2],edge[v2,v1]];
    edgeList = Join[edgeList,{edge[v1,v2],edge[v2,v1]}])

crossEdges[edge[v11_,v12_],edge[v21_,v22_]] :=
   If[ (v11 =!= v12 =!= v21 =!= v22) &&
       (Not[MemberQ[neighbourhood[v11],edge[v11,v22]]]) &&
       (Not[MemberQ[neighbourhood[v12],edge[v12,v21]]]),
       (removeEdge[edge[v11,v12]];
        removeEdge[edge[v21,v22]];
        addEdge[edge[v11,v22]];
        addEdge[edge[v12,v21]];)
      ,
       $crossFail]
       
tryToSwapEdge[e_edge,prob_] :=
   If[Random[Real] < prob,
      crossEdges[e,pickRandomEdge[]]]

randomizeLattice::usage = 
"use randomimzeLattice[prob_] to `cross' edges with a random other edge"
randomizeLattice[prob_:1] :=
   (Do[crossEdges[pickRandomEdge[],pickRandomEdge[]],
       {prob noOfEdges/2}];
    makeLatticeLists[];
    (*edgeList = Apply[Union, neighbourhoodList];*))
    
  

(*** utilities to display lattice ***)

reduceToSites[s_site] := s
$reduceToSites[edge[s1_, s2_]] := {s1, s2}
$reduceToSites[s_site] := {s}
$reduceToSites[l1_List,l2_List] := reduceToSites[Union[l1,l2]]
reduceToSites[edge[s1_, s2_]] := reduceToSites[edge[s1, s2]]
reduceToSites[edges_List] :=
    Apply[Union, Map[$reduceToSites, edges]]
expandSites[sl_List] :=
    Apply[Union, Map[neighbourhood, sl]]
expandRegion[s_site] := neighbourhood[s]
expandRegion[e_edge] := Union[neighbourhood[e]]
expandRegion[r_List] := expandSites[reduceToSites[r]]


expandRegion[midpoint_site,radius_] := 
   Nest[expandRegion, midpoint, radius]


(* to do with displaying the lattice *)

plotState[i_,pos_] :=
   {GrayLevel[0.8], Circle[pos,0.2]}

plotSite[x_site] :=
   plotState[sigma[x],pos[x]]

plotLattice[window_:plotRegion] :=
  Graphics[Map[plotSite, window]]

plotTime := 
   Graphics[If[plotTimeSwitch,
            {Text["t = " ~StringJoin~ ToString[N[time[]]],
                    {1,-1},{-1,0}]},{}]]

$rangeExpr := 
   ({{Min[#[[1]]] - 1, Max[#[[1]]] + 1},
     {Min[#[[2]]] - 2, Max[#[[2]]] + 1}} &) @ Transpose[pos /@ siteList]
     
$latticeRange=All
(* but $latticeRange=$rangeExpr is another solution *)

$latticeImageSize={400,400}
latticeImageSize=$latticeImageSize

showLattice[window_:plotRegion] :=
  Show[plotLattice[window],    
       plotTime,
       AspectRatio->Automatic,
       PlotRange:>$latticeRange,
       ImageSize->latticeImageSize,
       DisplayFunction->$DisplayFunction]

plotEdge[e_edge] :=
   Line[{pos[e[[1]]],pos[e[[2]]]}]

plotEdges[eList_:edgeList] :=
   Graphics[Map[plotEdge,eList]]

showEdges[edges_:edgeList] := 
  Show[plotEdges[edges], AspectRatio -> Automatic, 
    DisplayFunction -> $DisplayFunction]

plotRegion := siteList

siteMap[] :=
  Table[
   Table[site[{i,j}],
         {j,0,displayCols-1}],
   {i,0,displayRows-1}]

displayRules := 
  Thread[omega -> 
      Characters[StringTake[".abcdefghijklmnopqrstuvwxyzABC", Length[omega]]]]

displayLattice[] := 
   MatrixForm[Map[sigma, siteMap[], {2}]/.displayRules]

(* picking randomly... *)

pickRandomElement[from_List] :=
  from[[Random[Integer,{1,Length[from]}]]]
pickRandomSite[] :=
  site[Random[Integer, noOfSites-1]]
pickRandomEdge[] :=
  Module[{s=pickRandomSite[]},
         edge[s,pickRandomNeighbour[s]]]
pickRandomNeighbour[s_site] :=
  pickRandomElement[neighbourhood[s]][[2]]

pickFromDist[outcome_List, prob_List] :=
  outcome[[pickIndex[prob]]]

pickIndex[prob_List] :=
  Module[{u = Random[Real]},
    selectIndex[u, 1, prob]]

selectIndex[u_Real, i_Integer, prob_List] :=
  If[u >= First[prob], i,
    selectIndex[u - First[prob], i + 1, Rest[prob]]]  
   

cumulative[dist_] := FoldList[Plus, 0, dist]

pickFromDist[dist_] := 
   pickFromCumDist[cumulative[dist]]

pickFromCumDist[cd_] := 
   Module[{u=Random[Real],i=1},
     While[u>=cd[[i]],i++];
     i-1]
   
  
pickOutcomeFrom[outcome_List, prob_List] :=
  outcome[[pickIndex[prob]]]

pickIndex[prob_List] :=
  Module[{u = Random[Real]},
    selectIndex[u, 1, prob]]

selectIndex[u_Real, i_Integer, prob_List] :=
  If[u <= First[prob], i,
    selectIndex[u - First[prob], i + 1, Rest[prob]]]









(*** things to do with the state of the lattice ***)

omega::usage = "omega represents the set of possible states a site can be in"


state::usage = "state[object] represents the current state of an object (site, 
edge, configuration or set of sites"
state[s_site] := sigma[s]

state[e_edge] := {sigma[e[[1]]], sigma[e[[2]]]}

state[sites_List] := Map[state,sites]



setState[x_site, state_] := sigma[x] = state

setState[e_edge, {state1_, state2_}] :=
  {sigma[e[[1]]] = state1, sigma[e[[2]]] = state2}

config[s_site] := Map[state, neighbourhood[s]]
config[e_edge] := Map[state, neighbourhood[e], {2}]


omegaSquared :=
 Flatten[Outer[List,omega,omega],1]
omegaCubed :=
 Flatten[Outer[List, omega, omega, omega], 2]
      



(*** seeding lattices ***)

seedLatticeRandom[initRule_] :=
  Module[{dist,cdist},
     dist = Map[init,omega]//.initRule;
     cdist = cumulative[dist];
     Map[setState[#,pickState[cdist]]&,siteList];
     countLattice[];
    ]








swapStates[s1_,s2_]:=
  Block[{dummy},
    uncountBothWays/@state[neighbourhood[s1]];
    uncountBothWays/@state[neighbourhood[s2]];
    dummy=state[s1];
    setState[s1,state[s2]];
    setState[s2,dummy];
    countBothWays/@state[neighbourhood[s1]];
    countBothWays/@state[neighbourhood[s2]];
    ]

attemptSwapSites[goal_]:=
  Module[{hb,ha},
    s1=pickRandomSite[];
    s2=s1;
    While[s1\[Equal]s2,s2=pickRandomSite[]];
    hb=calcGoodness[goal];
    ha=hb;
    If[state[s1]=!=state[s2],
      (swapStates[s1,s2];
        ha=calcGoodness[goal];
        If[ha>hb,swapStates[s1,s2]];)
      ];
    ]



maxAdjustIter=10000;adjustPrecision=100;

modifyPairs[goal_]:=
  Module[{},
    For[adjustIter=0,
      (adjustIter<maxAdjustIter)&&((goodness=calcGoodness[goal])>
            adjustPrecision),
      adjustIter++,
      attemptSwapSites[goal]]]

















seedLattice[initRule_] :=
  Module[{dist,cdist},
     dist = Map[init,omega]//.initRule;
     cdist = cumulative[dist];
     Map[setState[#,pickState[cdist]]&,siteList];
     countLattice[];
    ]





(* things to do with the events affecting states *)

clearModel[] := Null

eventList::usage = "List of events affecting pairs in the form of {event[p1,p2]->rate,...} where rate can be specified in terms of
n[state,{x,y}] and n[{x,y},state]. The first event that applies to a given configuration is taken."

simulationRules::usage = "Use to specify numerical values of the parameters"

swap[{x_,y_}] := {y,x}

mirrorEvent[Rule[event[{a_, b_}, {c_, d_}], rate_]] :=
  Rule[event[{b, a}, {d, c}], (rate /.
       {nb[x_,{p1_,p2_}] -> nb[{p2,p1},x],
        nb[{p1_,p2_},x_] -> nb[x,{p2,p1}]})]

mirrorEvent[RuleDelayed[event[{a_, b_}, {c_, d_}], rate_]] :=
  RuleDelayed[event[{b, a}, {d, c}], (rate /.
       {nb[x_,{p1_,p2_}] :> nb[{p2,p1},x],
        nb[{p1_,p2_},x_] :> nb[x,{p2,p1}]})]


       

(* to find the number of neighbours in state in a given config *)
$neighbours[left,cfg_,pair_,state_] :=
  Count[cfg[[1]], {cfg[[1, 1, 1]], state}]
$neighbours[right,cfg_, pair_,state_] :=
  Count[cfg[[2]], {cfg[[2, 1, 1]], state}]




same[x_, y_] := If[x === y, 1, 0]
leftNb[x_, {y_, z_}] := nb[x, {y, z}] + same[x, z]
rightNb[{x_, y_}, z_] := nb[{x, y}, z] + same[x, z]
  
countNeighbours[cfg_] :=
   {nb[state_,{x_,y_}] :> $neighbours[left,cfg,{x,y},state],
    nb[{x_,y_},state_] :> $neighbours[right,cfg,{x,y},state]}
    

(* Things dealing with picking events *)




applyEventListOneWay[p1_?pairQ, p2_?pairQ] :=
  Plus @@ (Replace[event[p1,p2], #] & /@ eventList) /.
             event[_, _] -> 0 /. $focus->p1

applyEventListOtherWay[p1_?pairQ, p2_?pairQ] :=
  Plus @@ (Replace[event[swap[p1],swap[p2]], #]& /@ eventList) /.
             event[_, _] -> 0 /. nb[x_,y_]->nb[y,x] /.$focus->p1

applyEventList[p1_?pairQ, p2_?pairQ] :=
  applyEventListOneWay[p1,p2] + applyEventListOtherWay[p1,p2]


setupEventTable[] :=
  (TableOverSet[
    fullEvt[p1] =
    TableOverSet[
      evt[p1, p2] = 
        applyEventList[p1,p2] //.basicSimulationRules,
      {p2, omegaSquared}],
    {p1, omegaSquared}])
    

determineCfgDependence[] :=
   Map[(cfgDep[#]=!FreeQ[rates[#],nb])&,omegaSquared]

determineGlobDependence[] :=
   Map[(globDep[#]=!FreeQ[rates[#],p])&,omegaSquared]


compactEvtTable::usage = "exclude zero-rate events"
compactEvtTable[] :=
   Map[(zrPos = Position[fullEvt[#],0,1];
       evt[#] = Delete[fullEvt[#],zrPos];
       rates[#] = evt[#] dt/2;
       outcomes[#] = Delete[omegaSquared,zrPos];
       calculated[#] = False;
       )&,
       omegaSquared]




      

eventOccurs[sumRates_,sloppy] :=
   (Random[Real] < sumRates)
   
eventOccurs[sumRates_,precise] :=
   (Random[Real] < 1 - fExp[-sumRates])




calcRates[e_edge,pair_] :=
  Block[{},
    actualRates[pair] = rates[pair];
    If[ cfgDep[pair], 
        actualRates[pair] = actualRates[pair] /. countNeighbours[config[e]] ];
    If[ globDep[pair], 
        actualRates[pair] = actualRates[pair] /. p[x_] -> p[x,current] /. q[x_,y_] -> q[x,y,current] ];
    calculated[pair] = True;
    sumRates[pair] = Plus @@ actualRates[pair];
    If[sumRates[pair] > 0, relRates[pair] = actualRates[pair] / sumRates[pair], 
                           relRates[pair] = {} ];
    cumRelRates[pair] = Rest[FoldList[Plus,0,relRates[pair]]];
  ]

pickEventFor[pair_] :=
  Block[{u},
    u = Random[Real];
    pickEventFrom[u,cumRelRates[pair],outcomes[pair]] ]

pickEventFrom[u_,cr_,oc_] :=
  Block[{i=1},
    While[ u > cr[[i]], ++i];
    oc[[i]] ]



toCalc[pair_] :=  
  (!calculated[pair] || globDep[pair] || cfgDep[pair] ) 

pickEvent[e_edge,pair_] :=
  Module[{},
    If[ toCalc[pair], calcRates[e,pair] ];
    If[ sumRates[pair] > 0, 
           If[ eventOccurs[sumRates[pair],method],
                  pickEventFor[pair],
                  nop ]
        ]]
   



      
tryEvent[e_edge]:=
   Module[{pair},
      pair = state[e];
      pickEvent[e,pair]]




doEdgeEvent[] :=
   Module[{e=pickRandomEdge[],newPair},
      newPair = tryEvent[e];
      If[Not[AtomQ[newPair]], recordEvent[e,newPair]];
      {e,newPair}]









(*** registering events to do the bookkeeping ***)

eventRecord::usage = "gives all the events occurred in the simulation"

initCount[] :=
   (Map[(pairCount[#] = 0) &, omegaSquared];
    Map[(singleCount[#] = 0)&, omega];)

pairList[] := Map[state,edgeList]
singleList[] := Map[state,siteList]

countPair[i_] := ++pairCount[i]
countSingle[i_] := ++singleCount[i]

countBothWays[{i_,j_}] := (++pairCount[{i,j}]; ++pairCount[{j,i}])
uncountBothWays[{i_,j_}] := (--pairCount[{i,j}]; --pairCount[{j,i}])


countLattice::usage = "Counts all pairs on the lattice"
countLattice[] := 
   (initCount[];
    Map[countPair,pairList[]];
    Map[countSingle,singleList[]];)


listPairCounts[] :=
   Map[pairCount,omegaSquared]
listSingleCounts[] :=
   Map[singleCount,omega]
   
   

(*** counting triples on a lattice ***)

initTripleCount[] :=
   Map[(tripleCount[#] = 0)&, omegaCubed]
   
countTriple[i_] := ++tripleCount[i]

leftTriples[cfg_] :=
  Map[Prepend[#, cfg[[1, 1, 1]]] &, cfg[[2]]]
  
rightTriples[cfg_] :=
  Map[Append[#, cfg[[2, 1, 1]]] &, cfg[[1]]]
  
triples[cfg_] :=
  Join[leftTriples[cfg], rightTriples[cfg]]

countTriplesAroundEdge[edge_] :=
  Scan[countTriple, triples[state[neighbourhood[edge]]]]

listTripleCounts[] :=
  Map[tripleCount, omegaCubed]

noOfTriples :=
  Apply[Plus,listTripleCounts[]]

countTriples::usage = "Count triple frequencies on the lattice"
countTriples[] :=
(initTripleCount[];
  Scan[countTriplesAroundEdge, edgeList];
  listTripleCounts[])

tripletQ[x_] := MemberQ[omegaCubed,x]


pTriple[x_?tripletQ] := tripleCount[x]/noOfTriples
pTriple[x_?pairQ] := pairCount[x]/noOfEdges
pTriple[x_?singletQ] := singleCount[x]/noOfSites
qTriple[x_,y_?pairQ] := pTriple[Prepend[y,x]]/(ct[y][t]/noOfEdges)


pairCorrelation[{x_,y_}] :=
  pTriple[{x,y}]/(pTriple[x] pTriple[y])
tripleCorrelation[{x_,y_,z_}] :=
  pTriple[{x,y,z}]/
  (pTriple[x] pTriple[y] pTriple[z] pairCorrelation[{x,y}] pairCorrelation[{y,z}])

pairError[x_, {y_, z_}] :=
  qTriple[x, {y, z}] - q[x, y, t]





(*** bookkeeping ***)

pairQ[x_] := MemberQ[omegaSquared,x]
singletQ[x_] := MemberQ[omega,x]

historyList := 
  Join[{time}, Map[ct, omega], Map[ct, omegaSquared]]

ct[x_][current] := 
  Which[ pairQ[x], pairCount[x],
         singletQ[x], singleCount[x]
       ]


resetHistory[] := 
  Module[{}, 
   Map[(#[history]={})&, historyList];
   lastTimeAdded = -Infinity;
   Null]
  
addToHistory[] := 
  (If[time[] > lastTimeAdded,
     (lastTimeAdded = time[];
      Map[AppendTo[#[history], #[current]] &, historyList]), Print["huh"]];
   Null)

ct[x_][t_?NumberQ] := 
  If[t == time[], ct[x][current], 
                Interpolation[
                  Transpose[{time[history],ct[x][history]}],
                  InterpolationOrder -> 1][t]]
                    
     
     




(*** create a rules file that can be read by simPCA ***)


makeSimPCAparams::nonlin = "rate depends nonlinearly on local or global environment"

tr[p1_,p2_] := applyEventList[p1,p2] //.tempRules//.parameterValues

nstr[x_] := 
  (If[Not[NumberQ[x]],Message[makeSimPCAparams::nonlin]];
   ToString[CForm[N[x]]])
  

stateExpr[x_] :=
  (x /. displayRules) <> " " <> 
    nstr[init[x] //. parameterValues]


baseRate[p1_, p2_] := nstr[tr[p1, p2] /. nb[___] -> 0 /. p[___] -> 0]

(* SimPCA versions <= v.28

leftDep[p1_, p2_, x_] :=
  " " <> nstr[D[tr[p1, p2], nb[x, p1]]]

rightDep[p1_, p2_, x_] :=
    " " <> nstr[D[tr[p1, p2], nb[p1, x]]]
    
globDep[p1_, p2_, x_] :=
    " " <> nstr[D[tr[p1, p2], p[{x}]]]

lefties[p1_, p2_] := StringJoin @@ (leftDep[p1, p2, #] & /@ omega)

righties[p1_, p2_] := StringJoin @@ (rightDep[p1, p2, #] & /@ omega)

localEffects[p1_,p2_] :=
  If[FreeQ[tr[p1, p2], nb],
     " 0",
     " 1" <> lefties[p1,p2] <> righties[p1,p2] ]

globals[p1_,p2_] := StringJoin @@ (globDep[p1, p2, #] & /@ omega)

globalEffects[p1_,p2_] :=
  If[FreeQ[tr[p1, p2], p],
     " 0",
     " 1" <> globals[p1,p2] ]



ruleExpr[p1_, p2_] :=
    (p1 /. displayRules) <> 
      "->" <> (p2 /. displayRules) <> " " <> ToString[baseRate[p1, p2]] <> 
     localEffects[p1,p2] <> globalEffects[p1,p2]
*)





(* This is a hack to convert Mathematica expressions to strings
that are readable by Harald Helfgott's formulc package that is
now integrated into the simpca runtime *)


formulcForm[x_] := ToString[x]
formulcForm[x_Real] := "0"<>ToString[x]<>"0"
formulcForm[\[ExponentialE]^x_] := "exp("<>formulcForm[x]<>")"
formulcForm[E^x_] := "exp("<>formulcForm[x]<>")"
formulcForm[Log[y_]] := "ln("<>formulcForm[y]<>")"
formulcForm[Sin[y_]] := "sin("<>formulcForm[y]<>")"
formulcForm[Cos[y_]] := "cos("<>formulcForm[y]<>")"
formulcForm[Tan[y_]] := "tan("<>formulcForm[y]<>")"
formulcForm[ArcSin[y_]] := "asin("<>formulcForm[y]<>")"
formulcForm[ArcCos[y_]] := "acos("<>formulcForm[y]<>")"
formulcForm[ArcTan[y_]] := "atan("<>formulcForm[y]<>")"
formulcForm[Abs[y_]] := "abs("<>formulcForm[y]<>")"
formulcForm[Sqrt[y_]] := "sqrt("<>formulcForm[y]<>")"
formulcForm[x_^y_] := "("<>formulcForm[x]<>")^("<>formulcForm[y]<>")"
formulcForm[x_ y_] := "("<>formulcForm[x]<>")*("<>formulcForm[y]<>")"
formulcForm[x_ +y_] := formulcForm[x]<>"+"<>formulcForm[y]
formulcForm[p[{x_}]] := "glob("<>ToString[ordOmega[x]-1]<>")"
formulcForm[nb[x_,y_List]] := "lnbr("<>ToString[ordOmega[x]-1]<>")"
formulcForm[nb[y_List,x_]] := "rnbr("<>ToString[ordOmega[x]-1]<>")"

ruleExpr[p1_,p2_] :=
   (p1 /. displayRules) <> 
      "->" <> (p2 /. displayRules) <> " " <> 
      ToString[baseRate[p1, p2]] <> " " <>
      If[FreeQ[tr[p1,p2],nb],"0 ","1 "] <>
      If[FreeQ[tr[p1,p2], p],"0 ","1 "] <>
      formulcForm[N[tr[p1,p2]]]

(* end of formulc-related stuff *)


simPCAparamFn$ = "simpca.in"
simPCAparamFn := simPCAparamFn$ 

simPCAlatFn$ = ".lat"
simPCAlatFn := runName<>simPCAlatFn$
simPCAinitFn$ = ".sig"
simPCAinitFn := runName<>simPCAinitFn$

runName$ = "run"
runName = runName$

makeSimPCAparams[fname_String] := 
  Block[ {ruleFile,lstr},
        ruleFile=
          OpenWrite[fname, 
                    FormatType->OutputForm, PageWidth->Infinity];
         If[StringQ[simPCAlatFn],
            lstr="f "<>simPCAlatFn,
            (*else*)
            lstr = "rnd";
	        If[latticeType == triangular, lstr="tri"]];
	     Write[ruleFile, lstr];
	     Write[ruleFile, rows];
	     Write[ruleFile, cols];
	     Write[ruleFile,N[starttime]];
	     Write[ruleFile,N[stoptime]];
	     Write[ruleFile,N[dt]];
	     Write[ruleFile,N[dtout]];
	     Write[ruleFile,"sloppy"];
	     Write[ruleFile,"i "<>simPCAinitFn];
	     Write[ruleFile,runName];
	     Write[ruleFile, stateExpr[#]] & /@ omega;
	     TableOverSet[
	      TableOverSet[
	        If[tr[p1, p2] =!= 0, Write[ruleFile, ruleExpr[p1, p2]]], {p1, 
	          omegaSquared}], {p2, omegaSquared}];
	     Close[ruleFile];
  ]





(* make a lattice file readable by SimPCA *)

pcaNbr[strm_,edge[x_,site[y_]]]:=
  WriteString[strm," ",y] (* this is improvisation of the worst kind*)

pcaSite[strm_,site[i_]]:=
Block[{},
   WriteString[strm,
      i," ",
      loc[site[i]][[1]]," ",loc[site[i]][[2]]," ",
      pos[site[i]][[1]]," ",pos[site[i]][[2]]," ",
      Length[$nbh[site[i]]]];
   pcaNbr[strm,#]& /@ $nbh[site[i]];
   Write[strm," "]; 
  ]

makeLatticeFile[fname_String] :=
  Block[{strm},
     strm =
           OpenWrite[fname, 
                    FormatType->OutputForm, PageWidth->Infinity];
     pcaSite[strm,#]& /@ siteList;
     Close[strm];
   ]

makeInitFile[fname_String] :=
  Block[{strm},
    strm = OpenWrite[fname, 
                    FormatType->OutputForm, PageWidth->Infinity];
    WriteString[strm,sigma[#]/.displayRules]& /@ siteList;
    Close[strm];
  ]




  
prepExternal[] := 
  Block[ {ruleFile,lstr},
        Print["Creating parameter file "<>simPCAparamFn];
        makeSimPCAparams[simPCAparamFn];
        Print["Creating lattice file "<>simPCAlatFn];
        makeLatticeFile[simPCAlatFn];
        Print["Creating initial state file "<>simPCAinitFn];
        saveSigma[simPCAinitFn];
        If[FileType[runName]=!=Directory,CreateDirectory[runName]];
	    Print["... Ready"];
  ] 



readLattice[run_] :=
  Get[run<>".net"]

readLattice[] :=
  readLattice[runName]

readRun[run_] :=
  replayPCA[run <> ".dyn", 0, stoptime, dtout]    

readRun[] :=
  readRun[runName]

(*** and finally, running simulations ***)

(*setupSim::usage = "setupSim[] resets all counters"*)

setupTables[] :=
  Module[{},
    Clear[rates,cumRates,event,cfgDep,globDep,totalRate,outcome,calculated];
    setupEventTable[];
    compactEvtTable[];
    setupOutcome[];
    determineCfgDependence[];
    determineGlobDependence[];
   ]


setupSim[] :=
  Module[{},
    setupTables[];
    eventRecord = {};
    eventsOccurred = 0;
    initRandomState = $RandomState;
    nextOut = starttime + dtout;
    edgeStep = dt/noOfEdges;
   ]

pickState[cdist_] :=
  omega[[pickFromCumDist[cdist]]]




recordEventOld[e_edge,newPair_] := 
   (If[recordEvents,Write[theRecord,{time[],e,state[e],newPair}]];
    Decrement[pairCount[state[e]]];
    Decrement[pairCount[state[Reverse[e]]]];
    Map[(Decrement[pairCount[#]])&,config[e],{2}];
    Map[(Decrement[pairCount[Reverse[#]]])&,config[e],{2}];
    Map[Decrement[singleCount[#]]&,state[e]];
    setState[e,newPair];
    Map[Increment[singleCount[#]]&,newPair];
    Increment[pairCount[newPair]];
    Increment[pairCount[Reverse[newPair]]];
    Map[(Increment[pairCount[#]])&,config[e],{2}];
    Map[(Increment[pairCount[Reverse[#]]])&,config[e],{2}];
    ++eventsOccurred;
    newPair
   )

   
decPair[{x_,y_}] :=
  Module[{},
    --pairCount[{y,x}];
    --pairCount[{x,y}]
   ]
incPair[{x_,y_}] :=
  Module[{},
    ++pairCount[{y,x}];
    ++pairCount[{x,y}]
   ]

recordEvent[e_edge,newPair_] :=
  Module[{local},
    If[logEvents,Write[theRecord,{time[],e,state[e],newPair}]];
    local = Append[ Flatten[neighbourhood[e]],e ];
    (decPair[state[#]])& /@ local;
    Decrement[singleCount[#]]& /@ state[e];
    setState[e,newPair];
    Increment[singleCount[#]]& /@ newPair;
    (incPair[state[#]])& /@ local;
    ++eventsOccurred;
    newPair
   ]

recordEventFailSafe[e_edge,newPair_] :=
  Module[{},
    If[logEvents,Write[theRecord,{time[],e,state[e],newPair}]];
    setState[e,newPair];
    countLattice[];
    ++eventsOccurred;
    newPair
   ]

checkCounts[] :=
(*
   Module[{c1,c2},
      c1 = pairCount /@ omegaSquared;
      countLattice[];
      c2 = pairCount /@ omegaSquared;
      If[c1=!=c2, Print["HUH?\n",c1,"\n",c2]];
     ]*) Null
      

reportList := {Hold[time[]]}

reportState[] := 
  Map[Print[#," = ",N[ReleaseHold[#]]]&,reportList]

myNumberForm[x_] :=
" "<>ToString[IntegerPart[N[10000 x]]]







sigmaFileName[tim_] :=
  ToFileName[{Directory[],runName},"sigma"<>myNumberForm[tim]]
sigmaFileName[tim_,run_] :=
  ToFileName[{Directory[],run},"sigma"<>myNumberForm[tim]]

  (*in pascal use Write(x:12:5) *)
Off[NumberForm::"sigz"]

(* to reduce the size of the file --- significant 'cos there may 
be a great number of them --- only the ordinal number of the state is
saved *)

Clear[ordOmega]
ordOmega[x_:singletQ] := ordOmega[x] = Position[omega, x, 1][[1, 1]]
ordSigma[x_] := ordOmega[sigma[x]]


saveSigma[fileName_] := 
  Put[{N[time[],12], Map[ordSigma,siteList]}, fileName]

readSigma[tim_,run_] :=
   Block[{fn},
         fn = sigmaFileName[tim,run];
         If[FileType[fn] === File,
            theExpr = ReadList[fn];
           {time[], ordSigmaList} = theExpr[[1]];
           sigmaList = omega[[#]]& /@ ordSigmaList;
           MapThread[setState, {siteList, sigmaList}];
           True,
           (*or else*) False]]

readSigma[tim_] := readSigma[tim,runName]

appendState[fn_String] :=
  PutAppend[N[{time[], singleCount/@omega, pairCount/@omegaSquared}],fn]

readState[] :=
   Block[{sc,pc},
       {time[], sc, pc} = nextExpression;
        MapThread[(singleCount[#1] = #2) &, {omega, sc}];
        MapThread[(pairCount[#1] = #2) &, {omegaSquared, pc}];
        ]

  
simPCAdynOutFn[] := runName<>".dyn"
simPCAnetOutFn[] := runName<>".net"

report[] :=
   (If[logLattice, saveSigma[sigmaFileName[time[]]]];
    appendState[simPCAdynOutFn[]];
    reportState[];)

makeLogFile[] :=
  Block[{logfile},
     logfile = OpenWrite["simpca.log", 
                    FormatType->OutputForm, PageWidth->Infinity];
     Write[logfile,"simpca version 28 for 7-3-2004, Minus van Baalen"];
     Write[logfile,"mathematica output"];
     Close[logfile];
   ];

initSim[] :=
   (setupSim[];
    saveLattice[simPCAnetOutFn[]];
    seedLattice[basicSimulationRules];
    If[FileType[runName]=!=Directory,CreateDirectory[runName]];
    makeLogFile;
    time[]=starttime;
    resetHistory[];
    addToHistory[];
    If[FileType[simPCAdynOutFn[]]===File,
       DeleteFile[Evaluate[simPCAdynOutFn[]]]];
    report[];)

continueSim[newStoptime_:stoptime] :=
   (stoptime = newStoptime;
    stopSim = False;
    While[ (time[]<stoptime)&&Not[stopSim],
          (doEdgeEvent[];
           time[] += edgeStep;
           If[time[] >= nextOut,
              (addToHistory[];
               report[];
               nextOut += dtout)]) ];
    report[];)




prepLattice[] :=
  Block[{},
    Print["Creating New Lattice."];
    createLattice[latticeType, rows, cols];
    Print["Seeding Lattice."];
    seedLattice[simulationRules];
    time[] := starttime;
   ]

simPCA[] := 
  Block[{},
    If[FileType[runName]=!=Directory,CreateDirectory[runName]];
    initSim[];
    continueSim[];
    ]






replayPCA[outputFile_String, startReplay_, stopReplay_, dtReplay_] := 
  Block[{strm, nextExpression, sc, pc}, 
    If[FileType[outputFile] === File, 
      resetHistory[];
      strm = OpenRead[outputFile];
      time[] = -Infinity;
      nextExpression = Read[strm, Expression];
      While[(nextExpression =!= EndOfFile) && (time[] < stopReplay), 
        (*
        {time[], sc, pc} = nextExpression;
        MapThread[(singleCount[#1] = #2) &, {omega, sc}];
        MapThread[(pairCount[#1] = #2) &, {omegaSquared, 2 pc}];*)
        readState[];
        addToHistory[];
        nextExpression = Read[strm, Expression];
        ];
      Close[strm];]]
   


xrunThrough[toDo_,runName_, startReplay_, stopReplay_, dtReplay_] :=
  Block[{tim, fn}, 
    If[FileType[runName]===Directory,
       resetHistory[];
       For[tim = startReplay, tim <= stopReplay, tim += dtReplay,
         fn = sigmaFileName[tim];
         If[FileType[fn] == File,
            theExpr = ReadList[fn];
           {time[], ordSigmaList} = theExpr[[1]];
           sigmaList = omega[[#]]& /@ ordSigmaList;
           MapThread[setState, {siteList, sigmaList}];
           addToHistory[];
           toDo]];
      ]] 

runThrough[toDo_,runName_, startReplay_, stopReplay_, dtReplay_] :=
  Block[{tim, fn}, 
    If[FileType[runName]===Directory,
       resetHistory[];
       For[tim = startReplay, tim <= stopReplay, tim += dtReplay,
         If[readSigma[tim],
           addToHistory[];
           toDo]];
      ]] 
      
SetAttributes[runThrough,HoldFirst]

echo[x_] := (Print[x]; x)

fixform[n_, w_, d_] := 
  StringJoin@(ToString /@ cff[RealDigits[n, 10, w], w - d])

cff[{dl_, pt_}, apt_] := 
  Which[pt === apt, dl, 
        pt < apt, cff[{Prepend[Drop[dl, -1], 0], pt + 1}, apt],
        True, xxx]

$animationSize={700,700}

saveFrame[n_] :=
  Display[echo["frame "<>fixform[n,5,0]],
     Show[plotLattice[],    
       plotTime,
       AspectRatio->Automatic,
       PlotRange:>$latticeRange,
       DisplayFunction->Identity],"PICT",ImageSize->$animationSize]



makeAnimPicts[runName_, startReplay_, stopReplay_, dtReplay_] :=
  Block[{n=1},
    runThrough[saveFrame[n++],runName,startReplay,stopReplay,dtReplay]]

makeAnimation[runName_, startReplay_, stopReplay_, dtReplay_] :=
  runThrough[showLattice[],runName,startReplay,stopReplay,dtReplay]

rerun[runName_, startReplay_, stopReplay_, dtReplay_] :=
  runThrough[reportState[], runName, startReplay, stopReplay, dtReplay]

rerun[runName_] :=
  rerun[runName,starttime,stoptime,dtout]

(* saving results... *)

saveLattice[filename_] :=
  Module[{},
   If[ FileType[filename]===File, DeleteFile[filename]];
   Save[filename,
       {noOfSites,
        siteList,
        pos,
        loc,
        noOfEdges,
        edgeList,
        $nbh,
        $enbh,
        filename}]
       ]

saveState[filename_String] :=
  Module[{},
   If[ FileType[filename]===File, DeleteFile[filename]];
   Save[filename,
        {initRandomState,
         $RandomState,
         method,      (* calculation method *)
         starttime,
         stoptime,
         dt,          (* timestep *)
         dtout,
         time,        (* current time *)
         edgeStep,
         nextOut,
         omega,        (* set of states *)
         eventList,    (* events defining the model *)
         basicSimulationRules,   (* to obtain numerical values *)
         evt,          (* events applicable to pairs *)
         rates,
         outcomes,
         cfgDep,       (* says whether a rate depends on its config *)
         globDep,      (* says whether a rate depends on globals *)
         sigma,        (* current state of the lattice *)
         ct,           (* quantities that have been tracked *)
         eventsOccurred,
         lastTimeAdded,
         runName
         }]
        ]


(* deriving statistics *)

time[current] := time[]
theTime := time[history]

p[x_?pairQ,pca[t_]] := ct[x][t]/noOfEdges
p[x_?singletQ,pca[t_]] := ct[x][t]/noOfSites
p[{x_?singletQ},pca[t_]] := ct[x][t]/noOfSites

p[x_,current] := p[x,pca[time[]]]


p[x_, timeSeries] := Transpose[p[x, history]][[2]]

q[x_,y_,t_] := p[{x,y},t]/p[y,t]

p[x_||y_,t___] := p[x,t] + p[y,t]
p[{x_,y_||z_},t___] := p[{x,y},t] + p[{x,z},t]
p[{x_||y_,z_},t___] := p[{x,z},t] + p[{y,z},t]


pcaSeries[time] := theTime
pcaSeries[p[x_]] := p[x,pca[history]]
pcaSeries[q[x_,y_]] := q[x,y,pca[history]]

timeSeries[x_]:=Transpose[{theTime,pcaSeries[x]}]


PCAversion







