
      (*********************************)
      (** Improved Pair Approximation **)
      (**         methods             **)
      (**       Minus van Baalen      **)
      (**      (c)  1996 - 2004       **)
      (*********************************)

(*BeginPackage["ipa`"]*)

(* things to be defined by the user *)

IPAversion = "IPA.v16 June 23, 2004"

omega::usage = "omega should be the set of states"

phi::usage = "phi should be the connectedness parameter"

ihp::usage = "ihp is the complement of phi"

q::usage = "q[h,{i,j}] is the probability of finding an h 
next to an {i,j} pair"

p::usage = "p[x] gives the frequency of x (x may be singles, pairs or whatever)"

parameterList::usage = "parameterList, when defined as, {param -> default,...}
may lead to more tidy expressions"

eventList::usage = "eventList should be a list of {event[p1,p2] -> expr}"

ipaState::usage = "ipaState produces a vector with pair states"
ipaGrad::usage = "ipaGrad produces the differential equations for ipaState"

event::usage = "event[p1,p2] represents the mean rate of p1 -> p2
events"

latticeType::usage 
"specifies type of lattice"


theta
ateht 
ihp

(* definitions *)

defineModel::usage = "NOT IMPLEMENTED
defineModel[states,events,parameters] can be used to
define a model with other symbols than the default ones"


(* functions giving the results of the analysis *)

pairEvent::usage = "pairEvent[p1,p2] gives the mean rate of p1 -> p2 events"

transRate::usage = "transRate[p1,p2] gives the mean rate of p1 -> p2
events, including events in p1's neighbourhood"

transMat::usage = "transMat gives the entire transition matrix"

pairs::usage = "pairs gives the set of pairs"

represPairs::usage = "represPairs gives the set of representative pairs
(i.e. excludes mirror images)"

transMatRepr::usage = "transMatRepr gives the reduced transition matrix
(i.e. for the representative pairs)"

residentStates::usage = "residentStates[m] gives the set of states excluding
the mutant"

residentPairs::usage = "residentPairs[m] gives the set of states excluding
pairs with mutant m"

mutantPairs::usage = "mutantPairs[m] gives the set of states including
mutant m"

residentPairsRepr::usage = "residentPairsRepr[m] gives the set of states 
excluding pairs with mutant m and without mirror images"

mutantPairsRepr::usage = "mutantPairsRepr[m] gives the set of states 
excluding pairs with mutant m and without mirror images"

speciesDyn::usage = "speciesDyn[i] produces a set of differential equations for 
the pairs including species i"

stateVec::usage = "stateVec gives the full state vector of pairs"

stateVecRepr::usage = "stateVecRepr gives the state vector of representative pairs"

mutantVec::usage = "mutantVec[m] gives the state vector of mutant pairs"

mutantVecRepr::usage = "mutantVec[m] gives the state vector of 
representative mutant pairs"

gradFull::usage = "gradFull gives the differential equations for the full system"

gradRepr::usage = "gradRepr gives the differential equations for the full representative
system"

meanField::usage = "meanField produces the mean field equations"

speciesGrad::usage = "speciesGrad[i] gives the differential equations for species i"

mutantGrad::usage = "mutantGrad[i] gives the differential equations for species i
when rare"

invasionMatrix::usage = "invasionMat[i] gives the invasion matrix of species i"


(* rules to calculate p and q from pairs *)

pairApproximation::usage = 
"standard approximation of q[h,{i,j}]"

triangleApproximation::usage = 
"approximation of q[h,{i,j}] assuming a proportion theta of triangles"

thetaApproximation::usage = 
"approximation of q[h,{i,j}] assuming a proportion theta of triangles"

latticeTypes = {randomLattice, triangleLattice, thetaLattice}
   
qApproximation::usage = 
"qApproximation[latticeType] gives an expression approximating q[h,{i,j}]"

qExpand::usage = 
"expresses q[i,j] in terms of pair frequencies"

pExpand ::usage =  
"expresses p[{i}] in terms of pair frequencies"   

varExpand::usage =  
"expresses var[i,j] in terms of pair frequencies"   

pRules::usage = 
"gives a list of all pExpand"

qRules::usage = 
"gives a list of all qExpand"

qTripleRules::usage =
"gives a list for all q[h,{i,j}] "

pqRules::usage =
"combines all rules for p and q"

domesticRules::usage =
"required to pass on to Madonna"

neighbourApproximation::usage =
"gives the rules to approximate neighbours[state,{x,y}] (used in PCA simulations) 
in terms of local densities"

(* Auxiliary functions *)

checkUnity::usage = "checkUnity is an utility"

pExpansion::usage = "pExpansion gives rules to calculate p[i] from the pairs"

stateQ::usage = "stateQ[i] is true if i is a member of omega"

pairQ::usage = "pairQ[p] is true if p is a member of pairs"

qList::usage = "qList[qRule] gives the list of all expressions for q using qRule"

mirrorQ::usage = "mirrorQ[pair] yields true if the mirror image of pair is in
the set of representative pairs"

swap::usage = "swap[{i,j}] yields {j,i}"


(*** some additional functions ***)

TableOverSet[tosexpr_,{tosidx_,tosset_List}] := 
   Map[Function[tosidx,tosexpr],tosset]
   
SetAttributes[TableOverSet,HoldFirst]
   
SumOverSet[sosexpr_,{sosidx_,sosset_List}] := 
   Apply[Plus,TableOverSet[sosexpr,{sosidx,sosset}]]
   
myNDSolveAux[gradVec_, stateVecFun_, stateVec_, 
    iniCond_, {time_, start_, stop_}, ops___] := 
  NDSolve[Join[Thread[D[stateVecFun, time] == gradVec], 
      Thread[(stateVecFun /. time -> start) == iniCond]], 
    stateVec, {time, start, stop}, ops]
makefunc[svec_, t_] := Map[(# -> #[t]) &, svec]
myNDSolve[gradVec_, stateVec_, iniCond_, {time_, start_, stop_}, ops___] := 
  myNDSolveAux[gradVec /. makefunc[stateVec, time], 
      stateVec /. makefunc[stateVec, time],
      stateVec, iniCond, {time, start, stop}, ops][[1]]
myNDSolve::usage =
  "myNDSolve[{g1,g2,..},{s1,s2,...},{a1,a2,...},{t,t0,t1}] invokes NDSolve to solve the set of differential equations xi'[t] == gi subject to initial condition xi[t0] = ai"




(*Begin["`prv`"] no privacy!*)



defineModel[states_,events_,parameters_] :=
   Block[{}, omega = states;
             eventList = events;
             parameterList = parameters]


(*** some useful rules for p and q ***)

q[x_||y_,z_] := q[x,z] + q[y,z]
q[x_,y_||z_] := (p[{x,y}]+p[{x,z}])/(p[y]+p[z])


p[x_?stateQ,t___] := p[{x},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]







(* To construct the full transition matrix *)

transMat :=
  TableOverSet[
     TableOverSet[
        transElmt[kl,ij],
        {kl,pairs}],
     {ij,pairs}]

transElmt[p1_?pairQ,p1_?pairQ] :=
   - SumOverSet[transRate[p1,p2],{p2,pairs}]
   
transElmt[p1_?pairQ,p2_?pairQ] :=
   transRate[p1,p2]

transRate[{sti_,stj_},{sti_,stj_}] := 0

transRate[{sti_,stj_},{stk_,stj_}] :=
   Collect[
   checkPhi [ Expand [
   
   phi pairEvent[{sti,stj},{stk,stj},{sti,stj}] +
   
   ihp checkUnity[ 
                SumOverSet[ 
		           SumOverSet[
		              q[sth,{sti,stj}] pairEvent[{sth,sti},{sth2,stk},{sth,sti,stj} ],
		              {sth2, omega} ],
		           {sth,omega}]]]]
		           
		           ,q]

transRate[{stj_,sti_},{stj_,stk_}] :=
   transRate[{sti,stj},{stk,stj}]

transRate[{sti_,stj_},{stk_,stl_}] :=
   pairEvent[{sti,stj},{stk,stl},{sti,stj}]

(*
pairEvent[p1_?pairQ,p2_?pairQ,___] :=
   event[p1,p2]/phi /. eventList /.
                         swapEvent /.
                         eventList /.
                         event[_,_] -> 0 /.focus -> p1/.$focus->p1
*)                        

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



swapEvent :=
   event[{sti_,stj_},{stk_,stl_}] -> event[{stj,sti},{stl,stk}]












(* To construct the representative transition matrix *)

transMatReprOriginal :=
   TableOverSet[
      TableOverSet[
         transElmtRepr[stkl,stij],
         {stkl,represPairs}],
      {stij,represPairs}]

transElmtRepr[{stk_,stl_},{stk_,stl_}] :=
   - multiplicity[{stk,stl}] SumOverSet[transRateRepr[{stk,stl},stij],{stij, represPairs}] 
transElmtRepr[{stk_,stl_},{sti_,stj_}] := 
   multiplicity[{sti,stj}] transRateRepr[{stk,stl},{sti,stj}]

transRateRepr[{x_ ,y1_},{y2_,x_ }] := transRate[{x ,y1},{x ,y2}]
transRateRepr[{x1_,y_ },{y_ ,x2_}] := transRate[{x1,y },{x2,y }]
transRateRepr[{x1_,y1_},{x2_,y2_}] := transRate[{x1,y1},{x2,y2}]      
      

multiplicity[{x_,x_}] := 2
multiplicity[{x_,y_}] := 1

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











(* The following implements a faster way to construct the transition matrix *)

transElmtReprBasic[{stk_,stl_},{stk_,stl_}] := 0
transElmtReprBasic[{stk_,stl_},{sti_,stj_}] := transRateRepr[{stk,stl},{sti,stj}]

transMatReprBasic :=
   TableOverSet[
      TableOverSet[
         transElmtReprBasic[kl,ij],
         {kl, represPairs}],
      {ij, represPairs}]

setDiagonal[mat_] :=
   mat - IdentityMatrix[Length[mat]] Apply[Plus,mat]

setMultiplicity[mat_] :=
   mat Map[multiplicity,represPairs]

transMatRepr := 
   setMultiplicity[setDiagonal[transMatReprBasic]]












(* Constructing state vectors, etc. *)

pairs := pairList[omega]

represPairs := removeMirrors[pairs]

stateVec := Map[p,pairs]
stateVecRepr := Map[p,represPairs]

pairQ[p_] := MemberQ[pairs,p]

gradFull := transMat.stateVec

gradRepr := transMatRepr.stateVecRepr





(* To derive mean field equations *)
(* added March 8th 2003 *)

mfState := p /@ omega

mfGrad := 
  (mfDiffEqn /@ omega) 

mfDiffEqn[i_?stateQ] :=
  (Plus @@ ((allRates[{i,#}])& /@ omega)) /. p[{x_,y_}] -> p[{x}]p[{y}] /. 
     nb[x_,{y_,z_}] -> p[{x}]/phi

allRates[p1_?pairQ] :=
  Plus @@ ((- p[p1] pairEvent[p1,#] + p[#] pairEvent[#,p1])& /@ pairs) 

conservationLaw :=
  p[First[omega]] -> 1 - Plus @@ (p /@ Rest[omega])








(* mutant dynamics *)

mutantPairs[m_] := 
   Map[ findRepr,
   TableOverSet[{m,sti},{sti,omega}]
   ]

mutantVec[sti_] := Map[p,mutantPairs[sti]]
mutantVecRepr[sti_] := Map[p,mutantPairs[sti]]

residentPairs[m_] := 
   Map[ findRepr,
   pairList[residentStates[m]]
   ]
   
residentPairsRepr[m_] := removeMirrors[residentPairs[m]]


residentStates[m_List] := Complement[omega,m]
residentStates[m_] := Complement[omega,{m}]



findRepr[p_?mirrorQ] := swap[p]
findRepr[p_] := p














(* To derive invasion dynamics *)

speciesMat[sti_] :=
   TableOverSet[
      TableOverSet[
         transElmtRepr[stkl,stij],
         {stkl,represPairs}],
      {stij,mutantPairs[sti]}]
      
speciesGrad[sti_] :=
   TableOverSet[
      SumOverSet[
         transElmtRepr[stkl,stij] p[stkl],
         {stkl,represPairs}],
      {stij,mutantPairs[sti]}]
      
rewriteRes[m_] :=
   Flatten[ 
      TableOverSet[ 
         TableOverSet[ p[{iloc,jloc}] q[m,{iloc,jloc}] ->
                          p[{m,iloc}] q[jloc,{iloc,m}],
                       {iloc,residentStates[m]}],
         {jloc,residentStates[m]}]]

speciesRare[m_] := 
  Join[
    (q[m,#] -> 0)& /@ residentPairs[m] ,
    (q[m,#] -> 0)& /@ residentStates[m]
    ]

invasionMatrix[m_] := 
    Dj[ Expand[speciesGrad[m]]//.rewriteRes[m]//.speciesRare[m], mutantVec[m] ]


Dj[f_List,x_List] :=
   Table[Table[D[f[[dindexi]],x[[dindexj]]],{dindexj,Length[x]}],{dindexi,Length[x]}];






(* Dealing with conditional probabilities *)

pairApproximation := 
   q[h_,{i_,j_}] -> q[h,i]

corr[i_,j_] := p[{i,j}]/(p[{i}] p[{j}])


thetaApproximation := 
   {q[j_, {i_, j_}] :> 
      1 - SumOverSet[q[h, {i, j}], {h, Complement[omega, {j}]}], 
    q[h_, {i_, j_}] -> q[h, i] (ateht + theta q[h, j]/p[{h}])}


thetaApproximationNew := 
   {q[j_, {i_, j_}] :> 
      1 - SumOverSet[q[h, {i, j}], {h, Complement[omega, {j}]}], 
    q[h_, {i_, j_}] -> q[h, i] ateht + 
                       theta p[{h,i}] p[{h,j}]/(p[{h}] p[{i}] p[{j}])}

qApproximation[wrong] := 
   {q[h_, {i_, j_}] -> q[h, i] ateht + 
                       theta p[{h,i}] p[{h,j}]/(p[{h}] p[{i}] p[{j}])}



qApproximation[heuristic] :=
  q[i_,{a_,b_}] ->
     q[i,a] (1+theta (corr[i,b]-1))/
            (1+theta (Apply[Plus,Map[(q[#,a] corr[#,b])&,omega]]-1))

qApproximation[simplistic] :=
  q[i_,{a_,b_}] ->
     ateht q[i,a] + theta q[i,b] q[a,i]/
                    Apply[Plus,Map[(q[#,b] q[b,#])&,omega]]

qApproximation[reallySimplistic] :=
  q[i_,{a_,b_}] ->
     ateht q[i,a] + theta q[i,b]
     
     
 
    
qApproximation[randomLattice] := pairApproximation
qApproximation[triangleLattice] := thetaApproximation
qApproximation[thetaLattice] := thetaApproximation
qApproximation[triangular] := qApproximation[triangleLattice]
qApproximation[_] := qApproximation[randomLattice]


qExpand := 
   q[h_,i_] -> p[{h,i}]/p[{h}]

pExpand := 
   p[{i_}] -> SumOverSet[p[{i,j}],{j,omega}]
   
varExpand :=
   var[i_,j_] -> p[{i,j}]/(p[{i}] p[{j}])

makeRuleList[theRule,cases] :=
   cases //. theRule

pRules := 
    Flatten[
      TableOverSet[p[{x}] -> SumOverSet[p[{x,y}],{y,omega}],
                   {x, omega}]
           ]

qRules :=
    Flatten[
      TableOverSet[
         TableOverSet[q[x,y] -> p[{x,y}]/p[{y}], {x,omega}],
                   {y, omega}]
           ]

varExpansion :=
    Flatten[
      TableOverSet[
         TableOverSet[var[x,y] -> p[{x,y}]/(p[{x}] p[{y}]), {y,omega}],
                   {x, omega}]
           ]


qTripleRules[qRule_] :=
  Flatten[
    TableOverSet[ 
       TableOverSet[ temp[h,ij] -> q[h,ij], {h,omega}], 
       {ij, pairs}]] /. qRule //. temp -> q


pqRules[latticeType_] :=
   Join[ qTripleRules[qApproximation[latticeType]],
         qRules,
         pRules ]


 
domesticRules := {ihp -> 1-phi, ateht -> 1-theta}

tempRules := {}

simulationRules :=
  Join @@ {tempRules, parameterValues, domesticRules, pqRules[triangular]}
basicSimulationRules :=
  Join @@ {tempRules, parameterValues, domesticRules}

(* odds & ends *)


vec[x_List] := Map[p,x]

pairList[set_] := 
   Flatten[
      TableOverSet[
         TableOverSet[{x,y},
                      {y, set}], 
         {x, set}],
      1]

mirrorQ[x_?pairQ] := Not[MemberQ[represPairs,x]]

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

p[x_?mirrorQ] := p[swap[x]]

(*removeMirrors[pairs_] := Union[Map[Sort,pairs]]*)
(*this works, but doesn't preserve the order of pairs*)

removeMirrors[pairs_] :=
   rmAux[pairs,{}]

rmAux[{},pairList_] := pairList

rmAux[x_List,pairList_] :=
   If[ MemberQ[pairList,swap[First[x]]], 
          rmAux[Rest[x],pairList],
          rmAux[Rest[x],Append[pairList,First[x]]]]

        


checkUnity[expr_] := expr //.
   {SumOverSet[q[y,x_],{y,omega}] -> 1, 
    SumOverSet[c_ q[y,x_],{y,omega}] -> c}

checkPhi[expr_] := expr //.
   {c_ phi + c_ ihp -> c,
    -c_ phi - c_ ihp -> -c}


$na[x_,{y_,x_}] := ihp q[x,{y,x}]/phi + 1

$na[x_,{y_,z_}] := ihp q[x,{y,z}]/phi

neighbourApproximation :=
  {nb[x_, {y_, z_}] :> $na[x,{y,z}],
   nb[{z_, y_}, x_] :> $na[x,{y,z}]}

same[x_, y_] := If[x === y, 1, 0]

nbCalc[x_,{y_,z_}] := nb[x, {y, z}] + same[x, z]
nbCalc[{z_,y_},x_] := nb[{z, y}, x] + same[x, z]

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

ipaGrad := gradRepr //. neighbourApproximation

ipaState := stateVecRepr

ipaInit := init /@ represPairs

p[{x_?stateQ, y_?stateQ}, ipa[t_]] := p[{x, y}][t] /. IPAsol
p[{x_?stateQ}, ipa[t_]] := Plus @@ (p[{x, #}, ipa[t]] & /@ omega)

p[{x_?stateQ}, mf[t_]] := p[{x}][t] /. MFsol
p[{x_?stateQ, y_:stateQ}, mf[t_]] := p[{x}][t] p[{y}][t] /. MFsol

ipaTest := ipaGradN /. Thread[ipaState -> ipaInitN]
 
makeMadonnaFile[fName_String] :=
  makeMadonnaFile[fName,
     ipaGrad, {t,ipaState}, 
              {starttime,ipaInit},dt,(stoptime-starttime)/dt,simulationRules]


(* invasion analysis *)
(* some of the stuff I need probably can be found in Mathematica *)
(* already but I could not be bothered to look it up *)
(* problem is to find normalized Eigenvector associated with an Eigenvalue of 0 *)

firstCol[mat_] := Delete[Transpose[mat][[1]], 1]
subMat[mat_] := Transpose[Delete[Transpose[Delete[mat, 1]], 1]]
ev[mat_] := Prepend[Inverse[subMat[mat]].(-firstCol[mat]), 1]
norm[v_List] := v/Plus @@ v
qev[mat_] := norm[ev[mat]]



(*End[]*)
(*EndPackage[ ]*)

IPAversion

