program multibest(scanfeatures, multibestp, bestfeatures, lowerfeatures, output);
(* multibest: find features with the highest Ri in a given sequence range

  Elaine Bucheimer

modified by:
  Dr. Thomas D. Schneider
  toms@alum.mit.edu
  https://alum.mit.edu/www/toms/

*)

label 1; (* end of program *)

const
(* begin module version *)
   version = 1.37; (* of multibest.p 2022 Dec 19
2022 Dec 19: 1.37 check that prgmod modules are ok
2022 Dec 15: 1.36 prepare to go from O(n^2) to O(n)
2022 Dec 09: 1.35 fix output file header line
2022 Dec 09: 1.34 report counts of groups kept and removed
2022 Dec 09: 1.33 clean up
2022 Dec 09: 1.32 remove duplicate groups
2022 Dec 09: 1.31 tighten up speed Tricks
2022 Dec 09: 1.30 check for g <> h ...
2022 Dec 09: 1.29 document new logic, drop out of h loop, clean up old code.
2022 Dec 09: 1.28 clean up
2022 Dec 08: 1.27 extremely simplify filtering!!!
2022 Dec 08: 1.26 still infinite looop
2022 Dec 08: 1.25 solve case of last group removed
2022 Dec 08: 1.24 debug to handle triple sites - it works now!!
2022 Dec 08: 1.23 debug to handle triple sites
2022 Dec 08: 1.22 debug to handle triple sites
2022 Dec 08: 1.21 rebuild logic, solve infinite loops
2022 Dec 07: 1.20 detect infinite loops
2022 Dec 07: 1.19 g does NOT advance in list when deletion!, goes to preg!
2022 Dec 07: 1.18 work on 3 walker case at 583684 in E. coli
2022 Dec 07: 1.17 more tests
2022 Dec 07: 1.16 more tests - functional!!!!
2022 Dec 07: 1.15 it is possible ... infinite loops solved
2022 Dec 07: 1.14 debug g=preg = THE impossible PROGRAM ERROR
2022 Dec 07: 1.13 clean up
2022 Dec 07: 1.12 clean up
2022 Dec 07: 1.11 debug
2022 Dec 07: 1.10 rename best -> multibest
2022 Dec 07: 1.10 rename best -> multibest
2022 Dec 06: 1.09 handle comments in features
2022 Dec 06: 1.08 group filtering - bestfeatures is correct!
2022 Dec 06: 1.07 group filtering - write out
2022 Dec 02: 1.07 group filtering - functional!
2022 Dec 02: 1.06 group filtering
2022 Dec 02: 1.05 group coo and Ri table
2022 Dec 02: 1.04 showing groups
2022 Dec 02: 1.03 fresh code for multipart features
2022 Nov 29: 1.02 continue
2022 Nov 25: 1.01 buiding to multiscan
2022 Nov 24: 1.00 origin from localbest.p 1.25
*)
updateversion = 1.00; (* defines lowest acceptable current parameter file *)
versionupperbound = 9.00; (* highest version that will be accepted *)
(* end module version *)

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

(* begin module describe.multibest *)
(*
name
   multibest: find features with the highest Ri in a given sequence range

synopsis
   multibest(scanfeatures: in, multibestp: in,
             bestfeatures: out, lowerfeatures: out, output: out)

files
   scanfeatures:  A list of locations of binding sites and other features.  

   bestfeatures:  A list of the features from the original file which have the
      highest Ri in a given range

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

      window:  Features are discarded if the coordinate of another
         feature with a higher Ri is less than or equal to the window
         distance away from the one being looked at.  Note that if
         window is 100, the region accepted will be 2*100+1 = 201
         bases.

      asterisks: If character is 'A' then the filtered out features are shown
         as asterisks in a lowerfeatures file.

   lowerfeatures:  A list of the features from the original file which have the
      lower Ri values in a given range. These are the features which have been
      'discarded' by the program.

   output: messages to the user

description

   Scan and multiscan often produce a scanfeatures file with many
   sites picked up over a small range. This makes the map very
   difficult to look at, so this program attempts to filter the
   features and find the one with the highest Ri over a small range
   and only keep that one.

   This program, multibest, is designed for use with the output of the
   multiscan program.  multiscan creates 'groups' of features that
   must be handled together by multibest.  This is in contrast to the
   original scan program that works with single features which can be
   filtered using the localbest program.

   The user should be aware that this program discards valuable
   information in the form of possible binding sites.  For this
   reason, it is a good idea to not set the window at more than 1
   helix turn (10.6 bp) for DNA sites.  Sites at this distance apart
   are probably significant and should be kept in the scanfeatures
   file.

examples

documentation

see also
   {program that scans for single part features:} scan.p
   {program to use with scan.p:} localbest.p
   {scanp and localbest are used together.  multibest may not work
   with  scan.p output features.}

   {program to use with multibest; creates multi-part features:} multiscan.p

   {Program that uses any features file to create a postscript
   {display:} lister.p

author
   Elaine Bucheimer, Tom Schneider

bugs

technical notes

   Historical note:  multibest.p was derived from localbest.p to work
   with multipart features.

   The code uses a trick that makes it extremely simple at the cost of
   a small amount of memory.

   The trick is to first pre-mark all groups to be 'kept'.  Then pass
   through all the groups (variable g) and record whether to remove a
   particular group based on a pass through all groups again (varible
   h) to mark the ones to remove.  Then the showthegroups procedure
   (which replaces the original showgroups) merely writes or does not
   write the group.  This method costs ONE boolean of memory for each
   group.  Memory is cheap!!!

   It is a little slower because groups that have been removed are
   checked again.  But this should be fixable!!  Just add an if
   statement to have h skip work on the group g that is already to be
   removed!

   Also, once a g group has been marked false, the h loop can be
   terminated!

*)
(* end module describe.multibest *)

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

      dnamax = 3000;  (* length of dna arrays *)
      namelength = 200; (* maximum key name length *)
      linelength = 210; (* maximum line readable in book *)

(* end module book.const.for.lister version = 13.51; (@ of lister.p 2018 sep 19 *)

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

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

(* 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.51; (@ of lister.p 2018 sep 19 *)

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

   KeepDuplicates = false; (* If set to true, duplicate features
      are kept, if false, they are eliminated. *)

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.51; (@ of lister.p 2018 sep 19 *)

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

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

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

(* 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: real; (* 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;

(* begin module readdefinitionsandfeaturesVariables *)
readdefinitionsandfeaturesVariables = record
(* ALL LOCAL VARIABLES of procedure readdefinitionsandfeatures
ARE CARRIED IN readdefinitionsandfeaturesVariables.
*)
   verbose: boolean; (* give verbose output *)
   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 *)
   (* 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 *)
end;
(* end module readdefinitionsandfeaturesVariables *)

(* end module definition.feature.type version = 5.89; (@ of prgmod.p 2022 Dec 20 *)

var
(* begin module feature.globals *)
(* globals needed by lister.p: *)
   groupnumber: integer; (* counted group number *)
{zzzrrr}
   localRDAFvariables: readdefinitionsandfeaturesVariables;
   (* These variables are held for the readdefinitionsandfeatures
      procedure so that it can be called repeatedly and not forget
      the states of the variables.  This is generally not to be
      touched by the outside code - they are "local" to
      readdefinitionsandfeatures. *)
(* end module feature.globals version = 5.89; (@ of prgmod.p 2022 Dec 20 *)

   scanfeatures, (* file of features to be filtered *)
   multibestp, (* parameter file to control the program *)
   bestfeatures: text; (* output file of the features having the highest Ri
                          over the window *)
   lowerfeatures: text; (* output file of the features not having the highest
                           Ri over the window. They are written to appear as
                           asterisks *)


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

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

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

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

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

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

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

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

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

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

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

(* begin module readname *)
procedure readname(var afile: text; var id: name);
(* read a name from afile, code is similar to brname.
This routine is not needed now because it is replaced by
standard string reads. *)
var
   c: char; (* a character read *)
   i: integer; (* an index *)
begin
   skipblanks(afile);
   with id do begin
      length:=0;
      repeat
         length:=succ(length);
         read(afile,c);
         letters[length] := c
      until (eoln(afile)) or
            (length>=namelength) or
            (c=' ');
      if letters[length]=' ' then length:=length-1;
      if length<namelength
         then for i:=length+1 to namelength do letters[i]:=' '
   end;
end;
(* end module readname version = 13.51; (@ of lister.p 2018 sep 19 *)

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

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

(* 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 = 8.69; (@ of lister.p 1998 February 27 *)

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

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

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

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

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

(* 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;
            o := o^.next;
         end;
         o := mothers;
         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 has no MOTHERS.');
      end;

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

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

(* begin module writefeatures *)
procedure writefeatures(var afile: text; featurelist: featuretypeptr);
(* write the features in list featurelist *)
var
   f: featuretypeptr; (* current place in the featurelist *)
begin
   f := featurelist;
   while f <> nil do with f^ do begin
      writeafeature(afile, f);
      writeln(afile);
      f := f^.next;
   end;
end;
(* end module writefeatures version = 5.89; (@ of prgmod.p 2022 Dec 20 *)

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

(* begin module writedefinitions *)
procedure writedefinitions(var afile: text; definitionlist: definitiontypeptr);
(* write the definition in list definitionlist *)
var
   f: definitiontypeptr; (* current place in the definitionlist *)
begin
   f := definitionlist;
   while f <> nil do with f^ do begin
      writeadefinition(afile, f);
      f := f^.next;
   end;
end;
(* end module writedefinitions version = 5.89; (@ of prgmod.p 2022 Dec 20 *)

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

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

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

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

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

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

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

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

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

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

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

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

(* begin module readdefinitionsandfeatures *)
procedure readdefinitionsandfeatures(var afile: text;
                       var KeepDuplicates: boolean;
{zzzccc}
                       var comments: stringptr;
                       var definitionlist: definitiontypeptr;
                       var featurelist: featuretypeptr;
                       var groupnumber: integer;
                           windowsize: integer);
{
                       var localRDAFvariables: 
                           readdefinitionsandfeaturesVariables);
}
(* Read definitions and features from afile.

These are defined in the lister.p program.

* Keep duplicate features if KeepDuplicates is true.
* comments records comment strings
* comments records comment strings
* definitionlist: definitions of features
* featurelist: list of features
* groupnumber: number of the group of features we are working on
* window: size in bases to compare features

* localRDAFvariables: is a global that acts as local variables for
  this procedure.  ALL local variables are passed through this
  variable so that the procedure can remember them.
  This allows the afile to be read in chunks separated by the
  windowsize.

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

2022 Dec 16: if window > 0 make circular group list
????

*)

{ ALL LOCAL VARIABLES ARE CARRIED IN readdefinitionsandfeaturesVariables
var
   verbose: boolean; (* give verbose output *)
   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 *)
   (* 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 with localRDAFvariables do begin
{zzzrrr}
{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; end; (* addmotherlinks *)

begin with localRDAFvariables do 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;
         countmarks(groupTotalFeature,-1); (* final check *)

         (* 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
               (* 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"
                *)

               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,'SUPERDEBUG gap is ',gap:1);
writestring(output, leftstring);writeln(output,'<-  leftstring SUPERDEBUG');
writestring(output,rightstring);writeln(output,'<- rightstring SUPERDEBUG');
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;
                  ff := ff^.next;
               end;

if superdebug then begin
write(output,'SUPERDEB leftgap: '); writeafeature(output,leftgap); writeln(output);
write(output,'SUPERDEB 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,'SUPERDEB 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,'SUPERDEB: 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,'SUPERDEB 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);
               (* now alter aGapFeature so that it points to the
                  non-duplicate feature! *)
               (* 2016oct30: KEY CHANGE to fix hindi problem! *)
               if (not aGapFeature^.definition^.isgap)
               then aGapFeature^.others^.afeature := gg^.duplicateprevious;
            end;
            gg := gg^.next;
         end;

if superdebug then begin
writeln(output,'SUPERDEBUG: 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; (* end threadgroup *)

begin with localRDAFvariables do begin (* readdefinitionsandfeatures *)
{
   if (window <= 0) then begin
}

   reset(afile);

(* initialize readdefinitionsandfeaturesVariables *)
   verbose := false;
   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;
   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);}
         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
                  if equaldefinition(dtest, d)
                     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.');
                     barout;

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

                     uniquedefinitions := pred(uniquedefinitions);
                  end;
                  if not done then begin (* step forward *)
                     dtest := dtest^.next;
                  end;
               end;
            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);
               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 > 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;
               countpetals := succ(countpetals);
            end;
         end

         else if identifier = '@'
         then begin (* read feature *)
            countfeatures := succ(countfeatures);
            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 *)
(* 2022 dec 16: I moved the code here because the parts under the
(* following if statment were identical *)
               f^.next := newfeature;
               f := f^.next;
               f^.featurenumber := countfeatures;
               f^.featuregroup := groupnumber;
               f^.featureisfirstingroup := firstfeatureofgroup;

               if not featureinlist(featurelist, newfeature, oldfeature)
               then begin
                  uniquefeatures := succ(uniquefeatures);
               end
               else begin
                  if verbose then begin
                     writeln(output,'This feature was found more than once:');
                     writeafeature(output,newfeature);
                     writeln(output);
                  end;
                  if KeepDuplicates
                  then begin
                     (* Keep all duplicates in the list! *)
                     if verbose then writeln(output,'It is being kept.');
                     f^.duplicateprevious := nil;
                  end
                  else begin
                     if verbose then writeln(output,'It will be ignored.');
                     f^.duplicateprevious := oldfeature;
                  end

               end
            end;
            if firstfeatureofgroup then firstfeatureofgroup := false;

            threadgroup;
{ttt Is this the right place to examine the list or should it be later?}
         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;

   (* 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
            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');

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; (* readdefinitionsandfeatures *)
(* end module readdefinitionsandfeatures version = 5.89; (@ of prgmod.p 2022 Dec 20 *)

(* begin module embed.readparameters *)
procedure readparameters(var multibestp: text;
                         var window: integer;
                         var asterisk: boolean);
(* read the parameters from multibestp starting with the version
number *)
var
   parameterversion: real; (* parameter version number *)
   ch: char;
begin     
   reset(multibestp);

   readln(multibestp, parameterversion);
   if (parameterversion < updateversion) or
      (parameterversion > versionupperbound) then begin
      writeln(output, 'You have an old parameter file!');
      halt   
   end;  

   readln(multibestp, window);
   if (window < 2) then begin
      writeln(output, 'window parameter must be at least 2');
      halt;
   end;
   readln(multibestp, ch);
   if ch = 'A' then asterisk := true else asterisk := false;

   writeln(output);
   write(output, 'Any features which are within '); 
   write(output, window:1);
   writeln(output, ' bp of a feature');
   write(output, 'with a higher Ri will be ');
   
   if asterisk
   then writeln(output, 'written to the lowerfeatures file')
   else writeln(output, 'deleted');
end;
(* end module embed.readparameters *)

(* begin module filterfeatures *)
procedure filterfeatures(f: featuretypeptr;
                         var filteredfeature, discardfeature: featuretypeptr;
                         window: integer);
(* looks at each feature in f and decides whether it should be kept or not.
   filteredfeature is a list to work from that begins identical to f and ends
   containing only the filtered features. discardfeature is a list of discarded
   features, or the ones with the low Ri values which are rejected. window is
   the distance within which the maximum feature is to be found. *)
var 
   filteredfeature2: featuretypeptr; (* another reference to the filteredfeature list *)
   first: boolean;             (* has the first feature in filteredfeature been found *)
   firstf: featuretypeptr;     (* first feature in f *)
   firstfi: featuretypeptr;    (* first feature in filteredfeature *)
   found: boolean;             (* has a feature to be deleted been found *)

(* testing variables *)
c: integer; (* counter *)
t: featuretypeptr; (* location in discardfeature list *)

(* for testing integrity of the discardfeature list *)
procedure testit;
begin
c := 0;
 t := discardfeature;
 while t <> nil do begin
    c := succ(c);
    writeln(output,c:1);
    t := t^.next;
 end;
end;

begin

   firstf := f; (* sets the first element of the list of features *) 
halt;
{
procedure showfeatures(var afile: text; featurelist: featuretypeptr;
                       paging: boolean);
}
writeln(output,'========== showfeatures');
showfeatures(output, f, false);
halt;

   (* filteredfeature2 is one step ahead of filteredfeature. It is the feature
      which is actually compared to the one in f *)
   new(filteredfeature2); 
   filteredfeature2 := filteredfeature^.next; 

   discardfeature := nil;  (* this has to be initiated! 2005 May 26 *)

   (* find the first feature to be saved in the filteredfeature list. The
      program scans the list for the first feature which does not have a
      feature within the window range with a higher Ri *)
   first := false; 
   while (not first) do begin (* while the first feature has not been found *)

      first := true; (* assume this is the first feature *)
      f := firstf;
      found := false;

      (* scan the features in f until they have all been compared or a
         feature with a higher Ri than the filteredfeature within the
         window is found *)

      while (f <> nil) and (not found) do begin 
         (* if a feature with a higher Ri within the window is found, then it
            is not the first feature *)

{zzzqqq}
{procedure showafeature(var afile: text; f: featuretypeptr);}
writeln(output,'groupnumber: =============', groupnumber);
writeln(output,'feature: =============');
showafeature(output, f);
writeln(output,'======================');
halt;

{
if f^.definition^.istotal then begin
writeln(output,'istotal!');
end
else begin
writeln(output,'not istotal!');
end;
}


if f^.mothers <> nil then begin
   writeln(output,'this feature has a mommy!');
   if f^.mothers^.afeature <> nil then begin
      writeln(output,'f^.mothers^.afeature <> nil');
      if f^.mothers^.afeature^.definition^.istotal then begin
         writeln(output,'has total!');
      end
      else begin
         writeln(output,'no total!');
      end;
   end;
end;

{
writeln(output,
'f^.afeature^.definition^.istotal',
f^.afeature^.definition^.istotal
);
}


         if (abs(filteredfeature^.coordinate - f^.coordinate) <= window)
               and (filteredfeature^.Ri < f^.Ri) 
               and equalstring(filteredfeature^.nametag, f^.nametag) 
               and equalstring(f^.id, filteredfeature^.id) then begin

            (* bypass the first feature. Do not save in
               the filteredfeature list *)

            filteredfeature^.next := discardfeature;
            discardfeature := filteredfeature;
            filteredfeature := filteredfeature2;
            filteredfeature2 := filteredfeature^.next;

            found := true;
            first := false;
         end;
         f := f^.next;
      end;
   end;

writeln(output,'bubba halt');
halt;

   firstfi := filteredfeature; (* sets the first element in filteredfeature *)

   while filteredfeature2 <> nil do begin
      (* there are still features to be compared *)
      f := firstf; (* sets f to the beginning of the features list *)
      found := false; (* a feature to be deleted has not been found *)

      (* compare filteredfeature2 to each feature in f until the end of f is
         found or filteredfeature2 has a higher Ri than f and is within the
         window *)
      while (f <> nil) and (not found) do begin
         if (abs(filteredfeature2^.coordinate - f^.coordinate) <= window) 
            and (filteredfeature2^.Ri < f^.Ri) 
            and equalstring(filteredfeature2^.nametag, f^.nametag) 
            and equalstring(f^.id, filteredfeature2^.id) then begin

            (* if the Ri of filteredfeature is less than the Ri of f
               and is within the window, get rid of it *)

            found := true; (* a feature to be deleted has been found *)

            filteredfeature^.next := filteredfeature2^.next;
            filteredfeature2^.next := discardfeature;
            discardfeature := filteredfeature2;
            filteredfeature2 := filteredfeature^.next;

         end;
         f := f^.next;
      end;
     
      (* if a feature in f with higher Ri within the window is not found,
         simply step along the list to the next feature to be compared *)
      if not found then begin
         filteredfeature := filteredfeature^.next;
         filteredfeature2:= filteredfeature2^.next;
      end;
   end;

   filteredfeature := firstfi; (* set the filteredfeature list at the
                                  beginning to be written out *)
end;
(* end module filterfeatures *)

(* begin module writefilteredfeatures *)
procedure writefilteredfeatures(filteredfeature: featuretypeptr;
                                var features, bestfeatures: text;
                                var filterednumber: integer);
(* write the filtered features out, and count them. *)
var
   ch: char;
begin
   rewrite(bestfeatures);
   reset(features);
   filterednumber :=0;

   (* copy all the lines not beginning with '@' from features to
     bestfeatures.  This includes all the definitions *)
   while not eof(features) do begin
      if (features^ <> '@') then begin
         while not eoln(features) do begin
            read(features, ch);
            write(bestfeatures, ch);
         end;
         writeln(bestfeatures);
      end;
      readln(features);
   end;
   
   (* write the list of filteredfeatures to the bestfeatures file *)
   while filteredfeature <> nil do begin
      writeafeature(bestfeatures, filteredfeature);
      writeln(bestfeatures);
      filteredfeature := filteredfeature^.next;
      filterednumber := succ(filterednumber);
   end;

   writeln(output, 'number filtered:  ', filterednumber:1);

end;
(* end module writefilteredfeatures *)

(* begin module writediscardedfeatures *)
procedure writediscardedfeatures(var discardfeature: featuretypeptr;
                                     definitionlist: definitiontypeptr;
                                 var features, lowerfeatures: text;
                                 var discardednumber: integer);
(* write out the discarded features, count them. *)
var
   ch: char; (* a character to copy *)
   x: integer; (* index to strings *)
   firstd: definitiontypeptr; (* pointer to first definition *)
begin
   rewrite(lowerfeatures);
   reset(features);
   firstd := definitionlist;
   discardednumber :=0;

   while definitionlist <> nil do begin
      (* redefine the definition to only display an asterisk at
         coordinate zero.  Rename the features by adding '.filtered'
         to the end  *)
      with definitionlist^.nametag do begin
         letters[length+1] := '.';
         letters[length+2] := 'f';
         letters[length+3] := 'i';
         letters[length+4] := 'l';
         letters[length+5] := 't';
         letters[length+6] := 'e';
         letters[length+7] := 'r';
         letters[length+8] := 'e';
         letters[length+9] := 'd';
         length := length + 9;
      end;
      with definitionlist^.background do begin
         for x := 1 to length do if letters[x] = 'w' then letters[x] := ' ';
      end; 
      with definitionlist^.negparts do begin
         length := 1;
         letters[length] := '*'; 
      end;
      with definitionlist^.posparts do begin
         length := 1;
         letters[length] := '*'; 
      end;
      with definitionlist^ do begin
         locations[1] := 0;
      end;
      definitionlist^.marks := 1;
      definitionlist^.min := 0;
      definitionlist^.max := 0;

      definitionlist := definitionlist^.next;
   end;

   definitionlist := firstd;

   (* copy all the lines not beginning with '@' from features to
      lowerfeatures.  This includes all the definitions *)
   while not eof(features) do begin
      if (features^ <> '@') and (features^ <> 'd') then begin
         while not eoln(features) do begin
            read(features, ch);
            write(lowerfeatures, ch);
         end;
         writeln(lowerfeatures);
      end;
      if (features^ = 'd') then begin
         while definitionlist <> nil do begin
            writeadefinition(lowerfeatures, definitionlist);
            writeln(lowerfeatures);
            definitionlist := definitionlist^.next;
         end;
      end;
      readln(features);
   end;  

   (* write the list of discardfeature to the lowerfeatures file *)
   while discardfeature <> nil do begin

      (* change nametag of feature to match nametag of definition *)
      with discardfeature^.nametag do begin
         letters[length+1] := '.';
         letters[length+2] := 'f';
         letters[length+3] := 'i';
         letters[length+4] := 'l';
         letters[length+5] := 't';
         letters[length+6] := 'e';
         letters[length+7] := 'r';
         letters[length+8] := 'e';
         letters[length+9] := 'd';
         length := length + 9;
      end;

      (* delete quotestring so the only mark shown is an asterisk *)
      (* 2005 May 26 TDS: GEEPERS!!!!
         For some reason Elaine just DELETED the othertag!
         discardfeature^.othertag.length := 0;
         Well keep it!
      *)

      writeafeature(lowerfeatures, discardfeature);
      writeln(lowerfeatures);
      discardfeature := discardfeature^.next;
      discardednumber := succ(discardednumber);
   end;  

   writeln(output, 'number discarded: ', discardednumber:1);

end;
(* end module writediscardedfeatures *)

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

(* begin module showthegroups *)
procedure showthegroups(var afile: text;
                        var groupslist: featurecarrierptr;
                        keep: boolean;
                        var countgroups: integer);
(* Show the groups in grouplist to the afile.
If keep is true then show those where keepgroup is true
If keep is falsethen show those where keepgroup is false
countgroups counts the groups written out.
*)
var
   g: featurecarrierptr; (* current pointer into groupslist *)
   f: featuretypeptr; (* a feature in g *)
begin
   countgroups := 0;
   g := groupslist;
   (* Write definitions only for the FIRST group so as not to
      duplicate them! *)
   f := g^.afeature;
   while f <> nil do begin
      if f^.definition^.matrix <> nil then begin
         writeadefinition(afile,f^.definition);
         putriblmatrix(afile,f^.definition^.matrix);
      end;
      f := f^.next
   end;
   f := g^.afeature;
   while g <> nil do begin
      groupnumber := g^.afeature^.featuregroup;
      writeln(afile);
      if g^.keepgroup = keep then begin
         writeln(afile,'group ',groupnumber:1);
         countgroups := succ(countgroups);
         (* Do each feature separately to make sure only the
            relevant definitions are displayed! *)
         f := g^.afeature;
         while f <> nil do begin
            writeadefinition(afile,f^.definition); {writeln(afile);}
            writeafeature(afile,f); writeln(afile);
{zzzaaa}
            f := f^.next
         end;
      end
      else writeln(afile,'* group ',groupnumber:1,' removed by multibest');
      g := g^.next
   end;
end;
(* end module showthegroups *)

(* begin module showgroups *)
procedure showgroups(var afile: text;
                     var groupslist: featurecarrierptr);
(* Show the groups in grouplist to the afile *)
var
   g: featurecarrierptr; (* current pointer into groupslist *)
   f: featuretypeptr; (* a feature in g *)
begin
   g := groupslist;
   (* Write definitions only for the FIRST group so as not to
      duplicate them! *)
   f := g^.afeature;
   while f <> nil do begin
      if f^.definition^.matrix <> nil then begin
         writeadefinition(afile,f^.definition);
         putriblmatrix(afile,f^.definition^.matrix);
      end;
      f := f^.next
   end;
   f := g^.afeature;
   while g <> nil do begin
      groupnumber := g^.afeature^.featuregroup;
      writeln(afile);
      writeln(afile,'group ',groupnumber:1);
      (* Do each feature separately to make sure only the
         relevant definitions are displayed! *)
      f := g^.afeature;
      while f <> nil do begin
         writeadefinition(afile,f^.definition); {writeln(afile);}
         writeafeature(afile,f); writeln(afile);
         f := f^.next
      end;
      g := g^.next
   end;
end;
(* end module showgroups *)

(* begin module groupfeatures *)
procedure groupfeatures(var featurelist: featuretypeptr;
                        var groupslist: featurecarrierptr);
(* Group the features from featurelist together in group carriers
called the grouplist. *)
const debugging = false; (* set to true to see results *)
var
   f: featuretypeptr; (* one feature in f *)
   g: featurecarrierptr; (* current pointer into groupslist *)
   p: featuretypeptr; (* the previous feature in featurelist *)
   groupnumber: integer; (* group number *)
begin
{
writeln(output,'*** begin groupfeatures');
definitions for reference:
   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 *)
      next: featurecarrierptr; (* the next carriage *)
   end;
}
   if featurelist = nil then begin
      writeln(output,'groupfeatures was called with an empty feature list');
      halt
   end;

   (* initialize the groupslist *)
   new(groupslist);
   groupslist^.afeature := featurelist;
   groupslist^.next := nil;
   groupnumber := featurelist^.featuregroup;
   g := groupslist;

   f := featurelist;
   p := nil;
   while f <> nil do begin

      if debugging then begin
         writeln(output,'*********** for group ',groupnumber:1);
         writeln(output);
{
         showafeature(output,f);
}
      end;
{zzznnn}

      if f^.featuregroup <> groupnumber then begin
if debugging then writeln(output,'****** new group found');
         p^.next := nil; (* disconnect and finish previous group list *)
         groupnumber := f^.featuregroup;
         new(g^.next);
         g := g^.next;
         g^.afeature := f;
         g^.next := nil;
         (* complete the previous list if it exists *)
         if p <> nil then begin
            p^.next := nil;
         end;
      end
      else begin 
         if f^.definition^.istotal then begin
            (* This is the end of the group.
               Record the coo and Ri into the carrier *)
            g^.coo := f^.coordinate;
            g^.Ri  := f^.Aparam; (* since the bits is in Aparam, not Ri *)
if debugging then writeln(output,'*********** IS TOTAL');
{ this code checked the total
(* Note: the total value is not in f^.Ri but rather F^.Aparam *)
writeafeature(output,f);
writeln(output);
writeln(output,f^.Ri:1:infodecim,' bits');
writeln(output,'Ri=',f^.Ri:1:infodecim,' bits WTF????');
writeln(output,'Ri for total is in Aparam=',f^.Aparam:1:infodecim,' bits');
}
         end;
      end;
      p := f; (* hold the previous *)
      f := f^.next
   end;

   if debugging then begin
      writeln(output,'**** debugging groupfeatures');
      showgroups(output, groupslist);
   end;

if debugging then writeln(output,'*** done groupfeatures');
end;
(* end module groupfeatures *)

(* begin module multibest.themain *)
procedure themain(var scanfeatures, multibestp, bestfeatures, lowerfeatures: text);
(* the main procedure of the program *)
procedure myhalt; (* for debugging *)
begin
writeln(output,'SPECIAL HALT');
{
showgroups(bestfeatures, groupsoffeatures);
}
writeln(output, 'no goto');
{goto 1}
end;
var
   asterisk: boolean; (* should discarded features be written as asterisks
                         in the lowerfeatures file *)
   parameterversion: real; (* parameter version number *)
   definitionlist: definitiontypeptr; (* list of definitions *)
   discardfeature: featuretypeptr; (* list of discarded features *)
   discardednumber: integer; (* number of discarded features *)
   f: featuretypeptr; (* list of features from feature file *)
   filteredfeature: featuretypeptr; (* list of features,
                                       same as f, to work from *)
   filterednumber: integer; (* number of filtered features *)
   window: integer; (* distance within which the maximum feature
                    is to be found *)
   comments: stringptr; (* comment list from afile *)
   KeepDuplicates: boolean; (* keep duplicate features? *)

   groupsoffeatures: featurecarrierptr; (* carry features collected by groups *)
   g: featurecarrierptr; (* one group in groupsoffeatures *)
   preg: featurecarrierptr; (* the group previous to g *)
   h: featurecarrierptr; (* another group in groupsoffeatures *)

   debugging: boolean; (* turnon/off debugging *)
   Scountgroups: integer; (* number of groups counted for scanfeatures *)
   Bcountgroups: integer; (* number of groups counted for bestfeatures *)
   Lcountgroups: integer; (* number of groups counted for lowerfeatures *)
   Rcountgroups: integer; (* number of groups removed *)

begin
   writeln(output,'multibest ',version:4:2);

   (* swap these statements to turn on/off debugging: *)
   debugging := true;
   debugging := false;

   PatentMessage(output);

   reset(multibestp);
   readln(multibestp, parameterversion);
   if parameterversion < updateversion then begin
      writeln(output, 'You have an old parameter file!');
      halt
   end;

   readparameters(multibestp, window, asterisk);

   (* swap these statements to turn on/off debugging: *)
   KeepDuplicates := false;
   KeepDuplicates := true;

   reset(scanfeatures); (* required since window > 0 *)
writeln(output,'about to read defs with window = ',window:1);
   readdefinitionsandfeatures(scanfeatures, KeepDuplicates, comments,
                              definitionlist, f, groupnumber, window);
writeln(output,'done read defs');
{
halt;
}

   groupfeatures(f, groupsoffeatures);
{ use this to check groupfeatures:
showgroups(output, groupsoffeatures);
}
   writeln(output);

   rewrite(bestfeatures);
   writeln(bestfeatures,'* multibest ',version:4:2,' - bestfeatures');
   writestringlines(bestfeatures, comments);
   writeln(bestfeatures);

   rewrite(lowerfeatures);
   writeln(lowerfeatures,'* multibest ',version:4:2,' - lowerfeatures');
   writestringlines(lowerfeatures, comments);
   writeln(lowerfeatures);

   Scountgroups := 0;
   Rcountgroups := 0;
   g := groupsoffeatures;
   while g <> nil do begin
      g^.keepgroup := true; (* hold on to all groups initially *)
      g := g^.next;
   end;
   g := groupsoffeatures;
   while g <> nil do begin
      Scountgroups := succ(Scountgroups);
if debugging then writeln(output,'BUBBA 1 g=',ord(g):1, ' coo=',g^.coo:1);
      (* compare this group g to all other groups h *)
      h := groupsoffeatures;
      while (h <> nil) do begin
         (* Trick #3 don't look at equal cases! *)
         if g <> h then begin
if debugging then writeln(output,'BUBBA 2 h=',ord(h):1, ' coo=',h^.coo:1);
            (* Trick #2 don't look at ones already marked false! *)
            if g^.keepgroup then begin
               (* is h within the window to g? *)
               if abs(g^.coo - h^.coo) <= window then begin
                  (* is the information of g less than that of h? *)
                  if g^.Ri < h^.Ri then begin
                     g^.keepgroup := false; (* bu-bye! *)
                     (* Trick #1 stop looking at this g: already marked false! *)
                     h := nil; (* drop out of the loop!  Speed up!! *)
                  end
                  else begin (* look for duplicate features *)
                     if g^.coo = h^.coo then begin
                        if g^.Ri = h^.Ri then begin
                           (* g is NOT h and yet they have identical Ri! *)
                  { (* record comparison is dicey because they are binary,
                    so don't rely on it *)
                           if (g^.afeature == h^.afeature) then begin
                              writeln(output,'record comparson is =');
                              g^.keepgroup := false; (* bu-bye! *)
                           end
                           else begin
                              writeln(output,'record comparson is <>');
                           end;
                  }
                           if g^.keepgroup and h^.keepgroup then begin
                              (* only remove ONE of the two!! *)
                              g^.keepgroup := false; (* bu-bye! *)
                              h := nil; (* drop out of the loop!  Speed up!! *)
                              Rcountgroups := succ(Rcountgroups);
{
}
                           end;
                        end
                     end
                  end (* Ri matching and duplicate search *)
               end; (* coordinate matching *)
            end;
         end; (* g <> h test *)
         if h <> nil then h := h^.next; (* prevent crash when dropping out *)
      end; (* h loop *)
      g := g^.next;
   end; (* g loop *)

   showthegroups(bestfeatures,  groupsoffeatures, true,  Bcountgroups);
   showthegroups(lowerfeatures, groupsoffeatures, false, Lcountgroups);

   writeln(output, Scountgroups:3, ' groups were ... in scanfeatures');
   writeln(output, Rcountgroups:3, ' groups were ...... removed as duplicates');
   writeln(output, Bcountgroups:3, ' groups were put in bestfeatures');
   writeln(output, Lcountgroups:3, ' groups were put in lowerfeatures');

   if Scountgroups <> (Bcountgroups + Lcountgroups) then begin
      writeln(output, 'counts are wrong');
   end;

   writeln(output,'multibest ',version:4:2,' is done');

end;
(* end module multibest.themain *)

begin
   themain(scanfeatures, multibestp, bestfeatures, lowerfeatures);
1: end.
