program ri(inst, book, rsdata, values, rip, wave,
           rixyin, sequ, ribl, riinst, output);
(* ri: Rindividual is calculated for every site in the aligned book

  Thomas D. Schneider, Ph.D.
  Frederick National Laboratory for Cancer Research
  toms@alum.mit.edu
  https://alum.mit.edu/www/toms

module libraries required: delman, prgmod, delmod *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 2.83; (* of ri.p 2013 Jun 06
2013 Jun 06, 2.82: min/maxribl increased
2012 Jan 18, 2.81: maxribl increased
2012 Jan 18, 2.80: maximumrange increased by getting new delmod modules
2008 Apr 24, 2.79: reported value of neginfinity should have same wid/dec as table
2008 Apr 20, 2.78: niot must be real for precise computations
2005 Sep 28, 2.77: protect ln(x) from x <= 0.
2005 Sep 28, 2.76: lock unused modules
2005 May 23, 2.75: add patent message
2004 Dec 15, 2.74: update documentation of new parameters
2004 Dec 15, 2.73: prgmod module readheader upgraded.  cleanup.
2004 Dec 15, 2.72: crash on reading hold file using gpc compiler
2003 Aug 18, 2.71: niot is used when nxl < 1.
2003 Aug 18, 2.70: niot should accept -1000.
                   Ricalc: in this case compute appropriately.
2002 Aug 19, 2.69: add standard deviation computation to rixyin output EXPERIMENTAL
2002 Jun 21, 2.68: reported consensus range was wrong - corrected in prgmod.
2002 Apr 19, 2.67: documentation upgrade
2001 May 10, 2.65: BUG in Ri(bl) found and corrected!
2000 Dec 13, 2.64: upgrade alignment from delmod
2000 Dec 13, 2.63: upgrade See Also documentation
2000 Jul 30, 2.62: upgrade to GPC
2000 Jun 19, 2.59: fix symmetry detection in matrixsymmetry
1999 Oct  6, 2.55: now creates riinst file
1999 Jul  9, 2.53: automatic upgrade from rip 2.33 to 2.53
1999 Jul  5, 2.49: name and bounds in rip, standard routines from prgmod.p
1999 Jun  1, 2.47: upgrade to theline stuff for aligned list
1999 Mar 11, 2.44: attach wave definition to end of ribl 
1998 Oct 19, ----: range set to 1000
1997 Apr 17, 2.33: ribl now contains n(b,l) to allow detection of
    negative infinities (zero cases) by later programs (and whatever).
1995 Aug 29:  Standard deviation should be calculated from n-1.
origin 1990 Aug 15 from Dalvec 1.05 *)
updateversion = 2.53; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.ri *)
(*
name
   ri: Rindividual is calculated for every site in the aligned book

synopsis
   ri(inst: in, book: in, rsdata: in, values: in, rip: in, wave: in,
      rixyin: out, sequ: out, ribl: out, riinst: out, output: out)

files
   inst: delila instructions of the form 'get from 56 -5 to 56 +10;'
      (This file may be empty, in which case the sequences will be
      aligned by their 5' ends.)

   book: the book generated by delila using inst

   rsdata: data file from rseq program

   values: a file containing the values of the objects to which the Ri
      values are to be compared.  The file may be empty.  The file
      can contain some values and then terminate.  Further values
      will be assumed to be zero.  This allows one to copy a rixyin file
      produced by ri in to the values file.  Then one adds a new site
      runs delila etc and when ri runs again it will get the old values
      lined up to the new values, with zero for the one that was
      not there previously.

   rip:  Parameters to control the program.  The file must contain the
      following parameters, one per line:

      * parameterversion:  The version number of the program.  This allows
        the user to be warned if an old parameter file is used.

      * name (string in quotes): The name of the weight matrix.

      * thefrom, theto (integers): The FROM and TO over which to do the Ri
        calculation.  These must not exceed either of those in the inst/book
        or the rsdata.  The FROM value must be less than or equal to the TO
        value.

      * column (integer): defines the column of the values file to use.

      * lowerRi, upperRi (two reals): the lowest and highest Ri evaluation to
        report to rixyin, sequ and riinst.  If the first character of the
        line is 'a' then all evaluations are reported.  Otherwise two real
        numbers are expected.  Sequences within this range are printed to
        rixyin, sequ and riinst depending also on the fifth parameter.

      * lowerValue,upperValue (two reals): the lowest and highest Value
        evaluation to report to rixyin, sequ and riinst.  If the first
        character of the line is 'a' then all evaluations are reported.
        Otherwise two real numbers are expected.  Sequences within this range
        are printed to rixyin, sequ and riinst depending also on the fifth
        parameter.

      * printsequ (character):  determines whether or not to produce any raw
        sequences in the sequ file.  If the first character of the line is
        'p', sequences selected according to the third and fouth parameters
        are printed to sequ file.  (This is a complete on-off switch for the
        sequ file.)

      * printrixyin (character):  determines whether or not to print the
        sequence of the site being analyzed.  If the first character is 'p'
        then the sequence is printed to the rixyin file.

      * partials (character): determines whether or not to print sequences
        which have a partial site.  The problem is that if there is part of a
        site, then the Ri value is questionable, depending on where the
        deletion was.  The best analysis would not use a partial site, as it
        messes up the statistics.  If the first character is:

         n  Don't print the line at all.
         i  Keep the line, but force the Ri value to be -infinity.
            This allows the lines of rixyin to be correlated to the values
            still.
         -  (any other character): print as it is.

      * niot (character (one of 'slnb') and/or real): negative infinity or
         t:  determines what to do when f(b,l) = 0.  Positions for which
         f(b,l) = 0 will have negative infinity in the Ri(b,l) table.  The
         letter 's' means to use an "extended" Rodger Staden's method of
         giving f(b,l) = 1/(n+t), where t is a non-negative integer following
         the 's'.  When t = 0, this gives Staden's original method.  Using
         t = 1 would often give frequencies that are too low.  The value of
         t=2 corresponds to the Laplace Law of Succession.

         The ``law of succession of Laplace''
         \cite{Feller1968.Laplace,Papoulis1990} states that given $n$ trials
         in which there were $k$ results of one kind, the best estimate for
         the probability in another trial is $(k+1)/(n+2)$.  In the present
         case, we need the probability of the absence of a particular base
         when searching for {\em another\/} binding site, so $k = 0$ and the
         best estimate is $1/(n+2)$.  For this reason we set $t=2$ for most
         purposes.

         If instead of 's' for Staden's extended method, 'l' is used, the
         program will use the full Laplace Law of sucession

             f(b,l) = (k+1)/(n+2) = (n(b,l)+1)/(n(l)+2)

         as the estimate for the probability for k items out of n total.
         This OVERESTIMATES the information content and so is NOT to be used
         in general.  (This can be easily demonstrated by looking at the
         effect on the average of the distribution.)

         If the character is 'b' then Staden's method is used as above but
         Ehnb is not added.  This is ONLY to be used to simulate the bug
         corrected on 2001 May 10!  (See Technical Notes below.)

         If there is no 's', 'l' or 'b', (or when the line begins with 'n')
         then the program expects a number which the value for negative
         infinity.  It should be a value sufficiently below zero so that
         sites that are being excluded from the definition according to
         f(b,l) are separated from the true sites.  -1000 is a useful value,
         as it will always displace sites with exceptions far away from
         zero.

      * alignmenttype (character): 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.

       * Ribound (real): the Ri boundary for this definition, bits.
            The lower bound on Ri to find in the book.
            Use 0 to represent the Second Law of Thermodynamic bound.
            Use -500 or lower to locate all sites. 

       * Zbound (real): the Z boundary for this definition.
            The lower bound on Z (standard deviation) to find in the
            book.
            Use 100 to locate all sites.

       * Pbound (real): the probability boundary for this definition.
            The upper probability P to find in the book.
            Use 1.00 to locate all sites.

       Ribound, Zbound and Pbound are used in conjunction, so 0, 100,
       1.00 respectively would locate all sites greater than zero
       bits, with no other constraint.  Most of the time we adjust
       Ribound and do not play with the others.

   wave: Define a cosine wave.  See makelogo for the definition of this
         file.  These data are simply copied to the end of the ribl file.
         This allows the lister program to create a wave under the walkers If
         you use the same wave definition for makelogo, the same wave in the
         logo will be on the walker.

   rixyin: input to the xyplo program.  The file contains these columns of data:
         1 piece number
         2 piece name
         3 sequence of region analyzed 
         4 length of region analyzed on this piece
         5 aligning coordinate on the piece
         6 Rindividual for the piece
         7 value from the values file (or 0 if values is empty)
         8 standard deviation of Rindividual for that sequence.

zzzqqq describe Ri SD computation here.  This computation is EXPERIMENTAL.
The statistics is NOT PROVEN and the program is not even guaranteed
to compute this unproven value correctly!

********************************************************************************
2002 Jan 24

risd = Ri standard deviation

the test directory contains a functioning routine
that works - in pascal.  It needs to be boiled down more.

2002 Jan 29

Here is the comptuation, at a position l with n(l) sequences
and f(b,l) = n(b,l)/n(l).
If the frequency is zero, then use

f := 1/(2+n(l))
otherwise use
f := f(b,l)

The (corrected to account for n(l) variation) formula given in the perl
program is

variance = ((1-f)/f)/n(l)

or

variance = ((1/f-1)/n(l)

for most cases where n(b,l) > 0

For simplicity let
b = n(b,l)
n = n(l)

variance = ((1/f(b,l)-1)/n(l)
         = ((1/(n(b,l)/n(l))-1)/n
         = ((1/(b/n))-1)/n
         = (1/b)- (1/n)

which is amusingly simple.
********************************************************************************

   sequ: the raw sequences reported to rixyin if any selection is made
      (fourth line of rip file).  These end in periods, so they can be
      given to makebk to create a book.

   ribl: weight matrix Ri(b,l).  The file begins with the name of the weight
      matrix.  It contains the information content for each base b at
      each position l, in bits.  Lines that start with * are notes.  The next
      line contains the matrix FROM and TO coordinates (thefrom and theto).

      This is followed by the matrix in the order A, C, G, T from FROM to TO.

      At the end of each line of the matrix is the position (integer)
      followed by the numbers of bases (A, C, G, T) at that position.

      After the matrix, real numbers on individual lines report:
         Ri mean (Rsequence of selected region)
         Ri standard deviation
         Ri of consensus sequence
         Ri of anticonsensus sequence
         Ri average for random (equiprobable) sequence
         the value of negative infinity
      These are all for the given range.
      (Note: Although the mean Ri for the sites is Rsequence, to get a good
      estimate of this, it is better to use the value calculated by the rseq
      program because that is less sensitive to missing sequence data.)

      Then one integer is given; n the number of sequences.

      Then, the symmetry of the matrix is defined by the first character
      on the final line:
          a: asymmetric
          o: odd symmetry (the zero base is in the center)
          e: even symmetry (the zero base is left of the center)

      The file contains any number of wave definitions, copied
      from the wave file.

      The file ends with a period, ".".

   riinst: Delila instructions for the selected subset of sites according to
      the parameters lowerRi, upperRi, lowerValue, upperValue.  The gets are
      given in this form:


piece K01789; get from 283 -FROM to same +TO direction +; (@ 4.555974 bits @)

      This allows the user to substitute a range for the FROM and TO tags.
      For example, using the unix sed:

         cat riinst | sed -e "s/FROM/50/" | sed -e "s/TO/50/" > newinst

      gives results like this:

piece K01789; get from 283 -50 to same +50 direction +; (@ 4.555974 bits @)

      See bugs.
      (In the actual instructions the '@' is replaced with *; it
      can't be in these examples because they are comments in
      the source code.)

   output: messages to the user

description
   The program determines the individual informations of the sites in the book
   as aligned by the instructions, according to the frequency table given in
   the rsdata file.  The program calculates the Ri(b,l) table:

       Ri(b,l) := 2 - (- log2( f(b,l)))

   and sums this up for each sequence.  Ri is defined so that the average of
   the Ri's for a set of sequences is Rsequence.  However, if the sequences are
   incomplete, the average will probably be less than Rsequence.  The rixyin
   output is ready to read into the xyplo program for plotting and linear
   regression.  The ribl matrix is ready to be used to scan sequences with the
   scan program.

   The program can be used in subtle ways.  For example, one can analyze the
   individual information of the left half of a binding site.  This result can
   then be used in the values file to compare against the analysis of the right
   side of a binding site.

examples

rip:

2.54          version of ri that this parameter file is designed for.
"Fis"         name of the ribl matrix
-10 +10       From-to range to do the evaluation
1             column of the values file to copy to xyin
a 0 1000      lowest to highest Ri to put in xyin and sequ (a = any)
a -1000 +1000 lowest to highest Value to put in xyin and sequ (a = any)
n             p means print sequence to the sequ file
p             p means print sequence to the xyin file
-             -: accept all sites; n: no partials; i: partials -> -infinity
s 2           s: use Staden's Method, f(b,l)=1/(n+t); else negative infinity
i             alignmenttype first, book, instructions
-5000         Ribound: real; the Ri boundary for this definition
100           Zbound:  real; the Z boundary for this definition
1.00          Pbound:  real; the probability boundary for this definition

documentation

@article{Schneider.Ri,
author = "T. D. Schneider",
title = "Information Content of Individual Genetic Sequences",
journal = "J. Theor. Biol.",
volume = "189",
number = "4",
pages = "427-441",
note = "\htmladdnormallink
{https://alum.mit.edu/www/toms/paper/ri/}
{https://alum.mit.edu/www/toms/paper/ri/}",
comment = "indiv.tex",
comment = "Submitted, April 1997",
year = "1997"}

@article{Schneider.walker,
author = "T. D. Schneider",
title = "Sequence Walkers:
a graphical method to display how binding proteins
interact with {DNA} or {RNA} sequences",
journal = "Nucl. Acids Res.",
volume = "25",
comment = "walker.tex, November 1, issue 21",
note = "\htmladdnormallink
{https://alum.mit.edu/www/toms/paper/walker/}
{https://alum.mit.edu/www/toms/paper/walker/},
erratum: NAR 26(4): 1135, 1998",
pages = "4408-4415",
year = "1997"}

@article{Staden1984,
author = "R. Staden",
title = "Computer methods to locate signals in nucleic acid sequences",
journal = "Nucl. Acids Res.",
volume = "12",
pages = "505-519",
year = "1984"}

@inproceedings{Feller1968.Laplace,
author = "William Feller",
booktitle = "An Introduction to Probability Theory and Its Applications",
volume = "I",
edition = "3rd",
publisher = "John Wiley \& Sons, Inc.",
address = "New York",
isbn = "0-471-25708-7",
pages = "123-124",
year = "1968"}

@book{Papoulis1990,
author = "Athanasios Papoulis",
title = "Probability \& Statistics",
publisher = "Prentice Hall",
address = "Englewood Cliffs, NJ",
comment = "Englewood Cliffs NJ 07632
QA273.P197 1990 / ISBN 0-13-711698-5.
The Laplace law of succession is on page 173",
year = "1990"}

see also
   {program that generates the rsdata file: } rseq.p

   {program to plot the rixyin data file:   } xyplo.p

   {program that uses the ribl file for searches: } scan.p

   {program that uses the ribl to create an interactive
    display of sequence walkers: } makewalker.p

   {program that uses the output of the scan program
    to display multiple sequence walkers and other features
    and marks along a sequence: } lister.p

   {program that should be run to get a sequence logo
    before attempting to use this program to make walkers: } makelogo.p

   {modules used internally to this program are in: } prgmod.p

   {For people who have the ri program installed, example files are in your
   ii/delila directory under under files beginning with 'fis'.}

author

   Thomas D. Schneider

bugs

   In riinst the organism and chromosome names are not written because they
   are not available by the alignment method of reading.  For now the
   instructions are given as:
       organism ORG; chromosome CHR;
   Of course this only will work for a single organism, but that is useful!

technical notes

   2001 May 10:  Weight Matrix Bug corrected

The weight function is Ri = Enbh+log2(f) where f is the frequency of a
base in an aligned position, Enbh ~ 2 (2 with small sample
correction).  If f is zero the function goes to negative infinity. 
Since that is too extreme - another sequence might contain the base -
in the Staden-like method we use an estimate for f in that case,
1/(n+2) which is the best estimate of the frequency given the small
sample size.

Example of the bug.  With n = 18, (and putting Enhb=2 for simplicity)
we have Ri = 2 + log2(1/20) = -2.32.  Oh no!  The weights are -4.32! 
I looked into the code (the routine is Ricalc) and found that there is
a bug!  I forgot to add the 'Enbh' part.

Versions of Ri up to 2.64 incorrectly used Ri = log2(1/(n+2)) instead
of the correct function when the number of bases is zero.  This means
that the scan program has been missing some sites because they were
weighted too low.  To allow for comparison to previous versions of the
code, the niot parameter can be set to 'b'.  This should not be used
normally.

*)
(* end module describe.ri *)

{
more see also stuff:
   inst, book, rsdata, values, rip, rixyin, sequ, ribl
}

(* begin module info.const *)
   infofield = 10; (* size of field for printing information in bits *)
   infodecim = 6; (* number of decimal places for printing information *)
       (* these are used for conlist only *)
   nfield = 8; (* size of field for printing n, the number of sites *)
(* end module info.const version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module ribl.const *)
   minribl = -10001; (* lowest  ribl matrix from allowed *)
   maxribl = +10001; (* highest ribl matrix to   allowed *)
   defnegativeinfinity = -1000; (* default for negative infinity
                                   for the Ri(b,l) table *)
(* end module ribl.const version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module interact.const *)
(* begin module string.const *)
      maxstring = 2000; (* the maximum string *)
(* end module string.const version = 5.67; (@ of prgmod.p 2017 Sep 29 *)
(* end module interact.const version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module Ri.filler.const *)
      fillermax = 21; (* the size of the filler array for a string *)
(* end module Ri.filler.const *)

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

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

(* end module my.book.const version = 'delmod 6.60 91 Jan 11 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.75; {of delmod.p 2013 Jun 06} *)

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 *)
      previous: stringptr; (* the previous string in a series *)
      next: stringptr; (* the next string in a series *)
   end;
(* end module string.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* 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 = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* 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 = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module rs.type *)
   rstype = record (* types of data in the data table from rseq *)
      rstart, rstop: integer; (* range of the data *)
      l, (* position *)
      nal,ncl,ngl,ntl: integer; (* numbers of each base *)
      rsl, (* rsequence(l) *)
      rs,  (* running sum of rs *)
      varhnb, (* variance estimate for small sample correction *)
      sumvar: (* sum of varhnb *)
         real;
      nl: integer; (* = a+c+g+t *)
      ehnb: real; (* hg-e(n) *)
      flag: char; (* e = exact, a = approximate sampling statistics *)
   end;
(* end module rs.type version = 4.75; (@ of rsgra.p 1990 Oct 2 *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

      (* 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.75; {of delmod.p 2013 Jun 06} *)

(* begin module rsgra.type *)
(* begin module wave.type *)
      waveptr = ^waveparam; (* to link several wave definitions together *)
      waveparam = record (* parameters to define a cosine wave *)
         (* define a cosine wave: *)
         extreme: char; (* h or l, the high or low extreme to be defined *)
         wavelocation: real; (* the location in bases of the extreme *)
         wavebit: real; (* the location in bits of the extreme *)
         waveamplitude: real; (* the amplitude of the wave in bits *)
         wavelength: real; (* the wave length of the wave in bases *)
{ 2000 Nov 1.  This is difficult because the length of a curve is not
simply multiplied by a scale factor.
         BASEdashon: real; (* dashon interval in bases *)
         BASEdashoff: real; (* dashon interval in bases *)
         BASEdashoffset: real; (* dashon interval in bases *)
}
         dashon: real; (* the size of on dashes in cm.
                          dashon <= 0 means no dashes *)
         dashoff: real; (* the size of on dashes in cm *)
         dashoffset: real; (* the size of off dashes in cm *)
         thickness: real; (* thickness of the wave in cm.  <= 0 means default *)
         next: waveptr; (* the next wave *)
      end;
(* end module wave.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)
(* end module rsgra.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)
 
(* begin module scan.type *)
   (* define the ribl matrix.  Source library is prgmod.p
      REQUIRES ribl.const module
      REQUIRES base.type or book.type module
      REQUIRES wave.type module for wave definition! *)
   rblarray = array[a..t, minribl..maxribl] of real; (* real(B,L) *)
   ribltype = record
      riblname: string; (* name of the weight matrix *)
      riblheader: stringptr; (* misc header information, starting with '*' *)
      data: rblarray; (* real(B,L) *)
      numbers: array[a..t, minribl..maxribl] of integer; (* n(b,l) *)
      frombase, tobase: integer; (* range of the matrix *)
      mean, stdev: real; (* mean and standard deviation of distribution *)
      consensus: real; (* Ri value of consensus sequence *)
      anticonsensus: real; (* Ri value of anticonsensus sequence *)
      averageRi: real; (* average Ri for random sequence *)
      n: integer; (* number of sequences used to create matrix *)
      symmetry: char; (* the symmetry of the Riblmatrix *)
      waves: waveptr; (* cosine wave definitions for the matrix *)
      cmperbase, cmperbit: real; (* required for wave dash definition *)
      Ribound: real; (* the Ri boundary for this definition *) 
      Zbound: real; (* the Z boundary for this definition *) 
      Pbound: real; (* the probability boundary for this definition *) 
   end;
   ribltypeptr = ^ribltype; (* pointer to ribltype so it does
      not need to be created unless actually used *)
(* end module scan.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module datetime.type *)
   (* array for dates *)
   datetimearray = packed array[1..datetimearraylength] of char;
(* end module datetime.type version = 7.75; {of delmod.p 2013 Jun 06} *)

var
(* begin module Ri.var *)
   inst, (* the delila instructions required by the align procedures *)
   book, (* the book to be aligned *)
   rsdata, (* output of rseq program *)
   values, (* values of objects to correlate Ri with *)
   rip, (* parameters to control the program *)
   wave, (* parameters to control the program *)
   rixyin, (* output of Ri, columns of Ri and length data *)
   sequ, (* raw sequences if selection is being done *)
   ribl, (* output of Ri, Ri(b,l) weight matrix *)
   riinst: (* delila instructions for the sequences *)
      text;
(* end module Ri.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.75; {of delmod.p 2013 Jun 06} *)

(* 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 = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

GPC Functions:
function  GetUnixTime (var MicroSecond : Integer) : UnixTimeType;

http://agnes.dida.physik.uni-essen.de/~gnu-pascal/gpc_109.html#SEC109

7.10.8 Date And Time Routines 

procedure GetTimeStamp (var t : TimeStamp); 
function Date (t : TimeStamp) : packed array [1 .. DateLength] of Char; 
function Time (t : TimeStamp) : packed array [1 .. TimeLength] of Char; 

DateLength and TimeLength are implementation dependent constants. 

GetTimeStamp (t) fills the record `t' with values. If they are valid, the Boolean
flags are set to True. 

TimeStamp is a predefined type in the Extended Pascal standard. It may be
extended in an implementation, and is indeed extended in GPC. For the full
definition of `TimeStamp', see section 8.255 TimeStamp. 

*)
var
    t: TimeStamp; 

(* begin module pluckdigit *)
function pluckdigit(number, logplace:integer): char;
(* return the digit at the place value ('logplace') position of number.
example:
      pluckdigit(13625, 3) = 3
      pluckdigit(13625, 4) = 1
This routine was taken from module numberdigit in prgmod.p, but is
modified so as not to give the sign.  Instead it gives zeros
above the digits.  'myabsolute' replaced 'absolute', which
is apparently a keyword for GPC.  The name is kept for
to keep the code looking similar to its origin.
*)
var
      place: integer; (* the exponent of logplace *)
      count: integer; (* used to make place *)
      myabsolute: integer; (* the absolute value of number *)
      acharacter: char; (* the character to be returned *)
procedure digit;
(* extract a digit at the place position *)
var
      tenplace: integer; (* ten times place *)
      z: integer; (* an intermediate value *)
      d: integer; (* the digit extracted *)
begin (* digit *)
      tenplace:=10*place;
      z:=myabsolute-((myabsolute div tenplace)*tenplace);
      if place = 1
         then d:=z
         else d:= z div place;
      case d of
         0: acharacter:='0';
         1: acharacter:='1';
         2: acharacter:='2';
         3: acharacter:='3';
         4: acharacter:='4';
         5: acharacter:='5';
         6: acharacter:='6';
         7: acharacter:='7';
         8: acharacter:='8';
         9: acharacter:='9';
      end
end; (* digit *)
begin (* pluckdigit *)
      place:=1;
      for count:=1 to logplace do place:=10*place;
      if number=0 then begin
         acharacter:='0'
      end
      else begin
         myabsolute:=number;
         if myabsolute >= place
         then digit
         else acharacter := '0'
      end;
      pluckdigit:=acharacter
end; (* pluckdigit *)
(* end module pluckdigit *)

begin

 (* according to:
http://agnes.dida.physik.uni-essen.de/~gnu-pascal/gpc_109.html#SEC109
*)
   GetTimeStamp(t); 

(* Predefined time stamp:
http://agnes.dida.physik.uni-essen.de/~gnu-pascal/gpc_389.html#SEC389
  TimeStamp = {@@packed} record
      DateValid,
      TimeValid   : Boolean;
      Year        : Integer;
      Month       : 1 .. 12;
      Day         : 1 .. 31;
      DayOfWeek   : 0 .. 6;  { 0 means Sunday }
      Hour        : 0 .. 23;
      Minute      : 0 .. 59;
      Second      : 0 .. 61; { to allow for leap seconds }
      MicroSecond : 0 .. 999999
    end;
*)

  with t do begin
    if TimeValid then begin
{
       writeln(output,'valid time');
       writeln(output,'year   =',year:4);
       writeln(output,'month  =',month:2);
       writeln(output,'day    =',day:2);
       writeln(output,'hour   =',hour:2);
       writeln(output,'minute =',minute:2);
       writeln(output,'second =',second:2);
}
       adatetime := 'year/mm/dd hh:mm:ss';
       adatetime[ 1] := pluckdigit(Year,3);
       adatetime[ 2] := pluckdigit(Year,2);
       adatetime[ 3] := pluckdigit(Year,1);
       adatetime[ 4] := pluckdigit(Year,0);
       adatetime[ 6] := pluckdigit(Month,1);
       adatetime[ 7] := pluckdigit(Month,0);
       adatetime[ 9] := pluckdigit(Day,1);
       adatetime[10] := pluckdigit(Day,0);
       adatetime[12] := pluckdigit(Hour,1);
       adatetime[13] := pluckdigit(Hour,0);
       adatetime[15] := pluckdigit(Minute,1);
       adatetime[16] := pluckdigit(Minute,0);
       adatetime[18] := pluckdigit(Second,1);
       adatetime[19] := pluckdigit(Second,0);
    end
    else begin
       writeln(output,'getdatetime: invalid time!');
       halt;
    end
  end;
{
   Sun compiler method:
   date(adatetime,'%Y/%m/%d %H:%M:%S');
}
end;
(* end module getdatetime version = 1.17; (@ of timegpc.p 2002 Oct 9 *)

(* begin module readdatetime *)
procedure readdatetime (var thefile: text; var adatetime: datetimearray);
(* read the date and time from the file.  It must have this format:
123456789 123456789
         1         
1980/06/09 18:49:11
 *)
(*
2001 Jan 2: thefile is not written to!  Force 4 digit years.
2000 Oct 11: upgraded so that the p2c compiler does not object
to writing out the adatetime; added checks for the date.
*)
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,'You have an old datetime (only 2 year digits): ');
      for index:=1 to datetimearraylength do write(output,adatetime[index]);
      writeln(output);
      writeln(output,'Convert your database to 4 digit years.');
      halt;
   end;
   (* check the adatetime format.  Note that further checks for the
      other positions in the array could be done to be sure that
      they are numbers.  But this should be pretty good. *)
   if (adatetime[ 5]<>'/') or
      (adatetime[ 8]<>'/') or
      (adatetime[14]<>':') or
      (adatetime[17]<>':')
   then begin
      writeln(output,'readdatetime: bad date time read:');
      for index:=1 to datetimearraylength do write(output,adatetime[index]);
      writeln(output);
      halt
   end;
end;
(* end module readdatetime version = 1.17; (@ of timegpc.p 2002 Oct 9 *)

(* 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 = 1.17; (@ of timegpc.p 2002 Oct 9 *)

(* begin module timeseed *)    
(* 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.  This
has the nice property that the seed changes every second and does not repeat
for thousands of years!  *)
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));
}
   n := ord(c) - ord('0');
   if (n < 0) or (n > 9) then begin
      writeln(output,'timeseed: error in datetime');
      writeln(output,'it contains "',c,'" which is not a number.');
      writeln(output,'The getdatetime routine must be fixed.');
      halt;
   end;
   seed := seed + power*n
end; (* addtoseed *)

procedure makeseed(adatetime: datetimearray; var seed: real);
(* convert adatetime to a real number in seed, reversed order
Here is the standard adatetime format:
123456789 123456789
         1         
1980/06/09 18:49:11
*)
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[ 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 = 1.17; (@ of timegpc.p 2002 Oct 9 *)
 
(*[[*)
(* 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 = 1.17; (@ of timegpc.p 2002 Oct 9 *)
(*]]*)

(* ************************************************************************ *)
(* end module package.datetime version = 1.17; (@ of timegpc.p 2002 Oct 9 *)

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

(* begin module clearstring *)
(* These modules clear strings in various ways *)
(* ---- *)
procedure emptystring(var ribbon: string);
(* empty the contents of the string but do NOT
remove the pointer.  This is useful for clearing
one string within a linked list of them. *)
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; (* emptystring *)
(* ---- *)
procedure clearstring(var ribbon: string);
(* empty the string and remove the pointer.
As of 2010 Dec 04 there is a previous pointer that
allows proper removal of the linked list using disposestring.
This clearstring routine is for the initial clearing
of the variables. *)
begin (* clearstring *)
   with ribbon do begin
      emptystring(ribbon);
      previous := nil;
      next := nil;
   end
end; (* clearstring *)
(* ---- *)
procedure disposestring(var ribbon: stringptr);
(* remove the linked list from the string in revere order *)
var
   t: stringptr; (* a link in the string *)
begin
   t := ribbon;
   (* find the end *)
   while t^.next <> nil do t := t^.next;
   (* move backwards and dispose links *)
   while t^.previous <> nil do begin
      t := t^.previous;
      dispose(t^.next);
   end;
   dispose(t);
end; (* disposestring *)
(* ---- *)
procedure initializestring(var ribbon: string);
(* start the string with a nil pointer.   This routine
should be called before doing linked list work.  This allows
the standard string routines to clear the string without
killing the pointer.  This is now deprecated, do not use it
since clearstring still clears the next pointer. *)
begin (* initializestring *)
   writeln(output,'remove initializestring routine!');
   writeln(output,'replace it with clearstring routine!');
   halt; (* to force deprecation *)
   clearstring(ribbon);
   ribbon.next := nil;
end; (* initializestring *)
(* end module clearstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module 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 writestring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* 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 = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module filler.filltrigger *)
procedure filltrigger(var t: trigger;
                      a: filler);
(* fill the trigger t *)
begin (* filltrigger *)
         fillstring(t.seek,a)
end; (* filltrigger *)
(* end module filler.filltrigger version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* 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 = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module skipblanks *)

(* 2003 July 31: tab is considered a blank character *)

function isblank(c: char): boolean;
(* is the character c blank or tab? *)
const
   tab = 9; (* tab character *)
begin
   isblank := (c = ' ') or (ord(c) = tab)
end;

procedure skipblanks(var thefile: text);
(* skip over blanks until a non-blank, or end of line, is found *)
begin
   while isblank(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 (not isblank(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 = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* ************************************************************************ *)
(* end module package.trigger version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* 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 = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* ********************************************************************** *)
(* ********************************************************************** *)
(* ********************************************************************** *)

(* begin module readquotestring *)
procedure readquotestring(var afile: text; var s: string);
(* read the quoted string from afile, remove the quotes.  Trailing
blanks are converted to tab characters. *)
const
   tabcharacter = 9; (* the ordinal of the tab character *)
var
   c: char; (* a character read *)
   i: integer; (* index to the string s *)
begin
   emptystring(s);
   skipblanks(afile);
   with s do begin
      read(afile,c);
      if c <> '"' then begin
         writeln(output, 'quote string expected but " missing');
         writeln(output, 'character found instead:"',c,'"');
         halt
      end;
      if eoln(afile) then begin
         writeln(output, 'missing end of quote string " string');
         halt
      end;
      repeat
         read(afile,c);
         if c <> '"' then begin
            length := succ(length);
            letters[length] := c;
         end;
      until (c = '"')
            or eoln(afile) or (length = maxstring);
      if length = maxstring then begin
         writeln(output, 'A quote string exceeded ',maxstring:1,
                         ' characters.');
         writeln(output, 'Use a shorter string or increase constant',
                         ' maxstring.');
         halt
      end;
      if c <> '"' then begin
         writeln(output, 'The second " mark of a string is missing.');
         halt
      end;
      (* convert all blanks to tab characters *)
      for i := 1 to length do begin
         if letters[i] = ' '
         then letters[i] := chr(tabcharacter); (* set to tab character *)
      end;
   end;
end;
(* end module readquotestring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module grabtoken *)
procedure grabtoken(var thefile: text; var thestring: string);
(* skip any blanks and then grab the next token from the file *)
var
   c: char; (* a character in thefile *)
   done: boolean; (* done finding the name *)
begin
   skipblanks(thefile);
   done := false;
   with thestring do begin
      length := 0;
      while not done do begin
         if eoln(thefile)
         then done := true
         else begin
            read(thefile,c);
            if c = ' '
            then done := true
            else begin
               length := succ(length);
               letters[length] := c;
            end
         end
      end
   end
end;
(* end module grabtoken version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module detabstring *)
procedure detabstring(var s: string);
(* Tabs are converted to blanks.  This allows the tab to be a special "solid"
kind of blank, but the string written won't vary with position on the line
depending on the tabstops. *)
const
   tabcharacter = 9; (* the ordinal of the tab character *)
var
   c: char; (* a character write out *)
   i: integer; (* index to the string s *)
begin
   with s do for i := 1 to length do begin
      c := letters[i];
{
for testing:
      if ord(c) = tabcharacter then c := 'T'; (* tab converted to blank *)
}
      if ord(c) = tabcharacter then c := ' '; (* tab converted to blank *)
      letters[i] := c
   end;
end;
(* end module detabstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module writequotestring *)
procedure writequotestring(var afile: text; s: string);
(* Write the string s to afile, surrounded by double quote marks.  Tabs are
converted to blanks.  This allows the tab to be a special "solid" kind of
blank, but the string written won't vary with position on the line depending
on the tabstops. *)
begin
   detabstring(s);
   write(afile, '"');
   writestring(afile, s);
   write(afile, '"');
end;
(* end module writequotestring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module readstringline *)
procedure readstringline(var afile: text; var buffer: string);
(* get a line (as a string) from a file. *)
var   index: integer; (* of buffer *)
begin (* getstring *)
   clearstring(buffer);
   index := 0;
   while (not eoln(afile)) and (index < maxstring) do begin
      index := succ(index);
      read(afile, buffer.letters[index])
   end;
   if not eoln(afile) then begin
      writeln(output, 'readstringline: a line exceeds maximum string size (',
                      maxstring:1,')');
      halt
   end;
   buffer.length := index;
   buffer.current := 1;
   readln(afile);
end; (* readstringline *)
(* end module readstringline version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module writestringlines *)
procedure writestringlines(var afile: text; l: stringptr);
(* write lines of strings to a file. *)
var
   h: stringptr; (* pointer to l *)
begin (* getstring *)
   if l <> nil then begin
      h := l;
      while h <> nil do begin
         writestring(afile, h^);
         writeln(afile);
         h := h^.next;
      end;
   end
end;
(* end module writestringlines version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

{
(* LOCKED begin module readheader *)
procedure readheader(var afile: text; ch: char; var theheader: stringptr);
(* Read the header from a file.  Only read lines that begin with ch.
2004 Dec 15: Handle end of file (eof) gracefully. *)
var
   h: stringptr; (* pointer to the riblheader *)
begin
   if eof(afile) then begin
      clearstring(theheader^);
   end
   else begin
      if afile^ = ch then begin
         new(theheader);
         h := theheader;
         while afile^='*' do begin
            readstringline(afile, h^);
            if not eof(afile) then begin
               if afile^='*' then begin
                  new(h^.next);
                  h := h^.next;
               end
               else h^.next := nil;
            end
         end
      end
      else begin
         theheader := nil
      end;
   end;
end;
(* LOCKED end module readheader version = 5.16; (@ of prgmod.p 2005 May 23 *)
}

(* begin module readheader *)
procedure readheader(var afile: text; ch: char; var theheader: stringptr);
(* Read the header from a file.  Only read lines that begin with ch.
2004 Dec 15: Handle end of file (eof) gracefully. *)
var
   h: stringptr; (* pointer to the riblheader *)
begin
   if eof(afile) then begin
      clearstring(theheader^);
   end
   else begin
      if afile^ = ch then begin
         new(theheader);
         h := theheader;
         while afile^='*' do begin
            readstringline(afile, h^);
            if not eof(afile) then begin
               if afile^='*' then begin
                  new(h^.next);
                  h := h^.next;
               end
               else h^.next := nil;
            end
         end
      end
      else begin
         theheader := nil
      end;
   end;
end;
(* end module readheader version = 5.67; (@ of prgmod.p 2017 Sep 29 *)


(* ********************************************************************** *)
(* ********************************************************************** *)
(* ********************************************************************** *)

(* 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 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* ************************************************************************ *)
(* end module package.brpiece version = 7.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* ************************************************************************ *)
(* end module package.getpiece version = 7.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* 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 = 20000; (* 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.75; {of delmod.p 2013 Jun 06} *)

(* ************************************************************************ *)
(* end module package.align version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module align.maxminalignment *)
procedure maxminalignment(var inst, book: text; var theline: integer;
                          var fromparam, toparam: integer;
                          alignmenttype: char);
(* prescan the book to find the range over which the pieces of the
   book are spread, relative to the aligned base.  the procedure uses
   the same variables that align does (so it can call align itself), and
   it returns the range in fromparam and toparam.
   alignmenttype: 'f' means alignment by First internal coordinate
      base, 'b' means alignment by Book, 'i' means alignment by Instructions.
 *)
const
      maximumrange = 20000; (* the maximum size aligned piece;
                             this will presumably catch the alignment bug *)
var
      distance: integer; (* a distance to the aligned base *)
      pie:  pieceptr;
      length, alignedbase: integer;
begin
      new(pie);
      (* set an initial range for the two bounds *)
      fromparam:=+maxint;
        toparam:=-maxint;

      reset(book);
      reset(inst);
      while not eof(book) do begin
         case alignmenttype of
            'i': align(inst,book,theline,pie,length,alignedbase);
            'b','f':  begin
                 getpiece(book,theline,pie); (* read in the piece *)
                 length := piecelength(pie);
            end;
         end;

         if not eof(book) then begin
            case alignmenttype of
               'f': begin (* force alignment on first base *)
                    alignedbase := 0;
                    fromparam := 1;
                    distance:=length-alignedbase;
                    if toparam < distance then toparam:=distance;
               end;
               'i': begin (* use the alignedbase from the book *)
                    distance:=1-alignedbase;
                    if fromparam > distance then fromparam:=distance;

                    distance:=length-alignedbase;
                    if toparam < distance then toparam:=distance;
               end;
               'b': begin (* use the internal book *)
                    alignedbase := pietoint(0, pie);
                    distance:=1-alignedbase;
                    if fromparam > distance then fromparam:=distance;

                    distance:=length-alignedbase;
                    if toparam < distance then toparam:=distance;
               end;
            end;

            clearpiece(pie)
         end
      end;

      if toparam - fromparam > maximumrange then begin
         writeln(output,' in procedure maxminalignment:');
         writeln(output,' alignedbase = ',alignedbase:1);
         writeln(output,' fromparameter = ',fromparam:1);
         writeln(output,' toparameter = ',toparam:1);
         writeln(output,' this exceeds the maximum range allowed (',
                        maximumrange:1,')');
         writeln(output,' see notes in the procedure. ');
         halt
         (* notes:  if you desired this range, increase 'maximumrange'.
            otherwise, this may indicate a bug - either:
            1) locate the bug (and tell tom schneider, please...)
            2) reduce the size of the fragments, from one or the other end
               until the bombing is stopped. *)
      end;

      (* make the book readable again *)
      reset(book);
      reset(inst);
      dispose(pie)
end;
(* end module align.maxminalignment version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module align.withinalignment *)
function withinalignment(alignedposition, alignedbase, length: integer)
                         :boolean;
(* this function tells one if an aligned position, relative to an aligned
   base in a piece of some length is within the piece. *)
var   p: integer;(* the position on the piece *)
begin
      p := alignedposition + alignedbase;
      withinalignment := (p>0) and (p<=length)
end;
(* end module align.withinalignment version = 7.75; {of delmod.p 2013 Jun 06} *)

(* 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.75; {of delmod.p 2013 Jun 06} *)

(* ********************************************************************** *)
(* ********************************************************************** *)
(* ********************************************************************** *)

(* begin module readwaveparameters *)
procedure readawaveparameter(var afile: text; var wp: waveparam);
{
                                 cmperbase, cmperbit: real);
}
(* Read one wave parameter from a file. *)
{
2000 Nov 1:  The following experimental code is not implemented:
See the coscurve program
(https://alum.mit.edu/www/toms/delila/coscurve.html) for details
about the conversion from bases to cm.  The code now requires
cm per base and cmperbit.
const
   lover2pi = 1.2160067233; (* L/2pi for L being the length of
       the cosine function from 0 to 2pi *)
var
   multiplier: real; (* factor to multiply by to convert from bases to cm. *)
}
begin
   readln(afile,wp.extreme);
   readln(afile,wp.wavelocation);
   readln(afile,wp.wavebit);
   readln(afile,wp.waveamplitude);
   readln(afile,wp.wavelength);
   if wp.wavelength <= 0.0 then begin
      writeln(output,'wave parameters: wavelength must be positive');
      halt
   end;
   (* one may skip the dash and thickness parameters
      of the last wave definition.  This allows
      backwards compatability, but should not be
      used in general. *)
   if eof(afile)
   then begin
      wp.dashon := 0;
      wp.dashoff := 0;
      wp.dashoffset := 0;
      wp.thickness := 0
   end
   else begin
      (* implement new method to read dash information
         without breaking previous wave parameter files *)
      if afile^ <> 'd'
      then begin (* default to old method *)
         readln(afile,wp.dashon);
         wp.dashoff := wp.dashon;
         wp.dashoffset := 0;
      end
      else begin
         get(afile); (* skip the 'd' *)
         (* reading the parameters directly in cm would require
            that the user compute or guess the values.
            So read in as bases and compute cm.
            Unfortunately the length of a curve is not the simple
            scale factor.  Make the user guess for now.
         *)
         readln(afile,wp.dashon, wp.dashoff, wp.dashoffset)

{
experimental code to account for actual wave length:
         readln(afile, wp.BASEdashon, wp.BASEdashoff, wp.BASEdashoffset);
         multiplier    := 1;
         multiplier    := lover2pi
 * wp.waveamplitude (* in bits *)
 * cmperbit (* converts to cm *)
 / wp.wavelength (* 10.6 bp/turn *)
* cmperbase
;
multiplier    := lover2pi;
multiplier    := 1.0;
multiplier    := 1.0/lover2pi;
         multiplier := lover2pi
                       * wp.waveamplitude (* in bits *)
                       * cmperbit (* converts to cm *)
                       * wp.wavelength (* 10.6 bp/turn *)
                       * cmperbase;
         multiplier := 1.0
         wp.dashon     := multiplier * wp.BASEdashon;
         wp.dashoff    := multiplier * wp.BASEdashoff;
         wp.dashoffset := multiplier * wp.BASEdashoffset;
with wp do begin
   writeln(output, 'cmperbase  = ',cmperbase:5:2);
   writeln(output, 'cmperbit   = ',cmperbit:5:2);
   writeln(output, 'dashon     = ',dashon:5:2);
   writeln(output, 'dashoff    = ',dashoff:5:2);
   writeln(output, 'dashoffset = ',dashoffset:5:2);
   writeln(output, 'waveamplitude   = ',wp.waveamplitude:5:2);
   writeln(output, 'wavelength      = ',wp.wavelength:5:2);
end;
}

      end;
      if eof(afile)
      then wp.thickness := 0
      else readln(afile,wp.thickness);
   end;
end;

procedure readwaveparameters(var afile: text; var w: waveptr);
{
                                 cmperbase, cmperbit: real);
}
(* read from afile the wave parameters w.  The routine is done either at end
of file, when a period ('.') appears or at a blank line.  This makes it
compatable with different reading mechanisms. *)
var
   done: boolean; (* done reading *)
   morethanone: boolean; (* more than one wave? *)
   p: waveptr; (* one of the definitions of a cosine wave *)
procedure waystoend;
(* if end of file, a blank line or a period, consider the job done *)
var
   clear: boolean; (* not eof, so clear to test for comment *)
begin
   if eof(afile) then done := true;
   clear := false;
   if not done (* skip any number of comments *)
   then while (not clear) and (not done) do begin
      if eof(afile) then begin
         done := true;
         clear := true;
      end
      else if afile^ = '*'
           then readln(afile)
           else clear := true
   end;
   if not done then if eoln(afile) then done := true;
   if not done then if afile^='.' then begin
      readln(afile);
      done := true;
   end;
end;
begin
   done := false;
   morethanone := false;
   w := nil;
   waystoend;
   if not done then begin
      new(w);
      p := w;
      while not done do begin
         waystoend;
         if done
         then p^.next := nil
         else begin
            if afile^='*'
            then readln(afile) (* skip any number of comments *)
            else begin
               if morethanone then begin
                  new(p^.next);
                  p := p^.next
               end;
{
               readawaveparameter(afile, p^, cmperbase, cmperbit);
}
               readawaveparameter(afile, p^);
               (* if we ever come back again there will be more than one: *)
               morethanone := true;
            end
         end
      end
   end
   else w := nil;
end;
(* end module readwaveparameters version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module writewaveparameters *)
procedure writewaveparameters(var afile: text; w: waveptr);
(* write wave parameters to a file *)
const
   wid = 4; (* width of numbers *)
var
   i: waveptr; (* pointer to the waves in w *)
begin
   i := w;
   while i <> nil do with i^ do begin
      writeln(afile,extreme,
                    '  extreme: char; h or l, the high or low extreme to be defined');
      writeln(afile,wavelocation:wid:1,
                    '  wavelocation: real; the location in bases of the extreme');
      writeln(afile,wavebit:wid:1,
                    '  wavebit: real; the location in bits of the extreme');
      writeln(afile,waveamplitude:wid:1,
                    '  waveamplitude: real; the amplitude of the wave in bits');
      writeln(afile,wavelength:wid:1,
                    '  wavelength: real; the wave length of the wave in bases');
{
      writeln(afile,dash:wid:1,
                    '  dash: real; the size of dashes in cm.',
                    '  dash <= 0 means no dashes');
      writeln(afile,    BASEdashon:wid:1,
                   ' ', BASEdashoff:wid:1,
                   ' ', BASEdashoffset:wid:1,
                   ' dashon, dashoff, dashoffset (in bases)');
}
      writeln(afile,    dashon:wid:1,
                   ' ', dashoff:wid:1,
                   ' ', dashoffset:wid:1,
                   ' dashon, dashoff, dashoffset (in bases)');
      writeln(afile,thickness:wid:1,
                    '  thickness: real; the thickness of the cosine wave.',
                    '  <=0 means default');
      i := i^.next
   end;
   writeln(afile,'.');
end;
(* end module writewaveparameters version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module putriblmatrix *)
procedure putriblmatrix(var afile: text;
                    matrix: ribltypeptr);
(* write the ribl matrix to afile. *)
var
   lindex: integer; (* index to Riblmatrix, equivalent to l *)
begin
   if matrix <> nil then begin
      with matrix^ do begin
         write(afile,'ribl');
         write(afile,' ');
         writequotestring(afile, riblname);
         writeln(afile);
         (* write the header *)
         writestringlines(afile, riblheader);
         writeln(afile, frombase: nfield,' ', tobase: nfield,' frombase, tobase');
         for lindex := frombase to tobase do begin
            write(afile,' ',
                       ' ',data[a, lindex]:infofield:infodecim,
                       ' ',data[c, lindex]:infofield:infodecim,
                       ' ',data[g, lindex]:infofield:infodecim,
                       ' ',data[t, lindex]:infofield:infodecim,
                       ' ',lindex:nfield);
            write(afile,' ',numbers[a,lindex]:nfield,
                        ' ',numbers[c,lindex]:nfield,
                        ' ',numbers[g,lindex]:nfield,
                        ' ',numbers[t,lindex]:nfield);
            writeln(afile);
         end;
         writeln(afile,'*');
         writeln(afile,mean:infofield:infodecim, 
                       ' bits = mean (Rsequence of selected region)');
         writeln(afile,stdev:infofield:infodecim,
                       ' bits = standard deviation');
         writeln(afile,'*');
         writeln(afile,consensus:infofield:infodecim,
                       ' bits = Ri of consensus sequence',
                       ' from ',frombase:1,' to ',tobase:1);
         writeln(afile,anticonsensus:infofield:infodecim,
                       ' bits = Ri of anticonsensus sequence',
                       ' from ',frombase:1,' to ',tobase:1);
         writeln(afile,'*');
         writeln(afile,averageRi:infofield:infodecim,
                       ' bits = average Ri for random sequence',
                       ' from ',frombase:1,' to ',tobase:1);
         writeln(afile,'*');
         writeln(afile,n:1,
                       ' n, number of sequences used to create the matrix');
         writeln(afile,'*');
         case symmetry of
           'a': write(afile, 'asymmetric');
           'e': write(afile, 'even');
           'o': write(afile, 'odd');
         end;
         writeln(afile, ' symmetry of the matrix');
         writeln(afile,'*');
         writeln(afile,Ribound:infofield:infodecim,
                       ' Ri bound: lower bound on Ri');
         writeln(afile,Zbound:infofield:infodecim,
                       ' Z bound: lower bound on Z');
         writeln(afile,Pbound:infofield:infodecim,
                       ' P bound: upper probability');
         writeln(afile,'*');
         writewaveparameters(afile, waves);
         writeln(afile);
      end
   end
end;
(* end module putriblmatrix version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* ********************************************************************** *)
(* ********************************************************************** *)
(* ********************************************************************** *)

(* begin module Ri.Riheader *)
procedure Riheader(var infile, book: text; c: char; var outfile: text);
(* do the header of the plot, using c as the comment character *)
var
      index: integer; (* to lines in the infile *)
begin
   reset(infile);

   writeln(outfile,c,' Ri ',version:4:2);

   writeln(outfile,c);
   writeln(outfile,c,' Ri(b,l) table is from:');

   (* copy header lines to outfile *)
   for index := 1 to 3 do begin
      write(outfile,c);
      copyaline(infile,outfile);
   end;

   writeln(outfile,c);

   reset(book);
   writeln(outfile,c,' BOOK/INST sequences are from:');
   copyaline(book,outfile);
   write(outfile,c,' ');
   reset(inst);
   if not eof(inst) then copyaline(inst,outfile)
                    else writeln(outfile, '(no instructions)');
   writeln(outfile,c);
end;
(* end module Ri.Riheader *)

(* ********************************************************************** *)
(* ********************************************************************** *)

(* begin module rs.readrsrange *)
procedure readrsrange(var rsdata: text; var r: rstype);
(* read the range data from rsdata to r.  data is assumed to
be the rsdata file from program rseq. *)
var
   index: integer; (* for counting lines of rsdata *)
   skip: char; (* a character to skip the '*' on the line *)
begin
      for index := 1 to 11 do readln(rsdata);
      readln(rsdata, skip, r.rstart, r.rstop);
    (*  writeln(output, 'range: ',r.rstart:1,' ',r.rstop:1);*)
end;
(* end module rs.readrsrange version = 4.75; (@ of rsgra.p 1990 Oct 2 *)

(* begin module rs.getrsbegin *)
procedure getrsbegin(var infile: text);
(* skip to the beginning of the data in a data file from
the rseq program *)
var
      ch: char; (* a character read from infile *)  
      begindata: trigger; (* a trigger to locate the beginning of the data *)
begin
(*                                       1         2  *) 
(*                              123456789012345678901 *)
filltrigger(begindata         ,'l    a    c    g    t');   
      resettrigger(begindata);  

      reset(infile);
      while not begindata.found do begin
         if eoln(infile) then readln(infile);
         if eof(infile) then begin
            writeln(output,'beginning of data not found');
            halt;
         end;
         read(infile,ch);
         testfortrigger(ch,begindata); 
      end;
      readln(infile);
end;
(* end module rs.getrsbegin version = 4.75; (@ of rsgra.p 1990 Oct 2 *)

(* begin module rs.readrsdata *)
procedure readrsdata(var rsdata: text; var rdata: rstype);
(* read data from the data file of program rseq into the datatype *)
begin
   with rdata do begin
      read(rsdata,l,nal,ncl,ngl,ntl,rsl,rs,varhnb,sumvar,nl,ehnb);
      (* skip spaces to find the flag: *)
      while rsdata^=' ' do get(rsdata);
      readln(rsdata,flag);
(* writeln(output,'readrsdata: l a c g t flag = ',
l:1,' ',nal:1,' ',ncl:1,' ',ngl:1,' ',ntl:1,' ',flag); *)
   end
end;
(* end module rs.readrsdata version = 4.75; (@ of rsgra.p 1990 Oct 2 *)

(* ********************************************************************** *)
(* ********************************************************************** *)

(* begin module ri.ricalc *)
function Ricalc(ehnb: real; nxl, nl: integer;
                niot: real; ricalctype:char): real;
(* calculate the individual information Ri(b,l) for a base x having nxl
numbers out of a total of nl numbers total.  Ehnb is 2 - e(n), where
e(n) is the sampling correction.
If ricalctype is 'l', use the Laplace Law of sucession
             f(b,l) = (k+1)/(n+2) = (n(b,l)+1)/(n(l)+2)
as the estimate for the probability for k items out of n total.
Otherwise use the Staden method (ricalc = 's').

A special buggy mode is available (ricalctype = 'b') so that the program can
be used to create weight matricies generated prior to 2001 May 10 when the
matrix bug was corrected.

The bug is related to how the program is to make a weight matrix value when
the frequency of a base is zero.  Theoretically this is negative infinity: Ri
= ~2+log2(f) (where ~2 is 2 with small sample correction).  Obviously this is
overkill since the next sequence might have the base.  However one can
estimate the frequency as 1/(n+2), where n is the number of sequences.  In
this case Ri = ~2+log2(1/(n+2)).  The bug was that the ~2 was not added.

The effect is that if you have built a model in which some positions do not
have a base, the weight will be about 2 bits lower than it should be.  This
will exclude some sites that should have been included.  (Such sites will
have a black background on the letter when displayed with lister.) It means
that some sites will have been missed.

*)
begin

{
writeln(output, 'nl   = ',nl:1);
writeln(output, 'nxl  = ',nxl:1);
writeln(output, 'niot = ',niot:1);
writeln(output, 'ricalctype = ',ricalctype);
}

   case ricalctype of

{ removed 2003 Aug 18 - this code is incomplete because
it leads to a matrix filled with niot in all positions!
      'n': begin (* number *)
           Ricalc := niot
      end;
}

      'l': begin (* Laplace *)
           Ricalc := 2.0 + ln((nxl+1)/(nl+2))/ln(2);
      end;

      's','b','n': begin
         if nl <= 0 then begin
            writeln(output,' Ricalc: a position in the data has',
                           ' less than 1 example!',
                           ' ehnb = ',ehnb:8:5,
                           ' nxl = ',nxl:1,
                           ' nl = ',nl:1);
            halt
         end;
         if nxl < 1 then begin
            if ricalctype = 'b' (* BUG mode *)
            then Ricalc :=        ln(1/(nl+niot))/ln(2) (* BUG *) 
            else begin
               if niot > 1.0
               then Ricalc := ehnb + ln(1/(nl+niot))/ln(2) (* CORRECT *)
                   (* 2001 May 10: ehnb inserted above to correct bug *)
               else Ricalc := +niot;
            end
         end
         else Ricalc := ehnb + ln(nxl/nl)/ln(2);
      end
   end
end;
(* end module ri.ricalc *)

{
% original bug version, 2001 May 10
% (* LOCKED begin module ri.ricalc *)
% function Ricalc(ehnb: real; nxl, nl: integer;
%                 niot: integer; Staden: boolean; Laplace: boolean): real;
% (* calculate the individual information Ri(b,l) for a base x having nxl
% numbers out of a total of nl numbers total.  Ehnb is 2 - e(n), where
% e(n) is the sampling correction.
% If Laplace is true, use the Laplace Law of sucession
%              f(b,l) = (k+1)/(n+2) = (n(b,l)+1)/(n(l)+2)
% as the estimate for the probability for k items out of n total. *)
% begin
%    if Laplace then begin
%       Ricalc := 2.0 + ln((nxl+1)/(nl+2))/ln(2);
%    end
%    else begin
%       if nl <= 0 then begin
%          writeln(output,' Ricalc: a position in the data has',
%                         ' less than 1 example!',
%                         ' ehnb = ',ehnb:8:5,
%                         ' nxl = ',nxl:1,
%                         ' nl = ',nl:1);
%          halt
%       end;
%       if nxl < 1 
%       then if Staden
%            then Ricalc := ln(1/(nl+niot))/ln(2)
%            else Ricalc := niot
%       else Ricalc := ehnb + ln(nxl/nl)/ln(2);
%    end
% end;
% (* LOCKED end module ri.ricalc *)
}

(* begin module min *)
function min(a, b: real): real;
(* return the minimum of a and b *)
begin
   if a < b then min := a
            else min := b
end;
(* end module min *)

(* begin module max *)
function max(a, b: real): real;
(* return the maximum of a and b *)
begin
   if a > b then max := a
            else max := b
end;
(* end module max *)

(* begin module ri.computerandomav *)
function computerandomav(Riblmatrix: rblarray;
                         thefrom, theto: integer): real;
(* compute the average of the response of the matrix to equiprobable random
sequences to the file fout, for the range from thefrom to theto. *)
var
   average: real; (* running average *)
   position: integer; (* a location in the aligned sequence *)
   sum: real; (* running sum of Ribl for one position *) 
   count: integer; (* the number of non-infinite ribl values at the position *)
procedure addin(l: integer; b: base);
(* add the base b to the running sum and increment count
only if the Ri(b,l) value is not negative infinity. *)
begin
  if Riblmatrix[b,l] > defnegativeinfinity
  then begin
     count := succ(count);
     sum := sum + Riblmatrix[b,l]
  end;
end;
begin
   average := 0.0;
   for position := thefrom to theto do begin
     sum := 0.0;
     count := 0;
     addin(position,a);
     addin(position,c);
     addin(position,g);
     addin(position,t);

     if count > 0 then average := average + sum/count
   end;
   computerandomav := average;

{
   writeln(output,' ',average:infofield:infodecim,
                   ' bits = average Ri for random sequence',
                   ' from ',thefrom:1,' to ',theto:1,' *');
}
end;
(* end module ri.computerandomav *)

(* begin module scan.matrixsymmetry *)
function matrixsymmetry(matrix: rblarray;
                        fromwanted, towanted: integer): char;
(* determine if the matrix has a dyad axis of symmetry and the
kind.  The results are:  a asymmetric, o: odd, e: even.
NOTE: the computation is done over the range wanted on the rblarray,
which has both negative and positive index values. *)
var
   b: base; (* a base in the matrix *)
   lindex: integer; (* index to Riblmatrix *)
   length: integer; (* length of the used part of the matrix *)
   symmetric: boolean; (* true if no violation of matrix
      symmetry has been found yet *)
begin
   symmetric := true;
   length := towanted-fromwanted+1;
{
(* show matrix! *)
lindex := fromwanted;
while (lindex <= towanted) and symmetric do begin
  for b := a to t do begin
  write(output, 'matrix[',basetochar(b),',',lindex:1,']=',matrix[b,lindex]:5:2,' ');
  end;
writeln(output);
lindex := succ(lindex)
end;
}
   lindex := fromwanted;
{
writeln(output,'length     = ',length:3);
writeln(output,'fromwanted = ',fromwanted:3);
writeln(output,'towanted   = ',towanted:3);
}
   while (lindex <= towanted) and symmetric do begin
{
writeln(output,'lindex = ',lindex:1);
}
      for b := a to t do begin
        if    matrix[b,lindex]
           <> matrix[complement(b),towanted - lindex + fromwanted]
        then symmetric := false;
{
write(output, 'matrix[b,',lindex:1,']=',matrix[b,lindex]:5:2,' ');
write(output, 'matrix[comp(b),',towanted - lindex + 1:1,']=',
               matrix[complement(b),towanted - lindex + 1]:5:2,' ');
writeln(output,symmetric)
}
      end;
      lindex := succ(lindex)
   end;
   if symmetric
   then if (length mod 2) = 0 then matrixsymmetry := 'e'
                              else matrixsymmetry := 'o'
   else matrixsymmetry := 'a'
end;
(* end module scan.matrixsymmetry *)

(* begin module ri.readparameters *)
procedure readparameters(var rip: text;
                         var thename: string;
                         var thefrom,theto,column: integer;
                         var lowerRi: real;
                         var upperRi: real;
                         var lowerValue: real;
                         var upperValue: real;
                         var printsequ: boolean;
                         var printrixyin: boolean;
                         var partials: char;
                         var niot: real;
                         var ricalctype: char;
                         var alignmenttype: char;
                         var Ribound: real;
                         var Zbound: real;
                         var Pbound: real);
(* read the parameters *)
var
   parameterversion: real; (* parameter version number *)
   hold: text; (* text file for upgrading rip *)
procedure die;
begin
   writeln(output, 'missing parameter');
   halt
end;
procedure upgradefrom233tonow;
begin
writeln(output, 'FIXING THE PARAMETER file, rip!');
writeln(output, 'EDIT THE rip FILE TO GIVE YOUR MODEL A NEW NAME');
rewrite(hold);
reset(rip);
readln(rip); (* skip header line *)
while not eof(rip) do copyaline(rip, hold); (* grab it *)
rewrite(rip); (* here goes!! *)
writeln(rip,version:4:2,
'          version of ri that this parameter file is designed for.');
writeln(rip,
'"NAME_THIS_MODEL_IN_THE_rip_FILE_AND_RERUN_ri"        name of the ribl matrix');
reset(hold);
while not eof(hold) do copyaline(hold, rip); (* put it back *)
(* now add the new lines: *)
writeln(rip,
'-5000         Ribound: real; the Ri boundary for this definition');
writeln(rip,
'100           Zbound:  real; the Z boundary for this definition');
writeln(rip,
'1.00          Pbound:  real; the probability boundary for this definition');
(* now prepare for proper reading *)
reset(rip);
   readln(rip, parameterversion);
if parameterversion <> version then begin
  writeln(output,'utter failure of upgrade attempt!');
  halt
end
end;
begin
   reset(rip);

   readln(rip, parameterversion);
   if parameterversion < updateversion then begin
      writeln(output, 'You have an old parameter file!');
      if parameterversion = 2.33
      then upgradefrom233tonow
      else halt
   end;

   if eof(rip) then begin
      writeln(output,'missing From-To parameters');
      halt
   end;

   readquotestring(rip, thename);
   readln(rip);
   writequotestring(output,thename);
   writeln(output);

   readln(rip, thefrom, theto);

   if thefrom > theto then begin
      writeln(output,'thefrom (=',thefrom:1,
                     ') must be less than TO (=',
                     theto:1,') for Ribl');
      halt;
   end;

   if eof(rip) then begin
      writeln(output,'missing column parameter');
      halt
   end;

   readln(rip, column);
   if column < 1 then begin
      writeln(output,'column parameter must be positive');
      halt
   end;

   if eof(rip) then begin
      writeln(output,'You are missing the Ri bound parameters');
      halt
   end;

   if rip^='a' then begin
      lowerRi := -maxint;
      upperRi := +maxint;
      readln(rip)
   end
   else readln(rip, lowerRi, upperRi);
{
writeln(output, 'lowerRi = ', lowerRi:10:5);
writeln(output, 'upperRi = ', upperRi:10:5);
halt;
}

   if eof(rip) then begin
      writeln(output,'You are missing the Value bound parameters');
      halt
   end;

   if rip^='a' then begin
      lowerValue := -maxint;
      upperValue := +maxint;
      readln(rip)
   end
   else readln(rip, lowerValue, upperValue);

   if eof(rip) then begin
      writeln(output,'You are missing the selection parameter');
      halt
   end;

   if rip^ = 'p' then printsequ := true
                 else printsequ := false;
   readln(rip);

   if eof(rip) then die;
   if rip^ = 'p' then printrixyin := true
                 else printrixyin := false;
   readln(rip);

   if eof(rip) then die;
   readln(rip, partials);
   if (not (partials = 'n')) and (not (partials = 'i'))
   then partials := '-';

   if eof(rip) then die;
   if (rip^ = 's') or (rip^ = 'b') or (rip^ = 'n')
   then begin (* Staden's method, read t for niot *)
      read(rip, ricalctype);
      readln(rip, niot);
      if niot < 0.0 then begin
         writeln(output,'t must be non-negative');
         halt
      end;
   end
   else if rip^ = 'l'
   then begin
      readln(rip, ricalctype);
      niot := 0.0;
   end
   else begin (* read negative infinity for niot *)
      (* note: a '-' should preceed a negative number! *)
      if not (rip^ in ['0','1','2','3','4','5','6','7','8','9','n','-'])
      then begin
         writeln(output,'niot must be in "rlbn-" or a digit');
         halt
      end;
      if rip^ = 'n' then get(rip);  (* allow 'n' first' *)
      ricalctype := 'n';
      readln(rip, niot);
   end;

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

   if eof(rip) then die; readln(rip,Ribound);

   if eof(rip) then die; readln(rip,Zbound);

   if eof(rip) then die; readln(rip,Pbound);

end;
(* end module ri.readparameters *)

(* begin module ri.writeparameters *)
procedure writeparameters(var fout: text;
                          thefrom,theto,column: integer;
                          lowerRi: real;
                          upperRi: real;
                          lowerValue: real;
                          upperValue: real;
                          printsequ: boolean;
                          printrixyin: boolean;
                          partials: char;
                          niot: real;
                          ricalctype: char;
                          alignmenttype: char);
(* write the parameters to file fout *)
begin
   writeln(fout,'* PARAMETERS FOR Ri:');
   writeln(fout,'* ', thefrom:1,' ',theto:1,' from-to');
   writeln(fout,'* ',column:1,' column of value file');
   writeln(fout,'* ',lowerRi:infofield:infodecim,
                   ' ', upperRi:infofield:infodecim,
                   ' lowest to highest Ri selected');
   writeln(fout,'* ',lowerValue:infofield:infodecim,
                   ' ', upperValue:infofield:infodecim,
                   ' lowest to highest Value selected');
   write(fout,'* ');
   if not printsequ then write(fout,'not')
                    else write(fout,'   ');
   writeln(fout,' printing sequences to sequ');

   write(fout,'* ');
   if not printrixyin then write(fout,'not')
                    else write(fout,'   ');
   writeln(fout,' printing sequences to rixyin');

   write(fout,'* ');
   case partials of
     'n': writeln(fout,'n: no line printed when partial site');
     'i': writeln(fout,'i: keep line printed when partial site,',
                       ' but force Ri = -infinity');
     '-': writeln(fout,'-: whole line printed when partial site');
   end;

   write(fout,'* ');
   if (ricalctype = 's') or (ricalctype = 'b') then begin
      writeln(fout,' using Staden''s Method: when f(b,l) = 0,',
                   ' replace with f(b,l) = 1/(n+t), t = ',niot:2:1);
      if ricalctype = 'b' then begin
         writeln(fout,
         '**************************************************************');
         writeln(fout,
         '* WARNING: BUG version: Ri - Ehnb for cases where f(b,l) = 0 *');
         writeln(fout,
         '**************************************************************');
      end;
   end
   else if ricalctype = 'l'
        then writeln(fout, 'using Laplace''s Method',
             ' f(b,l) = (k+1)/(n+2) = (n(b,l)+1)/(n(l)+2)')
        else writeln(fout,niot:infofield:infodecim,
                     ' is the value of negative infinity');

   writeln(fout,'* ');
   writeln(fout,'* ',alignmenttype, ' alignmenttype first, book, instructions');
end;
(* end module ri.writeparameters *)

(* begin module PatentMessage *)
procedure PatentMessage(var f: text);
(* Write the patent message to file f
version = 3.02 of PatentMessage_module 2017 Jul 31
*)
procedure makebar;
begin
writeln(f,
'**********************************************************************');
end;
begin
makebar;
writeln(f,
'* By downloading this code you agree to the Source Code Use License: *');
writeln(f,
'* https://alum.mit.edu/www/toms/Source_Code_Use_License_Schneider.txt *');
writeln(f,
'* Contact: https://alum.mit.edu/www/toms/contacts.html                *');
makebar;
end;
(* end module PatentMessage version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module ri.themain *)
procedure themain(var inst, book, rsdata, values, rip, wave,
                      rixyin, sequ, ribl: text);
(* the main procedure of the program *)
var
   matrix: ribltypeptr; (* the generated ribl weight matrix *)

   anticonsensus: real; (* the value of the anticonsensus sequence *)
   b: base; (* a base in one of the sites *)
   basestdev: real; (* part of the computation of stdev *)
   column: integer; (* column of values file to use *)
   character: char; (* a character of the sequence to write out *)
   columnindex: integer; (* index for counting columns of values file *)
   consensus: real; (* the value of the consensus sequence *)
   data: rstype; (* data from rseq *)
   dontkillpartials: boolean; (* don't kill partial sites when partials <> n *)
   fromsign: char; (* sign for riinst from *)
   lengthanalyzed: integer; (* length of region analyzed by Value(b,l) *)
   ln2: real; (* ln(2) *)
   lowerRi: real; (* lowest Ri to report *)
   lowerValue: real; (* lowest Value to report *)
   mean: real; (* mean of Ritotal = Rsequence *)
   n: integer; (* number of sequences *)
   niot: real; (* negative infinity or t *)
   ricalctype: char; (* how to calculate Ri:
     'l' means use Laplace's method
     's' means use Staden's method
     'b' means use Staden's method with BUG (Ehnb not added in)
      *)
   position: integer; (* a location in the aligned sequence *)
   partials: char; (* n: no line, i: -infinity, - : keep anyway *)
   fullsite: boolean; (* true if the site is complete (not partial) *)
   printsequ: boolean; (* true if sequences are being printed to sequ *)
   printrixyin: boolean; (* true if sequences are being printed to rixyin *)
   Ritotal: real; (* the total of Ri(b,l) for a site *)
   stdev: real; (* standard deviation of Ritotal *)
   sumRi: real; (* running sum of Ritotal *)
   sumRi2: real; (* running sum of Ritotal squared *)
   thefrom: integer; (* the from base *)
   theto: integer; (* the to base *)
   tosign: char; (* sign for riinst to *)
   upperRi: real; (* highest Ri to report *)
   upperValue: real; (* highest Value to report *)
   value: real; (* a value to compare to ri *)
   valuesfull: boolean; (* true if there are data in the values file *)

{zzzqqq}
   (* variables for computing standard deviation *)
   nl: integer; (* number of bases at position l, n(l) *)
   f: real; (* frequency to use for variance *)
   variance: real; (* variance for a base *)
   sumVariance: real; (* variance for the site *)

   (* variables used by the align routines: *)
   apiece: pieceptr;
   length, alignedbase: integer;
   fromparam, toparam: integer;

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

   theline: integer; (* location in the book *)

   hold: text; (* header lines *)
begin
   writeln(output,'Ri ',version:4:2);
   PatentMessage(output);
   limitdate ('9','9','9','9','    /99/99 99:99:99');

   ln2 := ln(2);
   new(apiece);

   (* copy header stuff *)
   brinit(book,theline);
   rewrite(rixyin);
   (* by making the first line as follows, alist can use the file directly *)
   writeln(rixyin,'* bits');
   Riheader(rsdata,book,'*',rixyin);
   rewrite(ribl);
   rewrite(riinst);
   writeln(riinst,'title "WARNING: you need org and chr instructions";');
   writeln(riinst,'(* ');
   Riheader(rsdata,book,'*',riinst);
   writeln(riinst,' *)');
   writeln(riinst);
   writeln(riinst, 'organism ORG; chromosome CHR;');

   new(matrix);

   readparameters(rip,
                  matrix^.riblname,
                  thefrom,theto,column,
                  lowerRi,upperRi,
                  lowerValue,upperValue,
                  printsequ,printrixyin,partials,
                  niot,ricalctype, alignmenttype,
                  matrix^.Ribound,
                  matrix^.Zbound,
                  matrix^.Pbound);

   if printsequ then rewrite(sequ);
   dontkillpartials := (partials <> 'n');

   writeparameters(output,thefrom,theto,column,
                   lowerRi,upperRi,
                   lowerValue,upperValue,
                   printsequ,printrixyin,partials,
                   niot,ricalctype, alignmenttype);
   writeparameters(rixyin,thefrom,theto,column,
                   lowerRi,upperRi,
                   lowerValue,upperValue,
                   printsequ,printrixyin,partials,
                   niot,ricalctype, alignmenttype);

   (* Generate header lines in hold file: *)
   rewrite(hold);
   Riheader(rsdata,book,'*',hold);
   writeparameters(hold,thefrom,theto,column,
                   lowerRi,upperRi,
                   lowerValue,upperValue,
                   printsequ,printrixyin,partials,
                   niot,ricalctype, alignmenttype);
   writeln(hold,'*',
                 ' ','Ri(a,l)':infofield,
                 ' ','Ri(c,l)':infofield,
                 ' ','Ri(g,l)':infofield,
                 ' ','Ri(t,l)':infofield,
                 ' ','l':nfield,
                 ' ','a':nfield,
                 ' ','c':nfield,
                 ' ','g':nfield,
                 ' ','t':nfield);
   reset(hold);
   (* now slurp them into the data structure *)
   readheader(hold,'*',matrix^.riblheader);

   writeln(rixyin,'* Lengths file:');
   reset(values);
   if eof(values)
   then begin
      writeln(rixyin,'* empty');
      valuesfull := false;
      value := 0.0; (* this is the value that will be written out *)
   end
   else begin
      while values^='*' do copyaline(values,rixyin);
      valuesfull := true;
   end;
   writeln(rixyin,'*');

   (* Obtain the data for and create the Ri(b,l) table *)
   if eof(rsdata) then begin
      writeln(output,'empty rsdata file');
      halt
   end;

   (* find the range of the graph in bases *)
   reset(rsdata);
   readrsrange(rsdata,data);

{ 2013 Feb 15: wrong test!!! 
   if data.rstop - data.rstart + 1 > maxribl then begin
      writeln(output,' width of site ',data.rstop - data.rstart + 1:1,
zzqq
}

   if data.rstop > maxribl then begin
      writeln(output,' From-range of site ',data.rstop:1,
                     ' exceeds maxribl ',maxribl:1);
      halt
   end;

   if data.rstart < minribl then begin
      writeln(output,' To-range of site ',data.rstart:1,
                     ' exceeds minribl ',minribl:1);
      halt
   end;

   getrsbegin(rsdata);

   (* prepare the inst and book for reading *)
   maxminalignment(inst, book, theline, fromparam, toparam, alignmenttype);

   with data do writeln(rixyin, '* data are from ', rstart:1, ' to ', rstop:1);
   writeln(rixyin, '* book/inst alignment is from ',
                 fromparam:1, ' to ', toparam:1);
   writeln(rixyin,'* Ri analysis is from ',thefrom:1,' to ',theto:1);

   with data do begin
      if alignmenttype <> 'f'
      then if fromparam > thefrom then begin
         writeln(output, ' Aligned FROM = ',fromparam:1,
                         ' > requested FROM = ',thefrom:1);
         halt
      end;

      if rstart > thefrom then begin
         writeln(output, ' In file rsdata the FROM = ',rstart:1,
                         ' > requested FROM = ',thefrom:1);
         halt
      end;

      if alignmenttype <> 'f'
      then if toparam < theto then begin
         writeln(output, ' Aligned TO = ',toparam:1,
                         ' > requested TO = ',theto:1);
         halt
      end;

      if rstop < theto then begin
         writeln(output, ' In file rsdata the TO = ',rstop:1,
                         ' > requested TO = ',theto:1);
         halt
      end;
   end;

   (* prepare for reading the rsdata *)
   getrsbegin(rsdata);

   (* create the Ri(b,l) *)

   matrix^.frombase := thefrom;
   matrix^.tobase := theto;

   consensus := 0.0;
   anticonsensus := 0.0;
   for position := data.rstart to data.rstop do begin
      (* skip lines with an '*' *)
      if rsdata^ <> '*' then begin
         readrsdata(rsdata,data);
         with data do begin
            if position <> l
            then writeln(output, 'Warning: position should be',
                               ' ',position:1,', but is actually',
                               ' ',l:1);

            if (position >= thefrom) and (position <= theto) then begin
               matrix^.data[a, l] := Ricalc(ehnb,nal,nl,niot,ricalctype);
               matrix^.data[c, l] := Ricalc(ehnb,ncl,nl,niot,ricalctype);
               matrix^.data[g, l] := Ricalc(ehnb,ngl,nl,niot,ricalctype);
               matrix^.data[t, l] := Ricalc(ehnb,ntl,nl,niot,ricalctype);
               matrix^.numbers[a, l] := nal;
               matrix^.numbers[c, l] := ncl;
               matrix^.numbers[g, l] := ngl;
               matrix^.numbers[t, l] := ntl;

               consensus := max(
                                max(
                                    max(matrix^.data[a, l],
                                        matrix^.data[c, l]),
                                    matrix^.data[g, l]),
                                matrix^.data[t, l])
                            + consensus;

               anticonsensus := min(
                                    min(
                                        min(matrix^.data[a, l],
                                            matrix^.data[c, l]),
                                        matrix^.data[g, l]),
                                    matrix^.data[t, l])
                                + anticonsensus;
            end;
         end;
      end;
   end;

   (* Read the book using inst to align the pieces *)

   writeln(rixyin,'*');
   writeln(rixyin,'* Columns:');
   writeln(rixyin,'*  1 piece number');
   writeln(rixyin,'*  2 piece name');
   writeln(rixyin,'*  3 sequence region analyzed (if printed, - if not)');
   writeln(rixyin,'*  4 length of region analyzed on this piece ');
   writeln(rixyin,'*  5 aligning coordinate on piece');
   writeln(rixyin,'*  6 Rindividual for the piece');
   writeln(rixyin,'*  7 value from the values file');
   writeln(rixyin,'*  8 standard deviation of Rindividual for that sequence');
   writeln(rixyin,'*');

   sumRi := 0.0;
   sumRi2 := 0.0;
   n := 0;
   while not eof(book) do begin
      case alignmenttype of
         'i': align(inst,book, theline,apiece,length,alignedbase);
         'b','f': begin
              getpiece(book, theline,apiece); (* read in the piece *)
              length := piecelength(apiece);
         end;
      end;

      if not eof(book) then begin
         n := n + 1; (* count the sequences found *)

         case alignmenttype of
            'f': alignedbase := 1; (* force alignment on first base *)
            'i': ; (* use the alignedbase from the book *)

            'b': alignedbase := pietoint(0, apiece); (* use the internal book *)
         end;

         Ritotal := 0.0;
         sumVariance := 0.0;
         fullsite := true; (* innocent till proven guilty! *)
         for position:=thefrom to theto do begin
            if withinalignment(position, alignedbase, length) then begin
               b := getbase(position+alignedbase, apiece);
               Ritotal := Ritotal + matrix^.data[b, position];

{
If the frequency is zero, then use
 
f := 1/(2+n(l))
 
otherwise use
f := f(b,l)
The (corrected to account for n(l) variation) formula given in the perl
program is
   
variance = ((1-f)/f)/n(l)
} {zzzqqq}

               (* compute standard deviation *)
               (* compute frequency to use *)
               nl :=   matrix^.numbers[a, position]
                     + matrix^.numbers[c, position]
                     + matrix^.numbers[g, position]
                     + matrix^.numbers[t, position];
               if matrix^.numbers[b, position] = 0
               then f := 1/(2+nl)
               else f := matrix^.numbers[b, position] / nl;
               variance := ((1-f)/f)/nl;
               sumVariance := sumVariance + variance;

            end
            else fullsite := false (* proven guilty! *)
         end;
 
         (* obtain values from the values file *)
         if valuesfull then begin
            if not eof(values) then begin
               while values^='*' do readln(values);
               columnindex := 1;
               while columnindex < column do begin
                  skipblanks(values);
                  skipnonblanks(values);
                  columnindex := columnindex + 1
               end;
               if eoln(values) then begin
                  writeln(output,'Missing data column ',
                                  column:1,' in values file');
                  halt
               end;
               readln(values, value)
            end
            else begin
               value := 0.0;
               write(output,'Note: A zero was inserted into the rixyin file',
                            ' value column');
               if numbered then writeln(output,' for piece number ', number:1)
                           else writeln(output);
               writeln(output,'because the values file has finished.');
            end
         end;

         (* primary selection *)
         if (Ritotal >= lowerRi) and (Ritotal <= upperRi)
            and (fullsite or dontkillpartials)
         then begin

            (* secondary selection *)
            if (value >= lowerValue) and (value <= upperValue) then begin
               sumRi := sumRi + Ritotal;
               sumRi2 := sumRi2 + Ritotal*Ritotal;

               if numbered
               then write(rixyin, ' ', number:6, ' ')
               else write(rixyin, ' (no.#) ');

               (* name of the piece *)
               write(rixyin, ' ',apiece^.key.hea.keynam.letters);

               (* print the sequence and determine length of analyzed region *)
               if printrixyin then write(rixyin, ' ')
                              else write(rixyin, ' -');
               
               lengthanalyzed := 0;
               for position := thefrom to theto do begin
                  if withinalignment(position, alignedbase, length) then begin
                     character := basetochar(getbase(position+alignedbase,
                                                     apiece));
                     if printsequ then write(sequ,character);
                     if printrixyin then write(rixyin,character);
                     lengthanalyzed := lengthanalyzed + 1;
                  end
                  else if printrixyin then write(rixyin, '-');
               end;
               if printsequ then writeln(sequ,'.');

               (* length of the sequence analyzed *)
               write(rixyin, ' ', lengthanalyzed:nfield);

               (* coordinate of aligning base *)
               write(rixyin, ' ', inttopie(alignedbase,apiece):nfield);

               (* Ri *)
               if fullsite or (partials <> 'i')
               then write(rixyin, ' ',Ritotal:infofield:infodecim)
               (* I had niot here, but defnegativeinifinity is better
                  because it assures that partial sites are removed: *)
               else write(rixyin, ' ',defnegativeinfinity:infofield);

               write(rixyin, ' ', value:infofield:infodecim);

               write(rixyin, ' ', sqr(sumVariance):infofield:infodecim);

               writeln(rixyin);

               (* make partial riinst file: *)
               (* name of the piece *)
               write(riinst, 'piece ',apiece^.key.hea.keynam.letters:
                              apiece^.key.hea.keynam.length ,';');
               write(riinst,' get from ');
               write(riinst, inttopie(alignedbase,apiece):1);

               (* give the orientation *)
               case apiece^.key.piedir of
                 minus: begin
                    fromsign := '+';
                      tosign := '-';
                 end;
                 plus: begin
                    fromsign := '-';
                      tosign := '+';
                 end;
               end;

               write(riinst,' ',fromsign,'FROM to same ',
                                  tosign,'TO direction ',
                                  tosign);

{
               if apiece^.key.piedir = minus then write(riinst ,'-')
                                             else write(riinst ,'+');

}
               writeln(riinst,'; (* ', Ritotal:1:infodecim,' bits *)');
            end;
         end;

         clearpiece(apiece);
      end
   end;

   if n = 0 then begin
      writeln(output,'NO SEQUENCES FOUND IN BOOK!');
      rewrite(ribl);
      halt
   end;

   if n = 1 then begin
      writeln(output,'ONLY ONE SEQUENCE FOUND IN BOOK!');
      rewrite(ribl);
      halt
   end;

   if n = 2 then begin
      writeln(output,'ONLY TWO SEQUENCES FOUND IN BOOK!');
      rewrite(ribl);
      halt
   end;

   (* Calculation of the standard deviation must use n-1 since
      the mean is estimated.  See
         @book{Press1989,
         author = "W. H. Press
          and B. P. Flannery
          and S. A. Teukolsky
          and W. T. Vetterling",
         title = "Numerical Recipes in Pascal.
         The Art of Scientific Computing",
         publisher = "Cambridge University Press",
         address = "Cambridge",
         year = "1989"}
      page 500. *)

   mean := sumRi / n;

{
writeln(output,'n=',n:1);
writeln(output,'n-1=',n-1:1);
writeln(output,'sumRi=',sumRi:infofield:infodecim);
writeln(output,'mean=',mean:infofield:infodecim);
writeln(output,'mean*mean=',mean*mean:infofield:infodecim);
writeln(output,'sumRi2/n - (mean*mean)=',sumRi2/n - (mean*mean):infofield:infodecim);
writeln(output,'sqrt(sumRi2/n - (mean*mean))=',sqrt(sumRi2/n - (mean*mean)):infofield:infodecim);
}

   (* precompute this because the difference can be so close to
      zero that negative values are obtained.  Then the square root
      function objects *)
   basestdev := sumRi2/n - (mean*mean);

   if basestdev > 0
   then stdev := sqrt(n/(n-1))*sqrt(sumRi2/n - (mean*mean))
   else stdev := 0.0;

   matrix^.mean := mean;
   matrix^.stdev := stdev;
   matrix^.consensus := consensus;
   matrix^.anticonsensus := anticonsensus;
   matrix^.averageRi := computerandomav(matrix^.data,thefrom,theto);
   matrix^.n := n;
   matrix^.symmetry := matrixsymmetry(matrix^.data, thefrom, theto);

   reset(wave);
   readwaveparameters(wave, matrix^.waves);

   case matrix^.symmetry of
     'a': write(output,'asymmetric');
     'e': write(output,'even');
     'o': write(output,'odd');
   end;
   writeln(output, ' symmetry of the matrix');

   putriblmatrix(ribl, matrix);
end;
(* end module ri.themain *)

begin
   themain(inst, book, rsdata, values, rip, wave, rixyin, sequ, ribl);
1: end.
