(*--------------------------------------------------------------------------*)
(*----------------------------- INTERACT -----------------------------------*)
(*--------------------------------------------------------------------------*)

(* This package contains functions, examples of configurations and transi-  *)
(* tion rates for the simulation of spin system on two-dimensional grids    *)
(* The grids are supposed to be toric (periodic boundary conditions), and   *)
(* the flip rates at one site depend on the state (0 or 1) of the site as   *) 
(* well as the number of its neighbors in state 1.                          *) 


BeginPackage["UVW`Interact`"]
(* Version 2.0, 08/15/97 *)

(*----------------- Examples of basic configurations -----------------------*)
(* All configurations are two-dimensional arrays of 0's and 1's. The element*)
(* config[[x,y]] is interpreted as the state of site (x,y).                 *)  
(*--------------------------------------------------------------------------*)

Checkerboard::usage="
Checkerboard returns a 40-by-40 configuration of 0's and 1's arranged in
10 by 10 squares."

Diagonals::usage="
Diagonals returns a 40-by-40 configuration of 0's and 1's arranged in 
diagonal stripes."

RConfig::usage="
RConfig[p,width,height] returns a width-by-height array of independent 
random 0's and 1's, 1 being chosen with probability p."

(*------------------------ Examples of rates -------------------------------*)
(* All lists of rates are returned as a 2-by-5 list of reals. in such a     *)
(* list, rate[[i,j]] represents the rate at which the configuration will    *)
(* change at a site in state i (0 or 1) having j (from 0 to 4) neighbors in *)
(* state 1.                                                                 *)
(*--------------------------------------------------------------------------*)

Uniform::usage="
Uniform[lambda,mu] are the rates corresponding to the case where the        
configuration changes from 0 to 1 at rate lambda and from 1 to 0 at rate mu,
independently from the number of neighbors in state 1."

Ising::usage="
Ising[Alpha,Beta] returns the rates corresponding to the symmetric
Stochastic Ising Model, admitting a Gibbs measure as a reversible state.
Alpha is the potential of a site alone, Beta is the potential of a pair 
of neighboring sites."

Contact::usage="
Contact[lambda] returns the rates of the contact process. The transition
rates from 1 to 0 (curing) is constant. The rate of transition from
0 to 1 is proportional to the number of neighbors at 1 (infecting)."

Voter::usage="
Voter returns the rates of the Voter Model, where the transition rate from 0 
to 1 is proportional to the number of neighbors in state 1 and vice versa."

(*------------------- Treatment of a configuration -------------------------*)

Cyclic::usage="
Cyclic[n,bound] returns bound if n is 0, 1 if n is bound+1 , n in any other 
case. (Periodic boundary conditions are assumed for all configurations.)"

RepartConfig::usage="
RepartConfig[config] returns a 2-by-5 list of integers. Its element
[[i,j]] is the number of sites in state i (0 or 1) having j (from 0 to 4) 
neighbors in state 1."

Evolution::usage="
Evolution[initialconfig,rates,niter] simulates the evolution of a 
configuration according to the spin system corresponding to rates. niter
iteration are performed. One iteration consists of picking up a site at
random and decide to flip its state or not."

(*-------------------- Graphical representation ----------------------------*)

DrawConfig::usage="
DrawConfig[config,opts] plots config as a rectangular array of back and 
white squares."

(*--------------------------------------------------------------------------*)
(*--------------------------------------------------------------------------*)

Begin["`Private`"]

(*----------------- Examples of basic configurations -----------------------*)

Checkerboard :=
    Block[{i,j},
          Return[Table[If[ Xor[ Mod[Quotient[i,5],2]==0,
                                Mod[Quotient[j,5],2]==0 ]
                         ,1,0],
                       {i,0,39},{j,0,39}]
                ];
         ];

Diagonals :=
    Block[{i,j},
          Return[Table[If[ Mod[i+j,10]<5  
                         ,1,0],
                       {i,0,39},{j,0,39}]
                ];
         ];

RConfig[p_Real,height_Integer,width_Integer] :=
    Table[If[Random[]<p,1,0],{height},{width}];


(*------------------------ Examples of rates -------------------------------*)

Uniform[lambda_,mu_]:=
    Transpose[Table[{N[lambda],N[mu]},{5}]];

Ising[alpha_,beta_]:=
    Block[{i},
          Return[Transpose[
                   Table[{N[Exp[alpha+(2*i-4)*beta]],
                          N[Exp[-alpha+(4-2*i)*beta]]}
                         ,{i,0,4}]
                          ]
                ];
         ];

Contact[lambda_]:=
    Block[{i},
          Return[Transpose[Table[{N[lambda*i] , 1.} , {i,0,4}]]];
         ];

Voter :=
    Block[{i},
          Return[Transpose[Table[{N[i] , N[4-i]} , {i,0,4}]]];
         ];

(*------------------- Treatment of a configuration -------------------------*)

Cyclic[n_Integer,bound_Integer]:= 
    Switch[n,0,bound,bound+1,1,_,n];

RepartConfig[config_List]:=
    Block[{w,h,x,y,state,nei,dist},
                                       (* Dimensions of the grid------------*)
          w    = Length[Transpose[config]];
          h    = Length[config];

          dist = Table[0,{2},{5}];

          Do[(                         (* For all sites---------------------*)
                                       (* State of site (x,y)---------------*)
              state = config[[x,y]];
                                       (* Count its neighbors at 1----------*)
              nei = 1;
              nei = nei + config[[ Cyclic[(x-1),h],y]];
              nei = nei + config[[ Cyclic[(x+1),h],y]];
              nei = nei + config[[x,Cyclic[(y-1),w] ]];
              nei = nei + config[[x,Cyclic[(y+1),w] ]];
              state = state+1;
                                       (* Increment dist[[state]][[nei]]----*)
              dist[[state,nei]] = dist[[state,nei]]+1;
                                            
             ),{x,1,h},{y,1,w}];

          Return[dist];
         ];



Evolution[initialconfig_List,rates_List,niter_Integer]:=

    Block[{w,h,config,proba,x,y,nei},
                   (* INITIALIZATIONS *)
                                       (* Dimensions of the grid------------*)
          w      = Length[Transpose[initialconfig]];
          h      = Length[initialconfig];
                                       (* The curent configuration----------*)
          config = initialconfig;
                                       (* Probabilities to flip a site in a *)
                                       (* given state with n neighbors------*)
          proba = rates/Max[rates];

                      (* MAIN LOOP *)

          Do[
             (                         (* Choose a site at random-----------*)
              x = Random[Integer,{1,h}];
              y = Random[Integer,{1,w}];
                                       (* Count its neighbors at 1----------*)
              nei = 1;
              nei = nei + config[[ Cyclic[x-1,h],y]];
              nei = nei + config[[ Cyclic[x+1,h],y]];
              nei = nei + config[[x,Cyclic[y-1,w] ]];
              nei = nei + config[[x,Cyclic[y+1,w] ]];
                                       (* Decide to flip it or not----------*)
              If[Random[]<proba[[ 1+config[[x,y]],nei ]],
                                       (* Flip it---------------------------*)
                  config[[x,y]] = 1-config[[x,y]] 
                ];
             ),{niter}];

             Return[config];
         ];

(*-------------------- Graphical representation ----------------------------*)
DrawConfig[config_List,opts___]:=
    Show[Graphics[Raster[config]], opts,
                      Frame->True,
                      FrameTicks->None,
                      AspectRatio->
                        N[Length[config]/Length[Transpose[config]]]
        ]

End[ ]        (* End of Private Context-------------------------------------*)


EndPackage[ ]
(*--------------------------------------------------------------------------*)
