Dynamo 2x2 (v. 1.3)
Diagrams for Evolutionary Game Dynamics

When you have finished specifying parameters, click on the rightmost vertical bar and press Enter (or Shift-Enter) to run the program.

avoiderroronremove=1;
Unprotect[Xt];
Remove["Global`*"];
NotebookFind[SelectedNotebook[], "corecell", All, CellTags, AutoScroll -> False];
FrontEndExecute[FrontEndToken["RemoveFromEvaluationQueue"]];
NotebookFind[SelectedNotebook[], "initcell", All, CellTags, AutoScroll -> False];
FrontEndExecute[FrontEndToken["EvaluateCells"]];

NotebookFind[SelectedNotebook[], "corecell", All, CellTags, AutoScroll -> False];
FrontEndExecute[FrontEndToken["EvaluateCells"]];
Off[General::spell1];

User-defined Parameters

Choice of game

Functional form

Select the Functional form of the game. If the game is Linear, use the buttons to specify the structure of  F[1][x_,y_]  and F[2][x_,y_] .Then  fill the values for the payof matrices in the following section. 
To choose a nonlinear game, define F[1][x_,y_]  and F[2][x_,y_] as you like. In doing so, the components of x should be entered as x[[1]], x[[2]],  and the components of y should be entered similarly.  The default value of F[1][x_,y_] is A.y and F[2][x_,y_] is B’.x, the payoff vector field for the normal form game U. As an example, imagine we wanted to input the normal form structure as a nonlinear game. We could first define the payoff Matrix A={{1,0},{0,2}} and then write:
F[1][x_,y_]:={    A[[1,1]] y[[1]] + A[[1,2]] y[[2]]      ,    A[[2,1]] y[[1]] + A[[2,2]] y[[2]]   };
The most general way of defining payoffs is with a function F[x_] that assumes that x is a list of 2 vectors with 2 components each, so the components would be accessed as  x[[1,1]],x[[1,2]],x[[2,1]],x[[2,2]]; for an advanced example of this more general way of defining payoffs use the congestion game button.

matrixButtons[rowsize_]:=Module[{namesl,textsl,names ,texts,buttons,row,butrow,rown,rowt,name,text,
but,button},
namesl = {"Normal Form","General Linear","Congestion"};
textsl = {"F[1][x_,y_]:=A.y;
F[2][x_,y_]:=Transpose[B].x;
gamekind=NormalGame;",
"F[1][x_,y_]:=M[[1,1]].x+M[[1,2]].y;
F[2][x_,y_]:=M[[2,1]].x +M[[2,2]].y;
gamekind=LinearGame;","F[x_]:=Table[-Table[Sum[congestionCost[Φ[p][j][[k]]][u[Φ[j][k],x]],{k,Dimensions[Φ[p][j]][[1]]}],{j,dim[[p]]}],{p,pop}];
gamekind=CongestionGame;"};
names={};
texts={};
rown={};
rowt={};
row=1;
While[row<=Dimensions[namesl][[1]],
rown=Join[rown,{namesl[[row]]}];
rowt=Join[rowt,{textsl[[row]]}];
If[Or[IntegerQ[row/rowsize],row==Dimensions[namesl][[1]]],
names=Join[names,{rown}];
texts=Join[texts,{rowt}];
rown={};
rowt={};
];
row++;
];
buttons={};
row=1;
While[row<=Dimensions[names][[1]],
butrow={};
but=1;
While[ but<=Dimensions[names[[row]]][[1]],
button:=
Button[Column[{name},Alignment->{Center}],
(*NotebookLocate["unique_payofffunction"];*)
NotebookFind[SelectedNotebook[],"unique_payofffunction",Next,CellTags,AutoScroll->False];
SelectionMove[SelectedNotebook[],All,CellContents,AutoScroll->False];
NotebookWrite[SelectedNotebook[],RowBox[{text}],AutoScroll->False];
(*,ImageSize->{120,30}*)];
button=Style[button,ControlsRendering->"Generic"] /.text->texts[[row]][[but]]/.name->names[[row]][[but]];
butrow=Join[butrow,{button}];
(*butrow=Join[butrow,{Row[{" "}],button}];*)
but++;
];
buttons=Join[buttons,{butrow}];
row++;
];
TableForm [ buttons,TableAlignments->{Left,Top},TableSpacing->{1,1}]
(*Column[buttons]*)
];
matrixButtons[3]

F[1][x_,y_]:=A.y;
F[2][x_,y_]:=Transpose[B].x;
gamekind=NormalGame;

Payoff-related definitions

The multi population payoff

If[Or[
gamekind==NormalGame,
gamekind==LinearGame
],
Fmp[x_]:={F[1][x[[1]],x[[2]]],F[2][x[[1]],x[[2]]]};,
Fmp[x_]:=F[x];
];

The population's average payoff

Fbar[x_,F_]:=Table[Sum[x[[p,j]] F[x][[p,j]],{j,dim[[p]]}] ,{p,pop}];

Excess payoff vector

Fhat[x_,F_]:=Table[F[x][[p]] - Fbar[x,F][[p]],{p,pop}];

Vector of positive parts of excess payoffs

Fhatplus[x_,F_]:= Table[Table[Max[0,Fhat[x,F][[p,j]]],{j,dim[[p]]}], {p,pop}];

Projected payoff vector

PhiF[x_,F_]:= Table[projectionMatrix[{dim[[p]]}].F[x][[p]],{p,pop}];

Specification of payoff parameters

Specification of Normal Form Game

Choose a normal form game here by editing the payoff matrix below.  To use one of the built-in payoff matrices, click the corresponding button and the expression below where "U = ." appears will be updated.

Some standard normal form games

matrixButtons[rowsize_]:=Module[{namesl,matsl,names ,mats,buttons,row,butrow,rown,rowm,name,mat,but,button},
namesl = {"12 Coordination","Pure Coordination","Matching Pennies","Prisoner's Dilemma","Stag Hunt","Chain Store"};
matsl =  {
{{{{1,1}},{{0,0}}},{{{0,0}},{{2,2}}}},{{{{1,1}},{{0,0}}},{{{0,0}},{{1,1}}}},{{{{1,-1}},{{-1,1}}},{{{-1,1}},{{1,-1}}}},{{{{2,2}},{{0,3}}},{{{3,0}},{{1,1}}}},{{{{5,5}},{{0,4}}},{{{4,0}},{{2,2}}}},{{{{0,4}},{{0,4}}},{{{2,2}},{{-4,-4}}}}
};
names={};
mats={};
rown={};
rowm={};
row=1;
While[row<=Dimensions[namesl][[1]],
rown=Join[rown,{namesl[[row]]}];
rowm=Join[rowm,{matsl[[row]]}];
If[Or[IntegerQ[row/rowsize],row==Dimensions[namesl][[1]]],
names=Join[names,{rown}];
mats=Join[mats,{rowm}];
rown={};
rowm={};
];
row++;
];
buttons={};
row=1;
While[row<=Dimensions[names][[1]],
butrow={};
but=1;
While[ but<=Dimensions[names[[row]]][[1]],
button:=
Button[Column[{MatrixForm[mat],name},Alignment->{Center}],
(*NotebookLocate["unique_payoffmatrix"];*)
NotebookFind[SelectedNotebook[],"unique_payoffmatrix",Next,CellTags,AutoScroll->False];
SelectionMove[SelectedNotebook[],All,CellContents,AutoScroll->False];
NotebookWrite[SelectedNotebook[],RowBox[{"U=",ToBoxes[MatrixForm[mat]],";"}],AutoScroll->False];
(*,ImageSize->{150,80}*)];
button=Style[button,ControlsRendering->"Generic"] /.mat->mats[[row]][[but]]/.name->names[[row]][[but]];
butrow=Join[butrow,{button}];
(*butrow=Join[butrow,{Row[{" "}],button}];*)
but++;
];
buttons=Join[buttons,{butrow}];
row++;
];
Print[TableForm[buttons,TableAlignments->{Left,Top},TableSpacing->{1,1}],FromCharacterCode[10]]
(*Column[buttons]*)
];
matrixButtons[7]

U=(({{(({{1, 1}})), (({{0, 0}}))}, {(({{0, 0}})), (({{2, 2}}))}}));

A=(({{U[[1,1,1,1]], U[[1,2,1,1]]}, {U[[2,1,1,1]], U[[2,2,1,1]]}})) ;
B=(({{U[[1,1,1,2]], U[[1,2,1,2]]}, {U[[2,1,1,2]], U[[2,2,1,2]]}})) ;

Specification of General Linear Form Game

Choose a General Linear game here by specifying the values of the matrices of payoffs. The rows change depending on which population's payoff is being specified. The columns change depending on which population's state is being interacted with the parameters to produce the payoff. By default, the matrices are set to produce the 12 Coordination Normal Form Game.

M=(({{(({{0, 0}, {0, 0}})), (({{1, 0}, {0, 2}}))}, {(({{1, 0}, {0, 2}})), (({{0, 0}, {0, 0}}))}}));

Specification of nonlinear game

Definition of Congestion Game

Cost functions for facilities as functions of utilization levels.

congestionCost[1][u_] :=4+20u ; 
congestionCost[2][u_] :=2+30u^2;
congestionCost[3][u_] :=2+30u^2; 
congestionCost[4][u_] :=4+20u; 
congestionCost[5][u_] :=1+20u;
congestionCost[6][u_] :=0 ; 
congestionCost[7][u_] :=0 ;

Which strategies use which facilities.

Clear[Φ];
Φ[1][1]={1,2}; 
Φ[1][2]={3,4};
Φ[2][1]={2,3,5};
Φ[2][2]={1,4};
ρ[ϕ_]:=Table[DeleteCases[Table[
If[MemberQ[Φ[p][j],ϕ],j]
,{j,dim[[p]]}],Null],{p,pop}];
u[ϕ_,x_]:=Sum[Sum[
x[[p,ρ[ϕ][[p,j]]]]
,{j,Dimensions[ρ[ϕ][[p]]][[1]]}],{p,pop}] ;

Names of strategies

Give the strategies names to be printed on the phase diagram.  Vertices are ordered counterclockwise.

strategy[1][1]="U";
strategy[1][2]="D";
strategy[2][1]="L";
strategy[2][2]="R";
font="Palatino";

Computation of Nash equilibria

If the game you consider has a nontrivial component of Nash equilibria, the program will find a random subset of this component.  If you’d prefer not to compute the Nash equilibria, set findnashequilibria to 0. Similarly, you may skip finding the regular ESSs of the game by setting findESSs to 0. The option maxdenominator allows to control how the equilibria are displayed in the output. The output is presented in either decimal or rational form. If a fraction has a denominator that is bigger than maxdenominator, the number is presented in decimal form. This helps avoid having fractions with big denominators that are visually not appealing, but can be problematic if analytic results are desired. Setting the maximum denominator to “Infinity”, without quotes, makes the program display the exact rational numbers.

findnashequilibria=1 ;
findESSs=1;
maxdenominator=20;

Choice of dynamic

To specify a dynamic, press the button corresponding to the dynamic you desire.  Definitions of the dynamics can be found in the closed group at the end of this section.
Details:  If you use Logit[eta] or ILogit[eta], be sure to enter a value for the noise parameter eta in the closed group below.  The best response dynamic is actually the Logit[.001] dynamic, but with a different method of determining stability of rest points.  When using the SampleBR[k] dynamic, set an integer value of k > 1 in the closed group below. SampleBR is most noticeably different from a best response dynamic only for values of k below 5.

Some standard dynamics

matrixButtons[rowsize_]:=Module[{namesl,matsl,names ,buttons,row,butrow,rown,name,but,button},
namesl = {Replicator,BR,Defer[Logit[eta]],MSReplicator,TemperedBR[tQ],BNN,Defer[ILogit[eta]],Defer[SampleBR[k]],Smith};
names={};
rown={};
row=1;
While[row<=Dimensions[namesl][[1]],
rown=Join[rown,{namesl[[row]]}];
If[Or[IntegerQ[row/rowsize],row==Dimensions[namesl][[1]]],
names=Join[names,{rown}];
rown={};
];
row++;
];
buttons={};
row=1;
While[row<=Dimensions[names][[1]],
butrow={};
but=1;
While[ but<=Dimensions[names[[row]]][[1]],
button:=
Button[name,
(*NotebookLocate["unique_meandynamic"];*)
NotebookFind[SelectedNotebook[],"unique_meandynamic",Next,CellTags,AutoScroll->False];
SelectionMove[SelectedNotebook[],All,CellContents,AutoScroll->False];
NotebookWrite[SelectedNotebook[],RowBox[{"dyn:=",ToBoxes[name],";"}],AutoScroll->False];
(*,ImageSize->{100,30}*)];
button=Style[button,ControlsRendering->"Generic"]/.name->names[[row]][[but]];
(*butrow=Join[butrow,{Row[{" "}],button}];*)
If[ToString[names[[row]][[but]]]== " ",butrow=Join[butrow,{""}],butrow=Join[butrow,{button}]];(*butrow=Join[butrow,{button}];*)
but++;
];
buttons=Join[buttons,{butrow}];
row++;
];
TableForm[buttons,TableAlignments->{Left,Top},TableSpacing->{1,1}]
(*Column[buttons]*)
];
matrixButtons[3]

More dynamics

If you use SelMut[MM], enter the mutation matrix MM in its closed group.  ExcessPayoff , PairwiseComparison, Combined, and Other must be specified in greater detail below.

matrixButtons[rowsize_]:=Module[{namesl,matsl,names ,buttons,row,butrow,rown,name,but,button},
namesl = {Defer[SelMut[MM]],ExcessPayoff,PairwiseComparison,Project,Combined,Other};
names={};
rown={};
row=1;
While[row<=Dimensions[namesl][[1]],
rown=Join[rown,{namesl[[row]]}];
If[Or[IntegerQ[row/rowsize],row==Dimensions[namesl][[1]]],
names=Join[names,{rown}];
rown={};
];
row++;
];
buttons={};
row=1;
While[row<=Dimensions[names][[1]],
butrow={};
but=1;
While[ but<=Dimensions[names[[row]]][[1]],
button:=
Button[name,
(*NotebookLocate["unique_meandynamic"];*)
NotebookFind[SelectedNotebook[],"unique_meandynamic",Next,CellTags,AutoScroll->False];
SelectionMove[SelectedNotebook[],All,CellContents,AutoScroll->False];
NotebookWrite[SelectedNotebook[],RowBox[{"dyn:=",ToBoxes[name],";"}],AutoScroll->False];
(*,ImageSize->{100,30}*)];
button=Style[button,ControlsRendering->"Generic"]/.name->names[[row]][[but]];
(*butrow=Join[butrow,{Row[{" "}],button}];*)
butrow=Join[butrow,{button}];
but++;
];
buttons=Join[buttons,{butrow}];
row++;
];
TableForm[buttons,TableAlignments->{Left,Top},TableSpacing->{1,1}]
(*Column[buttons]*)
];
matrixButtons[6]

dyn := Replicator;

Noise level for logit dynamics

eta = .1;

Sample size for the sample best response dynamic

k = 2;

Transformation function for Tempered dynamics

tQ[x_]:=x;

Matrix of mutation rates for the selection-mutation dynamic

Entry ij of the array MM is the probability of a mutation from i to j ≠ i; entry ii is the probability of no mutation from i.   Note that the selection-mutation dynamic is only well defined for nonnegative payoffs, and that it is not invariant to the addition of a constant to all matrix entries.

MM={(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))};

dyn=dyn;

If the dynamic you consider has a nontrivial component of rest points, the program will compute and plot a random subset of this components.  If you set findrestpoints equal to 0, the computation and plotting is skipped.

findrestpoints=1;

Definition of combination of two dynamics

To define a combination of two dynamics, specify the two dynamics to be combined, the weights  on the dynamics, and the characterization of the rest points of the combined dynamic.  "Nash" is the default setting for the set or rest points.  This is correct, for example, with combinations of the replicator dynamic and an excess payoff dynamic like the BNN dynamic.  If the rest points aren't identical to the Nash equilibria for the dynamics you have specified, you can replace "Nash" with "Automatic" if the combined dynamic is continuous.  (If the combined dynamic isn't continuous, you're on your own!)  The last line of code formally defines the combined dynamic according to your specifications.

dyn[1]=Replicator;
dyn[2]=BNN;
weight[1]=.9;
weight[2]=.1;
RPCharacterization[Combined]="Nash";
SCharacterization[Combined]="Unknown";
Combined[x_,F_]:=weight[1]dyn[1][x,F]+weight[2] dyn[2][x,F];

Definitions of dynamics

Replicator Dynamic

Replicator[x_,F_]:=Table[Table[x[[p,j]]  Fhat[x,F][[p,j]],{j,dim[[p]]}],{p,pop}];

RPCharacterization[Replicator]="Automatic" ;
SCharacterization[Replicator]="Smooth";
AOCharacterization[Replicator]=False;

Maynard Smith Replicator Dynamic

MSReplicator[x_,F_]:=Table[(1/(Fbar[x,F][[p]]))Replicator[x,F][[p]],{p,pop}];

MSReplicator::negativepayoff = "Warning, Normal Form game with negative Payoff Matrix detected, the payoff matrix has been rescaled.";
If[dyn===MSReplicator,
Which[
And[gamekind==NormalGame,Min[U]<0],
U=U-Min[U]+1;
A=(({{U[[1,1,1,1]], U[[1,2,1,1]]}, {U[[2,1,1,1]], U[[2,2,1,1]]}})) ;
B=(({{U[[1,1,1,2]], U[[1,2,1,2]]}, {U[[2,1,1,2]], U[[2,2,1,2]]}})) ;
Message[MSReplicator::negativepayoff];
,
And[gamekind==LinearGame,Min[M]<0],
M=M-Min[M]+1;
Message[MSReplicator::negativepayoff];
];
If[Or[
gamekind==NormalGame,
gamekind==LinearGame
],
Module[{n,points},
n=10;
points=Tuples[Table[discSimplex[n,dim[[p]]],{p,pop}]];
Do[While[
Apply[And,Table[Fbar[points[[j]],Fmp][[p]]<=0,{p,pop}]],
U=U+1;
A=(({{U[[1,1,1,1]], U[[1,2,1,1]]}, {U[[2,1,1,1]], U[[2,2,1,1]]}})) ;
B=(({{U[[1,1,1,2]], U[[1,2,1,2]]}, {U[[2,1,1,2]], U[[2,2,1,2]]}})) ;
M=M+1;
],{j,Dimensions[points][[1]]}];
];
];
];

RPCharacterization[MSReplicator]="Automatic" ;
SCharacterization[MSReplicator]="HiperESS";
AOCharacterization[MSReplicator]=False;

I-Logit Dynamic

ILogit[eta_][x_,F_]:=Table[Table[((x[[p,j]] (E^((F[x][[p,j]])/eta)))/((∑_(k=1)^(dim[[p]]))x[[p,k]](E^((F[x][[p,k]])/eta)))) -x[[p,j]],{j,dim[[p]]}],{p,pop}] ;

RPCharacterization[ILogit[_]]="Automatic" ;
SCharacterization[ILogit[_]]="HiperESS";
AOCharacterization[ILogit[_]]=False;

Best Response dynamic

BR[x_,F_]:=Logit[.001][x,F];

RPCharacterization[BR]="Nash" ;
SCharacterization[BR]="ESS";
AOCharacterization[BR]=False;

Tempered Logit

TemperedLogit[eta_][Q_][x_, F_] := Module[{payoff,max},
payoff=F[x];
payoff=Table[Table[Q[payoff[[p,j]]],{j,dim[[p]]}],{p,pop}];
max=Table[Max[payoff[[p]]],{p,pop}];
Table[Table[
Exp[payoff[[p,j]]/eta]/Sum[Exp[payoff[[p,k]]/eta], {k,dim[[p]]}]Sum[x[[p,k]]*(max[[p]]-payoff[[p,k]]),{k,dim[[p]]}]
-x[[p,j]]*(max[[p]] -payoff[[p,j]])
,{j,dim[[p]]}],{p,pop}]
];

RPCharacterization[TemperedLogit[_][_]] = "Automatic" ;
SCharacterization[TemperedLogit[_][_]]="Unknown";
AOCharacterization[TemperedLogit[_][_]]=False;

Tempered Best Response

The tempered Best Response dynamic is an approximation through tempered logit of the tempered best response dynamic of Zusai (2011).

TemperedBR[Q_][x_, F_] :=TemperedLogit[.001][Q][x, F] ;

RPCharacterization[TemperedBR[_]] = "Automatic" ;
SCharacterization[TemperedBR[_]]="Unknown";
AOCharacterization[TemperedBR[_]]=False;

SampleBR Dynamic

The sample best response dynamic is introduced in Oyama, Sandholm, and Tercieux (2011).

sampleBRSimplex=Tuples[Table[discSimplex[k,dim[[p]]],{p,pop}]];
SampleBR[k_][x_,F_]:=Module[{Znk,Bk},
Znk=sampleBRSimplex;
Bk=Sum[Table[bestResponse[Znk[[l]],F][[p]]*Apply[Multinomial,k *Znk[[l,p]]]*Product[Limit[(x[[p,j]]+littleo)^(k *Znk[[l,p,j]]),littleo->0,Direction->-1],{j,dim[[p]]}]
,{p,pop}],{l,Dimensions[Znk][[1]]}];
Bk-x
];

RPCharacterization[SampleBR[_]]="Automatic" ;
SCharacterization[SampleBR[_]]="Smooth";
AOCharacterization[SampleBR[_]]=False;

Logit Dynamic

Logit[eta_][x_,F_]:=Table[Table[((E^((F[x][[p,j]])/eta))/((∑_(k=1)^(dim[[p]]))(E^((F[x][[p,k]])/eta))))-x[[p,j]],{j,dim[[p]]}],{p,pop}];

RPCharacterization[Logit[_]]="Automatic" ;
SCharacterization[Logit[_]]="Smooth";
AOCharacterization[Logit[_]]=False;

BNN dynamic

BNN[x_,F_]:=Table[Table[Fhatplus[x,F][[p,j]]-x[[p,j]] (∑_(k=1)^(dim[[p]]))Fhatplus[x,F][[p,k]],{j,dim[[p]]}],{p,pop}];

RPCharacterization[BNN]="Nash" ;
SCharacterization[BNN]="ESS";
AOCharacterization[BNN]=False;

Smith Dynamic

rho[x_,F_]:=Module[{Fx},
Fx=F[x];
Table[Table[Table[Max[Fx[[p,j]]-Fx[[p,k]],0],{j,dim[[p]]}],{k,dim[[p]]}],{p,pop}]
];

Smith[x_,F_]:=Table[Transpose[rho[x,F][[p]]]. x[[p]]-DiagonalMatrix[x[[p]]].rho[x,F][[p]].ConstantArray[1,dim[[p]]],{p,pop}];

RPCharacterization[Smith]="Nash" ;
SCharacterization[Smith]="ESS";
AOCharacterization[Smith]=False;

Excess payoff dynamic

Here one specifies the vector field sigmatilde that defines an excess payoff dynamic.  This definition should be stated in terms of the excess payoff vector Fhat[x,F], which is defined as Fhat[x,F] = F[x] - Fbar[x,F].  (Mathematica subtracts the scalar from the vector correctly.)  When computing rest points, the program assumes that sigmatilde is acute, so that the rest points of the dynamic are the Nash equilibria of the underlying game.

sigmatilde[x_,F_]:= Table[Table[(Max[0,Fhat[x,F][[p,j]]])^2,{j,dim[[p]]}],{p,pop}];

ExcessPayoff[x_,F_]:=Table[Table[sigmatilde[x,F][[p,j]]-x[[p,j]] (∑_(k=1)^(dim[[p]]))sigmatilde[x,F][[p,k]],{j,dim[[p]]}],{p,pop}];

RPCharacterization[ExcessPayoff]="Nash"  ;
SCharacterization[ExcessPayoff]="ESS";
AOCharacterization[ExcessPayoff]=False;

Pairwise Comparison Dynamic

Here one specifies the vector field rho that defines an excess payoff dynamic.  When computing rest points, the program assumes that rho is separable and sign-preserving, so that the rest points of the dynamic are the Nash equilibria of the underlying game. The leading example of the vector field is the one for the Smith dynamic defined above and it is the default.

rho1[x_,F_]:=rho[x,F];

PairwiseComparison[x_,F_]:=Table[Transpose[rho1[x,F][[p]]]. x[[p]]-DiagonalMatrix[x[[p]]].rho1[x,F][[p]].ConstantArray[1,dim[[p]]],{p,pop}];

RPCharacterization[PairwiseComparison]="Nash" ;
SCharacterization[PairwiseComparison]="ESS";
AOCharacterization[PairwiseComparison]=False;

Projection Dynamic

Project[x:pattern,F_]:=Module[{j,k,m,p,nc,Fx,boundaries,reduced,reducedx,npos,flag,project,result},
result=Table[ConstantArray[0,dim[[p]]],{p,pop}];
Fx=F[x];
For[p=1,p<=pop,p++,
If[Total[support[x[[p]]]]==dim[[p]],
result[[p]]=projectionMatrix[{dim[[p]]}].Fx[[p]];
,
nc=2^dim[[p]];
For[j=2,j<=nc,j++,
boundaries=numberSystem[dim[[p]],2,j];
If[Max[support[x[[p]],3]-boundaries]==1,Continue[],
npos=1;
reduced=ConstantArray[0,dim[[p]]];
reducedx=ConstantArray[0,dim[[p]]];
Do[
If[boundaries[[k]]==1,
reduced[[npos]]=Fx[[p,k]];
npos++;
];
,{k,dim[[p]]}];
npos--;
project=IdentityMatrix[npos]-1/npos*ConstantArray[1,{npos,npos}];
reduced=project.Take[reduced,npos];
flag=False;
Do[
If[And[boundaries[[k]]==0,Fx[[p,Position[boundaries,1][[1,1]]]]-Fx[[p,k]]-reduced[[1]]<0],(*result[[k]]=-x[[k]]*)
flag=True;
Break[];
];
,{k,dim[[p]]}];
If[flag,Continue[],
m=1;
Do[
If[boundaries[[k]]==1,
result[[p,k]]=reduced[[m]];
m++;
];
,{k,dim[[p]]}];
Break[];
];
];
];
];
];
result
];

RPCharacterization[Project]="Nash" ;
SCharacterization[Project]="ESS";
AOCharacterization[Project]=True;

Selection-mutation

SelMut[M_][x_,F_]:= Table[Transpose[M[[p]]] .DiagonalMatrix[x[[p]]].F[x][[p]]-(x[[p]].F[x][[p]])x[[p]],{p,pop}];

RPCharacterization[SelMut[___]]="Automatic" ;
SCharacterization[SelMut[___]]="Unknown";
AOCharacterization[SelMut[___]]=False;

Other

Define any other dynamic you want to use. When defining be sure you follow same style with previously defined dynamics. Otherwise, you may need to do major changes in the program.

Other[x_,F_]:=Table[Table[x[[p,j]]  Fhat[x,F][[p,j]],{j,dim[[p]]}],{p,pop}];

RPCharacterization[Other]="Automatic" ;
SCharacterization[Other]="NonSmooth";
AOCharacterization[Other]=False;

Choice of contour function

Choose a contour function here.  To use one of the built-in contour functions, press the button corresponding to the contour function you desire.  If you choose StableGameLyapunov[dyn], the program will use the Lyapunov function that is appropriate for the dynamic you selected in the previous section.  contourFunction is used in almost all figures including the background contour plots on phase diagrams (see the next section).

matrixButtons[]:=Module[{names ,buttons,col,colbut,name,but,button},
names ={{Speed,L1Speed,RSSpeed},{NormalFormPotential,CongestionPotential,NormalFormLogitPotential,CongestionLogitPotential},{Defer[StableGameLyapunov[dyn]]}};
buttons={};
col=1;
While[col<=Dimensions[names][[1]],
colbut={};
but=1;
While[ but<=Dimensions[names[[col]]][[1]],
button:=
Button[name,
(*NotebookLocate["unique_contourfunction"];*)
NotebookFind[SelectedNotebook[],"unique_contourfunction",Next,CellTags,AutoScroll->False];
SelectionMove[SelectedNotebook[],All,CellContents,AutoScroll->False];
NotebookWrite[SelectedNotebook[],RowBox[{"contourFunction=",ToBoxes[name],";"}],AutoScroll->False];
(*,ImageSize->{250,30}*)];
button=Style[button,ControlsRendering->"Generic"]/.name->names[[col]][[but]];
colbut=Join[colbut,{button}];
but++;
];
buttons=Join[buttons,{Column[colbut]}];
col++;
];
TableForm[{buttons},TableAlignments->{Center,Top}]
(*Row[buttons,Alignment->{Top}]*)
];
matrixButtons[]

contourFunction=Speed;

Definitions of contour functions

Speed functions

Speed is the standard (Euclidean) speed on the simplex.  L1Speed is self-explanatory.

Speed[x_]:= If[ dyn=== Project,stateNorm[PhiF[x,Fmp]],stateNorm[dyn[x,Fmp]]];

L1Speed[x_]:=stateNorm[dyn[x,Fmp],1];

RSSpeed[x_]:=stateNorm[Fhat[x,Fmp]];

Potential functions

perturbedLogitPotential[x_,eta_,m_:ConstantArray[1,pop]]:=Total[Table[m[[p]]eta Sum[x[[p,j]]/m [[p]]Log[Max[x[[p,j]]/m[[p]],10^(-7)]],{j,dim[[p]]}],{p,pop}]];

NormalFormPotential[x_]:=Total[Fbar[x,Fmp]]/Evaluate[If[pop==1,2,pop]];

NormalFormLogitPotential[x_,m_:ConstantArray[1,pop]]:=NormalFormPotential[x]-perturbedLogitPotential[x,eta];

CongestionPotential[x_]:=-Sum[Integrate[congestionCost[k][z],{z,0,u[k,x]}],{k,Max[Table[Table[Φ[p][j],{j,dim[[p]]}],{p,pop}]]}];

CongestionLogitPotential[x_]:=CongestionPotential[x]-perturbedLogitPotential[x,eta] ;

Lyapunov functions

StableGameLyapunov[Replicator][x_]:=- Product[Product[((x[[p,j]])^(NEa[[1,p,j]])) ,{j,dim[[p]]}],{p,pop}]  ;

StableGameLyapunov[Logit[eta]][x_]:= Sum[eta Log[(∑_(j=1)^(dim[[p]]))(E^((Fmp[x][[p,j]])/eta))]-x[[p]].Fmp[x][[p]],{p,pop}]-perturbedLogitPotential[x,eta];

StableGameLyapunov[BNN][x_]:=Sum[(∑_(j=1)^(dim[[p]]))((Fhatplus[x,Fmp][[p,j]])^2) ,{p,pop}] ;

StableGameLyapunov[Smith][x_]:=Sum[(∑_(k=1)^(dim[[p]]))x[[p,k]](∑_(j=1)^(dim[[p]]))(((Max[Fmp[x][[p,j]]-Fmp[x][[p,k]],0]))^2) ,{p,pop}] ;

StableGameLyapunov[Project][x_]:=Sum[Norm[x[[p]]-NEa[[1,p]]]^2,{p,pop}];

Other contour functions

Define your own contour function here.

Other[x_]:=Sum[ConstantArray[0,dim[[p]]],{p,pop}] ;

Specification of graphical output

Specifications for phase diagram

These parameters are used for contour function output.  If you do not include contour functions in your output, these parameters are not used.
drawphasediagram = 1 creates a phase diagram on the simplex.  The solution trajectories to be included are specified below.
ppoints determines how many points will be used when drawing the solution trajectories. 
pdcontourplot = 1 draws the phase diagram on top of a contour plot.
pdcolor = 1 generates color output; pdcolor = 0 generates black and white.  The scales used are drawn below.
conbrightness specifies how light the black and white coloring should be. 1 is the default, 1.5 will give a balanced overall lighter image. The option is ineffective when color is selected.
plotprecision specifies the number of points sampled in each dimension when drawing the contour plots.  Typically, this option is the main determinant of the program's running time.  Settings between 50 and 200 are accurate enough for most purposes.  For fast drafts a setting of 20 is adequate.
numberofcontours specifies the number of contour levels shown in the contour plots.
compressgraphic = 1 uses an unofficial patch to generate a contour plot that is both of higher quality and smaller size (see http://library.wolfram.com/infocenter/MathSource/7029/ ).  In particular, Mathematica’s routine for converting output into eps and pdf format works poorly unless this patch is used.  Using this patch can substantially increase the program’s running time (say, by a factor of 5 or more), so it is typically used only to create figures for publication. If compressgraphic = 0, the original Mathematica code, which since version 6 has created lower-quality output, is used.

drawphasediagram=1;
ppoints=50;
pdcontourplot=1;
pdcolor=1; 
conbrightness=1;
plotprecision=50;
numberofcontours=50;
compressgraphic =0;

Solution trajectories

Instructions for "automated" specification of solution trajectories

This section specifies the solution trajectories to be included in the phase diagram using built-in specifications.  These are of three types:
(i) initial conditions drawn from a grid, with the grid including points on the boundary ;  
(ii) initial conditions drawn from a grid, with the grid excluding points on the boundary ;
(iii) randomly chosen initial conditions.  
To use one of these options, press the appropriate button below.
endtime determines the time interval [0, endtime] of each solution.  
For the grid options, gridpoints1d specifies the number of initial conditions on each boundary of the grid.  
For the random option, randompoints specifies the number of random initial conditions.
numberofarrows specifies how many arrows to draw on each solution.  
arrowheadlength is self-explanatory; it can be set to zero if arrowheads are not desired.  
If numberofarrows = 1 is specified, the arrow appears at the end of each solution trajectory.  If numberofarrows = n, arrows are drawn at times {endtime/n,  2 endtime/n, ... , endtime}.  For more on how arrows are drawn, see “Adjusting positions of arrow heads and ‘tails’” below.

Instructions for "manual" specification of solution trajectories

The user can also specify which solution trajectories to graph directly.  To do this, trajectoryspecs must be defined as a table of the following form: {{{initial condition 1}, T1, {color1, thickness1, dashing1},{Arrow positions 1}, {Head lengths 1}}, {{initial condition 2}, T2, {color2, thickness2, dashing2}, {Arrow positions 2}, {Head lengths 2}}, ... }
Each initial condition is entered as a state (e.g.,{{1/2, 1/2},{2/3,1/3}}) which gives a list that contains for each population, one vector in the simplex.  Ti determines time interval [0, Ti] of the ith solution.  Arrow positions i and Head lengths i should be vectors of the same length, where this length is the number of arrowheads to be drawn on the solution trajectory.  Arrow positions are stated as points in the time interval [0, Ti].  For no arrows, set Arrow positions i and Head lengths i to {0}.  

{colori, thicknessi, dashingi} specifies the appearance of solutions trajectories.  You can specify any collection of graphics primitives you like (even ones beyond the three mentioned here), and you can place them in any order.  To use the default appearance, write {}. 
colori can be Hue[?], GrayLevel[?] or any other color command.
thicknessi can be Thickness[?] or AbsoluteThickness[?].
dashingi is specified as Dashing[{???}].

endtime=1;
 
trajgridpoints=15; 
trajrandompoints=6 ;

trajnumberofarrows=1 ;
trajarrowheadsize=.025 ;

matrixButtons[rowsize_]:=Module[{namesl,textsl,discard,names ,texts,buttons,row,butrow,rown,rowt,name,text,but,button},
namesl ={ "Grid initial conditions (boundary included)","Grid initial conditions (boundary excluded)","Random initial conditions"};
textsl = {
"trajectoryspecs =Module[{ds,arrows,arrowheads}, 
ds=Tuples[Table[discSimplex[trajgridpoints-1,dim[[p]]],{p,pop}]];
arrows=Table[((trajnumberofarrows - j + 1)*endtime)/trajnumberofarrows,{j,trajnumberofarrows}];
arrowheads=ConstantArray[trajarrowheadsize,trajnumberofarrows];
Table[{ds[[j]],endtime,{Black},arrows,arrowheads},{j,Dimensions[ds][[1]]}]
];"
,
"trajectoryspecs =Module[{ds,dsp,arrows,arrowheads}, 
ds=Tuples[Table[
dsp=discSimplex[trajgridpoints,dim[[p]]]+1/(dim[[p]]*trajgridpoints);
dsp=Table[Append[Drop[dsp[[j]],-1],dsp[[j,dim[[p]]]]-1/trajgridpoints],{j,Dimensions[dsp][[1]]}];
DeleteCases[Table[If[Min[dsp[[j]]]<0,Infinity,dsp[[j]]],{j,Dimensions[dsp][[1]]}],Infinity],{p,pop}]];
arrows=Table[((trajnumberofarrows - j + 1)*endtime)/trajnumberofarrows,{j,trajnumberofarrows}];
arrowheads=ConstantArray[trajarrowheadsize,trajnumberofarrows];
Table[{ds[[j]],endtime,{Black},arrows,arrowheads},{j,Dimensions[ds][[1]]}]
];"
,
"trajectoryspecs=Module[{ds,arrows,arrowheads}, 
ds=Tuples[Table[sampleAround[ConstantArray[1/dim[[p]],dim[[p]]],1,trajrandompoints],{p,pop}]];
arrows=Table[((trajnumberofarrows - j + 1)*endtime)/trajnumberofarrows,{j,trajnumberofarrows}];
arrowheads=ConstantArray[trajarrowheadsize,trajnumberofarrows];
Table[{ds[[j]],endtime,{Black},arrows,arrowheads},{j,Dimensions[ds][[1]]}]
];"
};
names={};
texts={};
rown={};
rowt={};
row=1;
While[row<=Dimensions[namesl][[1]],
rown=Join[rown,{namesl[[row]]}];
rowt=Join[rowt,{textsl[[row]]}];
If[Or[IntegerQ[row/rowsize],row==Dimensions[namesl][[1]]],
names=Join[names,{rown}];
texts=Join[texts,{rowt}];
rown={};
rowt={};
];
row++;
];
buttons={};
row=1;
While[row<=Dimensions[names][[1]],
butrow={};
but=1;
While[ but<=Dimensions[names[[row]]][[1]],
button:=
Button[Column[{name},Alignment->{Center}],
NotebookFind[SelectedNotebook[],"unique_trajectoryspecs",Next,CellTags,AutoScroll->False];
SelectionMove[SelectedNotebook[],All,CellContents,AutoScroll->False];
NotebookWrite[SelectedNotebook[],text,AutoScroll->False];
(*,ImageSize->{400,30}*)];
button=Style[button,ControlsRendering->"Generic"] /.text->texts[[row]][[but]]/.name->names[[row]][[but]];
butrow=Join[butrow,{button}];
(*butrow=Join[butrow,{Row[{" "}],button}];*)
but++;
];
buttons=Join[buttons,{butrow}];
row++;
];
TableForm [ buttons,TableAlignments->{Left,Top},TableSpacing->{1,1}]
(*Column[buttons]*)
];
matrixButtons[3]

trajectoryspecs =Module[{ds,dsp,arrows,arrowheads}, 
ds=Tuples[Table[
dsp=discSimplex[trajgridpoints,dim[[p]]]+1/(dim[[p]]*trajgridpoints);
dsp=Table[Append[Drop[dsp[[j]],-1],dsp[[j,dim[[p]]]]-1/trajgridpoints],{j,Dimensions[dsp][[1]]}];
DeleteCases[Table[If[Min[dsp[[j]]]<0,Infinity,dsp[[j]]],{j,Dimensions[dsp][[1]]}],Infinity],{p,pop}]];
arrows=Table[((trajnumberofarrows - j + 1)*endtime)/trajnumberofarrows,{j,trajnumberofarrows}];
arrowheads=ConstantArray[trajarrowheadsize,trajnumberofarrows];
Table[{ds[[j]],endtime,{Black},arrows,arrowheads},{j,Dimensions[ds][[1]]}]
];

Tuning rest point search

restPointSearchStartingValues contains the initial conditions used in the FindRoot command when the program searches for rest points.  If you suspect that the program has missed some rest points, add initial conditions in the regions where the missing rest points should be.  Each initial condition is written as a state.
restTolerance is a parameter in a filtering procedure used to eliminate false rest points.  If the speed of the dynamic at a prospective rest point is above (restTolerance) * (the maximum speed of the dynamic), the prospective rest point is thrown out.

restPointSearchStartingValues=Join[Tuples[Table[discSimplex[10,dim[[p]],True,.001],{p,pop}]],{
{{.3,.7},{.6,.4}},
{{.4,.6},{.7,.3}}
}];

restTolerance=.01 ;

Tuning stability tests

The first four parameters determine how rest points are checked for stability.  The parameter stabn specifies how many trajectories will be checked.  These initial conditions are random points at distance stabeps from the rest point.  Solutions from these points over the time interval [0, stabT] are computed.  If any of these solutions winds up more than stabdelta away from the rest point, the rest point is considered unstable; otherwise it is considered stable. Notice that stabeps should be large relative to stabdelta and at least four times as big for the algorithm to work properly.

Extra notes: if stabdelta is set too small, stable rest points wont be detected unless stabT is large giving time for the trajectories to converge. However, for large values of stabT the algorithm will tend to interpret saddles as sources, although this is a distinction of secondary importance in general, so that in principle a safe way to run the algorithm is with relatively large values of stabeps compared with stabdelta and a generous time span stabT. There is a way to differentiate sources from saddles and cycles. Notice that there will be an inevitable computational burden since this distinction will require large values for stabn so that even for a large stabT after enough sampling around the rest point at least one of the sampled surrounding starting conditions will be close enough to the direction that converges to the rest point and stay in it's proximity in spite of the large time span. To have the algorithm determine whether a rest point is a sink, a source, a cycle or a saddle the numericStabilityTest[] routine is available. Evaluate numericStabilityTest::notes to learn about its arguments, and ?numericStabilityTest for its general definition and output values.

stabn=10; 
stabeps=.01 ; 
stabdelta=.0025;
stabT=10 ;

Adjusting the position of arrows

To draw an arrow, the program must know where to put the tail of the arrow and where to put the head.  Even if one does not want the base-line of the arrow to appear in the picture, the position of its starting point must be specified so that the arrow points in some direction.  
The parameter relarrowstart below determines the positions of arrows' base-lines  when arrows are drawn on solution trajectories.  In particular, if an arrow is supposed to appear at a solution trajectory's position at time t, the program draws the end of the line at this point, placing the starting point (and therefore changing the direction of the arrow) at the solution trajectory's position at time t(1 –relarrowstart).
A value of .005 for relarrowstart works well in most cases:  the arrow line is not visible, and the arrow points in the right direction.  If after a diagram is produced the user notices that arrow base-lines are visible, a smaller value of relarrowstart should be specified.  If arrows point in the wrong directions, a larger value of relarrowstart could make the arrows look more aligned with the overall trajectory instead of its final portion.
The parameter arrowheadpos  specifies the location of the arrow-head within the base-line. The number is relative so it ranges from 0 to 1 as a proportion of the base-line.

relarrowstart=.005 ;
arrowheadpos=.8 ;

Adjusting Rest Points and NE dots

If the Rest Points and/or NE of the game have been calculated, this options control whether to show them and their appearence. The first two options, drawrestpoints and drawnashequilibra include the dots in the phase diagram if they are set to one. Notice that since for most dynamics these sets overlap, a hierachy exists. Dynamo priorizes the drawrestpoint option and only if it is disabled will the rest points which are also NE be drawn with the necolor. Since for imitative dynamics every NE is a rest point, the option drawnashequilibra is ineffective in that context unless the option drawrestpoints is set to 0. However,in principle it is possible that a point will emerge as a NE without being a rest point, this is why both option are enabled by default.
The option dotradius controls how big the dots are. The three options named srpcolor, urpcolor and necolor allow the user to modify the default grayscale color for stable rest points, unstable rest points and nash equilibria respectively. the defaults are black (0), white(1) and gray (.5).

drawrestpoints=1;
drawnashequilibria=1;
dotradius=.01;
srpcolor=0;
urpcolor=1;
necolor=.5;

Adding extra graphic elements

This section allows the user to add dots, lines, arrows, text, and color solution trajectories to a phase diagram/contour plot.

Enter the data for dots as follows:  customdots = {{center1, radius1, color1, open1}, {center2, radius2, color2, open2}, ...}. color1 should be a Mathematica-defined color function (Hue[.], GrayLevel[.], etc).  Set open1 = 0 for solid dot and = 1 for a dot with a black edge.  Centers should be entered as states like {{x1,x2},{y1,y2}}.
 
Enter the data for lines as follows:  customlines = {{startingpoint1, endpoint1, thickness1, dashing1, color1, arrowheadsize1}, {startingpoint2, endpoint2, thickness2, dashing2, color2, arrowheadsize2},...}.  If you want an arrow at the end of the line enter the value for arrowheadsize specifying the size of arrow head; if you do not want an arrow simply set it to 0. Thickness and Dashing are usual Mathematica commands. If you want to enter positions using simplex coordinates, you can directly enter the state (e.g.,{{0, 1},{1, 0}}).  Otherwise, 2d coordinates can be used with "planeToSimplex" command:  for example, the bottom-left vertex is at planeToSimplex[{0, 0}] ={{0,1},{1,0}}.

Enter data for text as follows:  customtext = {{"label1", position1, "fontfamily1","fontslant1", fontsize1}, {"label2", position2, "fontfamily2","fontslant2", fontsize2}, ...}. Remember that you are restricted to positions in the surroundings of the simplex. If you want to enter position using simplex coordinates, use the "simplexToPlane" command (e.g., simplexToPlane{,{{0, 1},{1, 0}}]).

Enter data for extra solution trajectories as follows:  customtrajectoryspecs = {{initial condition on simplex 1, endtime1 for NDSolve, beginningtime1 to be used in in the plot, {color1, thickness1}, {arrow position1}, {head lengths1}}, {initial condition on simplex 2, endtime2 for NDSolve, beginningtime2 to be used in in the plot, {color2, thickness2}, {arrow position2}, {head lengths2}},...}. 
The program computes a solution to the dynamic from time 0 through endtime?, and then plots the solution trajectory from the beginning time? through endtime.   color? should be a Mathematica-defined color function (Hue[.], GrayLevel[.], etc) (both color and thickness can be ignored, just write {} if you do not want to bother with these.  arrowpositon? must be a number strictly greater than beginingtime? and less than endtime (it must be strictly bigger than begining time).  Each trajectory can have multiple arrows, but even if there is only one arrow, be sure to put brackets around the arrow positions and head lengths.

Write customdots = {},  customlines = {}, customtext = {}, and customtrajectoryspecs = {} if you do not want to add anything.  (This is the default setting.)

Example:  This example generates one line, one dot, and one piece of text, one special solution trajectory (customtrajectoryspecs).
customdots={{{{1/2,1/2},{1/2,1/2}}, .01, GrayLevel[0], 1}};    [draws a black dot with white center at the center of the simplex]
customlines={{{{0,1},{1/2,1/2}},{{1,0},{1/2,1/2}}, Thickness[.004], Dashing[{}], GrayLevel[0], 0}};   [draws a horizontal line through at y1=1/2]
customtext={{"trial", simplexToPlane[{{1/2,1/2},{1/2,1/2}}], "Palatino", "Plain", 12}}; 
customtrajectoryspecs={{{{.7,.3},{.6,.4}}, 50, 45, {Hue[.83], Thickness[.005]}, {47}, {.03}}};

customdots={};
customlines={};
customtext={};
customtrajectoryspecs={};
customtrajectoryspecs={};

These are extra, less used parameters that are used for displaying accessory elements in the contour function.  If you do not include contour functions in your output, these parameters are not used.
colorbar = 1 to create a color scale (alternatively, a black and white scale) along with your contour plot.
showcontourformula = 1 displays the formula for the contour function in the output notebook.
drawbasinboundaries =1 displays the approximate separatrices of the basins of attraction. (This code is work in progress!)

colorbar=0;
showcontourformula=1;
drawbasinboundaries=0;

Specifications for other graphics

Vector field diagram

drawvectorfield = 1 draws a vector field diagram on the simplex.  The vector fields should take values in the tangent space of the simplex (i.e., values whose components sum to zero).  
vfcontourplot = 1 draws the vector field on top of a contour plot.

drawvectorfield=0;
vfcontourplot=1;

Instructions for specification of vector fields

Either one or two vector fields is drawn, according to whether drawvf[1], drawvf[2], or both are set to 1.
vectorfieldf[?] specifies the vector field to be drawn.  (Here and below, "?" equals 1 or 2.)  Common specifications for this are dyn (the dynamic specified in the "Choice of dynamic" section) and PhiF (the projection of the payoff vector field onto the tangent space of the simplex - see "Payoff-related definitions" in the "Choice of game"section).
vfscale[?] rescales the vectors, as is often needed to prevent the vector field arrows from overlapping or looking too small. A value of one should ensure that no arrow overlaps with the other.
vfarrowheadsize[?] is the relative size of the arrow heads; set it to zero if you do not want arrowheads.
vfcolor[?] determines the color of the vector field.  It should be set to Hue[t] or GrayLevel[t] using the charts below.

The base points for the vectors in the diagram can be chosen from a grid 
(i) with the grid including points on the boundary;
(ii) with the grid excluding points on the boundary;
Option (i) is the default setting.
To switch options press the appropriate palette button.
vfgridpoints specifies the number of initial conditions on each dimension of the grid.
arrowsfromrestpoints = 1 causes vector field arrows to be drawn from the rest points (in addition to all other arrows), conditional on rest points beeing drawn in the first place.

vectorfieldf[1]=PhiF ;
drawvf[1]=1;
vfscale[1]=1.5;
vfarrowheadsize[1]=.8;
vfcolor[1]= Hue[.64];

vectorfieldf[2]=dyn;
drawvf[2]=1;
vfscale[2]=1.5;
vfarrowheadsize[2]=.8;
vfcolor[2]= Hue[.84];

vfgridpoints= 25;
arrowsfromrestpoints=1;

matrixButtons[rowsize_]:=Module[{namesl,textsl,names ,texts,buttons,row,butrow,rown,rowt,name,text,but,button},
namesl ={ "Grid initial conditions (boundary included)","Grid initial conditions (boundary excluded)"};
textsl = {
"vectorfieldspecs=Tuples[Table[discSimplex[vfgridpoints,dim[[p]]],{p,pop}]];"
,
"vectorfieldspecs=Tuples[Table[Module[{ds},
ds=discSimplex[vfgridpoints,dim[[p]]];
ds=ds+1/(dim[[p]]*vfgridpoints);
ds=Table[Append[Drop[ds[[j]],-1],ds[[j,dim[[p]]]]-1/vfgridpoints],{j,Dimensions[ds][[1]]}];ds=DeleteCases[Table[If[Min[ds[[j]]]<0,Infinity,ds[[j]]],{j,Dimensions[ds][[1]]}],Infinity]
],{p,pop}]];"
};
names={};
texts={};
rown={};
rowt={};
row=1;
While[row<=Dimensions[namesl][[1]],
rown=Join[rown,{namesl[[row]]}];
rowt=Join[rowt,{textsl[[row]]}];
If[Or[IntegerQ[row/rowsize],row==Dimensions[namesl][[1]]],
names=Join[names,{rown}];
texts=Join[texts,{rowt}];
rown={};
rowt={};
];
row++;
];
buttons={};
row=1;
While[row<=Dimensions[names][[1]],
butrow={};
but=1;
While[ but<=Dimensions[names[[row]]][[1]],
button:=
Button[Column[{name},Alignment->{Center}],
NotebookFind[SelectedNotebook[],"unique_vectorfieldspecs",Next,CellTags,AutoScroll->False];
SelectionMove[SelectedNotebook[],All,CellContents,AutoScroll->False];
NotebookWrite[SelectedNotebook[],text,AutoScroll->False];
(*,ImageSize->{400,30}*)];
button=Style[button,ControlsRendering->"Generic"] /.text->texts[[row]][[but]]/.name->names[[row]][[but]];
butrow=Join[butrow,{button}];
(*butrow=Join[butrow,{Row[{" "}],button}];*)
but++;
];
buttons=Join[buttons,{butrow}];
row++;
];
TableForm [ buttons,TableAlignments->{Left,Top},TableSpacing->{1,1}]
(*Column[buttons]*)
];
matrixButtons[3]

vectorfieldspecs=Tuples[Table[discSimplex[vfgridpoints,dim[[p]]],{p,pop}]];

3D graph of the contour function

draw3dcontourfunction = 1 creates a 3d plot of the contour function. 
shiftpercentage specifies the distance above the simplex that the lowest point on the contour function will appear, as a percentage of the distance from the minimum value of the function to the maximum value of the function.
grayline = 1 draws a gray line on the edges of the 3d plot.
heightratio specifies the ratios of the three coordinates in the graphical output.  The default setting is {1, 1, 1.25}.
  
vpoint determines the value for the Mathematica command ViewPoint and hence the point from which the 3d plot is observed.  In Mathematica, ViewPoint -> {x, y, z} gives the position of the view point relative to the center of the three-dimensional box that contains the object being plotted.  The view point is given in a special scaled coordinate system in which the longest side of the bounding box has length 1.  (When heightratio is at its default setting of {1, 1, 1.25}, the longest side corresponds to the z coordinate.  The center of the bounding box is taken to have coordinates {0, 0, 0}. Choosing ViewPoint further away from the object reduces the distortion associated with perspective.

draw3dcontourfunction=0;
grayline=0 ;
shiftpercentage=.2 ;
heightratio={1,1,1.25} ;
vpoint={-2.730, -1.673, 1.096} ;

Program Core

time1=AbsoluteTime[] ;

Solving Numerically

Rest points and stability test

Module[{grid,step,accuracy},
step=.05;
accuracy=5;
grid=Chop[Tuples[Table[discSimplex[1/step,dim[[p]]],{p,pop}]],10^(-accuracy)];
tolLevel=restTolerance Max[Flatten[Table[Speed[grid[[j]]],{j,Dimensions[grid][[1]]}]]]
];

If[findrestpoints== 1,
RPs=restPoints[dyn,Fmp,restPointSearchStartingValues,tolLevel];
CRPs=stabilitySets[RPs,dyn,Fmp, Which[SCharacterization[dyn]=="Smooth","Smooth",SCharacterization[dyn]=="ESS","ESSNumeric",True,"Numeric"],stabn,stabT, stabdelta, stabeps];
SRPs=CRPs[[1]];
USRPs=CRPs[[2]];
If[SCharacterization[dyn]=="Smooth",
SRPvs=Table[Transpose[smoothStabilityValues[SRPs[[j]],dyn,Fmp]],{j,Dimensions[SRPs][[1]]}];
USRPvs=Table[Transpose[smoothStabilityValues[USRPs[[j]],dyn,Fmp]],{j,Dimensions[USRPs][[1]]}];
];
];
NEa=nashEquilibria[Fmp];
If[And[findnashequilibria==1,findESSs==1],
QEa= DeleteCases[Table[If[quasiStrictEqQ[NEa[[j]],Fmp] ,NEa[[j]],"discard"],{j,Dimensions[NEa][[1]]}],"discard"];
ESSs=stabilitySets[QEa,dyn,Fmp,"ESS",stabn,stabT,stabdelta,stabeps][[1]];
ESSs=regularESSs[NEa,Fmp][[1]];
];

Trajectories for Basins of Attraction

If[findrestpoints== 1 &&  drawbasinboundaries==1,
Module[{dta,timetopass,trajectories,trajectorySpecs},
trajectories=Join[
basinsTrajectories[USRPs,dyn,Fmp,AOCharacterization[dyn]],
basinsBoundaryTrajectories[USRPs,dyn,Fmp,AOCharacterization[dyn]]
];
trajectorySpecs =Table[
 {trajectories[[j,1]],timetopass, 0,{GrayLevel[0]}, {}, {}},{j,Dimensions[trajectories][[1]]}
];
customtrajectoryspecs=Join[customtrajectoryspecs,
eliminateDuplicates[.1,trajectorySpecs,Norm,Table[
Flatten[Join[trajectories[[j,1]],trajectories[[j,3]]]]
,{j,Dimensions[trajectories][[1]]}]];
];
];
];

Regular Trajectories

trajSol=Table[difEqSol[trajectoryspecs[[j,1]],dyn,Fmp,0,trajectoryspecs[[j,2]],"initial",AOCharacterization[dyn]],{j,Dimensions[trajectoryspecs][[1]]}];

Special Trajectories

cTrajSol=Table[difEqSol[customtrajectoryspecs[[j,1]],dyn,Fmp,customtrajectoryspecs[[j,3]],customtrajectoryspecs[[j,2]],"final",AOCharacterization[dyn]],{j,Dimensions[customtrajectoryspecs][[1]]}];

Drawing

If[drawvectorfield==0,
drawvf[1]=0;
drawvf[2]=0;
];
If[drawphasediagram==0,pdcontourplot=0] ;

Calculating color scale/gray scale

The approximate minimum and maximum values of the contour function are computed here.  These values are then used to assign colors/gray levels optimally.

If[pdcontourplot==1 ||vfcontourplot==1||draw3dcontourfunction==1 ,
Module[{contourFunctionValues,n},
n=10;
contourFunctionValues=Tuples[Table[discSimplex[n,dim[[p]]],{p,pop}]];
contourFunctionValues=
Table[Chop[contourFunction[contourFunctionValues[[j]]],(10^(-5))],{j,Dimensions[contourFunctionValues][[1]]}];
contourFunctionValues=DeleteCases[contourFunctionValues,Indeterminate];
If[contourFunction === Speed,
minAttained=0,
minAttained=N[Min[contourFunctionValues]]
];
maxAttained=N[Max[contourFunctionValues]];
rangeAttained=shiftpercentage (maxAttained-minAttained);
conColFunAux[Ψ_]=1-Max[0,Min[((Ψ-minAttained)/(maxAttained-minAttained)),1]] ;
If[pdcolor==0,conColFun[Ψ_]=GrayLevel[((.95)/1.45)(1-conColFunAux[conbrightness Ψ]+.45)];
,
conColFun[Ψ_]=Hue[ Which[conColFunAux[Ψ]≤(1/2),(2/5)conColFunAux[Ψ],(1/2)<conColFunAux[Ψ]≤(2/3), (6/5)conColFunAux[Ψ]-(2/5),(2/3)<conColFunAux[Ψ],(3/4)conColFunAux[Ψ]-(1/10)]];
] ;
If[(pdcontourplot==1 || vfcontourplot==1)&& colorbar==1,
colorBarPosition={1.05,1.1};
simplexColorBar={Graphics[(If[pdcolor==1,
{Hue[Which[1-#1 ≤(1/2),(2/5)(1-#1),(1/2)<(1-#1)≤(2/3), (6/5)(1-#1)-(2/5),(2/3)<(1-#1),(3/4)(1-#1)-(1/10)]],
(* ColorData["VisibleSpectrum"][450+#1*(650-450)], *)
Rectangle[{colorBarPosition[[1]],((#1+1.5)/3)},{colorBarPosition[[2]],((#1 +1.51)/3)}]}
,
{GrayLevel[((.95)/1.45)((#1)+.45)],Rectangle[{colorBarPosition[[1]],((#1+1.5)/3)},
{colorBarPosition[[2]],((#1 +1.51)/3)}]}
]&)/@Range[0,1,.01]], Graphics[{Text[PaddedForm[ minAttained,{3,2}],{colorBarPosition[[2]]+.05,.5}], Text[PaddedForm[maxAttained,{3,2}],{colorBarPosition[[2]]+.05,.83}]}]};
,
simplexColorBar=Graphics[{}];
];
If[pdcontourplot==1,pdColorBar=simplexColorBar,pdColorBar=Graphics[{}]] ;
If[vfcontourplot==1,vfColorBar=simplexColorBar,vfColorBar=Graphics[{}]] ;
]];

Drawing Frames and Dots

The cells that belong to this group are for drawing the simplex frame and strategy names.

simplex=Graphics[{EdgeForm[Thickness[.004]],FaceForm[],Polygon[planeCorners]}];

simplexStrategies=Flatten[Table[Table[
Text[strategy[p][j],simplexStrategiesPositions[[p,j]]]
,{j,dim[[p]]}],{p,pop}],1];

This part draws the dots for the rest points.  restPoints[] and stabilitySets[] should have been called in order to determine the rest points and their stability.  You may need to modify this part if  the stability test does not work for your choice of dynamic.

If[And[findrestpoints==1,drawrestpoints==1],
srpDots=Graphics[Table[{GrayLevel[srpcolor],Disk[simplexToPlane[SRPs[[j]]],dotradius]},{j,Dimensions[SRPs][[1]]}]] ;
urpDots=Graphics[Table[{GrayLevel[urpcolor],EdgeForm[Thickness[.4*dotradius]],Disk[simplexToPlane[USRPs[[j]]],.8*dotradius]},{j,Dimensions[USRPs][[1]]}]];
,
srpDots=Graphics[{}];
urpDots=Graphics[{}];
];
If[And[findnashequilibria== 1,drawnashequilibria==1],
neDots=
Graphics[Table[
{GrayLevel[necolor],EdgeForm[Thickness[.4*dotradius]],Disk[simplexToPlane[NEa[[j]]],.8*dotradius]}
,{j,Dimensions[NEa][[1]]}]] ;
,
neDots=Graphics[{}];
];

These cells draw custom text, lines and dots

cDots=If[customdots=={},Graphics[{}],
Graphics[DeleteCases[Table[If[customdots[[j,4]]==0,{customdots[[j,3]],Disk[simplexToPlane[customdots[[j,1]]], customdots[[j,2]]]}],{j,Dimensions[customdots][[1]]}],Null]]
];
cWDots=If[customdots=={},Graphics[{}],
Graphics[DeleteCases[Table[If[customdots[[j,4]]==1,{GrayLevel[1],EdgeForm[Thickness[.4*customdots[[j,2]]]],Disk[simplexToPlane[customdots[[j,1]]],.8 customdots[[j,2]]]}],{j,Dimensions[customdots][[1]]}],Null]]
];
cLines=If[customlines=={},Graphics[{}],
Graphics[DeleteCases[Table[{customlines[[j,3]],customlines[[j,4]],customlines[[j,5]],{Arrowheads[{{customlines[[j,6]], 4/5}}],Arrow[{simplexToPlane[customlines[[j,1]]],simplexToPlane[customlines[[j,2]]]}]}},{j,Dimensions[customlines][[1]]}],Null]]
];
cText=If[customtext=={},Graphics[{}],
Graphics[DeleteCases[Table[{
Text[Style[customtext[[j,1]],
FontFamily->customtext[[j,3]],
FontSlant->customtext[[j,4]],
FontSize->customtext[[j,5]]],N[customtext[[j,2]]]]
},{j,Dimensions[customtext][[1]]}],Null]]
];

These cells draw the frame, lines, dots and text for the 3D plot of the contour function

If[draw3dcontourfunction==1,
graphf[x_] := Append[simplexToPlane[x],contourFunction[x]];
graphmin[x_] := Append[simplexToPlane[x],minAttained-rangeAttained];
tdSimplex=Graphics3D[{FaceForm[],
Polygon[Table[Append[planeCorners[[j]],minAttained-rangeAttained],{j,Dimensions[planeCorners][[1]]}]]}
];
tdStrategies=Graphics3D[
Flatten[Table[Table[
Text[
Style[strategy[p][j],FontFamily->font,FontSize->10],
Append[simplexStrategiesPositions[[p,j]],minAttained-rangeAttained]
]
,{j,dim[[p]]}],{p,pop}],1]];
If[findrestpoints==1,
tdLines=Graphics3D[{Dashing[{.01,.01}],Flatten[Table[Line[{graphmin[RPs[[j]]],graphf[RPs[[j]]]}],{j,Dimensions[RPs][[1]],1}]]}];
tdPoints=Graphics3D[{PointSize[.02],Flatten[Table[Point[graphmin[SRPs[[j]]]],{j,Dimensions[SRPs][[1]]}]],GrayLevel[.5],Flatten[Table[Point[graphmin[USRPs[[j]]]],{j,Dimensions[USRPs][[1]]}]]}];
,
tdLines=Graphics3D[{}];
tdPoints=Graphics3D[{}];
];
,
tdSimplex=Graphics3D[{}];
tdStrategies=Graphics3D[{}];
tdLines=Graphics3D[{}];
tdPoints=Graphics3D[{}];
];

Drawing vector fields on the simplex

If[drawvectorfield==1 && (drawvf[1]==1 || drawvf[2]==1), 
Do[
If[drawvf[f]==1,
If[findrestpoints==1 && arrowsfromrestpoints==1,
vectorField[f]=customVectorField [ 
Join[
Table[simplexToPlane[vectorfieldspecs[[j]]],{j,Dimensions[vectorfieldspecs][[1]] }],
Table[simplexToPlane[RPs[[j]]],{j,Dimensions[RPs][[1]] }]
],
Join[
Table[simplexToPlane[vectorfieldf[f][vectorfieldspecs[[j]],Fmp]]-{1,0},
{j,Dimensions[vectorfieldspecs][[1]]}],
Table[simplexToPlane[vectorfieldf[f][RPs[[j]],Fmp]]-{1,0},
{j,Dimensions[RPs][[1]]}]
],vfcolor[f],vfscale[f],vfarrowheadsize[f],vfgridpoints];,
vectorField[f]=customVectorField [
Table[simplexToPlane[vectorfieldspecs[[j]]],{j,Dimensions[vectorfieldspecs][[1]] }],
Table[simplexToPlane[vectorfieldf[f][vectorfieldspecs[[j]],Fmp]]-{1,0},
{j,Dimensions[vectorfieldspecs][[1]]}],
vfcolor[f],vfscale[f],vfarrowheadsize[f],vfgridpoints];
];
,
vectorField[f]=Graphics[{}];
];
,{f,2}];
,
Do[vectorField[f]=Graphics[{}],{f,2}];
];

Drawing Solutions to Diff Eqs

trajLines=If[Dimensions[trajectoryspecs][[1]]==0,
Graphics[{}],Table[ParametricPlot[Evaluate[simplexToPlane[X[t]]/.trajSol[[j]]],{t,0,trajectoryspecs[[j,2]]}, DisplayFunction -> Identity,PlotPoints-> ppoints, Axes->{False,False},
PlotRange->simplexPlotRange,PlotStyle-> {trajectoryspecs[[j,3]]},AspectRatio->Automatic],{j,Dimensions[trajectoryspecs][[1]]}]];

cTrajLines=If[Dimensions[customtrajectoryspecs][[1]]==0,
Graphics[{}],
Table[ParametricPlot[Evaluate[simplexToPlane[X[t]]/.cTrajSol[[j]]],{t,customtrajectoryspecs[[j,3]],customtrajectoryspecs[[j,2]]}, DisplayFunction -> Identity,PlotPoints-> ppoints, Axes->{False,False},
PlotRange->simplexPlotRange,
PlotStyle-> {customtrajectoryspecs[[j,4]]},AspectRatio->Automatic],{j,Dimensions[customtrajectoryspecs][[1]]}]] ;

Drawing Arrows

invertarrows=0;
If[invertarrows==1,auxrelarrowstart=-1*relarrowstart,auxrelarrowstart=relarrowstart];
trajAr=If[Dimensions[trajectoryspecs][[1]]==0,
Graphics[{}],
Graphics[Flatten[Table[{
If[Dimensions[trajectoryspecs[[k,3]]][[1]]≠ 0, Join[Select[trajectoryspecs[[k,3]],(MatchQ[#1,Hue[__]] ||MatchQ[#1,GrayLevel[__]]||MatchQ[#1,RGBColor[__]]||MatchQ[#1,CMYKColor[__]])&],{GrayLevel[0]}][[1]],GrayLevel[0] ],
Table[{
Arrowheads[{{trajectoryspecs[[k,5,j]],arrowheadpos}}],
Arrow[{
Flatten[Evaluate[
simplexToPlane[X[trajectoryspecs[[k,4,j]]*(1-auxrelarrowstart)]]/.trajSol[[k]]]],
Flatten[Evaluate[
simplexToPlane[X[trajectoryspecs[[k,4,j]]]]/.trajSol[[k]]]]
}]
},{j,Dimensions[trajectoryspecs[[k,5]]][[1]]}]},{k,Dimensions[trajectoryspecs][[1]]}]]]
];

invertarrows=0;
If[invertarrows==1,auxrelarrowstart=-1*relarrowstart,auxrelarrowstart=relarrowstart];
cTrajAr=If[Dimensions[customtrajectoryspecs][[1]]==0,
Graphics[{}],
Graphics[Flatten[Table[{
If[Dimensions[customtrajectoryspecs[[k,4]]][[1]]≠ 0,Join[Select[customtrajectoryspecs[[k,4]],(MatchQ[#1,Hue[__]] ||MatchQ[#1,GrayLevel[__]]||MatchQ[#1,RGBColor[__]]||MatchQ[#1,CMYKColor[__]])&],{GrayLevel[0]}][[1]],GrayLevel[0]],
Table[{
Arrowheads[{{customtrajectoryspecs[[k,6,j]],arrowheadpos}}],Arrow[{
Flatten[Evaluate[
simplexToPlane[
X[customtrajectoryspecs[[k,5,j]]*(1-auxrelarrowstart )]]/.cTrajSol[[k]]]],Flatten[Evaluate[
simplexToPlane[
X[customtrajectoryspecs[[k,5,j]]]]/.cTrajSol[[k]]]]
}]
},{j,Dimensions[customtrajectoryspecs[[k,6]]][[1]]}]}
,{k,Dimensions[customtrajectoryspecs][[1]]}]]]];

Drawing Contour Plot

contourRanges={{0,1},{0,1}};

If[Or[And[drawphasediagram==1,pdcontourplot==1],And[drawvectorfield==1,vfcontourplot==1]],
Module[{f},
f[{x_,y_}]:=Evaluate[contourFunction[planeToSimplex[{x,y}]]];
con=ContourPlot[f[{x,y}],{x,contourRanges[[1,1]],contourRanges[[1,2]]},{y,contourRanges[[2,1]],contourRanges[[2,2]]},
DisplayFunction->Identity,
PlotPoints->plotprecision,
Contours->customContours[numberofcontours,50,f,contourRanges],
Frame->None,AspectRatio->Automatic,ContourStyle->None,Exclusions->None,ColorFunctionScaling->False,ColorFunction->(conColFun[#1]&)]
];
con=con/.contourPlotRule;
If[compressgraphic==1,con=fixPolygons[con]];
,
con=Graphics[{}];
];
If[pdcontourplot==1,pdCon=con,pdCon=Graphics[{}]];
If[vfcontourplot==1,vfCon=con,vfCon=Graphics[{}]];

Drawing 3D Plot of contour function

If[draw3dcontourfunction==1, 
Module[{restricted,f,tdEdgef},
restricted[x_]:=If[Apply[And,Table[simplexQ[x[[p]]],{p,pop}]],contourFunction[x],-Infinity];
f[{x_,y_}]:=Evaluate[restricted[planeToSimplex[{x,y}]]]; tdCon= Plot3D[ f[{x,y}],
{x,contourRanges[[1,1]],contourRanges[[1,2]]},{y,contourRanges[[2,1]],contourRanges[[2,2]]},
DisplayFunction -> Identity,
PlotPoints -> plotprecision, PlotRange -> {minAttained-rangeAttained,maxAttained+.05}, ViewPoint->vpoint,ColorFunctionScaling-> False , ColorFunction-> (conColFun[#3]&), 
Axes -> False, Boxed -> False, Mesh -> False, BoxRatios -> {1,1,1.25}] ;
If[grayline==1,
tdEdgef[x_,y_]:=ParametricPlot3D[ {x[t],y[t],contourFunction[planeToSimplex[{x[t],y[t]}]]},{t,0,1},
PlotStyle->{GrayLevel[.3],Thickness[.006]},
DisplayFunction->Identity,
PlotPoints->200, PlotRange->{Full,Full,{minAttained-rangeAttained,maxAttained+.05}}, ViewPoint->vpoint, 
Axes->False,Boxed->False,BoxRatios->{1,1,1.25}];
tdEdge={
tdEdgef[Function[t,0],Function[t,t]],tdEdgef[Function[t,t],Function[t,0]],
tdEdgef[Function[t,1],Function[t,t]],
tdEdgef[Function[t,t],Function[t,1]]
};
,
tdEdge=Graphics3D[{}];
];
];
,
tdCon=Graphics3D[{}];
tdEdge=Graphics3D[{}];
];

Generating the output notebook

Redirection

OutputToOneNotebook=0;
OutputClearNotebook=1;

If[$VersionNumber==8,
Module[{j,flag,str},
flag=True;
For[j=1,j<=Dimensions[ScheduledTasks[]][[1]],j++,
str=ToString[ScheduledTasks[][[j,2]]];
If[StringTake[str,Min[11,StringLength[str]]]=="dynamotask",
flag=False;
Break[]];
];
If[flag,RunScheduledTask[
dynamotask;
(* If[CurrentValue[Notebooks[][[1]],"NotebookFileName"]!="Messages", ...];*)
Quiet[Module[{nb,name},
nb=SelectedNotebook[];
name=CurrentValue[nb,"NotebookFileName"];
If[ToString[name]!="$Failed",
If[And[Or[StringTake[name,Min[13,StringLength[name]]]!="dynamo_output",StringTake[name,-Min[3,StringLength[name]]]==".nb"],ToString[Options[nb,WindowTitle][[1,2]]]≠"Automatic"],
SetOptions[nb,WindowTitle->Automatic]
(*NotebookSave[nb];*)
];
];
]];
]];
];
];
(* With the technology to keep track of all the dynamo_notebooks the following should be rewritten; only v8 and above. The use of Notebooks[] crashes mathematica when saving figures who knows why. *)

(*----------Redirecting output to a new notebook-------*)

Which[
Notebooks["dynamo_output"]=={} && OutputToOneNotebook==1,
$nb=CreateWindow[WindowSize-> {650,550}, WindowTitle->"dynamo_output"];
(*SetOptions[$nb,WindowTitle->Dynamic[If[ToString[Quiet[NotebookFileName[$nb]]]!="$Failed",Automatic,"dynamo_output"]]];
*)
,
Not[Notebooks["dynamo_output"]=={} ] && OutputToOneNotebook==1,
$nb=Notebooks["dynamo_output"][[1]];
SetSelectedNotebook[$nb];
If[OutputClearNotebook==1,
FrontEndExecute[FrontEndToken["SelectAll"]];
NotebookDelete[$nb];
(*NotebookPut[Notebook[{}],$nb]*)
];
,
OutputToOneNotebook==0,
SetNotebookCounter[]:=Module[{nb,NotebookCounter,notebook,namelength,tmpstring,tmpcounter},
NotebookCounter=1;
For[notebook=1,notebook<=Dimensions[Notebooks[]][[1]],notebook++,
nb=Notebooks[][[notebook]];
namelength=StringLength[CurrentValue[nb,"NotebookFileName"]];
If[namelength>14,
tmpstring=StringTake[CurrentValue[nb,"NotebookFileName"],14];
tmpcounter=ToExpression[StringTake[CurrentValue[nb,"NotebookFileName"],{15,namelength}]];
tmpcounter =If[IntegerQ[tmpcounter],tmpcounter,0];
(* Print[tmpstring,tmpcounter]; *)
If[ tmpstring=="dynamo_output_" && Not[tmpcounter <=0 ] && tmpcounter>NotebookCounter,NotebookCounter=tmpcounter];
];
];
NotebookCounter++;
If[NotebookCounter==2 && Notebooks["dynamo_output_1"]=={},NotebookCounter=1];
(*
While[
Not[Notebooks["dynamo_output_"~~ToString[NotebookCounter]]=={}],
NotebookCounter=NotebookCounter+1;
];
*)
NotebookCounter
];
$nb=CreateWindow[WindowSize-> {650,550},WindowTitle->"dynamo_output_"~~ToString[SetNotebookCounter[]]];
];
$DisplayFunction=NotebookWrite[$nb,Cell[BoxData[MakeBoxes[#,StandardForm]],"Output"]]& ;
(* 
$DisplayFunction=NotebookWrite[$nb, Cell[GraphicsData["PostScript", ExportString[#,"EPS"]],"Graphics"]]& ;
*)
$Post=If[#=!= Null,NotebookWrite[$nb,Cell[BoxData[MakeBoxes[#,StandardForm]],"Output"]]]& ;
Unprotect[Print] ;
Print[x_] := NotebookWrite[$nb,Cell[BoxData[MakeBoxes[x,StandardForm]]]] ;
Print[x__] :=Map[Print,{x}] ;
Protect[Print] ;
(*----------------Redirecting code ends----------------*)

SelectionMove[$nb,After,Notebook];
NotebookWrite[$nb,Cell[TextData["Dynamo 2x2 report"],"Title", FontFamily->"Palatino",FontColor-> GrayLevel[.7],Background-> Hue[RandomReal[],.4],TextAlignment-> Center,CellFrame-> False]] ;

Text output - the game

Print[Style[DisplayForm[GridBox[{{""},{"The Game"}}]],"Section"]];

FormatNumberSingle[x_,maxdenominator_:20,precision_:5]:=Module[{aux1,aux2},
aux1=N[x,precision];
aux2=Rationalize[aux1,10^(-precision)];
If[Denominator[aux2]>maxdenominator,aux1,aux2]
];
FormatNumber[x_,maxdenominator_:20,precision_:5]:=Module[{f,arg,d,aux},
f[arg_]:=If[Depth[arg]==1,FormatNumberSingle[arg,maxdenominator,precision],arg];
d=Depth[x];
aux=Apart[x];
While[d>0,
aux=Map[f,aux,{d}];
d--;
];
aux
];
FormatState[x_,maxdenominator_:20,precision_:5]:=If[Dimensions[x][[1]]==1,
FormatNumber[x[[1]],maxdenominator,precision],
FormatNumber[x,maxdenominator,precision]
];

If[gamekind==NormalGame,Print[Row[{
DisplayForm[RowBox[{"   Normal form game:  "}]],
DisplayForm[GridBox[{{""},{strategy[1][1]},{strategy[1][2]}}]],DisplayForm[GridBox[Join[{{strategy[2][1],strategy[2][2]}},U]]]}
]]];
Print[Row[{DisplayForm[RowBox[{"   Payoff vector field:  F"}]],MatrixForm[{x,y}], DisplayForm[RowBox[{" = "}]],MatrixForm[Flatten[Fmp[{{x1,x2},{y1,y2}}]]]}]] ;

If[findnashequilibria==1,
If[Dimensions[NEa][[1]]==1,Print[DisplayForm[RowBox[{"   Nash equilibrium:"}]]],Print[DisplayForm[RowBox[{"   Nash equilibria:"}]]]] ;
Print[Row[{
DisplayForm[RowBox[{"     "}]],
TableForm[Table[MatrixForm[{FormatState[NEa[[j]],maxdenominator]}],{j,Dimensions[NEa][[1]]}]]
}]]
If[findESSs==1,
Which[
Dimensions[ESSs][[1]]==0,Print[DisplayForm[RowBox[{"   No regular ESSs found"}]]],
Dimensions[ESSs][[1]]==1,Print[DisplayForm[RowBox[{"   Regular ESS:"}]]],
True,Print[DisplayForm[RowBox[{"   Regular ESSs:"}]]]
] ;
If[Dimensions[ESSs][[1]]>0,
Print[Row[{
DisplayForm[RowBox[{"     "}]],
TableForm[Table[MatrixForm[{FormatState[ESSs[[j]],maxdenominator]}],{j,Dimensions[ESSs][[1]]}]]
}]]
]
];
];

Text output - the dynamic

Print[Style[DisplayForm[GridBox[{{""},{"The Dynamic"}}]],"Section"]];

Print[Row[{DisplayForm[RowBox[{"   Name of dynamic:  "}]],  dyn}]];

Print[DisplayForm[RowBox[{"   Law of motion: "}]], Row[{
DisplayForm[RowBox[{"   (Subscript[V, F])"}]],MatrixForm[{x,y}], DisplayForm[RowBox[{" = "}]],MatrixForm[Simplify[Flatten[dyn[{{x1,x2},{y1,y2}},Fmp]]]] 
}]];

Text output - stability of rest points

If[findrestpoints==1,
If[Dimensions[SRPs][[1]]==1,
Print[DisplayForm[RowBox[{"   Stable rest point:"}]]],Print[DisplayForm[RowBox[{"   Stable rest points:"}]]]
];
If[Dimensions[SRPs][[1]]==0,Print[DisplayForm[RowBox[{"    None"}]]],
If[SCharacterization[dyn]=="Smooth",
Print[DisplayForm[RowBox[{"   (state x, relevant eigenvalues of (Subscript[DV, F])[x], relevant eigenvectors of (Subscript[DV, F])[x])"}]]];
Print[Row[{
DisplayForm[RowBox[{"   "}]],
TableForm[Table[{MatrixForm[{FormatState[SRPs[[j]],maxdenominator]}],MatrixForm[{FormatState[SRPvs[[j,1]],maxdenominator]}],MatrixForm[{FormatState[SRPvs[[j,2]],maxdenominator]}]},{j,Dimensions[SRPs][[1]]}],TableAlignments->Top]
} ]];
,
Print[Row[{
DisplayForm[RowBox[{"   "}]],
TableForm[Table[MatrixForm[{FormatState[SRPs[[j]],maxdenominator]}],{j,Dimensions[SRPs][[1]]}]]
} ]];
];
];
If[Dimensions[USRPs][[1]]==1,
Print[DisplayForm[RowBox[{"   Unstable rest point:"}]]],
Print[DisplayForm[RowBox[{"   Unstable rest points:"}]]]] ;

If[Dimensions[USRPs][[1]]==0,
Print[DisplayForm[RowBox[{"    None"}]]],
If[SCharacterization[dyn]=="Smooth",
Print[DisplayForm[RowBox[{"   (state x, relevant eigenvalues of (Subscript[DV, F])[x], relevant eigenvectors of (Subscript[DV, F])[x])"}]]];
Print[Row[{
DisplayForm[RowBox[{"   "}]],
TableForm[Table[{MatrixForm[{FormatState[USRPs[[j]],maxdenominator]}],MatrixForm[{FormatState[USRPvs[[j,1]],maxdenominator]}],MatrixForm[{FormatState[USRPvs[[j,2]],maxdenominator]}]},{j,Dimensions[USRPs][[1]]}],TableAlignments->Top]
} ]];
,
Print[Row[{
DisplayForm[RowBox[{"   "}]],
TableForm[Table[MatrixForm[{FormatState[USRPs[[j]],maxdenominator]}],{j,Dimensions[USRPs][[1]]}]]
}]];
];
];
]

Text output - the contour function

If[pdcontourplot==1||  draw3dcontourfunction==1,Print[Style[
DisplayForm[GridBox[{{""},{"Graphics"}}]],
"Section"]] 
];

If[pdcontourplot==1|| draw3dcontourfunction==1,Print[Row[{ 
DisplayForm[RowBox[{"   Name of Contour Function: "}]], 
contourFunction
}]]
]; 
If[showcontourformula==1,Print[
DisplayForm[RowBox[{"   Definition of Contour Function: "}]],
Row[{
DisplayForm[RowBox[{"   ψ"}]],
MatrixForm[ {x,y}],
DisplayForm[RowBox[{" = "}]],
contourFunction[{{x1,x2},{y1,y2}}] 
}]
]];
Print[Row[{
DisplayForm[RowBox[{"   Minimum Value = "}]],
N[minAttained]
}]];
Print[Row[{
DisplayForm[RowBox[{"   Maximum Value = "}]],
N[maxAttained]
}]] ;

Generating the vector field diagram

If[drawvectorfield==1, 
Do[Return[
Show[{vfCon,vfColorBar, simplex,vectorField[f],neDots,srpDots,urpDots},
ImageSize-> {500,500},
PlotRange->simplexPlotRange,
AspectRatio-> Automatic,BaseStyle->{FontFamily->font,FontSize->10},
Epilog->simplexStrategies]],{f,2}]]

Generating the phase diagram on simplex

If[drawphasediagram==1,
 Show[{pdCon,pdColorBar,simplex,trajLines,cTrajLines,cLines,trajAr,cTrajAr,neDots,srpDots,urpDots,cDots,cWDots,cText},
ImageSize-> {500,500},
PlotRange->simplexPlotRange,
AspectRatio-> Automatic,BaseStyle->{FontFamily->font,FontSize->10},
Epilog->simplexStrategies]]

Generating the 3D plot of the contour function

If[draw3dcontourfunction== 1, Show[{tdCon,tdSimplex,tdLines,tdPoints,tdEdge,tdStrategies},ImageSize->{500,500},ImagePadding->30]]

Printing time elapsed

PrintTimeElapsed[t_]:=Module[{dif},
dif=DateList[AbsoluteTime[]-t];
Print[Row[{
DisplayForm[RowBox[{"Runtime, "}]],
If[dif[[5]]<10,Row[{DisplayForm[RowBox[{"0"}]],dif[[5]]}],dif[[5]]],
DisplayForm[RowBox[{":"}]],
If[Floor[dif[[6]]]<10,Row[{DisplayForm[RowBox[{"0"}]],Floor[dif[[6]]]}],Floor[dif[[6]]]] 
}]
];
];
PrintTimeElapsed[time1]
Print[Row[{}]];

End redirection

(*----------Returning to original output production-------*)

$DisplayFunction=Identity;
$Post=If[#1=!=Null,CellPrint[Cell[BoxData[MakeBoxes[#1,StandardForm]],"Output"]]]&;
Unprotect[Print];
Print[x_]:=CellPrint[Cell[BoxData[MakeBoxes[x,StandardForm]]]];
Print[x__]:=Print/@{x};
Protect[Print];

(*----------------code ends----------------*)

Procedures; Maintenance

Procedures

Preliminary definitions and calculations

Norm

dim={2,2};
NormalGame=1;
LinearGame=2;
CongestionGame=3;

pop=Dimensions[dim][[1]];
stateNorm[x_,p_:2]:=Norm[Flatten[x],p];stateNorm[x_,p_:2]:=Norm[Table[Norm[Drop[x[[j]],-1],p],{j,Dimensions[x][[1]]}],p];
stateNorm[x_,p_:2]:=Norm[Table[Norm[x[[j]],p],{j,Dimensions[x][[1]]}],p];

General Purpose Routines

emptyVector[name_,dim_,var_:""]:=Module[{j,x},
x={};
For[j=1,j<=dim,j++,
(* x=Join[x,{ToExpression[ToString[Row[{name,j,var}]]]}]; *)
x=Join[x,{ToExpression[StringJoin[ToString[name],ToString[j],ToString[var]]]}];
];
x
];
(* Print[emptyVector[y,3,"_"]]; *)

vectorReplace[orig_,replace_,positions_]:=Module[{j,counter,result},
result={};
counter=1;
For[j=1,j<=Dimensions[orig][[1]],j++,
If[positions[[j]]==1,
result=Join[result,{replace[[counter]]}];
counter++;,
result=Join[result,{orig[[j]]}];
];
];
result
];
(* Print[vectorReplace[{1,1,1},{2,2},{1,0,1}]]; *)

eliminateDuplicates[delta_,reallist_,norm_:Norm,complist_:{}]:=Module[{j,k,reducedreallist,reducedcomplist},
If[complist=={},reducedcomplist=reallist,reducedcomplist=complist];
reducedreallist=reallist;
j=1;
While[j<=Dimensions[reducedcomplist][[1]],
k=1;
While[k<=Dimensions[reducedcomplist][[1]],
If[Not[j==k] && norm[reducedcomplist[[j]]-reducedcomplist[[k]]]<delta,
reducedcomplist=Drop[reducedcomplist,{k}];
reducedreallist=Drop[reducedreallist,{k}];
(* Print[j,k]; *)
Continue[];
];
k++;
];
j++
];
reducedreallist
];
(* eliminateDuplicates[2,{{6,6},{1,1},{1.1,.9},{5,7},{10,20},{1.1,.95}}] *)

subMatrix[matrix_,positions_]:=Module[{j,k,dim,result},
dim=Dimensions[positions][[1]];
result=matrix;
j=1;
result=Table[DeleteCases[Table[If[positions[[k]]==1,result[[j,k]],"discard"],{k,dim}],"discard"],{j,dim}];
result=DeleteCases[Table[If[positions[[j]]==1,result[[j]],"discard"],{j,dim}],"discard"];
result
];
(* subMatrix[{{1,3,4},{2,5,6},{7,9,8}},{1,0,1}] *)

simplexQ[x_,accuracy_:8]:=If[Max[x]>1||Min[x]<0||Abs[Total[x]-1]>10^(-accuracy),False,True];

checkBoundary[x_,delta_]:=Module[{j,boundaryvector},
boundaryvector={};
For[j=1,j<=Dimensions[x][[1]],j++,
If[Abs[x[[j]]]<delta,
boundaryvector=Append[boundaryvector,1];,
boundaryvector=Append[boundaryvector,0];
];
];
boundaryvector
];

sampleAround[x_,eps_,n_:1,variant_:"smallerthan"]:=Module[{validpoints,j,point},
validpoints={};
For[j=1,j<=n,j++,
point = RandomReal[1,Dimensions[x][[1]]];
point=point-Mean[point];
If[variant=="smallerthan",
point =RandomReal[{0,eps}]*point/Norm[point];
,
point =eps*point/Norm[point];
];
point =x+point;
If[simplexQ[point],
validpoints =Append[validpoints,point];
,
j--;
];
];
validpoints
];

numberSystem[npositions_,npossibilities_,j_]:=Module[{leftover,x},
leftover=j-1;
x=ConstantArray[0,npositions];
Do[
x[[k]]=IntegerPart[leftover/npossibilities^(npositions-k)];
leftover=leftover-x[[k]]*npossibilities^(npositions-k);
,{k,npositions}];
x
];

discSimplex[n_,dim_,perturbed_:False,eps_:Infinity]:=Module[{j,k,epsilon,x,grid,leftover,nc},
If[perturbed==True && eps>1/(2*dim),epsilon= 1/(2*dim),epsilon=eps];
grid = {};
nc=(n+1)^dim;
For[j=1,j<=nc,j+=n,
x=numberSystem[dim,n+1,j];
If[Total[x]== n,
x=x/n;
If[perturbed==True,x=sampleAround[x,epsilon][[1]]];
grid=Append[grid,x];
];
];
grid
];
(*
ListPointPlot3D[discSimplex[50,3]];
Dimensions[discSimplex[50,3]];
*)

support[x_,accuracy_:8]:=Module[{tolLevel},
tolLevel=10^(-accuracy);
Table[If[Abs[x[[j]]]>tolLevel,1,0],{j,Dimensions[x][[1]]}]
];

simplexSupport[x_,accuracy_:8]:=Module[{tolLevel},
tolLevel=10^(-accuracy);
Table[If[x[[j]]>tolLevel,1,0],{j,Dimensions[x][[1]]}]
];

simplexClosestPoint[x_]:=Module[{j,k,m,dim,nc,boundaries,reduced,npos,flag,project,result},
dim=Dimensions[x][[1]];
result=ConstantArray[0,dim];
 nc=2^dim;
For[j=2,j<=nc,j++,
boundaries=numberSystem[dim,2,j];
npos=1;
reduced=ConstantArray[0,dim];
Do[
If[boundaries[[k]]==1,
reduced[[npos]]=x[[k]];
npos++;
];
,{k,dim}];
npos--;
project=IdentityMatrix[npos]-1/npos*ConstantArray[1,{npos,npos}];
reduced=project.Take[reduced,npos]+1/npos*ConstantArray[1,npos];
If[Min[reduced]<0, Continue[],
flag=False;
Do[
If[And[boundaries[[k]]==0,x[[Position[boundaries,1][[1,1]]]]-x[[k]]-reduced[[1]]<0],
flag=True;
Break[];
];
,{k,dim}];
If[flag,Continue[],
m=1;
Do[
If[boundaries[[k]]==1,
result[[k]]=reduced[[m]];
m++;
];
,{k,dim}];
Break[];
];
];
];
result
];

projectionMatrix[dim_]:=Module[{pop},
pop=Dimensions[dim][[1]];
Flatten[Table[Transpose[Join[
ConstantArray[0,{If[p==1,0,Product[dim[[k]],{k,1,p-1}]],dim[[p]]}],
IdentityMatrix[dim[[p]]]-1/dim[[p]]*ConstantArray[1,{dim[[p]],dim[[p]]}],
ConstantArray[0,{If[p==pop,0,Product[dim[[k]],{k,p+1,pop}]],dim[[p]]}]
]],{p,pop}],1]
];

General Game Theory Routines

bestResponse::multip="Multiple Best Responses Found: Mean Dynamic is a Differential inclussion.";
bestResponse[x_,F_]:=Module[{Fx,y},
Fx=F[x];
Table[
y=ConstantArray[0,Dimensions[x[[p]]][[1]]];
y[[Ordering[Fx[[p]],-1]]]=1;
If[Dimensions[Position[Fx[[p]],Max[Fx[[p]]]]][[1]]>1,Message[bestResponse::multip]];
y,{p,Dimensions[x][[1]]}]
];

bestResponseSet[x_,F_,accuracy_:8]:=Module[{Fx,max},
Fx=Round[F[x],10^(-accuracy)];
Table[max=Max[Fx[[p]]];
Table[If[Fx[[p,j]]==max,1,0],{j,Dimensions[x[[p]]][[1]]}],{p,Dimensions[x][[1]]}]
];

quasiStrictEqQ[x_,F_,accuracy_:8]:=If[Table[Table[If[x[[p,j]]>0,1,0],{j,Dimensions[x[[p]]][[1]]}],{p,Dimensions[x][[1]]}]==bestResponseSet[x,F,accuracy],True,False];

nashEqQ[x_,F_,accuracy_:8]:=Module[{BRS},
BRS=bestResponseSet[x,F,accuracy];
If[Total[Table[If[Total[Table[If[x[[p,j]]>0,If[BRS[[p,j]]==1,1,0],1],{j,Dimensions[x[[p]]][[1]]}]]==Dimensions[x[[p]]][[1]],1,0],{p,Dimensions[x][[1]]}]]==Dimensions[x][[1]],True,False]
];

variationalIneqCheck[x_,support_,Fx_,variant_:"standard"]:=Module[{j,dim,flag},
flag=True;
dim=Dimensions[x][[1]];
For[j=1,j<=dim,j++,
If[support[[j]]==0,
If[variant=="quasistrict",
flag=(x-UnitVector[dim,j]).Fx>0;
,
flag=(x-UnitVector[dim,j]).Fx>=0;
];
If[Not[flag],Break[]];
];
];
flag
];

Differential Equations Routines

Protect[Xt];
pattern={{__?NumberQ}...};
X[t_]:=Table[Module[{result},
result=Table[Xt[p][j][t],{j,dim[[p]]-1}];Append[result,1-Total[result]]],{p,pop}];

scalarDyn[p_?NumericQ,j_?NumericQ,x:pattern,F_]:=dyn[x,F][[p,j]];
difEqSol[x_,dyn_,F_,t0_,T_,condition_:"initial",altopt_:False,accuracy_:10,maxsteps_:10000]:=Module[{cond,method,ms,solvedelayed,d},
Which[
condition=="initial" ,cond=t0;,
condition=="final" ,cond=T;,
True,cond=condition;
];
If[altopt,
method="Automatic";
(*method="ImplicitRungeKutta";*)
(*method="StiffnessSwitching";*)
(*method="ExplicitRungeKutta";*)
ms=10/maxsteps*(T-t0);
solvedelayed=False;
d=Table[Table[scalarDyn[p,j,X[t],F],{j,dim[[p]]}],{p,pop}];
,
method="Automatic";
ms=Automatic;
solvedelayed=False;
d=dyn[X[t],F];
];
NDSolve[
Flatten[Table[Table[{Xt[p][j]'[t]==d[[p,j]],Xt[p][j][cond]==x[[p,j]]},{j,dim[[p]]-1}],{p,pop}],1]
,Flatten[Table[Table[Xt[p][j],{j,dim[[p]]-1}],{p,pop}]],
{t,t0,T},
AccuracyGoal->accuracy,PrecisionGoal->accuracy,
Method->method,
MaxSteps->maxsteps,
MaxStepSize->ms,
(*MaxStepFraction->ms,*)
SolveDelayed->solvedelayed
]
];

difEqAdvance[x_,dyn_,F_,T_,altopt_:False,accuracy_:10,maxsteps_:10000,method_:"Automatic"]:=Module[{maxstepsize,result,r},
(*method="Automatic";*)
(*method="ImplicitRungeKutta";*)
(*method="StiffnessSwitching";*)
result=difEqSol[x,dyn,F,0,T,"initial",altopt,accuracy,maxsteps];
(*
result=NDSolve[
Module[{d},
d=dyn[X[t],F];
Flatten[Table[Table[{Xt[p][j]'[t]==d[[p,j]],Xt[p][j][0]==x[[p,j]]},{j,dim[[p]]-1}],{p,pop}],1]
]
,
Flatten[Table[Table[Xt[p][j],{j,dim[[p]]-1}],{p,pop}]]
,
{t,0,T},
AccuracyGoal->accuracy,PrecisionGoal->Infinity,
Method->method,
MaxSteps->maxsteps
];
*)
result=Evaluate[Flatten[Table[Drop[X[T][[p]],-1],{p,pop}]]/.result][[1]];
Table[
r=result[[
If[p==1,0,Product[dim[[j]]-1,{j,p-1}]]+1;;
Product[dim[[j]]-1,{j,p}]
]];
Append[r,1-Total[r]]
,{p,pop}]
];

Coordinate transformation (for plotting the phase diagram)

simplexToPlane maps points in the simplex to their 2D represetations for graphing.
planeToSimplex does the reverse.

simplexToPlane[x_]:={x[[2,2]],x[[1,1]]};
planeToSimplex[x_]:={{x[[2]],1-x[[2]]},{1-x[[1]],x[[1]]}};

simplexCorners=Tuples[Table[Table[UnitVector[dim[[p]],j],{j,dim[[p]]}],{p,pop}]];
planeCorners=Table[simplexToPlane[simplexCorners[[j]]],{j,Dimensions[simplexCorners][[1]]}];

planeCorners=planeCorners[[{1,2,4,3}]];
simplexPlotRange={{-.08,1.15},{-.08,1.08}};
simplexStrategiesPositions={{{-.025,1},{-.025,0}},{{0,1.025},{1,1.025}}};

Custom Graphics Routines

contourPlotRule=({EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[{r,Thickness[Small]}],r,i});
(*http://pages.uoregon.edu/noeckel/MathematicaGraphics.html*)

customContours[numberofcontours_,extremesgridsize_,function_,range_]:=Module[{values,extremes},
values=Flatten[Table[function[{range[[1,1]]+i/extremesgridsize*(range[[1,2]]-range[[1,1]]),range[[2,1]]+j/extremesgridsize*(range[[2,2]]-range[[2,1]])}],{j,0,extremesgridsize},{i,0,extremesgridsize}]];
extremes={Min[values],Max[values]};
Table[j,{j,extremes[[1]],extremes[[2]],(extremes[[2]]-extremes[[1]])/(numberofcontours -1)}]];

customVectorField[points_,vectors_,color_,magninorm_,magnihead_,gridpoints_]:=Module[{rescaler,maxnorm},
maxnorm=Max[Table[Norm[vectors[[j]] ],{j,Dimensions[vectors][[1]]}]];
rescaler=1/(gridpoints-1);
rescaler=magninorm rescaler/maxnorm;
Return[Table[
If[Norm[vectors[[j]]]>.01*maxnorm,
Graphics[{color,Arrowheads[{{magnihead /(2gridpoints),.99}}],
Arrow[{points[[j]],points[[j]]+vectors[[j]]rescaler}]
}],
Graphics[{color,
Disk[points[[j]],Scaled[0.005]]
}]
],
{j,Dimensions[vectors][[1]]}
]]
];

Rest Points and Equilibria

Calculation of rest points (continuous dynamics)

subFindRoot::inconsistency = "Warning, trivial root searched for a fixed point without any dimension to vary. Test if the isolated point is a root manually.";
subFindRoot[dyn_,F_,x0_,indicator_,accuracy_:8]:=
Module[{j,p,counter,pop,dim,subdim,subdimlast,x,xp,equations,restrictions,root},
pop=Dimensions[indicator][[1]];
dim=Table[Dimensions[indicator[[p]]][[1]],{p,pop}];
subdim=Table[Total[indicator[[p]]],{p,pop}];
subdimlast=Table[Max[Table[If[indicator[[p,j]]==1,j,0],{j,dim[[p]]}]],{p,pop}];
If[Total[subdim]>1,
x=emptyVector["x",Total[subdim]];
x=Table[vectorReplace[
x0[[p]],x[[Sum[subdim[[j]],{j,p-1}]+1;;Sum[subdim[[j]],{j,p}]]],
indicator[[p]]],{p,pop}];
For[p=1,p<=pop,p++,
x[[p,subdimlast[[p]]]]=0;
x[[p,subdimlast[[p]]]]=1-Total[x[[p]]];
];
equations={};
restrictions ={};
For[p=1,p<=pop,p++,
For[j=1,j<=dim[[p]],j++,
If[indicator[[p,j]]==1 && j!=subdimlast[[p]],
equations=Append[equations,dyn[x,F][[p,j]]];
restrictions = Append[restrictions,{x[[p,j]],x0[[p,j]]}];
];
]];
Off[FindRoot::jsing];
root=FindRoot[equations,restrictions, MaxIterations ->1000,AccuracyGoal->accuracy,PrecisionGoal->accuracy];
On[FindRoot::jsing];
root=Table[root[[j,2]],{j,Dimensions[root][[1]]}];root=Table[vectorReplace[x0[[p]],
Append[
root[[Sum[subdim[[j]]-1,{j,p-1}]+1;;Sum[subdim[[j]]-1,{j,p}]]]
,0],indicator[[p]]],{p,pop}];
Do[root[[p,subdimlast[[p]]]]=1-Total[root[[p]]];
,{p,pop}];
,
Message[subFindRoot::inconsistency];
If[dyn[x0,F]==Table[ConstantArray[0,dim[[p]]],{p,pop}],root={x0},root={}];
];
root
];

restPointsContinuous::usage="Calculates rest points of continuous dynamics, returning a list with them."  ;
restPointsContinuous::notes = "The arguments are the dynamic, the payoff function, a list of points to use as starting points for the calculation of the roots of the differential equation, a tolerance level for speed and two optional parameters. The first optional parameter is a proximity notion to eliminate duplicate rest points that are close to each other, and the second one gives a level of digits accuracy. There is an extra optional parameter that is used in low level programming that allows to search for rest points in a subsimplex by parsing an indicator vector that says through which dimensions should the rest point be looked for, ignoring the remaining ones.";
restPointsContinuous::inconsistency = "some of the roots of the differential equation have been adjusted beyond the tolerance level";
restPointsContinuous[dyn_,F_,startingValues_,tolLevel_:.00001,delta_:.001,accuracy_:8,sup_:{},dim_:dim]:= Module[{j,n,pop,support,grid,roots},
pop=Dimensions[dim][[1]]; 
If[sup=={},
support=Table[ConstantArray[1,dim[[p]]],{p,pop}],
support=sup];
n=Dimensions[startingValues][[1]];
(* Finds rest points by grid search on the subsimplex determined by support *) 
roots=Table[subFindRoot[dyn,F,startingValues[[j]],support,accuracy],{j,n}];
(* Eliminate rest points due to singularities in Jacobian*)
roots=DeleteCases[Table[If[And[startingValues[[j]]==roots[[j]],stateNorm[F[roots[[j]]]]>delta],"discard",roots[[j]]],{j,n}],"discard"];
(* Round to Accuracy *)
roots=Round[roots,10^(-accuracy)];
(* Eliminates double countings *)
roots=eliminateDuplicates[delta,roots,stateNorm];
(* Print[roots]; *)
n=Dimensions[roots][[1]];
(* clever code by emin, similar to Eliminate duplicates
reduce[points_]:=DeleteCases[Table[
If[(Abs[(x1 /. points[[1,1]])- (x1/. points[[j,1]])]>delta || Abs[(x2 /.  points[[1,2]])- (x2 /. points[[j,2]] )]>delta),points[[j]],0]
,{j,2,Dimensions[points][[1]]}],0];
roots=FixedPointList[reduce,roots];
roots=DeleteCases[roots,{}];
*)
(* Check points are in the simplex, eliminate them if not *)
roots = DeleteCases[Table[If[Total[Table[If[simplexQ[roots[[j,p]],accuracy],1,0],{p,pop}]]==pop,roots[[j]],"discard"],{j,n}],"discard"];
n=Dimensions[roots][[1]];
(* Check points to have "low speed", eliminate them if not. This makes no sense, since the points are roots of the equation, a double check after the roundoff I guess; look at the explanation in the startingValues definition above.*)
roots=DeleteCases[Table[If[Max[Abs[dyn[roots[[j]],F]]] > tolLevel, "discard",roots[[j]]],{j,n}],"discard"];
If[n>Dimensions[roots][[1]],Message[restPointsContinuous::inconsistency]];
roots
] ;

Calculation of restricted equilibria, Nash equilibria and Quasistrict equilibria

restrictedEquilibria::usage="is a procedure that finds the Restricted Nash equilibria of a game for payoff vector field F. It returns a list with the equilibria." ;
restrictedEquilibria::notes="The algorithm finds Restricted Equilibria by analyzing each possible subsimplex. For each subsimplex, it finds payoff indiference points or rest points of the Replicator Dynamic using root finding numerical methods. The arguments are the payoff fucntion and a number of optional computation parameters. The first parameter gives the method. The second optional parameter gives the notion of proximity used to discard repeated equilibria in the rest point method, the third one the level of decimal accuracy, the fourth one the size of the grid used to search for the equilibria and the fifth one a size of perturbation to build the grid.";
restrictedEquilibria[F_,method_:"indif",delta_:.01,accuracy_:8,gridsize_:10 ,eps_:.01,dim_:dim]:=Module[ {j,p,pop,subdim,subdimlast,x,fx, RNE,supports,restPoints,grid},
(* dim=Dimensions[F][[1]]; *)
pop=Dimensions[dim][[1]];
supports=Tuples[Table[Drop[Tuples[{0,1},dim[[p]]],1],{p,pop}]];
RNE={};
For[j=1,j<=Dimensions[supports][[1]],j++,
subdim=Table[Total[supports[[j,p]]],{p,pop}];
subdimlast=Table[Max[Table[If[supports[[j,p,k]]==1,k,0],{k,dim[[p]]}]],{p,pop}];
If[Total[subdim]>pop,
If[method=="indif",
x=emptyVector["x",Total[subdim]];
x=Table[x[[Sum[subdim[[k]],{k,p-1}]+1;;Sum[subdim[[k]],{k,p}]]],{p,pop}];
fx=F[Table[
vectorReplace[supports[[j,p]],x[[p]],supports[[j,p]]]
,{p,pop}]];
fx=Table[DeleteCases[Table[If[supports[[j,p,k]]==1,fx[[p,k]],"discard"],{k,dim[[p]]}],"discard"] ,{p,pop}];
restPoints=
Quiet[Solve[
Flatten[Table[Thread[
Append[
Drop[fx[[p]],-1]-ConstantArray[fx[[p,subdim[[p]]]],subdim[[p]]-1],
Total[x[[p]]]-1
]
==
ConstantArray[0,subdim[[p]]]
],{p,pop}],1]
,
Flatten[x,1]
]];
restPoints=Table[If[Dimensions[restPoints[[k]]][[1]]<Total[subdim],{},Table[restPoints[[k,l,2]],{l,Total[subdim]}]],{k,Dimensions[restPoints][[1]]}];
If[restPoints=={{}},
(*NON TRIVIAL COMPONENT DETECTED, PRINT WARNING MESSAGE*)
,
restPoints=Table[Table[restPoints[[k]][[Sum[subdim[[k]],{k,p-1}]+1;;Sum[subdim[[k]],{k,p}]]],{p,pop}],{k,Dimensions[restPoints][[1]]}];
restPoints=DeleteCases[Table[If[Total[Table[If[simplexQ[restPoints[[k,p]],accuracy],1,0],{p,pop}]]==pop,
Table[vectorReplace[supports[[j,p]],restPoints[[k,p]],supports[[j,p]]],{p,pop}]
,
"discard"
],{k,Dimensions[restPoints][[1]]}],"discard"];
];
,
grid=Tuples[Table[discSimplex[gridsize,subdim[[p]],True,eps],{p,pop}]];
grid = Table[Table[vectorReplace[ConstantArray[0,dim[[p]]],grid[[k,p]],supports[[j,p]]],{p,pop}],{k,Dimensions[grid][[1]]}];
(* Find the rest points for the Replicator dynamic without tolerance check and for the specific subsimplex *)
restPoints=restPointsContinuous[Replicator,F,grid,Infinity,delta,accuracy,supports[[j]]];
restPoints=DeleteCases[Table[
If[Total[dim]-Total[subdim]==Total[Table[Total[Table[If[restPoints[[k,p,l]]==0,1,0],{l,dim[[p]]}]],{p,pop}]],restPoints[[k]],"discard"]
,{k,Dimensions[restPoints][[1]]}],"discard"];
];
,
restPoints={supports[[j]]};
];
RNE=Append[RNE,restPoints];
];
RNE
];

nashEquilibria::usage="is a procedure that finds Nash equilibria of a game for payoff vector field F. It returns a list with the equilibria." ;
nashEquilibria::notes="The algorithm finds Nash Equilibria analitically by checking inequalities. Alternatively, it first finds the Restricted equilibria and then performs a test to see if the candidate equilibria respect the variational equation that defines a NE. The arguments are the payoff fucntion and a number of optional computation parameters. The first optional parameter chooses the method, while the second allows for a variant that calculates quasistrict equilibria instead of plain NE. The third option gives the notion of proximity used to discard repeated equilibria, the fourth one the level of decimal accuracy and the fifth one the size of the grid used to search for the equilibria. The next three parameters are used in case the sampling method is selected, making it active with the first of the three options, and specifying the size of the sample to check the variational equation and the radius of the sample.";
nashEquilibria[F_,method_:"indif",variant_:"standard",delta_:.001,accuracy_:8,gridsize_:10 ,n_:100,eps_:.01,dim_:dim]:=Module[ {j,k,l,pop, NE,RNE,supports,sample,Fx,flag},
pop=Dimensions[dim][[1]];
supports=Tuples[Table[Drop[Tuples[{0,1},dim[[p]]],1],{p,pop}]];
NE={};
RNE=restrictedEquilibria[F,method,delta,accuracy,gridsize];
For[j=1,j<=Dimensions[supports][[1]],j++,
(* Following lines eliminate the rest points that are not NE.*)
If[j<Dimensions[supports][[1]],
For[k=1,k<=Dimensions[RNE[[j]]][[1]],k++,
If[RNE[[j,k]]!={},
Fx=F[RNE[[j,k]]];
flag=True;
If[method=="indif",
(* this method will check against the pure strategies that are left out of the support, it's more efficient *)
flag=Apply[And,Table[variationalIneqCheck[RNE[[j,k,p]],supports[[j,p]],Fx[[p]],variant],{p,pop}]];
,
(* this method will check against a sample of points around the subsimplex rest point, it's less efficient *)
sample=Table[sampleAround[RNE[[j,k,p]],eps,n],{p,pop}];
sample=Table[Table[sample[[p,l]],{p,pop}],{l,n}];
For[l=1,l<= n,l++,
If[variant=="quasistrict",
flag=Apply[And,Table[(RNE[[j,k,p]]-sample[[l,p]]).Fx[[p]]>0,{p,pop}]];
,
flag=Apply[And,Table[(RNE[[j,k,p]]-sample[[l,p]]).Fx[[p]]>=0,{p,pop}]];
];
(* Print[i,RNE[[j,k]]]; *)
If[Not[flag],Break[]];
];
];
If[flag,NE=Append[NE,RNE[[j,k]]]];
];
];
,
If[RNE[[j]]!={{}},NE=Join[NE,RNE[[j]]]];
];
];
NE=eliminateDuplicates[delta,NE,stateNorm];
NE
];

Specification of rest points

restPoints::usage="Calculation of rest points of dynamics" ;
restPoints::notes = "The arguments are the dynamic, the payoff function, a list of points to use as starting points for the calculation of the roots of the differential equation, a tolerance level for speed and two optional parameters. The first optional parameter is a proximity notion to eliminate duplicate rest points that are close to each other, and the second one gives the number of digits of accuracy.";
restPoints[dyn_,F_,startingValues_,tolLevel_,delta_:.001,accuracy_:8]:=Which[
RPCharacterization[dyn]=== "Automatic",restPointsContinuous[dyn,F,startingValues,tolLevel,delta,accuracy],
RPCharacterization[dyn]=== "Nash",
nashEquilibria[F,"indif","standard",delta,accuracy],
RPCharacterization[dyn]==="RestrictedEquilibria",
restPointsContinuous[Replicator,F, startingValues,tolLevel,delta,accuracy],
RPCharacterization[dyn]==="RestrictedEquilibriaIndif",
restrictedEquilibria[F]
];

Stability

Stability Test Subroutine for Smooth Dynamics

smoothStabilityTest::usage="tests stability of a rest point for Smooth Dynamics through Linearization. The possible values that the function returns are 1 for a stable state, -1 for a source and 0 otherwise (saddle,cycle). This algorithm ignores 0 eigenvalues." ;
smoothStabilityTest::notes="Its arguments are the point, the name of the dynamic, the payoff vector field and optionally the method to use which is by default analytic and could be numeric, with a second optional parameter for regulating the precission to evaluate the sign of the eigenvalues which is useful under the numeric derivation method." ;
Needs["NumericalCalculus`"];
smoothStabilityValues[x_,dyn_,F_,method_:"Analytic",accuracy_:8]:=Module[{j,p,dim,pop,s,xv,xvp,DV,DVp,es,result},
pop=Dimensions[x][[1]];
dim=Table[Dimensions[x[[p]]][[1]],{p,pop}];
If[method=="Analytic",
xv=emptyVector["x",Total[dim]];
xvp=Table[xv[[Sum[dim[[k]],{k,p-1}]+1;;Sum[dim[[k]],{k,p}]]],{p,pop}];
DV=D[Flatten[dyn[xvp,F],1],{xv}]/.Thread[xv->Flatten[x,1]];
,
DV=Transpose[Flatten[Table[Table[
ND[Flatten[dyn[Table[Table[If[And[pt==p,jt==j],s,x[[pt,jt]]],{jt,dim[[p]]}],{pt,pop}],F],1],s,x[[p,j]]]
,{j,dim[[p]]}],{p,pop}],1]];
];
es=Eigensystem[DV.projectionMatrix[dim]];
(* notice the eigenvalues could be complex; Round[] works anyway. *)
(* not checking independence of the eigenvectors or that we have pop of them *)
DeleteCases[Table[If[Round[es[[1,j]],10^(-accuracy)]==0 &&Norm[es[[2,j]]-Flatten[Table[ConstantArray[es[[2,j,If[p==1,0,Product[dim[[k]],{k,p-1}]]+1]],dim[[p]]],{p,pop}]]]<10^(-accuracy)
,"discard",{es[[1,j]],es[[2,j]]}],{j,Dimensions[es[[1]]][[1]]}],"discard"]
];
smoothStabilityTest[x_,dyn_,F_,method_:"Analytic",accuracy_:8]:=Module[{values,es},
(*values=Flatten[smoothStabilityValues[x,dyn,F,method,accuracy],1];*)
values=smoothStabilityValues[x,dyn,F,method,accuracy];
es=Table[Re[values[[j,1]]],{j,Dimensions[values][[1]]}];
Which[
Round[Max[es],10^(-accuracy)]<=0,1,
Round[Min[es],10^(-accuracy)]>=0,-1,
True,0
]
];
(* Replicator, 1-2-3 coord
smoothStabilityTest[{1/3,1/3,1/3},dyn,F,"Numeric",10]
smoothStabilityTest[{0,0.6,0.4},dyn,F,"Analytic"]
smoothStabilityTest[{1,0,0},dyn,F]
*)

Numeric Stability Test Subroutine for non Smooth Dynamics

numericStabilityTest::usage="tests the stability of a rest point using a numeric approach which does not require dynamic differentiability, although it does require continuity. The test is imprecise and heuristic in nature. It returns a value of 1 if the point is stable, 0 if it's a cycle, -1 for a source and -2 if it's a saddle." ;
numericStabilityTest::notes="Its arguments are the point to be evaluated, the dynamic, the payoff function and a number of optional parameters. The first of these optional parameters gives the number of initial points for which tests will be conducted. The second one is the final time which will be used to see whether a solution trajectory is outside a neighborhood of the rest point. This neighborhood is defined by a ball of radious established by the following option, call it delta. The option that follows, say eps, defines another neighborhood of the rest point from which an initial condition for the DE is drawn. Notice that the choice of delta has to be smaller than that of eps, usually delta should be a tenth of epsilon and no more than a quarter of it." ;

numericStabilityTest[x_,dyn_,F_, n_:10,T_:10, delta_:.01, eps_:.05]:=Module[{j,k,dim,pop,sample,ns,trajectories,classifiedTrajectories},
pop=Dimensions[x][[1]];
dim=Table[Dimensions[x[[p]]][[1]],{p,pop}];
sample=Table[sampleAround[x[[p]],eps,n],{p,pop}];
sample=Table[Table[sample[[p,j]],{p,pop}],{j,n}];
(* the following code eliminated too many initial conditions that are close to each other *)
sample=DeleteCases[Table[If[stateNorm[x-sample[[j]]]<eps/3,"discard",sample[[j]]],{j,n}],"discard"];
ns=Dimensions[sample][[1]];
trajectories = Table[difEqAdvance[sample[[j]],dyn,F,T,AOCharacterization[dyn]],{j,ns}];(*- Checking Stability using delta neighborhood -*)
classifiedTrajectories =Table[
Which[
stateNorm[x-trajectories[[j]]]<delta,
1,
stateNorm[x-trajectories[[j]]]>eps+delta,
(*stateNorm[x-trajectories[[j]]]>stateNorm[x-sample[[j]]]+delta, *)
-1,
True,
0
(*
Abs[stateNorm[x-sample[[j]]]-stateNorm[x-trajectories[[j]]]]<delta,
0,
True,
-1
*)
(*
notice that the time elapsed has to be large if we want to capture a cycle that is spiriling outwards (heteroclinic cycle), so that the traj is suf further away from the rest point than the initial trajectory. However, the values returned will be either 0 or -1 but the point won't be confused with a stable rest point.
*)
]
,{j,ns}];
Which[
(* Min[classifiedTrajectories]==1, *)
And[Min[classifiedTrajectories]>=0,Max[classifiedTrajectories]==1],
1,
(* Max[classifiedTrajectories]==-1, *)
And[Max[classifiedTrajectories]<=0,Min[classifiedTrajectories]==-1],
-1,
And[Min[classifiedTrajectories]==0 ,Max[classifiedTrajectories]==0],
0,
And[Min[classifiedTrajectories]==-1,Max[classifiedTrajectories]==1],
-2,
True,
-3
]
];

Regular ESS

regularESSTest::usage="tests whether a supplied quasistrict equilibrium is a RTESS. The possible values that the function returns are 1 for a RTESS and 0 if not. In case all the eigenvalues have zero real part, the routine returns -1 implying it was unable to determine if it is a RTESS or not." ;
regularESSTest::notes="Its arguments are the point, the payoff vector field and optionally the method to use which is by default analytic and could be numeric, with a second optional parameter for regulating the precission to evaluate the sign of the eigenvalues which is useful under the numeric derivation method." ;
Needs["NumericalCalculus`"];
regularESSTest[x_,F_,method_:"Analytic",accuracy_:8]:=Module[{j,k,dim,pop,subdim,s,xv,xvp,DF,project,matrix,support,es}, 
pop=Dimensions[x][[1]];
dim=Table[Dimensions[x[[p]]][[1]],{p,pop}];
If[method=="Analytic",
xv=emptyVector["x",Total[dim]];
xvp=Table[xv[[Sum[dim[[k]],{k,p-1}]+1;;Sum[dim[[k]],{k,p}]]],{p,pop}];
DF=D[Flatten[F[xvp],1],{xv}]/.Thread[xv->Flatten[x,1]];
,
DF=Transpose[Flatten[Table[Table[
ND[Flatten[F[Table[Table[If[And[pt==p,jt==j],s,x[[pt,jt]]],{jt,dim[[p]]}],{pt,pop}]],1],s,x[[p,j]]]
,{j,dim[[p]]}],{p,pop}],1]];
];
matrix=DF+Transpose[DF];
support=Table[ConstantArray[1,dim[[p]]]-checkBoundary[x[[p]],10^(-accuracy)],{p,pop}];
matrix=subMatrix[matrix,Flatten[support,1]];
subdim=Table[Total[support[[p]]],{p,pop}];
project=projectionMatrix[subdim];
es=Eigensystem[project.matrix.project];
es=DeleteCases[Table[If[Round[es[[1,j]],10^(-accuracy)]==0 && es[[2,j]]==Flatten[Table[ConstantArray[es[[2,j,If[p==1,0,Product[subdim[[k]],{k,p-1}]]+1]],subdim[[p]]],{p,pop}]],"discard",{es[[1,j]],es[[2,j]]}],{j,Dimensions[es[[1]]][[1]]}],"discard"];
es=Table[Re[es[[j,1]]],{j,Dimensions[es][[1]]}];
If[Round[Max[es],10^(-accuracy)]<0,1,If[And[Max[es]==0,Min[es]==0],-1,0]]
];
(* Replicator, 1-2-3 coord
regularESSTest[{{1/3,1/3,1/3}},Replicator,Fmp,"Numeric",10]
regularESSTest[{{0,0.6,0.4}},Replicator,Fmp,"Analytic"]
regularESSTest[{{1,0,0}},Replicator,Fmp]
*)

regularESSs::usage="Finds the set of quasistrict equilibria which are RTESS points for a given Payoff Function. It retunrs two lists with the QSEa.";
regularESSs::notes="Its arguments are the rest points or NE to be evaluated, the payoff function and a number of optional parameters. The first of these optional parameters determines whether a numeric or a differential approach will be used to calculate the derivative function of the payoff function. The later is usually preferable. The last parameter indicates the number of digits of precission." ;
regularESSs[points_,F_,method_:"Analytic",accuracy_:8]:=Module[{j,stable,unstable},
stable={};
unstable={};
For[j=1,j<=Dimensions[points][[1]],j++,
If[quasiStrictEqQ[points[[j]],F],
If[regularESSTest[points[[j]],F,method,accuracy]==1,
stable=Append[stable,points[[j]]],
unstable=Append[unstable,points[[j]]]
];
];
];
{stable,unstable}
];

Stability Sets Routine

stabilitySets::usage="tests the stability of a set of rest points. It takes a list of them and returns two lists, the first one with the stable points and the second one with the unstable ones.";
stabilitySets::notes="Its arguments are the points to be evaluated, the dynamic, the payoff function and a number of optional parameters. The first of these optional parameters determines whether a numeric or a differential approach will be used. The later is preferable when possible since it does not rely on large sample measure properties and it is the default. The rest of the arguments parsed are usueful only under the numeric approach. The first argument after parsing the word numeric in quotes, indicates the number of initial points for which tests will be conducted. The second one is the final time which will be used to see whether a solution trajectory is outside a neighborhood of the rest point. This neighborhood is defined by a ball of radious established by the following option, call it delta. The option that follows, say eps, defines another neighborhood of the rest point from which an initial condition for the DE is drawn. Notice that the choice of delta has to be smaller than that of eps, usually delta should be a tenth of epsilon and no more than a quarter of it." ;

stabilitySets[points_,dyn_,F_, method_:"ESSNumeric", n_:10,T_:10, delta_:.01, eps_:.05]:=Module[{j,stablepoints,unstablepoints,flag},
stablepoints={};
unstablepoints={};
Which[
method=="Smooth",
For[j=1,j<=Dimensions[points][[1]],j++,
(* If[smoothStabilityTest[points[[j]],dyn,F]==1, *)
If[smoothStabilityTest[points[[j]],Replicator,F]==1,
stablepoints=Append[stablepoints,points[[j]]],
unstablepoints=Append[unstablepoints,points[[j]]]
];
];
,method=="Numeric",
For[j=1,j<=Dimensions[points][[1]],j++,
(* If[numericStabilityTest[points[[j]],dyn,F,n,T,delta,eps]==1, *)
If[numericStabilityTest[points[[j]],dyn,F,n,T,delta,eps]==1,
stablepoints=Append[stablepoints,points[[j]]],
unstablepoints=Append[unstablepoints,points[[j]]]
];
];
,method=="ESS",
For[j=1,j<=Dimensions[points][[1]],j++,
If[regularESSTest[points[[j]],F]==1,
stablepoints=Append[stablepoints,points[[j]]],
unstablepoints=Append[unstablepoints,points[[j]]]
];
];
,method=="ESSNumeric",
For[j=1,j<=Dimensions[points][[1]],j++,
If[quasiStrictEqQ[points[[j]],F],
flag=regularESSTest[points[[j]],F];
Which[
flag==1,
stablepoints=Append[stablepoints,points[[j]]],
(* the proofs in chapter 8 go in one direction only; if it's not an ESS we don't know if it's unstable
flag==0,
unstablepoints=Join[unstablepoints,{points[[j]]}],
*)
Or[flag==0,flag==-1],
If[numericStabilityTest[points[[j]],dyn,F,n,T,delta,eps]==1,
stablepoints=Append[stablepoints,points[[j]]],
unstablepoints=Append[unstablepoints,points[[j]]]
];
];
,
If[numericStabilityTest[points[[j]],dyn,F,n,T,delta,eps]==1,
stablepoints=Append[stablepoints,points[[j]]],
unstablepoints=Append[unstablepoints,points[[j]]]
];
];
];
];
{stablepoints,unstablepoints}
];

Basins

basinsTrajectories[unstableRestPoints_,dyn_,F_,altopt_:False,T_:20,n_:100,delta_:.05,eps_:.01,accuracy_:8]:=Module[{j,k,l,pop,sample,trajectories,solutions,nearTrajectories,selectedTrajectories},
pop=Dimensions[unstableRestPoints[[1]]][[1]];
selectedTrajectories={};
For[j=1,j<=Dimensions[unstableRestPoints][[1]],++j,
sample=Table[sampleAround[unstableRestPoints[[j,p]],eps,n],{p,pop}];
sample=Table[Table[sample[[p,k]],{p,pop}],{k,n}];
solutions=Table[difEqSol[sample[[k]],dyn,F,0,T,"final",altopt,accuracy],{k,n}];
trajectories=X[0] /. solutions;
nearTrajectories={};
k=1;
While[k<=Dimensions[unstableRestPoints][[1]],
If[k==j ,++k;Continue[]];
l=1;
While[l<=n,
(* Print[j,trajectories[[i,1]]]; *)
If[stateNorm[trajectories[[l,1]]-unstableRestPoints[[k]]]<delta,
nearTrajectories=Append[nearTrajectories,{sample[[l]],0,trajectories[[l,1]]}];
Break[];
];
++l
];
++k;
];
selectedTrajectories=Join[selectedTrajectories,nearTrajectories];
];
selectedTrajectories
];

basinsBoundaryTrajectories[unstableRestPoints_,dyn_,F_,altopt_:False,T_:20,n_:100,delta_:.05,eps_:.01,accuracy_:8]:=Module[{t0,j,k,l,pop,sample,trajectories,nearTrajectories,selectedTrajectories,solutions,point,possibleboundaries,boundaryvector1,boundaryvector2,doublebreak,nextpoint,repeated,ntrajectories},
selectedTrajectories={};
pop=Dimensions[unstableRestPoints[[1]]][[1]];
For[j=1,j<=Dimensions[unstableRestPoints][[1]],++j,
(* when the rest point is in the boundary, there's a trajectory that starts on the edge of the simplex and lands on it, find it*)
sample=Table[sampleAround[unstableRestPoints[[j,p]],eps,n],{p,pop}];
sample=Table[Table[sample[[p,k]],{p,pop}],{k,n}];
solutions=Table[difEqSol[sample[[k]],dyn,F,0,T,"final",altopt,accuracy],{k,n}];
trajectories=X[0] /. solutions;
(*Print[trajectories];*)
nearTrajectories={};
(* I was planning to use this to chck for connections to every possible boundary but then dropped the idea; the code could be useful in the future though; 
possibleboundaries =Take[Tuples[{0,1},Dimensions[unstableRestPoints[[j]]][[1]]],{2,-2}];
*)
boundaryvector1=Table[checkBoundary[unstableRestPoints[[j,p]],delta],{p,pop}];
(*Print[boundaryvector1];*)
l=1;
ntrajectories=0;
While[l<=n,
t0=T;
repeated=0;
While[t0>0,
point=X[t0] /. solutions[[l]];
point=point[[1]];
(* Print[point]; *)
boundaryvector2=Table[checkBoundary[point[[p]],delta],{p,pop}];
(* Print[{boundaryvector1,boundaryvector2}]; *)
For[k=1,k<=Dimensions[nearTrajectories][[1]],k++,
If[ stateNorm[nearTrajectories[[k,3]]-point]<delta,
repeated =1;
(* Print["k",k]; *)
Goto[nextpoint];
];
];
If[Total[boundaryvector2]>0 && Max[boundaryvector2-boundaryvector1]==1 && repeated==0,
nearTrajectories=Append[nearTrajectories,{sample[[l]],t0,point}];
(* Print[nearTrajectories]; *)
(* Print["t",t0]; *)
ntrajectories++;
Break[];
];
t0= t0-delta ;
];
(* Print["i",i]; *)
(* the following is usuful only in planar system where we know there can be at most 2 trajectories converging to a saddle*)
If[ntrajectories>=2,Break[]];
Label[nextpoint];
l++;
];
selectedTrajectories=Join[selectedTrajectories,nearTrajectories];
(*
Print[unstableRestPoints[[j]]];
Print[selectedTrajectories];
*)
];
selectedTrajectories
];

Lyapunov Exponents

MLE calculation

maximumLyapunovExponent::usage="calculates the maximal rate of expansion of the system in the atractor, as measured by the Maximum Lyapunov Exponent" ;
maximumLyapunovExponent::notes="Its arguments are the initial condition, the name of the dynamic, the payoff vector field, time elapsed to let the initial condition land in the attractor, the number of iterations, the number of initially ignored iterations, the size of the step at each iteration, the initial separation of the two trajectories, the distance criterion to assume the main trajectory is not moving and has converged, the accuracy of the ODE solver at each substep, the max steps of the ODE solver at each iteration. Notice that the effective precision of the ODE solver ep is obtained by substracting the order of magnitude of the maximum number of steps from the accuracy of the ODE solver at each substep. ep has to be larger than the initial separation of the two trajectories, and if the two trajectories at some point become closer than 10^(-ep) then the method fails." ;
maximumLyapunovExponent[x_,dyn_,F_,t0_,nsteps_:10000,nskip_:10000,step_:.01,separation_:.00000001,convergence_:.001,accuracy_:14,maxsteps_:100000]:=Module[{pop,main,neighbor,d,mains,ds},
pop=Dimensions[x][[1]];
main=difEqAdvance[x,dyn,F,t0,AOCharacterization[dyn],accuracy,maxsteps,"ExplicitRungeKutta"];
neighbor=Table[sampleAround[main[[p]],separation,1,"preciselyat"][[1]],{p,pop}];
mains={};
ds={};
Do[
main=difEqAdvance[main,dyn,F,step,AOCharacterization[dyn],accuracy,maxsteps,"ExplicitRungeKutta"];
neighbor=difEqAdvance[neighbor,dyn,F,step,AOCharacterization[dyn],accuracy,maxsteps,"ExplicitRungeKutta"];
d=stateNorm[main-neighbor];
neighbor=main+(neighbor-main)*(separation/d);
mains=Append[mains,main];
ds=Append[ds,d];
,
{nskip+nsteps}
];
{
If[Max[Table[stateNorm[mains[[j]]-mains[[nskip+1]]],{j,nskip+2,nskip+nsteps}]]<convergence,"convergent","non convergent"]
,
Total[Log[Drop[ds,nskip]/separation]]/nsteps/step
}
];

(* Replicator 9 9 9, 4 4 4, 1 1 1
maximumLyapunovExponent[{{.9,.05,.05}},dyn,Fmp,20,1000,100,.01,.001]
*)

FixPolygons

fixPolygons::usage =
"fixPolygons[ <graphic> , <options> ]:
  Joins all contiguous polygons into a single shape. 
  This can dramatically improve the image quality and 
  decrease the file size of a ContourPlot or RegionPlot.

Possible options (with defaults) are:
  earlyCleanupThreshold -> 1
";

earlyCleanupThreshold::usage = 
"At every merging iteration, all the vertices
shared by less than earlyCleanupThreshold polygons are 
deleted. Higher values means that cleanup is faster, but some 
edges might not be removed. The default of 1 leads calls for 
complete cleanup. Sometimes an overall speedup can be 
acheived by applying FixPloygons several times, reducing 
progressively the value of this option from approx. 6 to 1.";
Options[fixPolygons]={
  earlyCleanupThreshold -> 1
};
Begin["`Private`"];
(* 
fixPolygons function
Delve into GraphicsComplex and grab polygons to combine:
*)
fixPolygons[graf_,OptionsPattern[fixPolygons]] := 
  graf //. {
    Polygon[{a__},c___],
    Polygon[{b__},c___],d___} ->
      {Polygon[{a,b},c],d} /. 
        {Polygon[{ab__},c___],d___} :>
          {Polygon[cleanup[{ab},
             OptionValue[earlyCleanupThreshold]],c],d};
(* Cleanup polygons *)
Module[{vert,uvert,cvert,nvert,i,plist},
cleanup[polylist_,ect_]:=(
  plist = polylist;

(* Makes up a list of the occurrences of vertices in polygon list *)
  vert = Sort[Flatten[plist]]; 
  uvert = Union[vert]; nvert=Length@uvert;
(* Counts the occurrences of vertices in polygon list *)
  cvert = Length /@ Split@vert;
(* Preliminary cleanup based on earlyCleanupThreshold option *)
  vert = DeleteCases[
    Transpose@{uvert,cvert},
	{_,a_}/;a<=ect];
(* Cleans up the vertex list *)
  Do[
    If[Length@vert==0,Break[]];
    {plist,vert} = refine[{plist,vert},vert[[1,1]],ect];
  ,{i,1,nvert}];
(* Output and finish up, with a final cleanup from residual creaks *)
  plist //. {{c_,a_,b___,a_} :> {a,b},
	{c___,a_,b_,a_,d___} :> {c,a,d}}
);
];
(* Clean up vertices of polygons *)
Module[{ri,rj,lvp,nvl,ncl,ppoly,vpoly},
refine[{cl_,vl_},iv_,ect_]:=(
(* Gets the polygons having the shared vertex *)
  ppoly = Position[cl,{___,iv,___}]; 
  nvl = vl;
  vpoly = cl[[Flatten[ppoly]]];
  lvp = Length@vpoly;
(* Puts all the polygons in a position where the shared vertex is at position 1 *)
  vpoly = Table[
    RotateLeft[vpoly[[ri]],Position[vpoly[[ri]],iv][[1,1]]-1]
  ,{ri,1,lvp}];
(* Finds the neighbours of the polygon to be deleted *)
  ri = 1;
  rj = 2;
  While[ri <= lvp && rj <= lvp ,
    If[vpoly[[ri,2]] == vpoly[[rj,-1]],
      nvl = nvl /. {vpoly[[ri,2]],a_} :> {vpoly[[ri,2]],a-1};
      vpoly = Prepend[
        Delete[vpoly,{{ri},{rj}}],
		Join[vpoly[[rj]],Drop[vpoly[[ri]],2]]];
	  lvp--;
      ri = rj = 1
    ];
    rj++;
    If[rj==ri, rj++];
    If[rj>lvp, ri++; rj=1];
  ];
(* Cleans up polygons with "creaks" *)
  vpoly = Table[
    vpoly[[ri]] //. {{c_,a_,b___,a_} :> 
     (nvl = (nvl /. {a,e_} :> {a,e-1});
      {a,b}),
	{c___,a_,b_,a_,d___} :> 
     (nvl = (nvl /. {a,e_} :> {a,e-1});
     {c,a,d})}
   ,{ri,1,lvp}];
  ncl = Join[Delete[cl,ppoly],vpoly];
  nvl = DeleteCases[nvl,{iv,_}|{_,a_}/;a<=ect];
(* Output and finish: *)
  {ncl,nvl}
);
];
End[];

Maintenance

Reformat Open Group Cells

FrontEndExecute[FrontEndToken["SelectAll"]];
FrontEndExecute[FrontEndToken["SelectionCloseAllGroups"]];
NotebookFind[SelectedNotebook[],"opengroup", Next, CellTags, AutoScroll -> False];
While[If[Depth[Options[NotebookSelection[],CellTags]]==3,True,Not[ToString[Options[NotebookSelection[],CellTags][[1,2,1]]]=="finalcell"]],
FrontEndExecute[FrontEndToken["OpenCloseGroup"]];
NotebookFind[SelectedNotebook[],"opengroup", Next, CellTags, AutoScroll -> False];
];
NotebookFind[SelectedNotebook[],"opengroup", Next, CellTags, AutoScroll -> True];

Hide Cells

NotebookFind[SelectedNotebook[],"notopen", All, CellTags, AutoScroll -> False];
SetOptions[NotebookSelection[SelectedNotebook[]],CellOpen->False];
NotebookFind[SelectedNotebook[],"opengroup", Next, CellTags, AutoScroll -> True];
FrontEndExecute[FrontEndToken["SelectAll"]];
NotebookFind[SelectedNotebook[],"opengroup", Next, CellTags, AutoScroll -> True];

Buttons

FrontEndExecute[FrontEndToken["SelectAll"]];NotebookFind[SelectedNotebook[],"buttoncell", Next, CellTags, AutoScroll -> False];
While[If[Depth[Options[NotebookSelection[],CellTags]]==3,True,Not[ToString[Options[NotebookSelection[],CellTags][[1,2,1]]]=="finalcell"]],
If[ToString[Options[NotebookSelection[],GeneratedCell][[1,2]]]=="False",
SetOptions[NotebookSelection[SelectedNotebook[]],CellOpen->False,Evaluatable->False];
FrontEndExecute[FrontEndToken["OpenCloseGroup"]];
,
SetOptions[NotebookSelection[SelectedNotebook[]],CellOpen->True];
FrontEndExecute[FrontEndToken["OpenCloseGroup"]];
];
NotebookFind[SelectedNotebook[],"buttoncell", Next, CellTags, AutoScroll -> False];
];
NotebookFind[SelectedNotebook[],"opengroup", Next, CellTags, AutoScroll -> True];

Credits

Dynamo
Version 1.2  (8 July 2013)
by William H. Sandholm, Emin Dokumaci and Francisco Franchetti

mail:			Department of Economics
			University of Wisconsin
			1180 Observatory Drive
			Madison WI 53706

e-mail:			whs@ssc.wisc.edu
			edokumaci@wisc.edu
			franchetti@wisc.edu
	
websites:		http://www.ssc.wisc.edu/~whs/
			http://www.ssc.wisc.edu/~edokumac/
			https://mywebspace.wisc.edu/franchetti/web/
			
Dynamo website:	http://www.ssc.wisc.edu/~whs/dynamo
    
Financial support from the National Science Foundation under Grants SES-0092145, SES-0617753, SES-0851580, and SES-1155135 is gratefully acknowledged.