program orf(genomebook, scanfeatures, orfp, orffeatures, output);
(* orf: find ORFs for ribosome binding sites

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

 *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 1.36; (* of orf.p 2022 Dec 16
2022 Dec 16, 1.36: upgrade modules
2017 Oct 10, 1.35: continue debug, cite Hemm.Rudd2008
2017 Oct 02, 1.34: understand 'missing feature' [in caps] message
2017 Sep 29, 1.33: integrate with new multiscan groups
2016 Oct 09, 1.32: get new modules from lister.p; had to add groupnumber
2013 Nov 04, 1.31: clean up
2013 Oct 31, 1.30: clean up
2013 Oct 31, 1.29: clean up!
2013 Mar 08, 1.28: update address, can use book for genome.
2010 Jan 06, 1.27: don't put out any features if none in scanfeatures!
2010 Jan 06, 1.26: bombs in readdefinitionsandfeatures - just upgrade modules!
2008 Mar 25, 1.25: adjust documentation
2007 Apr 18, 1.24: clean up
2007 Apr 18, 1.23: bug when using complement extracted book for genomebook
2007 Apr 18, 1.21: orf orientation as <###] or [###>
2007 Apr 18, 1.20: correct direction of orf placement!
2007 Apr 18, 1.19: clean up
2007 Apr 17, 1.18: orientation bug: reversed sequence gets wrong orf
2007 Apr 07, 1.17: orffeatures must report total RBS info
2007 Apr 07, 1.16: upgrade modules
2007 Apr 07, 1.15: fix documentation
2007 Mar 26, 1.14: parameter to limit ors reported ??
2007 Mar 26, 1.13: total line records to the last base of stop codon.
                   feature, of course, to the first base.
2007 Mar 26, 1.12: small changes
2007 Mar 26, 1.11: color bar behind orf matches ir color - do in lister!
2007 Mar 26, 1.10: color bar behind orf!
2007 Mar 26, 1.09: put # codon into orf feature string
2007 Mar 25, 1.08: orf feature: write to orffeatures
2007 Mar 25, 1.07: orf feature: built!
2007 Mar 25, 1.06: clean up
2007 Mar 25, 1.05: use Bparam: orf in bases, Cparam: orf-end
2007 Mar 25, 1.04: use Aparam to store orf length (in codons)
2007 Mar 25, 1.03: reading the orf to the stop!
2007 Mar 25, 1.02: start to read features and compare to book piece
2007 Mar 24, 1.01: have parts to read genomebook and scanfeatures
2007 Mar 24, 1.00: origin
*)
updateversion = 1.00; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.orf *)
(*
name
   orf: find ORFs for ribosome binding sites

synopsis

   orf(genomebook:in, scanfeatures: in, orfp: in, orffeatures: out,
       output: out)

files

   genomebook:  a book from the Delila system containing a genome.
     If a genome is not available, a regular Delila book works ok.

   scanfeatures:  output of the multiscan program containing
     locations of Initiation Region (ir) features.

   orfp:  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.

      shortestorf: the shortest orf to report

      longestorf:  the  longest orf to report

   orffeatures:  ORF features reported as features for the lister program
      This file REPLACEs the scanfeatures in the features file.

   output: messages to the user

description

   The orf program reads a book containing a complete genome and,
   given the start points of translation (initiation regions, ir, the
   AUG, GUG or UUG starts) it finds the corresponding stop codons. 
   These are called 'orf' features.  They reside at the end of the
   open reading frame and report the length of the frame in codons. 
   The length includes the initiation and termination codons.

   Be sure to turn on 'predict peptides' in listerp!

   The procedure used by the program is:

   * Read in a whole genome book (or Delila book).

   * Read in a scanfeatures.

   * Look through the scanfeatures for ir features.  Determine the orf
   for each ir feature in the given orientation at that position in
   the genome.

   Note that the orffeatures REPLACE the scanfeatures in a features
   file and are given to the lister.p program for display.  This is
   because the scanfeatures are read in and then modified: the total
   lines corresponding to the ir get orf data added to them and an orf
   definition and features are inserted into the list.

   The orffeatures have special parameters for total features:

      The Aparam of the total feature is the information (in bits).

      The Bparam of the total feature is the orf length (in codons,
      including the stop).

      The Cparam of the total feature is the last base of the orf.

      The Dparam: see below.

   The parameters shortestorf and longestorf only affect one thing,
   the Dparam of the total feature of the ribosome binding site.  If
   the number of codons (including initiation codon and stop codon) is
   greater than or equal to shortestorf and less than or equal to
   longestorf then the total Dparam is set to '1', otherwise it is 0
   (zero).  This allows selection of the total features to generate
   tables.

examples

   Example parameters to use in listerp:
39    basesperline: number of bases per line in the listing
1     aastate: 0=no aa; 1=predict peptides; 2=translate all frames
7     frameallowed: binary; highest bit is highest frame on, etc.
1     codelength: 1 or 3 letters per amino acid

basesperline: must be a multiple of 3 for peptides
aastate:      1=predict peptides
frameallowed: 7 for all frames
codelength: 1 or 3 letters per amino acid could do either.

documentation

  We used the orf program to identify small proteins:

@article{Hemm.Rudd2008,
author = "M. R. Hemm
 and B. J. Paul
 and T. D. Schneider
 and G. Storz
 and K. E. Rudd",
title = "{Small membrane proteins found by comparative genomics and
ribosome binding site models}",
journal = "Mol. Microbiol.",
volume = "70",
pages = "1487--1501",
pmid = "19121005",
pmcid = "PMC2614699",
note = "\htmladdnormallink
{https://doi.org/10.1111/j.1365-2958.2008.06495.x}
{https://doi.org/10.1111/j.1365-2958.2008.06495.x}, \htmladdnormallink
{https://alum.mit.edu/www/toms/papers/smallproteins/}
{https://alum.mit.edu/www/toms/papers/smallproteins/}",
year = "2008"}

see also

   {program for display the results:} lister.p 

author

   Thomas Dana Schneider

bugs

   Orf is designed for a single bacterial genome at the moment, since
   it only handles one piece at a time.

   Can we use the ir color for the stop codon color bar?  No!  The
   color is available yet because mkpetal happens LATER - there is no
   way to get that color now!  This function can be done only in
   lister when there are colors available after mkpetals has been run.
   Method for lister: search the features for ir.  If the ir is
   followed by an orf, assign the color of the ir to the orf.

   When two ribosome binding sites are in frame, two orf stops are
   generated that are in the same place.  They only differ by the
   other string and so in the current lister they are considered to be
   duplicates.  Two solutions are (1) make lister check the 'other'
   string for identity and (2) do an 'other' string check only when
   the namestring is 'orf'.

technical notes

*)
(* end module describe.orf *)

globaldebug = false;

superdebug = false; (* allow special debugging mode for figuring out
                       vertical colored bars *)

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

(* begin module string.const *)
      maxstring = 2000; (* the maximum string *)
(* end module string.const version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(*[[*)
(* 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.85; (@ of prgmod.p 2022 Dec 08 *)
(*]]*)

(* begin module lister.const *)
   (* defaults for settable variables *)
   defbasesperline = 30;
   defaastate = 1;
   defframeallowed = 7;
   defcodelength = 3;
   defseqlines = 1;
   defpageaction = 'c'; (* page action *)
   defpagelength = 55; (* number of lines on a page *)
   defhopi = 'HOPI'; (* parts of display to print *)
   defnkfc = 'NKFC'; (* parts of piece information to print *)
   defndaf = 'NDAF'; (* default ndaf data controls *)
   defsof = 'SOF'; (* default sof feature controls *)
   prime = ''''; (* single quote mark *)

   shiftright = 3; (* number of characters to shift the features right
                      on the page *)
   tabcharacter = 9; (* the ordinal of the tab character *)
   maxdatacontrol = 4; (* maximum number of data controls *)
   maxsof = 3; (* maximum number of sof feature controls *)

   (* default bounding box definitions.  See the technical notes for how to
      set this for a different printer. *)


   (* The following bounding box is for the Canon Color Laser Copier 1150. *)
   defaultllx =   7.10999; (* default for llx, lower left x *)
   defaultlly =   7.01995; (* default for lly, lower left y *)
   defaulturx = 588.15;    (* default for urx, upper right x *)
   defaultury = 784.98;    (* default for ury, upper right y *)

{
   (* The following bounding box is for the Canon Color Laser Copier 800S. *)
   defaultllx =  10.08; (* lower left x *)
   defaultlly =   8.91; (* lower left y *)
   defaulturx = 588.06; (* upper left x *)
   defaultury = 779.85; (* upper left y *)
}

{
   (* The following bounding box is for the inkjet Tektronix Phaser 140. *)
   defaultllx =  14.4; (* lower left x *)
   defaultlly =  28.9; (* lower left y *)
   defaulturx = 596.8; (* upper left x *)
   defaultury = 777.7; (* upper left y *)
}

(* PostScript constants *)
   pwid = 8; (* width in character places to print PostScript numbers *)
   pdec = 5; (* decimal places to print PostScript numbers *)
   pdecolor = 4; (* decimal places for color descriptions (5 WILL CAUSE
                    NeWS 1.1 TO BOMB) *)

   defcharacterratio = 0.6; (* default for the ratio of the rightward
      horizontal motion when a character is drawn to the font height.  This is
      determined by the map file, but this program needs to know it also to
      do the mark computations. *)
   tolerateratio = 10000; (* the tolerance in decimal place to accept
      in comparing the value of defcharacterratio with that determined
      by the PostScript interpreter *)

   linemax = 1000; (* longest output line possible *)

    tenspacemark = '*'; (* the symbol to put every 10 bases 2012 Oct 25 *)
   fivespacemark = '.'; (* the symbol to put every 5 bases that is not
                          tenspacemark. 2012 Oct 25 *)

(* end module lister.const version = 13.24; (@ of lister.p 2016 Oct 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

type
(* 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.70; {of delmod.p 2005 Sep 15} *)

      (* 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 = 13.24; (@ of lister.p 2016 Oct 08 *)

(* begin module amino.type *)
      aminoacid = packed array[1..3] of char;
(* end module amino.type version = 13.24; (@ of lister.p 2016 Oct 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(*[[*)
(* 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.85; (@ of prgmod.p 2022 Dec 08 *)
(*]]*)

(* begin module definition.feature.type *)

   (* Features and their definitions for the lister program *)
   petaltype = record
      (* define parts of a colored rectangle around a walker *)
      el: integer; (* edgelinewidth: edge linewidth (integer) *)
      c : char;    (* color kind: how color is defined: r for RGB, h for HSB *)
      (* ---------------------------------------------- *)
      eh: real;    (* edgeh: edge hue        OR red     *)
      es: real;    (* edges: edge saturation OR blue    *)
      eb: real;    (* edgeb: edge brightness OR green   *)
      (* ---------------------------------------------- *)
      fh: real;    (* fillh: fill hue        OR red     *)
      fs: real;    (* fills: fill saturation OR blue    *)
      fb: real;    (* fillb: fill brightness OR green   *)
      (* ---------------------------------------------- *)
      spare: real;    (* spare variable *)
   end;

   commenttypeptr = ^commenttype; (* pointer to commettype *)
   commenttype = record (* store a comment *)
      name: string; (* name of the comment *)
      next: commenttypeptr; (* the next definition recorded *)
   end;

   definitiontypeptr = ^definitiontype; (* pointer to definitiontype *)
   definitiontype = record (* store a definition *)
(*
         define "Site" "-" "<v0*]" "[*0^>" -7 -3.5 0 +3.0 +5
*)
      nametag: string; (* name of the feature *)
      background: string; (* background padding characters *)
      negparts: string; (* parts of the feature display *)
      posparts: string; (* parts of the feature display *)
      locations: array [1..maxstring] of real; (* locations of the marks
         of the feature *)
      (* ---------------------------- *)
      (* computed values: *)
      marks: integer; (* the number of marks used in this feature *)
      min: real; (* one end of the feature (minimum location) *)
      max: real; (* the other end of the feature (maximum location) *)
      definitionnumber: integer; (* the number of the definition read in *)
      matrix: ribltypeptr; (* an ribl matrix if available *)
      (* The following boundaries are DIFFERENT from the originally defined
         boundaries of the matrix.  This allows the matrix to keep its bounds
         (eg for default values for scanning AND it allows the lister program
         to define separate boundaries for displaying green versus red on the
         zero bar.  *)
      RiboundFeature: real; (* the Ri boundary for this feature *)
      ZboundFeature: real; (* the Z boundary for this feature *)
      PboundFeature: real; (* the probability boundary for this feature *)
      petal: petaltype; (* color rectangle background *)
      isgap: boolean; (* the definition represents a gap between two
                         sequence walkers *)
      istotal: boolean; (* the definition represents a total feature for
                           a set of sequence walkers *)
      ismulti: boolean; (* the definition represents a sequence walker 
                           that is part of a set of sequence walkers *)
      next: definitiontypeptr; (* the next definition recorded *)
   end;

(* ---------------------------- *)
   featurecarrierptr = ^featurecarrier; (* pointer to featurecarrier *)
   featuretypeptr = ^featuretype; (* pointer to featuretype *)
   featurecarrier = record (* carries a list of features *)
      (* to make a list of features, we need a carrier that
         is not the original feature list. *)
      afeature: featuretypeptr; (* the feature being carried *)
      coo: integer; (* coordinate of total *)
      Ri: real; (* Ri of total *)
      keepgroup: boolean; (* keep this group when printing *)
      next: featurecarrierptr; (* the next carriage *)
   end;
   featuretype = record (* store a feature *)
(*
         @ K01789 270 +1 "Site" "Positive 270"
         @ K01789 270 -1 "Site" "Negative 270"
         @ K01789 280 +1 "Site" "Positive 280"
         @ K01789 229.0 -1 "dnaA" "+12.2 bits " 12.200338 -0.473212 0.318031
*)
      id: string; (* name of the sequence piece the feature is on *)
      featurenumber: integer; (* the number of the feature read in *)
      featuregroup: integer; (* the group the feature is a member of *)
      featureisfirstingroup: boolean; (* is this the first feature of a group? *)
      coordinate: real; (* location of zero base *)
      orientation: integer; (* -1 or +1 is in direction of piece *)
      nametag: string; (* the name tag - name of the feature *)
      othertag: string; (* the other tag - another tag for misc use *)
      (* ---------------------------- *)
      definition: definitiontypeptr; (* the definition of the feature *)
      (* ---------------------------- *)
      (* the following three values are only defined when there
         is a definition for the feature matrix *)
      Ri: real; (* individual information of this site *)
      Z: real; (* Z score for this site *)
      probability: real; (* probability for this site *)

      (* Three parameters replace Ri, Z and p when there is no matrix,
         Aparam, Bparam and Cparam.  If there is a matrix, Tparam
         will be blank.  Otherwise: *)
      Tparam: char; (* a parameter that defines the kind of parameter
                       that the remainder parameters are. *)
      Aparam: real; (* a parameter *)
      Bparam: real; (* a parameter *)
      Cparam: real; (* a parameter *)
      Dparam: real; (* a parameter *)
      relcoo: real; (* relative coordinate *)
         (* Since these parameters
         are now (2004 July 17) able to substitute for Ri, Z and
         probability, the Tparam cannot be placed in front of the
         three parameters in feature files.  Tparamvalues:
            h = hsb feature colors: hue, saturation, brightness
            r = rgb feature colors: red, green, blue
            b = bit range for the feature
            ' ' = defaults to sequence walker
         Summary:
         Tparam  Aparam     Bparam     Cparam     Dparam     For
          ' '    Ri         Z          P          (none)     walker
          'h'    hue        saturation brightness (none)     ascii lines
          'r'    red        green      blue       (none)     ascii lines
          'H'    hue        saturation brightness thickness  color bar
          'R'    red        green      blue       thickness  color bar
         In this scheme one can control the color of the bar as an
         ascii string OR one can control the colored rectangle behind
         the bar but NOT BOTH.
         *)
      (* ---------------------------- *)
      evencoordinate: real; (* location of zero base for asymmetric and odd
         symmetry sites, with 0.0 as the decimal.  For even symmetry sites
         it's the point of symmetry, with 0.5 as the decimal. *)
      (* ---------------------------- *)
      (* computed values: *)
      unsatisfied: boolean; (* true when the feature has not yet been used
                               during printing a line *)
      fromrange: real; (* the lower part of the feature *)
        torange: real; (* the higher part of the feature *)
      desiredline: integer; (* The line the feature wants to be printed on.
                            Zero means any value.  Line increases down the
                            page. *)
      (* There can be several features that this feature points to,
         From multiscan, Gap features have two and total features have
         at least two.  So others needs to be a list of features.
         This cannot be done on the feature list itself because that
         would wreck the feature list structure.  So the feature
         carriers are used.  Likewise, a single feature, such as
         an E. coli -35 promoter walker, can have several total features
         that point to it, representing several -10s.  So there can be
         several mother features. *)
      others:  featurecarrierptr; (* other features this feature points to! *)
      mothers: featurecarrierptr; (* The features that point to this feature *)
      (* ---------------------------- *)
      (* If this feature is duplicated, the next two pointers
      indicate the previous and next feature *)
      duplicateprevious: featuretypeptr;
{
      duplicatenext: featuretypeptr; (* making the chain go forward
        means we must search the entire list in featureinlist.
        This does not seem worthwhile yet. *)
}
{ppp}
      usedInGroup: boolean; (* feature has been used in the group *) 
      (* ---------------------------- *)
      next: featuretypeptr; (* the next feature recorded *)
   end;
(* end module definition.feature.type version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

var

   genomebook, scanfeatures, (* file used by this program *)
   orfp, (* file used by this program *)
   orffeatures: text; (* file used by this program *)

(* begin module lister.var *)
   groupnumber: integer; (* counted group number *)
(* end module lister.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 = 13.24; (@ of lister.p 2016 Oct 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* ************************************************************************ *)
(* end module package.brpiece version = 7.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* ************************************************************************ *)
(* end module package.getpiece version = 7.70; {of delmod.p 2005 Sep 15} *)

(* 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 = 13.24; (@ of lister.p 2016 Oct 08 *)

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

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module copystring *)
procedure copystring(a: string; var b: string);
(* copy string a to b.  The length, characters and
where the string points to are all copied. *)
var l: integer; (* index to the string *)
begin
   b.length := a.length;
   for l := 1 to a.length do b.letters[l] := a.letters[l];
   b.current := a.current;
   b.next := a.next;
end;
(* end module copystring version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module equalstring *)
function equalstring(a, b: string): boolean;
(* Test for equality between two strings at current positions.
NOTE:  A compiler bug results if one directly tests this way:
   if thedefinition^.nametag = aname
The reason is completely not clear!  I showed that the parts of the
strings were identical, but the whole was not by this test.
For this reason it is *critical to test strings with equalstring.
*)
var   index: integer; (* index to both strings *)
      equal: boolean; (* are letters in a and b the same? *)
begin (* equalstring *)
      if a.length = b.length
      then begin
         index := 1;
         repeat
            equal := (a.letters[index] = b.letters[index]);
            index := succ(index)
         until (not equal) or (index > a.length);
         equalstring := equal
      end
      else equalstring := false
end; (* equalstring *)
(* end module equalstring version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module writeNOTABstring *)
procedure writeNOTABstring(var tofile: text; var s: string);
(* write the string s to file tofile, no writeln, remove
all tab characters first. *)
var   notabs: string; (* no tab version of s *)
begin (* writeNOTABstring *)
   copystring(s,notabs);
   detabstring(notabs);
   writestring(tofile,notabs);
end; (* writeNOTABstring *)
(* end module writeNOTABstring version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

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

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module getriblmatrix *)
procedure getriblmatrix(var afile: text; var matrix: ribltypeptr);
(* Get the matrix from a file, with the defining coordinate limits,
followed by the mean and standard deviation and other data.
The entire matrix is read in. *)
const
   clearvalue = 5.0; (* If a large bit position is ever seen (for DNA), we
      KNOW there is a bug in the code! *)
var
   b: base; (* a base in the matrix *)
   l: integer; (* a coordinate in the matrix *)
   thel: integer; (* the l value in the file (should match l) *)

   teststring: string; (* string to test for 'ribl ' at matrix start *) 
   riblstart: string; (* word at the start of the matrix *)

procedure testforend;
(* test for end of file *)
begin
   if eof(afile) then begin
      writeln(output,'part or all of the ribl matrix is missing');
      halt;
   end
end;
procedure skip;
(* skip comment lines *)
begin
   testforend;
   while afile^='*' do begin
     readln(afile);
     testforend
   end
end;
begin
   new(matrix);
   with matrix^ do begin

      (* read name of the matrix *)
      grabtoken(afile, riblstart);
      teststring.letters[1] := 'r'; 
      teststring.letters[2] := 'i';
      teststring.letters[3] := 'b';
      teststring.letters[4] := 'l';
      teststring.length := 4; 
 
      if not equalstring(riblstart, teststring) then begin
         writeln(output,'The ribl matrix must begin with ''ribl "name"''');   
         writequotestring(output,riblstart);
         writeln(output, ' seen instead');
         halt
      end;
 
      skipblanks(afile); (* skip any spaces *)
      readquotestring(afile, matrix^.riblname);
      readln(afile);

      readheader(afile,'*',riblheader);

      readln(afile,frombase,tobase);

      if frombase < minribl then begin
         writeln(output,'Matrix is too large for program.');
         writeln(output,'frombase (=',frombase:1,')');
         writeln(output,'minribl (=',minribl:1,')');
         writeln(output,'Decrease minribl or increase frombase');
         halt
      end;

      if tobase > maxribl then begin
         writeln(output,'Matrix is too large for program.');
         writeln(output,'tobase (=',tobase:1,')');
         writeln(output,'maxribl (=',maxribl:1,')');
         writeln(output,'Increase maxribl or decrease tobase');
         halt
      end;

      (* clear unused matrix for safety *)
      for l := minribl to frombase-1
      do for b := a to t do matrix^.data[b,l] := clearvalue;
      for l := tobase+1 to maxribl
      do for b := a to t do matrix^.data[b,l] := clearvalue;

      for l := frombase to tobase do begin
         testforend;
         for b := a to t do read(afile,matrix^.data[b,l]);
         read(afile,thel);
         if thel <> l then begin
            writeln(output, 'positions in ribl do not match from-to range');
            halt;
         end;
         for b := a to t do read(afile,matrix^.numbers[b,l]);
         readln(afile)
      end;

      skip; readln(afile,mean);
      skip; readln(afile,stdev);
      skip; readln(afile,consensus);
      skip; readln(afile,anticonsensus);
      skip; readln(afile,averageRi);
      skip; readln(afile,n);
      skip; readln(afile,symmetry);
      skip; readln(afile,Ribound);
      skip; readln(afile,Zbound);
      skip; readln(afile,Pbound);
{this is a guess - set the converters to neutral for now
      cmperbase := 1.0;
      cmperbit := 1.0;
      readwaveparameters(afile, waves, cmperbase, cmperbit);
}
      readwaveparameters(afile, waves);
   end;
end;
(* end module getriblmatrix version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module locatefeature *)
procedure locatefeature(aname: string;
                        coo: real;
                        featurelist: featuretypeptr;
                        var thefeature: featuretypeptr;
                        StopAtTotal: boolean);

(* locate the feature corresponding to aname and coordinate coo in the
featurelist.  If the feature is found, it is pointed to by thefeature.
If it is not found, thefeature is nil.  If StopAtTotal is true, stop
when the first total feature is found.  This is used to confine the
search to a particular cluster of mutipart features.  Note that the
appropriate start point (ie, featurelist) must be determined by the
callig routines to begin at the start of the cluster. *)
var
   done: boolean; (* are we done searching for the feature? *)
begin
   thefeature := featurelist;
   done := false;
   repeat
      if thefeature = nil
      then done := true
      else begin
         (* is it the one we want? *)
   {
write(output,'locatefeature: thefeature^.nametag = ');
writequotestring(output,thefeature^.nametag);
write(output,' thefeature^.coordinate =', thefeature^.coordinate:1:1);
write(output,' coo =', coo:1:1);
writeln(output);
   }
         if equalstring(thefeature^.nametag, aname)
         then begin
            if thefeature^.coordinate = coo
            then done := true
            else thefeature := thefeature^.next
         end
         else if thefeature^.definition^.istotal and StopAtTotal   
         then begin
             (* we found the total and not the feature *)
             thefeature := nil
         end
         else thefeature := thefeature^.next
      end
   until done;
end;
{mmm}
(* end module locatefeature version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module writeadefinition *)
procedure writeadefinition(var afile: text; d: definitiontypeptr);
(* write the definition d to file afile *)
var i: integer; (* index to the string s *)
begin
   if d = nil then begin
      writeln(output,'writeadefinition was called with a nil definition!');
      halt;
   end;
   with d^ do begin
      write(afile,'define');
      write(afile,' '); writequotestring(afile,nametag);
      write(afile,' '); writequotestring(afile,background);
      write(afile,' '); writequotestring(afile,negparts);
      write(afile,' '); writequotestring(afile,posparts);
      for i := 1 to marks do write(afile,' ',locations[i]:1:1);
      writeln(afile);
{
      (* don't write the matrix this here - it's awkward if one wants
         to just write the parts above.  Just remember to write the
         matrix after the definition! *)
      writeln(afile);
      putriblmatrix(afile, matrix);
}
   end;
end;
(* end module writeadefinition version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module showadefinition *)
procedure showadefinition(var afile: text; d: definitiontypeptr);
(* show the definition d to file afile *)
const
   pwid = 8; (* width in character places to print PostScript numbers *)
   pdec = 5; (* decimal places to print PostScript numbers *)
var i: integer; (* index to the string s *)
begin
   writeln(afile,'DEFINITION: ');
   with d^ do begin
      write(afile,'nametag: ');
      writequotestring(afile,nametag);
      writeln(afile);

      write(afile,'background: ');
      writequotestring(afile,background);
      writeln(afile);

      write(afile,'negative parts: ');
      writequotestring(afile,negparts);
      writeln(afile);

      write(afile,'positive parts: ');
      writequotestring(afile,posparts);
      writeln(afile);

      write(afile,'locations: ');
      for i := 1 to marks do write(afile,' ',locations[i]:1:1);
      writeln(afile);

      writeln(afile,'marks: ',marks:1);
      writeln(afile,'min: ',min:5:1);
      writeln(afile,'max: ',max:5:1);

      write(afile,'petal: ');
      with petal do begin
         write(afile,el:1);
         write(afile,' ', c:1);
         write(afile,' ',eh:pwid:pdec);
         write(afile,' ',es:pwid:pdec);
         write(afile,' ',eb:pwid:pdec);
         write(afile,' ',fh:pwid:pdec);
         write(afile,' ',fs:pwid:pdec);
         write(afile,' ',fb:pwid:pdec);
         write(afile,' ',spare:pwid:pdec);
      end;
      writeln(afile);

{
      write(output,'THIS IS');
      if d^.istotal
      then writeln(output,      ' a total.')
      else if d^.isgap
           then writeln(output, ' a Gap.')
           else writeln(output, ' NEITHER a Gap nor a total.');
}
{sss}

   end;
end;
(* end module showadefinition version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module showdefinitions *)
procedure showdefinitions(var afile: text; definitionlist: definitiontypeptr);
(* show the definition in list definitionlist *)
var
   f: definitiontypeptr; (* current place in the definitionlist *)
begin
   f := definitionlist;
   while f <> nil do with f^ do begin
      showadefinition(afile, f);
      writeln(afile);
      f := f^.next;
   end;
end;
(* end module showdefinitions version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module writeafeature *)
procedure writeafeature(var afile: text; f: featuretypeptr);
(* write the feature f to file afile in the format:
@ K00000 36 +1 "nametag" "othertag" 12.200338 -0.473212 0.318031 -26
where the last three numbers are the Ri, Z and p
followed by the relative coordinate [2016 aug 18].
Do not carriage return.
*)
var
   o: featurecarrierptr; (* other or mother features that this one
      points to *)
begin
{writeln(output, '======> writeafeature zwwww');}
   if f = nil then begin
      writeln(output, 'writeafeature: nil feature pointer!');
      halt;
   end;
   with f^ do begin
      write(afile,'@ ');
      writeNOTABstring(afile,id);
      write(afile,' ',coordinate:1:1);
      write(afile,' ');
      if orientation > 0 then write(afile,'+',orientation:1)
                         else write(afile,    orientation:1);
      write(afile,' ');
      writequotestring(afile,nametag);
      write(afile,' ');
      writequotestring(afile,othertag);
{write(output, '[======> writeafeature zwwww did othertag]');}
      if definition = nil then begin
{
         writeln(output); (* break away from the previous line *)
         writeln(output, 'writeafeature: nil definition pointer');
         writeln(output, 'for feature:');
         writeNOTABstring(output,id);
         write(output,' ');
         writequotestring(output,nametag);
         write(output,' ');
         writequotestring(output,othertag);
         writeln(output);
}
         write(afile,' ',0.0:infofield:infodecim);
         write(afile,' ',0.0:infofield:infodecim);
         write(afile,' ',0.0:infofield:infodecim);
         write(afile,' ',0.0:3:1); (* relative coordinate *)
      end
      else begin
         if definition^.matrix <> nil then begin
{
            write(afile,' ',Ri:infofield:infodecim);
            write(afile,' ',Z:infofield:infodecim);
            write(afile,' ',probability:infofield:infodecim);
}
{write(output, '[======> writeafeature zwwww about to Ri]');}
            write(afile,' ',Ri:1:infodecim);
            write(afile,' ',Z:1:infodecim);
            write(afile,' ',probability:1:infodecim);
         end
         else begin (* parameters, new as of 2004 July 17 *)
            write(afile,' ',Tparam);
{
            write(afile,' ',Aparam:infofield:infodecim);
            write(afile,' ',Bparam:infofield:infodecim);
            write(afile,' ',Cparam:infofield:infodecim);
            write(afile,' ',Dparam:infofield:infodecim);
}
{write(output, '[======> writeafeature zwwww about to Aparam]');}
            write(afile,' ',Aparam:1:infodecim);
{write(output, '[======> writeafeature zwwww about to Bparam]');}
            write(afile,' ',Bparam:1:infodecim);
            write(afile,' ',Cparam:1:infodecim);
            write(afile,' ',Dparam:1:infodecim);
         end;
         write(afile,' ',relcoo:1:1); (* relative coordinate *)
{zzzqqq}
         write(afile,' #',featurenumber:1); (* feature number *)
         write(afile,' in group ',featuregroup:1); (* feature group *)
      end;
      if others <> nil then begin
         if others^.afeature <> nil then write(afile,' > ');
         o := others;
         while o <> nil do begin
            if o <> others then write(afile,', ');
            if o^.afeature <> nil then begin
               write(afile,'[');
               writeNOTABstring(afile,o^.afeature^.nametag);
               write(afile,']');
            end;
            o := o^.next;
         end;
      end;
   end
end;
(* originally from 9.29; (@ of lister.p 1999 July 5 *)
(* end module writeafeature version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module showafeature *)
procedure showafeature(var afile: text; f: featuretypeptr);
(* show the feature f to file afile *)
var
   o: featurecarrierptr; (* other features that this one points to *)
   totals: integer; (* totals pointing to f *)
   gaps: integer; (* gaps pointing to f *)
begin
   totals := 0;
   gaps := 0;
   with f^ do begin
      write(afile,'FEATURE #', featurenumber:1,
                  ' in group ',featuregroup:1);
{zzzfff}
      if featureisfirstingroup
      then write(afile,', FIRST in this group')
      else begin
         write(afile,', NOT FIRST in this group');
         if f^.definition^.istotal
         then write(afile,', TOTAL = end of this group')
         else write(afile,', MIDDLE feature');
      end;
      writeln(afile);
      write(afile,'id: ');
      writequotestring(afile,id);
      writeln(afile);

      writeln(afile,'coordinate: ',coordinate:1:1);

      write(afile,'orientation: ');
      if orientation > 0 then writeln(afile,'+',orientation:1)
                         else writeln(afile,    orientation:1);

      write(afile,'nametag: ');
      writequotestring(afile,nametag);
      writeln(afile);

      write(afile,'othertag: ');
      writequotestring(afile,othertag);
      writeln(afile);

      write(afile,'unsatisfied: ');
      if unsatisfied then writeln(afile,'true')
                     else writeln(afile,'false');

      writeln(afile,'fromrange: ',fromrange:5:1);
      writeln(afile,'  torange: ',  torange:5:1);
      writeln(afile,'desiredline: ',  desiredline:5);
      if duplicateprevious = nil
      then writeln(afile,'NOT a duplicate feature')
      else writeln(afile,'DUPLICATE of feature #',
                   duplicateprevious^.featurenumber:1);

      if others <> nil then begin
         writeln(afile,'This feature points at these OTHERS:');
         write(afile,'feature #',featurenumber:1,' -> ');
         o := others;
         while o <> nil do begin
{
write(afile,'-LOOP-');
if o = others then write(afile, '|o = others');
if o <> others then write(afile, '|o <> others');
if not (o = others) then write(afile, '|not (o = others)');
write(afile, '|');
}
{KKKKK}
            if o <> others then write(afile,', ');
            if o^.afeature <> nil then begin
               write(afile,'[');
               writeNOTABstring(afile,o^.afeature^.nametag);
               write(afile,']');
            end
            else write(afile,'[EMPTY]');
            o := o^.next;
         end;
         writeln(afile);
         o := others;
         while o <> nil do begin
            if o^.afeature <> nil then begin
               write(afile,'= ');
               writeafeature(afile,o^.afeature);
               writeln(afile);
            end;
            o := o^.next;
         end;
      end
      else begin
         writeln(afile,'This feature does not point to any OTHERS.');
      end;

      if mothers <> nil then begin
         writeln(afile,'The MOTHERS of this feature', 
                       ' (that point to this feature) are:');
         o := mothers;
         while o <> nil do begin
            if o^.afeature <> nil then begin
               write(afile,'feature #',featurenumber:1,' <- "');
               writeNOTABstring(afile,o^.afeature^.nametag);
               writeln(afile,'" #',o^.afeature^.featurenumber:1);
{
               writeln(afile,'"');
}
               if o^.afeature^.definition^.istotal
               then totals := succ(totals);
               if o^.afeature^.definition^.isgap
               then gaps := succ(gaps);
            end;
{writeln(output,'BUBBA: mothers part 1 about to "o := o^.next;" ');}
            o := o^.next;
         end;
         o := mothers;
{writeln(output,'BUBBA: mothers part 1 is done');}
         while o <> nil do begin
            if o^.afeature <> nil then begin
               write(afile,'= ');
               writeafeature(afile,o^.afeature);
               writeln(afile);
            end;
            o := o^.next;
         end;
{writeln(output,'BUBBA: mothers part 2 is done');}
      end
      else begin
         writeln(afile,'This feature has no MOTHERS.');
      end;
{writeln(output,'BUBBA: showafeature is done');}

   end;
   write(afile,gaps:1,' gap');
   if gaps <> 1 then write(afile,'s'); (* stolen from plural *)
   writeln(afile);
   write(afile,totals:1,' total');
   if totals <> 1 then write(afile,'s'); (* stolen from plural *)
   writeln(afile);
   if (totals > 1) then begin
      writeln(output,'showafeature: more than one total!');
      halt;
   end;
   if (gaps >1) and (totals < 1) then begin
      writeln(output,'showafeature: ',gaps:1,' gaps but no total!');
{NNN}
      halt;
   end;
end;
(* end module showafeature version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module showfeatures *)
procedure showfeatures(var afile: text; featurelist: featuretypeptr;
                       paging: boolean);
(* Show all of the features in list featurelist.
If paging is true, page between the features
*)
var
   f: featuretypeptr; (* current place in the featurelist *)
begin
   f := featurelist;
   while f <> nil do with f^ do begin
      showafeature(afile, f);
      f := f^.next;
      if f <> nil then begin
         (* don't add spacing or page for the last one! *)
         if paging
         then page(afile)
         else writeln(afile);
      end;
   end;
end;
(* end module showfeatures version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

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

(* begin module readadefinition *)
procedure readadefinition(var afile: text; var d: definitiontypeptr);
(* read a single definition from afile *)
var
   i: integer; (* an index *)
   location: real; (* a location on the sequence *)

procedure IsGapOrIsTotal(var d: definitiontypeptr);
(* Is this a definition which is for a Gap or a total connecting
multi-walkers?  If so, set isgap or istotal true (respectively).  *)
const
   tab = 9; (* tab characters separate name parts *)
   debugging = false; (* use to debug the code *)
var
   totstring: string; (* the string that indicates a IsGapOrIsTotal total *)
   gapstring: string; (* the string that indicates a IsGapOrIsTotal gap *)
   deflocation: integer; (* location in the nametag *)
   gaplocation: integer; (* location in the gapstring *)
   totlocation: integer; (* location in the total *)
   done: boolean; (* done comparing strings *)
begin

if debugging then begin
  writeln(output,'IsGapOrIsTotal test of this definition:'); {sss}
  writeadefinition(output,d); {sss}
end;

   clearstring(totstring);
   with totstring do begin
      letters[1] := chr(tab);
      letters[2] := 't';
      letters[3] := 'o';
      letters[4] := 't';
      letters[5] := 'a';
      letters[6] := 'l';
      length := 6;
   end;

   clearstring(gapstring);
   with gapstring do begin
      letters[1] := chr(tab);
      letters[2] := 'G';
      letters[3] := 'a';
      letters[4] := 'p';
      length := 4;
   end;

{
writeln(output, 'IsGapOrIsTotal: begin &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
}

   (* search for the Gap string *)
   done := false;
   d^.isgap := true;
   deflocation := d^.nametag.length;
   gaplocation := gapstring.length;
   while not done do begin
{
write(output,
'====>',
d^.nametag.letters[deflocation],
gapstring.letters[gaplocation]);
}
      if d^.nametag.letters[deflocation] <> gapstring.letters[gaplocation]
      then begin
         d^.isgap := false;
         done := true
      end;
      deflocation := deflocation - 1;
      gaplocation := gaplocation - 1;
{
write  (output,' deflocation = ',deflocation:1);
writeln(output,' gaplocation = ',gaplocation:1);
}
      if deflocation < 1 then done := true;
      if gaplocation < 1 then done := true;
   end;

(* just for debugging:
   if d^.isgap then begin
      if debugging then write(output, 'gap:   ');
      write(output, '====================');
      writestring(output, d^.nametag);
      writeln(output, '=======');
      writeadefinition(output, d);
      writeln(output);
      writeln(output, '====================');
{ok pointing to the features of this
has to happen with each feature of this type }
{      d^.isgap := true why set this???}
   end
   else begin
{
      if debugging then write(output, '       ');
      writestring(output, d^.nametag);
      writeln(output);
}
{      d^.isgap := false why set this???}
   end;
*)

   (* search for the total string *)
   d^.istotal := true;
   done := false;
   (* skip to the tab before the bits *)
   deflocation := d^.nametag.length;
   while not done do begin
{
write(output,d^.nametag.letters[deflocation]);
}
      if d^.nametag.letters[deflocation] = chr(tab)
      then done := true;
      deflocation := deflocation - 1;
      if deflocation < 1 then begin
         done := true;
         d^.istotal := false; (* fell off begin of string *)
      end;
   end;
{
writeln(output,'<-PRETOTAL');
}
   if d^.istotal = true then begin (* it still might be true *)
{
writeln(output,'...');
}
      totlocation := totstring.length;
      done := false;
{
this is defined above!
      deflocation := d^.nametag.length;
}
      while not done do begin
{
write(output,
'====>',
d^.nametag.letters[deflocation],
totstring.letters[totlocation]);
}
         if  d^.nametag.letters[deflocation]
         <> totstring.letters[totlocation]
         then begin
            d^.istotal := false;
            done := true
         end;
         deflocation := deflocation - 1;
         totlocation := totlocation - 1;
{
write  (output,' deflocation = ',deflocation:1);
writeln(output,' totlocation = ',totlocation:1);
}
         if deflocation < 1 then done := true;
         if totlocation < 1 then done := true;
      end;
   end;

   (* mark that this definition is part of a multiple set *)
   if d^.isgap or d^.istotal
   then d^.ismulti := true
   else d^.ismulti := false;

   if debugging then begin
{sss}
{
      writeadefinition(output, d);
}
      showadefinition(output, d);
{
if d^.istotal then begin
writeln(output, 'IsGapOrIsTotal: is total');
halt;
end;
}
      if d^.isgap
      then writeln(output,'Is a gap')
      else writeln(output,'Is NOT a gap');
      if d^.istotal
      then writeln(output,'Is a total')
      else writeln(output,'Is NOT a total');
      writeln(output, 'IsGapOrIsTotal: end &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
   end;

end; (* IsGapOrIsTotal *)

begin (* readadefinition *)
{writeln(output,'readadefinition');}
   skipnonblanks(afile);
   skipblanks(afile);
   with d^ do begin
      readquotestring(afile, nametag);

      (* read background character from a quotestring *)
      readquotestring(afile,background);
      if (background.length < 1) or
         (background.length > 2) then begin
         writeln(output,'Background string must contain only 1 or 2',
                        ' characters.');
         write(output,'Instead, it contained: ');
         writequotestring(output,background);
         writeln(output);
         halt
      end;

      readquotestring(afile,negparts);
      readquotestring(afile,posparts);
      if negparts.length <> posparts.length then begin
         write(output,'While reading definition of ');
         writequotestring(output,nametag);
         writeln(output);
         write(output,'the negative part ');
         writequotestring(output,negparts);
         writeln(output,' and');
         write(output,'the positive part ');
         writequotestring(output,posparts);
         writeln(output,' are not the same length.');
         writeln(output,'(',negparts.length:1,'<>',posparts.length:1,')');
         halt
      end;

      marks := negparts.length;
      if marks = 0 then begin
         (* multiscan uses empty strings sometimes ... *)
         marks := 2
      end;

      for i := 1 to marks do read(afile,locations[i]);

      matrix := nil;
      RiboundFeature := -maxint;
      ZboundFeature := +maxint;
      PboundFeature := 0.0;

      ismulti := false;

      next := nil;

      (* checks *)

      (* now compute the absolute coordinates that the feature resides in *)
      min := +maxint;
      max := -maxint;
      for i := 1 to marks do begin
         location := locations[i];
         if location < min then min := location;
         if location > max then max := location;
      end;
      if min = +maxint then begin
         writeln(output,'readadefinition: program error at min');
         halt
      end;
      if max = -maxint then begin
         writeln(output,'readadefinition: program error at max');
         halt
      end;

      (* determine if this is a Gap or total type feature *)
      IsGapOrIsTotal(d);

   end;
{
writeln(output,'-----');
showadefinition(list,d);
writeln(output,'-----');
}
end;
(* end module readadefinition version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module cleardefinition *)
procedure cleardefinition(var d: definitiontypeptr);
(* clear out the definition d *)
var
   l: integer; (* index to locations *)
begin
   with d^ do begin
      clearstring(nametag);
      clearstring(background);
      clearstring(negparts);
      clearstring(posparts);
      for l := 1 to maxstring do locations[l] := 0;
      marks := 0;
      min := 0.0;
      max := 0.0;
      definitionnumber := 0;
      matrix := nil;
      RiboundFeature := 0.0;
      ZboundFeature := 0.0;
      PboundFeature := 0.0;

      with petal do begin
         el := 0;
         c  := ' '; (* blank means NO PETAL *)
         eh := 0.0;
         es := 0.0;
         eb := 0.0;
         fh := 0.0;
         fs := 0.0;
         fb := 0.0;
         spare := 0.0;
      end;

      isgap := false;
      istotal := false;
      ismulti := false;
      next := nil;
   end
end;
(* end module cleardefinition version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module equaldefinition *)
function equaldefinition(a,b: definitiontypeptr): boolean;
(* compare the definitions *)
var
   same: boolean; (* are the two definitions EXACTLY the same? *)
   i: integer; (* index to marks *)
begin
   same := true;
   if equalstring(a^.nametag,        b^.nametag) then begin
      if equalstring(a^.background,  b^.background) then begin
         if equalstring(a^.negparts, b^.negparts) then begin
{I'm not sure the next two are good to do: }
            if a^.min = b^.min then begin
               if a^.max = b^.max then begin
   (* definitely not to be compared:
   if a^.number = b^.number then begin
   *)
                  if a^.marks = b^.marks then begin
                     for i := 1 to a^.marks do begin
                        if a^.locations[i] <> b^.locations[i]
                        then same := false
                     end
                  end
               end
               else same := false
            end
            else same := false
         end
         else same := false
      end
      else same := false
   end
   else same := false;
   equaldefinition := same
{
      for i := 1 to marks do write(afile,' ',locations[i]:1:1);
not worth the time savings:
      i := 1;
      while same and (a^.locations[i] = b^.locations[i]) do begin
         i := i + 1
      end;
}
end;
(* end module equaldefinition version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module clearfeature *)
procedure clearfeature(var f: featuretypeptr);
(* clear out the feature f *)
begin
   with f^ do begin
      clearstring(id);
      coordinate := 0;
      orientation := +1;
      clearstring(nametag);
      clearstring(othertag);
      definition := nil;
      duplicateprevious := nil;
      Ri := 0.0;
      Z := 0.0;
      probability := 0.0;
      Tparam := ' ';
      Aparam := 0.0;
      Bparam := 0.0;
      Cparam := 0.0;
      Dparam := 0.0;
      relcoo := 0.0;
      evencoordinate := 0.0;
      unsatisfied := true;
      fromrange := 0;
      torange := 0;
      featurenumber := 0;
      featuregroup := 0;
      featureisfirstingroup := false;
      desiredline := 0;
      others := nil;
      mothers := nil;
      usedInGroup := false; (* used in group *)
      next := nil;
   end
end;
(* end module clearfeature version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module readafeature *)
procedure readafeature(var afile: text;
                           definitionlist: definitiontypeptr;
                       var f: featuretypeptr);
(* Read a single feature from afile.  It has to match one of the
definitions in definitionlist.  The LAST definition read so far is the
one used for this feature, not any previous ones.
No carriage return is performed; this allows comments at the end of
the feature.
*)
var
   d: definitiontypeptr; (* current place in the definitionlist *)
procedure reportandhalt;
(* report a bad feature and halt *)
begin
   writeln(output, 'Feature details:');
   writeafeature(output, f);
   writeln(output);
   halt;
end;
begin
{writeln(output,'readafeature');}
   with f^ do begin
      skipnonblanks(afile);
      skipblanks(afile);
      grabtoken(afile, id);
      read(afile,coordinate);
      read(afile,orientation);
      readquotestring(afile,nametag);
      readquotestring(afile,f^.othertag);

      (* checks *)
      if not ((orientation = -1) or (orientation = +1)) then begin
         write(output,'while reading feature id ');
         writequotestring(output,id);
         writeln(output);
         writeln(output, 'feature orientation must be -1 or +1');
         writeln(output, 'It is ',orientation:1);
         reportandhalt
      end;

      (* Use the last definition read in before this feature *)
      (* 2008 Aug 06 *)
      d := definitionlist;
      while d <> nil do begin
         if equalstring(d^.nametag, f^.nametag) then begin
            definition := d;
         end;
         d := d^.next
      end;
      if definition = nil then begin
         write(output, 'No definition was found for feature ');
         writequotestring(output,nametag);
         writeln(output);
         reportandhalt
      end;
      d := definition;

      (* now we can read in the Ri, Z and probability *)
      if d^.matrix <> nil
      then begin
         read(afile, Ri, Z, probability);
         if not eoln(afile)
         then read(afile,relcoo) (* 2016aug19: new parameter *)
         else relcoo := 0.0; (* protect for reading old features *)
      end
      else begin
         (* give values that allow the site to always be found *)
         Ri := -maxint;
         Z  := maxint;
         probability := 1.0;
         (* read parameters for feature: *)

         (* To keep the previously defined features stable,
            read T param here by first making sure that the
            next object is not a number. *)

         skipblanks(afile);
         if eoln(afile) then begin (* handle blank line *)
            Tparam := ' ';
            Aparam := 0.0;
            Bparam := 0.0;
            Cparam := 0.0;
            Dparam := 0.0;
            relcoo  := 0.0;
         end
         else begin
            if (afile^ in
               ['-','+','0','1','2','3','4','5','6','7','8','9'])
            then begin
               (* Old style, just read it in *)
               read(afile, Aparam);
               read(afile, Bparam);
               read(afile, Cparam);
               Dparam := 0;
            end
            else begin
               read(afile, Tparam);
               read(afile, Aparam);
               read(afile, Bparam);
               read(afile, Cparam);
               read(afile, Dparam);
           end;
           if not eoln(afile)
           then read(afile,relcoo) (* 2016aug19: new parameter *)
           else relcoo := 0.0; (* protect for reading old features *)
        end;
      end;

(* 2004 July 17: I went ahead and removed this.  As far as I know
it is NOT used in feature files:
      if not eoln(afile) then read(afile, evencoordinate)
                         else evencoordinate := coordinate;
*)

      unsatisfied := true;
      fromrange := definition^.min;
        torange := definition^.max;

{
      featurenumber := globalfeaturenumber; (* number assigned *)
}
      desiredline := 0; (* initially it can go on any line *)
      next := nil; (* it is not connected to anything yet *)

      (*
      others:  featurecarrierptr; (@ other features this feature points to! @)
      mothers: featurecarrierptr; (@ The features that point to this feature @)
      *)
      others := nil;
      mothers := nil;

   end
end;
(* originally from 9.29; (@ of lister.p 1999 July 5 *)
(* end module readafeature version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module featureinlist *)
function featureinlist(featurelist: featuretypeptr;
                        newfeature: featuretypeptr;
                    var oldfeature: featuretypeptr): boolean;
(* determine if the new feature is found in the feature list.
If it was found, report the location in oldfeature.
Note that this is a SIDE EFFECT. *)
{2007 Mar 26 ppp}
var
   f: featuretypeptr; (* a feature in featurelist *)
{
   test: featuretype; (* a test feature.  We use this to
      avoid modifying newfeature *)
}
   done: boolean; (* done searching? *)
   found: boolean; (* the feature matches *)
procedure checkorpos;
(* check the orientation and position *)
begin
{
writeln(output,'checkorpos');
}
   found :=
      (newfeature^.orientation = f^.orientation)
         and
      (newfeature^.coordinate = f^.coordinate)
end; (* checkorpos *)
begin

{ just some checks for debugging
writeln(output);
writeln(output,'----------------------------------------------------');
showfeatures(output, featurelist, true);
writeln(output,'--------');
write(output,'new: ');
writeafeature(output,newfeature);
writeln(output);
}

   done := false;
   found := false;
   f := featurelist;
   oldfeature := nil;
{
write(output,'newfeature: ');
writeafeature(output,newfeature);
writeln(output);
}
   while not done do begin

{
write(output,'CHECK:   f: ');
writeafeature(output,f);
writeln(output);
}

      if equalstring(newfeature^.id, f^.id) and
         equalstring(newfeature^.nametag, f^.nametag)
      then if f^.definition = nil
         then checkorpos
         else if f^.definition^.matrix = nil
              then checkorpos
              else found :=

         (* Be paranoid about losing a walker!!
            They have to be the same evaluation: *)
         (newfeature^.Ri = f^.Ri) and
{
this would be overkill that probably would'nt help:
 and
         (newfeature^.Z = f^.Z) and
         (newfeature^.probability = f^.probability) and
}

         (* all three cases of matrix symmetry must be considered. *)

         (

           (* an odd symmetry will always match. *)
           (
             (f^.definition^.matrix^.symmetry = 'o')
               and
             (newfeature^.coordinate = f^.coordinate)
           )

           or

           (* an asymmetric symmetry will only match with = orientation. *)
           (
             (f^.definition^.matrix^.symmetry = 'a')
             and
             (newfeature^.orientation = f^.orientation)
             and
             (newfeature^.coordinate = f^.coordinate)
           )

           or
(* it doesn't seem possible to use the evencoordinate variable here!
Let it hang around though just in case.  The difficulty is that
at a coordinate jump evencoordinate will be different on the two sides.
*)

           (* an even symmetry will only match if one base away. *)
           (* THIS WILL FAIL AT A COORDINATE JUMP.  (fortunately rarely) *)
           (
             (f^.definition^.matrix^.symmetry = 'e')
             and
             (1 = abs(newfeature^.coordinate - f^.coordinate))
             and
             (newfeature^.orientation <> f^.orientation)
           )

           or

           (* an even symmetry will also match if identical. *)
           (
             (f^.definition^.matrix^.symmetry = 'e')
             and
             (newfeature^.coordinate = f^.coordinate)
             and
             (newfeature^.orientation = f^.orientation)
           )

         );

      if found then begin
         oldfeature := f;
         done := true;
{
;

write(output,'newfeature: ');
writeafeature(output,newfeature);
writeln(output);
write(output,'         f: ');
writeafeature(output,f);
writeln(output);

}
      end
      else if f^.next = nil
           then begin
              done := true;
              found := false
           end
           else f := f^.next
   end;
{
if found then
then writeln(output,'featureinlist: found = trew blew')
else writeln(output,'featureinlist: found = falseto');
}
   featureinlist := found
end;
(* end module featureinlist version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module plural *)
procedure plural(var thefile: text; number: integer);
(* if the number is not 1, return an s *)
begin (* plural *)
      if number <> 1 then write(thefile,'s')
end; (* plural *)
(* end module plural version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module locatedefinition *)
procedure locatedefinition(aname: string;
                           definitionlist: definitiontypeptr;
                           var thedefinition: definitiontypeptr);
(* locate the definition corresponding to aname
in the definitionlist. *)
var
   done: boolean; (* are we done searching for the definition? *)
begin
   thedefinition := definitionlist;
   done := false;
   repeat
      if thedefinition = nil
      then done := true
      else begin
         (* is it the one we want? *)
         if equalstring(thedefinition^.nametag, aname)
         then done := true
         else thedefinition := thedefinition^.next
      end
   until done;
end;
(* end module locatedefinition version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module readdefinitionsandfeatures *)
procedure readdefinitionsandfeatures(var afile: text;
                       var KeepDuplicates: boolean;
{zzzccc}
                       var comments: stringptr;
                       var definitionlist: definitiontypeptr;
                       var featurelist: featuretypeptr;
                       var groupnumber: integer);
(* Read definitions and features from afile.
These are defined in the lister.p program.
Keep duplicate features if KeepDuplicates is true.

The objects in the afile are identifiers.
The first character of lines in afile defines the identifier.
The legal identifier values are:

* comment
d definition of a feature
r ribl matrix
b boundary for site
c color rectangle
@ location of a feature
g beginning of a group that has gap connection lines and ends with a total line

2007 Mar 20.  One of the important tasks of readdefinitionsandfeatures
is to string together the features of multi-part features coming from
multiscan.  Traditionally this was done by searching for the related
parts while the list was being read and through all previously read
features.  However, this tends to lead to confusion in very complex
genetic structures where several sites of even the same kind can be in
the same position.  To avoid this the program (version 11.90 and on)
now keeps track of a new variable, cluster, which is a pointer into
the featurelist.  As each cluster is found, searches to find the parts
in markismulti are restricted to the cluster.  Also, all the threading
is done AFTER the features are read in, simplifying the process.
222

2007 Mar 26. The problem with actually removing duplicate features
from the featurelist is that the threading of mothers and others
fails.  The solution is to KEEP all duplicates, but to flag duplicates
(if KeepDuplicates is true) and then later to simply not count them
when displaying the results.

2016 Aug 12.  The problem with the above method is that it
relies on external coordinates to locate features.  This can
fail when going across coordinate bound of a circular piece.
New (tentative) method is to mark the multipart groups
explicitly with a 'group' line and then the threading should
just join the parts directly instead of clustering later.

2016 Sep 06: I removed all the code for markismulti

*)
var
   aname: string; (* a name of the ribl feature being read in *)
   commentparts: stringptr; (* comment parts *)
   countdefinitions: integer; (* count of the definitions read in *)
   countfeatures: integer; (* count of the features read in *)
   countribls: integer; (* count of the ribls read in *)
   countboundaries: integer; (* count of the boundaries read in *)
   countpetals: integer; (* count of the color rects read in *)
   currentcomment: stringptr; (* pointer to the currently read comment *)
   d: definitiontypeptr; (* current place in the definitionlist *)
   done: boolean; (* are we done searching for the definition? *)
   dprevious: definitiontypeptr; (* place in the definitionlist prior to d *)
   dtest: definitiontypeptr; (* place in the definitionlist for testing *)
   emptyafile: boolean; (* true if the afile is empty *)
   f: featuretypeptr; (* current place in the featurelist *)
   firstfeature: boolean; (* first feature of a feature list must be
                             read properly *)

(*gggGGG*)
   (* Group definitions 2016 Aug 12
      These variables are used to thread together a group of
      features and definitions with the total line.
      Most could be moved to procedure threadgroup.
    *)
   groupnumberinfile: integer; (* the group number read from the file *)
   mycharacter: char; (* for reading past the word 'group' in afile *)
   firstfeatureofgroup: boolean; (* for identifying the first feature
      of a group [this is not the firstfeature of a feature list!] *)
   ingroup: boolean; (* we are inside a group within the featurelist
      linked list of features *)
   groupfeature: featuretypeptr; (* the start of a group feature *)
   aGapFeature: featuretypeptr; (* the location of a gap feature *)
   groupTotalFeature: featuretypeptr; (* the location of the total feature *)
   groupGapCarrier: featurecarrierptr; (* the end of the Gap list *)
   groupTotalCarrier: featurecarrierptr; (* the end of the Total list *)
   leftgap:  featuretypeptr; (* left  side of a gap *)
   rightgap: featuretypeptr; (* right side of a gap *)
   ff: featuretypeptr; (* current place in the featurelist after the groupfeature:
                          then the feature we need to thread *)
   gg: featuretypeptr; (* current place in the featurelist after the groupfeature:
                          then the group we need to thread *)
   nameparts: stringptr; (* name parts of a total string *)
   anamepart: stringptr; (* one name in a total string *)
   leftstring  : string; (* left side of gap string *)
   gapstring   : string; (* middle gap string *)
   rightstring : string; (* right side of gap string *)
   gap: integer; (* gapstring converted to a number *)
   aletter: char; (* a letter in the gapstring *)

   identifier: char; (* identify the data type in afile *)
   matrix: ribltypeptr; (* ribl matrix being read in *)
   missingmatrix: boolean; (* has the matrix for a definition been read in? *)
   newfeature: featuretypeptr; (* a new feature.  If it is unique,
      we will add it to the f list. *)
   oldfeature: featuretypeptr; (* a feature that matches the newfeature *)
   thedefinition: definitiontypeptr; (* index to locate the definition
      corresponding to a ribl feature read in. *)
   uniquedefinitions: integer; (* count of the unique definitions read in *)
   uniquefeatures: integer; (* count of the unique features read in *)
procedure yell;
(* tell the user something is wrong and who we are *)
begin
   writeln(output,'ERROR in reading definitions and features');
end;

procedure barout;
(* just write a bar for warnings *)
begin
writeln(output, '***********************************',
                '***********************************');
end;

procedure splitstring(s: string; break, terminator: char;
                      var parts: stringptr);
(* split the string s by the break characters up to the terminator
character.  The resulting parts are linked strings.  This is for
use in splitting the total nametag into its parts.  Example:
"distal-proximal-p35-p10	4164377	total	15.0"
becomes:
"distal"
"proximal"
"p35"
"p10"
*)
const
   debugSplitstring = true; (* set to true to debug *)
{TTT}
var
   l: integer; (* location in the string s *)
   done: boolean; (* done finding parts *)
   spointer: stringptr; (* pointer to the string *)
begin
if superdebug then begin
if debugSplitstring then begin
writeln(output,'IN splitstring the string to split is:');
write(output,'"');
writestring(output,s);
writeln(output,'"');
writeln(output,'======> break character is: "',break,'"');
end;
end;
   l := 1;
   done := false;
   disposestring(parts); (* clear out any previous parts *)
   new(parts);
   spointer := parts;
   clearstring(spointer^);
   while not done do begin
      if s.letters[l] = terminator then done := true
      else begin
         if s.letters[l] = break then begin
if superdebug then begin
if debugSplitstring then begin
write(output,'debugSplitstring: "');
writestring(output,spointer^);
writeln(output,'"');
end;
end;
            new(spointer^.next);
            spointer := spointer^.next;
            clearstring(spointer^); {TTT}
         end
         else begin
            spointer^.length := spointer^.length + 1;
            spointer^.letters[spointer^.length] := s.letters[l];
if superdebug then begin
if debugSplitstring then begin
if (s.letters[l] = ' ') then begin
writeln(output,'debugSplitstring: splitstring write a space $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$');
end;
if (s.letters[l] = chr(9)) then begin
writeln(output,'debugSplitstring: splitstring write a TAB  $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$');
end;
if false then (* silence the verbose stuff for now *)
writeln(output,'WRITING: ord(', s.letters[l],')=', ord(s.letters[l]):1);
end;
end;
         end;
         l := succ(l);
         if l > maxstring then done := true
      end;
      if l = maxstring then begin
         writeln(output,'splitstring: maxed out! maxstring = l = ',l:1);
         halt;
      end;
   end;
if superdebug then begin
if debugSplitstring then begin
write(output,'debugSplitstring: splitstring result: "'); 
writestring(output,spointer^);
writeln(output,'"');
writeln(output,'OUT splitstring');
end;
end;
end; (* splitstring *)

procedure threadgroup;
(* Thread the current feature f and other features into a group - if
we find that we are at the end of the group. *)
var
   m: integer; (* index to marks of a total feature *)

(*-----*)
procedure countmarks(f:featuretypeptr; requiredcount: integer);
(* Count the others in feature f.
If the observed number is not the expected number (expectedcount),
then something is very wrong!
If requiredcount is -1, then the expected count is
the number of marks in the total.  Otherwise requiredcount is
the count expected.
*)
const
   countmarkdebug = false; (* display extra information when called *)
var
   expectedcount: integer; (* expected number of others *)
   observedcount: integer; (* count of others in a total feature *)
   otherptr: featurecarrierptr; (* location of an other *)
procedure countmarksprint; (* print the feature information *)
begin
   write(output,'++++++++ countmarks ++++++++');
   if requiredcount < 0
   then writeln(output,' final check')
   else writeln(output,' looking for ',requiredcount:1, ' marks');
   writeln(output,'++++++++ Counting marks for this total feature',
                  ' in group number: ',groupnumber:1,
                  ' FAILED: ');
   showadefinition(output,f^.definition);
   {writeafeature(output,f); writeln(output);}
   showafeature(output,f);
   writeln(output,f^.definition^.marks:1, ' marks are in this feature');

if false then begin (* final display for debugging *)
writeln(output);
writeln(output,'******* countmarks readdefinitionsandfeatures: *******');
writeln(output,'******* DEFINITIONS:');
if false then showdefinitions(output,definitionlist);
writeln(output);
writeln(output,'******* FEATURES:');
if true then showfeatures(output,featurelist, true);
writeln(output);
writeln(output,'******* countmarks readdefinitionsandfeatures: *******'); {sss}
{
halt;
}
end;

end; (* countmarksprint *)
begin
if countmarkdebug then writeln(output,'+++++++++++++++++++++ BEGIN countmarks');
   if (f^.definition^.istotal) then begin
      if f^.others <> nil then begin
         otherptr := f^.others;
         observedcount := 0;
         while otherptr <> nil do begin
if countmarkdebug then begin
writeln(output);
writeafeature(output,otherptr^.afeature);writeln(output);
writeln(output,'loop observedcount=',observedcount:1);
writeln(output,'isgap:   ',otherptr^.afeature^.definition^.isgap);
writeln(output,'istotal: ',otherptr^.afeature^.definition^.istotal);
end;
               if     (not otherptr^.afeature^.definition^.isgap)
                  and (not otherptr^.afeature^.definition^.istotal)
               then begin
                  observedcount := succ(observedcount);
               end;
            otherptr := otherptr^.next;
         end;
if countmarkdebug then writeln(output,'**** final observedcount=',observedcount:1);
{FFF}
         if requiredcount = -1
         then expectedcount := f^.definition^.marks
         else expectedcount := requiredcount;
         if observedcount <> expectedcount then begin
{fff}
            countmarksprint;
            writeln(output,'countmarks:',
                           ' expected: ', expectedcount:1,
                           ' observed: ', observedcount:1);
            writeln(output,'**** Wrong number of marks found',
                           '  in linked list for this total feature. ****');
            halt;
         end;
      end
      else begin
         if countmarkdebug then writeln(output,'countmarks no others');
      end
   end
   else begin
      writeln(output,'countmarks feature is not a total, cannot count marks');
      halt;
   end;
if countmarkdebug then writeln(output,'+++++++++++++++++++++ END countmarks');
end; (* countmarks *)
(*-----*)

procedure addmotherlinks(var agap: featuretypeptr;
                         var thefeature: featuretypeptr);
(* Thread the gap feature agap into the feature's mother list.
Don't add it if it is already in the list. *)
var
   mom: featurecarrierptr; (* pointer into the thefeature's mother list *)
   previousMom: featurecarrierptr; (* the previous p to allow adding
      to the end of the linked list *)
   gapthing: featurecarrierptr; (* holder for the new carrier *)
   done: boolean; (* done adding to the list *)
begin
{iii}
{
writeln(output,'@@@@@@@@@@@ in addmotherlinks @@@@@@@@@@@');
writeln(output,'------------ @@@@@ THREAD THIS gap:');
showafeature(output,agap);
writeln(output,'------------ @@@@@ INTO THIS feature:');
showafeature(output,thefeature);
}
   if thefeature^.mothers = nil then begin
if false then begin (* for debugging *)
      writeln(output,'************ ALL FEATURES: *************');
      writeln(output,'**** (also listed in the list file) ****');
      showfeatures(output,featurelist, false);
      writeln(output,'****************************************');
      writeln(output,'**** addmotherlinks: no mothers found for the feature:');
      showafeature(output, thefeature);
{ this is lister-specific code:
      rewrite(list);
      showfeatures(list,featurelist, true);
}
end;
{ halt; kills too soon }
{ttt}
      (* Probably this is wrong to do for gaps because we don't have
         the threading from total yet!! *)
      if not (    thefeature^.definition^.isgap
               or thefeature^.definition^.istotal)
      then halt;
   end
   else begin
{ddd}
      mom := thefeature^.mothers;
      previousMom := mom;
      done := false;
      while not done do begin
         if mom^.afeature = thefeature then begin
            done := true;
         end
         else begin
            if mom^.afeature^.definition^.istotal then begin
               (* We got to the end of the linked list,
                  so thread thefeature in front of the total *)
{writeln(output,'------------ addmotherlinks: THREADING IN');}
               new(gapthing);
               if mom = previousMom then begin
                  gapthing^.next := thefeature^.mothers;
                  thefeature^.mothers := gapthing;
               end
               else begin
                  previousMom^.next := gapthing; 
                  gapthing^.next := mom; 
               end;
               gapthing^.afeature := agap;  (* thread it in! *)
               done := true;
               if mom^.next <> nil then begin
                  writeln(output,'addmotherlinks: total points to something!');
                  halt;
               end;
            end
            else begin
               previousMom := mom;
               mom := mom^.next;
               if mom = nil then begin
                  writeln(output,'addmotherlinks: failed to find total');
                  halt;
               end
            end
         end;
      end;
{
writeln(output,'------------ @@@@@ adding more, result:');
showafeature(output,thefeature);
}
{
halt;
}
   end;
{
writeln(output,'------------ @@@@@ RESULT FOR THE feature:');
showafeature(output,thefeature);
}
end; (* addmotherlinks *)

begin (* threadgroup *)
   if ingroup then begin
      if groupfeature = nil then begin (* beginning of a group *)
{
writeln(output,'BEGIN OF GROUP ',groupnumber:1,' GGG');
writeafeature(output,f); writeln(output);
}
         (* this is the start of the group! *)
         groupfeature := f;
         (* this starts the group *)
      end
      else if f^.definition^.istotal then begin
         (* This total is the end of the group.
            At this point we have the entire group - since the total
            is supposed to be the end of the group.
            So the next thing to do is to thread the features
            into the total and the gap and total into the features. *)
         groupTotalFeature := f; (* use this for clarity *)
{
writeln(output);
writeafeature(output,f);writeln(output);
writeln(output);
writeln(output,'The total feature we are looking at: ========');
showafeature(output,f);
writeln(output,'The total feature''s definition: =============');
showadefinition(output,f^.definition);
}

{vvvv77}

{
writeln(output,'############## PRETHINGY');
countmarks(groupTotalFeature,-1);
writeln(output,'############## POSTTHINGY');
}
{
if groupTotalFeature^.others^.next^.next <> nil then begin
writeln(output,'***** PRE-BEFORE groupTotalFeature^.others^.next^.next <> nil');
showadefinition(output,groupTotalFeature^.definition);
showafeature(output,groupTotalFeature);
halt;
end;
}

         (* prepare to thread features into total. 
            gaps are done later based on the total. *)
         groupTotalCarrier := nil;

         (* Split out the names from the total 'nametag'. *)
{
writeln(output,'split out the total names from f^.nametag');
write(output,'"'); writestring(output,groupTotalFeature^.nametag);
writeln(output,'"');
}
         splitstring(groupTotalFeature^.nametag, '-', chr(9), nameparts);
                                                       (* 9 = tab *)

{
write(output,'***** fff ');
writeafeature(output,groupTotalFeature);
writeln(output);
}
         (* thread the features into the total *)
         anamepart := nameparts; (* anamepart steps through the nameparts *)
         for m := 1 to groupTotalFeature^.definition^.marks do begin
{
write(output,'***** fff mark ',m:1,' is ',location:1,' "');
writestring(output,anamepart^);
writeln(output,'"');
}
{
showafeature(output,groupTotalFeature);
}

            (* Find the feature(s) and the gap(s) for this mark *)
            gg := groupfeature; (* start at the first feature *)
{write(output,'===>TTT loop');}
            while gg <> nil do begin
{TTT}
{write(output,'===>'); writeafeature(output,gg); writeln(output);}

(*
if (gg^.number = 8)
or (gg^.number = 2) then begin
{it's failing to thread #8 for some reason - why? SOLVED}
writeln(output, 'TTT BUG POINT *******8888888888888888888888888888888 BEGIN #',
gg^.number:1);
write(output,'anamepart^:"');
writestring(output, anamepart^);writeln(output,'"');
write(output,'gg^.nametag:"');
writestring(output, gg^.nametag);writeln(output,'"');
if equalstring(anamepart^,gg^.nametag)
then writeln(output,    'equalstring(anamepart^,gg^.nametag)')
else writeln(output,'NOT equalstring(anamepart^,gg^.nametag)');
writeln(output, 'TTT BUG POINT *******8888888888888888888888888888888 END');
end;
*)

{
write(output,'   >'); writeadefinition(output,gg^.definition);
}
               if equalstring(anamepart^,gg^.nametag) then begin
                  if    gg^.definition^.isgap
                     or gg^.definition^.istotal then begin
                     (* If this happens something is VERY wrong! *)
                     writeln(output,'*** PROGRAM SHOULD NEVER BE HERE!');
                     writeln(output,'*** Total found matching gap or total');
                     halt;
                  end

                  else begin
{
writeln(output,'*** THREAD feature #',gg^.number:1,' INTO TOTAL TTT');
}
                     (* thread feature into total *)
                     if groupTotalCarrier = nil then begin
{
writeln(output,'*** if groupTotalCarrier = nil then begin');
}
                        new(groupTotalFeature^.others);
                        groupTotalCarrier := groupTotalFeature^.others;
                        groupTotalCarrier^.next := nil;
                     end
                     else begin
{
writeln(output,'*** if groupTotalCarrier <> nil then begin');
}
                        new(groupTotalCarrier^.next);
                        groupTotalCarrier := groupTotalCarrier^.next;
                        groupTotalCarrier^.next := nil;
                     end;
                     groupTotalCarrier^.afeature := gg;

                     (* thread total into feature! *)
                     if gg^.mothers = nil then begin
                        new(gg^.mothers);
                        gg^.mothers^.next := nil; (* NEW 2016aug30 *)
                        gg^.mothers^.afeature := groupTotalFeature;
{writeln(output,'*** threaded total into feature!');fff}
{IIIIIII}
                     end;
(*writeln(output,'ABOUT TO countmarks at mark number ',m:1);{TTT}*)
{fff}                countmarks(groupTotalFeature,m);
                  end;
               end;
               gg := gg^.next;
            end;

            anamepart := anamepart^.next;
         end; (* for threading marks into total loop *)
         if anamepart <> nil then begin
            writeln(output,'total names not all used');
            showafeature(output,f);
            showadefinition(output,f^.definition);
            halt;
         end;
{fff}    countmarks(groupTotalFeature,-1); (* final check *)

{iii}
{
writeln(output,'========= groupTotalFeature:');
showafeature(output,groupTotalFeature);
writeln(output,'============================');
}

         (* Now use the threaded total feature to thread the gaps!
            Thread all of the gaps in this group.
            Stop when we reach the total *)
         gg := groupfeature;
         while (not gg^.definition^.istotal) and (gg <> nil) do begin
            if gg^.definition^.isgap then begin
{
writeln(output);
writeln(output,'=== thread this gap: '); writeafeature(output,gg); writeln(output);
showafeature(output,gg);
}
{HIYA}
if superdebug then begin
writeln(output,'BUBBA gg^.number=', gg^.featurenumber:1,' ++++++++++++++++++++');
writeln(output,'BUBBA gg is: '); writeafeature(output,gg); writeln(output);
end;
               (* get the names and gap distance from the gap *)
               splitstring(gg^.nametag, '-', chr(9), nameparts); (* 9 = tab *)
               (* Example for a gap with nametag: "p35-(22)-p10 4164377 Gap"
                  nameparts^: "p35"
                  nameparts^.next^: "(22)"
                  nameparts^.next^.next^: "p10"
                *)

{
write(output,'* nameparts^: "');
writestring(output,nameparts^);
writeln(output,'"');
write(output,'* nameparts^.next^: "');
writestring(output,nameparts^.next^);
writeln(output,'"');
write(output,'* nameparts^.next^.next^: "');
writestring(output,nameparts^.next^.next^);
writeln(output,'"');
}

               leftstring  := nameparts^;
               gapstring   := nameparts^.next^;
                  (* gapstring is only needed for checking! *)
               rightstring := nameparts^.next^.next^;

               (* read the gap number: convert gapstring to a number *)
               gap := 0; (* start the number up *)
               if     (gapstring.letters[1] <> '(')
                  and (gapstring.letters[gapstring.length] <> ')')
               then begin
                  writeln(output,'Gap string found without surrounding "()"');
                  showafeature(output,gg);
                  halt;
               end
               else begin
                  gapstring.current := 2;
                  while gapstring.current <= gapstring.length - 1 do begin
                     aletter := gapstring.letters[gapstring.current];
                     if not (aletter in
                             ['0','1','2','3','4','5','6','7','8','9'])
                     then begin
                        writeln(output,'Gap contains non-digit: "',
                                        aletter,'"');
                        showafeature(output,gg);
                        halt;
                     end;
                           gap := (10*gap) + (ord(aletter) - ord('0'));
                     gapstring.current := succ(gapstring.current);
                  end;
               end;
               {if not((gap = abs(gg^.fromrange)) or (gap = abs(gg^.torange)))}
if superdebug then begin
writeln(output,'BUBBA gap is ',gap:1);
writestring(output, leftstring);writeln(output,'<-  leftstring BUBBA');
writestring(output,rightstring);writeln(output,'<- rightstring BUBBA');
end;
               if gap <> gg^.torange - gg^.fromrange
               then begin
                  writeln(output,'Gap in nametag does not match range:');
                  showafeature(output,gg);
                  halt;
               end;
               aGapFeature := gg; (* give the gap a good name *)

               (* Search the group features for ones that match this gap *)
               (* Use the names and gap distance to find the features *)
               leftgap := nil;
               rightgap := nil;
               ff := groupfeature;
               while     (ff <> groupTotalFeature)
                     and ((leftgap=nil) or (rightgap=nil))
               do begin
                  if not (ff^.definition^.isgap or ff^.definition^.istotal)
                  then begin
                     if round(ff^.relcoo) = aGapFeature^.fromrange then begin
{writeln(output,'fromrange coomatch');}
                        if equalstring(leftstring,ff^.nametag) then begin
{writeln(output,'equalstring(leftstring,ff^.nametag)');}
                           leftgap := ff;
                        end;
                     end;
                     if round(ff^.relcoo) = aGapFeature^.torange then begin
{writeln(output,'torange coomatch');}
                        if equalstring(rightstring,ff^.nametag) then begin
{writeln(output,'equalstring(rightstring,ff^.nametag)');}
{showafeature(output,ff);}
                           rightgap := ff;
                        end;
                     end;
                  end;
{writeln(output,'BUBBBA BEFORE ff := ff^.next; ');}
                  ff := ff^.next;
{writeln(output,'BUBBBA AFTER  ff := ff^.next; ');}
               end;

if superdebug then begin
write(output,'BUBBA leftgap: '); writeafeature(output,leftgap); writeln(output);
write(output,'BUBBA rightgap: '); writeafeature(output,rightgap); writeln(output);
end;

               (* Tests of the linked list so far: *)
               if ((leftgap=nil) or (rightgap=nil)) then begin
                  writeln(output,'gap ends not found:');
                  showafeature(output,aGapFeature);
                  halt;
               end;
               if leftgap^.mothers = nil then begin
                  writeln(output,'ERROR: missing total for leftgap');
                  halt;
               end;
               if rightgap^.mothers = nil then begin
                  writeln(output,'ERROR: missing total for rightgap');
                  halt;
               end;

               (* At this point we know the features that match the gap
                  at ff, they are leftgap and rightgap. *)

               (* Next, we thread the two gap features into the gap.
                  There are ONLY two. *)
               (* do the left side of the gap: *)
               if aGapFeature^.others <> nil then begin
                  writeln(output,'ERROR: aGapFeature not nil!');
                  halt;
               end;
               new(aGapFeature^.others);
               groupGapCarrier := aGapFeature^.others;
               groupGapCarrier^.next := nil;
               groupGapCarrier^.afeature := leftgap;
               (* do the right side of the gap: *)
               new(groupGapCarrier^.next);
               groupGapCarrier := groupGapCarrier^.next;
               groupGapCarrier^.next := nil;
               groupGapCarrier^.afeature := rightgap;

               (* Now thread the gap as a mother of each feature.
                  There can be multiple mothers for a feature. *)
{iii}
               addmotherlinks(aGapFeature, rightgap);
               addmotherlinks(aGapFeature,  leftgap);
if superdebug then begin
writeln(output,'BUBBA final aGapFeature is: ========');
showafeature(output,aGapFeature);
writeafeature(output,aGapFeature);writeln(output);
{
writeln(output,'mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm done');
writeln(output,'**********======= rightgap:');
showafeature(output,rightgap);
writeln(output,'**********======= leftgap:');
showafeature(output,leftgap);
}
end;

            end;
            gg := gg^.next

(* initial spot for if gg^.duplicateprevious <> nil then begin
addmotherlinks(aGapFeature, gg^.duplicateprevious);
end; - has pretty good effect. *)

         end; (* end of thread all of the gaps in this group *)

if superdebug then begin
writeln(output,'BUBBA: show AFTER threading:');
gg := groupfeature;
while gg <> nil do begin
writeafeature(output,gg);writeln(output);
gg := gg^.next
end;
end;

         (* Add mother links to the feature that has the vertical bar *)
{ddd}
{
writeln(output,'dddddddddddddddddddddddddddddddddddddd');
}
         gg := groupfeature; (* start at the first feature *)
         while gg <> nil do begin
{
writeln(output,'BUBBA before vertical thing gg is: ==='); writeafeature(output,gg);
}
{
            if gg^.duplicateprevious <> nil then begin
}
            (* 2016oct30 AVOID MODIFYING GAPS!!! *)
            if      (gg^.duplicateprevious <> nil)
                and (not gg^.definition^.isgap)
            then begin
               (* thread the gap into the feature: *)
               addmotherlinks(aGapFeature, gg^.duplicateprevious);
{
writeln(output,'gg^.duplicateprevious AFTER:');
showafeature(output,gg^.duplicateprevious);
writeln(output,'aGapFeature AFTER===================:');
showafeature(output,aGapFeature);
}
               (* now alter aGapFeature so that it points to the
                  non-duplicate feature! *)
{
writeln(output,'what aGapFeature points to:');
writeafeature(output,aGapFeature^.others^.afeature);writeln(output);
}
               (* 2016oct30: KEY CHANGE to fix hindi problem! *)
               if (not aGapFeature^.definition^.isgap)
               then aGapFeature^.others^.afeature := gg^.duplicateprevious;
{
writeln(output,'aGapFeature AFTER FIXED? ===================:');
writeln(output,'what aGapFeature points to AFTER FIX::');
writeafeature(output,aGapFeature^.others^.afeature);writeln(output);
writeln(output,'aGapFeature FINAL RESULT ===================:');
writeln(output,'what aGapFeature points to AFTER FIX::');
showafeature(output,aGapFeature);
}
            end;
{
writeln(output,'BUBBA after vertical thing gg is: ===='); writeafeature(output,gg);
}
            gg := gg^.next;
         end;

if superdebug then begin
writeln(output,'BUBBA: show AFTER motherlinks:========');
gg := groupfeature;
while gg <> nil do begin
writeafeature(output,gg);writeln(output);
gg := gg^.next
end;
end;

         groupTotalFeature := nil; (* reset *)
         ingroup := false;
{fff}
         groupfeature := nil; (* clear the group feature now that the
                                 total parts are threaded! *)
      end; (* total threading *)

{writeln(output,'END OF GROUP ',groupnumber:1,' GGG');}
   end; (* ingroup test *)
end; (* end threadgroup *)

begin (* readdefinitionsandfeatures *)
   reset(afile);
   emptyafile := eof(afile); (* defines if file contains anything *)
{zzzccc}
   comments := nil;
   definitionlist := nil;
   new(featurelist);
   clearfeature(featurelist);
   countdefinitions := 0;
   uniquedefinitions := 0;
   countfeatures := 0;
   uniquefeatures := 0;
   countribls := 0;
   countboundaries := 0;
   d := definitionlist;
   f := featurelist;
   countpetals := 0;
   groupnumber := 0;
   firstfeatureofgroup := false;

{gggg}
   groupfeature := nil;
   groupnumberinfile := 0;
   ingroup := false;
   new(nameparts);
   clearstring(nameparts^); (* start out clean *)
   firstfeature := true;
   while not eof(afile) do begin
      if eoln(afile) (* skip blank line *)
      then readln(afile)
      else begin
         skipblanks(afile);

         (* determine the identifier but don't read it: *)
         identifier := afile^;
{writeln(output,'identifier="',identifier,'"');}

         if identifier = '*'
         then begin (* read comment *)
            if (comments = nil) then begin
               new(comments);
               clearstring(comments^); (* start out clean *)
               currentcomment := comments;
            end
            else begin
               new(currentcomment^.next);
               currentcomment := currentcomment^.next;
            end;
            readstringline(afile,currentcomment^);
{writestringlines(output,currentcomment);}
{zzzccc}
         end

         else if identifier = 'd'
         then begin (* read definition *)
            if definitionlist = nil
            then begin
               new(definitionlist);
               cleardefinition(definitionlist);
               d := definitionlist;
               dprevious := d;
            end
            else begin
               dprevious := d;
               new(d^.next);
               d := d^.next;
               cleardefinition(d);
            end;
            readadefinition(afile, d);
{
write(output, 'read: '); writeadefinition(output,d); writeln(output);
}
            readln(afile); (* this allows comments at the end of the def *)
            countdefinitions := succ(countdefinitions);
            uniquedefinitions := succ(uniquedefinitions);
            d^.definitionnumber := uniquedefinitions;

            (* see if we can locate the definition, based on the nametag *)
            dtest := definitionlist;
            done := false;
            while not done do begin
               if dtest = nil
               then done := true
               else begin
{ add this 2008 Aug 06:}
                  if equaldefinition(dtest, d)
{ remove this 2008 Jul 22:
                  if equalstring(dtest^.nametag, d^.nametag)
}

{
                     and EqualDefLocations(dtest, d)
}
{rrr}
                     and (dtest^.definitionnumber <> d^.definitionnumber)
                  then begin
                     done := true;
                     (* We can't remove the duplicate because features
                        can already point to it!  So just warn them. *)
                     barout;
                     writeln(output, '*** WARNING: A duplicate definition',
                                     ' was found at definition ',
                                   countdefinitions:1,':');
                     write(output, '*** ');
                     writequotestring(output,d^.nametag);
                     writeln(output,'.');
                     writeln(output, '*** ',
                                     'The first definition (number ',
                                     dtest^.definitionnumber:1,
                                     ') will be used.');
{rrr}
{
                     if not EqualDefLocations(dtest, d) then begin
                        writeln(output, 'BUT WAIT!');
                        writeln(output, 'THE LOCATIONS ARE NOT EQUAL!');
                        writeln(output, 'THIS IS FATAL!');
writeln(output,'===');
                        showadefinition(output,dtest);
writeln(output,'===');
                        showadefinition(output,d);
writeln(output,'===');
                        writeadefinition(output,dtest);
writeln(output,'===');
                        writeadefinition(output,d);
writeln(output,'===');
                        barout;

  2006 Oct 14.  Can't be fatal: bars can have different sizes.
                        halt
                     end;
}
                     barout;

{
writeln(output,'definitionlist is:-------- BEFORE 0');
showdefinitions(output,definitionlist);
writeln(output,'-------------------------- BEFORE 1');

if dprevious^.next = d
then begin
writeln(output,'dprevious^.next = d as expected!');
end;

if d^.next = nil
then begin
writeln(output,'d^.next = nil as expected!');
end;
}

                     (* The following solved the duplicate definition bug *)
                     (* remove the item we just read from the list *)
                     dprevious^.next := nil;
                     dispose(d);
                     done := true; (* there is no more to do since there
                                      cannot be other duplicates *)
                     (* move d to the end of the list: *)
                     d := dprevious;
{
this is unnecessary, but worked:
                     d := definitionlist;
                     while d^.next <> nil do d := d^.next;
}

{
writeln(output,'definitionlist is:-------- AFTER 0');
showdefinitions(output,definitionlist);
writeln(output,'-------------------------- AFTER 1');
}

{
writeln(output, 'the dprevious defintion is:');
writeadefinition(output, dprevious); writeln(output);
}

                     uniquedefinitions := pred(uniquedefinitions);
                  end;
                  if not done then begin (* step forward *)
                     dtest := dtest^.next;
                  end;
               end;
            end;
(*
if ingroup then begin
   writeadefinition(output,d); {writeln(output); done by the writeadefinition}
end;
*)
         end

         else if identifier = 'r'
         then begin (* read ribl matrix *)
            (* read in the matrix *)
            getriblmatrix(afile, matrix);

            locatedefinition(matrix^.riblname,definitionlist,thedefinition);
            if thedefinition = nil then begin
               write(output,'A definition for Ribl matrix ');
               writequotestring(output,matrix^.riblname);
               writeln(output,' was not found in the features.');
               halt
            end;

            (* make sure that the definition is in the
               range of the weight matrix. *)
            if (thedefinition^.min < matrix^.frombase) or
               (thedefinition^.max > matrix^.tobase) then begin
               writeln(output, 'Definition:');
               writeadefinition(output,thedefinition);
               writeln(output);
               writeln(output, 'Violates The From-To Bounds',
                       ' of its Matrix:');
               writeln(output,'matrix frombase: ',matrix^.frombase:1);
               writeln(output,'matrix   tobase: ',matrix^.  tobase:1);
               halt
            end;

            (* slip the matrix into the data structure *)
            thedefinition^.matrix := matrix;
            countribls := succ(countribls);
         end

         else if identifier = 'b'
         then begin (* define boundaries *)
            (* read the boundary name *)
            skipnonblanks(afile);
            skipblanks(afile);
            readquotestring(afile, aname);
            locatedefinition(aname,definitionlist,thedefinition);
            if thedefinition = nil then begin
               write   (output,'WARNING: A definition for the boundaries: ');
               writequotestring(output,aname);
               writeln(output,'          was not found in the features.');
            end
            else with thedefinition^ do begin
               readln(afile, RiboundFeature, ZboundFeature, PboundFeature);
               countboundaries := succ(countboundaries);
            end;
         end

         else if identifier = 'p'
         then begin (* define color rectangle, a petal *)
         (* petal "ribl-name" el c eh es eb fh fs fb th *)

            (* read the color rectangle name *)
            skipnonblanks(afile);
            skipblanks(afile);
            readquotestring(afile, aname);
            locatedefinition(aname,definitionlist,thedefinition);
            if thedefinition = nil then begin
               write  (output,'WARNING: A definition',
                                      ' for the color rectangle: ');
               writequotestring(output,aname);
               writeln(output);
               writeln(output,'         was not found in the features.');
               readln(afile);
            end
            else with thedefinition^.petal do begin
               read(afile, el);
{
writeln(output,'the value of el is ',el:1);
}
               skipblanks(afile);
               read(afile, c);
               if not ((c='r') or (c='h')) then begin
                  yell;
                  writeln(output,' feature');
                  writequotestring(output,aname);
                  writeln(output,': color type must be r or h, not "', c,'"');
                  halt;
               end;
               readln(afile, eh, es, eb, fh, fs, fb, spare);

{
               if
                  (eh < 0) or (eh > 1) or
                  (es < 0) or (es > 1) or
                  (eb < 0) or (eb > 1) or
                  (fh < 0) or (fh > 1) or
                  (fs < 0) or (fs > 1) or
                  (fb < 0) or (fb > 1)
               then begin
                  yell;
                  writeln(output,' feature');
                  writequotestring(output,aname);
                  writeln(output,': eh, es, eb, fh, fs, and fb ');
                  writeln(output,'must have values between 0 and 1');
                  halt;
               end;
}

               if
                  (eh > 1) or
                  (es > 1) or
                  (eb > 1) or
                  (fh > 1) or
                  (fs > 1) or
                  (fb > 1)
               then begin
                  yell;
                  writeln(output,' feature');
                  writequotestring(output,aname);
                  writeln(output,': eh, es, eb, fh, fs, and fb ');
                  writeln(output,'must have values less than 1');
                  halt;
               end;

{
This is a spare variable for now so it has no restrictions
               if
                  (th < 0) or (th > 1)
               then begin
                  yell;
                  writeln(output,' feature');
                  writequotestring(output,aname);
                  writeln(output,': spare');
                  writeln(output,'must have values between 0 and 1');
                  halt;
               end;
}

               countpetals := succ(countpetals);
            end;
         end

         else if identifier = '@'
         then begin (* read feature *)
            countfeatures := succ(countfeatures);
{
            f^.featurenumber := countfeatures;
            f^.featuregroup := groupnumber;
            f^.featureisfirstingroup := firstfeature; (* mark into the feature *)
}

            if firstfeature then begin (* first feature of entire list *)
               (* the feature must be unique, just read it in *)
               readafeature(afile, definitionlist, f);
               readln(afile); (* this allows comments
                                 at the end of the feature *)
               f^.featurenumber := countfeatures;
               f^.featuregroup := groupnumber;
               f^.featureisfirstingroup := true;
               uniquefeatures := succ(uniquefeatures);
               firstfeature := false; (* set to false for later features *)
            end
            else begin
               (* determine if the feature is unique *)
               new(newfeature);
               clearfeature(newfeature);
               readafeature(afile, definitionlist, newfeature);
               readln(afile); (* this allows comments
                                 at the end of the feature *)
               if not featureinlist(featurelist, newfeature, oldfeature)
               then begin
                  (* since the feature is not in the list, add it in *)
                  f^.next := newfeature;
                  f := f^.next;
                  f^.featurenumber := countfeatures;
                  f^.featuregroup := groupnumber;
                  f^.featureisfirstingroup := firstfeatureofgroup;
{zzzfff}
                  uniquefeatures := succ(uniquefeatures);
               end
               else begin
                  writeln(output,'This feature was found more than once:');
                  writeafeature(output,newfeature);
                  writeln(output);

{2007 Mar 26 ppp}

                  (* Keep all duplicates in the list! *)
                  f^.next := newfeature;
                  f := f^.next;
                  f^.featurenumber := countfeatures;
                  f^.featuregroup := groupnumber;
{zzzfff}
                  f^.featureisfirstingroup := firstfeatureofgroup;

                  if KeepDuplicates
{ppp do something about orfs here }
                  then begin
                     writeln(output,'It is being kept.');
                     f^.duplicateprevious := nil;
                  end
                  else begin
                     writeln(output,'It will be ignored.');
                     f^.duplicateprevious := oldfeature;
                  end

               end
            end;
{zzzfff}
            if firstfeatureofgroup then firstfeatureofgroup := false;

            threadgroup;
{ttt Is this the right place to examine the list or should it be later?}
{
rewrite(list);
writeln(list,'after threadgroup showfeatures(list,featurelist, true);');
showfeatures(list,featurelist, true);
}
         end

         else if identifier = 'g'
         then begin (* start of group feature *)
            ingroup := true;
            groupnumber := groupnumber + 1;
{writeln(output, '--- groupnumber is now ',groupnumber:1);}
            mycharacter := '?';
{
writeln(output,'mycharacter is "',mycharacter,'"');
writeln(output,'ord(mycharacter) is "',ord(mycharacter):1,'"');
}
            while (mycharacter <> ' ') and (not eoln(afile)) do begin
               (* read past the word 'group' *)
               read(afile,mycharacter);
{writeln(output, '--- mycharacter "',mycharacter,'"');}
               if eoln(afile) then begin
                  writeln(output,'found end of line while reading "group"');
                  halt;
               end;
            end;
            readln(afile, groupnumberinfile);
{writeln(output, 'groupnumberinfile=', groupnumberinfile:1);}
{zzqq}
            if groupnumberinfile <> groupnumber then begin
               writeln(output, 'WARNING: counted groupnumber=',groupnumber:1,
                             ' does not match groupnumberinfile=',
                             groupnumberinfile:1);
{showfeatures(output,featurelist, true);}
            end;
            firstfeatureofgroup := true;
         end

         else begin
            writeln(output,'features file lines must be blank or',
                           ' begin with one of "*d@g', (* g=group *)
                           'rbp',
                           '", not "',identifier,'".');
            halt
         end
      end
   end;

{
showfeatures(output,featurelist, true);
writeln(output,'PASSPOINT');
halt;
}
{zzzqqqeee}

   (* 2008 Mar 11. There can be a non-empty features file
      which only contains definitions and no features.
      Call this an empty features file. *)
   if featurelist <> nil then begin
      if featurelist^.id.length = 0 then begin
         emptyafile := true;
      end
   end;

   if not emptyafile then begin
      (* check that any definition which claims to have a walker
         actually has a defined weight matrix *)
      if definitionlist <> nil then begin
         d := definitionlist;
         missingmatrix := false;
         while d <> nil do begin
{
write(output, 'testing:');
writeadefinition(output,d);
writeln(output);
}
            with d^ do if background.length = 2
            then if background.letters[2] in ['w','W','l','L','i','I']
            then if d^.matrix = nil then begin
               missingmatrix := true;
               writeln(output, 'Definition:');
               writeadefinition(output,d);
               writeln(output);
               writeln(output,'is missing a Ribl weight matrix');
            end;
            d := d^.next;
         end;
         if missingmatrix then begin
            writeln(output,'readdefinitionsandfeatures:',
                           ' missingmatrix is true');
            halt
         end;
      end;
   end; (* file check *)

   writeln(output,'**************************************************');
   if emptyafile then begin
      writeln(output,'** Empty Features File ***************************');
      writeln(output,'** The features file may have characters in it, **');
      writeln(output,'** but it contains no features.                 **');
      writeln(output,'**************************************************');
      (* make sure we don't give any empty features *)
      dispose(featurelist);
      featurelist := nil;
   end
   else begin
      write(output,countdefinitions:1,' definition');
      plural(output,countdefinitions);
      writeln(output,' read in');

      write(output,uniquedefinitions:1,' unique definition');
      plural(output,uniquedefinitions);
      writeln(output,' found');

      write(output,countribls:1,' definition');
      if countribls = 1
      then writeln(output,' has an Ribl matrix')
      else writeln(output,'s have Ribl matricies');

      write(output,countboundaries:1,' definition');
      if countboundaries = 1
      then writeln(output,' boundry defined')
      else writeln(output,' boundries defined');

      write(output,countpetals:1);
      if countpetals = 1
      then writeln(output,' petal (color rectangle behind a walker) defined')
      else writeln(output,' petals (color rectangles behind walkers) defined');

      write(output,countfeatures:1,' feature');
      plural(output,countfeatures);
      writeln(output,' read in');

      write(output,uniquefeatures:1,' unique feature');
      plural(output,uniquefeatures);
      writeln(output,' found');

{
(* show just the first feature *)
writeln(output,'~~~~~~~~~~~~~~~ p10 ~~~~~~~~~~~~~~~');
showafeature(output,featurelist);
writeln(output,'the feature should have 3 mothers but it is missing the p35!');
2016sep03fixed!
writeln(output,'this feature now has the correct 4 mothers!');
}
{halt;}
{iii debug}

if false then begin (* final display for debugging *)
writeln(output);
writeln(output,'******* @ end of readdefinitionsandfeatures: *******');
writeln(output,'******* DEFINITIONS:');
if false then showdefinitions(output,definitionlist);
writeln(output);
writeln(output,'******* FEATURES:');
if true then showfeatures(output,featurelist,true);
writeln(output);
writeln(output,'******* @ end of readdefinitionsandfeatures: *******'); {sss}
{
halt;
}
end;
{
writeln(output, 'end readdefinitionsandfeatures halt'); halt;
}

   end;
   writeln(output,'**************************************************');
end;
(* end module readdefinitionsandfeatures version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module isnamestring *)
function isnamestring(n: name; s: string): boolean;
(* is the delila-type name n the same as the string s? *)
var
   c: integer; (* index to n and s *)
   done: boolean; (* done checking *)
begin
(*
write(output,'isnamestring s(',s.length:1,'):');
writestring(output,s);
write(output,', n(',n.length:1,'):');
for c := 1 to n.length do write(output,n.letters[c]);
writeln(output);
*)
   if n.length = s.length then begin
      c := 1;
      done := false;
      while not done do begin
(*
writeln(output,'   ',s.letters[c],' ' ,n.letters[c]);
*)
         if s.letters[c] = n.letters[c]
         then begin
           c := succ(c);
           if c > s.length then begin
              isnamestring := true;
              done := true
           end
         end
         else begin
            isnamestring := false;
            done := true
         end
      end
   end
   else isnamestring := false;
end;
(* end module isnamestring version = 13.24; (@ of lister.p 2016 Oct 08 *)

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

(* 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.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module writename *)
procedure writename(var afile: text; s: name);
(* write the string s to afile *)
var i: integer; (* index to the string s *)
begin
   with s do for i := 1 to length do write(afile,letters[i]);
end;
(* end module writename version = 9.25; (@ of lister.p 1999 June 9 *)

(* begin module numberdigit *)
function numberdigit(number, logplace:integer): char;
(* return the digit at the place value ('logplace') position of number.
example:
      numberdigit(13625, 3) = 3
      numberdigit(13625, 4) = 1
2000 July 30 'myabsolute' replaced 'absolute', which
is apparently a keyword for GPC.  The name is kept
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 *)
procedure sign;
(* put a negative sign out or a positive sign *)
begin (* sign *)
      if number <0 then acharacter:='-'
                   else acharacter:='+'
end; (* sign *)
begin (* numberdigit *)
      place:=1;
      for count:=1 to logplace do place:=10*place;

      if number=0 then begin
         if place=1 then acharacter:='0'
                    else acharacter:=' '
      end
      else begin
         myabsolute:=abs(number);
         if myabsolute < (place div 10)
            then acharacter:=' '
            else if myabsolute >= place
                    then digit
                    else sign
      end;
      numberdigit:=acharacter
end; (* numberdigit *)
(* end module numberdigit version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module numbersize *)
function numbersize(n: integer):integer;
(* calculate amount of space to be reserved for the integer n *)
const ln10 = 2.30259; (* natural log of 10 - for conversion to log base 10 *)
      epsilon = 0.00001; (* a small number to correct log base 10 errors *)
var
   size: integer; (* intermediate result *)
begin (* numbersize *)
   if n = 0
   then numbersize:=1
   else begin
      size:=trunc(ln(abs(n))/ln10 + epsilon) + 1;
      (* the 1 is for the last digit *)
      (* the epsilon assures that we do not lose a place due to
         roundoff.  eg, sometimes log base 10 of 10 would be 0.9999
         instead of 1, and we would not do it right...  note: this will
         fail for very large numbers on the order of 1/epsilon. *)
      if n < 0 then size := succ(size);  (* account for minus sign *)
      numbersize := size;
   end
end; (* numbersize *)
(* end module numbersize version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module numberstring *)
procedure numberstring(anumber: integer; var astring: string);
(* Convert an integer number into a string.  Positive numbers do
not have a sign.  REQUIRES: numbersize and numberdigit. *)
var
   d: integer; (* the digit extracted *)
   size: integer; (* size of the number in digits *)
begin
   clearstring(astring);
   size := numbersize(anumber) - 1;
   for d := 0 to size do begin
      astring.letters[1+size-d] := numberdigit(anumber, d);
   end;
   astring.length := size+1;
end;
(* end module numberstring version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module concatstring *)
procedure concatstring(var a, b: string; var c: string);
(* concatenate a and b and put in c *)
var   abindex: integer; (* index to both strings *)
      cindex: integer; (* index to both strings *)
begin (* concatstring *)
   emptystring(c);
   c.length := a.length + b.length;
   if c.length > maxstring then begin
      writeln(output,'string too long in concatstring');
      halt
   end;

   cindex := 0;
   for abindex := 1 to a.length do begin
      cindex := succ(cindex);
      c.letters[abindex] := a.letters[abindex]
   end;

   for abindex := a.length+1 to a.length+b.length do begin
      cindex := succ(cindex);
      c.letters[abindex] := b.letters[abindex-a.length]
   end;
end; (* concatstring *)
(* end module concatstring version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

(* begin module attachstring *)
procedure attachstring(var a: string; b: string);
(* attach the string b onto the end of a *)
var
   hold: string; (* intermediate for doing the attachment *)
begin
   copystring(a, hold);
   concatstring(hold,b,a)
end;
(* end module attachstring version = 5.85; (@ of prgmod.p 2022 Dec 08 *)

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

(* begin module orf.themain *)
procedure themain(var genomebook, scanfeatures, orfp, orffeatures: text);
(* the main procedure of the program *)
var
   parameterversion: real; (* parameter version number *)
   debug: boolean; (* set true to debug *)

   apiece: pieceptr; (* a piece *)
   theline: integer; (* the current line number in the genomebook *)
   sequences: integer; (* number of the sequences in a genomebook *)
   orfs: integer; (* number of orfs reported *)
   definitionlist: definitiontypeptr; (* all definitions
                                         read in from scanfeatures *)
   featurelist: featuretypeptr; (* all features read in from scanfeatures *)
   KeepDuplicates: boolean; (* If set to true, duplicate features
      are kept, if false, they are eliminated.
      argument to readdefinitionsandfeatures *)
   comments: stringptr; (* argument to readdefinitionsandfeatures *)

   ir: string; (* the name of the initiation region *)
   f: featuretypeptr; (* a feature in the feature list *)
   b: base; (* a base found in the piece apiece *)
   b0,b1,b2: base; (* bases of a codon *)
   l: integer; (* position in a site *)
   pe: integer; (* position in the piece (external coordinates) *)
   pi: integer; (* position in the piece (internal coordinates) *)

   codon: integer; (* the codon we are at, initiation codon is zero. *)
   done: boolean; (* done reading the frame (or at end of sequence) *)
   stopcodon: boolean; (* true if we found a stop codon *)
   pielength: integer; (* length of the piece apiece *)
   totalfeature: featuretypeptr; (* the total feature for the ir feature *)

   d: definitiontypeptr; (* a definition for copying scanfeatures to
                            orffeatures *)
   matrix: ribltypeptr; (* a matrix read from scanfeatures for copying
                            to orffeatures *)

   orf: featuretypeptr; (* an orf feature constructed in this program *)
   orfdef: definitiontypeptr; (* orf definition *)
   fhold: featuretypeptr; (* a pointer for putting the orf into the
     featurelist *)
   dhold: definitiontypeptr; (* a pointer for putting the orfdef into the
     definitionlist *)

   codonword: string; (* the word ' condons' for use in displaying
                         the codon features *)
   orfname: string; (* the word 'orf' for use in displaying *)

   shortestorf: integer; (* shortest orf for total feature *)
   longestorf:  integer; (* shortest orf for total feature *)

   pieorientation: integer; (* orientation of the piece as an integer *)

newfeature: featuretypeptr; (* BUBBA *)

{aaa}

procedure report(var afile: text);
(* report results to a file *)
begin
   writeln(afile,'*');
   writeln(afile,'* ',sequences:1, ' sequences are in the genomebook');
   writeln(afile,'* ',orfs:1, ' ORFs are in the genomebook');
end;

function grabbase(l: integer): base;
(* grab base at position l relative to the zero coordinate
of the ir *)
begin
  pe := round(f^.coordinate+l*f^.orientation);
  pi := pietoint(pe, apiece);
  b := getbase(pi, apiece);
  grabbase := b;
{
  writeln(output, basetochar(b), ' ', pe:5);
}
end;

begin (* themain *)
   writeln(output,'orf ',version:4:2);
   reset(orfp);
   readln(orfp, parameterversion);
   if round(100*parameterversion) < round(100*updateversion) then begin
      writeln(output, 'You have an old parameter file!');
      halt
   end;

   readln(orfp, shortestorf);
   readln(orfp, longestorf);

   KeepDuplicates := true;
   readdefinitionsandfeatures(scanfeatures,
                       KeepDuplicates,
                       comments,
                       definitionlist,
                       featurelist);

{qqq}
{
if (featurelist = nil) then begin
   writeln(output,'nil list bubba *****************EWEEEEEEEEEEE');
end
else begin
   writeln(output,'non nil list bubba ****&&&&&&&&&&&&&&&&&&&&&&&&EE');
   showfeatures(output,featurelist);
   writeln(output,'non nil list bubba ****&&&&&&&&&&&&&&&&&&&&&&&&EE');
end;
}

   debug := false;
   clearstring(ir);
   with ir do begin
      letters[1] := 'i';
      letters[2] := 'r';
      length := 2
   end;

   (* build definition for orfs *)
   new(orfdef);
   cleardefinition(orfdef);
{yyy}
   with orfdef^ do begin
      with nametag do begin
         letters[1] := 'o';
         letters[2] := 'r';
         letters[3] := 'f';
         length := 3
      end;

      with background do begin
         letters[1] := '#';
         length := 1
      end;

      with negparts do begin
         letters[1] := '<';
         letters[2] := ']';
         length := 2
      end;

      with posparts do begin
         letters[1] := '[';
         letters[2] := '>';
         length := 2
      end;

      locations[1] :=  0.0;
      locations[2] :=  2.0;
      marks := 2;
      min := 0;
      max := 1;

   end;

   emptystring(codonword);
   with codonword do begin
      letters[1] := ' ';
      letters[2] := 'c';
      letters[3] := 'o';
      letters[4] := 'd';
      letters[5] := 'o';
      letters[6] := 'n';
      letters[7] := 's';
      length := 7
   end;

   emptystring(orfname);
   with orfname do begin
      letters[1] := 'o';
      letters[2] := 'r';
      letters[3] := 'f';
      length := 3
   end;

{
The list isn't written out and besides, the location has
to be before the first orf feature
   (* thread the new definition into the definitionlist *)
   dhold := definitionlist;
   while dhold^.next <> nil do begin
      dhold := dhold^.next;
writeln(output, 'definition');
writeadefinition(output,dhold);
   end;
if dhold^.next <> nil then begin
  writeln(output, 'dhold^.next <> nil');
end
else begin
  writeln(output, 'dhold^.next = nil');
end;
   dhold^.next := orfdef;
   dhold := dhold^.next;
writeln(output, 'last definition');
writeadefinition(output,dhold);
halt;
}

   brinit(genomebook, theline);
   rewrite(orffeatures);

   writeln(orffeatures,'* orf ',version:4:2);
   copyaline(genomebook,orffeatures);

   writeadefinition(orffeatures,orfdef);
writeadefinition(output,orfdef);
writeln(output);

   new(apiece);
   sequences := 0;
   orfs := 0;
   writeln(output,'reading genomebook');

   while not eof(genomebook) do begin
      getpiece(genomebook,theline,apiece);
      if not eof(genomebook) then begin
         sequences:=succ(sequences);

         write(output,'piece ');
         writename(output, apiece^.key.hea.keynam);
         writeln(output);

         f := featurelist;
         while f <> nil do begin
            (* is this feature on this piece of DNA? *)
            if isnamestring(apiece^.key.hea.keynam, f^.id) then begin
if debug then writeln(output,'ids MATCH!');

               (* is this feature an ir feature? *)
               if equalstring(f^.nametag, ir) then begin
                  writeln(output);
                  writeafeature(output,f);
                  writeln(output);

if debug then writeln(output,'names MATCH!');
{ they don't need to match orientation - just read in the right
orientation.
   if ((f^.orientation = +1) and
       (apiece^.key.piedir = plus))
      or
      ((f^.orientation = -1) and
       (apiece^.key.piedir = minus))
   then begin
writeln(output,'orientation matches');
halt;
   end;

for l := 0 to 10 do begin
  pe := round(f^.coordinate+l*f^.orientation);
  pi := pietoint(pe, apiece);
  b := getbase(pi, apiece);
  writeln(output, basetochar(b), ' ', pe:5);
end;
}

{
writeln(output,
'f^.orientation = ',
f^.orientation);
}


{
                  write(output,round(f^.coordinate):1,' ');
                  for l := 0 to 11 do begin

                     pe := round(f^.coordinate+l*f^.orientation);
                     pi := pietoint(pe, apiece);
                     b := getbase(pi, apiece);
                     if (l > 0) and (l mod 3 = 0) then write(output,'|');
                     case f^.orientation of
                        +1: begin
                           write(output, basetochar(b));
                        end;
                        -1: begin
                           write(output, basetochar(complement(b)));
                        end;
                     end;
                  end;
                  write(output, '...');
}

                  codon := 0;
                  stopcodon := false;

                  done := false;
                  pielength := piecelength(apiece);
                  case apiece^.key.piedir of
                     plus:  pieorientation := +1;
                     minus: pieorientation := -1;
                  end;
writeln(output);
                  write(output,round(f^.coordinate):1,' ');
                  while not done do begin

                     case f^.orientation*pieorientation of
                        +1: begin
                            b0 := grabbase(codon + 0);
                            b1 := grabbase(codon + 1);
                            b2 := grabbase(codon + 2);
                        end;
                        -1: begin
                            b0 := complement(grabbase(codon + 0));
                            b1 := complement(grabbase(codon + 1));
                            b2 := complement(grabbase(codon + 2));
                        end;
                     end;
                     write(output, basetochar(b0));
                     write(output, basetochar(b1));
                     write(output, basetochar(b2));
                     write(output,'|');
{zzz}

                     (* Identify a stop codon as taa, tga or tag: *)
                     if (b0 = t) and (
                           ((b1 = a) and (b2 = a))
                           or
                           ((b1 = g) and (b2 = a))
                           or
                           ((b1 = a) and (b2 = g))
                        )
                     then begin
                        done := true;
                        stopcodon := true;
                        orfs := orfs + 1;

if debug then begin
writeln(output,'stop codon! ',
basetochar(b0),
basetochar(b1),
basetochar(b2));
writeafeature(output,f);
writeln(output);
writeln(output,'----- this feature');
showafeature(output,f);
writeln(output,'----- the first mother of this feature GAP');
showafeature(output,f^.mothers^.afeature);
writeln(output,'----- the second mother of this feature TOTAL');
showafeature(output,f^.mothers^.next^.afeature);
writeln(output,'-----');
writeafeature(output,totalfeature);
{ writeafeature(output,f^.mothers^.next^.afeature); }
writeln(output);
writeln(output,'----- set Dparam to be codons');
end;

                        (* Sneak the orf info into the total feature:
                        Aparam: information in bits
                        Bparam: length in codons of the ORF
                        Cparam: end of orf, first base of stop codon
                        Dparam: 1 for inside, 0 for outside size of
                                orf range defined by shortestorf and
                                longestorf
                        *)

                        totalfeature := f^.mothers^.next^.afeature;
{
(* Is this really the totalfeature? *)
writeln(output,'BUBBA putative totalfeature: ===================++++++++++++');
writeafeature(output,totalfeature);
writeln(output); (* yup *)
}

                        with totalfeature^ do begin
                           (* Insert data into the total feature! *)
                           (* orf-codons *)
                           Bparam := (codon/3)+1; (* +1 => including stop *)


                           (* orf-end - account for piece orientation!*)
{aaa}
                           if f^.orientation*pieorientation = +1
                           then Cparam := 
                              inttopie(
                                 pietoint(
                                 round(f^.coordinate),apiece)+codon,
                                 apiece)
                           else Cparam := 
                              inttopie(
                                 pietoint(
                                 round(f^.coordinate),apiece)-codon,
                                 apiece);
                           (* For the negative orientation, we need to
                           switch to the OTHER end of the codon. *)

                           (* +codon   is the first base of the stop codon *)
                           (* +codon+2 is the  last base of the stop codon *)

                           (* implement cutoff in the total line *)
                           if (Bparam >= shortestorf) and
                              (Bparam <=  longestorf)
                           then Dparam := 1 (* THIS flags the orf cutoff! *)
                           else Dparam := 0;

                           writeln(output,
                           round(Cparam):1,' - ',
                           round(Bparam):1,' codon orf');
                        end;

                        (* build the orf feature *)

                        (* build the fixed parts of the orf features *)
                        new(orf);
                        clearfeature(orf);
                        copystring(orfname, orf^.nametag);

                        orf^.id := f^.id;
                        orf^.coordinate := totalfeature^.Cparam;

                        (* +codon   is the first base of the stop codon *)
                        (* +codon+2 is the  last base of the stop codon *)
                        (* Can't use totalfeature^.Cparam here since
                           total records the last base; features need
                           the first base. *)

                        orf^.orientation := totalfeature^.orientation;

                        (* put the initial coordinate , length,
                        "codons" here *)

                        numberstring(round(codon/3)+1, orf^.othertag);
                        attachstring(orf^.othertag, codonword);

                        (*  orf color is strong red: *)
                        orf^.Tparam := 'R'; (* use color bar *)
                        orf^.Aparam := 1.0; (* red *)
                        orf^.Bparam := 0.0; (* green *)
                        orf^.Cparam := 0.0; (* blue *)
                        orf^.Dparam := 1.0; (* thickness *)

                        (* orf color is color of ir (failed attempt): *)
{
(* Can we use the ir color for the stop codon color bar?
No!  The color is not here yet! - this gives black because mkpetal
happens LATER - there is no way to get that color now!
...
This can be done only in lister when there are colors available.
Method: search the features for ir.  If the ir is followed by an
orf, assign the color of the ir to the orf. *)
                        orf^.Tparam := 'H'; (* use color bar *)
                        orf^.Aparam := f^.Aparam; (* hue *)
                        orf^.Bparam := f^.Bparam; (* saturation *)
                        orf^.Cparam := f^.Cparam; (* brightness *)
                        orf^.Dparam := 1.0; (* thickness *)
}

{yyy}

                        orf^.definition := orfdef;
                        writeln(output,'orf:');
                        writeafeature(output,orf);
                        writeln(output);
{
                        writeln(output,'---- show the orf ---');
                        showafeature(output,orf);
                        writeln(output,'---------------------');
}
{
if f^.next = nil then begin
writeln(output,'THE END');
end
else begin
writeln(output,'THE MIDDLE');
end;
}
                        (* thread the orf into the features *)
                           fhold := f^.next;
                           f^.next := orf;
                           orf^.next := fhold;
{
writeln(output, 'totalfeature^.Aparam = "', totalfeature^.Aparam:1,'"');
writeln(output, 'totalfeature^.Tparam = "', totalfeature^.Tparam,'"');
}

if debug then begin
writeafeature(output,totalfeature);
writeln(output);
end;

                     end;
                     codon := codon + 3;
{zzz}
{
                     if pietoint(round(f^.coordinate),apiece)+codon
}
                     if pietoint(round(f^.coordinate),apiece)
                       + codon * f^.orientation*pieorientation
                       > pielength
                     then done := true;

                     (* 2013 Nov 04 I'm making this official code since
                        it terminates when the orf end is past the piece *)
                     if pietoint(round(f^.coordinate),apiece)+codon
                        > pielength
                     then begin
                        writeln(output,'- ORF IS NOT TERMINATED WITH A STOP');
                        writeln(output,'pielength ',
                                        pielength:1);
                        writeln(output,'codon ',
                                        codon:1);
                        writeln(output,'round(f^.coordinate) ',
                                        round(f^.coordinate):1);
                        writeln(output,'pietoint(round(f^.coordinate),apiece) ',
                                        pietoint(round(f^.coordinate),apiece):1);
                        writeln(output,'pietoint(round(f^.coordinate),apiece)',
                                       '+codon ',
                                        pietoint(round(f^.coordinate),apiece)
                                        +codon:1);
                        done := true;
                     end;

                  end;

               end;
            end;
            f := f^.next;
         end;
         clearpiece(apiece) (* clear the piece for reuse *)
      end
   end;

   writeln(output,'genomebook reading is done');

   (* Copy the scanfeatures to the orffeatures except
      when a feature is in the featurelist use the featurelist version.
      This substitutes changes to the featurelist (made above),
      namely the altered total features. *)
   f := featurelist;
   new(d);
   reset(scanfeatures);
   while not eof(scanfeatures) do begin
      if eoln(scanfeatures) then begin
         readln(scanfeatures);
         writeln(orffeatures);
      end
      else begin
if debug then begin
writeln(output,'scanfeatures^ = "', scanfeatures^,'"');
end;
         if not (scanfeatures^ in ['*','b','p',' ','g','d','r','@'])
         then begin
            writeln(output,'Illegal character in scanfeatures: "',
                            scanfeatures^,'"');
            writeln(output,'The first character of scanfeature lines');
            writeln(output,' can only contain "*bp dr@g"');
            halt;
         end;
         case scanfeatures^ of
            '*','b','p',' ','g': begin
               (* * = comment, b = boundaries
                  p = petal, g = group *)
               copyaline(scanfeatures, orffeatures);
               readln(scanfeatures);
            end;
            'd': begin
               readadefinition(scanfeatures, d);
               writeadefinition(orffeatures, d);
writeln(output,'BUBBAnow, read a definition:');
writeadefinition(output, d);
               readln(scanfeatures);
            end;
            'r': begin
               getriblmatrix(scanfeatures, matrix);
               putriblmatrix( orffeatures, matrix);
            end;
            '@': begin
writeln(output,'@ location');
               (* Substitute here.  The order has to be the
                  same as we originally read! *)

(* At this point we could
               f := f^.next;
when f is the orf feature!!
zzznnn
if equalstring(orfname, f^.nametag) then begin writeln(output,'HALTING');halt; end;
if equalstring(orfname, f^.nametag) then f := f^.next;
That failed - it never was true!
*)

new(newfeature); (* BUBBA required to avoid bus error *)
clearfeature(newfeature);
(* note: use old definitionlist to keep readafeature happy! *)
               readafeature(scanfeatures, definitionlist, newfeature);
               readln(scanfeatures); (* could do checks here!! *)
write(output,'current feature newfeature: ');
writeafeature(output, newfeature);
writeln(output);
{writeafeature(orffeatures, newfeature); writeln(orffeatures);}
               (* Do the switch to the modified feature.
                  We CANNOT use the newfeature because it does not
                  have the right links! *)
{
if newfeature^.definition^.istotal
then begin
}
               writeafeature(orffeatures, f); writeln(orffeatures);
{
end
else
writeafeature(orffeatures, newfeature); writeln(orffeatures);
}
if newfeature^.definition^.istotal then begin
writeln(output,'****** COMPARE newfeature TO f: *******');
writeafeature(output, newfeature); writeln(output,'=newfeature');
writeafeature(output, f         ); writeln(output,'=f');
writeln(output,'***************************************');
{WORKING HERE}
  writeln(output,'feature newfeature is a total: ');
  writeafeature(output, newfeature);
  writeln(output);
  writeln(output,'HERE is the spot to output the orf, just after the total');
  write(output,'feature f: ');
  writeafeature(output, f);
  writeln(output);
  (* use the current linked list: *)
{
  writeafeature(orffeatures, f);
  writeln(orffeatures);
}
end;
               f := f^.next;

               (* if there is a new orf next, write it out. *)
               if f <> nil then begin
                  if equalstring(orfname, f^.nametag) then begin
                     writeafeature(orffeatures, f);
                     writeln(orffeatures);
                     f := f^.next;
                  end;
               end;
            end;
         end;
      end;
   end;
{
if not eof(scanfeatures) then writeln(output, 'NOT EOF scanfeatures! BUBBA')
                         else writeln(output,  'IS EOF scanfeatures! BUBBA');
}

   (* Put out any additional features that were missed! *)
   if f <> nil then begin
      while f <> nil do begin
         writeln(output);
         writeln(output,'THIS FEATURE WAS NOT USED, writing it out:');
         writeafeature(output, f);
         writeln(output);
         writeafeature(orffeatures, f);
         writeln(orffeatures);
         f := f^.next;
      end;
   end;

   report(orffeatures);
   report(output);

end;
(* end module orf.themain *)

begin
   themain(genomebook, scanfeatures, orfp, orffeatures);
1: end.

