program frame(test,norm,result,output);
(* this program takes vectors from the file 'test' and does dotproducts to
   evaluate the potential reading frames of a sequence.
   by gary d. stormo
module libraries needed: delman, delmods, matmods *)

label 1;  (* end of program *)

const
(* begin module version *)
version = 1.17;  (* of frames 1994 sep 5
origin: july 28, 1983 *)
(* end module version *)

(* begin module describe.frame *)
(*
name
      frame: evaluator of potential reading frames

synopsis
      frame(test: in, norm: in, result: out, output: out)

files
      test: encoded vectors of the sequences to be tested for reading frames
      norm: encoded vector of the sequences used as the standard for testing
      result: the results of the tests; each sequence from test is evaluated
         for each of the three possible reading frames
      output: for messages to the user

description
      this calculates correlation coefficients between the standard and each
      of the three possible frames of the test sequences.  the sequences must
      be encoded so that each of the oligos (of whatever length is desired)
      are counted in each of the three frames.

examples
      the files framet and framen are examples of test and norm

documentation
      delman.use.frame

see also
      encode.p

author
      gary d. stormo

bugs
      none known
*)
(* end module describe.frame *)

(* begin module vector.const *)
      maxvecpart = 64; (* the number of elements in a 'part' of a vector *)
      vpagewidth = 64; (* the width (in characters) of the output vector
         file from procedure writevector *)
(* end module vector.const version = 'matmod 1.95 85 apr 18 tds/gds'; *)

type
(* begin module encode.type.param *)
(* the following types allow the user to specify parameters which will
   be used to encode the sequences in the book. *)
(* spaces are allowed between the encoded bases, and the number of
   bases to be skipped between each encoded pair of bases are kept
   in this linked list of integers *)
      spaceptr = ^spacelist;
      spacelist = record
         skips : integer; (* bases skipped to next coded base *)
         next : spaceptr; (* points to next spacing number *)
      end;

(* the encoding parameters for each region are stored in these records,
   and the records for each region are connected into a linked list of
   all the encoding parameters for the entire sequences *)
      endpoints = (start,stop);  (* end points of a coding region *)
      paramptr = ^parameter;
      parameter = record  (* these are the instructions for doing the coding *)
         range : array[endpoints] of integer; (* the bases to be coded by these
                                 instructions, relative to the alignedbase *)
         window : integer; (* the number of bases included in a coding vector *)
         wshift : integer;  (* movement of the window to the next site *)
         coding : integer; (* the 'level' at which the coding is being done,
                             i.e., 1: mononucleotides; 2: dinucleotides; ... *)
         spaces : spaceptr;  (* the spacing between the encoded bases *)
         cshift : integer;  (* shift to the next coding unit in the window *)

         (* values calculated at read-in time *)
         wvlength: integer; (* length of a window vector.  this is coding
            raised to the 4th power *)
         pvlength: integer; (* parameter vector length.  the vector made up
            of a series of window vectors. *)
         (* note: pvlength/wvlength is the number of windows in the
            parameter range *)

         next : paramptr;  (* set of instructions for coding another region *)
      end;

(* end module encode.type.param version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.type *)
(* vectors are useful for many applications, including the encoding of
   sequences.  this vector type is designed to be flexible enough to be
   used whenever one needs a vector.  it is designed as a linked list so
   it can contain as many elements as are ever needed.  the actual 'vector'
   is a record containing the length and a pointer to the first 'vectorpart'.
   that vectorpart is a record containing an array of the first 'maxvecpart'
   elements  (maxvecpart is a constant, from module vector.const, which must
   be included) and a pointer to the next 'vectorpart'.  the elements are
   of type real so that either integers or reals may be used.  *)
      partptr = ^vectorpart;
      vectorpart = record
         numbers : array[1..maxvecpart] of real;
         next : partptr;
      end;

      vector = record
         length : integer;
         part : partptr;
      end;
(* end module vector.type version = 'matmod 1.95 85 apr 18 tds/gds'; *)

var
      norm,  (* the standard vector for comparison to *)
      test,  (* the encoded vectors for testing the frames *)
      result: text; (* the results of the testing *)
      normparam,  (* the parameter record from file norm *)
      testparam: paramptr; (* the parameter record from file test *)
      normvector,  (* the vector from file norm *)
      testvector: vector; (* successively holds each vector from file test *)
      r,  (* the number of independently coded regions *)
      e: integer;  (* the end-of-sequence symbol *)

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

(* begin module unlimitln *)
procedure unlimitln(var afile: text);
(* this procedure removes a stupid system dependent limit on the number of
lines that one can write to a file.  you may remove it from the code if
your system does not want or need this.  suggested method: place comments
around the contents of the procedure. *)
begin
      linelimit(afile, maxint); (* set 'infinite' lines allowed for afile *)
end;
(* end module unlimitln version = 'delmod 6.51 85 apr 17 tds/gds' *)

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

(* begin module copylines *)
function copylines(var fin, fout: text; n: integer): integer;
(* copy n lines of file fin to file fout.
the actual number of lines copied is returned. *)
var
      index: integer; (* the current line number *)
begin (* copylines *)
      index := 0;
      while (not eof(fin)) and (index < n) do begin
         copyaline(fin, fout);
         index := succ(index)
      end;

      copylines := index
end; (* copylines *)
(* end module copylines version = 'delmod 6.51 85 apr 17 tds/gds' *)

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

(* begin module package.encode *)
(* *********************************************************************** *)
(* begin module encode.readencpar *)
(* read encode parameters *)

function vectorsize(param: paramptr): integer;
(* determine the size of the vector generated by a string of parameter
   records, begin with 'param' *)
var   size: integer;  (* of the vector, so far *)
begin
      size := 0;
      while (param <> nil) do begin
         size := size + param^.pvlength;
         param := param^.next;
      end;
      vectorsize := size;
end;

function paramsize(param: paramptr): integer;
(* determine the size of the vector generated by one parameter record *)
var   rangesize,  (* the number of bases in the range *)
      numwindows: integer; (* the number of windows which start in the
                  range, taking into account the shifting *)
begin
      with param^ do begin
         rangesize := range[stop] - range[start] + 1;
         numwindows := trunc((rangesize -1)/wshift + 1);
         paramsize := numwindows * wvlength
      end
end; (* paramsize *)


procedure readencpar(var thefile: text; var param: paramptr;
                    var regions, vectorlength: integer);
(* read the parameter list from 'thefile' and put the information into the
   linked parameter list, beginning with 'param'; the number of parameter
   records is returned in 'regions'.  the length of each vector is
   vectorlength *)
var   aparam,  (* from the parameter list *)
      newparam: paramptr;  (* the new parameter record from 'thefile' *)
      firstspaces,  (* in the spaces list *)
      newspaces: spaceptr; (* for the spaces list *)
      i,j: integer;  (* indices *)
      ch: char;  (* to check for proper positioning *)
begin
      reset(thefile);
      if eof(thefile) then begin
         writeln(output,' encoded sequence file is empty');
         halt
      end;
      new(aparam);
      if param = nil then new(param);
      aparam := param;
            (* skip the header *)
      for i := 1 to 2 do readln(thefile);

            (* skip a blank line in the header *)
      readln(thefile);

            (* read the number of regions *)
      readln(thefile,regions);

      for i := 1 to regions do begin
         with aparam^ do begin
               (* read the range *)
            read(thefile,range[start]);
            repeat read(thefile,ch) until (ch = 'o');
            readln(thefile,range[stop]);
               (* read the window size and shift *)
            readln(thefile,window);
            readln(thefile,wshift);
               (* read the coding level and arrangement *)
            read(thefile,coding);
            new(spaces);
            if (coding > 1) then begin
               repeat read(thefile,ch) until (ch = ':');
               new(firstspaces);
               spaces := firstspaces;
               read(thefile,spaces^.skips);
               for j := 2 to coding -1 do begin
                  new(newspaces);
                  spaces^.next := newspaces;
                  spaces := newspaces;
                  read(thefile,spaces^.skips);
               end;
               spaces^.next := nil;
               spaces := firstspaces
            end
            else spaces := nil;
            readln(thefile);
            readln(thefile,cshift);

            (* set up wvlength *)
            wvlength := round(exp(coding*ln(4)));
         end;

         (* set up pvlength *)
         aparam ^.pvlength := paramsize(aparam);

         if (i < regions) then begin
            new(newparam);
            aparam^.next := newparam;
            aparam := newparam;
         end;
      end;
      aparam^.next := nil;

      (* read the vector size *)
      readln(thefile,vectorlength);
      if vectorlength <> vectorsize(param)
      then begin
         writeln(output,' vector length in the encoded file');
         writeln(output,' does not correspond to the parameters');
         halt
      end;
end;

(* end module encode.readencpar version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module encode.functions *)
(* these functions are quite handy to help one know where one is
in an encoded vector.  they all take a parameter list as read by
readencpar, and the number of the element of the vector one is interested
in.  this element value runs from 0 to one less than the length of the
parameter vector.  examples are in evprint. *)

function beginpv(aparam: paramptr; element: integer): boolean;
(* are we at the beginning of a parameter vector? *)
begin (* beginpv *)
      if aparam <> nil (* is there a vector? *)
      then beginpv := (element = 0)
      else beginpv := true (* note circular use of aparam in evprint... *)
end; (* beginpv *)

function endpv(aparam: paramptr; element: integer): boolean;
(* are we at the end of a parameter vector?
the zero element is the first element of the parameter vector. *)
begin (* endpv *)
      if aparam <> nil
      then endpv := (element >= aparam^.pvlength)
      else endpv := true (* past the end of all vectors *)
end; (* endpv *)

function beginwv(aparam: paramptr; element: integer): boolean;
(* are we at the beginning of a window vector in the current parameter
vector?  the zero element is the first element of the parameter vector. *)
begin (* beginwv *)
      if aparam <> nil
      then with aparam^ do beginwv := ((element mod wvlength) = 0)
      else beginwv := false (* no further parameters *)
end; (* beginwv *)

function endwv(aparam: paramptr; element: integer): boolean;
(* are we at the end of a window vector in the current parameter vector?
the zero element is the first element of the parameter vector. *)
begin (* endwv *)
      if aparam <> nil
      then with aparam^ do endwv := (((element + 1) mod wvlength) = 0)
      else endwv := true (* no further parameters *)
(*;writeln(output,'element=',element:1,'   wvlength=',aparam^.wvlength:1); *)
end; (* endwv *)

function encposition(aparam: paramptr; element: integer): integer;
(* this returns the aligned position corresponding to the element in
a parmaeter list (aparam).
the zero element is the first element of the parameter vector. *)
begin (* encposition *)
      with aparam^
         do encposition := range[start] + (element div wvlength) * wshift
end; (* encposition *)
(* end module encode.functions version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module encode.evprint *)
(* procedures for encoded vector printing *)

procedure evln(var afile: text; element: integer; curparam: paramptr);
(* see evprint for definitions.  this procedure is used to provide
writelns to afile when the user wants to control the display better. *)
begin (* evln *)
      (* for end of windows: *)
      if beginwv(curparam,element)
      then writeln(afile)
      else if curparam = nil then writeln(afile); (* end of vector *)

      (* for end of parameters: *)
      if beginpv(curparam,element) then writeln(afile)
end; (* evln *)

procedure evheader(var afile: text;
                       posfield, valfield: integer;
                       firstparam: paramptr;
                   var curparam: paramptr);
(* see evprint for definitions of the variables and how to use evheader.
evheader must be called if you want to provide your own writelns
to the evprint display, ie, when ln is false *)
const
      codingmax = 4; (* the maximum coding level allowed *)
var
      bases: array[1..codingmax] of char; (* the bases *)
      i,j,n: integer; (* indices for loops *)
      size: integer; (* the number of characters in the encoded object.
                        a=1, aa=2, axa=3, etc *)
      s: spaceptr; (* for finding size *)
procedure space(i:integer);
(* put out i spaces to afile *)
var
      index: integer; (* for i *)
begin (* space *)
      for index := 1 to i do write(afile,' ');
end; (* space *)
procedure advance(level: integer);
(* advance the bases at the level.  aa -> ac -> ag -> ... -> tt *)
begin (* advance *)
      case bases[level] of
         'a': bases[level] := 'c';
         'c': bases[level] := 'g';
         'g': bases[level] := 't';
         't': begin
            bases[level] := 'a';
            advance(level - 1) (* aha.. *)
         end
      end
end; (* advance *)
begin (* evheader *)
      if curparam = nil then curparam := firstparam; (* start or recycle *)

      with curparam^ do begin
         if coding > codingmax then begin
            writeln(output,' evheader: coding levels must not be higher',
                           ' than ', codingmax:1);
            halt
         end;

         (* find the size of the encoded object *)
         size := coding;
         s := curparam^.spaces;
         while s <> nil do begin
            size := size + s^.skips;
            s := s^.next
         end;

         (* see if there is room for printing *)
         if abs(valfield) < size then begin
            writeln(output,' evheader: valfield (',valfield:1,
                           ') cannot be less than coding level and spaces (',
                           size:1,')');
            halt
         end;

         (* clear bases *)
         for i := 1 to coding do bases[i] := 'a';

         (* print the header *)
         (* room for space, ' to ' and two numbers *)
         space(1 + 4 + 2*posfield);
         for i := 1 to wvlength do begin
            space(1 + abs(valfield) - size);
            s := spaces;
            for n := 1 to coding do begin
               write(afile, bases[n]);
               if s <> nil then begin
                  for j := 1 to s^.skips
                     do write(afile,'x');
                  s := s^.next
               end
            end;

            (* advance the bases except for the last time *)
            if i <> wvlength then advance(coding)
         end
      end
end; (* evheader *)

procedure evprint(var afile: text;         (* the file to print to *)
                      posfield: integer;   (* the number of characters in
                                              which to print the position *)
                      val: real;           (* the value to print *)
                      valfield,            (* the number of characters in
                                              in which to print the value *)
                      valdecimal: integer; (* the number of characters
                                              devoted to the digits below the
                                              decimal point of val.  if <= 0
                                              val is printed as an integer *)
                      ln: boolean;         (* control writelns (see below) *)
                      firstparam: paramptr; (* the first parameter *)
                  var curparam: paramptr;  (* the current encoding parameter *)
                  var element: integer);   (* the position in the vector *)
(* evprint is an encoded vector printing routine.  it is designed to make
generation of display of encoded vectors or calculated vectors based on
encoded vectors easy.  the basic idea is that evprint will print all the
display associated with one vector value.  so at the beginning it prints
an appropriate header from the vector parameters, and for each line it
prints the beginning and end positions of the window.
1) all fields printed out are separated by one space.
2) the value arguments are to be set by the user before calling evprint.
evprint is then called sequentially with all the vector values placed in val.
this allows the data to come from any data structure; it is not forced
to come from a particular array or file type.
3) to initialize, read in firstparam, set the element to 0 and curparam to nil.
evprint will restart at firstparam after each pass across the vector.
4) once the printing of the display has begun,
do not modify the element or curparam yourself.  it is also not wise to
modify any value parameters except val and possibly valdecimal.  the element
argument is used internally by evprint to keep track of the display.
the value arguments, after being set by the user initially, control the
form of the display.

control of writelns:  if ln is true, then evprint will take care of
writelns.  otherwise the user must do the writeln, as:
      if beginpv(curparam,element) then begin
         evheader(afile,posfield,valfield,firstparam,curparam);
         writeln(afile,[your headers for each column]);
      end;
      evprint(afile,...false...)
      if beginwv(curparam,element)
         then write(afile,[your own material on the same line]);
      evln(afile,curparam,element);
note that evprint had moved the element to the next position, so the
test is for begin of wv, not end of wv.

if valfield is negative (but abs(valfield) is a reasonable size) then
val will not be printed, allowing the user to do so after evprint.
this allows more flexible displays.
note: ln must be false to use this feature.
 *)
procedure sider;
(* produce the side of the display *)
var
      position: integer; (* the location that a window begins
                            from the aligned base *)
begin (* sider *)
      with curparam^ do begin
         position :=encposition(curparam, element);
         write(afile,' ',position: posfield, ' to ',
                 (position + window - 1):posfield)
      end
end; (* sider *)
procedure middler;
(* produce the middle of the display *)
begin (* middler *)
      write(afile,' ');
      if valfield > 0  (* to display or not to display... *)
      then if valdecimal <= 0 (* chose display *)
           then write(afile, round(val):valfield)
           else write(afile, val:valfield:valdecimal)
      else if ln then begin
         writeln(output,' evprint: if valfield is negative,',
                        ' then ln must be false');
         halt
      end
end; (* middler *)
procedure step;
(* step to the next position *)
begin (* step *)
      element := succ(element); (* move along the vector *)

      if endpv(curparam, element) (* are we off the end of the vector? *)
      then begin (* yes *)
         curparam := curparam^.next;
         element := 0;
      end
end; (* step *)
begin (* evprint *)
      if ln then if beginpv(curparam,element) then begin
         evheader(afile, posfield, valfield, firstparam, curparam);
         writeln(afile)
      end;

      if beginwv(curparam, element) then sider;

      middler;

      step;

      if ln then evln(afile,element,curparam)
end; (* evprint *)
(* end module encode.evprint version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* *********************************************************************** *)
(* end module package.encode version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.functions *)
(* *********************************************************************** *)
(* begin module vector.io *)
(* *********************************************************************** *)
(* begin module vector.primitives *)
(* these functions and procedures do manipulations on vectors. *)

function vget(v:vector; pos:integer): real;
(* this returns from vector 'v' the value of the element at position 'pos' *)
var i: integer;
begin
      if (pos > v.length) or (pos < 1) then begin
         writeln(output,' error in call to function vget:',
            ' position ',pos:1,' is beyond the end of the vector');
         halt;
      end;

      (* move to the correct 'vectorpart' *)
      for i := 1 to ((pos -1) div maxvecpart) do  v.part := v.part^.next;

      (* get the proper array element from this part *)
      vget := v.part^.numbers[(pos -1) mod maxvecpart + 1];
end;

procedure vput(var v:vector; pos:integer; number: real);
(* this puts into vector 'v' the value of 'number' at position 'pos' *)
var i: integer;
      firstpart: partptr;  (* of the vector *)
begin
      if (pos > v.length) or (pos < 1) then begin
         writeln(output,' error in call to function vput:',
            ' position ',pos:1,' is beyond the end of the vector');
         halt;
      end;

      (* move to the correct 'vectorpart' *)
      firstpart := v.part;
      for i := 1 to ((pos -1) div maxvecpart) do  v.part := v.part^.next;

      (* put the 'number' into the proper array element for this part *)
      v.part^.numbers[(pos -1) mod maxvecpart + 1] := number;
      v.part := firstpart;
end;

procedure makevector(var v: vector; l: integer);
(* create the linked list of vector-parts which are needed for a vector
   of length 'l' *)
var   numparts,  (* number of parts needed for the vector *)
      i: integer; (* index to the parts of the vector *)
      firstpart,  (* of the vector *)
      newpart: partptr;  (* of the vector *)
begin
      if l < 1 then begin
         writeln(output,' makevector: positive length required');
         halt
      end;

      v.length := l;
      new(v.part);
      firstpart := v.part;
      numparts := (v.length - 1) div  maxvecpart + 1;
      for i := 1 to (numparts - 1) do begin
         new(newpart);
         v.part^.next := newpart;
         v.part := newpart;
      end;
      v.part^.next := nil;
      v.part := firstpart;
end;
(* end module vector.primitives version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.readvector *)
procedure readvector(var thefile:text; var v: vector);
(* read the elements of v into v from 'thefile'.  v must already be set
   up (all the parts created and linked together) before calling.
   'thefile' must contain, from the cursor point until the end of the
    vector, only numbers, either integers or reals; otherwise it will
    bomb.  the vector ends with an end of line so that
    end of file can be tested for. *)
var i,j: integer;  (* indecies *)
      numparts,  (* number of parts in the vector *)
      lastpart: integer; (* the number of elements in the last vector part *)
      firstpart: partptr;  (* pointer to the first vector part *)
begin
      firstpart := v.part;
      numparts := (v.length - 1) div maxvecpart + 1;
      lastpart := (v.length - 1) mod maxvecpart + 1;
      for i := 1 to (numparts - 1) do begin
         for j := 1 to maxvecpart do read(thefile,v.part^.numbers[j]);
         v.part := v.part^.next;
      end;
      for j := 1 to lastpart do read(thefile,v.part^.numbers[j]);
      readln(thefile);
      v.part := firstpart;
end;
(* end module vector.readvector version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.writevector *)
procedure writevector(var thefile: text; v: vector; y,z: integer);
(* writes the elements of 'v' to 'thefile'.  the integers y and z are:
    y: the size of the field for printing the number;
    z: the size of the decimal part of the field,
    if z = 0 then integers will be printed. *)
var   pos,   (* posititon in the vector *)
      i,j: integer;  (* indecies *)
      numparts,  (* number of parts in the vector *)
      lastpart: integer; (* the number of elements in the last vector part *)
      firstpart: partptr;  (* pointer to the first vector part *)
      x: integer; (* the number of elements to write on a line *)
begin
      x := trunc(vpagewidth/(y+1));
      pos := 0;
      firstpart := v.part;
      numparts := (v.length - 1) div maxvecpart + 1;
      lastpart := (v.length - 1) mod maxvecpart + 1;
      if (z = 0) then begin
         for i := 1 to (numparts - 1) do begin
            for j := 1 to maxvecpart do begin
               write(thefile,' ',round(v.part^.numbers[j]):y);
               pos := succ(pos);
               if (pos mod x = 0) then writeln(thefile);
            end;
            v.part := v.part^.next;
         end;
         for j := 1 to lastpart do begin
            write(thefile,' ',round(v.part^.numbers[j]):y);
            pos := succ(pos);
            if ((pos < v.length) and (pos mod x = 0)) then writeln(thefile);
         end
      end
      else begin
         for i := 1 to (numparts - 1) do begin
            for j := 1 to maxvecpart do begin
               write(thefile,' ',v.part^.numbers[j]:y:z);
               pos := succ(pos);
               if (pos mod x = 0) then writeln(thefile);
            end;
            v.part := v.part^.next;
         end;
         for j := 1 to lastpart do begin
            write(thefile,' ',v.part^.numbers[j]:y:z);
            pos := succ(pos);
            if ((pos < v.length) and (pos mod x = 0)) then writeln(thefile);
         end;
      end;
      writeln(thefile);
      v.part := firstpart;
end;
(* end module vector.writevector version = 'matmod 1.95 85 apr 18 tds/gds'; *)
(* *********************************************************************** *)
(* end module vector.io version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.dotproduct *)
function dotproduct(vectora,vectorb: vector): real;
(* this returns the dotproduct of 'vectora' and 'vectorb' *)
var i: integer;  (* an index *)
    j: real;  (* partial products *)
begin
      if vectora.length <> vectorb.length then begin
         writeln(output,' function dotproduct: vector lengths must',
            ' be equal');
         halt
      end;
      j := 0;
      for i := 1 to vectora.length do
         j := j + vget(vectora,i) * vget(vectorb,i);
      dotproduct := j;
end;
(* end module vector.dotproduct version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.magnitude *)
function magnitude(var v: vector): real;
(* find the magnitude (length) of the vector v *)
begin
      magnitude := sqrt(dotproduct(v,v));
end;
(* end module vector.magnitude version = 'matmod 1.95 85 apr 18 tds/gds'; *)

(* begin module vector.normalize *)
procedure normalize(var v: vector);
(* this replaces the vector v with the congruent vector of unit length,
   i.e.,  the resulting vector v has the property that v.v = 1   *)
var i: integer;  (* index *)
      length: real;  (* of the unnormalized vector *)
begin
      length := magnitude(v);
      for i := 1 to v.length do
         vput(v,i,vget(v,i)/length);
end;
(* end module vector.normalize version = 'matmod 1.95 85 apr 18 tds/gds'; *)
(* *********************************************************************** *)
(* end module vector.functions version = 'matmod 1.95 85 apr 18 tds/gds'; *)

procedure replacevector(var v: vector; n: integer);
(* each element of v is replaced by the difference between the element and
   the average value for all the elements of the same oligo.  'n' is the
   number of oligos for each frame (i.e., n = 4 if we are looking at monos,
   n = 16 if we are looking at dis, ...).
   for example, if we are looking at di-nucleotides we would replace the
   number of aa"s in each frame (vector positions 1, 17, 33) by their
   differences from the average value for aa"s.  the same is done for each
   oligo.  the resulting vector has the property that the sum of its elements
   is zero.  this vector is then normalized so that it also has the property
   that its length is one.  taking the dotproduct between two vectors which
   each have these properties is equivalent to calculating the correlation
   coefficient between the arrays which generated the vectors. *)
var  sum: real;  (* of the elements for a particular base *)
      i: integer;  (* an index *)
begin
      for i := 1 to n do begin
         sum := vget(v,i) + vget(v,i+n) + vget(v,i+2*n);
         vput(v,i,vget(v,i) - sum/3);
         vput(v,i+n,vget(v,i+n) - sum/3);
         vput(v,i+2*n,vget(v,i+2*n) - sum/3);
      end;
      normalize(v);
end;

procedure rearrange(var v: vector; n: integer);
(* rearrange the elements of v so as to have each oligo composition
   represent the in-frame case.  this is done by rotating the three parts of
   v into each of the three thirds of the vector.  n is the number of
   oligos in each part (i.e., n = 4 if we are looking at monos...) *)
var i: integer;
    r: real;
begin
      for i := 1 to n do begin
         r := vget(v,i);
         vput(v,i,vget(v,i+n));
         vput(v,i+n,vget(v,i+2*n));
         vput(v,i+2*n,r);
      end;
end;

procedure evaluate;
(* write the result of the dotproduct of the norm and test vectors for
   each of the three frames of the test vector *)
var i: integer;
begin
      write(result,' ':4,
            dotproduct(normvector,testvector):7:4);
      for i := 1 to 2 do begin
         rearrange(testvector,testparam^.wvlength);
         write(result,' ':4,
               dotproduct(normvector,testvector):7:4);
      end;
      writeln(result);
end;

procedure header;
begin
      writeln(result,' frame ',version:4:2,'; reading frame evaluation of');
      readln(test);
      copyaline(test,result);
      writeln(result,' by the standard of');
      readln(norm);
      copyaline(norm,result);
      writeln(result);
end;

begin  (* frames *)
      writeln(output,' frame ',version:4:2);
      reset(norm);
      reset(test);
      rewrite(result);
      header;
      new(normparam);
      new(testparam);
         (* read the parameter records from the input files *)
      readencpar(norm,normparam,r,e);
      readencpar(test,testparam,r,e);
      if (vectorsize(normparam) = vectorsize(testparam)) then begin
         makevector(normvector,vectorsize(normparam));
         makevector(testvector,vectorsize(testparam))
      end
      else begin
         writeln(output,' input vectors of unequal size');
         halt;
      end;
      readvector(norm,normvector);
      replacevector(normvector,normparam^.wvlength);
         (* finish the header on the result file *)
      writeln(result,' ',normparam^.coding:1,
         ' is the oligo length used for comparison');
      writeln(result);
      writeln(result,' ':4,'frame 0',' ':4,'frame 1',' ':4,'frame 2');
      writeln(result);
      readln(test);
      readln(test);
      while not eof(test) do begin
         readvector(test,testvector);
         replacevector(testvector,testparam^.wvlength);
         evaluate;
      end;
1: end.