{$N+}
{$R-}
{$define CGAMODE}
{ define BW}
{$define HISTORY}
{ define WILSON}
{$define WILSON2}
{ define WINTERBLIGHT}
{$define WILSON3} {$ifdef WILSON3} {$define WILSON2} {$endif} {Wilson3 is a variant of Wilson2}
{ define WEAK} {WEAK or STRONG refers to whether altruistic advantage is reflexive}
{ define SEX}
{ define DOMINANT} {means altruism is dominant; applies only if SEX is $defined}
{ define CORNERS}
program ECOOPVSC; {variant of Wilson's original evolution of cooperation on
                   a viscous grid}

uses GRAPH,CRT,DOS;

{type double=real;}
const mixed_random=0; top_bottom=1; solid_altr=2; solid_self=3; cancer=4; mixed_no_holes=5;
      initspec:byte=mixed_no_holes;  smartsearch:boolean=false;
      zero=true; nozero=false; half=0.5; after:boolean=false;
{$ifdef WEAK} fourorfive=4; {$else} fourorfive=5; {$endif}
      xmax=199; ymax=199; esc=#27; brk=#3; tiny=1E-8;
      empty=0; altr=1; self=2; hetz=3;
      central=0; neighbor=1; enhanced=2;
      {outfilename='BLITFACT.OUT'; logfilename='BLITFACT.LOG'; timelimit:longint=200000;}
      outfilename='OVRL-GEN.OUT'; logfilename='OVRL-GEN.LOG'; timelimit:longint=200000;
      lineacross='=================================================================';
      {parameters: }
{$ifdef WILSON2}
        cullingfreq:integer=00;    {how many generations between winters?}
        cullingfactor:double=0.05; {what percentage of the population survives winter?}
        b:double=0.01;  {benefit per site fr altruists}
        c:double=0.0107;  {cost=rel disadvantage of altruists}
        et:double=0.90; {number of lottery tickets assigned to an empty site}
        xt:double=0.0;
{$else WILSON2}
{$ifdef WILSON}
        f=0.60;  {base filling factor}
        b=0.016; {full benefit, for all altruists}
        c=0.01; {cost=rel disadvantage of altruists}
{$else WILSON}
        f1=0.68; {fitness of altr}
        f2=0.6870; {fitness of self}
{$ifdef DOMINANT}
        fz:double=f1;
{$else  DOMINANT}
        fz:double=f2;
{$endif DOMINANT}
{$ifdef WEAK}
       fcs:double=f1; {central selfish fitness is f1}
       fcz:double=f1;
{$else not WEAK}
       fcs:double=f2; {central selfish fitness is f2}
{$ifdef DOMINANT}
       fcz:double=f1;
{$else DOMINANT}
       fcz:double=f2;
{$endif DOMINANT}
{$endif WEAK}
         h=1.10; {fitness enhancement for neighbor of an altruist}
{$endif WILSON}
{$endif WILSON2}

{$ifdef CGAMODE}
{$ifdef BW}
      color:array[empty..hetz] of byte=(2,0,3,1);
{$else BW}
      color:array[empty..hetz] of byte=(0,1,2,3);
{$endif BW}
{$else CGAMODE}
      color:array[empty..hetz] of byte=(black,lightred,lightblue,magenta);
{$endif}

type
     occupant=empty..hetz;
     relation=central..enhanced;
     map=array[0..ymax,0..xmax] of byte;
     mapptr=^map;

var m,newm,p            :mapptr;
    tavghetz,
    neighborcorr,tavgncorr,
    tavgself,tavgaltr,
    tavgratio,tavgfill  :double;
    selfcum,selfsum,
    altrcum,altrsum,
    hetzcum,hetzsum,
    nsites,time       :longint;
    oneminusc,
    hf1,hf2,hfz,
    afill,sfill,sepn,
    w,z,
    rho1ss,rho2ss     :double;
    rstr,rstr2,
    rstr3,rstr4,rstr5 :string;
    east,west         :array[0..xmax] of integer;
    north,south       :array[0..ymax] of integer;
    probaltr,probfill :array[0..5,0..5] of double;
    fin,fout          :text;
    headerdone        :boolean;
    ki                :char;

function TwoStr(x :word; zero :boolean):string;
         var ws :string;
         begin
         Str(x:2,ws);
         if (ws[1]=' ') then
           if (zero) then ws[1]:='0' else ws:=ws[2];
         TwoStr:=ws;
         end;

function TimeStr:string;
         var hr,mn,sc,s100,yr,mo,da,dow :word;
         begin
         GetTime(hr,mn,sc,s100);
         GetDate(yr,mo,da,dow);
         TimeStr:=twostr(mo,nozero)+'/'+twostr(da,nozero)+'/'+twostr(yr mod 100,zero)
                 +' at '+twostr(hr,nozero)+':'+twostr(mn,zero);
         end;

procedure VGAMode;
          var driver, mode: integer;
          begin
{$ifdef CGAMODE}
          driver:=CGA; mode:=CGAC1;
{$else}
          driver:=VGA; mode:=VGAHI;
{$endif}
          InitGraph(driver,mode,'C:\BP\BGI');
          if (grerror<>0) then writeln(grerror,': ',GraphErrorMsg(grerror));
          if (xmax>=150) then SetTextStyle(smallfont,horizdir,4);
          end;

procedure InitProbArrays;
          var altrcount,selfcount :integer;
              norm,pb             :double;
          begin
          for altrcount:=0 to 5 do begin
            pb:=b*altrcount;
            for selfcount:=0 to (5-altrcount) do begin
              if (altrcount=0) and (selfcount=0) then begin
                probaltr[altrcount,selfcount]:=0;
                probfill[altrcount,selfcount]:=0;
                end
              else begin
                {pb:=5*b*altrcount/(altrcount+selfcount);}
                norm:=altrcount*(1+pb-c) + selfcount*(1+pb) + (5-altrcount-selfcount)*et + xt;
                probaltr[altrcount,selfcount]:=altrcount*(1+pb-c)/norm;
                probfill[altrcount,selfcount]:=probaltr[altrcount,selfcount] + selfcount*(1+pb)/norm;
                end;
              end;
            end;
          end;

procedure InitPopulation;
          var y,x :integer;
          begin
          case initspec of
            mixed_random : begin
                           for y:=0 to ymax do
                             for x:=0 to xmax do
                               m^[y,x]:=random(3);
                           end;
          mixed_no_holes : begin
                           for y:=0 to ymax do
                             for x:=0 to xmax do
                               m^[y,x]:=1+random(2);
                           end;
            top_bottom : begin
                         for y:=0 to ymax div 2 do
                           for x:=0 to xmax do
                             if (random(4)=0) then m^[y,x]:=empty else m^[y,x]:=altr;
                         for y:=succ(ymax div 2) to ymax do
                           for x:=0 to xmax do
                             if (random(4)=0) then m^[y,x]:=empty else m^[y,x]:=self;
                             {if (random<0.5) then m^[y,x]:=empty else m^[y,x]:=altr;}
                         end;
           solid_self : begin
                        for y:=0 to ymax do
                           for x:=0 to xmax do
                             if (random<0.3) then m^[y,x]:=0 else m^[y,x]:=self;
                        end;
    cancer,solid_altr : begin
                        for y:=0 to ymax do
                           for x:=0 to xmax do
                             if (random<0.3) then m^[y,x]:=0 else m^[y,x]:=altr;
                        if (initspec=cancer) then
                           for y:=ymax div 2 - 5 to ymax div 2 + 5 do
                             for x:=xmax div 2 - 5 to xmax div 2 + 5 do
                               m^[y,x]:=self;
                        end;
            end; {case initspec}
          end;

procedure OpenFileAndWriteHeader(fname :string);
          begin
          assign(fout,fname); {$I-} append(fout); {$I+}
          if (ioresult>0) then rewrite(fout)
          else writeln(fout,lineacross);
          writeln(fout,succ(xmax),' * ',succ(ymax),' grid.  '+TimeStr);
          {$ifdef WILSON} write(fout,'WILSON '); {$endif}
          {$ifdef WILSON3} write(fout,'WILSON3 '); {$else}
          {$ifdef WILSON2} write(fout,'WILSON2 '); {$endif} {$endif}
          {$ifdef WEAK} write(fout,'WEAK '); {$endif}
          {$ifdef SEX} write(fout,'SEX '); {$endif}
          {$ifdef DOMINANT} write(fout,'DOMINANT '); {$endif}
          {$ifdef CORNERS} write(fout,'CORNERS '); {$endif}
          write(fout,'Grid Initialization: ');
          case initspec of
             mixed_random : write(fout,'mixed random');
               top_bottom : write(fout,'top/bottom');
               solid_self : write(fout,'solid SELF');
               solid_altr : write(fout,'solid ALTR');
                   cancer : write(fout,'ALTR w/ SELF cancer');
           mixed_no_holes : write(fout,'mixed, no holes');
              end; {case initspec}
          writeln(fout);
{$ifdef WILSON2}
          if (cullingfreq=0) then writeln('No winters')
          else begin
{$ifdef WINTERBLIGHT}
            write(fout,'Block ');
{$endif WINTERBLIGHT}
            writeln(fout,'Culling: ',cullingfactor:5:3,' every ',cullingfreq,' generations.');
            end;
          writeln(fout,'b=',b:7:5,'  c=',c:7:5,'  et=',et:5:3,'  xt=',xt:5:3);
{$endif WILSON2}
          writeln(fout,'Theo SS filling:  Altr=',rho1ss:5:3,'   Self=',rho2ss:5:3);
          writeln(fout);
          write(fout,'   TIME    ALTR    SELF');
{$ifdef SEX}
          write(fout,'    HETZ');
{$endif SEX}
          if (smartsearch) then writeln(fout,'  b TA(ln(ratio)) TA(fill) Rho1SS  g  TA(Ncorr)')
          else writeln(fout,'    NCorr');
          end;

procedure InitInit;
          begin
          VGAMode;
          Randomize;
          new(m); new(newm); new(p);
          nsites:=succ(xmax)*succ(ymax);
          end;

procedure Init;
          var x,y :integer;
              ws  :string;
              aa,bb,cc :double;
          begin
          ki:=#0;
{$ifdef WILSON2}
          InitProbArrays;
          oneminusc:=1-c;
          Str(c:7:5,ws); rstr:='c: '+ws;
          Str(b:7:5,ws); rstr:=rstr+' b: '+ws;
          Str(et:5:3,ws);
          rstr2:='Fitness of void: '+ws;
          Str(cullingfactor:5:3,ws);
          Str(xt:5:3,ws);
          if (xt<>0) then
          rstr3:='Invisible holes: '+ws
          else rstr3:='';
          if (cullingfreq=0) then begin
            rstr4:='No winters'; rstr5:='';
            aa:=4*b; {In this new version, every A is guaranteed its own b}
            bb:=(1-c+b-et-4*b)*half;  {Total is (1+4*rho)*b}
            cc:=et+xt/5-1+c-b;
            {aa:=fourorfive*b;  OLD VERSION - SLIGHTLY DIFFERENT.
            bb:=(1-c-fourorfive*b-et)/2;    Total is (5*rho)*b
            cc:=c-1+et+xt/fourorfive;}
            if (sqr(bb)-aa*cc>=0) and (aa<>0) then
              rho1ss:=(-bb+sqrt(sqr(bb)-aa*cc))/aa
            else rho1ss:=0;
            if (et>=1) then rho2ss:=0
            else rho2ss:=1-(xt/fourorfive)/(1-et);
            end
          else begin
            Str(cullingfactor:5:3,ws);
{$ifdef WINTERBLIGHT} rstr4:='Block Culling: '+ws; {$else} rstr4:='Culling: '+ws; {$endif}
            Str(cullingfreq,ws);
            rstr5:='every '+ws+' genrns';
            end;
{$else WILSON2}
{$ifdef WILSON}
          oneminusc:=1-c;
          rho1ss:=f*(1+b);
          rho2ss:=f;
          Str(c:7:5,ws); rstr:='c: '+ws;
          Str(b:7:5,ws); rstr:=rstr+' b: '+ws;
          Str(f:5:3,ws);
          rstr2:='Base filling f: '+ws;
{$else WILSON}
          hf1:=h*f1;
          hf2:=h*f2;
          hfz:=h*fz;
          rho1ss:=1/f1 + 1/hf1 - 1/(f1*hf1);
          rho2ss:=2/f2 - 1/sqr(f2);
          Str(f1:5:3,ws); rstr:='f1: '+ws;
          Str(f2:5:3,ws); rstr:=rstr+'  f2: '+ws;
          Str(h:5:3,ws);
          rstr2:='';
{$ifdef CGAMODE}
          if (xmax>150) then rstr2:='h: '+ws else
{$endif CGAMODE}
          rstr:=rstr+'  h: '+ws;
{$endif WILSON}
{$endif WILSON2}
          time:=0;
          for y:=0 to ymax do begin
              north[y]:=pred(y);
              south[y]:=succ(y);
              end;
          north[0]:=ymax; south[ymax]:=0;
          for x:=0 to xmax do begin
              west[x]:=pred(x);
              east[x]:=succ(x);
              end;
          west[0]:=xmax; east[xmax]:=0;
          InitPopulation;
{$ifdef HISTORY}
          OpenFileAndWriteHeader(outfilename);
{$endif}
          selfcum:=0; altrcum:=0; hetzcum:=0;
          end;

procedure ComputeStats;
          var x,y,eastx,southy  :integer;
              sssum,aasum,assum,
              essum,aesum,eesum :longint;
              denom,
              xavg,xxavg,xyavg,
              rhobar,newratio   :double;
          begin
          altrsum:=0; selfsum:=0; hetzsum:=0;
          aasum:=0; eesum:=0; sssum:=0; aesum:=0; assum:=0; essum:=0;
          for y:=0 to ymax do begin
            southy:=south[y];
            for x:=0 to xmax do begin
              eastx:=east[x];
              case m^[y,x] of self : begin
                                     inc(selfsum);
                                     case m^[southy,x] of
                                        self : inc(sssum);
                                        altr : inc(assum);
                                       empty : inc(essum);
                                         end; {case}
                                     case m^[y,eastx] of
                                        self : inc(sssum);
                                        altr : inc(assum);
                                       empty : inc(essum);
                                         end; {case}
                                     end;
                              altr : begin
                                     inc(altrsum);
                                     case m^[southy,x] of
                                        self : inc(assum);
                                        altr : inc(aasum);
                                       empty : inc(aesum);
                                         end; {case}
                                     case m^[y,eastx] of
                                        self : inc(assum);
                                        altr : inc(aasum);
                                       empty : inc(aesum);
                                         end; {case}
                                     end;
                              hetz : inc(hetzsum);
                             empty : begin
                                     case m^[southy,x] of
                                        self : inc(essum);
                                        altr : inc(aesum);
                                       empty : inc(eesum);
                                         end; {case}
                                     case m^[y,eastx] of
                                        self : inc(essum);
                                        altr : inc(aesum);
                                       empty : inc(eesum);
                                         end; {case}
                                     end;
                              end; {case}
              end; {for x}
            end; {for y}
          if (aasum+aesum=0) then afill:=0
          else afill:=aasum/(aasum+aesum);
          if (sssum+essum=0) then sfill:=0
          else sfill:=sssum/(sssum+essum);
          {neighborcorr= (<xy>-<x><y>)/(<xý>-<x>ý) where <x>=<y>}
          xyavg:=(sssum+aasum-assum)/(2*nsites);
          xxavg:=(altrsum+selfsum)/nsites;
          xavg:=(altrsum-selfsum)/nsites;
          denom:=xxavg-sqr(xavg);
          if (abs(denom)<tiny) then neighborcorr:=-1
          else neighborcorr:=(xyavg-sqr(xavg)) / denom;
          if (altrsum*selfsum=0) or (aasum+sssum+assum=0) then sepn:=-1
          else begin
            w:=(aasum+sssum)/(aasum+sssum+assum);
            w:=w * sqr(altrsum+selfsum)/(2*altrsum*selfsum);
            sepn:=(w-1)/(w+1);
            if (assum>0) then w:=0.5*(aasum+sssum)/assum
            else w:=-1;
            rhobar:=0.5*(rho1ss+rho2ss);
            if (rhobar<1) then z:=eesum/(2*nsites*sqr(1-rhobar))
            else z:=-1
            end;
          rhobar:=(altrsum+selfsum {$ifdef SEX}+hetzsum{$endif})/nsites;
          if (selfsum>0) and (altrsum>0) then newratio:=ln(altrsum/selfsum)
          else if (altrsum>0) then newratio:=nsites
          else newratio:=-nsites;
          if (time<=100) then begin
            tavgratio:=newratio;
            tavgncorr:=neighborcorr;
            tavgfill:=rhobar;
            tavgself:=selfsum; tavgaltr:=altrsum; {$ifdef SEX}tavghetz:=hetzsum;{$endif}
            end
          else begin
            tavgratio:=0.98*tavgratio+0.02*newratio;
            tavgncorr:=0.98*tavgncorr+0.02*neighborcorr;
            tavgfill:=0.98*tavgfill+0.02*rhobar;
            tavgself:=0.98*tavgself+0.02*selfsum; tavgaltr:=0.98*tavgaltr+0.02*altrsum;
            {$ifdef SEX} tavghetz:=0.98*tavghetz+0.02*hetzsum;{$endif}
            end;
          end; {ComputeStats}

procedure Accumulate;
          begin
          inc(selfcum,selfsum);
          inc(altrcum,altrsum);
          end;

{$ifdef WILSON2}
{$else WILSON2}
{$ifdef WILSON}
{$else WILSON}
{$ifdef SEX}
function Seed(source :occupant; rel :relation):occupant;
         begin
         if (source=empty) then Seed:=empty
         else case rel of
           central : case source of
               altr : begin
                      if (random<f1) then Seed:=altr else Seed:=empty;
                      end;
               hetz : begin
                      if (random<fcz) then Seed:=hetz else Seed:=empty;
                      end;
               self : begin
                      if (random<fcs) then Seed:=self else Seed:=empty;
                      end;
             end; {central}
           neighbor : case source of
               altr : begin
                      if (random<f1) then Seed:=altr else Seed:=empty;
                      end;
               hetz : begin
                      if (random<fz) then Seed:=hetz else Seed:=empty;
                      end;
               self : begin
                      if (random<f2) then Seed:=self else Seed:=empty;
                      end;
             end; {neighbor}
           enhanced : case source of
               altr : begin
                      if (random<hf1) then Seed:=altr else Seed:=empty;
                      end;
               hetz : begin
                      if (random<hfz) then Seed:=hetz else Seed:=empty;
                      end;
               self : begin
                      if (random<hf2) then Seed:=self else Seed:=empty;
                      end;
             end; {enhanced}
           end; {case rel}
         end;
{$else SEX}
function Seed(source :occupant; rel :relation):occupant;
         begin
         if (source=empty) then Seed:=empty
         else case rel of
           central : begin
             if (source=altr) then begin
                if (random<f1) then Seed:=altr else Seed:=empty;
                end
             else {source=self} begin
                if (random<fcs) then Seed:=self else Seed:=empty;
                end;
             end; {central}
           neighbor : begin
             if (source=altr) then begin
                if (random<f1) then Seed:=altr else Seed:=empty;
                end
             else {source=self} begin
                if (random<f2) then Seed:=self else Seed:=empty;
                end;
             end; {neighbor}
           enhanced : begin
             if (source=altr) then begin
                if (random<hf1) then Seed:=altr else Seed:=empty;
                end
             else {source=self} begin
                if (random<hf2) then Seed:=self else Seed:=empty;
                end;
             end; {enhanced}
           end; {case}
         end;
{$endif SEX}
{$endif WILSON else}
{$endif WILSON2 else}

(*
{$ifdef CORNERS}
function RandomCorner(y,x :integer):occupant;
         var x1,y1 :integer;
         begin
         case random(4) of
            0 : begin
                if x=0 then x1:=xmax else x1:=pred(x);
                if y=0 then y1:=ymax else y1:=pred(y);
                end;
            1 : begin
                if x=0 then x1:=xmax else x1:=pred(x);
                if y=ymax then y1:=0 else y1:=succ(y);
                end;
            2 : begin
                if x=xmax then x1:=0 else x1:=succ(x);
                if y=0 then y1:=ymax else y1:=pred(y);
                end;
           else {3} begin
                if x=xmax then x1:=0 else x1:=succ(x);
                if y=ymax then y1:=0 else y1:=succ(y);
                end;
           end; {case}
         RandomCorner:=Seed(y1,x1,(m^[y,x]=altr));
         end;
{$endif CORNERS}
*)
{$ifdef WILSON}
procedure IncrCount(y,x :integer; var altrcount,selfcount :integer);
          begin
          case m^[y,x] of
             self : inc(selfcount);
             altr : inc(altrcount);
             end;
          end;

function SomeoneSeeds(y,x :integer):boolean;
         var altrcount,selfcount,totalcount :integer;
             prob                           :double;
         begin
         altrcount:=0; selfcount:=0;
         IncrCount(y,x,altrcount,selfcount);
         IncrCount(y,east[x],altrcount,selfcount);
         IncrCount(y,west[x],altrcount,selfcount);
         IncrCount(north[y],x,altrcount,selfcount);
         IncrCount(south[y],x,altrcount,selfcount);
         totalcount:=altrcount+selfcount;
         if (totalcount=0) then SomeoneSeeds:=false
         else begin
           prob:=f*(1+b*altrcount/totalcount);
           SomeoneSeeds:=(random<prob);
           end;
         end;

procedure IncrWeightedCount(y,x :integer; var altrcount,selfcount :double);
          begin
          case m^[y,x] of
             self : selfcount:=selfcount+1;
             altr : altrcount:=altrcount+oneminusc;
             end;
          end;

function RandomNeighbor(y,x :integer):occupant;
         var altrcount,selfcount,totalcount :double;
         begin
         altrcount:=0; selfcount:=0;
         IncrWeightedCount(y,x,altrcount,selfcount);
         IncrWeightedCount(y,east[x],altrcount,selfcount);
         IncrWeightedCount(y,west[x],altrcount,selfcount);
         IncrWeightedCount(north[y],x,altrcount,selfcount);
         IncrWeightedCount(south[y],x,altrcount,selfcount);
         totalcount:=altrcount+selfcount;
         {if (totalcount=0) then RandomNeighbor:=empty    never happens
         else}if (selfcount=0) then RandomNeighbor:=altr
         else if (altrcount=0) then RandomNeighbor:=self
         else if (altrcount>totalcount*random) then RandomNeighbor:=altr
         else RandomNeighbor:=self;
         end;

procedure NextTimeStep;
          var x,y    :integer;
              swap   :mapptr;
              rel    :relation;
          begin
          for y:=0 to ymax do
            for x:=0 to xmax do begin
              if SomeoneSeeds(y,x)
              then newm^[y,x]:=RandomNeighbor(y,x)
              else newm^[y,x]:=empty;
              end;
          swap:=m; m:=newm; newm:=swap; inc(time);
{$ifdef HISTORY}
          Accumulate;
{$endif HISTORY}
          end;
{$else WILSON}
{$ifdef WILSON3}
procedure ComputeFitnessArray;
          var x,y,pp          :integer;
          begin
          for y:=0 to ymax do
            for x:=0 to xmax do begin
{$ifndef WEAK}if (m^[y,x]=altr) then pp:=1 else {$endif} pp:=0;
              if (m^[north[y],x]=altr) then inc(pp);
              if (m^[south[y],x]=altr) then inc(pp);
              if (m^[y,east[x]]=altr) then inc(pp);
              if (m^[y,west[x]]=altr) then inc(pp);
{$ifdef CORNERS}
              if (m^[north[y],east[x]]=altr) then inc(pp);
              if (m^[south[y],west[x]]=altr) then inc(pp);
              if (m^[south[y],east[x]]=altr) then inc(pp);
              if (m^[north[y],west[x]]=altr) then inc(pp);
{$endif}
              p^[y,x]:=pp;
              end;
          end;
{$endif WILSON3}

{$ifdef WILSON2}
(*
function RandomNeighbor(y,x :integer):occupant;
         {this version holds a weighted lottery - unnecessarily slow}
         var r,norm       :double;
             pb           :double; {partial benefit}
             orderedneighbor :array[0..4] of occupant; {0,1,2,3,4=C,N,S,E,W}
             tix,cumtix   :array[0..4] of double;       {0,1,2,3,4=C,N,S,E,W}
             i,altrcount,selfcount,totalcount  :integer;

   function TixFunc(o :occupant):double;
            begin
            case o of
              empty : TixFunc:=et;
               altr : TixFunc:=1-c+pb;
               self : TixFunc:=1+pb;
               end;
            end;

         begin
         orderedneighbor[0]:=m^[y,x];
         orderedneighbor[1]:=m^[north[y],x];
         orderedneighbor[2]:=m^[south[y],x];
         orderedneighbor[3]:=m^[y,east[x]];
         orderedneighbor[4]:=m^[y,west[x]];
         altrcount:=0; selfcount:=0;
         for i:=0 to 4 do
           case orderedneighbor[i] of
              self : inc(selfcount);
              altr : inc(altrcount);
              end;
         totalcount:=altrcount+selfcount;
         if (totalcount=0) then RandomNeighbor:=empty
         else begin
           pb:=5*b*altrcount/totalcount;
           {or pb:=b*altrcount;}
           for i:=0 to 4 do
             tix[i]:=TixFunc(orderedneighbor[i]);
           cumtix[0]:=tix[0];
           for i:=1 to 4 do cumtix[i]:=cumtix[pred(i)]+tix[i];
           norm:=cumtix[4];
           r:=random*norm;
           for i:=0 to 3 do
             if (r<cumtix[i]) then begin
                RandomNeighbor:=orderedneighbor[i];
                exit;
                end;
           RandomNeighbor:=orderedneighbor[4];
           end;{totalcount>0}
         end;
*)
{$ifdef WILSON3}
function SimulatedLottery(y,x :integer):occupant;
         var altrcount,selfcount,emptycount,
             r,sum                           :double;
         begin
         altrcount:=0; selfcount:=0; emptycount:=xt;
         case m^[y,x] of
           self : selfcount:=selfcount+1.1+p^[y,x]*b;
           altr : altrcount:=altrcount+1.1-c+p^[y,x]*b;
           else emptycount:=emptycount+et;
           end;
         case m^[north[y],x] of
           self : selfcount:=selfcount+1+p^[north[y],x]*b;
           altr : altrcount:=altrcount+1-c+p^[north[y],x]*b;
           else emptycount:=emptycount+et;
           end;
         case m^[south[y],x] of
           self : selfcount:=selfcount+1+p^[south[y],x]*b;
           altr : altrcount:=altrcount+1-c+p^[south[y],x]*b;
           else emptycount:=emptycount+et;
           end;
         case m^[y,east[x]] of
           self : selfcount:=selfcount+1+p^[y,east[x]]*b;
           altr : altrcount:=altrcount+1-c+p^[y,east[x]]*b;
           else emptycount:=emptycount+et;
           end;
         case m^[y,west[x]] of
           self : selfcount:=selfcount+1+p^[y,west[x]]*b;
           altr : altrcount:=altrcount+1-c+p^[y,west[x]]*b;
           else emptycount:=emptycount+et;
           end;
{$ifdef CORNERS}
         case m^[north[y],east[x]] of
           self : selfcount:=selfcount+1+p^[north[y],x]*b;
           altr : altrcount:=altrcount+1-c+p^[north[y],x]*b;
           else emptycount:=emptycount+et;
           end;
         case m^[south[y],west[x]] of
           self : selfcount:=selfcount+1+p^[south[y],x]*b;
           altr : altrcount:=altrcount+1-c+p^[south[y],x]*b;
           else emptycount:=emptycount+et;
           end;
         case m^[south[y],east[x]] of
           self : selfcount:=selfcount+1+p^[y,east[x]]*b;
           altr : altrcount:=altrcount+1-c+p^[y,east[x]]*b;
           else emptycount:=emptycount+et;
           end;
         case m^[north[y],west[x]] of
           self : selfcount:=selfcount+1+p^[y,west[x]]*b;
           altr : altrcount:=altrcount+1-c+p^[y,west[x]]*b;
           else emptycount:=emptycount+et;
           end;
{$endif CORNERS}
         sum:=altrcount+selfcount+emptycount;
         r:=random;
         if (r<altrcount/sum) then SimulatedLottery:=altr
         else if (r<(altrcount+selfcount)/sum) then SimulatedLottery:=self
         else SimulatedLottery:=empty;
         end;
{$else WILSON3}
function SimulatedLottery(y,x :integer):occupant;
         var i,altrcount,selfcount  :integer;
             r                      :double;
         begin
         altrcount:=0; selfcount:=0;
         case m^[y,x] of
           self : inc(selfcount);
           altr : inc(altrcount);
           end;
         case m^[north[y],x] of
           self : inc(selfcount);
           altr : inc(altrcount);
           end;
         case m^[south[y],x] of
           self : inc(selfcount);
           altr : inc(altrcount);
           end;
         case m^[y,east[x]] of
           self : inc(selfcount);
           altr : inc(altrcount);
           end;
         case m^[y,west[x]] of
           self : inc(selfcount);
           altr : inc(altrcount);
           end;
         r:=random;
         if (r<probaltr[altrcount,selfcount]) then SimulatedLottery:=altr
         else if (r<probfill[altrcount,selfcount]) then SimulatedLottery:=self
         else SimulatedLottery:=empty;
         end;
{$endif WILSON3}
{$else WILSON2}
function RandomNeighbor(y,x :integer; rel :relation):occupant;
         begin
         case random(4) of
           0 : RandomNeighbor:=Seed(m^[y,east[x]],rel);
           1 : RandomNeighbor:=Seed(m^[y,west[x]],rel);
           2 : RandomNeighbor:=Seed(m^[north[y],x],rel)
           else {3} RandomNeighbor:=Seed(m^[south[y],x],rel);
           end; {case}
         end;
(*
function RandomNeighbor(y,x :integer; rel :relation):occupant;
         var x1,y1 :integer;
         begin
         case random(4) of
           0 : begin
               y1:=y;
               if x=0 then x1:=xmax else x1:=pred(x);
               end;
           1 : begin
               y1:=y;
               if x=xmax then x1:=0 else x1:=succ(x);
               end;
           2 : begin
               x1:=x;
               if y=0 then y1:=ymax else y1:=pred(y);
               end;
           else {3} begin
               x1:=x;
               if y=ymax then y1:=0 else y1:=succ(y);
               end;
            end; {case}
         if (m^[y,x]=altr) then RandomNeighbor:=Seed(m^[y1,x1],enhanced)
         else RandomNeighbor:=Seed(m^[y1,x1],neighbor);
         end;
*)
{$endif WILSON2}




{$ifdef SEX}
procedure TwoRandomNeighbors(y,x :integer; rel :relation; var parent1,parent2 :occupant);
          var x1,y1 :integer;
          begin
          case random(6) of
            0 : begin
                parent1:=Seed(m^[y,east[x]],rel);
                parent2:=Seed(m^[y,west[x]],rel);
                end;
            1 : begin
                parent1:=Seed(m^[y,east[x]],rel);
                parent2:=Seed(m^[north[y],x],rel);
                end;
            2 : begin
                parent1:=Seed(m^[y,east[x]],rel);
                parent2:=Seed(m^[south[y],x],rel);
                end;
            3 : begin
                parent1:=Seed(m^[y,west[x]],rel);
                parent2:=Seed(m^[north[y],x],rel);
                end;
            4 : begin
                parent1:=Seed(m^[y,west[x]],rel);
                parent2:=Seed(m^[south[y],x],rel);
                end;
            5 : begin
                parent1:=Seed(m^[south[y],x],rel);
                parent2:=Seed(m^[north[y],x],rel);
                end;
            end; {case}
          if (parent1=empty) then parent1:=RandomNeighbor(y,x,rel);
          if (parent2=empty) then parent2:=RandomNeighbor(y,x,rel);
          end; {TwoRandomNeighbors}

function Cross(a,b :occupant):occupant;
{empty=0; altr=1; self=2; hetz=3;}
         begin
         case a*10+b of
           0,1,2,3,
           10,20,30 : Cross:=empty;
                 11 : Cross:=altr;
                 22 : Cross:=self;
                 33 : case random(4) of
                        0 : Cross:=self;
                        1 : Cross:=altr;
                        else {2,3} Cross:=hetz;
                        end; {case}
              12,21 : Cross:=hetz;
              13,31 : case random(2) of
                        0 : Cross:=hetz;
                        1 : Cross:=altr;
                        end;
              23,32 : case random(2) of
                        0 : Cross:=hetz;
                        1 : Cross:=self;
                        end;
            end; {case}
         end;

procedure NextTimeStep;
          var x,y             :integer;
              swap            :mapptr;
              rel             :relation;
              parent1,parent2 :occupant;
          begin
          for y:=0 to ymax do
            for x:=0 to xmax do begin
              newm^[y,x]:=Cross(Seed(m^[y,x],central),RandomNeighbor(y,x,central));
              if (newm^[y,x]=empty) then begin
                if (m^[y,x]=altr) {$ifdef DOMINANT} or (m^[y,x]=hetz) {$endif} then rel:=enhanced
                else rel:=neighbor;
                TwoRandomNeighbors(y,x,rel,parent1,parent2);
                newm^[y,x]:=Cross(parent1,parent2);
                end;
{$ifdef CORNERS}
              if (newm^[y,x]=empty) then
                if (random<cc) then newm^[y,x]:=RandomCorner(y,x);
{$endif}
              end;
          swap:=m; m:=newm; newm:=swap; inc(time);
{$ifdef HISTORY}
          Accumulate;
{$endif HISTORY}
          end;

{$else SEX, no WILSON}

{$ifdef WILSON2}
procedure Blight;
          var y0,x0,y1,x1,y2,x2,y,x :integer;
          begin
          y1:=random(succ(ymax));  x1:=random(succ(xmax));
          y2:=y1+round(sqrt(1-cullingfactor)*ymax);  x2:=x1+round(sqrt(1-cullingfactor)*xmax);
          for y0:=y1 to y2 do begin
            y:=y0 mod succ(ymax);
            for x0:=x1 to x2 do begin
              x:=x0 mod succ(xmax);
              m^[y,x]:=empty;
              end;
            end;
          end;

procedure Winter;
          var y,x :integer;
          begin
          for y:=0 to ymax do
            for x:=0 to xmax do
              if (random>cullingfactor) then m^[y,x]:=empty;
          end;

procedure NextTimeStep;
          var x,y    :integer;
              swap   :mapptr;
              rel    :relation;
          begin
          if (cullingfreq>0) and (time mod cullingfreq = 0) then
{$ifdef WINTERBLIGHT} Blight; {$else} Winter; {$endif}
{$ifdef WILSON3}
          ComputeFitnessArray;
{$endif WILSON3}
          for y:=0 to ymax do
            for x:=0 to xmax do
              newm^[y,x]:=SimulatedLottery(y,x); {RandomNeighbor is equivalent but slower}
              {newm^[y,x]:=RandomNeighbor(y,x);}
          swap:=m; m:=newm; newm:=swap; inc(time);
{$ifdef HISTORY}
          Accumulate;
{$endif HISTORY}
          end;

{$else WILSON2}

procedure NextTimeStep;
          var x,y    :integer;
              swap   :mapptr;
              rel    :relation;
          begin
          for y:=0 to ymax do
            for x:=0 to xmax do begin
              newm^[y,x]:=Seed(m^[y,x],central);
              if (newm^[y,x]=empty) then begin
                if (m^[y,x]=altr) then rel:=enhanced else rel:=neighbor;
                newm^[y,x]:=RandomNeighbor(y,x,rel);
                end;
{$ifdef CORNERS}
              if (newm^[y,x]=empty) then
                if (random<cc) then newm^[y,x]:=RandomCorner(y,x);
{$endif}
              end;
          swap:=m; m:=newm; newm:=swap; inc(time);
{$ifdef HISTORY}
          Accumulate;
{$endif HISTORY}
          end;
{$endif WILSON2 (else)}
{$endif SEX (else)}
{$endif WILSON (else)}

procedure OutputToDisk;
          begin
          append(fout);
          writeln(fout,time:7,altrcum:10,selfcum:10,{$ifdef SEX}hetzcum:10,{$endif}neighborcorr:10:3);
          altrcum:=0; selfcum:=0;
          close(fout);
          end;

procedure Display;
          var x,y,ht :integer;
              ws     :string;
          begin
          ComputeStats;
          ClearViewPort;
(*
          SetBkColor(white); SetColor(0);
*)
          for y:=0 to ymax do
            for x:=0 to xmax do
               PutPixel(x,y,color[m^[y,x]]);
          Str(time:5,ws);
          ht:=round(TextHeight(ws)*1.4);
          y:=3;  x:=xmax+3;
          ws:='Time:         '+ws;
          OutTextXY(x,y,ws);   inc(y,ht);
          Str(100*(selfsum+altrsum+hetzsum)/nsites:4:0,ws);
          ws:='Occupancy rate:'+ws+'%';
          OutTextXY(x,y,ws);   inc(y,ht);
          Str(altrsum:5,ws);
          ws:='Tot altruists: '+ws;
          OutTextXY(x,y,ws);   inc(y,ht);
          Str(selfsum:5,ws);
          ws:='Tot defectors: '+ws;
          OutTextXY(x,y,ws);   inc(y,ht);
{$ifdef SEX}
          Str(hetzsum:5,ws);
          ws:='Tot heterozyg: '+ws;
          OutTextXY(x,y,ws);   inc(y,ht);
{$endif SEX}
{$ifndef WILSON2}
          Str(w:5:3,ws);
          ws:='w: '+ws;
          OutTextXY(x,y,ws);   inc(y,ht);
          Str(z:5:3,ws);
          ws:='z: '+ws;
          OutTextXY(x,y,ws);   inc(y,ht);
{$endif}
          if (neighborcorr<=-1) then ws:='undef'
          else Str(neighborcorr:5:3,ws);
          ws:='Neighb corr: '+ws;
          OutTextXY(x,y,ws);   inc(y,ht);
          Str(afill:5:3,ws);
          ws:='Altr filling: '+ws;
          OutTextXY(x,y,ws);   inc(y,ht);
          Str(sfill:5:3,ws);
          ws:='Self filling: '+ws;
          OutTextXY(x,y,ws);   inc(y,2*ht);
          OutTextXY(x,y,rstr); inc(y,ht);
          OutTextXY(x,y,rstr2); inc(y,ht);
{$ifdef WILSON2}
          OutTextXY(x,y,rstr3); if (rstr3>'') then inc(y,2*ht);
          OutTextXY(x,y,rstr4); inc(y,ht);
          OutTextXY(x,y,rstr5); if (rstr5>'') then inc(y,ht);
{$else WILSON2}
{$ifndef WILSON}
{$ifdef WEAK}
          ws:='Weak';
{$else}
          ws:='Strong'
{$endif WEAK else}
          ws:=ws+' selfish adv';
          OutTextXY(x,y,ws);   inc(y,ht);
{$endif not WILSON}
{$endif not WILSON2}
          if (cullingfreq=0) then begin
            Str(100*rho1ss:5:1,ws);
            ws:='Theo SSaltr: '+ws+'%';
            OutTextXY(x,y,ws);   inc(y,ht);
            Str(100*rho2ss:5:1,ws);
            ws:='Theo SSself: '+ws+'%';
            OutTextXY(x,y,ws);   inc(y,ht);
            end;
if after then delay(500);
          end; {Display}

procedure LogFile;
          begin
          if (smartsearch) and (headerdone) then begin
            assign(fout,logfilename); append(fout);
            end
          else OpenFileAndWriteHeader(logfilename);
          headerdone:=true;
          write(fout,time:7,altrsum:8,selfsum:8);
          if (smartsearch) then writeln(fout,'  b=',b:6:4,tavgratio:8:3,tavgfill:8:4,rho1ss:8:4,4*b/(c-b):8:2,tavgncorr:6:3)
          else writeln(fout);
          {assign(fout,logfilename); append(fout);
          case initspec of
            solid_altr :
              writeln(fout,time:7,b:8:3,c:8:3,et:8:3,xt:8:3,altrsum:8,rho1ss:8:3,altrsum/nsites:8:3);
            solid_self :
              writeln(fout,time:7,et:8:3,xt:8:3,selfsum:8,rho2ss:8:3,selfsum/nsites:8:3);
            end;}
          close(fout);
          end;

procedure Shuffle; {zero-viscosity case - destroy all geographic info.  Keep neighbor corr=0}
          var i                 :longint;
              y1,x1,y2,x2       :word;
              swap              :byte;
          begin
          for i:=1 to 2*nsites do begin
            x1:=random(succ(xmax));
            y1:=random(succ(ymax));
            x2:=random(succ(xmax));
            y2:=random(succ(ymax));
            swap:=m^[y1,x1];
            m^[y1,x1]:=m^[y2,x2];
            m^[y2,x2]:=swap;
            end; {for i}
          end;

procedure RunOneCase;
          begin;
          Init;
          repeat
            if (time mod 5=0) then Display;
            {Shuffle;}
{$ifdef HISTORY}
            if (time mod 5=0) then OutputToDisk;
{$endif}
            NextTimeStep;
            if (keypressed) then ki:=readkey;
          until (ki in [esc,brk]) or ((altrsum=0) or (selfsum=0)) or (time>=timelimit);
          {$ifdef HISTORY} OutputToDisk; {$endif}
          LogFile;
          end;

function UseBatchFile:boolean;
         var header :string;
         begin
         UseBatchFile:=false;
         if (paramcount<>1) then exit;
         assign(fin,paramstr(1)); {$I-} reset(fin); {$I+}
         if (ioresult<>0) then exit;
         UseBatchFile:=true;
         readln(fin,header);
         end;

procedure ReadInputLine;
          begin
          readln(fin,b,c,et,xt,cullingfreq,cullingfactor,initspec,timelimit);
          if (b<0) then begin
            b:=-b;
            smartsearch:=true;
            end
          else smartsearch:=false;
          if (b>c) then b:=b/5;
          end;

procedure RegressionCalc(n :integer;xsum,ysum,xysum,x2sum,y2sum :double; var R2,intercept :double);
          var xvar,yvar,A,B :double;
          begin
          xsum:=xsum/n; ysum:=ysum/n; xysum:=xysum/n;
          x2sum:=x2sum/n; y2sum:=y2sum/n;
          xvar:=x2sum-sqr(xsum); yvar:=y2sum-sqr(ysum);
          R2:=sqr(xysum-xsum*ysum)/(xvar*yvar);
          A:=(xysum-xsum*ysum)/xvar;
          B:=ysum-A*xsum;
          intercept:=-B/A; {value of x=g that makes y=ln(ratio)=0}
          end;

procedure HuntForCrossoverC;
          var maxtimelimit :longint;
              incr,g       :double;
              count,n      :integer;
              R2,stalemate_g,
              xsum,ysum,xysum,
              x2sum,y2sum  :double;

   procedure RegisterForRegression;
             var x,y :double;
             begin
             if (abs(tavgratio)>10) then exit;
             inc(n);
             x:=g; y:=tavgratio/time;
             xsum:=xsum+x; ysum:=ysum+y; xysum:=xysum+x*y;
             x2sum:=x2sum+sqr(x); y2sum:=y2sum+sqr(y);
             end;

          begin
          maxtimelimit:=timelimit;
          headerdone:=false;
          timelimit:=500;
          count:=0; 
          n:=0; xsum:=0; ysum:=0; x2sum:=0; y2sum:=0; xysum:=0;
          repeat
            inc(count);
            RunOneCase;
            timelimit:=round(1.2*timelimit);
            if (selfsum=0) then incr:=-0.2*5000/time
            else if (altrsum=0) then incr:=0.2*5000/time
            else begin
              timelimit:=round(1.2*timelimit);
              {incr:=ln(selfsum/altrsum)*500/time;}
              incr:=-tavgratio*1500/time;
              end;
            if (timelimit>maxtimelimit) then timelimit:=maxtimelimit;
            if (incr>0.2) then incr:=0.2
            else if (incr<-0.2) then incr:=-0.2;
            g:=4*b/(c-b) * (1+sqrt(random)*(incr));
            RegisterForRegression;
            if (g<1) then g:=1;
            b:=g*c/(4+g);
            if (timelimit=maxtimelimit) and (altrsum>0) and (selfsum>0) and (count<20) then count:=20;
          until (count=22);
          {assign(fout,logfilename); append(fout); writeln(fout,'ÍÍ> *b:c=',4*b/(c-b):6:2); close(fout);}
          RegressionCalc(n,xsum,ysum,xysum,x2sum,y2sum,R2,stalemate_g);
          assign(fout,logfilename); append(fout);
            writeln(fout,'ÍÍ> stalemate g = ',stalemate_g:7:3,' (Rý=',R2:5:3,')');
          close(fout);
          end;

begin
InitInit;
if (UseBatchFile) then begin
   while not eof(fin) and (b<>0) do begin
     ReadInputLine;
     if (smartsearch) then HuntForCrossoverC
     else RunOneCase;
     end; {while not eof}
   end {if batch}
else RunOneCase;
CloseGraph;
end.

repeat
  if (time mod 5=0) then Display;
{$ifdef HISTORY}
  if (time mod 50=0) then OutputToDisk;
{$endif}
  NextTimeStep;
  if (keypressed) then after:=true;
until ((keypressed) and (ReadKey in [esc,brk])) or
      (altrsum=0) or (selfsum=0);
{$ifdef HISTORY} OutputToDisk; {$endif}
LogFile;
Delay(1000); CloseGraph;
end.