program encode(inst,book,encseq,encodep,output);
(* encode a book of sequences into strings of integers.
  by Gary Stormo.

modified by:

  Dr. Thomas D. Schneider
  toms@alum.mit.edu
  permanent email: toms@alum.mit.edu (use only if first address fails)
  https://alum.mit.edu/www/toms/

modules needed: delmods, delman, matmods *)

label 1;  (* the end of the program *)

const
(* begin module version *)
   version = 1.42;  (* of encode.p 2007 Jun 22
2007 Jun 22, 1.42: handle mutated names; load from delmod
2002 Nov  1, 1.41: fix output message
2000 Nov 20, 1.40: upgrade documentation and to gpc
2000 Nov 17, 1.38: upgrade to new delmod for reading name inst.
1999 Jun 11, 1.37: upgrade to theline in alignedlist stuff
1995 Jan 29: empty inst implies book alignment
origin july 7, 1983 *)
(* end module version *)

(* begin module describe.encode *)
(*
name
   encode: encodes a book of sequences into strings of integers

synopsis
   encode(inst: in, book: in, encseq: out, encodep: in, output: out)

files
   inst: the instructions generating the book; for aligning the sequences
       If the inst file is empty, then the sequences are aligned by
       the zero coordinate of the book (this allows the use of the
       "default coordinate zero" option of Delila) or by the first
       base of the piece, as defined by the first parameter.

   book: the sequences to be encoded

   encseq: the encoded sequences

   encodep: parameter file for describing how the sequences are to be
         encoded.

   The first parameter, the first character on the first line, defines how
   to align the pieces.  See the alist program for the detailed logic.
   There are three choices, as in alist:

      'f' (for 'first') then the sequences are always aligned by their
      first base.

      'i' then the sequences are aligned by the delila instructions.  If
      the inst file is empty, alignment is forced to the 'b' mode.

      'b' (for 'internal') then the alignment is on the internal zero of
      the book's sequence.  This option is to be used when "default
      coordinate zero" is used in the Delila instructions.

   The remaining parameters are stored as a list of parameter records, of
   which there may be any number.  Each parameter record has five lines of
   information which it must include (all i's and j's are integers):

   1.  i j specify the nucleotides, relative to the aligned base,
           over which this parameter record is to operate; these may
           be any integers, but i <= j is required;
   2.  i   is the size of the windows to be encoded; within the window
           the number of each oligonucleotide of length 'coding' are
           determined and printed as part of the total sequence vector;
   3.  i   is the shift to the next window to be encoded;
   4.  i : j1 j2 j3 ...  is the 'coding'-level and arrangement; the
           'coding'-level, i, is the number of nucleotides in the oligos we
           are counting, i.e., 1 means monos, 2 means dis, ...;  if i > 1
           then we can also skip bases between the ones we are encoding;
           if the i is followed next by a colon, there must be i-1 integers
           (j1..j(i-1)) which specify the number of bases to be skipped
           between the ones which are encoded; for example, if we have the
           sequence xyz and we are interested in the di-nucleotides we can
           get the xy by the parameter '2 : 0', or we could get the xz by
           parameter '2 : 1'; if there is no colon all the skips are
           assumed to be zero;
   5.  i   is the shift to the next coding site within the window;
           this allows us to encode only some of the oligos within a window,
           such as only those that are in-frame;
   multiple parameter records can be concatenated in the encodep file
   and then each sequence in the book will be encoded according to each
   parameter record into a single vector of integers.

   output: for messages to the user

description

   This program is used to encode a book of sequences into a string of
   integers.  Each sequence in the book is encoded into a single string of
   integers (ended by an 'end of sequence' symbol) according to the user
   specified parameters, which are in the file 'encodep'.

examples

documentation

@article{Schneider1984,
author = "T. D. Schneider
 and G. D. Stormo
 and M. A. Yarus
 and L. Gold",
title = "Delila system tools",
journal = "Nucl. Acids Res.",
volume = "12",
pages = "129-140",
year = "1984"}

see also

   {Example parameter file:} encodep

   {delman.use.encode: }
   https://alum.mit.edu/www/toms/delman1.html#delman.use.encode.1

   {delman.use.aligned.books:}
   https://alum.mit.edu/www/toms/delman1.html#delman.use.aligned.books

   {Before using encode, one should always check the sequences
    by looking at them as an aligned list with} alist.p

   {The output of this program is used by} rseq.p

author

   Gary Stormo

bugs

      none known

technical notes

*)
(* end module describe.encode *)

(* begin module book.const *)
(* constants needed for book manipulations *)
      dnamax = 1024;  (* length of dna arrays *)
      namelength = 100; (* maximum key name length *)
      linelength = 200; (* maximum line readable in book *)
(* end module book.const version = 7.72; {of delmod.p 2007 Jul 23} *)

(* begin module vector.const *)
      maxvecpart = 64; (* the number of elements in a 'part' of a vector *)
      vpagewidth = 64; (* the width (in characters) of the output vector
         file from procedure writevector *)
(* end module vector.const version = 'matmod 1.95 85 apr 18 tds/gds'; *)

      pagewidth = 32; (* number of encoding elements written per line *)
      codingmax = 6; (* the maximum coding-level allowed *)

(* begin module interact.const *)
(* begin module string.const *)
      maxstring = 2000; (* the maximum string *)
(* end module string.const version = 4.86; (@ of prgmod.p 2004 Sep 8 *)
(* end module interact.const version = 7.72; {of delmod.p 2007 Jul 23} *)
 
(* begin module filler.const *)
      fillermax = 50; (* the size of the filler array for a string *)
(* end module filler.const version = 4.18; (@ of prgmod.p 1996 September 12 *)

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.86; (@ of prgmod.p 2004 Sep 8 *)
(* end module interact.type version = 7.72; {of delmod.p 2007 Jul 23} *)

(* 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.18; (@ of prgmod.p 1996 September 12 *)

(* 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.18; (@ of prgmod.p 1996 September 12 *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

      (* 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.72; {of delmod.p 2007 Jul 23} *)

(* begin module encode.type.param *)
(* the following types allow the user to specify parameters which will
   be used to encode the sequences in the book. *)
(* spaces are allowed between the encoded bases, and the number of
   bases to be skipped between each encoded pair of bases are kept
   in this linked list of integers *)
      spaceptr = ^spacelist;
      spacelist = record
         skips : integer; (* bases skipped to next coded base *)
         next : spaceptr; (* points to next spacing number *)
      end;

(* the encoding parameters for each region are stored in these records,
   and the records for each region are connected into a linked list of
   all the encoding parameters for the entire sequences *)
      endpoints = (start,stop);  (* end points of a coding region *)
      paramptr = ^parameter;
      parameter = record  (* these are the instructions for doing the coding *)
         range : array[endpoints] of integer; (* the bases to be coded by these
                                 instructions, relative to the alignedbase *)
         window : integer; (* the number of bases included in a coding vector *)
         wshift : integer;  (* movement of the window to the next site *)
         coding : integer; (* the 'level' at which the coding is being done,
                             i.e., 1: mononucleotides; 2: dinucleotides; ... *)
         spaces : spaceptr;  (* the spacing between the encoded bases *)
         cshift : integer;  (* shift to the next coding unit in the window *)

         (* values calculated at read-in time *)
         wvlength: integer; (* length of a window vector.  this is coding
            raised to the 4th power *)
         pvlength: integer; (* parameter vector length.  the vector made up
            of a series of window vectors. *)
         (* note: pvlength/wvlength is the number of windows in the
            parameter range *)

         next : paramptr;  (* set of instructions for coding another region *)
      end;

(* end module encode.type.param version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.type *)
(* vectors are useful for many applications, including the encoding of
   sequences.  this vector type is designed to be flexible enough to be
   used whenever one needs a vector.  it is designed as a linked list so
   it can contain as many elements as are ever needed.  the actual 'vector'
   is a record containing the length and a pointer to the first 'vectorpart'.
   that vectorpart is a record containing an array of the first 'maxvecpart'
   elements  (maxvecpart is a constant, from module vector.const, which must
   be included) and a pointer to the next 'vectorpart'.  the elements are
   of type real so that either integers or reals may be used.  *)
      partptr = ^vectorpart;
      vectorpart = record
         numbers : array[1..maxvecpart] of real;
         next : partptr;
      end;

      vector = record
         length : integer;
         part : partptr;
      end;
(* end module vector.type version = 'matmod 1.95 85 apr 18 tds/gds'; *)

var
      inst,  (* the instructions that generated the book, for alignment of
               the sequences *)
      book,  (* the book of sequences *)
      encseq,  (* the encoded sequences, which are arrays of integers *)
      encodep: text; (* the file of parameters which specifies how to
                     encode the sequences *)
      regions, (* the number of independently coded regions *)
      length,  (* the length of the sequence pointed to *)
      alignedbase: integer;  (* base in each sequence to be aligned by *)
      apiece: pieceptr; (* a pointer the the piece currently being encoded *)
      firstparam: paramptr; (* the beginning of the list of parameters *)
      fourpowers: array[0..codingmax] of integer; (* for storing powers of 4 *)
      seqvector: vector;  (* the encoded sequence vector *)
      noinst: boolean; (* true if the inst file is empty *)
   theline: integer; (* location in book *)

      alignmenttype: char; (* 'f' means alignment by First internal coordinate
         base, 'b' means alignment by Book, 'i' means alignment by Instructions *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* ************************************************************************ *)
(* end module package.primitive version = 7.72; {of delmod.p 2007 Jul 23} *)

(* begin module vector.functions *)
(* *********************************************************************** *)
(* begin module vector.io *)
(* *********************************************************************** *)
(* begin module vector.primitives *)
(* these functions and procedures do manipulations on vectors. *)

function vget(v:vector; pos:integer): real;
(* this returns from vector 'v' the value of the element at position 'pos' *)
var i: integer;
begin
      if (pos > v.length) or (pos < 1) then begin
         writeln(output,' error in call to function vget:',
            ' position ',pos:1,' is beyond the end of the vector');
         halt;
      end;

      (* move to the correct 'vectorpart' *)
      for i := 1 to ((pos -1) div maxvecpart) do  v.part := v.part^.next;

      (* get the proper array element from this part *)
      vget := v.part^.numbers[(pos -1) mod maxvecpart + 1];
end;

procedure vput(var v:vector; pos:integer; number: real);
(* this puts into vector 'v' the value of 'number' at position 'pos' *)
var i: integer;
      firstpart: partptr;  (* of the vector *)
begin
      if (pos > v.length) or (pos < 1) then begin
         writeln(output,' error in call to function vput:',
            ' position ',pos:1,' is beyond the end of the vector');
         halt;
      end;

      (* move to the correct 'vectorpart' *)
      firstpart := v.part;
      for i := 1 to ((pos -1) div maxvecpart) do  v.part := v.part^.next;

      (* put the 'number' into the proper array element for this part *)
      v.part^.numbers[(pos -1) mod maxvecpart + 1] := number;
      v.part := firstpart;
end;

procedure makevector(var v: vector; l: integer);
(* create the linked list of vector-parts which are needed for a vector
   of length 'l' *)
var   numparts,  (* number of parts needed for the vector *)
      i: integer; (* index to the parts of the vector *)
      firstpart,  (* of the vector *)
      newpart: partptr;  (* of the vector *)
begin
      if l < 1 then begin
         writeln(output,' makevector: positive length required');
         halt
      end;

      v.length := l;
      new(v.part);
      firstpart := v.part;
      numparts := (v.length - 1) div  maxvecpart + 1;
      for i := 1 to (numparts - 1) do begin
         new(newpart);
         v.part^.next := newpart;
         v.part := newpart;
      end;
      v.part^.next := nil;
      v.part := firstpart;
end;
(* end module vector.primitives version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.readvector *)
procedure readvector(var thefile:text; var v: vector);
(* read the elements of v into v from 'thefile'.  v must already be set
   up (all the parts created and linked together) before calling.
   'thefile' must contain, from the cursor point until the end of the
    vector, only numbers, either integers or reals; otherwise it will
    bomb.  the vector ends with an end of line so that
    end of file can be tested for. *)
var i,j: integer;  (* indecies *)
      numparts,  (* number of parts in the vector *)
      lastpart: integer; (* the number of elements in the last vector part *)
      firstpart: partptr;  (* pointer to the first vector part *)
begin
      firstpart := v.part;
      numparts := (v.length - 1) div maxvecpart + 1;
      lastpart := (v.length - 1) mod maxvecpart + 1;
      for i := 1 to (numparts - 1) do begin
         for j := 1 to maxvecpart do read(thefile,v.part^.numbers[j]);
         v.part := v.part^.next;
      end;
      for j := 1 to lastpart do read(thefile,v.part^.numbers[j]);
      readln(thefile);
      v.part := firstpart;
end;
(* end module vector.readvector version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.writevector *)
procedure writevector(var thefile: text; v: vector; y,z: integer);
(* writes the elements of 'v' to 'thefile'.  the integers y and z are:
    y: the size of the field for printing the number;
    z: the size of the decimal part of the field,
    if z = 0 then integers will be printed. *)
var   pos,   (* posititon in the vector *)
      i,j: integer;  (* indecies *)
      numparts,  (* number of parts in the vector *)
      lastpart: integer; (* the number of elements in the last vector part *)
      firstpart: partptr;  (* pointer to the first vector part *)
      x: integer; (* the number of elements to write on a line *)
begin
      x := trunc(vpagewidth/(y+1));
      pos := 0;
      firstpart := v.part;
      numparts := (v.length - 1) div maxvecpart + 1;
      lastpart := (v.length - 1) mod maxvecpart + 1;
      if (z = 0) then begin
         for i := 1 to (numparts - 1) do begin
            for j := 1 to maxvecpart do begin
               write(thefile,' ',round(v.part^.numbers[j]):y);
               pos := succ(pos);
               if (pos mod x = 0) then writeln(thefile);
            end;
            v.part := v.part^.next;
         end;
         for j := 1 to lastpart do begin
            write(thefile,' ',round(v.part^.numbers[j]):y);
            pos := succ(pos);
            if ((pos < v.length) and (pos mod x = 0)) then writeln(thefile);
         end
      end
      else begin
         for i := 1 to (numparts - 1) do begin
            for j := 1 to maxvecpart do begin
               write(thefile,' ',v.part^.numbers[j]:y:z);
               pos := succ(pos);
               if (pos mod x = 0) then writeln(thefile);
            end;
            v.part := v.part^.next;
         end;
         for j := 1 to lastpart do begin
            write(thefile,' ',v.part^.numbers[j]:y:z);
            pos := succ(pos);
            if ((pos < v.length) and (pos mod x = 0)) then writeln(thefile);
         end;
      end;
      writeln(thefile);
      v.part := firstpart;
end;
(* end module vector.writevector version = 'matmod 1.95 85 apr 18 tds/gds'; *)
(* *********************************************************************** *)
(* end module vector.io version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.dotproduct *)
function dotproduct(vectora,vectorb: vector): real;
(* this returns the dotproduct of 'vectora' and 'vectorb' *)
var i: integer;  (* an index *)
    j: real;  (* partial products *)
begin
      if vectora.length <> vectorb.length then begin
         writeln(output,' function dotproduct: vector lengths must',
            ' be equal');
         halt
      end;
      j := 0;
      for i := 1 to vectora.length do
         j := j + vget(vectora,i) * vget(vectorb,i);
      dotproduct := j;
end;
(* end module vector.dotproduct version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.magnitude *)
function magnitude(var v: vector): real;
(* find the magnitude (length) of the vector v *)
begin
      magnitude := sqrt(dotproduct(v,v));
end;
(* end module vector.magnitude version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.normalize *)
procedure normalize(var v: vector);
(* this replaces the vector v with the congruent vector of unit length,
   i.e.,  the resulting vector v has the property that v.v = 1   *)
var i: integer;  (* index *)
      length: real;  (* of the unnormalized vector *)
begin
      length := magnitude(v);
      for i := 1 to v.length do
         vput(v,i,vget(v,i)/length);
end;
(* end module vector.normalize version = 'matmod 1.95 85 apr 18 tds/gds'; *)
(* *********************************************************************** *)
(* end module vector.functions version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module encode.vectorsize *)
function vectorsize(param: paramptr): integer;
(* determine the size of the vector generated by a string of parameter
   records, begin with 'param' *)
var   size: integer;  (* of the vector, so far *)
begin
      size := 0;
      while (param <> nil) do begin
         size := size + param^.pvlength;
         param := param^.next;
      end;
      vectorsize := size;
end;

(* end module encode.vectorsize version = 'matmod 1.65 85 apr 9 tds/gds'; *)

(* 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 writeline(var afile: text; l: lineptr;
                        carriagereturn: boolean);
(* write a line to a file, with carriage return if
carriagereturn is true. *)
var
   index: integer; (* index to characters in l *)
begin
   with l^ do begin
      for index := 1 to length
      do write(afile, letters[index]);
   end;
   if carriagereturn then writeln(afile);
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);
{
This is illegal according to gpc because one cannot
write a pointer to a text file.  It can be unearthed
for debugging.
      write(output, ', pointer id: ',l:1);
}
      writeln(output);
      l := l^.next
   end;
end;

procedure cleardna(var l: dnaptr);
(* clear the dna strutures to the free list *)
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 *)
(* Note: the dirhomologous and dircomplement are treated as
plus and minus directions, which MIGHT NOT BE RIGHT! *)
var   i: integer; (* an intermediate value *)
begin
      with pie^.key do begin
         case piedir of
            dirhomologous,
            plus:  if p>=piebeg
                   then i:=p-piebeg+1
                   else i:=(p-coobeg)+(cooend-piebeg)+2;
            dircomplement,
            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 *)
(* Note: the dirhomologous and dircomplement are treated as
plus and minus directions, which MIGHT NOT BE RIGHT! *)
var   p: integer; (* an intermediate value *)
begin
      with pie^.key do begin
         case piedir of
            dirhomologous,
            plus: begin
               p:=piebeg+(i-1);
               if p>cooend
                  then if coocon=circular
                     then p:=p-(cooend-coobeg+1)
            end;
            dircomplement,
            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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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;
      (* protection added 2005 Sep 15 *)
      while (not eoln(thefile)) and (i < linelength)  do begin
         i:=succ(i);
         read(thefile,acharacter);
         l^.letters[i]:=acharacter
      end;
      (* protection added 2005 Sep 15 *)
      if not eoln(thefile) then begin
         writeln(output,'***********************************************');
         writeln(output,'* WARNING: brline: book line length exceeded');
         writeln(output,'* linelength > ',linelength:1,' characters');
         writeln(output,'* Only ',linelength:1,' characters read from book');
         writeln(output,'***********************************************');
      end;
      l^.length:=i;
      l^.next:=nil;
      readln(thefile);
      theline := succ(theline)
end;
(* end module book.brline version = 7.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* ************************************************************************ *)
(* end module package.brpiece version = 7.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* ************************************************************************ *)
(* end module package.getpiece version = 7.72; {of delmod.p 2007 Jul 23} *)

(* 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.18; (@ of prgmod.p 1996 September 12 *)

(* 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.18; (@ of prgmod.p 1996 September 12 *)

(* 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.18; (@ of prgmod.p 1996 September 12 *)

(* 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.18; (@ of prgmod.p 1996 September 12 *)

(* 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.18; (@ of prgmod.p 1996 September 12 *)

(* 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.18; (@ of prgmod.p 1996 September 12 *)

(* ************************************************************************ *)
(* end module package.trigger version = 4.18; (@ of prgmod.p 1996 September 12 *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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.72; {of delmod.p 2007 Jul 23} *)

(* 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 asterisks) 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. *)

   gettrigger: trigger; (* trigger to find 'get' *)
   defaulttrigger: trigger; (* trigger to find 'default' *)
   nametrigger: trigger; (* trigger to find 'name' *)
   piecetrigger: trigger; (* trigger to find 'piece' *)
   settrigger: trigger; (* trigger to find 'set' *)

   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: '}' *)

   quote1trigger: trigger; (* trigger to find single quote ' *)
   quote2trigger: trigger; (* trigger to find double quote " *)

   dotteddone: boolean; (* a dot '.' has been found in the name - ignore
      the rest of the name - for comparisons with mutations. *)

{
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;

procedure skipquote(quote: trigger);
(* skip an entire quote of either the ' or " persuasion *)
var
   kind: char; (* the kind of quote, ' or " *)
begin
   kind := quote.seek.letters[1];
{writeln(output,'skipquote ',kind);}
   repeat
      findnonblank(inst,ch); (* get to the quote *)
   until (ch = kind) or eof(inst);
   if ch <> kind then begin
      writeln(output,'end of quote starting with ',kind,' not found');
      halt;
   end;
end;

begin
 filltrigger(defaulttrigger,'default');
     filltrigger(gettrigger,'get    ');
    filltrigger(nametrigger,'name   ');
   filltrigger(piecetrigger,'piece  ');
     filltrigger(settrigger,'set    ');
   filltrigger(begincomment,'(*     ');
     filltrigger(endcomment,'*)     ');
     filltrigger(begincurly,'{      ');
       filltrigger(endcurly,'}      ');
  filltrigger(quote1trigger,'''      ');
  filltrigger(quote2trigger,'"      ');

  resettrigger(defaulttrigger);
  resettrigger(gettrigger);
  resettrigger(nametrigger);
  resettrigger(piecetrigger);
  resettrigger(settrigger);

  resettrigger(begincomment);
  resettrigger(begincurly);
  resettrigger(quote1trigger);
  resettrigger(quote2trigger);

   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 (* we are not inside a comment *)

                     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, quote1trigger);
                     if quote1trigger.found then begin
                        skipquote(quote1trigger);
                     end;

                     testfortrigger(ch, quote2trigger);
                     if quote2trigger.found then begin
                        skipquote(quote2trigger);
                     end;

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

                     testfortrigger(ch, settrigger);
                     if settrigger.found then begin
                        indefault := true;
                        resettrigger(settrigger)
                     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
}
                           (* 2007 Jun 22: replace loop with while
                              so that we can drop out when dotted
                              names are detected. *)
                           p := 1;
                           dotteddone := false;
                           while not dotteddone do begin
                              if eoln(inst)
                              then dotteddone := true
                              else begin
                                 read(inst,ch);
                                 (* ignore names after a dot *)
{
if ch = '.' then writeln(output,'inst dotteddone');
}
                                 if ch         = '.' then dotteddone := true;
                                 if letters[p] = '.' then dotteddone := true;
{
if ch = '.' then writeln(output,'book dotteddone');
writeln(output,'BUBBa ch = ',ch,'     ',p:1);
}
{zzz}
                                 if (letters[p] <> ch) and (not dotteddone) and
                                    (ch <> ';')
                                 then begin
                                    writeln(output,
                                           'The piece name in the book: ');
                                    writeln(output,letters:length);
                                    writeln(output,'does not match',
                                                 ' the inst file piece name:');
                                    (* write the letters that matched: *)
                                    for p1 := 1 to p-1 
                                        do write(output,letters[p1]);
                                    (* 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;
                                 p := p + 1;
                                 if p > length then dotteddone := true;
                              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.72; {of delmod.p 2007 Jul 23} *)

(* begin module book.getbase *)
function getbase(position: integer; pie: pieceptr):base;
(* Get a base from the position (internal coordinates) of the
   piece.  Protection is made against positions outside the piece.
In the case of circles it would be convenient to wrap around when
requests are off the end.  So the routine will do a modular
wrap for positions outside the range 1 to the length.
This is a new feature as of 2000 March 22.
*)
var
   workdna: dnaptr; (* pointer to the dna part of pie *)
   p: integer; (* current count of bases into the workdna *)
   spot: integer; (* the last base of the dna part *)
   thelength: integer; (* the length of the piece *)
begin
{
writeln(output,'NEW getbase: position=',position:1,'^^^^^^^^^^^^^^^^^^^^');
}
   (* handle cases of position out of range by circular wrapping *)
   thelength := piecelength(pie);
   while position < 1         do position := position + thelength;
   while position > thelength do position := position - thelength;

   workdna:=pie^.dna;
   p:=workdna^.length;
   while position > p do begin
{
writeln(output,' workdna^.length=',workdna^.length:1);
}
      workdna := workdna^.next;
      if workdna = nil then begin
        writeln(output,'error in function getbase!');
        halt
      end;
      p := p + workdna^.length;
   end;
{
writeln(output,'p=',p:1);
}
   if workdna = nil then begin
      writeln(output,'error in getbase: request off end of piece');
      halt
   end
   else begin
      spot := workdna^.length - (p-position);
{
writeln(output,'spot=',spot:1);
showdnasegment(output,workdna, spot);
}
      if (spot <= 0) then begin
         writeln(output,'error in getbase, spot (= ',spot:1, ') must be positive');
         halt
      end;
      if (spot > workdna^.length) then begin
         writeln(output,'error in getbase, spot (=',spot:1,
                        ') must be less than length (=',workdna^.length:1,')');
         halt
      end;
{
writeln(output,'base = ', workdna^.part[spot]);
}
      getbase:=workdna^.part[spot]
   end
end;
(* end module book.getbase version = 7.72; {of delmod.p 2007 Jul 23} *)

procedure initialize;
var i: integer;  (* an index *)
begin
      writeln(output,' encode, version ',version:4:2);
      reset(inst);
      if eof(inst) then noinst := true
                   else noinst := false;
      brinit(book, theline);
      rewrite(encseq);
      new(firstparam);
      new(apiece);
      fourpowers[0] := 1;  (*  4**0 = 1  *)
      for i := 1 to codingmax do fourpowers[i] := 4 * fourpowers[i-1];
end;

procedure setparam;
(* get the parameter list from the file encodep and put the information
   into the linked parameter list *)
var   newparam,  (* the new parameter list from 'encodep' *)
      param: paramptr;  (* a parameter record being read in *)
      firstspaces,  (* in the spaces list *)
      newspaces: spaceptr; (* for the spaces list *)
      i: integer;  (* an index *)
      ch: char;  (* to check for ':' *)

begin
      reset(encodep);
      if not eof(encodep) then begin
         (* read the alignment type *)
         readln(encodep, alignmenttype);
         if not (alignmenttype in ['f','i','b']) then begin
            writeln(output,'alignment type inst must be f, b, or i');
            halt
         end;

         param := firstparam;
         regions := 0;
         repeat
            regions := succ(regions);
            with param^ do begin
               readln(encodep,range[start],range[stop]);
               if (range[start] > range[stop]) then begin
                  writeln(output,' error in parameter file:',
                     ' end of range cannot be less than beginning of range');
                  halt;
               end;
               readln(encodep,window);
               if (window <= 0) then begin
                  writeln(output,' error in parameter file:',
                     ' window size must be positive');
                  halt;
               end;
               readln(encodep,wshift);
               if (wshift <= 0) then begin
                  writeln(output,' error in parameter file:',
                     ' window shift must be positive');
                  halt;
               end;
               read(encodep,coding);
               if (coding <= 0) then begin
                  writeln(output,' error in parameter file:',
                     ' coding must be positive');
                  halt;
               end;
               if coding > codingmax then begin
                  writeln(output,' error in parameter file:',
                     ' coding can not be greater than ',codingmax:1);
                  halt
               end;
               if (coding > codingmax) then begin
                  writeln(output,' error in coding file:',
                     ' requested coding level too large');
                  halt;
               end;
               new(spaces);
               ch := ' ';
               while (not eoln(encodep)) and (ch = ' ') do read(encodep,ch);
               if (ch = ':') then begin
                  new(firstspaces);
                  spaces := firstspaces;
                  read(encodep,spaces^.skips);
                  for i := 2 to coding -1 do begin
                     new(newspaces);
                     spaces^.next := newspaces;
                     spaces := newspaces;
                     read(encodep,spaces^.skips);
                  end;
                  spaces^.next := nil;
                  spaces := firstspaces
               end
               else spaces := nil;
               readln(encodep);
               readln(encodep,cshift);
               if (cshift <= 0) then begin
                  writeln(output,' error in parameter file:',
                     ' coding shift must be positive');
                  halt;
               end;
                  (* wvlength is 4**coding *)
               wvlength := round(exp(coding*ln(4)));
            (* pvlength is wvlength * (the number of windows in the range) *)
               pvlength := wvlength *
                        trunc((range[stop] - range[start])/wshift + 1);
            end;
            if eof(encodep) then param^.next := nil
            else begin
               new(newparam);
               param^.next := newparam;
               param := newparam;
            end;
         until eof(encodep)
      end
      else begin
         writeln(output,' error: empty parameter file (encodep)');
         halt;
      end;
end;

procedure encheader;
(* write the header information to encseq *)
var
    param: paramptr;  (* points to a parameter record *)
    aspace: spaceptr;  (* in the spaces list *)
    i: integer;  (* an index *)

begin
      writeln(encseq,' encode ',version:4:2,'; encoding of sequences in');
      write(encseq,' ');
      copyaline(book,encseq);
      writeln(encseq);
      writeln(encseq,' ',regions:1,' independently coded regions');
      writeln(encseq);

      param := firstparam;
      while param <> nil do with param^ do begin
         write(encseq,' ',range[start]:1,' to ',range[stop]:1);
         writeln(encseq,' is encoded as:');
         writeln(encseq,' ':6,window:1,' long windows');
         writeln(encseq,' ':6,wshift:1,' bases shift to new window');
         write(encseq,' ':6,coding:1);
         if (coding > 1) then begin
            write(encseq,' :');
            if (spaces = nil) then for i := 1 to coding - 1 do
               write(encseq,' 0')
            else begin
               aspace := spaces;
               for i := 1 to coding - 1 do begin
                  write(encseq,' ',aspace^.skips:1);
                  aspace := aspace^.next;
               end;
            end;
         end;
         writeln(encseq,' are the coding units');
         writeln(encseq,' ':6,cshift:1,' bases shift to new coding site');
         param := param^.next;
         writeln(encseq);
      end;
      writeln(encseq,' ',vectorsize(firstparam):1,' is the vector length');
      writeln(encseq);
end;

procedure codeit(apiece: pieceptr; alignedbase: integer; {aparam: paramptr;}
                     var v: vector);
(* this procedure takes one sequence (which has been aligned) and encodes
   it into a string of integers according to the parameters the user has
   specified.  each window of the sequence causes a vector to be printed,
   the size of which is determined by the coding level specified.  for
   example, if the coding level is 1 (monos are being counted) then the
   vector has four elements.  if the coding level is 2 (dis are being
   counted) then the vector has 16 elements.  the vectors for each window
   are concatenated to give the vector for each sequence, which is ended
   with the 'end of sequence' symbol.
   the procedure takes advantage of the fact that the type 'base' is an
   ordered set of the nucleotides, with a,c,g,t being assigned internally
   the values 0,1,2,3.  the total number of elements of a vector is
   4**coding-level, and these correspond to the number of oligos of each
   type from all a's to all t's.  this makes it easy to increase the
   vector element for the oligo which exists at the coding site.  for
   example, there are 64 tri-nucleotides and the vector which encodes
   them has elements 0 (for the oligo aaa) through 63 (for the oligo ttt).
   the number of the oligo cgt would be stored in element 27, as seen from
       16 * 1   (16 for 4**(coding-level - 1) and 1 for the c)
     +  4 * 2   ( 4 for 4**(coding-level - 2) and 2 for the g)
     +  1 * 3   ( 1 for 4**(coding-level - 3) and 3 for the t)
     ---------
         27
   another example is finding the position of the oligo taac in the vector
   of all the 256 tetramers, from position 0 to 255.
       64 * 3   (64 = 4**(coding-level - 1) and 3 for the t)
    +  16 * 0   (16 = 4**(coding-level - 2) and 0 for the a)
    +   4 * 0   ( 4 = 4**(coding-level - 3) and 0 for the a)
    +   1 * 1   ( 1 = 4**(coding-level - 4) and 1 for the c)
    ---------
        193    ( is the vector element corresponding to the oligo taac )
*)

var   param: paramptr;  (* list of parameters for encoding *)
      aspace: spaceptr;  (* in the spaces list *)
      startsite,  (* of a window on the sequence *)
      firstpos,  (* of a coding site within a window *)
      pos,  (* of a particular base to be encoded *)
      sitesize,  (* total number of bases in a coding site, both those that
                  are encoded and those that are skipped *)
      element,  (* of the encoding vector, from 1 to v.length *)
      welement,  (* of the window vector, from 0 to 4**coding - 1 *)
      length,  (* of the sequence *)
      i: integer;  (* an index *)

begin
      length := piecelength(apiece);
      param := firstparam;
                  (* clear the vector *)
      for i := 1 to v.length do vput(v,i,0);
         (* start at the beginning fo the vector *)
      element := 1;

      repeat  (* for each parameter record in the list *)
         with param^ do begin
               (* determine the number of bases in a coding site *)
            aspace := spaces;
            sitesize := coding;
            while (aspace <> nil) do begin
               sitesize := sitesize + aspace^.skips;
               aspace := aspace^.next;
            end;

               (* set the window beginning at the range start *)
            startsite := alignedbase + range[start];
            repeat  (* move the window across the range *)
                  (* set coding site beginning at window beginning *)
               firstpos := startsite;
               repeat
                        (* if the whole coding site is on the piece *)
                  if (firstpos > 0) then
                  if (firstpos + sitesize -1 <= length) or
                     ((apiece^.key.piecon = circular) and (firstpos <= length))
                  then begin
                     pos := firstpos;
                     welement := 0;
                     aspace := spaces;
                     for i := coding downto 1 do begin
                        welement := welement + fourpowers[i-1] *
                           ord(getbase(pos,apiece));
                        if (aspace <> nil) then begin
                           pos := pos +1 + aspace^.skips;
                           aspace := aspace^.next
                        end
                        else pos := pos +1;
                           (* for circular sequences *)
                        if (pos > length) then pos := pos - length;
                     end;
                  (* increment the vector at position 'element + welement' *)
                     vput(v,element + welement,vget(v,element + welement) + 1);

                  end;
                  firstpos := firstpos + cshift;
               until (firstpos >= startsite + window);

               startsite := startsite + wshift;
                  (* advance element beyond window *)
               element := element + wvlength;
            until (startsite > alignedbase + range[stop]);
         end;
         param := param^.next;
      until (param = nil);
end;  (* codeit *)

begin  (* encode *)
      initialize;
      setparam;
      encheader;
      seqvector.length := vectorsize(firstparam);
      makevector(seqvector,seqvector.length);
      if noinst
      then if alignmenttype = 'i' then begin
         alignmenttype := 'b';
         writeln(output,'There are no instructions so alignment type',
                        ' is forced to b mode');
      end;

      while not eof(book) do begin
         if alignmenttype <> 'i' then begin
            getpiece(book,theline,apiece); (* read in the piece *)
            if not eof(book) then begin
               length := piecelength(apiece);
               case alignmenttype of
                  'i': begin writeln(output,'prgm error'); halt end;
                  'b': begin
                         alignedbase := pietoint(0, apiece);
                  end;
                  'f': begin
                         alignedbase := 1;
                  end;
               end;
            end;
         end
         else align(inst,book,theline,apiece,length,alignedbase);
         if not eof(book) then begin
            codeit(apiece,alignedbase,{firstparam,}seqvector);
            writevector(encseq,seqvector,3,0);
            clearpiece(apiece);
         end;
      end;
1: end. (* encode *)
