program diffribl(ribla, riblb, diffriblp, posdiff, drxyin, output);
(* diffribl: calculate the difference between two ribls

  by Rye Shultzaberger

modified by:

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

and

  Zehua Chen

 *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 1.28; (* of diffribl.p 2005 Jun 3
2005 Jun  3, 1.24: documentation upgrade, comment out Zhang 4D mode
2005 Jun  2, 1.23: documentation upgrade
2005 Jun  1, 1.22: documentation upgrade
2005 May 31, 1.21: documentation upgrade
2005 May 18, 1.20: z and y modes changed, by Zehua Chen 
2005 Apr  8, 1.19: y mode: documentation
2005 Apr  8, 1.18: y mode: is sqrt(8)L
2005 Apr  5, 1.17: y mode: like z but normalized
2004 Sep  9, 1.16: upgrade for GPC
2004 Feb 12, 1.15: Zhang 3D and 4D distances were added.
2004 Jan 22, 1.14: documentation: Zhang method
2004 Jan 22, 1.13: bug: in r mode userfrom and userto were not set!!
2002 Dec  3, 1.12: documentation: flexrbs
2002 Jan 23, 1.08: Introduce average scan mode for calctype
2001 Sep 12, 1.07: Euclidean distance, 4 dimensional is now in bits
                   at each position, as is reasonable
                   (ie square root of sum of differences squared) TDS
2001 Sep 12, 1.06: output includes parameters TDS
2001 May  8, 1.05: spelling corrections
2000 Aug 21, 1.04: compiles under gpcc
2000 Jun 15, 1.03: added drxyin
2000 May 30, 1.00 origin *)
updateversion = 1.00; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.diffribl *)
(*
name
   diffribl: calculate the difference between two ribls

synopsis
   diffribl(ribla: in, riblb: in, diffriblp: in,
            posdiff: out, drxyin: out, output: out)

files

   ribla:  The ribl output of the Ri program for the first of 2
   compared ribls

   riblb:  The ribl output of the Ri program for the second of 2
   compared ribls

   posdiff: An output file which can be used with the xyplo program. This
            file lists all the difference or distance values for each
            position. The columns are as follows:

            1. This coordinate is the relative position of ribl A.

            2. This value is the sum of differences or distances at
            that position

           This file is only useful for the non-scrolling calculations.

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

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

      range of calculation (char OR char with integers): 
         This allows for the user to specify the range of the matrix.  The
         user can use the range of the matrix by using an 'r'.  The user can
         specify their own range by using a 'u' and then the desired range.
         The user range must be smaller than the range of the ribl.

      scrolling function (char OR char with integers):
         The matrix can be scrolled over itself by using 'v' and the range of
         the scroll.  To not use the scrolling function, use 'n'.

      calctype: type of calculation (char):
         The user can use one of several types of calculations. 
         In many cases, the units reported are in bits.

        (e) The first, specified by "e" is a measurement of the Euclidean
            distance between two positions in two matrices.  This is done
            with the following equation:

            Positional distance = square root( (A1 - A2)^2 + (C1 - C2)^2 +
                                               (G1 - G2)^2 + (T1 - T2)^2 )

            This positional distance is then summed for all positions giving
            the total sum of positional distances.

            When "e" is used with the scrolling function, it calculates only
            for the overlapping part of the matrices.  This feature can be
            used with both symmetric and asymmetric models.


	(o) The second, specified by "o", is a measurement of the Euclidean
	    distance between two matrices.   As opposed to the calculation
	    done in "e", this treats each matrix as a point in 4^(l)
	    dimensional space.

            Since there is only one point, in this case there are no positional
            differences and so the values in posdiff are reported the
            same as for the "e" option.

	(d) The third, specified by "d", is a measurement of difference
	    between the two matrices.  This is done with the following
	    equation:

            Positional difference = (A1 - A2) + (C1 - C2) +
                                    (G1 - G2) + (T1 - T2)

	    This positional difference is then summed for all positions
	    giving the total sum of positional differences.

	    When "d" is used with the scrolling function, it calculates only
	    for the overlapping part of the matrices.  This feature cannot be
	    used with an asymmetric model.

	(s) The fourth, specified by "s", is a measurement of the
            average response the ribla should make to passing across
            sites in riblb.  This is computed as:
               $\sum_l \sum_b f_b(b,l-offset)*Ri_a(b,l)$

            (This is LaTeX typesetting notation, \sum is sum;
            "_" means subscript.)

            NOTE: the frequency f is computed from the number of bases
            at the given position.

        (z) The fifth, specified by "z", is a measurement of the three
            dimentional distance between the two matrices,
            following Zhang.Zhang1991a and Zhang.Zhang1991b.

            Base frequencies are computed from the ribl data file.
            Then each set of frequencies, a, c, g, and t, for which,
            by definition,

            a + c + g + t = 1                                      (1)

            can be represented in three dimensions as:

            x = (a+g) - (t+c) = 2(a+g) - 1                         (2)
            y = (a+c) - (t+g) = 2(a+c) - 1
            z = (a+t) - (g+c) = 2(a+t) - 1 

            These are three independent variables defined by Zhang.
            They map into a tetrahedron in three dimensions.

            Zhang consideres the above to be a 'reduced'
            coordinate system.  The non-reduced system is:

            X = [sqrt(3)/4] x                                      (3)
            Y = [sqrt(3)/4] y
            Z = [sqrt(3)/4] z

            The positional distance is calculated as:

            positional distance = sqrt (   (X2-X1)^2               (4)
                                         + (Y2-Y1)^2 
                                         + (Z2-z1)^2 )

            where sqrt is the square root.

            From Zhang.Zhang1991a (page 46 and 47), this simplifies
            to:

            positional distance = [sqrt(3)/2]                      (5)

                                  * sqrt( (a2-a1)^2 + (g2-g1)^2 +
                                          (c2-c1)^2 + (t2-t1)^2  )

            where the (a1, g1, c1, t1), (a2, g2, c2 and t2) are
            probabilities of two different matrices at one position.

            This positional distance is then summed for all positions giving
            the total sum of positional distances.

            Because frequencies sum to 1, there really are only three
            independent degrees of freedom and therefore only three
            dimensions.  So equations for three and four dimensions
            give the same results.

        (y) The seventh, specified by "y", is computed as "z" and then
            it is normalized by the maximum possible distance between
            points in the tetrahedrons.

            With frequencies, the largest distance in the Zhang
            tetrahedron is sqrt(3/2), along the edge.  For L
            positions, the largest possible distance is therefore
            sqrt(3/2)L.  All values are divided by this maximum.

            The following shows that the maximum distance in the
            non-reduced coordinate system is sqrt(3/2).  Using
            equations (2) and (3), for the case of all G, the point is
            at

                (X, Y, Z) = (sqrt(3)/4, -sqrt(3)/4, -sqrt(3)/4)

            while for all C, the point is

                (X, Y, Z) = (-sqrt(3)/4, sqrt(3)/4, -sqrt(3)/4)

            The distance between these points is sqrt(3/2).

  drxyin:  This gives the total sum of distance or difference, depending on
	   which calculation function is being used.  When the scrolling
	   function is being used, it will report the total sum value along
	   with the position of the scroll.  The position of the scroll is
	   the distance between the zero coordinates of the two matrices.

   output: messages to the user

description

   This program looks at the differences in two individual information
   weight matricies (ribls) by finding the difference in information
   at each position, for each base, and then summing the differences. 
   Then all of the differences at each position are summed to express
   a diffribl value.

   Actually, the program now has a number of other ways of comparing
   the ribls, depending on a user parameter.

examples

examples of diffriblp

1.24        version of diffribl that this parameter file is designed for.
r u -10 +10 r use the from/to coords from ribl, u means use user specified
n v -21 +21 v and coords=move riblB across riblA for the range, n=no
eodszy      e:Euclid, o:Euclid4d, d:difference, s:scan, z:Zhang, y:znorm

documentation

@article{Shultzaberger.Schneider2001,
author = "R. K. Shultzaberger
 and R. E. Bucheimer
 and K. E. Rudd
 and T. D. Schneider",
title = "{Anatomy of \emph{Escherichia coli}
Ribosome Binding Sites}",
journal = "J. Mol. Biol.",
volume = "313",
pages = "215-228",
comment = "Shultzaberger.Schneider.flexrbs",
note = "\htmladdnormallink
{https://alum.mit.edu/www/toms/paper/flexrbs/}
{https://alum.mit.edu/www/toms/paper/flexrbs/}",
year = "2001"}

see also

   {example parameter file:} diffriblp

   {Description of use is in Shultzaberger.Schneider2001:}
   https://alum.mit.edu/www/toms/paper/flexrbs/

   {program that generates ribls: } ri.p
   {program that uses ribls to find sites:} scan.p
   {graphics program for xyin:} xyplo.p
   {source of program modules:} lister.p

author

   Ryan Shultzaberger
   Thomas D. Schneider
   Zehua Chen

bugs

   There is a problem with comparing different sized ribls.  I (Ryan?)
   need to fix this.  For now, only use this program with same sized
   ribls.  The result will be wrong if done otherwise.

   Comparisons in 4 dimensional space are not appropriate because the
   4 probabilities are not independent.  To avoid this, one can
   replace the 4 dimensional space with a 3 dimensional one according
   to Zhang's methods:

@article{Zhang.Zhang1991a,
author = "C.-T. Zhang
 and R. Zhang",
title = "Diagrammatic representation of the distribution of {DNA}
bases and its applications",
journal = "Int. J. Biol. Macromol.",
volume = "13",
pages = "45-49",
note = "tetrahedron method",
year = "1991"}

@article{Zhang.Zhang1991b,
author = "C.-T. Zhang
 and R. Zhang",
title = "Analysis of distribution of bases in the coding sequences
by a diagrammatic technique",
journal = "Nucleic Acids Res.",
volume = "19",
pages = "6313-6317",
note = "tetrahedron method",
year = "1991"}

@article{Zhang1997,
author = "C.-T. Zhang",
title = "A Symmetrical Theory of {DNA} Sequences and Its Applications",
journal = "J. Theor. Biol.",
volume = "187",
pages = "297-306",
year = "1997"}

   To do this, the Ri can be converted to probabilities according to
   Ri = 2 + log2(Pi).  Then the probabilities are converted to the
   Zhang XYZ space.  Distances are then measured in that XYZ space.
   However it is better to use the Pi directly from the ribl file.

technical notes

*)
(* end module describe.diffribl *)

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

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

type  

(* begin module discan.types *)
   box = record (* values related to distribution *)
      range: real;    (* the numbers that make up the range of the dist *)
      occur: real;    (* the number of occurrences at each position,
                         This can be a real if one is making up a
                         distribution. *)
      probdist: real; (* the log2 (prob for each position) *)
      isitzero: boolean;    (* checks to see if the occurrence is 0 *)
   end;  
 
   map = array[1..1000] of box; (* a set of boxes *)
(* end module discan.types *)

(* begin module base.type *)
   (* define the four nucleotide bases *)
   base = (a,c,g,t);
(* end module base.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

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

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

var
   ribla: text; (* the ribl of model A *)
   riblb: text; (* the ribl of model B *)
   diffriblp: text; (* parameter file for diffribl *)
   posdiff: text; (* a xyplo input file *)
   drxyin: text; (* a xyplo input file *)

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

(* begin module skipblanks *)

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

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

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

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

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

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

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

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

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

(* begin module readstringline *)
procedure readstringline(var afile: text; var buffer: string);
(* get a line (as a string) from a file. *)
var   index: integer; (* of buffer *)
begin (* getstring *)
   clearstring(buffer);
   index := 0;
   while (not eoln(afile)) and (index < maxstring) do begin
      index := succ(index);
      read(afile, buffer.letters[index])
   end;
   if not eoln(afile) then begin
      writeln(output, 'readstringline: a line exceeds maximum string size (',
                      maxstring:1,')');
      halt
   end;
   buffer.length := index;
   buffer.current := 1;
   readln(afile);
end; (* readstringline *)
(* end module readstringline version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

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

(* begin module readquotestring *)
procedure readquotestring(var afile: text; var s: string);
(* read the quoted string from afile, remove the quotes.  Trailing
blanks are converted to tab characters. *)
const
   tabcharacter = 9; (* the ordinal of the tab character *)
var
   c: char; (* a character read *)
   i: integer; (* index to the string s *)
begin
   emptystring(s);
   skipblanks(afile);
   with s do begin
      read(afile,c);
      if c <> '"' then begin
         writeln(output, 'quote string expected but " missing');
         writeln(output, 'character found instead:"',c,'"');
         halt
      end;
      if eoln(afile) then begin
         writeln(output, 'missing end of quote string " string');
         halt
      end;
      repeat
         read(afile,c);
         if c <> '"' then begin
            length := succ(length);
            letters[length] := c;
         end;
      until (c = '"')
            or eoln(afile) or (length = maxstring);
      if length = maxstring then begin
         writeln(output, 'A quote string exceeded ',maxstring:1,
                         ' characters.');
         writeln(output, 'Use a shorter string or increase constant',
                         ' maxstring.');
         halt
      end;
      if c <> '"' then begin
         writeln(output, 'The second " mark of a string is missing.');
         halt
      end;
      (* convert all blanks to tab characters *)
      for i := 1 to length do begin
         if letters[i] = ' '
         then letters[i] := chr(tabcharacter); (* set to tab character *)
      end;
   end;
end;
(* end module readquotestring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module readheader *)
procedure readheader(var afile: text; ch: char; var theheader: stringptr);
(* Read the header from a file.  Only read lines that begin with ch.
2004 Dec 15: Handle end of file (eof) gracefully. *)
var
   h: stringptr; (* pointer to the riblheader *)
begin
   if eof(afile) then begin
      clearstring(theheader^);
   end
   else begin
      if afile^ = ch then begin
         new(theheader);
         h := theheader;
         while afile^='*' do begin
            readstringline(afile, h^);
            if not eof(afile) then begin
               if afile^='*' then begin
                  new(h^.next);
                  h := h^.next;
               end
               else h^.next := nil;
            end
         end
      end
      else begin
         theheader := nil
      end;
   end;
end;
(* end module readheader version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module readwaveparameters *)
procedure readawaveparameter(var afile: text; var wp: waveparam);
{
                                 cmperbase, cmperbit: real);
}
(* Read one wave parameter from a file. *)
{
2000 Nov 1:  The following experimental code is not implemented:
See the coscurve program
(https://alum.mit.edu/www/toms/delila/coscurve.html) for details
about the conversion from bases to cm.  The code now requires
cm per base and cmperbit.
const
   lover2pi = 1.2160067233; (* L/2pi for L being the length of
       the cosine function from 0 to 2pi *)
var
   multiplier: real; (* factor to multiply by to convert from bases to cm. *)
}
begin
   readln(afile,wp.extreme);
   readln(afile,wp.wavelocation);
   readln(afile,wp.wavebit);
   readln(afile,wp.waveamplitude);
   readln(afile,wp.wavelength);
   if wp.wavelength <= 0.0 then begin
      writeln(output,'wave parameters: wavelength must be positive');
      halt
   end;
   (* one may skip the dash and thickness parameters
      of the last wave definition.  This allows
      backwards compatability, but should not be
      used in general. *)
   if eof(afile)
   then begin
      wp.dashon := 0;
      wp.dashoff := 0;
      wp.dashoffset := 0;
      wp.thickness := 0
   end
   else begin
      (* implement new method to read dash information
         without breaking previous wave parameter files *)
      if afile^ <> 'd'
      then begin (* default to old method *)
         readln(afile,wp.dashon);
         wp.dashoff := wp.dashon;
         wp.dashoffset := 0;
      end
      else begin
         get(afile); (* skip the 'd' *)
         (* reading the parameters directly in cm would require
            that the user compute or guess the values.
            So read in as bases and compute cm.
            Unfortunately the length of a curve is not the simple
            scale factor.  Make the user guess for now.
         *)
         readln(afile,wp.dashon, wp.dashoff, wp.dashoffset)

{
experimental code to account for actual wave length:
         readln(afile, wp.BASEdashon, wp.BASEdashoff, wp.BASEdashoffset);
         multiplier    := 1;
         multiplier    := lover2pi
 * wp.waveamplitude (* in bits *)
 * cmperbit (* converts to cm *)
 / wp.wavelength (* 10.6 bp/turn *)
* cmperbase
;
multiplier    := lover2pi;
multiplier    := 1.0;
multiplier    := 1.0/lover2pi;
         multiplier := lover2pi
                       * wp.waveamplitude (* in bits *)
                       * cmperbit (* converts to cm *)
                       * wp.wavelength (* 10.6 bp/turn *)
                       * cmperbase;
         multiplier := 1.0
         wp.dashon     := multiplier * wp.BASEdashon;
         wp.dashoff    := multiplier * wp.BASEdashoff;
         wp.dashoffset := multiplier * wp.BASEdashoffset;
with wp do begin
   writeln(output, 'cmperbase  = ',cmperbase:5:2);
   writeln(output, 'cmperbit   = ',cmperbit:5:2);
   writeln(output, 'dashon     = ',dashon:5:2);
   writeln(output, 'dashoff    = ',dashoff:5:2);
   writeln(output, 'dashoffset = ',dashoffset:5:2);
   writeln(output, 'waveamplitude   = ',wp.waveamplitude:5:2);
   writeln(output, 'wavelength      = ',wp.wavelength:5:2);
end;
}

      end;
      if eof(afile)
      then wp.thickness := 0
      else readln(afile,wp.thickness);
   end;
end;

procedure readwaveparameters(var afile: text; var w: waveptr);
{
                                 cmperbase, cmperbit: real);
}
(* read from afile the wave parameters w.  The routine is done either at end
of file, when a period ('.') appears or at a blank line.  This makes it
compatable with different reading mechanisms. *)
var
   done: boolean; (* done reading *)
   morethanone: boolean; (* more than one wave? *)
   p: waveptr; (* one of the definitions of a cosine wave *)
procedure waystoend;
(* if end of file, a blank line or a period, consider the job done *)
var
   clear: boolean; (* not eof, so clear to test for comment *)
begin
   if eof(afile) then done := true;
   clear := false;
   if not done (* skip any number of comments *)
   then while (not clear) and (not done) do begin
      if eof(afile) then begin
         done := true;
         clear := true;
      end
      else if afile^ = '*'
           then readln(afile)
           else clear := true
   end;
   if not done then if eoln(afile) then done := true;
   if not done then if afile^='.' then begin
      readln(afile);
      done := true;
   end;
end;
begin
   done := false;
   morethanone := false;
   w := nil;
   waystoend;
   if not done then begin
      new(w);
      p := w;
      while not done do begin
         waystoend;
         if done
         then p^.next := nil
         else begin
            if afile^='*'
            then readln(afile) (* skip any number of comments *)
            else begin
               if morethanone then begin
                  new(p^.next);
                  p := p^.next
               end;
{
               readawaveparameter(afile, p^, cmperbase, cmperbit);
}
               readawaveparameter(afile, p^);
               (* if we ever come back again there will be more than one: *)
               morethanone := true;
            end
         end
      end
   end
   else w := nil;
end;
(* end module readwaveparameters version = 5.67; (@ of prgmod.p 2017 Sep 29 *)



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

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

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

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

      readheader(afile,'*',riblheader);

      readln(afile,frombase,tobase);

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

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

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

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

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

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

(* begin module cleanmatrix *)
(* places zero's in the bases outside of the specified matrix *)
procedure cleanmatrix(var matrix: ribltypeptr;
                      frombase, tobase: integer);
var
   x: integer; (* generic variable *)

begin

   x:= minribl; 

   while x < frombase do begin
      matrix^.data[a,x]:= 0; 
      matrix^.data[c,x]:= 0; 
      matrix^.data[g,x]:= 0; 
      matrix^.data[t,x]:= 0; 
      x:= x + 1;
   end;

   while (x >= frombase) and (x<= tobase) do begin
      x:=x+1;
   end;

   while (x > tobase) and (x <= maxribl) do begin
      matrix^.data[a,x]:= 0;
      matrix^.data[c,x]:= 0;
      matrix^.data[g,x]:= 0;
      matrix^.data[t,x]:= 0;
      x:= x + 1; 
   end;

end;

(* end module cleanmatrix *)


(* begin module subribl *)
procedure subribl(matrixA, matrixB: ribltypeptr; var matrixC: ribltypeptr;
                  var posdiff, drxyin: text; space: integer;
                  rangebool: boolean; userfrom, userto: integer);
(* calculate the difference between points *)

var
   x: integer; (* counter *)
   sum: real; (* sum of differences at position *)
   totalsum: real; (*the sum of all the differences for all the positions *)

begin
   totalsum:= 0;

   if rangebool = false then begin
      for x:= (matrixA^.frombase - abs(space)) to (matrixA^.tobase + abs(space)) do begin
         matrixC^.data[a,x]:= (matrixA^.data[a,x] - matrixB^.data[a,x - space]);
         matrixC^.data[c,x]:= (matrixA^.data[c,x] - matrixB^.data[c,x - space]);
         matrixC^.data[g,x]:= (matrixA^.data[g,x] - matrixB^.data[g,x - space]);
         matrixC^.data[t,x]:= (matrixA^.data[t,x] - matrixB^.data[t,x - space]);
         sum:= (matrixC^.data[a,x] + matrixC^.data[c,x] + matrixC^.data[g,x]
         + matrixC^.data[t,x]);
         write(posdiff,x:10);
         writeln(posdiff,sum:15:6);
 
         totalsum:= totalsum + sum;
      end;
      writeln(posdiff, '*');
      writeln(posdiff,'* the sum of all the differences is: ', totalsum:15:6,
                      ' bits');

      writeln(drxyin,space:10,totalsum:15:6);

      if space = 0 then begin
         writeln(output,'The sum of all distances is: ',totalsum:10:6,
                        ' bits');
      end;

   end; 

   if rangebool then begin
      for x:= (userfrom - abs(space)) to (userto + abs(space)) do begin
         matrixC^.data[a,x]:= (matrixA^.data[a,x] - matrixB^.data[a,x - space]);
         matrixC^.data[c,x]:= (matrixA^.data[c,x] - matrixB^.data[c,x - space]);
         matrixC^.data[g,x]:= (matrixA^.data[g,x] - matrixB^.data[g,x - space]);
         matrixC^.data[t,x]:= (matrixA^.data[t,x] - matrixB^.data[t,x - space]);
         sum:= (matrixC^.data[a,x] + matrixC^.data[c,x] + matrixC^.data[g,x]
         + matrixC^.data[t,x]);
         write(posdiff,x:10);
         writeln(posdiff,sum:15:6);
 
         totalsum:= totalsum + sum;
      end;
      writeln(posdiff, '*');
      writeln(posdiff,'* the sum of all the differences is: ',
                      totalsum:15:6, 'bits');
     
      writeln(drxyin,space:10,totalsum:15:6);

      if space = 0 then begin
         writeln(output,'The sum of all distances is: ',
                        totalsum:10:6, 'bits');

      end;

   end;

end;

(* end module subribl *)

{
zzz
}

(* begin module rifreq *)
function rifreq(ri: real): real;
(* compute the frequency by inverting the Ri value.
NO SMALL SAMPLE CORRECTION is used.
    Ri = 2 + log_2 (f)
so
    f = 2^(Ri-2)
or
    f = (Ri-2)  I THINK THIS IS WRONG
note that
2^x = exp(ln(2^x)) = exp(x*ln(2))
THIS CAN BE EASILY AVOIDED BY USING THE ACTUAL NUMBERS!
*)
begin
   rifreq := exp((ri-2.0)*ln(2))
end;
(* end module rifreq *)

(* begin module scanribl *)
procedure scanribl(matrixA, matrixB: ribltypeptr; var matrixC: ribltypeptr;
                  var posdiff, drxyin: text; space: integer;
                  rangebool: boolean; userfrom, userto: integer);
(* calculate the scan average between matrixA and the frequencies of matrixB.
by TDS *)

var
   b: base; (* index to bases *)
   x: integer; (* counter *)
   sum: real; (* sum of averages at position *)
   totalsum: real; (* the sum of all the differences for all the positions *)
   thefrom, theto: integer; (* range to use *)
   nl: integer; (* numbers of bases at a position l *)

begin

{
writeln(output,'scan average at spacing: ',space:1);
}

   totalsum := 0.0;

   if rangebool then begin
      thefrom := matrixA^.frombase;
      theto   := matrixA^.tobase;
   end
   else begin
      thefrom := userfrom;
      theto   := userto;
   end;

{
writeln(output,'   thefrom: ',   thefrom:1);
writeln(output,'     theto: ',     theto:1);
}

   for x:= (thefrom - abs(space)) to (theto + abs(space)) do begin
      nl := 0;
      for b := a to t do begin
         nl := nl + matrixB^.numbers[b,x - space];
{
if space = 0 then writeln(output,'matrixB^.numbers[b,x - space] = ',
         matrixB^.numbers[b,x - space]:1);
}
      end;

      sum := 0.0;
      for b := a to t do begin
{
         matrixC^.data[b,x] :=   matrixA^.data[b,x]
                               * rifreq(matrixB^.data[b,x - space]);
}
         if matrixB^.numbers[b,x - space] <> 0 then begin

            matrixC^.data[b,x] := matrixA^.data[b,x]
                                  * matrixB^.numbers[b,x - space]/ nl;

            (* one gets the same result if the matricies are reversed:
            matrixC^.data[b,x] := matrixB^.data[b,x]
                                  * matrixA^.numbers[b,x - space]/ nl;
            *)

            sum := sum + matrixC^.data[b,x];
         end;
      end;
{
writeln(output,'* at an x of ',x:3, ', the sum is: ', sum:15:6, ' bits');
}

      write(posdiff, x:10);
      writeln(posdiff, sum:15:6);

      totalsum:= totalsum + sum;
   end;
   writeln(posdiff,'* at a spacing of ',space:1,
                   ', the average is: ', totalsum:15:6, ' bits');
   writeln(posdiff, '*');

   writeln(drxyin,space:10, totalsum:15:6);

   if space = 0 then begin
      writeln(output,'At zero spacing, the average is: ',
                     totalsum:10:6, ' bits');
   end;

end;

(* end module scanribl *)


{
zzz
}

(* begin module euribl *)
(* calculates the Euclidean distance between points *)
procedure euribl(matrixA, matrixB: ribltypeptr; var matrixC: ribltypeptr;
                  var posdiff, drxyin: text; space: integer;
                  rangebool: boolean; userfrom, userto: integer);

var
   x: integer; (* counter *)
   sum: real; (* sum of differences at position *)
   totalsum: real; (*the sum of all the differences for all the positions *)

begin

   totalsum:= 0;

   if rangebool = false then begin
      for x:= (matrixA^.frombase - abs(space)) to (matrixA^.tobase + abs(space)) do begin
         matrixC^.data[a,x]:= sqr(matrixA^.data[a,x] - matrixB^.data[a,x - space]);
         matrixC^.data[c,x]:= sqr(matrixA^.data[c,x] - matrixB^.data[c,x - space]);
         matrixC^.data[g,x]:= sqr(matrixA^.data[g,x] - matrixB^.data[g,x - space]);
         matrixC^.data[t,x]:= sqr(matrixA^.data[t,x] - matrixB^.data[t,x - space]);
         sum:= sqrt(matrixC^.data[a,x] + matrixC^.data[c,x] + matrixC^.data[g,x]
         + matrixC^.data[t,x]);
         write(posdiff,x:10);
         writeln(posdiff,sum:15:6);

         totalsum:= totalsum + sum;

      end;
      writeln(posdiff, '*');
      writeln(posdiff,'* the sum of all the Euclidean distances is: ',
                             totalsum:15:6, 'bits');

      writeln(drxyin,space:10,totalsum:15:6);

      if space = 0 then begin
         writeln(posdiff, '*');
         writeln(output,'The sum of all distances is: ',totalsum:10:6, 'bits');
      end;

   end;

   if rangebool then begin
      for x:= (userfrom - abs(space)) to (userto + abs(space)) do begin
         matrixC^.data[a,x]:= sqr(matrixA^.data[a,x] - matrixB^.data[a,x - space]);
         matrixC^.data[c,x]:= sqr(matrixA^.data[c,x] - matrixB^.data[c,x - space]);
         matrixC^.data[g,x]:= sqr(matrixA^.data[g,x] - matrixB^.data[g,x - space]);
         matrixC^.data[t,x]:= sqr(matrixA^.data[t,x] - matrixB^.data[t,x - space]);
         sum:= sqrt(matrixC^.data[a,x] + matrixC^.data[c,x] + matrixC^.data[g,x]
         + matrixC^.data[t,x]);
         write(posdiff,x:10);
         writeln(posdiff,sum:15:6);
 
         totalsum:= totalsum + sum;
      end;  
      writeln(posdiff, '*');
      writeln(posdiff,'* the sum of all the Euclidean distances is: ',
                      totalsum:15:6, ' bits');
         
      writeln(drxyin,space:10,totalsum:15:6);

      if space = 0 then begin
         writeln(output,'The sum of all distances is: ',
                         totalsum:10:6, ' bits');
      end;
   end;

end;


(* end module euribl *)

(* begin module onepointribl *)
(* calculates the Euclidean distance between points *)
procedure onepointribl(matrixA, matrixB: ribltypeptr; var matrixC: ribltypeptr;
                  var posdiff, drxyin: text; space: integer;
                  rangebool: boolean; userfrom, userto: integer);

var
   x: integer; (* counter *)
   sum: real; (* sum of differences at position *)
   totalsum: real; (*the sum of all the differences for all the positions *)
   onedist: real; (* the distance between 2 ribls *)
   ln2: real; (* the ln (2) *)

begin

   totalsum:= 0;
   ln2:= ln(2);

   if rangebool = false then begin
      for x:= (matrixA^.frombase - abs(space)) to (matrixA^.tobase + abs(space)) do begin
         matrixC^.data[a,x]:= sqr(matrixA^.data[a,x] - matrixB^.data[a,x - space]);
         matrixC^.data[g,x]:= sqr(matrixA^.data[g,x] - matrixB^.data[g,x - space]);
         matrixC^.data[t,x]:= sqr(matrixA^.data[t,x] - matrixB^.data[t,x - space]);
         matrixC^.data[c,x]:= sqr(matrixA^.data[c,x] - matrixB^.data[c,x - space]);
         sum:= (matrixC^.data[a,x] + matrixC^.data[c,x] + matrixC^.data[g,x]
         + matrixC^.data[t,x]);
         write(posdiff,x:10);
         writeln(posdiff,sqrt(sum):15:6);

         totalsum:= totalsum + sum;

      end;
      onedist:= sqrt(totalsum);
      writeln(posdiff, '*');
      writeln(posdiff,'* the sum of all the Euclidean distances is: ',
                       onedist:15:6, 'bits');

      writeln(drxyin,space:10,onedist:15:6);

      if space = 0 then begin
         writeln(output,'The sum of all distances is: ',onedist:10:6, ' bits');
      end;

   end;

   if rangebool then begin
      for x:= (userfrom - abs(space)) to (userto + abs(space)) do begin
         matrixC^.data[a,x]:= sqr(matrixA^.data[a,x] - matrixB^.data[a,x - space]);
         matrixC^.data[c,x]:= sqr(matrixA^.data[c,x] - matrixB^.data[c,x - space]);
         matrixC^.data[g,x]:= sqr(matrixA^.data[g,x] - matrixB^.data[g,x - space]);
         matrixC^.data[t,x]:= sqr(matrixA^.data[t,x] - matrixB^.data[t,x - space]);
         sum:= (matrixC^.data[a,x] + matrixC^.data[c,x] + matrixC^.data[g,x]
         + matrixC^.data[t,x]);
         write(posdiff,x:10);
         writeln(posdiff,sum:15:6);
 
         totalsum:= totalsum + sum;
      end;  
      onedist:= sqrt(totalsum);
      writeln(posdiff, '*');
      writeln(posdiff,'* the sum of all the Euclidean distances is: ',
                       onedist:15:6, ' bits');
         
      writeln(drxyin,space:10,onedist:15:6);

      if space = 0 then begin
         writeln(output,'The sum of all distances is: ',onedist:10:6,' bits');
      end;
   end;
end;
(* end module onepointribl *)



(* begin module overlapribl *)
(* calculates the Euclidean distance for only the overlapping matrices *)
procedure overlapribl(matrixA, matrixB: ribltypeptr; var matrixC: ribltypeptr;
                  var posdiff, drxyin: text; space: integer;
                  rangebool: boolean; userfrom, userto: integer);
 
var
   x: integer; (* counter *)
   sum: real; (* sum of differences at position *)
   totalsum: real; (*the sum of all the differences for all the positions *)
   newfrom: integer; (* the new from range value of the overlap *)
   newto: integer; (* the new to range value of the overlap *)
 
begin
 
   totalsum:= 0; 

   if space <= 0 then begin
      if rangebool = false then begin
         newfrom:= matrixA^.frombase;
         newto:=  matrixA^.tobase + space;
      end;

      if rangebool then begin
         newfrom:= userfrom;
         newto:= userto + space;
      end;

   end;

   if space > 0 then begin
      if rangebool = false then begin
         newfrom:= matrixA^.frombase + space;
         newto:=  matrixA^.tobase;
      end;

      if rangebool then begin
         newfrom:= userfrom + space;
         newto:= userto;
      end;

   end;
      
 
   for x:= (newfrom) to (newto) do begin
      matrixC^.data[a,x]:= sqr(matrixA^.data[a,x] - matrixB^.data[a,x - space]);
      matrixC^.data[c,x]:= sqr(matrixA^.data[c,x] - matrixB^.data[c,x - space]);
      matrixC^.data[g,x]:= sqr(matrixA^.data[g,x] - matrixB^.data[g,x - space]);
      matrixC^.data[t,x]:= sqr(matrixA^.data[t,x] - matrixB^.data[t,x - space]);
      sum:= sqrt(matrixC^.data[a,x] + matrixC^.data[c,x] + matrixC^.data[g,x]
            + matrixC^.data[t,x]);
      write(posdiff,x:10);
      writeln(posdiff,sum:15:6);

      totalsum:= totalsum + sum;
   end;

   writeln(posdiff, '*');
   writeln(posdiff,'* the sum of all the Euclidean distances is: ',
                   totalsum:15:6, ' bits');

   writeln(drxyin,space:10,totalsum:15:6); end;

(* end module overlapribl *)


(* begin module zhang3ribl *)
(* calculates the distance in Zhang XYZ space *)
procedure zhang3ribl(matrixA, matrixB: ribltypeptr; var matrixC: ribltypeptr;
                  var posdiff, drxyin: text; space: integer;
                  rangebool: boolean; userfrom, userto: integer;
                  normalize: boolean);

var
   x: integer; (* counter *)
   sum: real; (* sum of differences at position *)
   totalsum: real; (*the sum of all the differences for all the positions *)
   newfrom: integer; (* the new from range value of the overlap *)
   newto: integer; (* the new to range value of the overlap *)
   maxfactor: real; (* the maximum value possible *)
   numA: integer; (* number of sites of matrix A *)
   numB: integer; (* number of sites of matrix B *)
   b: base; (* A, G, C, T *)
   piAa: real; (* probability of A of matrix A *)
   piAg: real; (* probability of G of matrix A *)
   piAc: real; (* probability of C of matrix A *)
   piAt: real; (* probability of T of matrix A *)
   piBa: real; (* probability of A of matrix B *)
   piBg: real; (* probability of G of matrix B *)
   piBc: real; (* probability of C of matrix B *)
   piBt: real; (* probability of T of matrix B *)
   piCa: real; (* sqr(piAa - piBa) *)
   piCc: real; (* sqr(piAc - piBc) *)
   piCg: real; (* sqr(piAg - piBg) *)
   piCt: real; (* sqr(piAt - piBt) *)

begin
   
   totalsum:= 0;
   piAa := 0; piAc := 0; piAg := 0; piAt := 0;
   piBa := 0; piBc := 0; piBg := 0; piBt := 0;
   piCa := 0; piCc := 0; piCg := 0; piCt := 0;
  
   if space <= 0 then begin
      if rangebool = false then begin
         newfrom:= matrixA^.frombase;
         newto:=  matrixA^.tobase + space;
      end;
 
      if rangebool then begin
         newfrom:= userfrom;
         newto:= userto + space;
      end;
 
   end;
 
   if space > 0 then begin
      if rangebool = false then begin
         newfrom:= matrixA^.frombase + space;
         newto:=  matrixA^.tobase;
      end;
 
      if rangebool then begin
         newfrom:= userfrom + space;
         newto:= userto;
      end;
 
   end;
 
   if normalize then begin
      maxfactor := (newto - newfrom -1) * sqrt(3/2)
   end
   else begin
      maxfactor := 1.00;
   end;
 
   for x:= (newfrom) to (newto) do begin
       numA := 0;
       numB := 0;
       for b := a to t do begin
                 numA := numA + matrixA^.numbers[b,x];
                 numB := numB + matrixB^.numbers[b,x];
       end;

       piAa := matrixA^.numbers[a,x]/numA;
       piAc := matrixA^.numbers[c,x]/numA;
       piAg := matrixA^.numbers[g,x]/numA;
       piAt := matrixA^.numbers[t,x]/numA;
       piBa := matrixB^.numbers[a,x]/numB;
       piBc := matrixB^.numbers[c,x]/numB;
       piBg := matrixB^.numbers[g,x]/numB;
       piBt := matrixB^.numbers[t,x]/numB;
       piCa := sqr(piAa - piBa);
       piCc := sqr(piAc - piBc);
       piCg := sqr(piAg - piBg);
       piCt := sqr(piAt - piBt);
       
       sum:= sqrt(3*(piCa + piCc + piCg + piCt))/2;
       write(posdiff,x:10);
       if normalize then begin
           writeln(posdiff,sum/maxfactor:15:6);
       end
       else begin
           writeln(posdiff,sum:15:6);
       end;
 
       totalsum:= totalsum + sum;
   end;
 
   writeln(posdiff, '*');
 
   if normalize then begin
     writeln(posdiff,'* NORMALIZED the sum of all the Euclidean distances is: ',
                        totalsum/maxfactor:15:6);
   end
   else begin
     writeln(posdiff,'* the sum of all the Euclidean distances is: ',
                        totalsum:15:6);
   end;
 
   if normalize then begin
      writeln(drxyin,space:10,totalsum/maxfactor:15:6);
   end
   else begin
      writeln(drxyin,space:10,totalsum:15:6);
   end;
 
end;
(* end module zhang3ribl *)


{ comment out this old unused method:
(* begin module zhang4ribl *)
(* calculates the zhang 4D distance *)
procedure zhang4ribl(matrixA, matrixB: ribltypeptr; var matrixC: ribltypeptr;
                  var posdiff, drxyin: text; space: integer;
                  rangebool: boolean; userfrom, userto: integer);
(* This is the old computation from the Ribl matrix in four
dimensions. *) 
var
   x: integer; (* counter *)
   sum: real; (* sum of differences at position *)
   totalsum: real; (*the sum of all the differences for all the positions *)
   newfrom: integer; (* the new from range value of the overlap *)
   newto: integer; (* the new to range value of the overlap *)
begin
 
   totalsum:= 0;
 
   if space <= 0 then begin
      if rangebool = false then begin
         newfrom:= matrixA^.frombase;
         newto:=  matrixA^.tobase + space;
      end;
 
      if rangebool then begin
         newfrom:= userfrom;
         newto:= userto + space;
      end;
 
   end;

   if space > 0 then begin
      if rangebool = false then begin
         newfrom:= matrixA^.frombase + space;
         newto:=  matrixA^.tobase;
      end;

      if rangebool then begin
         newfrom:= userfrom + space;
         newto:= userto;
      end;

   end;

   for x:= (newfrom) to (newto) do begin
       matrixA^.data[a,x]:= exp((matrixA^.data[a,x] - 2)*ln(2));
       matrixA^.data[c,x]:= exp((matrixA^.data[c,x] - 2)*ln(2)); 
       matrixA^.data[g,x]:= exp((matrixA^.data[g,x] - 2)*ln(2));
       matrixA^.data[t,x]:= exp((matrixA^.data[t,x] - 2)*ln(2));
       matrixB^.data[a,x]:= exp((matrixB^.data[a,x] - 2)*ln(2));
       matrixB^.data[c,x]:= exp((matrixB^.data[c,x] - 2)*ln(2));
       matrixB^.data[g,x]:= exp((matrixB^.data[g,x] - 2)*ln(2));
       matrixB^.data[t,x]:= exp((matrixB^.data[t,x] - 2)*ln(2));
       matrixC^.data[a,x]:= sqr(matrixA^.data[a,x] - matrixB^.data[a,x - space]);
       matrixC^.data[c,x]:= sqr(matrixA^.data[c,x] - matrixB^.data[c,x - space]);
       matrixC^.data[g,x]:= sqr(matrixA^.data[g,x] - matrixB^.data[g,x - space]); 
       matrixC^.data[t,x]:= sqr(matrixA^.data[t,x] - matrixB^.data[t,x - space]);
       sum:= sqrt(3*(matrixC^.data[a,x] + matrixC^.data[c,x] + matrixC^.data[g,x]
                   + matrixC^.data[t,x]));
       sum:= sum/2;
       write(posdiff,x:10);
       writeln(posdiff,sum:15:6);

       totalsum:= totalsum + sum;
   end;

   writeln(posdiff, '*');
   writeln(posdiff,'* the sum of all the Euclidean distances is: ',
                      totalsum:15:6, ' bits');

   writeln(drxyin,space:10,totalsum:15:6); end;

(* end module zhang4ribl *)
}


(* begin module dlapribl *)
(* calculates the difference for only the overlapping matrices *)
procedure dlapribl(matrixA, matrixB: ribltypeptr; var matrixC: ribltypeptr;
                  var posdiff, drxyin: text; space: integer;
                  rangebool: boolean; userfrom, userto: integer);
 
var
   x: integer; (* counter *)
   sum: real; (* sum of differences at position *)
   totalsum: real; (*the sum of all the differences for all the positions *)
   newfrom: integer; (* the new from range value of the overlap *)
   newto: integer; (* the new to range value of the overlap *)

begin
 
   totalsum:= 0;
 
   if space <= 0 then begin
      if rangebool = false then begin
         newfrom:= matrixA^.frombase;
         newto:=  matrixA^.tobase + space;
      end;
 
      if rangebool then begin
         newfrom:= userfrom;
         newto:= userto + space;
      end;
 
   end;
 
   if space > 0 then begin
      if rangebool = false then begin
         newfrom:= matrixA^.frombase + space;
         newto:=  matrixA^.tobase;
      end;
 
      if rangebool then begin
         newfrom:= userfrom + space;
         newto:= userto;
      end;
 
   end;
 

   for x:= (newfrom) to (newto) do begin
      matrixC^.data[a,x]:= (matrixA^.data[a,x] - matrixB^.data[a,x - space]);
      matrixC^.data[c,x]:= (matrixA^.data[c,x] - matrixB^.data[c,x - space]);
      matrixC^.data[g,x]:= (matrixA^.data[g,x] - matrixB^.data[g,x - space]);
      matrixC^.data[t,x]:= (matrixA^.data[t,x] - matrixB^.data[t,x - space]);
      sum:= (matrixC^.data[a,x] + matrixC^.data[c,x] + matrixC^.data[g,x]
      + matrixC^.data[t,x]);
      write(posdiff,x:10);
      writeln(posdiff,sum:15:6);
 
      totalsum:= totalsum + sum;
   end;
   writeln(posdiff, '*');
   writeln(posdiff,'* the sum of all the Euclidean distances is: ', totalsum:15:6);
 
   writeln(drxyin,space:10,totalsum:15:6);
end;
 
(* end module dlapribl *)

(* begin module readdrparam *)
(* read in all the parameter files *)
procedure readdrparam(var diffriblp: text;
                      var parameterversion: real;
                      var rangechar, spacechar: char;
                      var userfrom, userto, spacefrom, spaceto: integer;
                      var rangebool, spacebool: boolean;
                      var calctype: char);
(* read diffribl (dr) parameters. *)
begin

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

   read(diffriblp,rangechar); 
   rangebool:= false;
   spacebool:= false;
 
   if rangechar = 'r' then begin
      readln(diffriblp);   
   end;  
 
   if rangechar = 'u' then begin
      rangebool:= true;
      read(diffriblp, userfrom); 
      readln(diffriblp, userto);
   end;  
 
   if (rangechar <> 'r') and (rangechar <> 'u') then begin 
      writeln(output, 'the range needs to specified with an "r" or "u"');
      halt;
   end;

   read(diffriblp, spacechar);

   if spacechar = 'n' then begin
      readln(diffriblp);
   end;

   if spacechar = 'v' then begin
      spacebool:= true;
      read(diffriblp, spacefrom);
      readln(diffriblp, spaceto);
   end;

   if (spacechar <> 'n') and (spacechar <> 'v') then begin  
      writeln(output, 'the variability range needs to specified',
                      ' with an "n" or "v"'); 
      halt; 
   end;

   read(diffriblp, calctype);

   readln(diffriblp);

   if not (calctype in ['e','o','d','s','z','y']) then begin
      writeln(output, 'calculation type must be specified',
                      ' by one of "eodszy"');
      halt;
   end;

end;
(* end module readdrparam *)

(* begin module writedrparam *)
procedure writedrparam(var f: text;
                       parameterversion: real;
                       rangechar, spacechar: char;
                       userfrom, userto, spacefrom, spaceto: integer;
                       rangebool, spacebool: boolean;
                       calctype: char);
(* write the diffribl (dr) parameters to file f *)
begin
   (* for now only report one more parameter *)
   write(f, '* ', calctype,' ');

   case calctype of
     'e': writeln(f, 'Euclidean distance');
     'o': writeln(f, 'Euclidean distance, 4 dimensional');
     'd': writeln(f, 'Difference distance');
     's': writeln(f, 'scan average');
     'z': writeln(f, 'Zhang, position-by-position');
     'Z': writeln(f, 'Zhang 4 dimensional');
     'y': writeln(f, 'Zhang, position-by-position, normalized by sqrt(3/2)L');
   end;

end;
(* end module writedrparam *)

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

(* begin module diffribl.themain *)
(* the main procedure of the program *)
procedure themain(var ribla, riblb, diffriblp, posdiff, drxyin: text);
(* the main procedure of the program *)

var
   parameterversion: real; (* parameter version number *)
   rangechar: char; (* this char tells whether to use from to range or not *)
   spacechar: char; (* switches on variable overlap test *)
   userfrom: integer; (* the user specified from range *)
   userto: integer; (* the user specified to range *)
   spacefrom: integer; (* the from range for the spacing function *)
   spaceto: integer; (* the to range for the spacing function *)

   rangebool: boolean; (* if true, use user's range *)
   spacebool: boolean; (* if true, use user's variability range *)

   calctype: char; (* The type of calculation *)

   ribltypeptrA: ribltypeptr; (* ribla *)
   ribltypeptrB: ribltypeptr; (* riblb *)
   ribltypeptrC: ribltypeptr; (* new ribl that contains dist or diff values *)

   space: integer; (* distance from zero position *)
   
begin
   writeln(output,'diffribl',version:4:2);

   patentmessage(output);

   (* reads in all of the parameters *)
   readdrparam(diffriblp,
               parameterversion, rangechar, spacechar, userfrom, userto,
               spacefrom, spaceto, rangebool, spacebool,
               calctype);

   reset(ribla);
   reset(riblb);

   (* reads in the two matrixes *)
   getriblmatrix(ribla, ribltypeptrA);
   getriblmatrix(riblb, ribltypeptrB);

   (* makes a fake matrices to be filled with distance or difference values *)
   reset(ribla);
   getriblmatrix(ribla, ribltypeptrC);

   if rangechar = 'r' then begin
      (* take the range from the matrix (which one?????) *)
      userfrom := ribltypeptrA^.frombase;
      userto   := ribltypeptrA^.tobase;
      writeln(output,'range is taken from the A matrix as:',
                     ' ', userfrom:1,
                     ' ', userto:1);
   end;

   (* check the userfrom / userto range *)
   if rangebool then begin
      if (userfrom < ribltypeptrA^.frombase) or
                   (userto > ribltypeptrA^.tobase) then begin
         writeln(output,'userfrom userto range must be smaller than ribl size');
         halt;
       end;
    end;

   (* places zeroes in positions outside of the specified matrix *)
   if rangebool = false then begin
      cleanmatrix(ribltypeptrA, ribltypeptrA^.frombase, ribltypeptrA^.tobase);
      cleanmatrix(ribltypeptrB, ribltypeptrB^.frombase, ribltypeptrB^.tobase);
   end;

   if rangebool then begin
      cleanmatrix(ribltypeptrA, userfrom, userto);
      cleanmatrix(ribltypeptrB, userfrom, userto);
   end;

   rewrite(drxyin);

   rewrite(posdiff);
   writeln(posdiff,'* diffribl ',version:4:2);

   writedrparam(posdiff, 
                parameterversion, rangechar, spacechar, userfrom, userto,
                spacefrom, spaceto, rangebool, spacebool,
                calctype);

   write(posdiff, '* Matrix A: ');
   writestring(posdiff,ribltypeptrA^.riblname);
   writeln(posdiff);

   write(posdiff, '* Matrix B: ');
   writestring(posdiff,ribltypeptrB^.riblname);
   writeln(posdiff);

   writeln(posdiff, '*');
   writeln(posdiff,'* column 1: position');
   writeln(posdiff,'* column 2: square root of the sum of squared',
                  ' differences or distances');
   writeln(posdiff, '*');


   (* do the calculations *)
   if spacebool = false then begin
      spacefrom := 0;
      spaceto   := 0;
   end;

   for space := spacefrom to spaceto do begin

      case calctype of
        'e': (* Euclidean distance *)
            overlapribl(ribltypeptrA, ribltypeptrB, ribltypeptrC,
                        posdiff, drxyin,
                        space, rangebool, userfrom, userto);

        'o': (* Euclidean distance, 4 dimensional *)
            onepointribl(ribltypeptrA, ribltypeptrB, ribltypeptrC,
                         posdiff, drxyin, space, 
                         rangebool, userfrom, userto);

        'd': (* Difference distance *)
            dlapribl(ribltypeptrA, ribltypeptrB, ribltypeptrC,
                     posdiff, drxyin,
                     space, rangebool, userfrom, userto);

        's': (* scan average *)
             scanribl(ribltypeptrA, ribltypeptrB, ribltypeptrC,
                      posdiff, drxyin,
                      space, rangebool, userfrom, userto);

        'z': (* Zhang - corresponds to Euclidean distance *)
             zhang3ribl(ribltypeptrA, ribltypeptrB, ribltypeptrC,
                      posdiff, drxyin,
                      space, rangebool, userfrom, userto,
                      false);

        'y': (* Zhang - "z" normalized *)
             zhang3ribl(ribltypeptrA, ribltypeptrB, ribltypeptrC,
                      posdiff, drxyin,
                      space, rangebool, userfrom, userto,
                      true);

{ Block this old method:
        (Z) The sixth, specified by "Z", is the direct measurement of
            four dimentional distance between the two matrices.  This
            is done with the following equations:
.
            First, the probabilities (Pi) are use as in (z). 
            Positional distance is calculated according to Zhang's
            method:
.
            positional distance = sqrt(3)/2 * sqrt[(a2-a1)^2 + (c2-c1)^2
                                                 + (g2-g1)^2 + (t2-t1)^2]
.
            This positional distance is then summed for all positions
            giving the total sum of positional distances.
            It should give the same result as z.
(except that z is now directly from the probabilities)
        'Z': (* Zhang 4D - corresponds to Euclidean 4 dimensional *)
             zhang4ribl(ribltypeptrA, ribltypeptrB, ribltypeptrC,
                       posdiff, drxyin,
                       space, rangebool, userfrom, userto);
}

      end;
   end;

end;
(* end module diffribl.themain *)

begin
   themain(ribla, riblb, diffriblp, posdiff, drxyin);
1: end.
