program mfoldfea(cootab, {mfoldfeap,} book, mfoldfeatures, output);
(* mfoldfea: convert mfold RNA structure to lister feature

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

 *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 1.35; (* of mfoldfea.p 2007 mar 15
2007 mar 15, 1.34; segmentation fault - needed prgmod module upgrade
2005 May 25, 1.33; add patent message, this DOES use ii
2005 May 23, 1.32: add patent message? No, this doesn't use ii
2000 Oct 12, 1.28: fix bug: more than 10 sequences handled properly
2000 Oct  6, 1.25: only handles last sequence bug
2000 Oct  1, 1.24: piece orientation bug fixed
2000 Jul 30, 1.23: upgrade to GPC
2000 Jun  2, 1.22: debug
1999 Aug  8, 1.21: functional version
origin 1999 July 4 *)
updateversion = 1.00; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.mfoldfea *)
(*
name
   mfoldfea: convert mfold RNA structure to lister feature

synopsis
   mfoldfea(cootab: in, book: in,
            mfoldfeatures: out, output: out)

files
   cootab:  Concatenation of the .ct files from mfold that correspond
      to structures in the book.

   mfoldfeatures:  features for the lister program.

   book:  the Delila book used by mfoldseq to make the mfoldsequ, which was
       then used by mfold to make the cootab (coordinate table) file.

   output: messages to the user

description

   Convert the mfold RNA structures to features for the lister program.

   Mfold will only fold one sequence at a time.  So the mfoldseq program is
   used a series of times inside the run.mfold script to get all the folds in
   a book.  These are gathered together by run.mfold into the cootab file.
   This program then takes the cootab file along with the original book to
   generate features.  It uses the original book to get the coordinate
   system, which was lost in converting to the mfold file format.

examples

   Starting from a book, run mfoldseq to get a sequ file that contains one
   sequence.  For the second sequence in the mfold.book example, mfoldseq gives
   this mfoldsequ file:

;
mfold.demo.2#2
AAAACCAGGCGCAAATTTTTTTTGCGCCGGG1

Mfold generates the mfoldsequ.ct

   31 ENERGY =   -14.5    mfold.demo.2#2                
    1 A       0    2    0    1
    2 A       1    3    0    2
    3 A       2    4    0    3
    4 A       3    5    0    4
    5 C       4    6   30    5
    6 C       5    7   29    6
    7 A       6    8    0    7
    8 G       7    9   28    8
    9 G       8   10   27    9
   10 C       9   11   26   10
   11 G      10   12   25   11
   12 C      11   13   24   12
   13 A      12   14   23   13
   14 A      13   15   22   14
   15 A      14   16   21   15
   16 T      15   17    0   16
   17 T      16   18    0   17
   18 T      17   19    0   18
   19 T      18   20    0   19
   20 T      19   21    0   20
   21 T      20   22   15   21
   22 T      21   23   14   22
   23 T      22   24   13   23
   24 G      23   25   12   24
   25 C      24   26   11   25
   26 G      25   27   10   26
   27 C      26   28    9   27
   28 C      27   29    8   28
   29 G      28   30    6   29
   30 G      29   31    5   30
   31 G      30    0    0   31

which is input to this program in the cootab file, along with the original
book.

The structure, from mfoldsequ.out.html, is:

 FOLDING BASES    1 TO   31 OF mfold.demo.2#2                
 ENERGY  =     -14.5
 
        10        
AAAA| A        TT 
    CC GGCGCAAA   
    GG CCGCGTTT  T
---G^ -        TT 
   30         20  

When the mfoldfeatures are given to lister, the resulting list is:

        *10       *         *20       *         *30       *       
 5' a a a a c c a g g c g c a a a t t t t t t t t g c g c c g g g 3' 
            (-(---(-(-(-(-(-(-(-(-----------)-)-)-)-)-)-)-)-)-) helix.mfold.demo.2#2[8.30]-14.5

The parenthesis show bases that are paired together.

documentation

see also

   {The script that runs mfoldseq and mfoldfea:} run.mfold
   {The program that pulls sequences from the book for run.mfold:} mfoldseq.p
   {The program that uses the mfoldfeatures result:} lister.p
   {Small book for demonstrating folding:} mfold.book

   {More documentation on mfold is at Michael Zuker's Web site:}
   http://www.rpi.edu/~zukerm/

author
   Thomas Dana Schneider

bugs

technical notes

   The mfold program has no clue about coordinate systems as in the Delila
   system, it works with from 1 to n numbering like GenBank.  As a result the
   folded coordinates need to be mapped to the book coordinates so that the
   features can be written properly.  This is a little tricky because there
   can be multiple folds per sequence.  So the mfoldseq program extracts a
   single sequence and numbers it like EG12#18.  This name passes through
   mfold unscathed.  The mfoldfea program then separates apart the EG12 part
   as the name from the 18 as the book number.  When mfoldfea is writing the
   features, both the name and the sequence number must match to the book.
   The number is just counted, it is independent of the number that a user
   can assign with the 'set numbering' command in Delila (<NUMBERING
   DEFAULT>).  The coordinates then are taken from the book.

*)
(* end module describe.mfoldfea *)

(* begin module ct.documentation *)
(*
From zuker@snark.wustl.edu Thu Jul  8 19:52 EDT 1999
Date: Thu, 8 Jul 1999 18:51:38 -0500 (CDT)
From: Michael Zuker <zuker@snark.wustl.edu>
Message-Id: <199907082351.SAA03602@snark.wustl.edu>
To: Tom Schneider <toms@alum.mit.edu>
Subject: Re: stuck
Cc: +tom@snark.wustl.edu

Tom:

OK, I get it now. You really want to work with the 'ct' file. It
gives you a reasonable set of foldings that you can work with.

'ct' stands for 'coordinate table'. I got the name and format from
Richard Feldmann. The format is fairly clear, but just in case:

1st record: number of bases, followed by a label. The energy is
usually given, but this is not strictly necessary. My software
assumes that the 'equal sign', '=', is always present, and uses
this information to chop up the file into individual foldings.

subsequent records:

  143 C     142  144  201  143
  144 G     143  145  200  144
  145 C     144  146  199  145
  146 G     145  147    0  146
  147 A     146  148  197  147

  i seq(i) 5'(i) 3'(i) bp(i) hst_num(i)

The ith record after the first must have i as the first number. This
is an error check. seq(i) is an ASCII character that represents the
base. There a no restrictions on what characters may be used. 5'(i) is
the upstream covalently connected base. Usually i-1, but may be 0 to
indicate that there is no connection. 3'(i) is the downstream
covalently connected base. It's usually i+1, but may be 0 to indicate
that there is no connection. Thus the 'ct' file format automatically
handles multi-molecular folding. It also handles circular molecules,
since you can set 5'(1) = n and 3'(n) = 1, where 'n' is the number of
bases in the molecule. bp(i) and i are base paired, except when bp(i)
is 0. In this case, i is single-stranded. Note that base triples
cannot be represented, but that pseudoknots are no problem. hst_hum(i)
is the historical numbering for base i. This can be anything. It could
represent genomic DNA numbering before processing. Negative integers
are allowed. The numbering may decrease. You could have something
like: "1, 2, 3, ... 10, 1, 2, ... 10" that might represent 2 10-mers
that are hybridized.

Michael
*)
(* end module ct.documentation *)

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

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

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

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

type

(* begin module string.type *)
   stringptr = ^string; (* pointer to a string *)
   string = record (* a string of characters *)
      letters: array[1..maxstring] of char; (* the letters in the string *)
      length: integer; (* the number of characters in the string *)
      current: integer; (* the letter we are working on *)
      previous: stringptr; (* the previous string in a series *)
      next: stringptr; (* the next string in a series *)
   end;
(* end module string.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module 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.52; {of delmod.p 2000 Jul 30} *)

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


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

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

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


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

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

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

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

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

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

(* end module book.type version = 7.52; {of delmod.p 2000 Jul 30} *)


(* begin module wave.type *)
      waveptr = ^waveparam; (* to link several wave definitions together *)
      waveparam = record (* parameters to define a cosine wave *)
         (* define a cosine wave: *)
         extreme: char; (* h or l, the high or low extreme to be defined *)
         wavelocation: real; (* the location in bases of the extreme *)
         wavebit: real; (* the location in bits of the extreme *)
         waveamplitude: real; (* the amplitude of the wave in bits *)
         wavelength: real; (* the wave length of the wave in bases *)
{ 2000 Nov 1.  This is difficult because the length of a curve is not
simply multiplied by a scale factor.
         BASEdashon: real; (* dashon interval in bases *)
         BASEdashoff: real; (* dashon interval in bases *)
         BASEdashoffset: real; (* dashon interval in bases *)
}
         dashon: real; (* the size of on dashes in cm.
                          dashon <= 0 means no dashes *)
         dashoff: real; (* the size of on dashes in cm *)
         dashoffset: real; (* the size of off dashes in cm *)
         thickness: real; (* thickness of the wave in cm.  <= 0 means default *)
         next: waveptr; (* the next wave *)
      end;
(* end module wave.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module scan.type *)
   (* define the ribl matrix.  Source library is prgmod.p
      REQUIRES ribl.const module
      REQUIRES base.type or book.type module
      REQUIRES wave.type module for wave definition! *)
   rblarray = array[a..t, minribl..maxribl] of real; (* real(B,L) *)
   ribltype = record
      riblname: string; (* name of the weight matrix *)
      riblheader: stringptr; (* misc header information, starting with '*' *)
      data: rblarray; (* real(B,L) *)
      numbers: array[a..t, minribl..maxribl] of integer; (* n(b,l) *)
      frombase, tobase: integer; (* range of the matrix *)
      mean, stdev: real; (* mean and standard deviation of distribution *)
      consensus: real; (* Ri value of consensus sequence *)
      anticonsensus: real; (* Ri value of anticonsensus sequence *)
      averageRi: real; (* average Ri for random sequence *)
      n: integer; (* number of sequences used to create matrix *)
      symmetry: char; (* the symmetry of the Riblmatrix *)
      waves: waveptr; (* cosine wave definitions for the matrix *)
      cmperbase, cmperbit: real; (* required for wave dash definition *)
      Ribound: real; (* the Ri boundary for this definition *) 
      Zbound: real; (* the Z boundary for this definition *) 
      Pbound: real; (* the probability boundary for this definition *) 
   end;
   ribltypeptr = ^ribltype; (* pointer to ribltype so it does
      not need to be created unless actually used *)
(* end module scan.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

   definitiontypeptr = ^definetype; (* pointer to definetype *)
   definetype = 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 *)
      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 *)
      coordinate: real; (* location of zero base *)
      orientation: integer; (* -1 or +1 is in direction of piece *)
      nametag: string; (* the name tag - name of the feature *)
      othertag: string; (* the other tag - another tag for misc use *)
      (* ---------------------------- *)
      definition: definitiontypeptr; (* the definition of the feature *)
      (* ---------------------------- *)
      (* the following three values are only defined when there
         is a definition for the feature matrix *)
      Ri: real; (* individual information of this site *)
      Z: real; (* Z score for this site *)
      probability: real; (* probability for this site *)

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


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

(* begin module mfold.type *)
   pairingptr = ^pairing; (* pointer to a ctdata record *)
   pairing = record (*  Define one coordinate table record.
      See module ct.documentation for details.
         147 A     146  148  197  147
         i seq(i) 5'(i) 3'(i) bp(i) hst_num(i)
      *)
      i: integer; (* position in the sequence *)
      thebase: char; (* the base at i *)
      five: integer; (* the 5' connection to base i *)
      three: integer; (* the 3' connection to base i *)
      basepaired: integer; (* the base paired to i *)
      hstnum: integer; (* historical numbering ("external coordinates") *)
      next: pairingptr; (* the next ctdata record *)
   end;
   structureptr = ^structure; (* pointer to a single RNA structure *)
   structure = record (* a single RNA structure *)
      idnumber: string; (* the sequence identifier for this structure
          with a number on the end, '#5' as read from cootab *)
      sequnumber: integer; (* the number '#5' as a number *)
      id: string; (* the sequence identifier for this structure *)
      nametag: string; (* the identifier for this structure *)
      length: integer; (* the length of this structure *)
      energy: real; (* the energy for this folding, in kcal/mole *)
      pairs: pairingptr; (* the ctdata records *)
      orientation: integer; (* orientation of the sequence*)
      next: structureptr; (* the ctdata records *)
   end;
(* end module mfold.type *)

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

var
   cootab, (* file used by this program *)
   book, (* file used by this program *)
   mfoldfeatures: text; (* file used by this program *)

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

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

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

(* ************************************************************************ *)
(* end module book.var version = 7.52; {of delmod.p 2000 Jul 30} *)


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

(* begin module crash *)
procedure crash;
(* Crash the program by trying to open a nonexistant file.  This allows
tracing by the dbx program.  To use:  insert call into the halt program
or whereever a traceable stop is desired. *)
var
   bogus: text; (* boghous internal file *)
begin
  writeln(output,' program crash.');
  reset(bogus);
end;
(* end module crash version = 7.52; {of delmod.p 2000 Jul 30} *)

(* begin module halt *)
procedure halt;
(* stop the program.  the procedure performs a goto to the end of the
   program.  you must have a label:
      label 1;
   declared, and also the end of the program must have this label:
      1: end.
   examples are in the module libraries.
   this is the only goto in the delila system. *)
begin
      writeln(output,' program halt.');
      goto 1
end;
(* end module halt version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

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

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

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

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

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

procedure writeline(var afile: text; l: lineptr;
                        carriagereturn: boolean);
(* write a line to a file, with carriage return if
carriagereturn is true. *)
var
   index: integer; (* index to characters in l *)
begin
   with l^ do begin
      for index := 1 to length
      do write(afile, letters[index]);
   end;
   if carriagereturn then writeln(afile);
end;

procedure showfreedna;
(* show the freedna list *)
var
   counter: integer; (* count of freedna list *)
   l: dnaptr; (* pointer into freedna list *)
begin
   l := freedna;
   counter := 0;
   while l <> nil do begin
      counter := succ(counter);
      write(output,counter:1);
      write(output, ', length = ',l^.length:1);
{
This is illegal according to gpc because one cannot
write a pointer to a text file.  It can be unearthed
for debugging.
      write(output, ', pointer id: ',l:1);
}
      writeln(output);
      l := l^.next
   end;
end;

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

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

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

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

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

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

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

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

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

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

(* end module book.basis version = 7.52; {of delmod.p 2000 Jul 30} *)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(* ************************************************************************ *)
(* end module package.brpiece version = 7.52; {of delmod.p 2000 Jul 30} *)

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

(* ************************************************************************ *)
(* end module package.getpiece version = 7.52; {of delmod.p 2000 Jul 30} *)

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

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

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

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

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

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

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

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

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

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

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

(* begin module attachcharacterstring *)
procedure attachcharacterstring(c: char; var a: string);
(* attach the character c to the end of a *)
begin
   if a.length = maxstring then begin
      writeln(output,'attachcharacterstring: maxstring exceeded');
      halt
   end;
   with a do begin
      length := succ(length);
      letters[length] := c
   end
end;
(* end module attachcharacterstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

(* begin module numberdigit *)
function numberdigit(number, logplace:integer): char;
(* return the digit at the place value ('logplace') position of number.
example:
      numberdigit(13625, 3) = 3
      numberdigit(13625, 4) = 1
2000 July 30 'myabsolute' replaced 'absolute', which
is apparently a keyword for GPC.  The name is kept
to keep the code looking similar to its origin.
*)
var
      place: integer; (* the exponent of logplace *)
      count: integer; (* used to make place *)
      myabsolute: integer; (* the absolute value of number *)
      acharacter: char; (* the character to be returned *)
procedure digit;
(* extract a digit at the place position *)
var
      tenplace: integer; (* ten times place *)
      z: integer; (* an intermediate value *)
      d: integer; (* the digit extracted *)
begin (* digit *)
      tenplace:=10*place;
      z:=myabsolute-((myabsolute div tenplace)*tenplace);
      if place = 1
         then d:=z
         else d:= z div place;
      case d of
         0: acharacter:='0';
         1: acharacter:='1';
         2: acharacter:='2';
         3: acharacter:='3';
         4: acharacter:='4';
         5: acharacter:='5';
         6: acharacter:='6';
         7: acharacter:='7';
         8: acharacter:='8';
         9: acharacter:='9';
      end
end; (* digit *)
procedure sign;
(* put a negative sign out or a positive sign *)
begin (* sign *)
      if number <0 then acharacter:='-'
                   else acharacter:='+'
end; (* sign *)
begin (* numberdigit *)
      place:=1;
      for count:=1 to logplace do place:=10*place;

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

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

(* begin module numberrealstring *)
procedure numberrealstring(var s: string;
                            n: real;
                            thewidth: integer;
                            thedecimal: integer);
(* Put the real number n into the string s taking up thewidth characters and
using thedecimal decimal places.  Note that the procedure multiplies n by a
perfectly integer power (which is quite clean) and then rounds that.  This
ensures that the last decimal place is rounded but no others.  Then the
program picks apart that clean number using numberdigit.  NOTE: If thewidth
is too small to fit the number then the size is increased.
Note:  The general rule is to throw away leading zeros EXCEPT in the
special case like 0.1!  This code does not account for that case.
2005 May 30: That's not acceptable.  No leading spaces are allowed
unless the string is shorter than thewidth.  Add a bit of code to do that.
*)
{
const dndebugpoint = -20; updebugpoint = -0;
const dndebugpoint = -1.1; updebugpoint = 100;
}
var
   decimalplace: integer; (* amount of space for a decimal place *)
   d: integer; (* one of the digit places *)
   digits: integer; (* the number n is multiplied by powers of 10 until
      all relevant digits are above zero.  Then it is rounded to produce
      'digits'.  This guarantees that the rounding is correct over the
      entire number. *)
   i: integer; (* index to the digits of the number *)
   power: integer; (* powers of 10 *)
   minus: boolean; (* whether we made a - *)
   newwidth: integer; (* the larger width to fit the number into *)
   sign: integer; (* amount of space for a sign *)
   symbol: char; (* a symbol to put into s *)
begin
   (* numberrealstring MUST NOT CALL clearstring(s); because
      clearstring WIPES OUT the rest of a linked list *)
   emptystring(s);
   if thewidth > 0 then with s do begin
      power := 1;
      for i := 1 to thedecimal do begin
         power := power * 10;
      end;
      digits := round(power*abs(n));
      if thedecimal > 0 then decimalplace := 1
                        else decimalplace := 0;
      if n < 0 then sign := 1
               else sign := 0;
      newwidth := numbersize(digits)+thedecimal+sign;

      (* add space for leading zero for numbers below abs(1).
      2005 May 30: This puts an extra space
      at the start of 0.97 when using only one decimal place.
      This must be removed later in some cases *)
      if abs(n) < 1.0 then newwidth := succ(newwidth);

      (* use the requested width if it is bigger *)
      if thewidth > newwidth
      then newwidth := thewidth;

{
if (n > dndebugpoint) and (n < updebugpoint) then begin
writeln(output);
write(output,'n = ',n:1:1);
write(output,'  n = ',n:thewidth:(thedecimal+3));
write(output,'  digits = ',digits:1);
write(output,'  trunc(n) = ',trunc(abs(n)):1);
write(output,'  numbersize(trunc(abs(n))) = ',numbersize(trunc(abs(n))):1);
writeln(output);
write(output,'  newwidth = ',newwidth:1);
write(output,'  thewidth = ',thewidth:1);
write(output,'  thedecimal = ',thedecimal:1);
writeln(output);
end;
}

      for d := newwidth-1-decimalplace downto 0 do begin
          symbol := numberdigit(digits,d);
{
if (n > dndebugpoint) and (n < updebugpoint) then begin
write(output,'        d = ',d:1:1);
write(output,'        symbol = ',symbol:1);
writeln(output);
end;
}
          if symbol = '-' then symbol := '0';
          if symbol = '+' then symbol := '0';
          if symbol = ' ' then symbol := '0';
          length := succ(length);
          letters[length] := symbol;
          if d = thedecimal then if thedecimal > 0 then begin
{
if (n > dndebugpoint) and (n < updebugpoint) then begin
write(output,'            at the decimal    = .');
writeln(output);
end;
}
             length := succ(length);
             letters[length] := '.';
          end;
      end;
{
if (n > dndebugpoint) and (n < updebugpoint) then begin
writequotestring(output,s); writeln(output);
end;
}

      (* convert leading zero's to blanks *)
      current := 1;
      while (letters[current] = '0') and (current < length) do begin
         if letters[current+1] <> '.' then letters[current] := ' ';
         current := succ(current);
{
if (n > dndebugpoint) and (n < updebugpoint) then begin
writequotestring(output,s); writeln(output);
end;
}

      end;
      (* now put the minus sign back *)
      if n < 0 then begin
         current := 1;
         minus := false;
         while letters[current] = ' ' do begin
            if letters[current+1] <> ' ' then begin
               letters[current] := '-';
               minus := true;
            end;
            current := succ(current);
            if current >= maxstring then begin
               write  (output,'numberrealstring:');
               writeln(output,' cannot place minus sign in string');
               writeln(output,'placing number at ',current:1);
               writeln(output,'the number is n =',n:newwidth:thedecimal);
               writestring(output,s);
               halt;
            end;
         end;
         if not minus then begin
            writeln(output,'numberrealstring ERROR:',
                           ' missing minus sign for number ',
                           n:newwidth:thedecimal);
            writeln(output,'Increase one of the number widths.');
            halt
         end;
      end;
   end;

   (* 2005 May 30: At this point the number 0.97 with one decimal
   place comes out as: " 1.0".  This is not acceptable, so remove
   any leading blanks. *)
   if s.length > thewidth
   then with s do while letters[1] = ' ' do begin
      for i := 1 to length - 1
      do letters[i] := letters[i + 1];
      length := length - 1;
   end;

{
if (n > dndebugpoint) and (n < updebugpoint) then begin
writeln(output,'s.length = ',s.length:1);
writeln(output,'n=',n:thewidth:thedecimal);
writequotestring(output,s); writeln(output);
end;
}

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

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

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

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

(* 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;
      desiredline := 0;
      others := nil;
      mothers := nil;
      usedInGroup := false; (* used in group *)
      next := nil;
   end
end;
(* end module clearfeature version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

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

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

(* begin module writeadefinition *)
procedure writeadefinition(var afile: text; d: definitiontypeptr);
(* write the definition d to file afile *)
var i: integer; (* index to the string s *)
begin
   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 do this here - it's awkward if one wants
         to just write the parts above.  Just rember to write the
      matrix after the definition! *)
      writeln(afile);
      putriblmatrix(afile, matrix);
}
   end;
end;
(* end module writeadefinition version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

(* 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
   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);
      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(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(afile,' ',Aparam:1:infodecim);
            write(afile,' ',Bparam:1:infodecim);
            write(afile,' ',Cparam:1:infodecim);
            write(afile,' ',Dparam:1:infodecim);
         end;
         write(afile,' ',relcoo:1:1); (* relative coordinate *)
         write(afile,' #',featurenumber:1); (* feature number *)
      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.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module skipblanks *)

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

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

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

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

procedure skipcolumn(var thefile: text);
(* skip over a data column *)
begin
   skipblanks(thefile); skipnonblanks(thefile)
end;
(* end module skipblanks version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

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

(* begin module unnumber *)
procedure unnumber(idnumber: string; var id: string; var n: integer);
(* remove the number id (eg '#5') from the end of idnumber
and put into id. *)
var
   c: char; (* a character in the idnumber *)
   d: integer; (* the digit corresponding to c *)
   power: integer; (* the power of 10 for the digit of the number *)
begin
{
   write(output,'unnumber, idnumber is');
   writestring(output,idnumber);
   writeln(output);
}
   n := 0;
   power := 1;
   copystring(idnumber, id);
   while (id.length > 0) and (id.letters[id.length] <> '#') do begin
      (* suck off the number *)
      c := id.letters[id.length];
      d := ord(c) - ord('0');
      if (d < 0) or (d > 9) then begin
         writeln(output,'unnumber: the id string has a bad number after the "#":');
         write(output,'The entire id string is "');
         writestring(output,id);
         writeln(output,'"');
         halt;
      end;
      n := n + d*power; (* tack the last digit on *)
      power := 10*power;
{writeln(output,'c=',c, ' d= ',d:1, ' n = ',n:1);  }
      (* move back one position in the string: *)
      id.length := pred(id.length);
   end;
   id.length := pred(id.length);
   if id.length <= 0 then begin
      writeln(output,'unnumber: id does not have a "#"');
      writestring(output,idnumber);
      writeln(output);
      halt;
   end;
{writeln(output,'THE sequnumber is ',n:1);}
end; (* unnumnber *)
(* end module unnumber *)

(* begin module readcoordinatetable *)
procedure readcoordinatetable(var cootab: text; var st: structureptr);
(* read the coordinate table ct from file cootab into the st datastructure.
Example table:
   27 ENERGY =   -15.8    lac                           
    1 C       0    2   26    1
    2 C       1    3   25    2
    3 A       2    4    0    3
    ...
   22 G      21   23    6   22
   23 C      22   24    5   23
   24 C      23   25    4   24
   25 G      24   26    2   25
   26 G      25   27    1   26
   27 G      26    0    0   27
See module ct.documentation for details.
*)
var
{
   equals: char; (* the "=" in the cootab *)
}
   h: structureptr; (* pointer into the st *)
   c: pairingptr; (* pointer into one structure *)
   toldversion: boolean; (* tell the version of mfold being read *)
   word: string; (*  a word in the table *)

begin
   writeln(output,'Reading in structures:');
   reset(cootab);
   if eof(cootab) then begin
      writeln(output,'cootab is empty');
      halt
   end;

   new(st);
   st^.next := nil;
   h := st;
   toldversion := false;

   while not eof(cootab) do begin (* read in all the strucures *)
      with h^ do begin
         (* read the header *)
         read(cootab,length);
{
writeln(output,'cootab length = ',length:1);
}
         skipcolumn(cootab); (* skip "ENERGY" *)
         skipblanks(cootab);
{
         read(cootab, equals); (* = *)
}
         grabtoken(cootab,word);
         if word.letters[1] <> '=' then begin
            writeln(output,'readcoordinatetable: "=" not found');
{ handle this format, probably a bug in mfold 3.1
   51   [initially    -2.8]    EG12#1
}
            writeln(output,'checking for "["');
            if word.letters[word.length] <> ']' then begin
               writeln(output, 'readcoordinatetabe: could not find "]"');
               writestring(output, word);
               writeln(output);
               halt
            end;
            energy := 0.0;
            grabtoken(cootab,idnumber);
         end
         else begin
            read(cootab,energy);
            grabtoken(cootab,idnumber);
         end;


         if idnumber.letters[1] = '[' then begin
            (* handle version 3.1 cootab format:
   27   dG = -15.1  [initially   -15.1]    mfold.demo.1#1                   
            To do this, skip to the next ] and then continue.
            *)
            grabtoken(cootab,idnumber);
            grabtoken(cootab,idnumber);
            unnumber(idnumber, id, sequnumber);
            if not toldversion then begin
               toldversion := true;
               writeln(output,'reading mfold 3.1 format')
            end;
         end
         else begin
            (* handle version 3.0 cootab format:
   51 ENERGY =    -8.1    EG12#1                        
            *)
            unnumber(idnumber, id, sequnumber);
            if not toldversion then begin
               toldversion := true;
               writeln(output,'reading mfold 3.1 format')
            end;
         end;

         readln(cootab);

         write(output,'  ');
         writequotestring(output,id);
         write(output,' length = ',length:2);
         write(output,' energy = ',energy:5:1);
         write(output,' idnumber = ');
         writestring(output,idnumber);
         write(output,' sequnumber = ',sequnumber:1);
         writeln(output);

         (* read in the pairs that define this structure *)
         new(pairs);
         c := pairs;
         c^.next := nil;
         repeat
            with c^ do begin
               read(cootab, i);
{
write(output,' i=',i:2);
}
               skipblanks(cootab);
               read(cootab, thebase); (* the base *)
               read(cootab, five); (* a number *)
               read(cootab, three); (* a number *)
               read(cootab, basepaired);
               read(cootab, hstnum);
               readln(cootab);
{
write(output,' basepaired=',basepaired:2);
}

               if eof(cootab) then begin
{
writeln(output,' eof of cootab');
}
                  if i <> h^.length then begin
                     writeln(output,'table ended at record ',i:1,
                                    ' This is not the reported length (',
                                    h^.length:1,')');
                     halt;
                  end;
               end
               else begin
{
writeln(output,' NOT eof of cootab');
}
                  if i <> h^.length then begin
                     new(c^.next);
{
writeln(output,' new link');
writeln(output,'i = ',i:1, ' h^.length =',h^.length:1);
}
                     c := c^.next;
                     c^.next := nil;
                  end
{
else writeln(output,'final link');
}
               end

            end;
         until c^.i = h^.length
      end;

      if eof(cootab) then begin
{
writeln(output,' eof of cootab, wesa done');
}
      end
      else begin
{
writeln(output,' NOT eof of cootab more structure to go');
}
{could string links together? no - separate loop needed}
         new(h^.next);
         h := h^.next;
         h^.next := nil;
      end
   end;
end;
(* end module readcoordinatetable *)

(* begin module writecoordinatetable *)
procedure writecoordinatetable(var f: text; st: structureptr);
(* write the coordinate table out from structure st *)
var
   s: structureptr; (* a pointer to st *)
   p: pairingptr; (* a pointer to s^.pairs *)
begin
   s := st;
   while s <> nil do with s^ do begin
      write(f, length:5);
      write(f, ' ENERGY = ');
      write(f, energy:7:1,'    ');
      writestring(f, id);
      writeln(f);
      p := pairs;
      while p <> nil do with p^ do begin
         write(f, i:5);
         write(f, ' ',thebase);
         write(f, ' ',five:7);
         write(f, ' ',three:4);
         write(f, ' ',basepaired:4);
         write(f, ' ',hstnum:4);
         writeln(f);
         p := p^.next
      end;
      s := s^.next
   end
end;
(* end module writecoordinatetable *)

(* begin module findlink *)
function findlink(n: integer; pa: pairingptr): pairingptr;
(* locate the link corresponding to n in the structure pa *)
var
   s: pairingptr; (* index to the pairs *)
   done: boolean; (* done with looking? *)
begin
   s := pa;
   done := false;
   while not done do begin
{
writeln(output,'s^.i=',s^.i:1);
}
      if s^.i = n then begin
         findlink := s;
         done := true;
      end
      else begin
         if s^.next = nil then begin
            writeln(output,'linkct: can''t locate pairing ', n:1);
            halt;
         end;
         s := s^.next;
      end
   end
end;
(* end module findlink *)

(* begin module linkct *)
procedure linkct(var st: structureptr);
(* make cross links in the structure st *)
var
   a, b: pairingptr; (* points in structure st *)
begin
   with st^ do begin
       a := pairs;
       b := a;
{ ... is this really needed yet? }
   end;
end;
(* end module linkct *)

(* begin module buildnames *)
procedure buildnames(var st: structureptr);
(* build the name of each structure *)
const
   dot = '.'; (* name part separater *)
   thedecimal = 1; (* decimal places for energy *)
   thewidth = 5; (* width in places for energy *)
var
   anumber: string; (* a string with a number in it *)
   helix: string; (* string with "helix" in it *)
   last: pairingptr; (* the last coordinate of s *)
   s: structureptr; (* pointer into st *)
begin
   writeln(output,'Building names of structures:');

   (* start the strings *)
   clearstring(helix);
   helix.letters[1] := 'h';
   helix.letters[2] := 'e';
   helix.letters[3] := 'l';
   helix.letters[4] := 'i';
   helix.letters[5] := 'x';
   helix.length := 5;

   s := st;
   while s <> nil do with s^ do begin
      (* the name starts with the string given.  Attach to it *)
{
      copystring(id, nametag);

      (* add to the name tag the number in the book
         so that the name becomes unique *)
      attachcharacterstring('#',nametag);
      numberstring(sequnumber, anumber); (* first coordinate *)
      attachstring(nametag,anumber);
      attachcharacterstring(dot,nametag);

      attachstring(nametag, helix); (* attach "helix" *)
}

      copystring(helix, nametag); (* start name with "helix" *)
      attachcharacterstring(dot,nametag);
      attachstring(nametag, id); (* attach id *)

      (* add to the name tag the number in the book
         so that the name becomes unique *)
      attachcharacterstring('#',nametag);
      numberstring(sequnumber, anumber); (* first coordinate *)
      attachstring(nametag,anumber);

      attachcharacterstring('[',nametag);

      numberstring(pairs^.hstnum, anumber); (* first coordinate *)
      attachstring(nametag,anumber);
      attachcharacterstring(dot,nametag);

      last := findlink(st^.length, st^.pairs);
{
writeln(output,'last^.i = ',last^.i:1);
writeln(output,'last^.hstnum = ',last^.hstnum:1);
}
      numberstring(last^.hstnum, anumber); (* last coordinate *)
      attachstring(nametag,anumber);
      attachcharacterstring(']',nametag);

      (* energy *)
      anumber.length := thewidth;
      numberrealstring(anumber,energy,1,thedecimal);

      attachstring(nametag,anumber);

      write(output,'  ');
      writequotestring(output,nametag);
      writeln(output);

      s := s^.next
   end;
end;
(* end module buildnames *)

(* begin module buildcoordinates *)
procedure buildcoordinates(var book: text; var st: structureptr);
(* Build the coordinate (hstnum) of each structure, using the information
in the book. *)
var
   apiece: pieceptr; (* pointer a piece of dna *)
   c: pairingptr; (* pointer into st pairs *)
   done: boolean; (* done going through structures for a piece? *)
   s: structureptr; (* pointer into st *)
   theline: integer; (* line number of the book *)
   thenumber: integer; (* the number of the sequence in the book *)
begin
   writeln(output,'Building coordinate systems of structures:');

   (* start the book *)
   brinit(book, theline);
   new(apiece);
   thenumber := 0;

   while (not eof(book)) do begin
      getpiece(book, theline, apiece);
      if not eof(book) then begin

         (* locate structures found on this piece *)
         s := st;
         done := false;
         thenumber := succ(thenumber);
         while not done do begin
{
writeln(output,'build loop eof(book) = ',eof(book));
if s=nil then writeln(output,'s=nil')
         else writeln(output,'s<>nil');
writequotestring(output,s^.id);
writeln(output,': try to match this to the book name');
writeln(output,apiece^.key.hea.keynam.letters:apiece^.key.hea.keynam.length);
}

            (* the names must match and the number has to match also! *)
            if ((isnamestring(apiece^.key.hea.keynam, s^.id)) and
                (thenumber = s^.sequnumber))
            then begin

{
writequotestring(output,s^.id);
writeln(output,': structure match to book name');
}

               with apiece^.key do begin
                  with hea.keynam do write(output, letters:length);
                  write(output,' ',piebeg:1);
                  write(output,' ',pieend:1);
                  case piedir of
                    minus: write(output,' ',-1:2);
                    plus : write(output,' ',+1:2);
                  end;
write(output,' ');
writestring(output,s^.idnumber);
write(output,' sequence number: ', thenumber:1);
                  writeln(output);

                  (* put the coordinate system into the structure *)
                  if s^.length <> piecelength(apiece) then begin
                     writeln(output, 'piece length in book (',
                                     piecelength(apiece):1,
                                     ') does not match',
                                     ' piece length of structure (',
                                     s^.length:1,
                                     ')');
                     halt
                  end;
                  (* go through the structure and insert the coordinates *)
                  c := s^.pairs;
                  while c <> nil do begin
                     c^.hstnum := inttopie(c^.i, apiece);
                     case piedir of
                       minus: s^.orientation := -1;
                       plus : s^.orientation := +1;
                     end;
{
write  (output,'c^.i = ',c^.i:1);
writeln(output,'  c^.hstnum = ',c^.hstnum:1);
}
                     c := c^.next
                  end;

               end;
            end;
            s := s^.next;
            if s = nil then begin
               done := true;
            end
         end;

         clearpiece(apiece); (* clear the piece for reuse *)
      end
      else if s <> nil then begin
         writeln(output,'buildnames: end of book found before end of structures');
         halt
      end;
   end
end;
(* end module buildcoordinates *)

(* begin module buildfeature *)
procedure buildfeatures(var thefile: text; st: structureptr);
(* build the feature corresponding to structure s.  Write the results to
thefile. *)
var
   c: pairingptr; (* pointer into st pairs *)
   d: definitiontypeptr; (* a definition being generated *)
   f: featuretypeptr; (* a feature being generated *)
   s: structureptr; (* pointer into st *)
begin
   writeln(output,'Writing features');

   (* set up the permanent parts of the definition *)
   new(d);
   cleardefinition(d);
   d^.background.letters[1] := '-';
   d^.background.length := 1;

   (* set up the permanent parts of the feature *)
   new(f);
   clearfeature(f);

   s := st;
   while s <> nil do begin
      copystring(s^.nametag, d^.nametag);
      c := s^.pairs;

      with d^ do begin
         clearstring(negparts);
         clearstring(posparts);
         marks := 0;
         while c <> nil do begin
            if (c^.basepaired > 0) then begin
               if c^.i < c^.basepaired
               then attachcharacterstring('(',negparts)
               else attachcharacterstring(')',negparts);
               marks := succ(marks);
               locations[marks] := c^.i-1; (* so zero base is start *)
            end;
            c := c^.next
         end;
         copystring(negparts, posparts);
      end;

      copystring(s^.id, f^.id);
      copystring(s^.nametag, f^.nametag);
      with f^ do begin
         coordinate := s^.pairs^.hstnum;
         orientation := s^.orientation;
      end;

      (* write out the results from d and f *)
      writeadefinition(thefile,d);
      writeafeature(thefile,f);
      writeln(thefile);
      writeln(thefile);

      s := s^.next;
   end;
end;
(* end module buildfeature *)

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

(* begin module mfoldfea.themain *)
procedure themain(var cootab, {mfoldfeap,} book, mfoldfeatures: text);
(* the main procedure of the program *)
var
{
   parameterversion: real; (* parameter version number *)
}
   st: structureptr; (* linked list of structures *)
begin
   writeln(output,'mfoldfea ',version:4:2);
   PatentMessage(output);
{
Not in use yet:
   reset(mfoldfeap);
   readln(mfoldfeap, parameterversion);
   if parameterversion < updateversion then begin
      writeln(output, 'You have an old parameter file!');
      halt
   end;
}

   reset(cootab);
   rewrite(mfoldfeatures);
   writeln(mfoldfeatures,'* mfoldfea ',version:4:2);

   readcoordinatetable(cootab, st);
   buildcoordinates(book, st);
{
   for testing purposes:
   writecoordinatetable(output, st);
halt;
}

   buildnames(st);
   buildfeatures(mfoldfeatures, st);

   writeln(output,'mfoldfea is done.');

end;
(* end module mfoldfea.themain *)

begin
   themain(cootab, {mfoldfeap,} book, mfoldfeatures);
1: end.










