program derivative(xyin, derivativep, data, output);
(* derivative: take the derivative of x versus y in xyin

  Tom Schneider
  NCI/FCRDC Bldg 469. Room 144
  P.O. Box B
  Frederick, MD  21702-1201
  (301) 846-5581 (-5532 for messages)
  permanent email: toms@alum.mit.edu
  toms@ncifcrf.gov
  http://www.lecb.ncifcrf.gov/~toms/

  National Cancer Institute
  Laboratory of Experimental and Computational Biology

 *)

label 1; (* end of program *)

const
(* begin module version *)
version = 1.03; (* of derivative.p 1999 September 24
origin 1999 Sep 24 *)
updateversion = 1.00; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.derivative *)
(*
name
   derivative: take the derivative of x versus y in xyin

synopsis
   derivative(xyin: in, derivativep: in, data: out, output: out)

files
   xyin:  a data file consisting of two columns, the first is called x
      and the second is y.  Comment lines can begin with "*", and empty
      lines are allowed.

      output columns are:
       1: slope
       2: x
       3: y
       4: trigger value (see trigger parameter)

   data:  output of this program, a single data file.  The file may be
      passed to genhis for plotting.

   derivativep:  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.

      outwid, outdec: two integers that determine the width in characters
         of output numbers and their number of decimal places.

      scale: if the first character is 'l' then it is a log-log scale
         and log base 2 will be take of both x and y before computing the
         slope.

      trigger: when the slope is below the trigger, the trigger value
         in xyin will be 0, otherwise it will be 1.

   output: messages to the user

description

   Two data columns in the xyin file are read as x and y.
   The difference between each x value and its previous one is deltax.
   The difference between each y value and its previous one is deltay.
   The computed slope is deltay/deltax.
   The program will halt if deltax is zero.

examples

documentation

see also
   genhis.p, derivativep

author
   Thomas Dana Schneider

bugs

technical notes

*)
(* end module describe.derivative *)

var
   xyin, (* file used by this program *)
   derivativep, (* file used by this program *)
   data: text; (* file used by this program *)

(* 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.16 84 mar 12 tds/gds'; *)

(* begin module getxy *)
procedure getxy(var xyin: text; var x, y: real; var gotten: boolean);
(* get two numbers, x and y from the xyin file.  If end of file
was found, gotten is true, otherwise false. *)
var
  done: boolean; (* done searching for the data? *)
begin
{writeln(output,'in getxy');}
   gotten := false; (* so pesimistic! *)
   done := false;
   while not done do begin
      if eof(xyin)
      then done := true
      else if eoln(xyin)
           then readln(xyin)
           else if xyin^ = '*'
                then readln(xyin)
                else begin
                   readln(xyin, x, y);
                   gotten := true;
                   done := true;
{
write  (output,'x = ',x:1:1);
write  (output,' y = ',y:1:1);
writeln(output);
}
           end;
   end;
end;
(* end module getxy *)

(* begin module derivative.themain *)
procedure themain(var xyin, derivativep, data: text);
(* the main procedure of the program *)
var
   deltax: real; (* change in x *)
   deltay: real; (* change in y *)
   gotten: boolean; (* have we gotten new data? *)
   ln2: real; (* ln(2) *)
   parameterversion: real; (* parameter version number *)
   outwid: integer; (* output width in places *)
   outdec: integer; (* output decimal places *)
   previousx: real; (* previous value of x *)
   previousy: real; (* previous value of y *)
   scale: char; (* l means log-log scale *)
   slope: real; (* slope of curve *)
   trigger: real; (* trigger value *)
   x: real; (* current x *)
   y: real; (* current y *)
function log2(x:real): real;
(* compute log base 2 of x *)
begin
   log2 := ln(x)/ln2;
end; (* log2 *)
begin
   writeln(output,'derivative ',version:4:2);
   reset(derivativep);
   readln(derivativep, parameterversion);
   if parameterversion < updateversion then begin
      writeln(output, 'You have an old parameter file!');
      halt
   end;

   readln(derivativep, outwid, outdec); 
   readln(derivativep, scale); 
   readln(derivativep, trigger); 

   rewrite(data);
   writeln(data,'* derivative ',version:4:2);

   ln2 := ln(2);
   reset(xyin);

   getxy(xyin, previousx, previousy, gotten);

   if not gotten then begin
      writeln(output,'NO DATA');
      halt
   end;

   while not eof(xyin) do begin
      getxy(xyin, x, y, gotten);
      if gotten then begin
         if scale = 'l' then begin
            deltax := log2(x) - log2(previousx);
            deltay := log2(y) - log2(previousy);
         end
         else begin
            deltax := x - previousx;
            deltay := y - previousy;
         end;
         if deltax = 0 then begin
            writeln(output,'DIVISION BY ZERO');
            writeln(output,'DIVISION BY ZERO');
            write  (output,'x = ',x:1:1);
            write  (output,'previousx = ',previousx:1:1);
            writeln(output);
            halt
         end;
         slope := deltay/deltax;
         write(data, slope:outwid:outdec);
         write(data, ' ',x:outwid:outdec);
         write(data, ' ',y:outwid:outdec);
         if slope < trigger
         then write(data, ' 0')
         else write(data, ' 1');
         writeln(data);
         previousx := x;
         previousy := y;
      end;
   end;

end;
(* end module derivative.themain *)

begin
   themain(xyin, derivativep, data);
1: end.