program comp(book,cmp,compp,output);
(* determine the oligonucleotide composition of a book.
   by Gary Stormo and Tom Schneider
module libraries required: delman, delmods, auxmods *)

label 1; (* end of program *)

const
(* begin module version *)
version = 5.27; (* of comp.p 1999 Oct 13
1999 Oct 13: 5.27: added E. coli composition example.
previous changes: 1994 Sep 5
origin before 1988 oct 10 *)
(* end module version *)

(* begin module describe.comp *)
(*
name
      comp: determine the composition of a book.

synopsis
      comp(book: in, cmp: out, compp: in, output: out)

files
      book: the sequences;
      cmp: the composition, determined for mononucleotides up to
            oligonucleotides of length "compmax", see file compp;
      compp: parameter file used to set the length of the oligonucleotides for
            which the composition is to be determined ("compmax");  that number
            must be the first thing in the file; if the file is empty
            compmax is set by default to the constant "defcompmax";
      output: for messages to the user.

description
      Comp counts the number of each oligonucleotide (from length 1 to
      compmax) in the book and prints that to file "cmp".  The output is
      printed in order of increasing length of oligonucleotide (i.e., first
      the monos, then the dis, ...).  If there are no occurences of an
      oligonucleotide, but its one-shorter parent did occur, it will be given
      a zero.  None of its descendants will be printed in the composition
      file.

examples

   As an example of the output format, the composition to depth 3 of E. coli
   (U00096, 16-OCT-1997) is:

 comp 5.27: composition of 
 * 1999/05/04 14:41:13, 1999/05/04 14:38:08, dbbk 3.33
 3 is the longest oligo counted

 *
 0-long oligos (the total number of bases)
          4639221
 *
 1-long oligos
      a   1142136      c   1179433      g   1176775      t   1140877
 *
 2-long oligos
      aa   337835      ac   256658      ag   237851      at   309792
      ca   325118      cc   271649      cg   346636      ct   236029
      ga   267234      gc   383865      gg   270083      gt   255593
      ta   211948      tc   267261      tg   322205      tt   339463
 *
 3-long oligos
      aaa  108901      aac   82578      aag   63364      aat   82992
      aca   58633      acc   74899      acg   73263      act   49863
      aga   56618      agc   80848      agg   50611      agt   49774
      ata   63692      atc   86476      atg   76229      att   83395
      caa   76607      cac   66752      cag  104785      cat   76974
      cca   86442      ccc   47764      ccg   87031      cct   50412
      cga   70934      cgc  115673      cgg   86870      cgt   73159
      cta   26762      ctc   42714      ctg  102900      ctt   63653
      gaa   83490      gac   54737      gag   42460      gat   86547
      gca   96010      gcc   92961      gcg  114609      gct   80285
      gga   56199      ggc   92123      ggg   47470      ggt   74291
      gta   52670      gtc   54225      gtg   66108      gtt   82590
      taa   68837      tac   52591      tag   27241      tat   63279
      tca   84033      tcc   56025      tcg   71733      tct   55469
      tga   83483      tgc   95221      tgg   85132      tgt   58369
      tta   68824      ttc   83846      ttg   76968      ttt  109825

see also
      compan.p, histan.p, markov.p

authors
      Gary Stormo and Tom Schneider

bugs
      none known

technical note
      The algorithm is an interesting application of linked lists.  The
      composition is stored as a tree, and a number of "spiders" climb the
      tree during its construction.

*)
(* end module describe.comp *)

(* begin module comp.const *)
      defcompmax = 2;  (* default value for variable compmax *)
(* end module comp.const *)

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

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

(* end module book.const version = 'delmod 6.51 85 apr 17 tds/gds' *)

type
(* begin module book.type *)
(* types needed for book manipulations *)

      chset = set of 'a'..'z';

      (* types defined in book definition *)

      alpha = packed array[1..namelength] of char; (* this is not alfa *)

      (* name is a left justified string with blanks following the
         characters *)
      name = record
         letters: alpha;
         length: 0..namelength (* zero means an unspecified structure *)
      end;

      lineptr = ^line;
      line = record (* a line of characters *)
         letters: packed array [1..linelength] of char;
         length: 0..linelength;
         next: lineptr
      end;

      direction = (plus, minus, dircomplement, dirhomologous);
      configuration = (linear, circular);
      state = (on, off);


      header = record (* header of key *)
         keynam: name; (* key name of structure *)
         fulnam: lineptr; (* full name of structure *)
         note:   lineptr (* note key *)
      end;


      (* base types *)
      base = (a,c,g,t);
      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 = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module comp.type *)
            (* the composition is stored in a tree of these nodes *)
      compnodeptr = ^compnode; (* points to a node of the tree *)
      compnode = record    (* a node of the composition tree *)
         count : integer; (* the number of oligos for this node *)
         son : array[base] of compnodeptr; (* the pointers to 'descendants'
                                           of this node of the tree *)
      end;

            (* spiders are used to make the composition tree *)
      spiderptr = ^spider; (* points to a 'spider' *)
      spider = record  (* a spider climbs the composition tree, its path
                        determined by the sequences, and increments 'count' at
                        all the nodes it passes, thereby determining the
                        composition *)
         depth : integer; (* the level of the node now at *)
         place : compnodeptr; (* a pointer to the current node *)
         next : spiderptr; (* the next spider in the collection *)
      end;

            (* the total number of composition entries at a given level
               is stored in the linked list of type comptotal *)
      compptr = ^comptotal;
      comptotal = record
         count: integer;  (* the number at a given level *)
         next: compptr;  (* pointer to the next level totals *)
      end;

           (* 'path' is used in printing the tree *)
      pathptr = ^path; (* pointer into a path on the tree *)
      path = record  (* the path of bases to get to a particular node *)
         bas : base;
         next : pathptr;
      end;
(* end module comp.type version = 'auxmod 1.37 85 apr 4 gds/tds'; *)


var
      book, cmp, compp: text;
      compmax : integer;  (* the longest oligonucleotide for which the
                                 composition is determined *)
      pie : pieceptr; (* the piece currently being used *)
      root : compnodeptr;  (* the 'root' of the composition tree *)
      firstspider : spiderptr; (* the first in the string of spiders *)

(* 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 = 'delmod 6.51 85 apr 17 tds/gds' *)

(* 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 = 'auxmod 1.37 85 apr 4 gds/tds'; *)

(* begin module copyaline *)
procedure copyaline(var fin, fout: text);
(* copy a line from file fin to file fout *)
begin (* copyaline *)
      while not eoln(fin) do begin
         fout^ := fin^;
         put(fout);
         get(fin)
      end;
      readln(fin);
      writeln(fout);
end; (* copyaline *)
(* end module copyaline version = 'delmod 6.51 85 apr 17 tds/gds' *)

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

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

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

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

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

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

procedure 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 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 = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.getto *)
function getto(var thefile: text; ch: chset): char;
(* search the file for a character in the first line which is a
   member of the set ch. *)
var achar: char;
begin
      achar:=' ';
      while (not(achar in ch)) and (not eof(thefile))
         do readln(thefile,achar);
      if (achar in ch) then getto:=achar
                       else getto:=' '
end;
(* end module book.getto version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.skipstar *)
procedure skipstar(var thefile:text);
(* skip start of line (or star = '*'). *)
begin (* skipstar *)
      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; (* skipstar *)
(* end module book.skipstar version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brreanum *)
procedure brreanum(var thefile: text; var reanum: real);
(* read a real number from the file *)
begin
      skipstar(thefile);
      readln(thefile,reanum);
end;
(* end module book.brreanum version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brnumber *)
procedure brnumber(var thefile: text; var num: integer);
(* read a number from the file *)
begin
      skipstar(thefile);
      readln(thefile,num)
end;
(* end module book.brnumber version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brname *)
procedure brname(var thefile: text; 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)
end; (* brname *)
(* end module book.brname version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brline *)
procedure brline(var thefile: text; var l: lineptr);
(* read a line from the file *)
var
      i,j: 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;
      if i<l^.length then for j:=i+1 to l^.length do l^.letters[j]:=' ';
      l^.length:=i;
      l^.next:=nil;
      readln(thefile)
end;
(* end module book.brline version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brdirect *)
procedure brdirect(var thefile: text; var direct: direction);
(* read a direction *)
var   ch: char;
begin
      skipstar(thefile);
      readln(thefile,ch);
      if ch='+' then direct:=plus
                else direct:=minus
end;
(* end module book.brdirect version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brconfig *)
procedure brconfig(var thefile: text; var config: configuration);
(* read a configuration *)
var   ch: char;
begin
      skipstar(thefile);
      readln(thefile,ch);
      if ch='l' then config:=linear
                else config:=circular
end;
(* end module book.brconfig version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brnotenumber *)
procedure brnotenumber(var thefile: text; 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);
         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)
            until thefile^ = 'n';
            readln(thefile)
         end
         else readln(thefile)
      end
end; (* brnotenumber *)
(* end module book.brnotenumber version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brnote *)
procedure brnote(var thefile: text; 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);
         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,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)
         end
         else readln(thefile)
      end
end; (* brnote *)
(* end module book.brnote version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brheader *)
procedure brheader(var thefile: text; var hea: header);
(* read the header of a key. *)
begin
      with hea do begin
         (* read key name *)
         brname(thefile,keynam);

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

         (* read note key *)
         if readnumber then brnotenumber(thefile,note)
                       else brnote(thefile,note)
      end
end;
(* end module book.brheader version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brpiekey *)
procedure brpiekey(var thefile: text; var pie: piekey);
(* read piece key *)
begin
      with pie do begin
         brheader(thefile,hea);
         brreanum(thefile,mapbeg);
         brconfig(thefile,coocon);
         brdirect(thefile,coodir);
         brnumber(thefile,coobeg);
         brnumber(thefile,cooend);
         brconfig(thefile,piecon);
         brdirect(thefile,piedir);
         brnumber(thefile,piebeg);
         brnumber(thefile,pieend);
      end
end;
(* end module book.brpiekey version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brdna *)
procedure brdna(var thefile: text; var dna: dnaptr);
(* read in dna from thefile *)
(* 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,['d']);
      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 *)
         read(thefile,ch); (* ch is either '*' or 'd' *)
      end;
      readln(thefile)
end;
(* end module book.brdna version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brpiece *)
procedure brpiece(var thefile: text; var pie: pieceptr);
(* read in a piece *)
begin
      brpiekey(thefile,pie^.key);
      if numbered or (not skipunnum)
         then brdna(thefile,pie^.dna)
end;
(* end module book.brpiece version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brinit *)
procedure brinit(var book: text);
(* 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;
end; (* brinit *)
(* end module book.brinit version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* ************************************************************************ *)
(* end module package.brpiece version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.getpiece *)
procedure getpiece(var thefile: text; var pie: pieceptr);
(* move to and read in the next piece in the book *)
var   ch: char;
begin
      ch:=getto(thefile,['p']); (* get to the next p(iece) in the book *)
      if ch<>' ' then begin
         brpiece(thefile,pie);
         ch:=getto(thefile,['p']); (* read past closing p *)
      end
end;
(* end module book.getpiece version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* ************************************************************************ *)
(* end module package.getpiece version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.stepbase *)
function stepbase(startdna: dnaptr; var dna: dnaptr; var d: dnarange): base;
(* advance d by one base in dna and then return the base at the new d.
   (this means that one should initialize d to zero)
   if we go past the last base, we restart at startdna.
   note: d is not the number of the base... it is used as a
   record for stepbase.  do not mess with it, and do not use it to find
   out what base you are on.  use a separate counter.  *)
begin
      if (d=dnamax) or (d=dna^.length) then begin
         d:=1;
         dna:=dna^.next;
         if dna=nil then dna:=startdna
      end
      else d:=succ(d);
      stepbase:=dna^.part[d]
end;
(* end module book.stepbase version = 'delmod 6.51 85 apr 17 tds/gds' *)

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

procedure makespiders;  (* make as many spiders as there are levels to the
                        composition tree, which is 'compmax', and leave
                        firstspider pointing to the beginning of the string *)
var i: integer;  (* an index of the 'for' loop *)
    aspider : spiderptr;
begin
      new(firstspider);
      firstspider^.depth := 0;
      firstspider^.place := root;
      firstspider^.next := nil;
      for i := 1 to (compmax-1) do begin
         new(aspider);
         aspider^.next := firstspider;
         firstspider := aspider;
      end;
end;

procedure initialize;
var ba : base;
begin (* initialize *)
      writeln(output,'comp ',version:4:2);
      brinit(book);
      rewrite(cmp);
      writeln(cmp,' comp ',version:4:2,': composition of ');
      write(cmp,' ');
      copyaline(book,cmp);

         (* get the parameter 'compmax' from compp or set to default *)
      reset(compp);
      if not eof(compp)
         then readln(compp,compmax)
         else compmax := defcompmax;
      if compmax <= 0 then begin
         writeln(output,' composition maximum must be positive');
         halt
      end;
      writeln(cmp,' ',compmax:1,' is the longest oligo counted');
      writeln(cmp);

         (* build the root of the tree *)
      new(root);
      root^.count := 0;
      for ba := a to t do root^.son[ba] := nil;

      makespiders;
      new(pie);
end;

procedure increment(var s: spiderptr; var ba: base);
(* the 'ba'-son of the node pointed to by the spider s has its count
   incremented, or is created if it does not exist *)
var bai : base;
begin
      with s^ do begin
         if (depth = 0) then place^.count := succ(place^.count);
         depth := succ(depth);
         if (place^.son[ba] = nil) then begin
            new(place^.son[ba]);
            place := place^.son[ba];
            with place^ do begin
               count := 1;
               for bai := a to t do son[bai] := nil;
            end
         end
         else begin
            place := place^.son[ba];
            place^.count := succ(place^.count);
         end;
         if (depth = compmax) then begin
            depth := 0;
            place := root;
         end;
      end;
end;

procedure maketree(pie: pieceptr; root: compnodeptr; firstspider: spiderptr);
(* increment the counts of the nodes of the tree (starting at the root)
   according to the sequence in piece pie, using the spiders *)
var active : spiderptr;
    i,j,   (* indices for loops *)
    length : integer; (* length of the sequence *)
    pos : dnarange;  (* position of previous base in the sequence *)
    thedna : dnaptr;  (* the internal current dna string *)
    ba : base;
begin
      pos := 0;
      thedna := pie^.dna;
      active := firstspider;
         (* point all spiders to root *)
      for i := 1 to compmax do begin
         active^.depth := 0;
         active^.place := root;
         active := active^.next;
      end;
      length := piecelength(pie);
      if (length >= compmax) then begin (* this is the usual case *)
            (* get started *)
         for i := 1 to compmax do begin
            active := firstspider;
            ba := stepbase(pie^.dna,thedna,pos);
            for j := 1 to i do begin
               increment(active,ba);
               active := active^.next;
            end;
         end;
            (* keep rolling through the sequence *)
         for i := compmax + 1 to length do begin
            active := firstspider;
            ba := stepbase(pie^.dna,thedna,pos);
            for j := 1 to compmax do begin
               increment(active,ba);
               active := active^.next;
            end;
         end;
         if (pie^.key.piecon = circular) then
               (* handle the case of a circular sequence *)
            for i := (length + 1) to (length + compmax -1) do begin
               active := firstspider;
               ba := stepbase(pie^.dna,thedna,pos);
               for j := 1 to compmax do begin
                  if (active^.depth >= (i - length)) then increment(active,ba);
                  active := active^.next
               end
            end
      end
      else begin   (* this is the unusual case of compmax > piecelength *)
         for i := 1 to length do begin
            active := firstspider;
            ba := stepbase(pie^.dna,thedna,pos);
            for j := 1 to i do begin
               increment(active,ba);
               active := active^.next;
            end;
         end;
         if (pie^.key.piecon = circular) then
            for i := length + 1 to compmax do begin
               active := firstspider;
               ba := stepbase(pie^.dna,thedna,pos);
               for j := 1 to compmax do begin
                  if (active^.depth >= (i - length)) then increment(active,ba);
                  active := active^.next;
               end;
            end;
      end;
end;

(* begin module comp.writecomp *)

procedure writetree(var thefile: text; root : compnodeptr);
(* start at the root of the tree and write out the paths to all the nodes
   and all the node counts.  write a zero to the empty nodes that occur
   before the depth of compmax *)
var start,  (* of the paths through the tree *)
    point : pathptr;  (* some point on the tree path *)
    printdepth,  (* the depth of the tree currently being printed *)
    i : integer;


procedure advancetonode(depth: integer; thenode : compnodeptr; point : pathptr);
(* this procedure recursively calls itself until it gets to the depth that
   is currently being printed, and writes the count at that node to the
   composition file.  if a nil branch is found at the printing level a zero
   is the output.  *)
var ba : base;


procedure writecomp(point : pathptr; count : integer);
(* write the path to the node, recorded in the string beginning with start *)
var i : integer;
begin
      write(thefile,' ':6);
            (* write the oligo for this node *)
      write(thefile,basetochar(point^.bas));
      for i := 2 to depth do begin
         point := point^.next;
         write(thefile,basetochar(point^.bas));
      end;
      write(thefile,' ',count:10-depth);
      if (point^.bas = t) then writeln(thefile);
end;

begin (* advancetonode *)
      if (depth = printdepth) then
         if (thenode <> nil) then writecomp(start,thenode^.count)
         else writecomp(start,0)
      else if (thenode <> nil) then
         for ba := a to t do begin
            point^.bas := ba;
            case ba of
               a : advancetonode(depth+1,thenode^.son[a],point^.next);
               c : advancetonode(depth+1,thenode^.son[c],point^.next);
               g : advancetonode(depth+1,thenode^.son[g],point^.next);
               t : advancetonode(depth+1,thenode^.son[t],point^.next);
            end;
         end;
end;

begin (* writetree *)
         (* make the string of path pointers *)
      new(start);
      start^.next := nil;
      for i := 1 to compmax -1 do begin
         new(point);
         point^.next := start;
         start := point;
      end;
         (* write the total number of bases *)
      writeln(thefile,' *');
      writeln(thefile,' 0-long oligos (the total number of bases)');
      writeln(thefile,' ',root^.count:16);
         (* print out each node of the tree, breadth first *)
      for printdepth := 1 to compmax do begin
         writeln(thefile,' *');
         writeln(thefile,' ',printdepth:1,'-long oligos');
         advancetonode(0,root,start);
      end;
end;
(* end module comp.writecomp *)

begin  (* comp *)
      initialize;
      while not eof(book) do begin
         getpiece(book,pie);
         if not eof(book) then begin
            maketree(pie,root,firstspider);
            clearpiece(pie);
         end;
      end;
      writetree(cmp,root);
1: end.