unit TCProg;

interface

uses TCG;

const
  MaxOps = 10;
  MaxTrace = 60;

type
  Ops = array [0..MaxOps-1] of TCG.Gate;

  ProgP = ^ Prog;
  Prog = record
    wid: integer;  { Width of gates in bits. }
    len: integer;  { Number of ops in op[]. }
    op: Ops;
  end;

  Trace = string[MaxTrace];

function mismatch(
    var x, y: Prog;   { Programs to compare. }
    ix, iy: integer;  { Positions to compare. }
    wd: integer       { Window width. }
  ): integer;
  { Compares x[ix+i] with y[iy+i] (modulo MaxOps),
    for -wd <= i <= wd. Returns weighted difference. }

function best_y_match(
    var x, y: Prog;  { Programs to compare }
    jx, jy: integer; { Positions to compare. }
    maxoff: integer; { Offsets to compare. }
    wd: integer      { Comparison window radius. }
  ): integer;
  { Returns the index ky = (jy + offy) mod py.len such that
    -maxoff <= offy <= maxoff and y.op[ky] best matches
    x.op[jx] in a window of width wd.
  }

procedure mate_and_mutate(px, py: ProgP; var z: Prog; var trc: Trace);
  { Puts in z a combination of px^ and py^ with mutations. }

implementation

uses Crt;

function mismatch(
    var x, y: Prog;   { Programs to compare. }
    ix, iy: integer;  { Positions to compare. }
    wd: integer       { Window width. }
  ): integer;
  var d, w, sum: integer;

  function diff1(kx, ky: integer): integer;
    var xg, yg: TCG.Gate;
    begin
      if (kx >= 0) and (kx < x.len) then
        xg := x.op[kx]
      else
        TCG.clear(xg);
      if (ky >= 0) and (ky < y.len) then
        yg := y.op[ky]
      else
        TCG.clear(yg);
      diff1 := TCG.diff(xg, yg)
    end;

  begin
    w := 1 shl wd;
    sum := w * diff1(ix, iy);
    { Compare center positions: }
    for d := 1 to wd do
      begin
        w := w shr 1;
        sum := sum + w * diff1(ix + d, iy + d);
        sum := sum + w * diff1(ix - d, iy - d);
      end;
    mismatch := sum
  end;

function best_y_match(
    var x, y: Prog;  { Programs to compare }
    jx, jy: integer; { Positions to compare. }
    maxoff: integer; { Offsets to compare. }
    wd: integer      { Comparison window radius. }
  ): integer;
  var aoff: integer;
      adiff, best_diff: integer;
      ix, iy: integer;
  begin
    if (y.len = 0) then
      best_y_match := 0
    else
      begin
        ix := jx;
        iy := jy mod y.len;
        best_y_match := iy;
        best_diff := mismatch(x, y, ix, iy, wd);
        for aoff := 1 to maxoff do
          begin
            iy := (jy + aoff) mod y.len;
            adiff := mismatch(x, y, ix, iy, wd);
            if adiff < best_diff then
              begin best_diff := adiff; best_y_match := iy end;
            iy := (jy + y.len - aoff) mod y.len;
            if iy < 0 then iy := iy + y.len;
            adiff := mismatch(x, y, ix, iy, wd);
            if adiff < best_diff then
              begin best_diff := adiff; best_y_match := iy end;
          end;
      end;
  end;


procedure mate_and_mutate(px, py: ProgP; var z: Prog; var trc: Trace);
  var jx, jy, jz: integer;
      nit, i, nmut, nzest: integer;
      mutp: real;

  function toss(p: real): boolean;
    begin
      toss := (Random < p)
    end;

  procedure ptrace(c: char);
    var ln: integer;
    begin
      ln := ord(trc[0]);
      if ln < MaxTrace then
        begin
          inc(ln);
          trc[ln] := c;
          trc[0] := chr(ln)
        end
    end;

  procedure copyx;
    begin
      if (jx < px^.len) and (jz < MaxOps) then
        begin
          z.op[jz] := px^.op[jx];
          inc(jz);
          inc(jx);
          if (jy < py^.len) then inc(jy);
          { Let's not count it as a mutation for now: }
          ptrace('x');
        end;
    end;

  procedure copyxy;
    begin
      if (jx < px^.len) and (jy < py^.len)
      and (jz < MaxOps) then
        begin
          TCG.merge(px^.op[jx], py^.op[jy], z.op[jz]);
          if TCG.eq(px^.op[jx], py^.op[jy]) then
            ptrace('=')
          else if TCG.eq(z.op[jz], px^.op[jx]) then
            ptrace('X')
          else if TCG.eq(z.op[jz], py^.op[jy]) then
            ptrace('Y')
          else
            begin ptrace('#'); inc(nmut) end;
          inc(jz);
          inc(jx);
          inc(jy);
        end;
    end;

  procedure skipx;
    begin
      if (jx < px^.len) then
        begin
          inc(jx);
          inc(nmut);
          ptrace('.')
        end
    end;

  procedure insertz;
    var n1, n2, nm, nv: integer;
    begin
      if (jz < MaxOps) then
        begin
          n1 := random(z.wid);
          n2 := random(z.wid);
          if n1 > n2 then nm := n1 else nm := n2;
          n1 := 1 + random(z.wid - nm);
          n2 := 1 + random(z.wid - nm);
          if n1 < n2 then nv := n1 else nv := n2;
          TCG.toss(z.wid, nm, nv, z.op[jz]);
          inc(jz);
          inc(nmut);
          ptrace('+')
        end
    end;

  procedure deletez;
    begin
      if (jz > 0) then
        begin
          dec(jz);
          inc(nmut);
          ptrace('-')
        end
    end;

  procedure swapz;
    var g: TCG.Gate;
    begin
      if (jz >= 2) then
        begin
          g := z.op[jz-1];
          z.op[jz-1] := z.op[jz-2];
          z.op[jz-2] := g;
          inc(nmut);
          ptrace('~')
        end
    end;

  procedure mutatez;
    begin
      if (jz > 0) then
        begin
          TCG.mutate(z.wid, z.op[jz-1]);
          inc(nmut);
          ptrace('!')
        end
    end;

  procedure crunchz;
    begin
      if (jz > 1) then
        begin
          TCG.merge(z.op[jz], z.op[jz-1], z.op[jz-1]);
          dec(jz);
          inc(nmut);
          ptrace('c')
        end
    end;

  procedure switchxy;
    var jt: integer;
        pt: ProgP;
    begin
      pt := px; px := py; py := pt;
      jt := jx; jx := jy; jy := jt;
      { not necessarily a mutation...}
      ptrace('*');
    end;

  procedure syncy(maxoff: integer);
    var jynew: integer;
    begin
      if (jx < px^.len) then
        begin
          jynew := best_y_match(px^, py^, jx, jy, maxoff, 2);
          if jynew < jy then ptrace('<');
          if jynew > jy then ptrace('>');
          { Not necessarily a mutation; may undo them... }
          jy := jynew
        end
    end;

  procedure jumpx;
    var dx: integer;
    begin
      if px^.len > 0 then
        begin
          dx := 1;
          while toss(0.50) do
            begin
              dx := dx shl 1;
              if toss(0.50) then dx := dx or 1
            end;
          if toss(0.50) then dx := -dx;
          jx := (jx + dx) mod px^.len;
          if jx < 0 then jx := jx + px^.len;
          syncy(py^.len div 2);
          inc(nmut);
          ptrace('j')
        end
    end;

  procedure sort_and_simplify;
    var i, j, nsimp: integer;
        g1, g2: TCG.Gate;
        e1, e2: boolean;
        ok: array [0..MaxOps] of boolean;
          { True if we know that (z[i-1], z[i]) is simplified. }
    begin
      if jz > 0 then
        begin
          ok[0] := true;
          for i := 1 to jz-1 do ok[i] := false;
          ok[jz] := true;
          TCG.simp1(z.op[0]);
          i := 1;
          nsimp := 0;
          while (i < jz) do
            begin
              { At this point, ok[1..i-1] is all true. }
              if ok[i] then
                inc(i)
              else
                begin
                  g1 := z.op[i-1];
                  g2 := z.op[i];
                  TCG.simp2(z.op[i-1], z.op[i]);
                  e1 := TCG.eq(g1, z.op[i-1]);
                  e2 := TCG.eq(g2, z.op[i]);
                  if not (e1 and e2) then
                    begin
                      inc(nsimp);
                      if TCG.eq(g1, z.op[i]) and TCG.eq(g2, z.op[i-1]) then
                        ptrace('$')
                      else
                        ptrace('%')
                    end;
                  if (i >= 2) and not e1 then
                    ok[i-1] := false;
                  ok[i] := true;
                  if (i <= jz-2) and not e2 then
                    ok[i+1] := false;
                  if not ok[i-1] then dec(i) else inc(i)
                end;
            end;
          { Eliminate no-op gates: }
          i := 0; j := 0;
          while j < jz do
            begin
              if (z.op[j].val and not z.op[j].msk) <> 0 then
                begin z.op[i] := z.op[j]; inc(i) end
              else
                ptrace('0');
              inc(j)
            end;
          jz := i;
        end;
    end;

  begin
    jx := 0;
    jy := 0;
    jz := 0;
    z.len := 0;
    if px^.wid > py^.wid then z.wid := px^.wid else z.wid := py^.wid;
    nit := 0;
    nmut := 0;
    trc := '';
    repeat
      { Estimate number of gates still to be copied: }
      nzest := (px^.len - jx + py^.len - jy) div 2;
      if nzest < 1 then nzest := 1;
      mutp := (1.0 - nmut)/nzest;
      if toss(mutp) then
        begin
          { Mutate the copy just made: }
          case random(100) of
            00..59:  mutatez;
            60..69:  swapz;
            70..79:  if (jz + nzest) < MaxOps then insertz else mutatez;
            80..89:  if toss((1.0 + jz + nzest) / MaxOps)
                       then deletez else insertz;
            90..92:  crunchz;
            93..95:  jumpx;
            96..99:  switchxy;
          end;
        end
      else if (jx < px^.len) then
        begin
          { Copy one strand, or both: }
          syncy(5);
          if (jy < py^.len)
          and TCG.eq(px^.op[jx], py^.op[jy]) then
            begin
              copyxy;
              if toss(0.20) then switchxy;
            end
          else
            begin
              copyx
            end;
        end;
      { Count iterations: }
      inc(nit);
    until (nit > 10*MaxOps)
    or (((jx = px^.len) or (jz >= MaxOps)) and toss(1.00 - mutp));
    sort_and_simplify;
    z.len := jz;
    for i := length(trc) + 1 to MaxTrace do
      ptrace(' ');
  end;

end.