program sebo(inst, book, target, sebop, instout, list, output);
(* sebo: search for aligned book sequences in another book

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

module libraries:  delmod
 *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 1.47; (* of sebo.p 2010 Oct 28
2010 Oct 28: 1.47: give name instruction outside comment!
2009 Mar 03: 1.46: correct bug from Blake
2008 Jan 16: 1.45: add space before 'times in the target'
2006 Oct 26: 1.44: correct reading params
2006 Oct 26: 1.43: implement rangecontrol
2006 Oct 26: 1.42: allow output instructions to match original fragment
2004 Sep  9: 1.41: upgrade for GPC
2004 Apr 30: 1.40: spelling correction, guarantee
2003 Aug 25: 1.39: bug: output at every base for book alignment
2002 Aug 27: 1.38: Do not halt if probe is larger than target sequence
2001 May  7: 1.37: report total number of probes found.
2001 Mar 26: 1.34: note in bugs on shift parameter
2001 Mar 23: 1.33: count targets and probes
2001 Mar 23: 1.32: clean up after memory leak solved.
2001 Mar 23: 1.31: hunt for a memory leak
2000 Nov  8: 1.25: upgrade documentation and fix writeinstructions bflip
2000 jun 27: 1.20: names of instout default to previous names
2000 jun 27: 1.18: add examples
2000 jun 27: 1.17: fully functional
2000 jun 26: 1.14: alert parameter
2000 jun 26: 1.12: allow mismatches
2000 jun 26: 1.11: fully functional: found that Ecoli ahpC is next to dsbG!
2000 jun 26: 1.04: barely functional
2000 jun 24: 1.01: reads target and book.
2000 jun 24: 1.00: origin *)
updateversion = 1.21; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.sebo *)
(*
name
   sebo: search for aligned book sequences in another book

synopsis
   sebo(inst: in, book: in, target: in, sebop: in,
        instout: out, list: out, output: out)

files

   inst:  Delila instructions to create a book.

   book:  The Delila book created with inst, containing fragments to search
          in target.

   target:  Target book to search, usually a whole genome.  The sequence is
      treated as a circle to guarantee correct results on circular genomes.
      The incorrect results on linear chromosomes will be quite rare because
      of the telomeres.

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

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

      2. If the first character of the second line is
         'f' (for 'first') then the sequences are always aligned by their
         first base.

         'i' then the sequences are aligned by the delila instructions.  If
         the inst file is empty, alignment is forced to the 'b' mode.

         'b' (for 'internal') then the alignment is on the internal zero of
         the book's sequence.  This option is to be used when "default
         coordinate zero" is used in the Delila instructions.

      3. fromdo, todo: (integer) range of aligned book to use in the search.


      4. rangecontrol, frominst, toinst: (char, integer, integer):
         Range of Delila instructions output to instout.

         If range control, the first character on the line is:

         'a' then the instructions will give the ORIGINAL fragment in
         absolute coordinates.

         'r' then the instructions will give the ORIGINAL fragment in
         relative coordinates.

         For any other initial character, the instructions will give
         the range frominst to toinst.

      5. mismatches: (integer) number of mismatches allowed.

      6. alertinterval: (integer) Many genomes are large and a search may
         take a while, making the program seem to be doing nothing.  To
         observe that the program is running, set this parameter to a number
         of bases.  Sebo will write out the location it is currently
         searching every alertinterval bases.  Nonpositive values disable the
         feature.

   instout:  revised inst giving Delila instructions for finding the
      fragments specified in inst/book in the target.  Full names of pieces
      from the book (such as those assigned by Delila instructions of the
      form "name 'abinddingsite';") will default to previous full names if
      there are any.  This follows the convention of blank name on the alist
      being the previous name.

   list: display of aligned locations found.

   output: messages to the user

description

   Nowadays it is often the case that old GenBank entries become fused into a
   single genome.  It would be nice to be able to automatically convert old
   delila instructions to the new genome.

examples

********************************************************************************
                              Example 1

----- sebo.inst:

title "An OxyR site";
{Example instruction file for sebo }
{ number only the pieces, starting at 1 }
  default numbering piece;
  default numbering 1;
  default out-of-range reduce-range; 
organism E.coli;
chromosome E.coli;
name "ahpC";
piece D13187; { ECOD13187 };
get from 116 -200 to 116 +200 direction +;

----- sebop:

1.21       version of sebo that this parameter file is designed for.
i          f: first base, i: inst, b: book alignment
-30 30     fromdo todo: range of aligned book to use in the search.
-200 200   frominst toinst: range of delila instructions output to instout
2          number of mismatches allowed
500000     alertinterval (bases)

----- target:

E. coli genome U00096

----- instout:

title "sebo 1.21 search with book";
{ sebop parameters:
i: use the alignedbase from the book
       -30         30 range of book alignment used for search
      -200        200 range of resulting output instructions
         2            mismatches allowed
    500000            alert interval
}
organism E.coli; chromosome E.coli;
piece U00096; { the target sequence is 4639221 bases }

{ original source probe: }
{ name ""; piece D13187; get from 1 to 316 direction +; }
{ final target instructions: }
name "ahpC"; get from 638089 -200 to same +200 direction +;

----- list:

sebo 1.21 search with book";
sebop parameters:
i: use the alignedbase from the book
       -30         30 range of book alignment used for search
      -200        200 range of resulting output instructions
         2            mismatches allowed
    500000            alert interval

Target sequence: 
piece U00096; get from 1 to 4639221 direction +; { 4639221 bases }
   Searching target with probe: name "ahpC"; piece D13187; get from 1 to 316 direction +;

   ---------------------                   +++++++++++++++++++++
   322222222221111111111--------- +++++++++111111111122222222223
   0987654321098765432109876543210123456789012345678901234567890
   ttacgaaggttgtaaggtaaaacttatcgatttgataatggaaacgcattaccggaatcgg probe
                                                      X X        differences
   ttacgaaggttgtaaggtaaaacttatcgatttgataatggaaacgcattagccgaatcgg target
                                 ^ found at 638089 on the target

********************************************************************************
                              Example 2

Suppose that you have some sequences in a file and one would like to locate
them on the chromosome of an organism in the form of Delila instructions.
Here is how to create the input files for sebo:

1.  Convert the sequences into a Delila book using one of the programs
rawbk, makebk, dbbk, or mkdb (see the web page mentioned below).  This gives
you the file 'book' in directory 1.

2.  In another directory (2) build a Delila library for your organism (this
is also described on that web page).

3.  Use delila to extract both the sequence and the complement, for example:

   title "Delila instructions for Bacteriophage T4 and its complement";
   organism c.T4; chromosome c.T4; piece AF158101;
   get all piece;
   get all piece direction -;

4. Move or link the resulting book into the file 'target' in Directory 1.

5. Then set up the sebop for searching without Delila instructions (since no
delila instructions were used to make the book):

1.21       version of sebo that this parameter file is designed for.
f          f: first base, i: inst, b: book alignment
0  40      fromdo todo: range of aligned book to use in the search.
0  40      frominst toinst: range of delila instructions output to instout
10         number of mismatches allowed
500000     alertinterval (bases)

6. Finally use sebo to create the instout.  You may have to allow for some
mismatches.  These instructions should allow you to extract the sequences
from the organism database

7. Sebo reports its results by adding lots of comments to instout in the form
{}.  If you want to remove them, (so you have a clean set of instructions)
you can use the nocom program.  You can remove the extra blanks with
noblank.  In unix:

nocom < instout | sed -e 's/{}//g' | noblank > instclean

You can also remove redundant delila instructions for organism, chromosome
and piece by hand, especially if there is only one piece in the target.

********************************************************************************

documentation


see also

   {Sequence extraction program:} delila.p
   {Simple search engine:} search.p
   {More details on alignment types:} alist.p

   {example sebop parameter file:} sebop
   {example inst file:} sebo.inst

   {Programs for conversion of raw sequences into a Delila Book:}
   rawbk.p makebk.p dbbk.p mkdb.p

   {Conversion of a book of sequences into a Delila Library:}
   https://alum.mit.edu/www/toms/delilalibraries.html

   {Remove the comments created by sebo:} nocom.p
   {Remove the blank lines created by sebo:} noblank.p

   {Shift the instructions to whereever one wants:} instshift.p

author

   Thomas Dana Schneider

bugs

technical notes

*  The org/chr/pie names written to delila instructions come from the current
   target sequence.

*  Delila instructions for sites frequently come in pairs:

get from 163 -200 to 163 +200 direction +;
get from 163 +200 to 163 -200 direction -;

   So the target will be searched in both orientations.

*  constant dnamax is set to a large value to increase
   the speed of the search.

* a nice parameter to have would allow the user to shift the coordinate of
the zero into any position they want in the piece.  Fortunately this can be
done with the instshift program, so this is not high priority.  The base that
one specifies would become the zero base of the Delila instructions.

* 2001 April 17: if the fromdo/todo range given for the book is outside of
the piece in the book (in aligned searching), then the book sequence will not
be found in the target.

*)
(* end module describe.sebo *)

(* 2001 Mar 25
* memory leak information:
The memory leak can be seen with the SIZE parameter of the top program.
The value of SIZE continuously increases as the program runs.
Memory leak occurs under gpc and sun compilers.
The freedna variable always has 0 links during the loop!
The cleardna routine is never called, no halt there.
The clearpiece routine IS called, halt there halts.
The call to clearpiece is from getocp.
Why are no links going onto the freedna list at that point?

I have demonstrated that there is indeed a memory leak!  By cranking up the
size of the dna segments from 5000000 (5e6) to 100000000 (1e8) I got sebo to
find a bunch of targets and THEN to crash:

Ran out of memory
Trace/Breakpoint Trap

Since it should reuse the arrays, this should not happen unless the arrays
are being wasted (memory is allocated but lost and then the next sequence
uses more memory).

Also, top shows it going up more rapidly.

      dnamax = 100000000; { causes crash}

Now to find out why it is loosing the dna links and not reusing them.

It was rather subtle.  I was looking at a routine called getocp (get
organism, chromosome, piece) which DOES clear the DNA with each use.  To my
surprise sebo was not using getocp.  But the program can also use alignments,
depending on the sebop.  Indeed, you had set the sebop to use delila
instructions, but had left the inst file EMPTY.  This caused the program to
use the alignement routines to read the sequences.  But that routine does not
clear the DNA and so there was an overflow.

So I have put a halt in; 
 
sebo 1.31 
i: use the alignedbase from the book 
        75        125 range of book alignment used for search 
      -100        100 range of resulting output instructions 
         2            mismatches allowed 
     50000            alert interval 
The sebop requests instructions to be used 
but the inst file is empty! 
 program halt. 

However these routines are called:

  align
  getpiece
  brpiece
  brdna
  getdna

So getdna is probably guilty!  Here is the original code:

      if freedna<>nil
         then begin 
IT NEVER REACHES HERE
            l:=freedna;
            freedna:=freedna^.next
         end
         else new(l);
      l^.length:=0;
      l^.next:=nil

Fortunately the bug was simple.  I used getocp which automatically takes care
of memory.  But for the align routine, it can't do that and so the sebo
program has to.  I forgot to clearpiece(Bpie); Adding this solved the
problem.

*)

(* LOCK begin module book.const *)
(* constants needed for book manipulations *)
{
      dnamax = 1024;  (* length of dna arrays *)
}
      dnamax = 5000000;  (* length of dna arrays
                            INCREASE TO GENOME SIZE FOR SPEED *)

      namelength = 100; (* maximum key name length *)
      linelength = 80; (* maximum line readable in book *)
(* LOCK end module book.const version = 7.49; {of delmod.p 2000 June 26} *)

(* begin module string.const *)
      maxstring = 2000; (* the maximum string *)
(* end module string.const version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* begin module filler.const *)
      fillermax = 50; (* the size of the filler array for a string *)
(* end module filler.const version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

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

(* begin module filler.type *)
      (* the following is an array used to fill a string.
      it is convenient to have it much shorter than the maxstring, so that
      it is easy to fill the string using procedure fillstring.
      the user must declare the value of constant fillermax. *)
      filler = packed array[1..fillermax] of char;
(* end module filler.type version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* begin module trigger.type *)
      trigger = record (* an object to be searched for *)
         seek: string; (* the characters looked for *)
         state: integer; (* how close to triggering we are *)
         skip: boolean; (* trigger not found- skip the line *)
         found: boolean (* the trigger was found *)
      end;
(* end module trigger.type version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

      (* 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.67; {of delmod.p 2004 Sep 8} *)

var
   inst, book, (* file used by this program *)
   target, (* file used by this program *)
   sebop, (* file used by this program *)
   list, (* file used by this program *)
   instout: 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.67; {of delmod.p 2004 Sep 8} *)

(* 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 = 7.67; {of delmod.p 2004 Sep 8} *)

(* 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 = 7.67; {of delmod.p 2004 Sep 8} *)

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

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

(* begin module clearstring *)
procedure clearstring(var ribbon: string);
(* empty the string *)
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; (* clearstring *)
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. *)
begin (* initializestring *)
   clearstring(ribbon);
   ribbon.next := nil;
end; (* initializestring *)
(* end module clearstring version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* 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 = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* begin module filler.fillstring *)
procedure fillstring(var s: string; a: filler);
(* this procedure makes it reasonably easy to fill the string s with
characters.  one calls the procedure as: *)
(*                           1         2         3         4         5 *)
(*                  12345678901234567890123456789012345678901234567890 *)
(*   fillstring(s, 'this-is-the-string                                ');
the two comments make it easy to line the characters up. also, for this
example, it was assumed that the length of filler as defined by the
constant fillermax was 50. *)
var
      length: integer; (* of the string without trailing blanks *)
      index: integer; (* of s *)
begin (* fillstring *)
      clearstring(s);
      length := fillermax;
      while (length > 1) and (a[length] = ' ') do length := pred(length);
      if (length = 1) and (a[length] = ' ') then begin
         writeln(output, 'fillstring: the string is empty');
         halt
      end;

      for index := 1 to length do s.letters[index] := a[index];
      s.length := length;
      s.current := 1
end; (* fillstring *)
(* end module filler.fillstring version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

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

(* begin module trigger.proc *)
(* this module allows one to scan a series of characters, as from
an array or a file, and to "trigger" or detect a simple string
in the series.  the advantage of the trigger is that several triggers
can "observe" a stream of characters at once, each looking for a
different thing.
some other modules required: interact.const, interact.type *)

procedure resettrigger(var t: trigger);
(* reset the trigger to ground state *)
begin (* resettrigger *)
      with t do begin
         state := 0;
         skip := false;
         found := false
      end
end; (* resettrigger *)

procedure testfortrigger(ch: char; var t: trigger);
(* look at the character ch.
   if it is part of the trigger (at the current trigger state),
       then the trigger state goes higher.
   if it is not part of the trigger then the trigger state is reset,
      skip is true and one should skip onward to find the trigger.
   if the trigger is found, found is true.
1996 Sep 12: Bug found!  In the case of a trigger "ab", the program
used to miss it for situations like "aab".  This was because at the
first a it would step up.  Then it would see the second a and recognize
that was not part of ab.  It would fail to realize that it could be
the start of a new one.  The code now accounts for that possibility. *)
begin (* testfortrigger *)
      with t do begin
         state := succ(state);
{
writestring(list,seek);
writeln(list,'testfortrigger seek.letters[',state:1,']:',
seek.letters[state],' ch:',ch);
}
         if seek.letters[state] = ch
         then begin
            skip := false;
            if state = seek.length then found := true
                                   else found := false
         end
         else begin
            (* it failed.  But wait!  It could be the beginning
            of a NEW trigger string! *)
            if seek.letters[1] = ch then begin
               state := 1;
               skip := false;
               found := false
            end
            else begin (* reset trigger *)
               state := 0;
               skip := true;
               found := false
            end
         end
      end
end; (* testfortrigger *)
(* end module trigger.proc version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* begin module skipblanks *)

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

function isblank(c: char): boolean;
(* is the character c blank or tab? *)
begin
   isblank := (c = ' ') or (c = chr(9))
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 = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* ************************************************************************ *)
(* end module package.trigger version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

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

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

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

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

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

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

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

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

procedure cleardna(var l: dnaptr);
(* clear the dna strutures to the free list *)
var   lptr: dnaptr;
begin
      if l<>nil then begin
         lptr:=l;
         l:=l^.next;
         lptr^.next:=freedna;
         freedna:=lptr
      end
end;

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

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

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

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

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

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

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

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

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

(* end module book.basis version = 7.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* ************************************************************************ *)
(* end module package.brpiece version = 7.67; {of delmod.p 2004 Sep 8} *)

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* ************************************************************************ *)
(* end module package.getpiece version = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module findblank *)
procedure findblank(var afile: text);
(* read a file to find the next blank character *)
var   ch: char;
begin
      repeat read(afile,ch) until ch = ' '
end;
(* end module findblank version = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module findnonblank *)
procedure findnonblank(var afile: text; var ch: char);
(* find the next non blank character in a file, return it in ch. *)
begin
      ch:=' ';
      while (not eof(afile)) and
            (ch = ' ')
      do begin
         read(afile,ch);
         if eoln(afile) then readln(afile)
      end
end;
(* end module findnonblank version = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module align.align *)
procedure align(var inst, book: text; var theline: integer;
                var pie:  pieceptr;
                var length, alignedbase: integer);
(* documentation on align is in module info.align and
delman.use.aligned.books.
1996 Sep 12: The routine now uses the trigger functions found
in prgmod.  The bug in the oldalign routine (that it misses the
end of comments that end in a series of astrisks) has been fixed.
It now checks that the piece corresponds to the book. *)
const maximumrange = 10000; (* if the alignment point is more than this
   distance from the piece ends, the program halts in an attempt to catch
   the alignment bug... 1991 Jan 11 It appears that the rewrite of the
   code has removed the bug, but the check will be kept. *)
   semicolon = ';'; (* end of delila instruction *)
var
   ch: char; (* a character in inst *)
   p: integer; (* index to a piece name *)
   p1: integer; (* another index to a piece name *)
   done: boolean; (* done finding an aligning get *)
   thebase: integer; (* the base read in *)
   indefault: boolean; (* true when within a default statement.  These
      can contain the word 'piece', which must be ignored. *)

   gettrigger: trigger; (* trigger to find 'get' *)
   defaulttrigger: trigger; (* trigger to find 'default' *)
   nametrigger: trigger; (* trigger to find 'name' *)
   piecetrigger: trigger; (* trigger to find 'piece' *)
   settrigger: trigger; (* trigger to find 'set' *)

   begincomment: trigger; (* trigger to find '(-*' (ignore the dash!) *)
   endcomment: trigger; (* trigger to find '*-)'  (ignore the dash!) *)
   begincurly: trigger; (* trigger to find comments: '{' *)
   endcurly: trigger; (* trigger to find comments: '}' *)

   quote1trigger: trigger; (* trigger to find single quote ' *)
   quote2trigger: trigger; (* trigger to find double quote " *)

{
procedure rd(var f: text; var ch: char);
(* read ch from f allowing inspection of the result *)
begin
   read(f,ch);
   write(output,ch);
   write(list,ch);
   write(output,'<',ch,'>');
end;

procedure rdln(var f: text);
(* readln f allowing inspection of the result *)
begin
   readln(f);
   writeln(output);
   writeln(list);
end;
}

procedure skipcomment(var f: text);
(* skip an entire comment *)
var
   comment: boolean; (* true means we are inside a comment *)
begin
   (* skip to end of comment *)
   resettrigger(endcomment);
   comment := true;
   while comment do begin
      if eof(f) then begin
         writeln(output,'A comment does not end!');
         halt
      end;
      if eoln(f) 
      then readln(f)
{ rdln(f) }
      else begin
{write(output,'<');
rd(f,ch);
write(output,'>');}
         read(f,ch);
         testfortrigger(ch, endcomment);
         if endcomment.found
         then comment := false;
      end
   end
end;

procedure skipcurly(var f: text);
(* skip an entire comment made by {}*)
var
   comment: boolean; (* true means we are inside a comment *)
begin
   (* skip to end of comment *)
   resettrigger(endcurly);
   comment := true;
   while comment do begin
      if eof(f) then begin
         writeln(output,'A comment does not end!');
         halt
      end;
      if eoln(f) 
      then readln(f)
{ rdln(f) }
      else begin
{write(output,'<');
rd(f,ch);
write(output,'>');}
         read(f,ch);
         testfortrigger(ch, endcurly);
         if endcurly.found
         then comment := false;
      end
   end
end;

procedure skipquote(quote: trigger);
(* skip an entire quote of either the ' or " persuasion *)
var
   kind: char; (* the kind of quote, ' or " *)
begin
   kind := quote.seek.letters[1];
{writeln(output,'skipquote ',kind);}
   repeat
      findnonblank(inst,ch); (* get to the quote *)
   until (ch = kind) or eof(inst);
   if ch <> kind then begin
      writeln(output,'end of quote starting with ',kind,' not found');
      halt;
   end;
end;

begin
 filltrigger(defaulttrigger,'default');
     filltrigger(gettrigger,'get    ');
    filltrigger(nametrigger,'name   ');
   filltrigger(piecetrigger,'piece  ');
     filltrigger(settrigger,'set    ');
   filltrigger(begincomment,'(*     ');
     filltrigger(endcomment,'*)     ');
     filltrigger(begincurly,'{      ');
       filltrigger(endcurly,'}      ');
  filltrigger(quote1trigger,'''      ');
  filltrigger(quote2trigger,'"      ');

  resettrigger(defaulttrigger);
  resettrigger(gettrigger);
  resettrigger(nametrigger);
  resettrigger(piecetrigger);
  resettrigger(settrigger);

  resettrigger(begincomment);
  resettrigger(begincurly);
  resettrigger(quote1trigger);
  resettrigger(quote2trigger);

   indefault := false;

   if not eof(book) then begin (* if there is still more to the book ... *)
      getpiece(book,theline,pie); (* read in the piece *)
      if not eof(book) then begin (* if we found a piece ... *)
         length:=pietoint(pie^.key.pieend,pie); (* calculate piece length *)

         (* now find in inst the next occurance of 'get' *)
         done := false;
         while not done do begin
            if eof(inst) then begin (* no instructions? *)
               alignedbase := 1; (* simply align by the first base *)
               done := true
            end
            else begin
               if eoln(inst)
               then readln(inst)
{then rdln(inst)}
               else begin
{rd(inst,ch);}
                  read(inst,ch);

                  testfortrigger(ch, begincomment);
                  testfortrigger(ch, begincurly);
                  if begincomment.found or begincurly.found
                  then begin
                     if ch = '*' then begin
                        skipcomment(inst);
                        resettrigger(begincomment);
                     end
                     else begin
                        resettrigger(begincurly);
                        skipcurly(inst);
                     end
                  end
                  else begin (* we are not inside a comment *)

                     testfortrigger(ch, gettrigger);
                     if gettrigger.found then begin
                        findnonblank(inst,ch); (* get to "from" *)
                        findblank(inst); (* get past "from" *)
                        read(inst,thebase); (* read in the alignedbase *)
{writeln(output);writeln(output,'thebase = ',thebase:1);}
                        alignedbase:=pietoint(thebase,pie);
{writeln(output,'alignedbase=',alignedbase:1);}
                        done := true
                     end;

                     testfortrigger(ch, quote1trigger);
                     if quote1trigger.found then begin
                        skipquote(quote1trigger);
                     end;

                     testfortrigger(ch, quote2trigger);
                     if quote2trigger.found then begin
                        skipquote(quote2trigger);
                     end;

                     testfortrigger(ch, defaulttrigger);
                     if defaulttrigger.found then begin
                        indefault := true;
                        resettrigger(defaulttrigger)
                     end;
                     if ch = semicolon then indefault := false;

                     testfortrigger(ch, settrigger);
                     if settrigger.found then begin
                        indefault := true;
                        resettrigger(settrigger)
                     end;
                     if ch = semicolon then indefault := false;

                     (* check that piece names are correct *)
                     testfortrigger(ch, piecetrigger);
                     if not indefault then if piecetrigger.found then begin
                        skipblanks(inst); (* get to name *)
                        with pie^.key.hea.keynam do begin
                           for p := 1 to length do begin
                              read(inst,ch);
                              if letters[p] <> ch then begin
                                 writeln(output,'The piece name in the book: ');
                                 writeln(output,letters:length);
                                 writeln(output,'does not match',
                                              ' the inst file piece name:');
                                 (* write the letters that matched: *)
                                 for p1 := 1 to p-1 
                                     do write(output,letters[p1]);
                                 (* write the offending letter: *)
                                 write(output, ch);
                                 (* get the rest of the name and show it: *)
                                 done := eoln(inst);
                                 while not done do begin
                                   done := eoln(inst);
                                   if not done then begin
                                     read(inst,ch); 
                                     if (ch = ' ') or (ch = ';') 
                                     then done := true;
                                     if not done then write(output,ch);
                                   end; 
                                 end;
                                 writeln(output);
                                 (* mark the first letter
                                    that does not match: *)
                                 for p1 := 1 to p-1 do write(output,' ');
                                 write(output,'^'); 
                                 writeln(output);
                                 halt
                              end;
                           end
                        end;
                     end;

                  end
               end
            end
         end;

         if (alignedbase <= -maximumrange) or
            (alignedbase > length + maximumrange) then begin
            writeln(output,' In procedure align:');
            writeln(output,' read in base was ',thebase:1);
            writeln(output,' in internal coordinates: ',alignedbase:1);
            writeln(output,' maximum range was ',maximumrange:1);
            writeln(output,' piece length was ',length:1);
            with pie^.key.hea.keynam
               do writeln(output,' piece name: ',letters:length);
            writeln(output,' piece number: ',number:1);
            writeln(output,' aligned base is too far away... see the code');
            halt
         end
      end
   end
end;
(* end module align.align version = 7.67; {of delmod.p 2004 Sep 8} *)

(* ************************************************************************ *)
(* end module package.align version = 7.67; {of delmod.p 2004 Sep 8} *)

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

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

(* 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.67; {of delmod.p 2004 Sep 8} *)

(* begin module book.name *)
procedure clearname(var n: name);
(* clear the name n *)
var
  i: integer; (* index to piece name *)
begin
  n.length := 0;
  for i := 1 to namelength do n.letters[i] := ' ';
end;

procedure writename(var f: text; n: name);
(* write the name n to file f *)
var
  i: integer; (* index to piece name *)
begin
   for i := 1 to n.length do
       write(f,n.letters[i]);
end;

procedure copyname(a: name; var b: name);
(* copy name a to name b *)
var
  i: integer; (* index to piece name *)
begin
   for i := 1 to a.length
   do b.letters[i] := a.letters[i];
   b.length := a.length;
end;

function equalname(a,b: name): boolean;
(* is name a equal to name b? *)
var
  i: integer; (* index to piece name *)
  same: boolean; (* temporary variable to hold the answer *)
begin
   same := true; (* what optimism!! *)
   if b.length = a.length
   then begin
      i := 1;
      while same and (i <= a.length) do begin
         same := (b.letters[i] = a.letters[i]);
         i := succ(i);
      end
   end
   else same := false;
   equalname := same;
end;
(* end module book.name version = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module book.brorgkey *)
procedure brorgkey(var thefile: text; var theline: integer; var org: orgkey);
(* read organism key *)
begin
      with org do begin
{bbb}
         brheader(thefile,theline,hea);
         getline(mapunit);
         brline(thefile,theline,mapunit);
      end
end;
(* end module book.brorgkey version = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module book.brchrkey *)
procedure brchrkey(var thefile: text; var theline: integer; var chr: chrkey);
(* read chromosome key *)
begin
      with chr do begin
{bbb}
         brheader(thefile,theline,hea);
         brreanum(thefile,theline,mapbeg);
         brreanum(thefile,theline,mapend);
      end
end;
(* end module book.brchrkey version = 7.67; {of delmod.p 2004 Sep 8} *)
 
(* begin module book.getocp *)
procedure getocp(var thefile: text; var theline: integer;
                 var org: orgkey;   var orgchange, orgopen: boolean;
                 var chr: chrkey;   var chrchange, chropen: boolean;
                 var pie: pieceptr; var piechange, pieopen: boolean);
(* Get the next piece and its organism and chromosome keys.  The three change
variables indicate whether or not a new organism, chromosome or piece name
was found.  If a piece is not found, then pieopen will be false.  orgopen,
chropen and pieopen are used by getocp to tell when it has entered an
organism, chromosome or piece.  All booleans should be set to false
initially.  There should be one triplet for each book read.

It is important to initialize ALL variables, including pie:

   orgchange := false; orgopen := false;
   chrchange := false; chropen := false;
   piechange := false; pieopen := false;
   pie := nil;
   theline := 0;

1999 June 2  The book reading routines now treat data objects more
precisely.  Rather than test for eof, the endo of book occurs when pieopen is
returned as false.  A book reading loop now looks like this:

   repeat
      getocp(book, theline,
                   org, orgchange, orgopen,
                   chr, chrchange, chropen,
                   pie, piechange, pieopen);
      writeln(output,'pieopen: ',pieopen);
      if pieopen then begin
         writeln(output,'piece at line: ',theline:1);
      end;
   until not pieopen;
*)

var
      ch: char;
      newchr: chrkey;
      neworg: orgkey;
      newpie: pieceptr;
begin
      ch:='a';
      while not (ch in [' ','p']) do begin
         ch:=getto(thefile,theline,['o','c','p']);
         if ch <> ' ' then begin
            case ch of
               'o': if orgopen
                  then begin
                     readln(thefile);
                     (* move past the word 'organism'
                        - new definition 1999 Mar 13 *)
                     orgopen:=false (* close organism *)
                  end
                  else begin
                     brorgkey(thefile,theline,neworg);
                     if (neworg.hea.keynam.letters <> org.hea.keynam.letters) and
                        (neworg.hea.keynam.length <> org.hea.keynam.length)
                     then begin
{
writeln(output,'--------orgchanged!');
write  (output,'--------old org:"',   org.hea.keynam.letters);
writeln(output,              '" ',    org.hea.keynam.length:1);
write  (output,'--------new org:"',neworg.hea.keynam.letters);
writeln(output,              '" ',neworg.hea.keynam.length:1);
}
(*ccc*)
                        orgchange:=true;
                        copyheader(neworg.hea,org.hea);
                        (* move the mapunit over to the org! *)
                        org.mapunit := neworg.mapunit;
                        clearline(neworg.mapunit);

                     end
                     else orgchange:=false;
                     orgopen:=true;
                  end;
               'c': if chropen
                  then begin
                     readln(thefile);
                     (* move past the word 'chromosome'
                        - new definition 1999 Mar 13 *)
                     chropen:=false (* close chromosome *)
                  end
                  else begin
                     brchrkey(thefile,theline,newchr);
                     if (newchr.hea.keynam.letters <> chr.hea.keynam.letters) and
                        (newchr.hea.keynam.length <> chr.hea.keynam.length)
                     then begin
{
writeln(output,'--------chrchanged!');
write  (output,'--------old chr:"',   chr.hea.keynam.letters);
writeln(output,              '" ',    chr.hea.keynam.length:1);
write  (output,'--------new chr:"',newchr.hea.keynam.letters);
writeln(output,              '" ',newchr.hea.keynam.length:1);
}
                        chrchange:=true;
                        copyheader(newchr.hea,chr.hea);
                        (* move the map range over to the chr! *)
                        chr.mapbeg := newchr.mapbeg;
                        chr.mapend := newchr.mapend;
                     end
                     else chrchange:=false;
                     chropen:=true;
                  end;
               'p': if pieopen
                  then begin
                     pieopen:=false; (* close last piece *)
                     ch:='a' (* prevent falling out of the loop *)
                  end
                  else begin
                     new(newpie);
                     brpiece(thefile,theline,newpie);
                     if pie = nil
                     then piechange := true
                     else begin
                        if (newpie^.key.hea.keynam.letters
                         <> pie^.key.hea.keynam.letters) and
                           (newpie^.key.hea.keynam.length
                         <> pie^.key.hea.keynam.length)
                        then begin
                           piechange:=true;
                        end
                        else piechange:=false;
                     end;
                     pieopen:=true;
                     (* we always have to switch over to the new piece,
                        because although the name may be the same,
                        the DNA sequence could be different.  That is,
                        the book may contain two pieces with the same
                        name, and we want to be sure to search the
                        new one, not the old one. *)
                     if pie <> nil then begin
                        clearpiece(pie); (* save the links *)
                        dispose(pie); (* close up shop *)
                     end;
                     pie := newpie;
                  end
            end
         end
         else begin
            pieopen := false
         end
      end
end;
(* origin: search version = 6.39 *)
(* end module book.getocp version = 7.67; {of delmod.p 2004 Sep 8} *)

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

(* begin module book.iwname *)
procedure iwname(var thefile: text; thename: name);
(* write the name to the file *)
var   c: integer;
begin
      for c:=1 to thename.length
         do write(thefile, thename.letters[c])
end;
(* end module book.iwname version = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module book.iworg *)
procedure iworg(var afile: text; org: orgkey);
(* write an organism specification.  no writeln is done to allow
   write orgchr to do this. *)
begin
      write(afile,'organism ');
      iwname(afile,org.hea.keynam);
      write(afile,';');
end;
(* end module book.iworg version = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module book.iwchr *)
procedure iwchr(var afile: text; chr: chrkey);
(* write an chromosome specification.  no writeln is done to allow
   write orgchr to do this. *)
begin
      write(afile,'chromosome ');
      iwname(afile,chr.hea.keynam);
      write(afile,';');
end;
(* end module book.iwchr version = 7.67; {of delmod.p 2004 Sep 8} *)
 
(* begin module book.iwpie *)
procedure iwpie(var afile: text; pie: piekey);
(* write a piece specification *)
begin
      write(afile,'piece ');
      iwname(afile,pie.hea.keynam);
      writeln(afile,';');
end;
(* end module book.iwpie version = 7.67; {of delmod.p 2004 Sep 8} *)
 
(* begin module book.iworgchr *)
procedure iworgchr(var afile: text;
                       org: orgkey; orgchange, orgopen: boolean;
                       chr: chrkey; chrchange, chropen: boolean);
(* write both organism and chromosome specifications, based
on whether the organism or chromosome changed (orgchange and chrchange)
and whether they are currently open (orgopen, chropen).
See getocp in the br routines. *)
begin
   if orgchange and orgopen then iworg(afile,org);
   if orgchange and chrchange and orgopen and chropen then write(afile,' ');
   if chrchange and chropen then iwchr(afile,chr);
   if (orgchange and orgopen) or (chrchange and chropen) then writeln(afile)
end;
(* end module book.iworgchr version = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module book.iwget *)
procedure iwget(var afile: text; pie: pieceptr;
                    fromplace, pieceplace, toplace: integer;
                    flip: boolean; insttype: integer;
                    carriagereturn: boolean);
(* print a get delila instruction in the orientation of pie,
   from fromplace to toplace pieceplace. +/- direction;
If flip is false, the piece direction is as on the piece,
if it is true, the it is the opposite direction.
insttype: instruction type.
          insttype=1 means the form   get from p  -/+f to same +/-t dir +/-;
          insttype=2 means the form   get from p1 -/+f to p2   +/-t dir +/-;
          where p, p1 and p2 are locations
carriagereturn: if true, add a carriage return to the end of the line.
 *)
procedure iwposition(relative: integer; sameallowed: boolean);
(* write the <POSITION> *)
procedure iwrelative(relative: integer);
begin
           if relative>=0 then write(afile,' +', relative:1)
      else if relative< 0 then write(afile,' ',relative:1)
end;
begin (* iwposition *)
      if (insttype = 1) and sameallowed
      then write(afile,' same')
      else write(afile,' ',pieceplace:1);
      case pie^.key.piedir of
         plus:  iwrelative(+relative);
         minus: iwrelative(-relative);
      end
end;
begin (* iwget *)
      write(afile,'get from');
      iwposition(fromplace, false);
      write(afile,' to');
      iwposition(toplace, true);
      write(afile,' direction');
      case pie^.key.piedir of
                  dirhomologous, (* handle case, may not be right *)
         plus:  case flip of
                  false: write(afile,' +');
                   true: write(afile,' -');
                end;
                  dircomplement, (* handle case, may not be right *)
         minus: case flip of
                  false: write(afile,' -');
                   true: write(afile,' +');
                end
      end;         
      write(afile,';');
      if carriagereturn then writeln(afile);
end;
(* end module book.iwget version = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module book.iwgetsimple *)
procedure iwgetsimple(var afile: text; pie: pieceptr;
                      fromplace, toplace: integer;
                      flip: boolean;
                      carriagereturn: boolean);
(* print a get delila instruction in the orientation of pie,
   from fromplace to toplace pieceplace. +/- direction;
If flip is false, the piece direction is as on the piece,
if it is true, the it is the opposite direction.
Format:
   get from [pieceplace+fromplace]
         to [pieceplace+  toplace] +/- direction;
where [] means to precompute the value.
carriagereturn: if true, add a carriage return to the end of the line.
 *)
begin (* iwget *)
      write(afile,'get from');
      write(afile,' ',fromplace:1);
      write(afile,' to');
      write(afile,' ',  toplace:1);
      write(afile,' direction');
      case pie^.key.piedir of
                  dirhomologous, (* handle case, may not be right *)
         plus:  case flip of
                  false: write(afile,' +');
                   true: write(afile,' -');
                end;
                  dircomplement, (* handle case, may not be right *)
         minus: case flip of
                  false: write(afile,' -');
                   true: write(afile,' +');
                end
      end;         
      write(afile,';');
      if carriagereturn then writeln(afile);
end;
(* end module book.iwgetsimple version = 7.67; {of delmod.p 2004 Sep 8} *)

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

(* begin module package.numbar *)
(* ************************************************************************ *)
(* 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 = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* 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 = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* begin module numberbar *)
function firstlastmax(firstnumber, lastnumber: integer): integer;
(* compute the sizes of firstnumber and lastnumber (including + or -
sign) and then determine which number is larger *)
var
   firstlines: integer; (* number of lines needed for firstumber *)
   lastlines: integer; (* number of lines needed for lastnumber *)
begin
   firstlines := numbersize(firstnumber);
   if firstnumber > 0
   then firstlines := firstlines + 1; (* add one more for + sign *)
   lastlines := numbersize(lastnumber);
   if lastnumber > 0
   then lastlines := lastlines + 1; (* add one more for + sign *)
   if firstlines > lastlines
   then firstlastmax := firstlines
   else firstlastmax := lastlines;
end;

procedure numberbar(var afile: text; spaces, firstnumber, lastnumber: integer;
                    var linesused: integer);
(* write a bar of numbers to a file, with several spaces before.
   the number of lines used is returned *)
var
   logplace: integer; (* the log of the digit being looked at *)
   number: integer; (* the current number being written *)
   spacecount: integer; (* count of spaces *)
begin
{ 2000 June 24:  This code was not sufficient to deal with the sign
correctly.  The numbersize routine now does *not* give the + sign
(which makes it useful for other purposes) so we have to account
for it here now.
   if abs(firstnumber) > abs(lastnumber)
      then linesused:= numbersize(firstnumber)
      else linesused:= numbersize(lastnumber);
}
   linesused := firstlastmax(firstnumber, lastnumber);
{
writeln(output,'numberbar says linesused = ',linesused:1);
      for logplace:=linesused-1 downto 0 do begin
1999 July 15: this is changed to linesused since numbersize now
accounts for the sign:
   for logplace := linesused downto 0 do begin
2000 June 24:  back to the old code ...
}
   for logplace := linesused-1 downto 0 do begin
      for spacecount:=1 to spaces do write(afile,' ');
      for number:=firstnumber to lastnumber
         do write(afile,numberdigit(number,logplace));
      writeln(afile)
   end
end;
(* end module numberbar version = 4.86; (@ of prgmod.p 2004 Sep 8 *)
(* ************************************************************************ *)
(* end module package.numbar version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* begin module space *)
procedure space(var afile: text; spaces: integer);
(* space over some amount to file afile. *)
var
   spacecount: integer; (* count of spaces *)
begin
   for spacecount:=1 to spaces do write(afile,' ');
end;
(* end module space *)

(* begin module sebo.readparameters *)
procedure readparameters(var sebop: text; (* the parameter file to read *)
                         var parameterversion: real;
                         var alignmenttype: char;
                         var fromdo, todo: integer;
                         var rangecontrol: char;
                         var frominst, toinst: integer;
                         var mismatches: integer;
                         var alertinterval: integer);
(* read the parameters as defined in the main routine. *)
begin
   reset(sebop);
   readln(sebop, parameterversion);
{
writeln(output,'parameterversion:', parameterversion:40:38);
writeln(output,'updateversion:',    updateversion:40:38);
}
   if round(100*parameterversion) < round(100*updateversion) then begin
      writeln(output, 'You have an old parameter file!');
      halt
   end;

   readln(sebop, alignmenttype);

   readln(sebop, fromdo, todo);
   if fromdo > todo then begin
      writeln(output,'fromdo (',fromdo:1,
                     ') must be less than todo (',todo:1,')');
      halt
   end;

   if sebop^ in ['a','r'] then begin
      readln(sebop, rangecontrol);
      frominst := 0;
      toinst := 0;
   end
   else begin
      readln(sebop, frominst, toinst);
      if frominst > toinst then begin
         writeln(output,'frominst (',frominst:1,
                        ') must be less than toinst (',toinst:1,')');
         halt
      end;
   end;

   readln(sebop, mismatches);
   if mismatches < 0 then begin
      writeln(output,'mismatches must be non-negative');
      halt
   end;

   readln(sebop, alertinterval);

end;
(* end module sebo.readparameters *)


(* begin module sebo.writeparameters *)
procedure writeparameters(var afile: text;
                              alignmenttype: char;
                              fromdo, todo: integer;
                              rangecontrol: char;
                              frominst, toinst: integer;
                              mismatches: integer;
                              alertinterval: integer);
(* write the parameters to afile. *)
const
   wid = 10; (* width of output numbers *)
begin
   case alignmenttype of
      'f': writeln(afile,'f: force alignment on first base');
      'i': writeln(afile,'i: use the alignedbase from the book');
      'b': writeln(afile,'b: use the internal book');
   end;
   writeln(afile, fromdo:wid, ' ',todo:wid,
                  ' range of book alignment used for search');
   if not (rangecontrol in ['a','r']) then begin
      writeln(afile, frominst:wid, ' ',toinst:wid,
                    ' range of resulting output instructions');
   end
   else begin
      write  (afile, rangecontrol, ' ':wid, ' ':wid);
      if rangecontrol = 'r' then write(afile,' RELATIVE');
      if rangecontrol = 'a' then write(afile,' ABSOLUTE');
      writeln(afile,' instructions');
   end;
   writeln(afile, mismatches:wid,   ' ', ' ':wid, ' mismatches allowed');
   writeln(afile, alertinterval:wid,' ', ' ':wid, ' alert interval');
end;
(* end module sebo.writeparameters *)

(* begin module identifypiece *)
procedure identifypiece(var afile: text; var Ppie: pieceptr;
                            carriagereturn: boolean);
(* identify the piece Ppie by delila instructions;
write result to afile.
carriagereturn: if true, add a carriage return to the end of the line.
example:
piece U00096; get from 637551 to 638151 direction +;
*)
var
   Pflip: boolean; (* is the piece flipped? *)
begin
   with Ppie^.key.hea do begin
      if fulnam <> nil then if fulnam^.length > 0 then begin
         write(afile, 'name "');
         writeline(afile, fulnam, false);
         write(afile, '"; ');
      end;
   end;
   Pflip := (Ppie^.key.piedir <> Ppie^.key.coodir);
   write(afile, 'piece ');
   iwname(afile, Ppie^.key.hea.keynam);
{zz}
   write(afile, '; ');
   if Pflip
   then iwgetsimple(afile, Ppie,
                    inttopie(1, Ppie),
                    inttopie(piecelength(Ppie), Ppie),
                    not Pflip, false)  (* with cr *)
   else iwgetsimple(afile, Ppie,
                    inttopie(1, Ppie),
                    inttopie(piecelength(Ppie), Ppie),
                    Pflip, false); (* with cr *)
   if carriagereturn then writeln(afile);
end;
(* end module identifypiece *)

(* begin module sebo.themain *)
procedure themain(var inst, book, target, sebop, instout, list: text);
(* the main procedure of the program *)
var

   alertinterval: integer; (* interval in bases to tell the user we are still
                              hard at work *)
   alignedbase: integer; (* the aligning base of Bpie *)
   alignmenttype: char; (* 'f' means alignment by First internal coordinate
      base, 'b' means alignment by Book, 'i' means alignment by Instructions *)
   done: boolean; (* done searching this Bpie against this Tpie? *)
   parameterversion: real; (* parameter version number *)
   Bflip: boolean; (* probe direction = target direction? *)
   Tflip: boolean; (* target pie = target coo? *)
   found: boolean; (* the sequence was found in the target? *)
   fromdo, todo: integer; (* user defined range to use for searching *)
   frominst, toinst: integer; (* range of resulting instout instructions *)
   hits: integer; (* number of probes found in target *)
   indent: integer; (* amount to indent book pieces on display *)
   insttype: integer; (* type of instructions to write out; see iwget *)
   l: integer; (* aligned position *)
   linesused: integer; (* line count from numberbar *)
   mismatches: integer; (* number of mismatches allowed *)
   misses: integer; (* number of mismatches found so far *)
   probenumber: integer; (* count of probes *)
   rangecontrol: char; (* control of output range.  'a': absolute
      coordinates; 'r': relative; other: frominst to toinst *)
   targetnumber: integer; (* count of targets *)
   thefulnam: lineptr; (* the current name to use in instout.  It will
                     default to the previous name. *)
   totalprobes: integer; (* total number of probes found in target *)

   (* book pieces of DNA to search with *)
   Btheline: integer; (* location in the book *)
   Borg: orgkey;   Borgchange, Borgopen: boolean;
   Bchr: chrkey;   Bchrchange, Bchropen: boolean;
   Bpie: pieceptr; Bpiechange, Bpieopen: boolean;
   Blength: integer; (* length of Bpie *)
   Bbase: base; (* a base in the book sequence *)
   Bposition: integer; (* position on Bpie *)

   (* target piece to search *)
   Ttheline: integer; (* location in the target book *)
   Torg: orgkey;   Torgchange, Torgopen: boolean;
   Tchr: chrkey;   Tchrchange, Tchropen: boolean;
   Tpie: pieceptr; Tpiechange, Tpieopen: boolean;
   Tlength: integer; (* length of Tpie *)
   Tbase: base; (* a base in the target book sequence *)
   Tposition: integer; (* position on Tpie *)

procedure showalignment(var fout: text; t: integer);
(* show the alignment to file fout for coordinate t *)
var
   index: integer; (* index to the book sequence *)
begin

   (* show the probe sequence *)
   writeln(fout);
   numberbar(fout, indent, fromdo, todo, linesused);
   space(fout, indent);
   for index := fromdo to todo do begin
      if withinalignment(index, alignedbase, Blength)
      then begin
         Bbase := getbase(index+alignedbase, Bpie);
         write(fout, basetochar(Bbase));
      end
      else begin
         write(fout, ' ');
      end;
   end;
   writeln(fout, ' probe');

   (* show the probe-target differences *)
   if misses > 0 then begin
      space(fout, indent);
      for index := fromdo to todo do begin
         if withinalignment(index, alignedbase, Blength)
            and
            withinalignment(index, t, Tlength)
         then begin
            if getbase(index+alignedbase, Bpie)
               =
               getbase(index+t, Tpie)
            then write(fout, ' ')
            else write(fout, 'X');
         end
         else begin
            write(fout, ' ');
         end;
      end;
      writeln(fout, ' differences');
   end;

   (* show the target sequence if it differs from the probe *)
   if misses > 0 then begin
      space(fout, indent);
      for index := fromdo to todo do begin
         if withinalignment(index, t, Tlength)
         then begin
            Tbase := getbase(index+t, Tpie);
            write(fout, basetochar(Tbase));
         end
         else begin
            write(fout, ' ');
         end;
      end;
      writeln(fout, ' target');
   end;

   (* show the zero base location *)
   space(fout, indent);
   space(fout, abs(fromdo));
   write(fout,'^ found at ', inttopie(t, Tpie):1);
   writeln(fout,' on the target');

   writeln(fout);

end; (* showalignment *)

procedure copyline(a: lineptr; var b: lineptr);
(* copy the contents of line a to b *)
var
   index: integer; (* index to characters in l *)
begin
   for index := 1 to a^.length
   do b^.letters[index] := a^.letters[index];
   b^.length := a^.length;
end;

procedure writeinstructions(var instout: text; t: integer);
(* write instructions out for coordinate t *)
begin

   writeln(instout);

(* 2009 Mar 3: this code gave entirely wrong coordinates -
8065429 on the E. coli genome! - replace with identifypiece below.
   writeln(instout,'{ BUBBA original source probe: }');
*){
   write(instout,'{');
   if Bpie^.key.hea.fulnam^.length <> 0 then begin
      write(instout,' name "');
      writeline(instout, Bpie^.key.hea.fulnam, false);
      write(instout,'";');
   end;
   write(instout,' piece ');
   iwname(instout, Bpie^.key.hea.keynam);
   write(instout,'; ');
   if Bflip
   then iwgetsimple(instout, Bpie,
                    inttopie(Bpie^.key.pieend,Bpie),
                    inttopie(Bpie^.key.piebeg,Bpie),
                    Bflip, false)  (* with cr *)
   else iwgetsimple(instout, Bpie,
                    inttopie(Bpie^.key.piebeg,Bpie),
                    inttopie(Bpie^.key.pieend,Bpie),
                    Bflip, false); (* with cr *)
} (*   writeln(instout,' }'); *)

   (* 2009 Mar 3 search target as reported to output *)
   writeln(instout, '{Searching target ',targetnumber:1,
   ' with probe ',probenumber:1,': ');
   identifypiece(instout, Bpie, false);
   writeln(instout, ' }');

   writeln(instout,'{ final target instructions: }');

   if Bpie^.key.hea.fulnam^.length <> 0 then begin
      write(instout,'name "');
      writeline(instout, Bpie^.key.hea.fulnam, false);
{zz}
      write(instout,'"; ');
   end
   else begin (* give name instruction from search sequence piece name: *)
      write(instout,'name "');
      iwname(instout, Bpie^.key.hea.keynam);
      write(instout,'"; ');
   end;

{original:
   iwget(instout, Tpie, frominst, inttopie(t, Tpie),
                  toinst, Tflip, insttype, false);
This gave the wrong instruction in the case that the
target was flipped.  Force it to be correct:
}

   if not (rangecontrol in ['a','r']) then begin
      iwget(instout, Tpie, frominst, inttopie(t, Tpie),
                     toinst, false, insttype, true);
   end
   else begin
      case rangecontrol of
      'a': begin (* absolute instructions *)
           iwget(instout, Tpie,
                     frominst,
                               inttopie(t, Tpie),
                     frominst + (piecelength(Bpie)-1),
{ may need orientation stuff here }
                     false, 2, true);
           end;
      'r': begin (* relative instructions *)
           iwget(instout, Tpie,
                     0, (* relative is from the given coordinate *)
                     inttopie(t, Tpie),
                     piecelength(Bpie)-1, (* base 0 counts as one *)
                     false, insttype, true);
           end;
      end;
   end;
(*
writeln(instout);
write(instout,' {');
writeln(instout,'frominst: ', frominst:1);
writeln(instout,'  toinst: ',   toinst:1);

if Tflip
then writeln(instout,'Tflip is true')
else writeln(instout,'Tflip is false');

if Bflip
then writeln(instout,' Bflip is true')
else writeln(instout,' Bflip is false');
write(instout,' }');
*)

end; (* writeinstructions *)

procedure announceTarget(var fout: text);
(* announce the Target sequence to file fout *)
begin
   writeln(fout);
   writeln(fout,'Target sequence: ');
   identifypiece(fout, Tpie, false);
   writeln(fout,' { ',Tlength:1,' bases }');
end; (* announceTarget *)

procedure announceProbe(var fout: text);
(* announce the Probe sequence to file fout *)
begin
   space(fout, indent);
   writeln(fout, 'Searching target ',targetnumber:1,
                 ' with probe ',probenumber:1,': ');
{zz}
   identifypiece(fout, Bpie, true);
end; (* announceProbe *)

procedure usingpreviousname(var fout: text);
(* make an announcement to fout about the probe name *)
begin
   write(fout,'{  Using a previous probe full name: "');
   writeline(fout,thefulnam, false);
   writeln(fout,'" as the full name for the next probe. }');
end; (* usingpreviousname *)

procedure reporthits(var fout: text);
(* make an announcement to fout of the number of probes found *)
begin
   if hits = 0 then begin
      write(fout,'{ probe');

      if thefulnam^.length <> 0 then begin
         write(fout,' "');
         writeline(fout,thefulnam, false);
         write(fout,'"');
      end;

      writeln(fout,' was not found in the target. }');
   end
   else if hits > 1 then begin
      write(fout,'{  probe "');
      writeline(fout,thefulnam, false);
      writeln(fout,'" was found ',hits:1, ' times in the target. }');
   end
end; (* reporthits *)

procedure dosearch;
(* we now have the book piece and the target, so we can search *)
var
   t: integer; (* position on Tpie that we are checking against Bposition *)
begin
   case alignmenttype of
      'f': alignedbase := 1; (* force alignment on first base *)
      'i': ; (* use the alignedbase from the inst *)
      'b': alignedbase := pietoint(0, Bpie); (* use the internal book *)
   end;

   hits := 0;
{
writeln(output,'Blength = ',Blength:1,' =================');
writeln(output,'fromdo = ',fromdo:1,' =================');
writeln(output,'todo = ',todo:1,' =================');
}
   for t := 1 to Tlength do begin
      if alertinterval > 0 then if (t mod alertinterval) = 0
           then writeln(output, '... ',inttopie(t,Tpie):1);

      done := false;
      found := true;
      l := fromdo;
      misses := 0;
      while (l <= todo) and found do begin
         (* location on B of the base to compare: *)
         Bposition := l + alignedbase;

         (* make sure we are using only the book sequence! *)
         if (Bposition >= 1) and (Bposition <= Blength) then begin

            (* location on T of the base to compare: *)
            Tposition := t + l;

            (* circular wrap for target: *)
            if Tposition < 1       then Tposition := Tposition + Tlength;
            if Tposition > Tlength then Tposition := Tposition - Tlength;

            Tbase := getbase(Tposition, Tpie);
            Bbase := getbase(Bposition, Bpie);

            if Bbase <> Tbase then begin
               misses := succ(misses);
               if misses > mismatches then found := false;
            end;
         end;

         l := succ(l);
      end;


      if found then begin
         showalignment(output,t);
         showalignment(list,t);
         writeinstructions(instout,t);
         hits := succ(hits);
         totalprobes := totalprobes + hits;
      end;

   end;
   reporthits(output);
   reporthits(list);
   reporthits(instout);

end; (* dosearch *)

begin (* themain *)
   writeln(output,'sebo ',version:4:2);
   readparameters(sebop, parameterversion,
                         alignmenttype,
                         fromdo, todo,
                         rangecontrol,
                         frominst, toinst,
                         mismatches,
                         alertinterval);
   writeparameters(output,
                         alignmenttype,
                         fromdo, todo,
                         rangecontrol,
                         frominst, toinst,
                         mismatches,
                         alertinterval);

   (* solution to memory leak 2001 Mar 23: *)
   if alignmenttype = 'i' then begin
      reset(inst);
      if eof(inst) then begin
         writeln(output, 'The sebop requests instructions to be used');
         writeln(output, 'but the inst file is empty!');
         halt;
      end;
   end;

   rewrite(instout);
   writeln(instout,'title "sebo ',version:4:2,' search with book";');
   writeln(instout,'{ sebop parameters:');
   writeparameters(instout,
                         alignmenttype,
                         fromdo, todo,
                         rangecontrol,
                         frominst, toinst,
                         mismatches,
                         alertinterval);
   writeln(instout,'}');
   insttype := 1;

   rewrite(list);
   writeln(list,'sebo ',version:4:2,' search with book";');
   writeln(list,'sebop parameters:');
   writeparameters(list,
                         alignmenttype,
                         fromdo, todo,
                         rangecontrol,
                         frominst, toinst,
                         mismatches,
                         alertinterval);

   brinit(target, Ttheline);

   new(Tpie);
   Torgchange := false; Tchrchange := false; Tpiechange := false;
   Torgopen   := false; Tchropen   := false; Tpieopen   := false;

   new(Bpie);
   Borgchange := false; Bchrchange := false; Bpiechange := false;
   Borgopen   := false; Bchropen   := false; Bpieopen   := false;

   linesused := 0;
   indent := 3;
   getline(thefulnam);
   targetnumber := 0;

   totalprobes := 0;

   while not eof(target) do begin

      (* read in the target piece *)
      writeln(output);
      writeln(output, 'reading in target book ...');
      getocp(target, Ttheline,
                     Torg, Torgchange, Torgopen,
                     Tchr, Tchrchange, Tchropen,
                     Tpie, Tpiechange, Tpieopen);
      writeln(output, '... done');
      iworgchr(instout, Torg, Torgchange, Torgopen,
                        Tchr, Tchrchange, Tchropen);

      if not eof(target) then begin
         targetnumber := succ(targetnumber);
         Tlength := piecelength(Tpie);

         write(instout,'piece ');
         writename(instout, Tpie^.key.hea.keynam);
         writeln(instout, '; { the target sequence is ',
                          Tlength:1,' bases }');

         announceTarget(output);
         announceTarget(list);

         reset(inst);
         brinit(book, Btheline);
         probenumber := 0;
         while not eof(book) do begin

            case alignmenttype of
               'i': align(inst, book, Btheline, Bpie, Blength, alignedbase);
               'b','f': begin

                    getocp(book, Btheline,
                           Borg, Borgchange, Borgopen,
                           Bchr, Bchrchange, Bchropen,
                           Bpie, Bpiechange, Bpieopen);
                    Blength := piecelength(Bpie);
               end;
            end;

            if not eof(book) then begin
               probenumber := succ(probenumber);

               (* Determine relative orientation (flip).  We want to report
                  instructions in instout that are how to grab the given
                  fragment from the Target, irrespective of how it was obtained
                  from the probe library.  So we only care if the target pie
                  is the complement of the coo. *)
               Tflip := (Tpie^.key.piedir <> Tpie^.key.coodir);
               Bflip := (Bpie^.key.piedir <> Bpie^.key.coodir);

{
if Tflip
then writeln(output,'Tflip is true')
else writeln(output,'Tflip is false');

if Bflip
then writeln(output,'Bflip is true')
else writeln(output,'Bflip is false');
}

{
write(output,'Bpie^.key.hea.fulnam= "');
writeline(output,Bpie^.key.hea.fulnam, false);
writeln(output,'"');
write(output,'thefulnam         = "');
writeline(output,thefulnam, false);
writeln(output,'"');
}

               (* this is the memory of the non-empty previous probe name: *)
               if Bpie^.key.hea.fulnam^.length <> 0
               then copyline(Bpie^.key.hea.fulnam, thefulnam)
               else if thefulnam^.length <> 0
                    then begin (* retrieve the previous memory *)
                       usingpreviousname(output);
                       usingpreviousname(list);
                       usingpreviousname(instout);
                       copyline(thefulnam, Bpie^.key.hea.fulnam);
                    end;

               (* announce the book *)
               announceProbe(output);
               announceProbe(list);

               if Tlength < Blength then begin
                  writeln(output,'Target length is less than probe in book.');
{
                  halt;
}
               end
               else dosearch;

               (* solution to memory leak 2001 Mar 23: *)
               clearpiece(Bpie);

            end;
         end;
      end
      else begin
         writeln(output,'No further targets.');
      end;
   end;

   writeln(output);
   writeln( output,'probes:       ', probenumber:1);
   writeln( output,'total probes found: ', totalprobes:1);
   writeln(   list,'total probes found: ', totalprobes:1);
   writeln(instout,'{total probes found: ', totalprobes:1,'}');
   writeln(output,'targets: ', targetnumber:1);
   writeln(output,'sebo search completed');

end;
(* end module sebo.themain *)

begin
   themain(inst, book, target, sebop, instout, list);
1: end.
























