{Recombinant reproduction, haploid individuals. There are three loci. The Coop
locus governs whether the individual cooperates, and has three alleles: C, D,
and S. The Punish locus has two alleles, P and NP. The
MinPun locus causes Cooperators to defect if the level of defection on the
previous period was greater than MinPun.}
unit StrongRDiploidNoLink;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  StdCtrls, XYGraphu, HelpU, Histogram, RAverageU, Menus;

const
  Vision = 100;
  PeriodsPerPlot = 100;
type
  TForm1 = class(TForm)
    RunBtn: TButton;
    Label3: TLabel;
    TotalPeriodsBox: TEdit;
    Label4: TLabel;
    CurrentPeriodBox: TEdit;
    HaltBtn: TButton;
    NumGroupsBox: TEdit;
    Label5: TLabel;
    BBox: TEdit;
    CpBox: TEdit;
    CBox: TEdit;
    GammaBox: TEdit;
    MuBox: TEdit;
    Label8: TLabel;
    Label10: TLabel;
    Label11: TLabel;
    Label12: TLabel;
    Label13: TLabel;
    Phi0Box: TEdit;
    Label14: TLabel;
    QuitBtn: TButton;
    Label15: TLabel;
    MutationRateBox: TEdit;
    Label6: TLabel;
    AgentsPerGroupBox: TEdit;
    Rho0CBox: TEdit;
    Label17: TLabel;
    OpenDialog1: TOpenDialog;
    MinGroupSizeBox: TEdit;
    Label20: TLabel;
    MinSSeedBox: TEdit;
    MaxSSeedBox: TEdit;
    HelpButton: TButton;
    Label19: TLabel;
    Label21: TLabel;
    MinMutatedSBox: TEdit;
    MaxMutatedSBox: TEdit;
    ExecuteButton: TButton;
    ExecutionCycleBox: TEdit;
    Label24: TLabel;
    Label1: TLabel;
    Label9: TLabel;
    Rho0DBox: TEdit;
    Rho0PBox: TEdit;
    Label25: TLabel;
    Label26: TLabel;
    Label22: TLabel;
    Rho0SBox: TEdit;
    PopupMenu1: TPopupMenu;
    procedure RunBtnClick(Sender: TObject);
    procedure HaltBtnClick(Sender: TObject);
    procedure FormClose(Sender: TObject; var Action: TCloseAction);
    procedure QuitBtnClick(Sender: TObject);
    procedure SetParameters;
    procedure HaltBtnMouseDown(Sender: TObject; Button: TMouseButton;
      Shift: TShiftState; X, Y: Integer);
    procedure QuitBtnMouseDown(Sender: TObject; Button: TMouseButton;
      Shift: TShiftState; X, Y: Integer);
    procedure HelpButtonClick(Sender: TObject);
    procedure ExecuteButtonClick(Sender: TObject);
  private
  public
  end;
type
	AgentCoop = (Selfish,Cooperator,Defector);
  AgentStatus = (InGroup,InPool);
	Agent = class(TObject)
		ACoop : AgentCoop;
    APunish,Cooperated,MakeCopy : Boolean;
    AStatus : AgentStatus;
    {Cost of being ostracized used by selfish agent to decide
     how much to shirk}
    SCost,Score,MinPun : Real;
    GMember,Mate,NumCopies : Integer;
		Kill  : Boolean;
    Constructor Init; virtual;
    Destructor Destruct;
    Procedure Copy(A : Agent);
    {Shirking rate as a function of the number of reciprocators
     and the agent's SCost}
    Function ASigma(FractionReciprocator : Real): Real;
    Procedure Mutate;
	end;
  Group = class(TObject)
    GNumber : Integer;
    Fs,Fc,Fd,Fp : Real;
    LastGSigma,GSigma : Real;
    GSize : Integer;
    Constructor Init; virtual;
    Procedure AddMember;
  end;
var
  XYPopComposition,PopTypes : TXYGraph;
  NumAgents,NextAgent : Integer;
  Agents : Array of Agent;
  Groups : Array of Group;
  MaxSCost,MinSCost,MaxMutatedSCost,MinMutatedSCost : Real;
	Periods : Longint;
	NumGroups : Integer;
  TotalPeriods,NextFamily : Longint;
  HaltSimulation : Boolean;
	TotalCount : Longint;
  Rho0C,Rho0D,Rho0P,Rho0S,B,Cp,C,Gamma,Mu,RhoC,RhoS,RhoD,RhoP: Real;
  APerGroup,PopSize,MinGroupSize,TotalInPool : Integer;
  ExecutionCycle,TotalAgents: Integer;
  Phi0 : Real;
  Execute : Boolean;
  Form1: TForm1;
  MutationRate : Real;
  ShirkVision,CoopVision,SelfishVision,DefectorVision,
    SelfishPunishVision,CooperatePunishVision,PunisherVision,
    FracInPool,AveScore,AveMinPun : TRAverage;
  TotalSR,TotalExecutions : Integer;
  SFrac,PFrac,CFrac : Real;
  AgentA,AgentB : Agent;

implementation

{$R *.DFM}

Function GetAgent(FracFromPool : Extended) : Integer; forward;
Function GetGroupAgent : Integer; forward;

Constructor Agent.Init;
begin
	inherited Create;
  AStatus := InPool;
  Cooperated := True;
end;

Destructor Agent.Destruct;
begin
  inherited;
end;

Procedure Agent.Copy(A : Agent);
begin
	ACoop := A.ACoop;
	APunish := A.APunish;
  AStatus := A.AStatus;
  Cooperated := A.Cooperated;
  Score := A.Score;
  SCost := A.SCost;
  MinPun := A.MinPun;
  GMember := A.GMember;
  Kill := False;
  NumCopies := 0;
end;

Function Agent.ASigma(FractionReciprocator: Real) : Real;
var
  R : Real;
begin
  R := 1 - FractionReciprocator*SCost/(2*C);
  if R < 0 then R := 0;
  ASigma := R;
end;

Constructor Group.Init;
begin
end;

Procedure Group.AddMember;
var
	I : Integer;
  R : Real;
  Around : Boolean;
begin
	Around := False;
	I := 1 + Random(NextAgent-1);
	while Agents[I].AStatus = InGroup do begin
  	Inc(I);
  	if I >= NextAgent then begin
    	if Around then begin
      	ShowMessage('Overflow in Group.Addmember!');
        Halt;
      end;
    	I := 1;
      Around := True;
    end;
  end;
  Agents[I].GMember := GNumber;
  Agents[I].AStatus := InGroup;
  Agents[I].APunish := Random < PFrac;
  R := Random;
  if R < CFrac then Agents[I].ACoop := Cooperator
  else if R < SFrac then Agents[I].ACoop := Selfish
  else Agents[I].ACoop := Defector;
end;

Procedure InitXYPopComposition;
const
	GraphHeight = 300;
  GraphWidth = 550;
begin
	with XYPopComposition do begin
  	Initialize(GraphWidth,GraphHeight,400,GraphHeight+20);
    Caption := 'Population Composition';
    ShowGraph;
    SetVLabel('');
    SetHLabel('Period');
    SetXAfterDecimalPoint(0);
    SetYAfterDecimalPoint(2);
    SetXYLimits(0,TotalPeriods,0,1);
    SetLegend('% Selfish');
    SetLegend('% Cooperator');
    SetLegend('% Defector');
    SetLegend('% Punisher');
    SetLegend('Shirking Level');
    SetLegend('% Pop in Pool');
  end;
end;

Procedure InitPopTypes;
const
	GraphHeight = 300;
  GraphWidth = 550;
begin
	with PopTypes do begin
  	Initialize(GraphWidth,GraphHeight,400,GraphHeight-300);
    Caption := 'Population Composition';
    ShowGraph;
    SetVLabel('');
    SetHLabel('Period');
    SetXAfterDecimalPoint(0);
    SetYAfterDecimalPoint(2);
    SetXYLimits(0,TotalPeriods,0,1);
    SetLegend('% Coop & Punish');
    SetLegend('% Selfish & Punish');
    SetLegend('% Cooperator MinPun');
  end;
end;

Procedure FillGroups;
var
	I,J: Integer;
begin
	for I := 1 to NumGroups do begin
    Groups[I].GNumber := I;
    for J := 1 to APerGroup do Groups[I].AddMember;
    Groups[I].GSize := APerGroup;
  end;
end;

Procedure CalcGSize;
var
	I: Integer;
begin
	for I := 1 to NumGroups do Groups[I].GSize := 0;
  for I := 1 to NextAgent-1 do
  	if Agents[I].AStatus = InGroup then
      Inc(Groups[Agents[I].GMember].GSize);
end;

Procedure CalculatePayoffs;
var
	I,G : Integer;
begin
	for I := 1 to NumGroups do begin
    Groups[I].Fs := 0;
    Groups[I].Fc := 0;
    Groups[I].Fd := 0;
    Groups[I].Fp := 0;
    Groups[I].LastGSigma := Groups[I].GSigma;
    Groups[I].GSigma := 0;
  end;
  CalcGSize;
  for I := 1 to NextAgent-1 do begin
    if Agents[I].AStatus = InGroup then begin
      G := Agents[I].GMember;
      case Agents[I].ACoop of
        Selfish : Groups[G].Fs := Groups[G].Fs + 1;
        Cooperator : Groups[G].Fc := Groups[G].Fc + 1;
        Defector : Groups[G].Fd := Groups[G].Fd + 1;
      end;
      if Agents[I].APunish then Groups[G].Fp := Groups[G].Fp + 1;
    end
  end;
  for I := 1 to NumGroups do begin
    if Groups[I].GSize > 0 then begin
      Groups[I].Fs := Groups[I].Fs/Groups[I].GSize;
      Groups[I].Fc := Groups[I].Fc/Groups[I].GSize;
      Groups[I].Fd := Groups[I].Fd/Groups[I].GSize;
      Groups[I].Fp := Groups[I].Fp/Groups[I].GSize;
    end;
    Groups[I].GSigma:=0;
  end;
  for I := 1 to NextAgent-1 do if Agents[I].AStatus = InGroup then begin
    G := Agents[I].GMember;
    case Agents[I].ACoop of
      Selfish : Groups[G].GSigma :=
        Groups[G].GSigma + Agents[I].ASigma(Groups[G].Fp);
      Defector : Groups[G].GSigma := Groups[G].GSigma + 1;
     Cooperator : if Groups[G].LastGSigma > Agents[I].MinPun then begin
        Groups[G].GSigma := Groups[G].GSigma + 1;
        Agents[I].Cooperated := False;
      end
      else Agents[I].Cooperated := True;
     end;
  end;
  for I := 1 to NumGroups do if Groups[I].GSize > 0 then
    Groups[I].GSigma := Groups[I].GSigma/Groups[I].GSize;
end;

Procedure Agent.Mutate;
begin
  case Random(3) of
    0 : ACoop := Defector;
    1 : ACoop := Selfish;
    2 : ACoop := Cooperator;
  end;
  APunish := Random(2) = 0;
  SCost := MinMutatedSCost + Random*(MaxMutatedSCost-MinMutatedSCost);
  if Random(2) = 0 then MinPun := MinPun*0.97
  else MinPun := MinPun/0.97;
end;

// Agents from groups who want to emigrate but are not admitted to
// another group remain where they started.
Procedure Immigrate;
var
	I,K,AInGroups,AInPool : Integer;
  FracFromPool, Admit : Extended;
begin
  AInGroups := 0;
  AInPool := 0;
  for I := 1 to NextAgent-1 do begin
    if (Agents[I].AStatus = InGroup) then Inc(AInGroups)
    else Inc(AInPool);
  end;
  if AInPool = 0 then FracFromPool := 0
  else FracFromPool := AInPool/(Gamma*AInGroups + AInPool);
  for I := 1 to NumGroups do begin
    Admit := Mu*AInGroups/NumGroups;
    while Random < Admit do begin
      K := GetAgent(FracFromPool);
      if K > 0 then begin
        Inc(Groups[I].GSize);
        Agents[K].AStatus := InGroup;
        Agents[K].GMember := I;
      end;
      Admit := Admit - 1;
    end;
  end;
end;

Function LargestGroup : Integer;
var
  I,Max,GMax : Integer;
begin
  Max := 0;
  GMax := 0;
  for I := 1 to NumGroups do begin
    if Groups[I].GSize > Max then begin
      Max := Groups[I].GSize;
      GMax := I;
    end;
  end;
  if Max > APerGroup then LargestGroup := GMax
  else LargestGroup := 0;
end;

Procedure RePopulateGroups;
var
	I,J,K,L,LL : Integer;
begin
	for I := 1 to NumGroups do begin
    CalcGSize;
  	if Groups[I].GSize < MinGroupSize then begin
      // Send rest of agents to pool
      for J := 1 to NextAgent-1 do
        if Agents[J].GMember = I then Agents[J].AStatus := InPool;
      Groups[I].GSize := 0;
      // Repopulation the group
      while Groups[I].GSize < APerGroup do begin
        // Recruit from largest group
      	J := LargestGroup;
        if J = 0 then Break; // No oversized groups left
        L := Groups[J].GSize - APerGroup;
        if L > APerGroup-Groups[I].GSize then L := APerGroup-Groups[I].GSize;
        LL := 0;
        K := 1;
        while LL < L do begin
          if (Agents[K].AStatus = InGroup) and (Agents[K].GMember = J) then begin
            Dec(Groups[J].GSize);
            Inc(Groups[I].GSize);
            Agents[K].GMember := I;
            Inc(LL);
          end;
          Inc(K);
          if K = NextAgent then Break;
        end;
      end;
      // Now recruit from pool. This will very rarely happen
      while Groups[I].GSize < APerGroup do begin
       	K := 1 + Random(NextAgent-1);
        while (Agents[K].AStatus = InGroup) and (K < NextAgent-1) do Inc(K);
        if K = NextAgent then begin
          K := 1;
          while (Agents[K].AStatus = InGroup) and (K < NextAgent) do Inc(K);
          if K = NextAgent then Break; // pool must be empty!
        end;
        Agents[K].AStatus := InGroup;
        Agents[K].GMember := I;
        Inc(Groups[I].GSize);
      end;
      // Finally, recruit from small groups. Almost never get here.
      while Groups[I].GSize < APerGroup do begin
      	K := 1 + Random(NextAgent-1);
        if Agents[K].AStatus = InGroup then Dec(Groups[Agents[K].GMember].GSize);
        Agents[K].AStatus := InGroup;
        Agents[K].GMember := I;
        Inc(Groups[I].GSize);
      end;
    end;
  end;
end;

Procedure Ostracize;
var
	I,G : Integer;
  FPunish : Real;
begin
	for I := 1 to NextAgent-1 do begin
    G := Agents[I].GMember;
    FPunish := Groups[G].Fp;
  	if Agents[I].AStatus = InGroup then
      case Agents[I].ACoop of
        Selfish : if Random < Agents[I].ASigma(FPunish)*FPunish then begin
            Dec(Groups[G].GSize);
            Agents[I].AStatus := InPool;
          end;
        Defector : if Random < FPunish then begin
          Dec(Groups[G].GSize);
          Agents[I].AStatus := InPool;
        end;
        Cooperator: if (not Agents[I].Cooperated) and (Random < FPunish) then begin
          Dec(Groups[G].GSize);
          Agents[I].AStatus := InPool;
        end;
      end;
  end;
end;

Procedure CullPopulation;
var
	I : Integer;
begin
  if (NextAgent < 10) or (NextAgent < TotalAgents div 10) then begin
    ShowMessage('There are now only '+IntToStr(NextAgent)+' agents!');
    Halt;
  end;
  while NextAgent >= TotalAgents+1 do begin
    I := 1 + Random(NextAgent-2);
    if I < NextAgent-1 then Agents[I].Copy(Agents[NextAgent-1]);
    Dec(NextAgent);
  end;
  CalcGSize;
end;

Procedure InitializeAgents;
var
	I : Integer;
begin
 	for I := 1 to NumAgents do begin
  	Agents[I].Score := 0;
    Agents[I].SCost := MinSCost + Random*(MaxSCost-MinSCost);
  	Agents[I].ACoop := Selfish;
    Agents[I].APunish := True;
    Agents[I].AStatus := InPool;
    Agents[I].Cooperated := True;
    Agents[I].MinPun := 0.40*Random;
	end;
end;

{The Score of an agent is its rate of growth, not including the possibility
of being ostracized, which is Fr*ASigma. After reproduction, this fraction
of selfish agents must be transferred to the pool.}
Procedure CalculateGlobalPayoffs;
var
	I,G : Integer;
  GSigma,Fs,Fp,ASigma : Real;
begin
	for I := 1 to NextAgent-1 do begin
    G := Agents[I].GMember;
    GSigma := Groups[G].GSigma;
    Fs := Groups[G].Fs;
    Fp := Groups[G].Fp;
    ASigma := Agents[I].ASigma(Fp);
		if Agents[I].AStatus = InPool then Agents[I].Score := Phi0
    else begin
      case Agents[I].ACoop of
        Cooperator: if not Agents[I].Cooperated then
          Agents[I].Score := B*(1-GSigma)
          else Agents[I].Score := B*(1-GSigma) - C;
        Selfish: Agents[I].Score := B*(1-GSigma) - C*(1-ASigma)*(1-ASigma);
        Defector : Agents[I].Score := B*(1-GSigma);
      end;
      if Agents[I].APunish then
        Agents[I].Score := Agents[I].Score - Cp*Fs*GSigma;
    end;
  end;
end;

Procedure Recombine(I : Integer);
var
  J,K : Integer;
begin
  J := Agents[I].Mate;
  if J = 0 then Exit;
  Agents[I].MakeCopy := False;
  Agents[J].MakeCopy := False;
  AgentA.SCost := Agents[I].SCost;
  AgentB.SCost := Agents[J].SCost;
  AgentA.AStatus := InGroup;
  AgentB.AStatus := InGroup;
  AgentA.MakeCopy := False;
  AgentB.MakeCopy := False;
  AgentA.Cooperated := True;
  AgentB.Cooperated := True;
  if Random(2) = 0 then AgentA.APunish := Agents[I].APunish
  else AgentA.APunish := Agents[J].APunish;
  if Random(2) = 0 then AgentB.APunish := Agents[I].APunish
  else AgentB.APunish := Agents[J].APunish;
  if Random(2) = 0 then AgentA.ACoop:= Agents[I].ACoop
  else AgentA.ACoop := Agents[J].ACoop;
  if Random(2) = 0 then AgentB.ACoop:= Agents[I].ACoop
  else AgentB.ACoop := Agents[J].ACoop;
  if AgentA.ACoop = Cooperator then begin
    if Random(2) = 0 then begin
      if Agents[I].ACoop = Cooperator then AgentA.MinPun:= Agents[I].MinPun
      else AgentA.MinPun := Agents[J].MinPun;
    end
    else begin
      if Agents[J].ACoop = Cooperator then AgentA.MinPun:= Agents[J].MinPun
      else AgentA.MinPun := Agents[I].MinPun;
    end;
  end;
  if AgentB.ACoop = Cooperator then begin
    if Random(2) = 0 then begin
      if Agents[I].ACoop = Cooperator then AgentB.MinPun:= Agents[I].MinPun
      else AgentB.MinPun := Agents[J].MinPun;
    end
    else begin
      if Agents[J].ACoop = Cooperator then AgentB.MinPun:= Agents[J].MinPun
      else AgentB.MinPun := Agents[I].MinPun;
    end;
  end;
  Agents[NextAgent].Copy(AgentA);
  if Random < MutationRate then Agents[NextAgent].Mutate;
  K := Agents[I].GMember;
  Agents[NextAgent].GMember := K;
  Inc(Groups[K].GSize);
  Inc(NextAgent);
  if NextAgent > NumAgents then begin
    ShowMessage('Population Exploded!');
    Halt;
  end;
  Agents[NextAgent].Copy(AgentB);
  if Random < MutationRate then Agents[NextAgent].Mutate;
  K := Agents[I].GMember;
  Inc(Groups[K].GSize);
  Agents[NextAgent].GMember := K;
  Inc(NextAgent);
  if NextAgent > NumAgents then begin
    ShowMessage('Population Exploded!');
    Halt;
  end;
end;

Procedure Reproduce;
var
	I,J,NumCopies : Integer;
begin
	for I := 1 to NextAgent-1 do begin
  	Agents[I].Kill := False;
    Agents[I].MakeCopy := False;
    Agents[I].Mate := 0;
  end;
  NumCopies := 0;
 	for I := 1 to NextAgent-1 do begin
  	if Agents[I].AStatus = InPool then begin
      if Phi0 < 0 then begin
        if Random < -Phi0 then Agents[I].Kill := True
      end
      else if Random < Phi0 then begin
        Agents[I].MakeCopy := True;
        Inc(NumCopies);
      end;
    end
    else if Agents[I].Score < 0 then begin
    	if Random < -Agents[I].Score then Agents[I].Kill := True;
    end
    else if Random < Agents[I].Score then begin
      Agents[I].MakeCopy := True;
      Inc(NumCopies);
    end;
  end;
  while NumCopies > 1 do begin
    I := 1+Random(NextAgent-1);
    repeat
      if Agents[I].MakeCopy and (Agents[I].Mate = 0) then Break;
      Inc(I);
      if I = NextAgent then I := 1;
    until False;
    J := 1+Random(NextAgent-1);
    repeat
      if (J <> I) and Agents[J].MakeCopy and (Agents[J].Mate = 0) then Break;
      Inc(J);
      if J = NextAgent then J := 1;
    until False;
    Agents[I].Mate := J;
    Agents[J].Mate := I;
    Dec(NumCopies,2);
  end;
  for I := 1 to NextAgent-1 do begin
		if Agents[I].Kill then begin
      if Agents[I].AStatus = InGroup then Dec(Groups[Agents[I].GMember].GSize);
      Dec(NextAgent);
  		Agents[I].Copy(Agents[NextAgent]);
    end;
  end;
  for I := 1 to NextAgent-1 do
    if Agents[I].MakeCopy then Recombine(I);
end;

Function AveSigma : Real;
var
	I,TotalJ : Integer;
  TotalR : Real;
begin
  TotalJ := 0;
  TotalR := 0;
	for I := 1 to NextAgent-1 do begin
    if Agents[I].AStatus = InGroup then begin
      Inc(TotalJ);
      case Agents[I].ACoop  of
        Selfish : TotalR := TotalR +
          Agents[I].ASigma(Groups[Agents[I].GMember].Fp);
        Defector : TotalR := TotalR + 1;
      end;
    end;
  end;
  if TotalJ > 0 then ShirkVision.Add(TotalR/TotalJ);
  AveSigma := ShirkVision.Ave;
end;

Procedure Averages;
var
  I,TotalSelfish,TotalCooperator,TotalDefector,SP,CP,
    TotalPunisher,SInPool,CInPool,DInPool : Integer;
  MinPunAve : Extended;
begin
  MinPunAve := 0;
  TotalSelfish := 0;
  TotalCooperator := 0;
  TotalDefector := 0;
  TotalPunisher := 0;
  TotalInPool := 0;
  SInPool := 0;
  CInPool := 0;
  DInPool := 0;
  CP := 0;
  SP := 0;
	for I := 1 to NextAgent-1 do begin
    case Agents[I].ACoop of
      Selfish : begin
        Inc(TotalSelfish);
        if Agents[I].APunish then Inc(SP);
      end;
      Defector : Inc(TotalDefector);
      Cooperator : begin
        Inc(TotalCooperator);
        if Agents[I].APunish then Inc(CP);
        MinPunAve := MinPunAve + Agents[I].MinPun;
      end;
    end;
    if Agents[I].APunish  then Inc(TotalPunisher);
    if Agents[I].AStatus = InPool then begin
      case Agents[I].ACoop of
        Selfish : Inc(SInPool);
        Defector : Inc(DInPool);
        Cooperator : Inc(CInPool);
      end;
    end;
  end;
  CoopVision.Add(TotalCooperator/NextAgent);
  SelfishVision.Add(TotalSelfish/NextAgent);
  DefectorVision.Add(TotalDefector/NextAgent);
  PunisherVision.Add(TotalPunisher/NextAgent);
  SelfishPunishVision.Add(SP/NextAgent);
  CooperatePunishVision.Add(CP/NextAgent);
  TotalInPool := DInPool + CInPool + SInPool;
  FracInPool.Add(TotalInPool/NextAgent);
  if TotalCooperator > 0 then AveMinPun.Add(MinPunAve/TotalCooperator);
end;

Procedure CalcStats;
var
  I : Integer;
begin
  Popsize := 0;
  RhoS := 0;
  RhoC := 0;
  RhoD := 0;
  RhoP := 0;
  for I := 1 to NumGroups do Inc(Popsize,Groups[I].GSize);
  for I := 1 to NextAgent-1 do begin
  	Case Agents[I].ACoop of
      Selfish : RhoS := RhoS+1;
      Cooperator : RhoC := RhoC+1;
  	  Defector : RhoD := RhoD+1;
    end;
    if Agents[I].APunish then RhoP := RhoP+1;
  end;
  RhoC := RhoC/(NextAgent-1);
  RhoD := RhoD/(NextAgent-1);
  RhoP := RhoP/(NextAgent-1);
  RhoS := RhoS/(NextAgent-1);
end;

// Return zero if the pool is empty
Function GetPoolAgent : Integer;
var
  I : Integer;
  Around : Boolean;
begin
  GetPoolAgent := 0;
  Around := False;
  I := 1+ Random(NextAgent-1);
  while Agents[I].AStatus <> InPool do begin
    Inc(I);
    if I = NextAgent then begin
      if Around then Exit;
      I := 1;
      Around := True;
    end;
  end;
  GetPoolAgent := I;
end;

Function GetAgent(FracFromPool : Extended) : Integer;
begin
  if Random < FracFromPool then GetAgent := GetPoolAgent
  else GetAgent := GetGroupAgent;
end;

Function GetGroupAgent : Integer;
var
  I : Integer;
  Around : Boolean;
begin
  GetGroupAgent := 0;
  Around := False;
  I := 1 + Random(NextAgent-1);
  while Agents[I].AStatus <> InGroup do begin
    Inc(I);
    if I = NextAgent then begin
      if Around then Exit;
      I := 1;
      Around := True;
    end;
  end;
  GetGroupAgent := I;
end;

Procedure GetTotalInPool;
var
  I : Integer;
begin
  TotalInPool := 0;
  for I := 1 to NextAgent-1 do
    if Agents[I].AStatus = InPool then Inc(TotalInPool);
end;

Procedure Play;
begin
  GetTotalInPool;
	CalculatePayoffs;
  CalculateGlobalPayoffs;
  Reproduce;
  Ostracize;
  Immigrate;
  CullPopulation;
  RePopulateGroups;
end;

Procedure SetSeedParameters;
var
  R : Real;
begin
  R := Rho0C+Rho0S+Rho0D;
  CFrac := Rho0C/R;
  SFrac := CFrac + Rho0S/R;
  PFrac := Rho0P;
end;

Procedure TForm1.SetParameters;
var
  I : Integer;
begin
  APerGroup := StrToInt(AgentsPerGroupBox.Text);
  NumGroups := StrToInt(NumGroupsBox.Text);
  if Length(Groups) < NumGroups+1 then begin
    SetLength(Groups,NumGroups+1);
    for I := 1 to NumGroups do
      if not Assigned(Groups[I]) then Groups[I] := Group.Init;
  end;
  TotalAgents := APerGroup*NumGroups;
  NextAgent := TotalAgents+1;
  NumAgents := Round(TotalAgents*1.5);
  if Length(Agents) < 1+NumAgents then begin
    SetLength(Agents,1+NumAgents);
    for I := 1 to NumAgents do
      if not Assigned(Agents[I]) then Agents[I] := Agent.Init;
  end;
  if not Assigned(AgentA) then AgentA := Agent.Init;
  if not Assigned(AgentB) then AgentB := Agent.Init;
  Rho0C := StrToFloat(Rho0CBox.Text);
  Rho0D := StrToFloat(Rho0DBox.Text);
  Rho0P := StrToFloat(Rho0PBox.Text);
  Rho0S := StrToFloat(Rho0SBox.Text);
  SetSeedParameters;
  C := StrToFloat(CBox.Text);
  Cp := StrToFloat(CPBox.Text);
  B := StrToFloat(BBox.Text);
  Gamma := StrToFloat(GammaBox.Text);
  Mu := StrToFloat(MuBox.Text);
  Phi0 := StrToFloat(Phi0Box.Text);
  MutationRate := StrToFloat(MutationRateBox.Text);
  MinGroupSize := StrToInt(MinGroupSizeBox.Text);
  TotalPeriods := Round(StrToFloat(TotalPeriodsBox.Text));
  if TotalPeriods < PeriodsPerPlot then begin
    ShowMessage('Total Periods must be greater than '+IntToStr(PeriodsPerPlot));
    Halt;
  end;
  MaxSCost := StrToFloat(MaxSSeedBox.Text);
  MinSCost := StrToFloat(MinSSeedBox.Text);
  MaxMutatedSCost := StrToFloat(MaxMutatedSBox.Text);
  MinMutatedSCost := StrToFloat(MinMutatedSBox.Text);
end;

Procedure CalcScore;
var
  I : Integer;
  TScore : Real;
begin
  TScore := 0;
  for I := 1 to NextAgent-1 do
      TScore := TScore + Agents[I].Score;
  AveScore.Add(TScore/NextAgent);
end;

Procedure TForm1.RunBtnClick(Sender: TObject);
var
	I : Integer;
  Stop : Boolean;
begin
  Stop := False;
  ShirkVision := TRAverage.Init(Vision);
  CoopVision := TRAverage.Init(Vision);
  SelfishVision := TRAverage.Init(Vision);
  DefectorVision := TRAverage.Init(Vision);
  PunisherVision := TRAverage.Init(Vision);
  SelfishPunishVision := TRAverage.Init(Vision);
  CooperatePunishVision := TRAverage.Init(Vision);
  AveScore := TRAverage.Init(Vision);
  AveMinPun:= TRAverage.Init(Vision);
  FracInPool := TRAverage.Init(Vision);
  SetParameters;
  if not Assigned(XYPopComposition) then
    XYPopComposition := TXYGraph.Create(Self);
  InitXYPopComposition;
  if not Assigned(PopTypes) then
    PopTypes := TXYGraph.Create(Self);
  InitPopTypes;
  InitializeAgents;
  FillGroups;
	TotalCount := 0;
  Periods := 0;
  HaltSimulation := False;
  HaltBtn.SetFocus;
	for I := 1 to TotalPeriods do begin
    if Stop then Break;
    Inc(Periods);
		Play;
    CalcStats;
    if Periods mod PeriodsPerPlot = 0 then begin
      Application.ProcessMessages;
			CurrentPeriodBox.Text := Format('%d',[Periods]);
      CurrentPeriodBox.Update;
      Application.ProcessMessages;
      if HaltSimulation then Break;
      Averages;
      with XYPopComposition do begin
      	Addpoint(1,Periods,SelfishVision.Ave);
      	Addpoint(2,Periods,CoopVision.Ave);
      	Addpoint(3,Periods,DefectorVision.Ave);
      	Addpoint(4,Periods,PunisherVision.Ave);
        Addpoint(5,Periods,AveSigma); 
        AddPoint(6,Periods,FracInPool.Ave);
      end;
      with PopTypes do begin
      	Addpoint(1,Periods,CooperatePunishVision.Ave);
      	Addpoint(2,Periods,SelfishPunishVision.Ave);
      	Addpoint(3,Periods,AveMinPun.Ave);
      end;
		end;
	end;
  RunBtn.SetFocus;
end;

Procedure TForm1.HaltBtnClick(Sender: TObject);
begin
	HaltSimulation := True;
end;

procedure TForm1.FormClose(Sender: TObject; var Action: TCloseAction);
begin
	HaltSimulation := True;
  Application.Terminate;
end;

procedure TForm1.QuitBtnClick(Sender: TObject);
begin
  Halt;
end;

procedure TForm1.HaltBtnMouseDown(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
begin
  HaltBtnClick(Sender);
end;

procedure TForm1.QuitBtnMouseDown(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
begin
  HaltBtnClick(Sender);
  Halt;
end;

procedure TForm1.HelpButtonClick(Sender: TObject);
begin
  HelpForm.ShowModal;
end;

procedure TForm1.ExecuteButtonClick(Sender: TObject);
begin
  Execute := True;
  TotalExecutions := StrToInt(ExecutionCycleBox.Text);
  TotalSR := 0;
  ExecutionCycle := 1;
  while ExecutionCycle <= TotalExecutions do begin
    RunBtnClick(Sender);
    ExecutionCycleBox.Text := IntToStr(TotalExecutions-ExecutionCycle);
    ExecutionCycleBox.Update;
    Application.ProcessMessages;
    if HaltSimulation then begin
      Execute := False;
      ExecutionCycleBox.Text := IntToStr(TotalExecutions);
      ExecutionCycleBox.Update;
      Exit;
    end;
    Inc(ExecutionCycle);
  end;
  Execute := False;
  ExecutionCycleBox.Text := IntToStr(TotalExecutions);
  ExecutionCycleBox.Update;
end;

Function GetFloat(var S : String) : Real;
var
  T : String;
begin
  T := '';
  while (S<>'') and (S[1] in ['0'..'9','.',' ']) do begin
    if S[1]<>' ' then T := T + S[1];
    Delete(S,1,1);
  end;
  if S <> '' then Delete(S,1,1);
  GetFloat := StrToFloat(T);
end;

initialization
  XYPopComposition := Nil;
  Execute := False;
end.