leopardSpots28IV00.htm

file:C:/MAPLE/ leopardsSpot28april00.mws

L E O P A R D S ' S P O T S B Y E N T R O P Y R E D U C T I O N

algorithm and program by W. Richard Stark ,

stark@math.usf.edu,

Department of Mathematics,

University of South Florida,

Tampa, FL 33620 .

0. Initialization, basic constants .

> restart():
Digits := 4:
with(plots):
with(plottools):

Constants N, sYn, avgD, exR, inR, States.

> N := 50;
sYn := floor(1.5*sqrt(N));
avgD := 7:
exR := evalf(80/sqrt(N)):
inR := evalf(90*sqrt(avgD/(Pi*N))):
'[exR,inR]' = [exR,inR];
States := [Y1,Y2,B1,B2,B3,R1,R2,R3]:

1. The tissue simulator.

> run := proc (maxRun,interval) local step;
global A,AA,AFreq,C,Cxy,CO,COxy,Colors,H,N,Q,S,States,sYn,T,
G,GO,inR,exR, edgeCount,eLabels,entHistory,Spectrum,SS,
EDGES,BCs,GCs,PTO,RCs,RYB,YCs,SEA,STY,VWd;

makeS();
initialize();

print(oOo);
print([step, normalizedEdgeEntropy]);

for step from 0 to maxRun do
entropy(step);
if type(step/interval,integer)
then print(entHistory[-1]);
disPlay(false);
if halting()
then step:=maxRun+1; print(HALT); fi;
elif step<interval
then print(entHistory[-1]);
disPlay(false); fi;
nextS(); od;

print(oOo);
print("-------------- final global state ------------------");
disPlay(true);

print(oOo);
print("------------ normalized edge entropy ---------------");
print('maximumEdgeLabelEntropy(inBits)' = maxEntropy);
print(plot(entHistory[3..-1]));
print([initialEntropy, finalEntropy] = [entHistory[3][2],entHistory[-1][2]])
end:

Here are functions which make basic structures.

C is the set of cells, Q is the set of cell-states, S:C-->Q is a global state,

A is a random subset of C of cells which will be active at the next step.

Their activity changes their cell-states and so changes the global state from S to the nextS.

> makeS := proc () local c; global C,N,States,S;
S := array(1..N);
for c in C do S[c]:=randEle(States[1..-3]); od;
evalm(S);
end:
nextS := proc () local c; global A,AA,C,G,N,S,States,sYn,T;
(`AA,C,S,sYn,T must be defined`);
makeA();
for c in A do T[c] := nextQ(c); od;
for c in A do S[c] := T[c]; od;
evalm(S);
end:
makeA := proc () local c,CA,nCA; global A,AA,C,G,N,S,sYn;
(`AA, and sYn must altrady exits`);
CA := subtract(C,[op(AA[-1]),op(AA[-2])]);
nCA := nops(CA);
A :=[];
for c in CA do
if evalb(rand(nCA+1)()<sYn) then
A := [op(A),c]; fi;od;
AA := [op(AA),A];
A;
end:

Initialize sets up certain structures.

Entropy computes the [Shannon] entropy of the set of state pairs defined by the global state.

DisPlay graphs the tissue-structures -- eg figure 1.

Halting determines when the pattern of color assignments has reached a fixed-point .

> initialize := proc () local c; global AA,C,entHistory,maxEntropy,T;
(`this construction preceeds the first makeA() in the first makeS()`):
AA := [[],[]];
T := array(1..N);
entHistory := [];
maxEntropy := evalf(log[2](36));
for c in C do T[c]:=Y4; od;
end:
entropy := proc (step) local entr,i,ns,p;
global edgeCount,entHistory,G,lE,lEmatrix,maxEntropy,nS,S,
Spectrum,States;
labelEO();
makeSpectrum();
entr := 0;
ns := nops(States);
for i from 1 to ns*(1+ns)/2 do
p := Spectrum[i];
if 0<p then
entr := entr - p*log[2](p); fi; od;
entHistory := [op(entHistory),[step,entr/maxEntropy]]:
end:
disPlay := proc (bv) local c;
global C,S,Colors,BRY,BCs,RCs,WCs,YCs,PTO,EDGES,SEA,STY,VWd;
makeBRY(bv):
BCs := map(c->cellC(c,blue), BRY[1]):
RCs := map(c->cellC(c,red), BRY[2]):
YCs := map(c->cellC(c,yellow), BRY[3]):
WCs := map(c->cellC(c,white), BRY[4]):
print(PLOT(PTO,EDGES,op(BCs),op(RCs),op(YCs),op(WCs),SEA,STY,VWd)):
end:
halting := proc () local c,stability;
global C,G,N,S,States,sYn,T;
stability := 0;
for c in C do
if color(S[c],true)=color(nextQ(c),true)
then stability:=stability+1; fi; od;
evalb(stability=N);
end:

The following functions support disPlay.

> makePconsts := proc () local c,col;
global torusSeams,PTS,AST,VEW,EOxy,Coxy,inR,
SEA,EDGES,PTO,STY,VWd,C,S,BRY;
torusSeams := [[0,0],[0,100],[100,100],[100,0],[0,0]]:
PTS := POINTS(op(Cxy),SYMBOL(CIRCLE),COLOUR(RGB,1.0,0.3,0.3)):
AST := AXESSTYLE(BOX):
VEW := VIEW(-inR..100+inR, -inR..100+inR):
SEA := POLYGONS(torusSeams, COLOR(RGB,0.95,0.98,1.0)):
EDGES := POLYGONS(op(EOxy)):
PTO := POINTS(op(COxy), SYMBOL(CIRCLE), COLOR(RGB,1.0,0.4,0.4)):
STY := AXESSTYLE(BOX):
VWd := VIEW(-inR..100+inR, -inR..100+inR):
end:
makeBRY := proc (bv) local c,col,colorPartC; global C,Colors,BRY,S;
BRY := [[],[],[],[]]:
for c in C do
col := color(S[c],bv);
if (col=blue) then BRY[1]:= [op(BRY[1]),c];
elif (col=red) then BRY[2]:= [op(BRY[2]),c];
elif (col=yellow) then BRY[3]:= [op(BRY[3]),c];
elif (col=white) then BRY[4]:= [op(BRY[4]),c]; fi;od:
BRY:
end:
color := proc (q,bv)
if q=Y1 then yellow;
elif (q=Y2 and bv) then yellow;
elif q=B1 then blue;
elif q=B2 then blue;
elif q=B3 then blue;
elif q=R1 then red;
elif q=R2 then red;
elif q=R3 then red;
else white; fi;
end:
cellC := proc (c,col) global Cxy,S;
circle(Cxy[c],5,thickness=7,color=col);
end:

2. Cell-state transition rules.

An active cell's next cell-state is determined by nextQ.

In other words this procedure determine when an edge in the cell-state transition graph will be reaversed.

The subprocedure Y1Y2 describes when an active cell c in state S[c]=Y1 will change to state Y2.

> nextQ := proc (c) local nsp,s; global C,G,S,States;
nsp:=nSpec(c); s:=S[c];
if (Y1Y1(s,nsp) or Y2Y1(s,nsp) or B1Y1(s,nsp) or
B2Y1(s,nsp) or B3Y1(s,nsp) or R1Y1(s,nsp) or
R2Y1(s,nsp) or R3Y1(s,nsp) )
then Y1;
elif (Y1Y2(s,nsp) or Y2Y2(s,nsp) or B1Y2(s,nsp) or
B2Y2(s,nsp) or B3Y2(s,nsp) or R1Y2(s,nsp) or
R2Y2(s,nsp) or R3Y2(s,nsp) )
then Y2;
elif (B1B1(s,nsp) or B3B1(s,nsp) or Y2B1(s,nsp))
then B1;
elif (B1B2(s,nsp) or B2B2(s,nsp))
then B2;
elif (B2B3(s,nsp) or B3B3(s,nsp))
then B3;
elif (Y2R1(s,nsp) or R1R1(s,nsp) or R3R1(s,nsp))
then R1;
elif (R1R2(s,nsp) or R2R2(s,nsp))
then R2;
elif (R2R3(s,nsp) or R3R3(s,nsp))
then R3;
else Y1; fi;
end:

Y1Y1 := proc (s,nsp)
evalb(s=Y1 and (0<bL(nsp) or 1<rD(nsp)));
end:
Y1Y2 := proc (s,nsp)
evalb(s=Y1 and 0=bL(nsp)+rD(nsp));
end:

Y2Y2 := proc (s,nsp)
evalb(s=Y2 and 0=bL(nsp)+rD(nsp) and 0<nsp[1]);
end:
Y2R1 := proc (s,nsp)
evalb(s=Y2 and aLL(nsp)=nsp[2]);
end:
Y2B1 := proc (s,nsp)
evalb(s=Y2 and 1=rD(nsp));
end:
Y2Y1 := proc (s,nsp)
evalb(s=Y2 and ((0=rD(nsp) and 0<bL(nsp)) or 1<rD(nsp)));
end:

B1B1 := proc (s,nsp)
evalb(s=B1 and 1=nsp[6] and 0=nsp[7]+nsp[8]);
end:
B1B2 := proc (s,nsp)
evalb(s=B1 and 1=nsp[7] and 0=nsp[6]+nsp[8]+nsp[5]);
end:
B1Y1 := proc (s,nsp)
evalb(s=B1 and ((0<nsp[4] and 0<nsp[6]) or
(0<nsp[5] and 0<nsp[7]) or
(0<nsp[3] and 0<nsp[8]) or 1<>rD(nsp)));
end:
B1Y2 := proc (s,nsp)
evalb(s=B1 and aLL(nsp)=yW(nsp));
end:

B2B2 := proc (s,nsp)
evalb(s=B2 and 1=nsp[7] and 0=nsp[6]+nsp[8]);
end:
B2B3 := proc (s,nsp)
evalb(s=B2 and 1=nsp[8] and 0=nsp[6]+nsp[7]+nsp[3]);
end:
B2Y1 := proc (s,nsp)
evalb(s=B2 and ((0<nsp[4] and 0<nsp[6]) or
(0<nsp[5] and 0<nsp[7]) or
(0<nsp[3] and 0<nsp[8]) or 1<>rD(nsp)));
end:
B2Y2 := proc (s,nsp)
evalb(s=B2 and aLL(nsp)=yW(nsp));
end:

B3B3 := proc (s,nsp)
evalb(s=B3 and 1=nsp[8] and 0=nsp[6]+nsp[7]);
end:
B3B1 := proc (s,nsp)
evalb(s=B3 and 1=nsp[6] and 0=nsp[7]+nsp[8]+nsp[4]);
end:
B3Y1 := proc (s,nsp)
evalb(s=B3 and ((0<nsp[4] and 0<nsp[6]) or
(0<nsp[5] and 0<nsp[7]) or
(0<nsp[3] and 0<nsp[8]) or 1<>rD(nsp)));
end:
B3Y2 := proc (s,nsp)
evalb(s=B3 and aLL(nsp)=yW(nsp));
end:

R1R1 := proc (s,nsp)
evalb(s=R1 and nsp[2]+nsp[3]+nsp[5]=aLL(nsp) and 0<nsp[2]+nsp[5]);
end:
R1R2 := proc (s,nsp)
evalb(s=R1 and aLL(nsp)=nsp[3]);
end:
R1Y1 := proc (s,nsp)
evalb(s=R1 and 0<nsp[1]+rD(nsp));
end:
R1Y2 := proc (s,nsp)
evalb(s=R1 and aLL(nsp)=yW(nsp));
end:

R2R2 := proc (s,nsp)
evalb(s=R2 and aLL(nsp)=nsp[3]+nsp[4] and 0<nsp[3]);
end:
R2R3 := proc (s,nsp)
evalb(s=R2 and aLL(nsp)=nsp[4]);
end:
R2Y1 := proc (s,nsp)
evalb(s=R2 and 0<nsp[1]+rD(nsp));
end:
R2Y2 := proc (s,nsp)
evalb(s=R2 and aLL(nsp)=yW(nsp));
end:

R3R3 := proc (s,nsp)
evalb(s=R3 and aLL(nsp)=nsp[4]+nsp[5] and 0<nsp[4]);
end:
R3R1 := proc (s,nsp)
evalb(s=R3 and aLL(nsp)=nsp[5]);
end:
R3Y1 := proc (s,nsp)
evalb(s=R3 and 0<nsp[1]+rD(nsp));
end:
R3Y2 := proc (s,nsp)
evalb(s=R3 and aLL(nsp)=yW(nsp));
end:

XnY1 := proc (s,nsp)
evalb(1<rD(nsp));
end:

bL := proc (nsp) nsp[3]+nsp[4]+nsp[5]; end:
rD := proc (nsp) nsp[6]+nsp[7]+nsp[8]; end:
yW := proc (nsp) nsp[1]+nsp[2]; end:
aLL := proc (nsp) yW(nsp)+bL(nsp)+rD(nsp); end:
sINK := proc (nsp) nsp[1]; end:

3. Building structures and collecting data.

These functions support the recording of the global states' edge-entropy.

> makeSpectrum := proc () local i,j;
global edgeCount,eLabels,G,lEmatrix,nS,S,States,Spectrum;
makeElabels();
makeEmatrix();
Spectrum :=
[seq(seq(evalf(lEmatrix[i,j]/edgeCount),j=i..nS),i=1..nS)];
end:
makeEmatrix := proc () local e,iE,ll,mv,zeros;
global G,eLabels,lEmatrix,nS,S,States;
nS := nops(States);
iE := map(sort, map(e->[ndex(e[1]),ndex(e[2])],eLabels));
zeros := [seq(0, mv=1..nS^2)];
lEmatrix := matrix(nS,nS,zeros);
for ll in iE do
lEmatrix[ll[1],ll[2]] := 1+lEmatrix[ll[1],ll[2]]; od;
lEmatrix;
end:
makeElabels := proc () local e; global eLabels,E,S;
eLabels := map(e->[S[e[1]],S[e[2]]], E);
end:

The following functions construct the (somewhat) random tissue

as a set of cells (ie., vertices) communicating with negihbors.

Cells are randomly selected points on the 100x100 torus

subject to the condition that no two are less than exR apart.

When cells c,d satisfy exR<distance(c,d)<inR they are connected by an edge.

G is the list of neighborhoods -- G[c] is the list of cells connected to c by an edge.

E is the set of edges.

CCxy makes the set C of N cells so that no two are with in exR of each other,

and the locations Cxy of these cells.

> makeG := proc () global exR,inR,N,C,Cxy,G;
G := map(neighbors,C);
end:
makeE := proc () local c,d; global C,E,edgeCount,G;
E := [];
for c in C do
for d in G[c] do
if c<d then E:=[op(E),[c,d]]; fi;od;od;
edgeCount := nops(E);
E;
end:
makeS := proc () local c; global C,N,States,S;
S := array(1..N);
for c in C do S[c]:=randEle(States[1..-3]); od;
evalm(S);
end:
neighbors := proc (p::posint) local neigh,q; global exR,inR,C,Cxy;
neigh := [];
for q in C do
if edge(p,q) then neigh:=[op(neigh),q]; fi;od;
neigh;
end:
edge := proc (p::posint,q::posint) local pqD; global Cxy,exR,inR;
pqD := torusDist(Cxy[p],Cxy[q]);
(exR<pqD and pqD<=inR);
end:
makeCCxy := proc () local i,n,xy; global exR,N,C,Cxy;
C := [seq(i,i=1..N)]:
Cxy := [randLoc()];
while evalb(nops(Cxy)<N) do
xy := randLoc();
if farEnough(xy) then
Cxy:=[op(Cxy),xy]; fi;od;
Cxy;
end:

GO, CO, COxy, and EOxy are the corresponding structures with the edge-overlap seen in the diagrams.

This overlap aids in understanding the influences of cells near the seams of the torus.

> makeCCOxy := proc () local c,cd,n,xy; global CO,Cxy,COxy,inR,overLap;
c := 100; cd := c-inR;
overLap := [];
for xy in Cxy do
if xy[1]<inR
then if xy[2]<inR then
overLap:=[op(overLap), [xy[1]+c, xy[2]+c]];
elif cd<xy[2] then
overLap:=[op(overLap), [xy[1]+c, xy[2]-c]];
else
overLap:=[op(overLap), [xy[1]+c, xy[2] ]];fi;
elif cd<=xy[1] then
if xy[2]<inR then
overLap:=[op(overLap),[xy[1]-c, xy[2]+c]];
elif cd<=xy[2] then
overLap:=[op(overLap),[xy[1]-c, xy[2]-c]];
else
overLap:=[op(overLap),[xy[1]-c, xy[2] ]];fi;
elif xy[2]<inR then
overLap:=[op(overLap), [xy[1], xy[2]+c]];
elif cd<=xy[2] then
overLap:=[op(overLap), [xy[1], xy[2]-c]];fi;
od;
COxy := e2sort([op(Cxy), op(overLap)]);
CO := map(cellName,COxy);
COxy;
end:
makeGO := proc () global CO,COxy,GO,inR;
GO := map(sort,map(reduce,map(neighborsCO,CO)));
end:
makeEOxy := proc () local e,nO,p,q; global COxy,Cxy,EO,EOxy,inR,eLabels;
EOxy := []; nO := nops(COxy);
for p from 1 to nO-1 do
for q from p+1 to nO do
if edgeCO(p,q)
then
EOxy := [op(EOxy),[COxy[p],COxy[q]]];fi;
od;od;
EO := map(e->[cellName(e[1]),cellName(e[2])],EOxy);
end:
neighborsCO := proc (p) local neigh,q; global CO,COxy,inR;
neigh := [];
for q in CO do
if edgeCO(p,q) then neigh:=[op(neigh),q]; fi;od;
neigh;
end:
edgeCO := proc (p,q) local pqD; global COxy,inR;
pqD := eDist(COxy[p],COxy[q]);
(0<pqD and pqD<=inR);
end:

These low-level functions support those above.

> cellName := proc (xy) local name,uv; global Cxy;
uv := [xy[1] mod 100, xy[2] mod 100];
member(uv,Cxy,'name');
name;
end:
farEnough := proc (p) local c,bv; global exR,Cxy;
bv := true;
for c from 1 to nops(Cxy) do
bv := (bv and evalb(exR<torusDist(Cxy[c],p)));od;
bv;
end:
torusDist := proc (c1,c2)
evalf(eDist([0,0], C2Cvector(c1,c2)));
end:
eDist := proc (c1,c2)
evalf(sqrt((c1[1]-c2[1])^2 + (c1[2]-c2[2])^2));
end:
C2Cvector := proc (c1,c2)
[min(abs(c1[1]-c2[1]), 100-abs(c1[1]-c2[1])),
min(abs(c1[2]-c2[2]), 100-abs(c1[2]-c2[2]))];
end:

> subtract := proc (s1,s2) local c,s3;
s3 := [];
for c in s1 do
if not(member(c,s2)) then s3 := [op(s3),c]; fi;od;
s3;
end:
reduce := proc (lst) local i,p,rlist;
rlist := [];
for i from 1 to nops(lst) do
p := lst[i];
if not(member(p,lst[i+1..-1]))
then rlist:=[op(rlist),p]; fi;od;
rlist;
end:

randLoc := () -> [rand(100)(), rand(100)()]:
randEle := proc (lst) local nlst;
nlst := nops(lst);
lst[rand(nlst)()+1]:
end:

e2sort := proc (xyList)
sort(xyList,e2ineq);
end:
e2ineq := proc (uv,xy)
evalb(uv[1]<xy[1] or (uv[1]=xy[1] and uv[2]<=xy[2]));
end:

nSpec := proc (c) local ndx,ns,s; global G,S,States;
ns := [seq(0,ndx=1..nops(States))];
for s in SoG(c) do
member(s,States,'ndx');
ns[ndx] := 1+ns[ndx]; od;
ns;
end:
ndex := proc (s) local ndx;global States;
member(s,States,'ndx');
ndx;
end:

SoG := proc (c) local nc; global G,S;
map(nc->S[nc], G[c]);
end:

makeGG := proc () local c,d,Gc,GGc; global G,GG;
GG := [];
for c in C do
Gc := G[c];
GGc := reduce(map(op, map(d->G[d],Gc)));
GG := [op(GG), sort(subtract(GGc,[op(Gc),c]))]; od;
GG;
end: