program sga (input, output, param, out);

{ A Simple Genetic Algorithm - SGA }

const maxpop    = 100;
      maxstring = 30;

type allele     = boolean;
     chromosome = array[1..maxstring] of allele;
     individual = record
                   chrom: chromosome;
                   fitness: real;
                   parent1, parent2, xsite: integer;
                  end;
     population = array [1..maxpop] of individual;

var oldpop, newpop: population;
    popsize, lchrom, gen, maxgen: integer;
    pcross, pmutation, sumfitness, randomseed: real;
    nmutation, ncross: integer;
    avg, max, min: real;
    param, out: text;
    linelength: integer;


{in interfac.p}
function objfunc (chrom: chromosome; lbits: integer): real; external;

{in utility.p}
procedure page (var out: text);
{ Issue form feed to device of file }
begin writeln (out); writeln (out); end;

procedure repchar (var out: text; ch: char; repcount: integer);
{ repeatedly write a character to an output device }
var j: integer;
begin for j:= 1 to repcount do write (out, ch) end;

procedure skip (var out: text; skipcount: integer);
{ Skip skipcount lines on dvice out }
var j: integer;
begin for j := 1 to skipcount do writeln (out) end;

function power (x, y: real): real;
{ Raise x to the yth power }
begin power := exp ( y * ln(x) ) end;

{in random.p}
{ random.p: contains random number generator and related utilities
            including advance_arndom, warmup_random, random, randomize,
            flip, rnd }

{ Global variables - Don't use these names in other code}
var oldrand: array[1..55] of real; { Array of 55 random numbers }
    jrand: integer;                { current random }

procedure advance_random;
{ Create next batch of 55 random numbers }
var j1: integer;
    new_random: real;
begin
 for j1 := 1 to 24 do
  begin
   new_random := oldrand[j1] - oldrand[j1+31];
   if (new_random < 0.0) then new_random := new_random + 1.0;
   oldrand[j1] := new_random;
  end;
 for j1 := 25 to 55 do
  begin
   new_random := oldrand[j1] - oldrand[j1-24];
   if (new_random < 0.0) then new_random := new_random + 1.0;
   oldrand [j1] := new_random;
  end;
end; {advance_random}

procedure warmup_random (random_seed: real);
{ Get random off and runnin }
var j1, ii: integer;
    new_random, prev_random: real;
begin
 oldrand[55] := random_seed;
 new_random := 1.0e-9;
 prev_random := random_seed;
 for j1 := 1 to 54 do
  begin
   ii := 21*j1 mod 55;
   oldrand[ii] := new_random;
   new_random := prev_random - new_random;
   if (new_random < 0.0) then new_random := new_random + 1.0;
   prev_random := oldrand[ii];
  end;
  advance_random; advance_random; advance_random;
  jrand := 0;
end; { warmup_random }

function random: real;
{ Fetch a single random random number between 0.0 and 1.0 - Substractive }
begin
 jrand := jrand + 1;
 if (jrand > 55) then
  begin jrand := 1; advance_random end;
 random := oldrand[jrand];
end; { random }

function flip (probability: real): boolean;
{ Flip a biased coin - true if heads }
begin
 if probability = 1.0 then flip := true
    else flip := (random <= probability);
end; { flip }

function rnd (low, high: integer): integer;
{ Pick a random integer between low and high }
var i: integer;
begin
 if low >= high then i := low
    else begin
     i := trunc( random * (high - low + 1) + low );
     if i > high then i := high;
    end;
 rnd := i;
end; { rnd }


procedure statistics (popsize: integer;
                      var max, avg, min, sumfitness: real;
                      var pop: population);
{ Calculate population statistics }
var j: integer;
begin
 { Initialize }
 sumfitness := pop[1].fitness;
 min        := pop[1].fitness;
 max        := pop[1].fitness;
 { Loop for max, min, sumfitness }
 for j := 2 to popsize do with pop[j] do begin
  sumfitness := sumfitness + fitness;   { Accumulate fitness sum }
  if fitness > max then max := fitness; { New max }
  if fitness < min then min := fitness; { New min }
 end;
 { Calculate average }
 avg := sumfitness / popsize;
end; { statistics }

{in initial.p}
procedure initdata;
{ Interactive data inquiry and set up }
{var ch: char; j: integer;}
begin
 rewrite(out); { Set up for list device }
 write (out, chr(15));
{ clrscr;}       { Clear screen }
 writeln;
 repchar (output, ' ', 25); writeln ('---------------------------------');
 repchar (output, ' ', 25); writeln ('A Simple Genetic Algorithm - SGA ');
 repchar (output, ' ', 25); writeln ('---------------------------------');
 writeln ('* Initialization of SGA parameters from "param" ...');
 reset(param);
 readln (param, popsize); { population size }
 readln (param, lchrom); { cromosome length }
 readln (param, maxgen); { max generations }
 readln (param, pcross); { cpossover probability }
 readln (param, pmutation); { mutation probalility }
 readln (param, randomseed); { random seed }
 close (param);
 writeln ('is complete.');
 { Initialize random number generator }
 warmup_random (randomseed);
 { Initilaize counters }
 nmutation := 0;
 ncross := 0;
 linelength := 2*lchrom + 50;
end; { initdata }

procedure initreport;
{ Initial report }
begin
 writeln (out, '-------------------------------------');
 writeln (out, '   A Simple Genetic Algorithm - SGA  ');
 writeln (out, '-------------------------------------');
 skip (out, 5);
 writeln (out, '   SGA parameters ');
 writeln (out, '   -------------- ');
 writeln (out);
 writeln (out, ' Population size (popsize)        = ', popsize);
 writeln (out, ' Chromosome length (lchrom)       = ', lchrom);
 writeln (out, ' Maximum # of generation (maxgen) = ', maxgen);
 writeln (out, ' Crossover probability (pcross)   = ', pcross);
 writeln (out, ' Mutation probability (pmutation) = ', pmutation);
 skip (out, 8);
 writeln (out, ' Initial population maximum fitness = ', max);
 writeln (out, ' Initial population average fitness = ', avg);
 writeln (out, ' Initial population minimum fitness = ', min);
 writeln (out, ' Initial population sum of fitness  = ', sumfitness);
 page (out); { New page }
end; { initreport }

procedure initpop;
{ Initialize a population at random }
var j, j1: integer;
begin
 for j := 1 to popsize do with oldpop[j] do begin
  for j1 := 1 to lchrom do chrom[j1] := flip(0.5); { A fair coin toss }
  fitness := objfunc (chrom, lchrom); { Evaluate initial fitness }
  parent1 := 0; parent2 := 0; xsite := 0; { Initialize printout vars }
 end;
end; { initpop }

procedure initialize;
{ Initialization Coordinator }
begin
 initdata;
 initpop;
 statistics (popsize, max, avg, min, sumfitness, oldpop);;
 initreport;
end; { initialize }

{in triops.p}
{ 3-operators: Reproduction (select), Crossover (crossover),
               & Mutation (mutation)                        }

function select (popsize: integer; sumfitness: real;
                 var pop: population): integer;
{ Select a single individual via roulette wheel selection }
var rand, partsum: real; { random point on wheel, selection }
    j: integer;          { population index }
begin
 partsum := 0.0; j := 0; { Zero out counter and accumulator }
 rand := random * sumfitness; { Wheel point calc. uses random number [0,1]}
 repeat { Find wheel slot }
  j := j + 1;
  partsum := partsum + pop[j].fitness;
 until (partsum >= rand) or (j = popsize);
 { Return individual number }
 select := j;
end; { select }

function mutation (alleleval: allele; pmutation: real;
                   var nmutation: integer): allele;
{ Mutate an allele w/ pmutation, count number of mutations }
var mutate: boolean;
begin
 mutate := flip(pmutation); {flip the biased coin }
 if mutate then begin
  nmutation := nmutation + 1;
  mutation  := not alleleval;  { change bit value }
 end else
  mutation := alleleval;       { no change }
end; { mutation }

procedure crossover (var parent1, parent2, child1, child2: chromosome;
                     var lchrom, ncross, nmutation, jcross: integer;
                     var pcross, pmutation: real);
{ Cross 2 parent strings, place in 2 child strings }
var j: integer;
begin
 if flip (pcross) then begin   { Do crossover with p(cross) }
  jcross := rnd (1, lchrom-1); { Cross between 1 and l-1}
  ncross := ncross + 1;        { Increment crossover counter }
 end else { Otherwise set cross site to force mutation }
 jcross := lchrom;
 { 1st exchange 1 to 1 and 2 to 2 }
 for j := 1 to jcross do begin
  child1[j] := mutation (parent1[j], pmutation, nmutation);
  child2[j] := mutation (parent2[j], pmutation, nmutation);
 end;
 { 2nd exchange 1 to 2 and 2 to 1 }
 if jcross <> lchrom then { Skip if cross site is lchrom -- no crossover }
  for j := jcross + 1 to lchrom do begin
   child1[j] := mutation(parent2[j], pmutation, nmutation);
   child2[j] := mutation(parent1[j], pmutation, nmutation);
  end;
end; { crossoever }

{in generate.p}
procedure generation;
{ Create a new generation through select, crossover, and mutation }
{ Note: generation assumes and even-numbered popsize              }
var j, mate1, mate2, jcross: integer;
begin
 j := 1;
 repeat  { select, crossover, and mutation until newpop is filled }
  mate1 := select (popsize, sumfitness, oldpop); {pick pair of mates}
  mate2 := select (popsize, sumfitness, oldpop);
  { Crossover and mutation - mutation embedded within crossover }
  crossover (oldpop[mate1].chrom, oldpop[mate2].chrom,
             newpop[j    ].chrom, newpop[j + 1].chrom,
             lchrom, ncross, nmutation, jcross, pcross, pmutation);
  { Decode string, evaluate fitness, record parantage date on both children}
  with newpop[j] do begin
   fitness := objfunc (chrom, lchrom);
   parent1 := mate1;
   parent2 := mate2;
   xsite   := jcross;
  end;
  with newpop[j+1] do begin
   fitness := objfunc (chrom, lchrom);
   parent1 := mate1;
   parent2 := mate2;
   xsite   := jcross;
  end;
  { Increment population index }
  j := j + 2;
 until j > popsize
end; { generate.gsa }

procedure writechrom (var out: text; chrom: chromosome; lchrom: integer);
{ Write a chromosome as a string of 1's (true's) and 0's (false's) }
var j: integer;
begin
 for j := lchrom downto 1 do
  if chrom[j] then write (out, '1')
     else write (out, '0');
end; { writechrom }

procedure report (gen: integer);
{ Write the population report }
var j: integer;
begin
 repchar (out, '-', linelength); writeln(out);
 repchar (out, ' ', 50); writeln (out, 'Population Report');
 repchar (out, ' ', 23); write (out, 'Generation ', gen-1: 2);
 repchar (out, ' ', 57); writeln (out, 'Generation ', gen: 2);
 writeln (out);
 write (out, '  #   ', 'string':round(lchrom/2)+3, ' ':round(lchrom/2)-3,
	' fitness ');
 write (out, '        #  parents  xsite');
 writeln (out, '             string              fitness');
 repchar (out, '-', linelength); writeln (out);
 for j := 1 to popsize do begin
  write (out, j:3, ')  ');
  { Old string }
  with oldpop[j] do begin
   writechrom (out, chrom, lchrom);
   write (out, ' ', ' ', fitness:6:4, '    |');
  end;
  { New string }
  with newpop[j] do begin
   write (out, '   ', j:2, ') (', parent1:2, ',', parent2:2, ')   ',
               xsite:2, '   ');
   writechrom (out, chrom, lchrom);
   writeln (out, ' ', ' ', fitness:6:4);
  end;
 end;
 repchar (out, '-', linelength); writeln (out);
 { Generation statistics and accumulated values }
 writeln (out, 'Gen: ', gen:2, '  Accumulated Statistics: '
             , ' max=', max:6:4, ', min=', min:6:4, ', avg=', avg:6:4
             , ', sum=', sumfitness:6:4, 
	     ', nmutation=', nmutation:3
             , ', ncross=', ncross:3);
 repchar (out, '-', linelength); writeln(out);
 page(out);
end; { report }

begin
 gen := 0;
 initialize;
 repeat
  gen := gen + 1;
  generation;
  statistics (popsize, max, avg, min, sumfitness, newpop);
  report (gen);
  oldpop := newpop;
  until (gen >= maxgen)
end. {sga}