program patlrn (funcbook,funcinst,nfuncbook,nfuncinst,pattern,start,
                 minmax,ignore,patlrnp,output);
(* patlrn - pattern learning program
      by gary stormo (modified by tom schneider)

  Thomas D. Schneider, Ph.D.
  toms@alum.mit.edu
  https://alum.mit.edu/www/toms 

modules needed: delman, delmods, auxmods  *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 3.32; (* of patlrn.p 2016 Jan 24
2016 Jan 24, 3.32: remove nonstandard pascal:
       write(output, ', pointer id: ',l:1);
1999 Dec 13, 3.31: y2k upgrade
1998 July 9: upgrade to new alist
1995 Aug 30: probably conversion to unix
origin july 9, 1982 *)
(* end module version *)

(* begin module describe.patlrn *)
(*
name
      patlrn: pattern learning

synopsis
      patlrn(funcbook: in, funcinst: in, nfuncbook: in, nfuncinst: in,
            pattern: out, start: in, minmax: in, ignore: in, patlrnp: in,
            output: out)

files
      funcbook: the book of sequences belonging to the functional class;
      funcinst: the instructions for funcbook, for aligning the sequences;
      nfuncbook: the book of sequences for the nonfunctional class;
      nfuncinst: the instructions for nfuncbook, for aligning the seqs;
      pattern: the resulting wmatrix which separates the classes;
      start: a matrix for initializing wmatrix to. it is initialized to
            all 0's if this file is empty;
      minmax: to set the values of funcmin (the minimum value for a functional
            sequence) and nfuncmax (the maximum value for a nonfunctional
            sequence). if this file is empty they are set to 1 and 0,
            respectively, and vary along with the matrix;
      ignore: a file specifying regions of the sequences which are to be
            ignored in the learning process; the maximum number of regions
            which can be ignored is set by the constant 'maxignore';  the
            file must contain two integers per line, the first specifying
            the 5' end and the second the 3' end of the region to be ignored.
      patlrnp: parameter file for setting maxtimes, the number of times
            through all the sequences before stopping without a solution;
      output: for messages to the user.

description
      patlrn uses the 'perceptron' algorithm to find a weighting function
      (a 'wmatrix') which serves to distinguish the sequences in the two
      classes from one another.  our paper, stormo et.al., nar 10, 2995 (1982),
      describes the algorithm in detail and gives an example of its use.

see also
      patlst.p, patana.p, patser.p, patval.p, sepa.p

author
      gary d. stormo (modified by tom schneider)

bugs
      the section of code for ignoring regions of the sequences in the
      learning process (i.e., when the file 'ignore' is not empty) has
      been overlayed over the rest of the code, rather than worked into
      it, and consequently, using this feature can be quite inefficient.

technical note
      the program will be more efficient if the constant 'dnamax' in the
      module 'book.const' is made to be the size of the sequences used
      by the program.  for instance, setting it to whatever 'maxmatrix'
      is would be a good idea.
*)
(* end module describe.patlrn *)


(* LOCK begin module book.const *)
(* constants needed for book manipulations *)

(* THIS MODULE IS LOCKED SO THAT THE DNAMAX IS RESTRICTED FOR CONSERVATION OF
MEMORY *)

      dnamax = 100;  (* length of dna arrays *)
      namelength = 100; (* maximum key name length *)
      linelength = 80; (* maximum line readable in book *)

(* LOCK end module book.const version = 'delmod 6.87 98 Feb 27 tds/gds' *)

(* begin module datetime.const *)
   datetimearraylength = 19; (* length of dataarray for dates,
       It is just long enough to include the 4 digit year - solving the
       year 2000 problem:
1980/06/09 18:49:11
123456789 123456789 
         1         2
*)
(* end module datetime.const version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module matrix.const *)
      maxmatrix = 120; (* the maximum size of a pattern matrix *)
(* end module matrix.const version = 'auxmod 1.40 94 Sep 5 gds/tds'; *)
      maxignore = 10; (* maximum number of regions to ignore *)

(* begin module interact.const *)
(* begin module string.const *)
      maxstring = 150; (* the maximum string *)
(* end module string.const version = 4.39; (@ of prgmod.p 1999 November 28 *)
(* end module interact.const version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module filler.const *)
      fillermax = 50; (* the size of the filler array for a string *)
(* end module filler.const version = 4.21; (@ of prgmod.p 1997 October 22 *)

type

(* begin module interact.type *)
(* begin module string.type *)
   stringptr = ^string; (* pointer to a string *)
   string = record (* a string of characters *)
      letters: array[1..maxstring] of char; (* the letters in the string *)
      length: integer; (* the number of characters in the string *)
      current: integer; (* the letter we are working on *)
      next: stringptr; (* the next string in a series *)
   end;
(* end module string.type version = 4.39; (@ of prgmod.p 1999 November 28 *)
(* end module interact.type version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module filler.type *)
      (* the following is an array used to fill a string.
      it is convenient to have it much shorter than the maxstring, so that
      it is easy to fill the string using procedure fillstring.
      the user must declare the value of constant fillermax. *)
      filler = packed array[1..fillermax] of char;
(* end module filler.type version = 4.21; (@ of prgmod.p 1997 October 22 *)

(* begin module trigger.type *)
      trigger = record (* an object to be searched for *)
         seek: string; (* the characters looked for *)
         state: integer; (* how close to triggering we are *)
         skip: boolean; (* trigger not found- skip the line *)
         found: boolean (* the trigger was found *)
      end;
(* end module trigger.type version = 4.21; (@ of prgmod.p 1997 October 22 *)

(* begin module datetime.type *)
   (* array for dates *)
   datetimearray = packed array[1..datetimearraylength] of char;
(* end module datetime.type version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.type *)
(* types needed for book manipulations *)

      chset = set of 'a'..'z';

      (* types defined in book definition *)

      alpha = packed array[1..namelength] of char; (* this is not alfa *)

      (* name is a left justified string with blanks following the
         characters *)
      name = record
         letters: alpha;
         length: 0..namelength (* zero means an unspecified structure *)
      end;

      lineptr = ^line;
      line = record (* a line of characters *)
         letters: packed array [1..linelength] of char;
         length: 0..linelength;
         next: lineptr
      end;

      direction = (plus, minus, dircomplement, dirhomologous);
      configuration = (linear, circular);
      state = (on, off);

      header = record (* header of key *)
         keynam: name; (* key name of structure *)
         fulnam: lineptr; (* full name of structure *)
         note:   lineptr (* note key *)
      end;

(* begin module base.type *)
   (* define the four nucleotide bases *)
   base = (a,c,g,t);
(* end module base.type version = 7.40; {of delmod.p 2000 Feb 18} *)

      (* sequence types *)
      dnaptr = ^dnastring;
      dnarange = 0..dnamax;
      seq = packed array[1..dnamax] of base;
      dnastring = record
         part: seq;
         length: dnarange;
         next: dnaptr
      end;


      orgkey = record (* organism key *)
         hea: header;
         mapunit: lineptr (* genetic map units *)
      end;

      chrkey = record (* chromosome key *)
         hea: header;
         mapbeg: real; (* number of genetic map beginning *)
         mapend: real (* number of genetic map ending *)
      end;

      pieceptr = ^piece;
      piekey = record (* piece key *)
         hea: header;
         mapbeg: real; (* genetic map beginning *)
         coocon: configuration; (* configruation (circular/linear) *)
         coodir: direction; (* direction (+/-) relative to genetic map *)
         coobeg: integer; (* beginning nucleotide *)
         cooend: integer; (* ending nucleotide *)
         piecon: configuration; (* configruation (circular/linear) *)
         piedir: direction; (* direction (+/-) relative to coordinates *)
         piebeg: integer; (* beginning nucleotide *)
         pieend: integer; (* ending nucleotide *)
      end;


      piece = record
         key: piekey;
         dna: dnaptr
      end;

      reference = record
         pienam : name;  (* name of piece referred to *)
         mapbeg : real;  (* genetic map beginning *)
         refdir : direction;  (* direction relative to coordinates *)
         refbeg : integer;  (* beginning nucleotide *)
         refend : integer;  (* ending nucleotide *)
      end;

      genkey = record  (* gene key *)
         hea : header;
         ref : reference;
      end;

      trakey = record  (* transcript key *)
         hea : header;
         ref : reference;
      end;

      markerptr = ^marker;
      markey = record (* marker key *)
         hea : header;
         ref : reference;
         sta : state;
         phenotype : lineptr;
         next : markerptr;
      end;

      marker = record
         key : markey;
         dna : dnaptr;
      end;

(* end module book.type version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module matrix.type *)
      matrix = array[base,1..maxmatrix] of integer; (* a pattern matrix *)
(* end module matrix.type version = 'auxmod 1.40 94 Sep 5 gds/tds'; *)

      alpieceptr = ^alpiece;
      alpiece = record
         length : integer;
         alignedbase : integer;
         dna : dnaptr;
         next : alpieceptr
      end;
      sequenceclass = boolean;
      ends = (fiveprime,threeprime); (* endpoints of ignoring regions *)

var
(* begin module book.var *)
(* ************************************************************************ *)
(* global variables needed for book manipulations *)

      (* free storage: *)
      freeline: lineptr; (* unused lines *)
      freedna: dnaptr; (* unused dnas *)

      readnumber: boolean; (* whether to read a number from the notes, or
                              to read in the notes *)
      number: integer; (* the number of the item just read *)
      numbered: boolean; (* true when the item just read is numbered *)
      skipunnum: boolean; (* a control variable to allow skipping of
                             un-numbered items in the book *)

(* ************************************************************************ *)
(* end module book.var version = 7.40; {of delmod.p 2000 Feb 18} *)

      patlrnp, (* parameter file, to specify 'maxtimes', the number of
              times through the sequences before stopping without a solution *)
      pattern, (* the output solution weighting function *)
      ignore,  (* file of integers, two per line, specifying the endpoints
                 of regions to be ignored (relative to aligned base) *)
      minmax, (* file of two integers, the first to set funcmin and the
              second to set nfuncmax.  if this file is empty they
              will start at 0 and -1 and vary with the matrix. *)
      start, (* file used to initialize wmatrix. if this file is empty
              it will be initialized as all 0's.  the format of
              the file must be the same as the output of this program. *)
      funcbook,  (* book of functional  sequences  *)
      nfuncbook,  (* book of nonfunctional sequences *)
      funcinst, (* inst for aligning functional pieces *)
      nfuncinst : text;  (* inst for aligning  nonfunctional pieces *)
      maxtimes, (* the number of times through the sequences before stopping
                  without a solution *)
      maxaligned,  (* maximum bases up to (and including)  aligned  base *)
      maxafter : integer; (* maximum bases after  aligned base *)
      width : integer;  (* width of seqmatrix (= maxaligned + maxafter) *)
      funcpie, (* pointer to list  of functional pieces *)
      nfuncpie : alpieceptr; (* pointer to list of nonfunctional pieces *)
      first : array[sequenceclass] of alpieceptr; (* pointer to first pieces  *)
      apiece : pieceptr; (* used by align routine to store piece *)
      wmatrix : matrix; (* weight function as a matrix *)
      changes,  (* number of times wmatrix has changed during this run *)
      funcmin,    (* minimum 'valu' for functional sequences and *)
      nfuncmax : integer; (* maximum 'valu' for nonfunctional sequences. *)
         (* see procedure evaluate for details of how 'valu' is determined *)
      varthresh, (* true if program determines funcmin and nfuncmax,
                     false if they are set by minmax *)
      done : boolean; (* true if program is to terminate *)
      seqclass : sequenceclass; (* the type of sequence being used *)
      mostseqs : sequenceclass; (* the class with the most sequences *)
      functional,  (* functional = true, nonfunctional = false *)
      nonfunctional : boolean; (* constants for seqclass variables *)
      adatetime : datetimearray;  (* for dating output *)
      timesthru : integer; (* number of times through all the sequences; must
                            be <= maxtimes *)
      ignorebases : array[ends,1..maxignore] of integer;  (* array of endpoints
                   of regions to ignore, read from input file ignore *)
      theline : integer; (* line number of the book *)


(* begin module package.primitive *)
(* ************************************************************************ *)
(* begin module halt *)
procedure halt;
(* stop the program.  the procedure performs a goto to the end of the
   program.  you must have a label:
      label 1;
   declared, and also the end of the program must have this label:
      1: end.
   examples are in the module libraries.
   this is the only goto in the delila system. *)
begin
      writeln(output,' program halt.');
      goto 1
end;
(* end module halt version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module copyaline *)
procedure copyaline(var fin, fout: text);
(* copy a line from file fin to file fout *)
begin (* copyaline *)
      while not eoln(fin) do begin
         fout^ := fin^;
         put(fout);
         get(fin)
      end;
      readln(fin);
      writeln(fout);
end; (* copyaline *)
(* end module copyaline version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module copylines *)
function copylines(var fin, fout: text; n: integer): integer;
(* copy n lines of file fin to file fout.
the actual number of lines copied is returned. *)
var
      index: integer; (* the current line number *)
begin (* copylines *)
      index := 0;
      while (not eof(fin)) and (index < n) do begin
         copyaline(fin, fout);
         index := succ(index)
      end;

      copylines := index
end; (* copylines *)
(* end module copylines version = 7.40; {of delmod.p 2000 Feb 18} *)

(* ************************************************************************ *)
(* end module package.primitive version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module package.trigger *)
(* ************************************************************************ *)

(* begin module interact.clearstring *)
procedure clearstring(var ribbon: string);
(* empty the string *)
var   index: integer; (* to the ribbon *)
begin (* clearstring *)
      with ribbon do begin
         for index := 1 to maxstring do letters[index] := ' ';
         length := 0;
         current := 0;
      end
end; (* clearstring *)
(* end module interact.clearstring version = 4.21; (@ of prgmod.p 1997 October 22 *)

(* begin module interact.writestring *)
procedure writestring(var tofile: text; var s: string);
(* write the string s to file tofile, no writeln *)
var   i: integer; (* index to s *)
begin (* writestring *)
      with s do for i := 1 to length do write(tofile, letters[i])
end; (* writestring *)
(* end module interact.writestring version = 4.21; (@ of prgmod.p 1997 October 22 *)

(* begin module filler.fillstring *)
procedure fillstring(var s: string; a: filler);
(* this procedure makes it reasonably easy to fill the string s with
characters.  one calls the procedure as: *)
(*                           1         2         3         4         5 *)
(*                  12345678901234567890123456789012345678901234567890 *)
(*   fillstring(s, 'this-is-the-string                                ');
the two comments make it easy to line the characters up. also, for this
example, it was assumed that the length of filler as defined by the
constant fillermax was 50. *)
var
      length: integer; (* of the string without trailing blanks *)
      index: integer; (* of s *)
begin (* fillstring *)
      clearstring(s);
      length := fillermax;
      while (length > 1) and (a[length] = ' ') do length := pred(length);
      if (length = 1) and (a[length] = ' ') then begin
         writeln(output, 'fillstring: the string is empty');
         halt
      end;

      for index := 1 to length do s.letters[index] := a[index];
      s.length := length;
      s.current := 1
end; (* fillstring *)
(* end module filler.fillstring version = 4.21; (@ of prgmod.p 1997 October 22 *)

(* begin module filler.filltrigger *)
procedure filltrigger(var t: trigger;
                      a: filler);
(* fill the trigger t *)
begin (* filltrigger *)
         fillstring(t.seek,a)
end; (* fillstring *)
(* end module filler.filltrigger version = 4.21; (@ of prgmod.p 1997 October 22 *)

(* begin module trigger.proc *)
(* this module allows one to scan a series of characters, as from
an array or a file, and to "trigger" or detect a simple string
in the series.  the advantage of the trigger is that several triggers
can "observe" a stream of characters at once, each looking for a
different thing.
some other modules required: interact.const, interact.type *)

procedure resettrigger(var t: trigger);
(* reset the trigger to ground state *)
begin (* resettrigger *)
      with t do begin
         state := 0;
         skip := false;
         found := false
      end
end; (* resettrigger *)

procedure testfortrigger(ch: char; var t: trigger);
(* look at the character ch.
   if it is part of the trigger (at the current trigger state),
       then the trigger state goes higher.
   if it is not part of the trigger then the trigger state is reset,
      skip is true and one should skip onward to find the trigger.
   if the trigger is found, found is true.
1996 Sep 12: Bug found!  In the case of a trigger "ab", the program
used to miss it for situations like "aab".  This was because at the
first a it would step up.  Then it would see the second a and recognize
that was not part of ab.  It would fail to realize that it could be
the start of a new one.  The code now accounts for that possibility. *)
begin (* testfortrigger *)
      with t do begin
         state := succ(state);
{
writestring(list,seek);
writeln(list,'testfortrigger seek.letters[',state:1,']:',
seek.letters[state],' ch:',ch);
}
         if seek.letters[state] = ch
         then begin
            skip := false;
            if state = seek.length then found := true
                                   else found := false
         end
         else begin
            (* it failed.  But wait!  It could be the beginning
            of a NEW trigger string! *)
            if seek.letters[1] = ch then begin
               state := 1;
               skip := false;
               found := false
            end
            else begin (* reset trigger *)
               state := 0;
               skip := true;
               found := false
            end
         end
      end
end; (* testfortrigger *)
(* end module trigger.proc version = 4.21; (@ of prgmod.p 1997 October 22 *)

(* begin module skipblanks *)
procedure skipblanks(var thefile: text);
(* skip over blanks until a non-blank, or end of line, is found *)
begin
      while (thefile^ = ' ') and not eoln(thefile) do get(thefile);
end;

procedure skipnonblanks(var thefile: text);
(* skip over nonblanks until a blank, or end of line, is found *)
begin
      while (thefile^ <> ' ') and not eoln(thefile) do get(thefile);
end;

procedure skipcolumn(var thefile: text);
(* skip over a data column *)
begin
   skipblanks(thefile); skipnonblanks(thefile)
end;
(* end module skipblanks version = 4.21; (@ of prgmod.p 1997 October 22 *)

(* ************************************************************************ *)
(* end module package.trigger version = 4.21; (@ of prgmod.p 1997 October 22 *)

(* begin module package.align *)
(* ************************************************************************ *)
(* begin module package.getpiece *)
(* ************************************************************************ *)
(* begin module package.brpiece *)
(* ************************************************************************ *)
(* begin module book.basis *)
(* procedures needed for book manipulations *)

(* get procedures should be used for all linked lists of records *)

procedure getline(var l: lineptr);
(* obtain a line from the free line list or by making a new one *)
begin
      if freeline<>nil
         then begin
            l:=freeline;
            freeline:=freeline^.next
         end
         else new(l);
      l^.length:=0;
      l^.next:=nil
end;

procedure getdna(var l: dnaptr);
begin
      if freedna<>nil
         then begin
            l:=freedna;
            freedna:=freedna^.next
         end
         else new(l);
      l^.length:=0;
      l^.next:=nil
end;

(* clear procedures should be called each time the records are no longer needed
         failure to do this may result in a stack overflow. *)

procedure clearline(var l: lineptr);
(* return a line to the free line list *)
var   lptr: lineptr;
begin
      if l<>nil then begin
         lptr:=l;
         l:=l^.next;
         lptr^.next:=freeline;
         freeline:=lptr
      end
end;

procedure showfreedna;
(* show the freedna list *)
var
   counter: integer; (* count of freedna list *)
   l: dnaptr; (* pointer into freedna list *)
begin
   l := freedna;
   counter := 0;
   while l <> nil do begin
      counter := succ(counter);
      write(output,counter:1);
      write(output, ', length = ',l^.length:1);
(* 2016 Jan 24 this is probably not standard Pascal, so hide it:
      write(output, ', pointer id: ',l:1);
*)
      writeln(output);
      l := l^.next
   end;
end;

procedure cleardna(var l: dnaptr);
var   lptr: dnaptr;
begin
      if l<>nil then begin
         lptr:=l;
         l:=l^.next;
         lptr^.next:=freedna;
         freedna:=lptr
      end
end;

procedure clearheader(var h: header);
(* clear the header h (remove lines to free storage) *)
begin
      with h do begin
         clearline(fulnam);
         while note<>nil do clearline(note)
      end
end;

procedure clearpiece(var p: pieceptr);
(* clear the dna of the piece *)
begin
      while p^.dna<>nil do cleardna(p^.dna);
      clearheader(p^.key.hea)
end;

function chartobase(ch:char):base;
(* convert a character into a base *)
begin
      case ch of
         'a': chartobase:=a;
         'c': chartobase:=c;
         'g': chartobase:=g;
         't': chartobase:=t
      end
end;

function basetochar(ba:base):char;
(* convert a base into a character *)
begin
      case ba of
         a: basetochar:='a';
         c: basetochar:='c';
         g: basetochar:='g';
         t: basetochar:='t';
      end
end;

function complement(ba:base):base;
(* take the complement of ba *)
begin
      case ba of
         a: complement:=t;
         c: complement:=g;
         g: complement:=c;
         t: complement:=a;
      end
end;

function chomplement(b: char): char;
(* create the character complement of base b.  I must be getting hungry! *)
begin
  chomplement := basetochar(complement(chartobase(b)));
end;

function pietoint(p: integer; pie: pieceptr): integer;
(* p is a coordinate on the piece.
   we want to transform p into a number
   from 1 to n: an internal coordinate system for
   easy manipulation of piece coordinates *)
var   i: integer; (* an intermediate value *)
begin
      with pie^.key do begin
         case piedir of
            plus:  if p>=piebeg
                   then i:=p-piebeg+1
                   else i:=(p-coobeg)+(cooend-piebeg)+2;
            minus: if p<=piebeg
                   then i:=piebeg-p+1
                   else i:=(cooend-p)+(piebeg-coobeg)+2
         end;
         pietoint:=i
      end
end;

function inttopie(i: integer; pie: pieceptr):integer;
(* i is in the range 1 to some maximum.  it is an internal coordinate
   system for the program.  we want to do a
   coordinate transformation to obtain
   a value in the range of the piece called pie:
   i=1 corresponds to piebeg and
   i=its maximum corresponds to pieend *)
var   p: integer; (* an intermediate value *)
begin
      with pie^.key do begin
         case piedir of
            plus: begin
               p:=piebeg+(i-1);
               if p>cooend
                  then if coocon=circular
                     then p:=p-(cooend-coobeg+1)
            end;
            minus: begin
               p:=piebeg-(i-1);
               if p<coobeg
                  then if coocon=circular
                     then p:=p+(cooend-coobeg+1)
            end
         end;
         inttopie:=p
      end
end;

function piecelength(pie: pieceptr): integer;
(* return the length of the dna in pie *)
begin
      piecelength:=pietoint(pie^.key.pieend,pie)
end;

(* end module book.basis version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.getto *)
function getto(var thefile: text; var theline: integer; ch: chset): char;
(* search the file for a character in the first line which is a
   member of the set ch.
Note:  on 1999 March 10 the definition of this function was cleaned
up.  Instead of putting thefile on the line AFTER the charcter ch
has been found, it puts thefile ON the line.  Other routines like
brdna and brpiece have to move to the next line themselves.  This makes
getto give the OBJECT.  *)
var
   achar: char; (* a character in thefile *)
   done: boolean; (* done finding achar *)
begin
   achar:=' ';
   done := false;
   while not done do begin
      if eof(thefile)
      then done := true
      else begin
         achar := thefile^;
         if achar in ch
         then done := true
         else begin
            readln(thefile);
            theline := succ(theline)
         end
      end
   end;
   if (achar in ch) then getto:=achar
                    else getto:=' '
{ The old method -
      while (not(achar in ch)) and (not eof(thefile))
         do begin
            readln(thefile,achar);
            theline := succ(theline)
         end;
   if (achar in ch) then getto:=achar
                    else getto:=' '
}
end;
(* end module book.getto version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.skipstar *)
procedure skipstar(var thefile: text);
(* skip start of line (or star = '*'). *)
begin (* skipstar *)
   if eof(thefile) then begin
      writeln(output,' procedure skipstar: end of book found');
      halt
   end
   else begin
      if thefile^ <> '*' then begin
         writeln(output,' procedure skipstar: bad book');
         writeln(output,' "*" expected as first character on the line, but "',
                        thefile^,'" was found');
         halt
      end;
      get(thefile); (* skip the star *)
      if thefile^ <> ' ' then begin
         writeln(output,' procedure skipstar: bad book');
         writeln(output,' "* " expected on a line but "*',
                        thefile^,'" was found');
         halt
      end;
      get(thefile) (* skip the blank *)
   end
end; (* skipstar *)
(* end module book.skipstar version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brreanum *)
procedure brreanum(var thefile: text; var theline: integer; var reanum: real);
(* read a real number from the file *)
begin
      skipstar(thefile);
      readln(thefile,reanum);
      theline := succ(theline)
end;
(* end module book.brreanum version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brnumber *)
procedure brnumber(var thefile: text; var theline: integer; var num: integer);
(* read a number from the file *)
begin
      skipstar(thefile);
      readln(thefile,num);
      theline := succ(theline)
end;
(* end module book.brnumber version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brname *)
procedure brname(var thefile: text; var theline: integer; var nam: name);
(* read a name from the file *)
var   i: integer; (* an index to the name *)
      c: char; (* a character read *)
begin (* brname *)
      skipstar(thefile);
      with nam do begin
         length:=0;
         repeat
            length:=succ(length);
            read(thefile,c);
            letters[length] := c
         until (eoln(thefile)) or
               (length>=namelength) or
               (letters[length]=' ');
         if letters[length]=' ' then length:=length-1;
         if length<namelength
            then for i:=length+1 to namelength do letters[i]:=' '
      end;
      readln(thefile);
      theline := succ(theline)
end; (* brname *)
(* end module book.brname version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brline *)
procedure brline(var thefile: text; var theline: integer; var l: lineptr);
(* read a line from the file *)
var
      i: integer;
      acharacter: char;
begin
      skipstar(thefile);
      i:=0;
      while (not eoln(thefile)) do begin
         i:=succ(i);
         read(thefile,acharacter);
         l^.letters[i]:=acharacter
      end;
      l^.length:=i;
      l^.next:=nil;
      readln(thefile);
      theline := succ(theline)
end;
(* end module book.brline version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brdirect *)
procedure brdirect(var thefile: text; var theline: integer; var direct: direction);
(* read a direction *)
var   ch: char;
begin
      skipstar(thefile);
      readln(thefile,ch);
      theline := succ(theline);
      if ch='+' then direct:=plus
                else direct:=minus
end;
(* end module book.brdirect version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brconfig *)
procedure brconfig(var thefile: text; var theline: integer; var config: configuration);
(* read a configuration *)
var   ch: char;
begin
      skipstar(thefile);
      readln(thefile,ch);
      theline := succ(theline);
      if ch='l' then config:=linear
                else config:=circular
end;
(* end module book.brconfig version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brnotenumber *)
procedure brnotenumber(var thefile: text; var theline: integer; var note: lineptr);
(* book note reading to obtain the number of the object.
the procedure returns the value of the number as a global.
(this is not such a good practice, but we are stuck with it for now.) *)
begin (* brnotenumber *)
      note:=nil;
      numbered := false;
      number := 0; (* force number to zero if there
         is no number at all *)
      (* the next character is n or * depending on whether there are notes *)
      if thefile^ = 'n' then begin
         readln(thefile);
         theline := succ(theline);
         if thefile^ <> 'n' then begin
            skipstar(thefile);
            if not eoln(thefile) then begin
               if thefile^ = '#' then begin
                  numbered := true;
                  get(thefile); (* move past the number symbol *)
                  read(thefile,number);
               end
            end;
            repeat
               readln(thefile);
               theline := succ(theline)
            until thefile^ = 'n';
            readln(thefile);
            theline := succ(theline)
         end
         else begin
            readln(thefile);
            theline := succ(theline)
         end
      end
end; (* brnotenumber *)
(* end module book.brnotenumber version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brnote *)
procedure brnote(var thefile: text; var theline: integer; var note: lineptr);
(* read note key *)
var
      newnote: lineptr; (* the new note *)
      previousnote: lineptr; (* the last line of the notes *)
begin (* brnote *)
      note:=nil;
      if thefile^ = 'n' then begin (* enter note *)
         readln(thefile);
         theline := succ(theline);
         if thefile^ <> 'n' then begin (* abort null note (n/n) *)
            getline(note);
            newnote:=note;
            while thefile^ <> 'n' do begin (* wait until end of note *)
               brline(thefile,theline,newnote);
               previousnote:=newnote;
               (* get next note *)
               getline(newnote^.next);
               newnote:=newnote^.next;
            end;
            (* last note was not used, so: *)
            clearline(newnote);
            previousnote^.next:=nil;
            readln(thefile);
            theline := succ(theline);
         end
         else begin
            readln(thefile);
            theline := succ(theline);
         end;
      end
end; (* brnote *)
(* end module book.brnote version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brheader *)
procedure brheader(var thefile: text; var theline: integer; var hea: header);
(* read the header of a key. *)
begin
      with hea do begin
         readln(thefile); (* move past the object name - new definition 1999 Mar 13 *)
         theline := succ(theline);
{bbb}
         (* read key name *)
         brname(thefile,theline,keynam);

         (* read full name *)
         getline(fulnam);
         brline(thefile,theline,fulnam);

         (* read note key *)
         if readnumber then brnotenumber(thefile,theline,note)
                       else brnote(thefile,theline,note)
      end
end;
(* end module book.brheader version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.copyheader *)
procedure copyheader(fromhea: header; var tohea: header);
(* copy the header fromhea into tohea.  Note that the
linked objects are NOT copied, but merely pointed to. *)
begin
   tohea.keynam.letters := fromhea.keynam.letters;
   tohea.keynam.length  := fromhea.keynam.length;
   tohea.note           := fromhea.note;
   tohea.fulnam         := fromhea.fulnam;
end;
(* end module book.copyheader version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brpiekey *)
procedure brpiekey(var thefile: text; var theline: integer; var pie: piekey);
(* read piece key, track the line number *)
begin
      with pie do begin
         brheader(thefile,theline,hea);
         brreanum(thefile,theline,mapbeg);
         brconfig(thefile,theline,coocon);
         brdirect(thefile,theline,coodir);
         brnumber(thefile,theline,coobeg);
         brnumber(thefile,theline,cooend);
         brconfig(thefile,theline,piecon);
         brdirect(thefile,theline,piedir);
         brnumber(thefile,theline,piebeg);
         brnumber(thefile,theline,pieend);
      end
end;
(* end module book.brpiekey version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brdna *)
procedure brdna(var thefile: text; var theline: integer; var dna: dnaptr);
(* read in dna from thefile, track the line *)
(* note: if the dna were circularized, by linking the last dnastring
to the first, then the cleardna routine could not clear properly,
and would loop forever... there is no reason to do that, since a simple
mod function will allow one to access the circle. *)
var
      ch: char;
      workdna: dnaptr;
begin
      getdna(dna);
      workdna:=dna;
      ch:=getto(thefile,theline,['d']);
      readln(thefile);
      theline := succ(theline);
      read(thefile,ch); (* skipstar *)
      while (ch = '*') do
      begin
         read(thefile,ch); (* skip blank *)
         repeat
            read(thefile,ch);
            if ch in ['a','c','g','t'] then begin
               if workdna^.length=dnamax then begin
                  getdna(workdna^.next);
                  workdna:=workdna^.next
               end;
               workdna^.length:=succ(workdna^.length);
               workdna^.part[workdna^.length]:=chartobase(ch)
            end
         until eoln(thefile);
         readln(thefile); (* go to next line *)
         theline := succ(theline);
         read(thefile,ch); (* ch is either '*' or 'd' *)
      end;
      readln(thefile); (* read past the d *)
      theline := succ(theline);
end;
(* end module book.brdna version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brpiece *)
procedure brpiece(var thefile: text; var theline: integer; var pie: pieceptr);
(* read in a piece, change theline to reflect the lines traversed *)
begin
{
   readln(thefile); (* move past the word 'piece' - new definition 1999 Mar 13 *)
   theline := succ(theline); (* BUG: was below! *)
bbb}
   brpiekey(thefile,theline,pie^.key);
   if numbered or (not skipunnum)
   then brdna(thefile,theline,pie^.dna);
   readln(thefile); (* move past the word 'piece' - new definition 1999 Mar 13 *)
   theline := succ(theline);
end;
(* end module book.brpiece version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.brinit *)
procedure brinit(var book: text; var theline: integer);
(* check that the book is ok to read, and
   set up the global variables for br routines *)
begin (* brinit *)
      (* halt if the book is bad (first word is 'halt') or the first
         character is not * *)
      reset(book);
      if not eof(book) then begin
         (* check for the date line *)
         if book^ <> '*' then begin
            if book^ <> 'h'
            then writeln(output, ' this is not the first line of a book:')
            else writeln(output, ' bad book:');
            write(output, ' ');

            while not (eoln(book) or eof(book)) do begin
               write(output, book^);
               get(book)
            end;
            writeln(output);
            halt
         end
      end
      else begin
         writeln(output, ' book is empty');
         halt
      end;

      (* initialize free storage *)
      freeline:=nil;
      freedna:=nil;

      readnumber:=true; (* usually we read in numbers for items *)
      number:=0; (* arbitrary value *)
      numbered:=false; (* the piece has no number (none yet read in) *)
      skipunnum:=false;
   theline := 1;
end; (* brinit *)
(* end module book.brinit version = 7.40; {of delmod.p 2000 Feb 18} *)

(* ************************************************************************ *)
(* end module package.brpiece version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module book.getpiece *)
procedure getpiece(var thefile: text; var theline: integer; var pie: pieceptr);
(* move to and read in the next piece in the book *)
var   ch: char;
begin
   ch:=getto(thefile,theline,['p']); (* get to the next p(iece) in the book *)
   if ch<>' ' then begin
      brpiece(thefile,theline,pie);
{
1999 june 2: removed this:
      ch:=getto(thefile,theline,['p']); (* read to end of p *)
}
{
bbb - now done in brpiece
      readln(thefile); (* read past piece *)
      theline := succ(theline);
}
   end
   else clearpiece(pie);
end;
(* end module book.getpiece version = 7.40; {of delmod.p 2000 Feb 18} *)

(* ************************************************************************ *)
(* end module package.getpiece version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module findblank *)
procedure findblank(var afile: text);
(* read a file to find the next blank character *)
var   ch: char;
begin
      repeat read(afile,ch) until ch = ' '
end;
(* end module findblank version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module findnonblank *)
procedure findnonblank(var afile: text; var ch: char);
(* find the next non blank character in a file, return it in ch. *)
begin
      ch:=' ';
      while (not eof(afile)) and
            (ch = ' ')
      do begin
         read(afile,ch);
         if eoln(afile) then readln(afile)
      end
end;
(* end module findnonblank version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module align.align *)
procedure align(var inst, book: text; var theline: integer;
                var pie:  pieceptr;
                var length, alignedbase: integer);
(* documentation on align is in module info.align and
delman.use.aligned.books.
1996 Sep 12: The routine now uses the trigger functions found
in prgmod.  The bug in the oldalign routine (that it misses the
end of comments that end in a series of astrisks) has been fixed.
It now checks that the piece corresponds to the book. *)
const maximumrange = 10000; (* if the alignment point is more than this
   distance from the piece ends, the program halts in an attempt to catch
   the alignment bug... 1991 Jan 11 It appears that the rewrite of the
   code has removed the bug, but the check will be kept. *)
   semicolon = ';'; (* end of delila instruction *)
var
   ch: char; (* a character in inst *)
   p: integer; (* index to a piece name *)
   p1: integer; (* another index to a piece name *)
   done: boolean; (* done finding an aligning get *)
   thebase: integer; (* the base read in *)
   indefault: boolean; (* true when within a default statement.  These
      can contain the word 'piece', which must be ignored. *)

   defaulttrigger: trigger; (* trigger to find 'default' *)
   gettrigger: trigger; (* trigger to find 'get' *)
   piecetrigger: trigger; (* trigger to find 'piece' *)
   begincomment: trigger; (* trigger to find '(-*' (ignore the dash!) *)
   endcomment: trigger; (* trigger to find '*-)'  (ignore the dash!) *)
   begincurly: trigger; (* trigger to find comments: '{' *)
   endcurly: trigger; (* trigger to find comments: '}' *)

{
procedure rd(var f: text; var ch: char);
(* read ch from f allowing inspection of the result *)
begin
   read(f,ch);
   write(output,ch);
   write(list,ch);
   write(output,'<',ch,'>');
end;

procedure rdln(var f: text);
(* readln f allowing inspection of the result *)
begin
   readln(f);
   writeln(output);
   writeln(list);
end;
}

procedure skipcomment(var f: text);
(* skip an entire comment *)
var
   comment: boolean; (* true means we are inside a comment *)
begin
   (* skip to end of comment *)
   resettrigger(endcomment);
   comment := true;
   while comment do begin
      if eof(f) then begin
         writeln(output,'A comment does not end!');
         halt
      end;
      if eoln(f) 
      then readln(f)
{ rdln(f) }
      else begin
{write(output,'<');
rd(f,ch);
write(output,'>');}
         read(f,ch);
         testfortrigger(ch, endcomment);
         if endcomment.found
         then comment := false;
      end
   end
end;

procedure skipcurly(var f: text);
(* skip an entire comment made by {}*)
var
   comment: boolean; (* true means we are inside a comment *)
begin
   (* skip to end of comment *)
   resettrigger(endcurly);
   comment := true;
   while comment do begin
      if eof(f) then begin
         writeln(output,'A comment does not end!');
         halt
      end;
      if eoln(f) 
      then readln(f)
{ rdln(f) }
      else begin
{write(output,'<');
rd(f,ch);
write(output,'>');}
         read(f,ch);
         testfortrigger(ch, endcurly);
         if endcurly.found
         then comment := false;
      end
   end
end;

begin
     filltrigger(gettrigger,'get    ');
   filltrigger(piecetrigger,'piece  ');
 filltrigger(defaulttrigger,'default');
   filltrigger(begincomment,'(*     ');
     filltrigger(endcomment,'*)     ');
     filltrigger(begincurly,'{      ');
       filltrigger(endcurly,'}      ');
    resettrigger(gettrigger);
  resettrigger(piecetrigger);
resettrigger(defaulttrigger);
  resettrigger(begincomment);
  resettrigger(begincurly);
   indefault := false;

   if not eof(book) then begin (* if there is still more to the book ... *)
      getpiece(book,theline,pie); (* read in the piece *)
      if not eof(book) then begin (* if we found a piece ... *)
         length:=pietoint(pie^.key.pieend,pie); (* calculate piece length *)

         (* now find in inst the next occurance of 'get' *)
         done := false;
         while not done do begin
            if eof(inst) then begin (* no instructions? *)
               alignedbase := 1; (* simply align by the first base *)
               done := true
            end
            else begin
               if eoln(inst)
               then readln(inst)
{then rdln(inst)}
               else begin
{rd(inst,ch);}
                  read(inst,ch);

                  testfortrigger(ch, begincomment);
                  testfortrigger(ch, begincurly);
                  if begincomment.found or begincurly.found
                  then begin
                     if ch = '*' then begin
                        skipcomment(inst);
                        resettrigger(begincomment);
                     end
                     else begin
                        resettrigger(begincurly);
                        skipcurly(inst);
                     end
                  end
                  else begin

                     testfortrigger(ch, gettrigger);
                     if gettrigger.found then begin
                        findnonblank(inst,ch); (* get to "from" *)
                        findblank(inst); (* get past "from" *)
                        read(inst,thebase); (* read in the alignedbase *)
{writeln(output);writeln(output,'thebase = ',thebase:1);}
                        alignedbase:=pietoint(thebase,pie);
{writeln(output,'alignedbase=',alignedbase:1);}
                        done := true
                     end;

                     testfortrigger(ch, defaulttrigger);
                     if defaulttrigger.found then begin
                        indefault := true;
                        resettrigger(defaulttrigger)
                     end;
                     if ch = semicolon then indefault := false;

                     (* check that piece names are correct *)
                     testfortrigger(ch, piecetrigger);
                     if not indefault then if piecetrigger.found then begin
                        skipblanks(inst); (* get to name *)
                        with pie^.key.hea.keynam do begin
                           for p := 1 to length do begin
                              read(inst,ch);
                              if letters[p] <> ch then begin
                                 writeln(output,'Piece name in book: ');
                                 writeln(output,letters:length);
                                 writeln(output,'does not match',
                                              ' the inst file name:');
                                 (* write the letters that matched: *)
                                 for p1 := 1 to p-1 do write(output,letters[p]);
                                 (* write the offending letter: *)
                                 write(output, ch);
                                 (* get the rest of the name and show it: *)
                                 done := eoln(inst);
                                 while not done do begin
                                   done := eoln(inst);
                                   if not done then begin
                                     read(inst,ch);
                                     if (ch = ' ') or (ch = ';') then done := true;
                                     if not done then write(output,ch);
                                   end;
                                 end;
                                 writeln(output);
                                 (* mark the first letter that does not match: *)
                                 for p1 := 1 to p-1 do write(output,' ');
                                 write(output,'^');
                                 writeln(output);
                                 halt
                              end;
                           end
                        end;
                     end;

                  end
               end
            end
         end;

         if (alignedbase <= -maximumrange) or
            (alignedbase > length + maximumrange) then begin
            writeln(output,' In procedure align:');
            writeln(output,' read in base was ',thebase:1);
            writeln(output,' in internal coordinates: ',alignedbase:1);
            writeln(output,' maximum range was ',maximumrange:1);
            writeln(output,' piece length was ',length:1);
            with pie^.key.hea.keynam
               do writeln(output,' piece name: ',letters:length);
            writeln(output,' piece number: ',number:1);
            writeln(output,' aligned base is too far away... see the code');
            halt
         end
      end
   end
end;
(* end module align.align version = 7.40; {of delmod.p 2000 Feb 18} *)

(* ************************************************************************ *)
(* end module package.align version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module package.datetime *)
(* ************************************************************************ *)
(* begin module getdatetime *)
procedure getdatetime(var adatetime: datetimearray);
(* get the date and time into a single array from the system clock.
   adatetime contains the date:
       1980/06/09 18:49:11
         ye mo da ho mi se
(year, month, day, hour, minute, second).
As of 2000 February 18, the Sun Pascal compiler requires a formatting
statement.  This statement allows the date to be generated in this
standard Delila format in a single call.  Information about the
formatting statement is available on the manual page for date in Unix.
If a computer does not have this method, see the 'oldgetdatetime' routine
in delmod.p (https://alum.mit.edu/www/toms/delila/delmod.html)
for some conversion code.
*)
begin
   date(adatetime,'%Y/%m/%d %H:%M:%S');
end;
(* end module getdatetime version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module readdatetime *)
procedure readdatetime (var thefile: text; var adatetime: datetimearray);
(* read the date and time from the file *)
var
      index: integer; (* to the udatetime *)
      (* the following is an unpacked date time array, to avoid reading
      into a packed array.  reading into a packed array is not transportable *)
      udatetime: array[1..datetimearraylength] of char;
begin
      for index:=1 to datetimearraylength do read(thefile,udatetime[index]);
      pack(udatetime, 1, adatetime);
      if (adatetime[3]='/') and (adatetime[12]=':') then begin
         writeln(output,' old datetime (only 2 year digits) read: ',
                          adatetime:datetimearraylength);
      end;
end;
(* end module readdatetime version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module writedatetime *)
procedure writedatetime(var thefile: text; adatetime: datetimearray);
(* expand the date and time out and print in the file *)
var
      index: integer; (* index of datetime *)
begin
      for index:=1 to datetimearraylength do write(thefile,adatetime[index])
end;
(* end module writedatetime version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module timeseed *)
procedure addtoseed(var seed, power: real; c: char);
(* add the digit represented by c to the seed at the power position *)
var
      n: integer; (* the character represented by c *)
begin (* addtoseed *)
      power := power/10;
{
writeln(output,'addtoseed, c = ',c);
writeln(output,'addtoseed, ord(c) = ',ord(c));
}
      case c of
         ' ': begin
                 writeln(output,'timeseed: error in datetime');
                 halt;
              end;
         '0': n := 0;
         '1': n := 1;
         '2': n := 2;
         '3': n := 3;
         '4': n := 4;
         '5': n := 5;
         '6': n := 6;
         '7': n := 7;
         '8': n := 8;
         '9': n := 9
      end;
(*writeln(output,'timeseed number is [',n:1,']'); (@ debug *)
      seed := seed + power*n
end; (* addtoseed *)

procedure makeseed(adatetime: datetimearray; var seed: real);
(* convert adatetime to a real number in seed, reversed order *)
var
   power: real; (* a digit of the seed such as 0.01 *)
begin
   seed := 0.0;
   power := 1.0;
   addtoseed(seed, power, adatetime[19]);
   addtoseed(seed, power, adatetime[18]);
   (* : *)
   addtoseed(seed, power, adatetime[16]);
   addtoseed(seed, power, adatetime[15]);
   (* : *)
   addtoseed(seed, power, adatetime[13]);
   addtoseed(seed, power, adatetime[12]);
   (*   *)
   addtoseed(seed, power, adatetime[10]);
   addtoseed(seed, power, adatetime[9]);
   (* / *)
   addtoseed(seed, power, adatetime[7]);
   addtoseed(seed, power, adatetime[6]);
   (* / *)
   addtoseed(seed, power, adatetime[4]);
   addtoseed(seed, power, adatetime[3]);
   addtoseed(seed, power, adatetime[2]);
   addtoseed(seed, power, adatetime[1]);
end;

procedure orderseed(adatetime: datetimearray; var seed: real);
(* convert adatetime to a real number in seed, normal order *)
var
   power: real; (* a digit of the seed such as 0.01 *)
begin
   seed := 0.0;
   power := 1.0;
   addtoseed(seed, power, adatetime[1]);
   addtoseed(seed, power, adatetime[2]);
   addtoseed(seed, power, adatetime[3]);
   addtoseed(seed, power, adatetime[4]);
   (* / *)
   addtoseed(seed, power, adatetime[6]);
   addtoseed(seed, power, adatetime[7]);
   (* / *)
   addtoseed(seed, power, adatetime[9]);
   addtoseed(seed, power, adatetime[10]);
   (*   *)
   addtoseed(seed, power, adatetime[12]);
   addtoseed(seed, power, adatetime[13]);
   (* : *)
   addtoseed(seed, power, adatetime[15]);
   addtoseed(seed, power, adatetime[16]);
   (* : *)
   addtoseed(seed, power, adatetime[18]);
   addtoseed(seed, power, adatetime[19]);
end;

procedure timeseed(var seed: real);
(* read the computer date and time.  reverse the order of the digits
and put a decimal point in front.  this gives a fraction between
zero and one that varies quite quickly, and is always unique (if the
computer has sufficient accuracy).  it is to be used as a seed to
a random number generator. *)
var
   adatetime: datetimearray; (* a date and time *)
begin (* timeseed *)
   getdatetime(adatetime);
{
writeln(output,'timeseed: adatetime: ',adatetime);
}
   makeseed(adatetime, seed);
end; (* timeseed *)
(* end module timeseed version = 7.40; {of delmod.p 2000 Feb 18} *)
 
(*[[*)
(* begin module limitdate *) 
procedure limitdate(a,b,c,d: char; limitdatetime: datetimearray);
(* test whether the current time is before the limit.
If it is later, halt the program *)
var
   adatetime: datetimearray; (* a date and time *)
   Dday: real; (* the critical day *)
   now: real;  (* this very moment *)
begin
   getdatetime(adatetime);
{
writeln(output,'adatetime:',adatetime);
writeln(output,'adatetime[1]: ', adatetime[1]);
writeln(output,'adatetime[2]: ', adatetime[2]);
writeln(output,'adatetime[3]: ', adatetime[3]);
writeln(output,'adatetime[4]: ', adatetime[4]);
writeln(output,'adatetime[5]: ', adatetime[5]);
}

   orderseed(adatetime, now);

   if 
      (limitdatetime[1] <> ' ') or
      (limitdatetime[2] <> ' ') or
      (limitdatetime[3] <> ' ') or
      (limitdatetime[4] <> ' ') then halt;

   limitdatetime[1] := a;
   limitdatetime[2] := b;
   limitdatetime[3] := c;
   limitdatetime[4] := d;

   orderseed(limitdatetime, Dday);

{
   writeln(output,'now:  ',now:20:10);
   writeln(output,'Dday: ',Dday:20:10);
}
   if now > Dday then begin
       writeln(output,'This program expired on ',limitdatetime);
       writeln(output,'See: https://alum.mit.edu/www/toms/walker/contacts.html');
       halt
   end
end;
(* end module limitdate version = 7.40; {of delmod.p 2000 Feb 18} *)

(*]]*)
(* ************************************************************************ *)
(* end module package.datetime version = 7.40; {of delmod.p 2000 Feb 18} *)

function getalbase(position: integer; pie: alpieceptr):base;
(* get a base from the nth position (internal coordinates) of the
   piece.  no protection is made against positions outside the piece.
   this procedure does the same thing as the standard getbase procedure
   from delmods, except that it uses for the variable 'pie' one of type
   alpieceptr instead of one of type pieceptr.  the program could be
   changed to the more standard procedure, and it probably will be ... *)
var
      workdna: dnaptr;
      p: integer; (* the last base of the dna part *)
begin
      workdna:=pie^.dna;
      p:=dnamax;
      while position>p do begin
         p:=p+dnamax;
         workdna:=workdna^.next
      end;
      getalbase:=workdna^.part[position-(p-dnamax)]
end;

(* begin module matrix.readmatrix *)
procedure findlearnend(var pattern: text);
(* locate the end of the cyclic learning part of the pattern file *)
const
      learnmax = 16; (* the size of learnend *)
var
      (* this array contains the pattern that defines the end *)
      learnend: packed array[1..learnmax] of char;
      state: integer; (* how close we are to finding learnend *)
begin (* findlearnend *)
      reset(pattern);
      learnend := 'end of learning.';
      state := 1;
      while not eof(pattern) and (state < learnmax) do begin
         if learnend[state] = pattern^
         then state := succ(state)
         else state :=1;
         if eoln(pattern) then readln(pattern)
                          else get(pattern)
      end;
      if eof(pattern) then begin
         write  (output,'pattern matrix does not contain "');
         writeln(output,learnend,'" signal');
         halt
      end
end; (* findlearnend *)
procedure getcolon(var f: text);
(* move to the next colon in file f *)
begin (* getcolon *)
      while f^<>':' do begin
         get(f);
         if eof(f) then begin
            writeln(output,'pattern is missing colons');
            halt
         end
      end;
      get(f); (* move past the colon *)
end; (* getcolon *)

procedure readmatrix(var thefile: text; var wmatrix: matrix;
                     var beginning, width: integer);
(* read a pattern matrix wmatrix from the file thefile, returning
the aligned position for the first position (beginning) and the width *)
var   i : integer;  (* an index *)
      ba : base;  (* an index *)
begin (* readmatrix *)
      findlearnend(thefile);

      (* skip min and max values *)
      getcolon(thefile); getcolon(thefile);

      getcolon(thefile); (* skip to range *)
      read(thefile,beginning); (* read first position of the matrix *)

      getcolon(thefile); (* skip to width *)
      read(thefile,width); (* read width of the matrix *)

      if (width > maxmatrix) then begin
         writeln(output,'input matrix too large');
         halt;
      end;

      for i := 1 to width do   (* read in matrix *)
      begin
         getcolon(thefile); (* skip position label *)
         for ba := a to t do read(thefile,wmatrix[ba,i]);
         readln(thefile);
      end;
            (* set all the unused elements of wmatrix to 0 *)
      for i := (width + 1) to maxmatrix do
         for ba := a to t do wmatrix[ba,i] := 0;
end; (* readmatrix *)
(* end module matrix.readmatrix version = 'auxmod 1.40 94 Sep 5 gds/tds'; *)

procedure initialize; (* initializes all variables, including wmatrix *)
var ba  : base;
    beginning, (* of the starting matrix *)
    i : integer;
begin (* initialize *)
      writeln(output,'patlrn ',version:4:2);
      brinit(funcbook, theline);
      brinit(nfuncbook, theline);
      reset(funcinst);
      reset(nfuncinst);
      reset(minmax);
      reset(start);
      rewrite(pattern);

      functional := true; (* these are constants *)
      nonfunctional :=  false;

      beginning := 0;
      maxafter  := 0;
      first[functional] := nil;
      first[nonfunctional] := nil;
      changes := 0;
      timesthru := 0;
      done := false;
            (* wmatirx initialization *)
      if not eof(start) then readmatrix(start,wmatrix,beginning,width)
      else for ba := a to t do
         for i := 1 to maxmatrix do wmatrix[ba,i] := 0;
      new(apiece);
      maxaligned := 1 - beginning;
end;

procedure patheader;
(* write the pattern header to file pattern *)
begin (* patheader *)
      writeln(pattern,' patlrn ',version:4:2,' pattern learning program');

      getdatetime(adatetime);
      write(pattern,' ');
      writedatetime(pattern,adatetime);
      writeln(pattern);

      writeln(pattern,' functional sequences are from:');
      write(pattern,' ');
      copyaline(funcbook,pattern);

      writeln(pattern,' nonfunctional sequences are from:');
      write(pattern,' ');
      copyaline(nfuncbook,pattern);

      writeln(pattern);
      reset(start);
      if not eof(start) then
      begin
         readln(start);
         write(pattern,' matrix initialized as ');
         copyaline(start,pattern)
      end
      else writeln(pattern,' matrix initialized as all 0''s');
      writeln(pattern);
end;

procedure setparam; (* read from file 'patlrnp' the value of 'maxtimes',
                       the number of times through the sequences before
                       stopping without a solution *)
begin
      reset(patlrnp);
      if not eof(patlrnp) then readln(patlrnp,maxtimes)
      else maxtimes := 100;
end;

procedure readignore;  (* this reads from the input ignore the regions
               which are not to be considered in the pattern leaarning.
               the format of the ignore file is two integers per line,
               up to maxignore lines, where the first integer is the
               beginning of the region to be ignored and the second is
               the end (relative to the aligned base). the last can be
               equal to the first but not less than it.  *)
var   index : integer; (* the number of regions to be ignored *)

begin
      reset(ignore);
      for index := 1 to maxignore do ignorebases[fiveprime,index] := -10000;
      index := 1;
      if not eof(ignore) then
         writeln(pattern,' the following regions are ignored by',
            ' setting wmatrix to 0 :');
      while (not eof(ignore) and (index <= maxignore)) do
      begin
         read(ignore,ignorebases[fiveprime,index]);
         readln(ignore,ignorebases[threeprime,index]);
         if ignorebases[fiveprime,index] > ignorebases[threeprime,index] then
         begin
            writeln(output,' ignore file format error - fiveprime greater',
               ' than threeprime');
            halt
         end;
         writeln(pattern,ignorebases[fiveprime,index]:6,' to ',
            ignorebases[threeprime,index]);
         index := succ(index);
      end;
      if ((index > maxignore) and not eof(ignore)) then
      begin
         writeln(output,' ignore file has too many lines; increase the',
            ' constant "maxignore" and recompile');
         halt
      end;
      writeln(pattern,' ');
end;

procedure ignoring;  (* sets the regions specified by the file ignore
                           to zero in wmatrix *)
var   i,index : integer;
      ba : base;

begin
      index := 1;
      while ((ignorebases[fiveprime,index] <> -10000) and
            (index <= maxignore)) do
      begin
         for i := ignorebases[fiveprime,index] to ignorebases[threeprime,index]
            do for ba := a to t do wmatrix[ba,maxaligned + i]  := 0;
         index := succ(index);
      end;
end;

procedure makealignedpieces;
      (* we make a linked list of the functional sequences and their
         alignments, and then do the same for the nonfunctional  sequences.
         the maximum sequence length is determined by procedure matrixwidth *)
         (* funcmin and nfuncmax vary with changes to wmatrix (increasing
            or decreasing as wmatrix is added to or subtracted from)
            or they can be set by the file minmax, which has two
            integers, the first for funcmin and the second for
            nfuncmax. *)

var  seqnumb : array[sequenceclass] of integer;

procedure alignedpiece(alpie : alpieceptr; var book,inst : text);

procedure matrixwidth; (* this determines maxaligned and
               maxafter, which determine the pattern width *)
begin
      if maxaligned < alpie^.alignedbase then
         maxaligned := alpie^.alignedbase;
      if maxafter < (alpie^.length - alpie^.alignedbase) then
         maxafter := (alpie^.length - alpie^.alignedbase);
      if ((maxaligned + maxafter) > maxmatrix) then
      begin
         write(output,'constant maxmatrix is too small');
         halt
      end;
end;

begin (* alignedpiece *)
      seqnumb[seqclass] := 0;
      while not eof(book) do
      begin
         new(alpie);
         align(inst,book,theline,apiece,alpie^.length,alpie^.alignedbase);
         if not eof(book) then
         begin
            seqnumb[seqclass] := succ(seqnumb[seqclass]);
            alpie^.dna := apiece^.dna;
            apiece^.dna := nil;
            alpie^.next  := first[seqclass];
            first[seqclass]  := alpie;
            clearpiece(apiece);
            matrixwidth;
         end;
      end;
end;

begin   (* makealignedpieces *)
      seqclass := functional;
      alignedpiece(funcpie,funcbook,funcinst);
      seqclass := nonfunctional;
      alignedpiece(nfuncpie,nfuncbook,nfuncinst);
      if (seqnumb[functional] > seqnumb[nonfunctional]) then
         mostseqs := functional
      else mostseqs := nonfunctional;
      width := maxaligned + maxafter;
      if eof(minmax) then begin
         funcmin := 0;
         nfuncmax := -1;
         varthresh := true;
      end
      else begin
         read(minmax,funcmin);
         read(minmax,nfuncmax);
         varthresh := false;
      end;
end;

procedure evaluatew; (* this chooses a sequence and evaluates
               it by wmatrix.  if the
               required condition (given by funcmin and nfuncmax)
               is not met wmatrix is changed accordingly *)


var   thispie : alpieceptr;

procedure chooseseq;

procedure arewedone; (* we are done if (1) the  sequence list we  just
               came to the end of  (seqclass) is the longest
               list  (mostseqs) and (2) the number of changes
               is equal to 0 or we have been through all the sequences the
               maximum allowed times (maxtimes) *)

begin
      if (mostseqs = seqclass) then
      begin
         writeln(pattern,changes:8);
         if ((changes = 0) or (timesthru = maxtimes)) then done := true
         else begin
            changes := 0;
            timesthru := succ(timesthru);
         end;
      if (changes <> 0) then
         writeln(pattern,' maximum allowed times through sequences : ',
            timesthru:2);
      end;
end;

begin  (* chooseseq *)
      if (seqclass = functional) then
      begin
         thispie := funcpie;
         if (thispie = nil) then
         begin
            arewedone;
            if not done then
            begin
               funcpie := first[functional];
               thispie := funcpie
            end
         end
         else funcpie := funcpie^.next
      end
      else begin
         thispie := nfuncpie;
         if (thispie = nil) then
         begin
            arewedone;
            if not done then
            begin
               nfuncpie := first[nonfunctional];
               thispie := nfuncpie
            end
         end
         else nfuncpie := nfuncpie^.next;
      end;
end;

procedure evaluate; (* wmatrix is evaluated as follows: the corresponding
               elements of seqmatrix and wmatrix are multiplied and all the
               resulting products are summed to give 'valu'. if the sequence
               is functional and valu is < funcmin, wmatrix is changed by
               adding seqmatrix to it. if the sequence is nonfunctional and
               valu is > nfuncmax, wmatrix is changed by  subtracting seqmatrix
               from it. after each test of wmatrix by a sequence, a sequence
               of the other class is chosen to test it next and we repeat
               this  until done = true.  *)
var   missing :  integer; (* the number fewer bases preceding the aligned
               base in this sequence than in the maximum *)
      valu : integer;
{
variable ba is neither used nor set
      ba : base;
}
      i : integer;

begin
      missing := (maxaligned - thispie^.alignedbase);
      valu := 0;
      for i := 1 to thispie^.length do
         valu := valu + (wmatrix[getalbase(i,thispie),missing + i]);
      if ((seqclass = functional) and (valu < funcmin)) then
      begin
         changes := succ(changes);
         for i := 1 to thispie^.length do
            wmatrix[getalbase(i,thispie),missing + i] :=
               succ(wmatrix[getalbase(i,thispie),missing + i]);
         if varthresh then begin
            funcmin := pred(funcmin);
            nfuncmax := pred(nfuncmax);
         end;
         ignoring;
      end;
      if ((seqclass = nonfunctional) and (valu > nfuncmax)) then
      begin
         changes := succ(changes);
         for i := 1 to thispie^.length do
            wmatrix[getalbase(i,thispie),missing + i] :=
               pred(wmatrix[getalbase(i,thispie),missing + i]);
         if varthresh then begin
            funcmin := succ(funcmin);
            nfuncmax := succ(nfuncmax);
         end;
         ignoring;
      end;
      seqclass := not(seqclass);
end;

begin (* evaluatew *)
      chooseseq;
      if not done then evaluate;
end;

procedure printoutw;
(* print out the w matrix to the pattern file *)
var   ba : base;
      i : integer;
begin (* printoutw *)
      writeln(pattern,' end of learning.');
      writeln(pattern);
      writeln(pattern,' minimum value for a functional sequence is: ',
         funcmin:1);
      writeln(pattern,' maximum value for a nonfunctional sequence is: ',
         nfuncmax:1);
      writeln(pattern);

      write(pattern,' range: ',(1-maxaligned):1,' to ',(width-maxaligned):1);
      writeln(pattern,'  width: ',width:1);
      writeln(pattern);
         (* write the matrix vertically to avoid being too wide to transport *)
            (* label the columns *)
      writeln(pattern,'  pos. ','a':5,'c':5,'g':5,'t':5);
            (* make the upper matrix-border *)
      write(pattern,' ':6);
      for i := 1 to 23 do write(pattern,'-');
      writeln(pattern);
            (* write the matrix positions and elements *)
      for i := 1 to width do begin
         write(pattern,' ',(i-maxaligned):4,' :');
         for ba := a to t do write(pattern,' ',wmatrix[ba,i]:4);
         writeln(pattern,' :');
      end;
            (* make the lower matrix-border *)
      write(pattern,' ':6);
      for i := 1 to 23 do write(pattern,'-');
      writeln(pattern);
end; (* printoutw *)

begin (* patlrn *)
      initialize;
      setparam;
      patheader;
      readignore;

      writeln(pattern,' the number of changes to wmatrix in each round');
      makealignedpieces;
      funcpie := first[functional];
      nfuncpie := first[nonfunctional];
      seqclass := mostseqs; (* we begin with a sequence from the most
                              abundant class *)
      while not done do evaluatew;

      printoutw;
1: end. (* patlrn *)
