program denri(data, denrip, den, xyplop, output);
(* denri: compute density for delta Ri

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

  National Cancer Institute
  Laboratory of Mathematical Biology

 *)

label 1; (* end of program *)

const
(* begin module version *)
version = 1.09; (* of denri.p 1999 June 28
1.09: 1999 June 28: documentation brush up.
1.08: 1997 April 1: major use
origin 1997 March 31 *)
updateversion = 1.00; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.denri *)
(*
name
   denri: compute density for delta Ri

synopsis
   denri(data: in, denrip: in, den: out, xyplop: out, output: out)

files
   data:  output of scan program

   den:  output of this program.  Use it as the xyin input to xyplo.

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

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

      1. numbertoprocess: number of data items to process (negative = all)

      2. wantfrom, wantto: position from-to range for computing

      3. wantminRi, wantmaxRi: bit from-to range for computing

   xyplop: control file for xyplo

   output: messages to the user

description

   Take the delta Ri values from the scan program and compute the density at
   each position and bit level.  Bit values from the scan data file are
   rounded and counts are placed in bins.  The output can then be plotted
   with xyplo.

examples

   denrip:

1.01      version of denri that this parameter file is designed for.
1000      number of data items to process (negative = all)
-20 20    wantfrom, wantto: position from-to range for computing
-70 20    wantminRi, wantmaxRi: bit from-to range for computing

documentation

see also
   scan.p, xyplo.p

author
   Thomas Dana Schneider

bugs

technical notes
   Four constants in module denri.const define the largest
   position range and bit range the program can handle.

*)
(* end module describe.denri *)

(* begin module denri.const *)
   rangefrom = -500; (* lowest position allowed *)
   rangeto   = +500; (* highest position allowed *)
   minRi = -70; (* lowest Ri allowed *)
   maxRi = +20; (* highest Ri allowed *)
(* end module denri.const *)

var
   data, (* file used by this program *)
   denrip, (* file used by this program *)
   den, (* file used by this program *)
   xyplop: 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 = 4.20; (@ of prgmod.p 1997 March 15 *)

(* begin module skipblanks *)
procedure skipblanks(var thefile: text);
(* skip over blanks until a non-blank, or end of line, is found *)
begin
      while (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 (thefile^ <> ' ') and not eoln(thefile) do get(thefile);
end;

procedure skipcolumn(var thefile: text);
(* skip over a data column *)
begin
   skipblanks(thefile); skipnonblanks(thefile)
end;
(* end module skipblanks version = 4.20; (@ of prgmod.p 1997 March 15 *)

(* begin module denri.themain *)
procedure themain(var data, denrip, den, xyplop: text);
(* the main procedure of the program *)
var
   bins: array[rangefrom..rangeto, minRi..maxRi] of integer;
      (* bins to store the data *)
   c: integer; (* count number of data items to fill in *)
   d: integer; (* index to delta Ri values *)
   deltaRi: real; (* delta Ri values *)
   numbertoprocess: integer; (* number of data items to process *)
   parameterversion: real; (* parameter version number *)
   position: integer; (* position in the sequences *)
   wantfrom, wantto: integer; (* range to process *)
   wantminRi, wantmaxRi: integer; (* bit range to process *)

procedure readparameters(var denrip: text);
(* read user defined parameters *)
begin
   reset(denrip);
   readln(denrip, parameterversion);
   if parameterversion < updateversion then begin
      writeln(output, 'You have an old parameter file!');
      halt
   end;

   readln(denrip, numbertoprocess);
   readln(denrip, wantfrom, wantto);
   readln(denrip, wantminRi, wantmaxRi);

   if wantfrom < rangefrom then begin
      writeln(output,'wantfrom cannot be less than ',rangefrom:1);
      halt;
   end;

   if wantto > rangeto then begin
      writeln(output,'wantto cannot be more than ',rangeto:1);
      halt;
   end;

   if wantfrom > wantto then begin
      writeln(output,'wantfrom cannot be more than wantto');
      halt;
   end;

   if wantminRi < minRi then begin
      writeln(output,'wantminRi cannot be less than ',minRi:1);
      halt;
   end;

   if wantmaxRi > maxRi then begin
      writeln(output,'maxRi cannot be more than ',maxRi:1);
      halt;
   end;

   if wantminRi > wantmaxRi then begin
      writeln(output,'wantminRi cannot be more than wantmaxRi');
      halt;
   end;

end;

procedure writeparameters(var out: text);
begin
   writeln(out,'* denri ',version:4:2);
   writeln(out,'* ', parameterversion:9:2,' parameterversion for denri');
   writeln(out,'* ', numbertoprocess:9,' numbertoprocess');
   writeln(out,'* ', wantfrom:4, ' ', wantto:4 ,' wantfrom, wantto');
   writeln(out,'* ', wantminRi:4, ' ',wantmaxRi:4,' wantminRi, wantmaxRi');
end;

procedure clearbins;
(* clear the bins out *)
begin
   writeln(output, 'clearing bins ...');
   for position := wantfrom to wantto
   do for d := wantminRi to wantmaxRi
   do bins[position,d] := 0;
   writeln(output, 'cleared');
end;

procedure fillonedatum;
(* read and store one data item *)
begin
   if data^ = '*'
   then readln(data)
   else begin
      skipcolumn(data);
      skipcolumn(data);
      skipcolumn(data);
      read(data,position);
      skipcolumn(data);
      skipcolumn(data);
      skipcolumn(data);
      skipcolumn(data);
      read(data,deltaRi);
      readln(data);
      d := round(deltaRi);

{
      if (d < minRi) or (d > maxRi)
         or
         (position < rangefrom) or (position > rangeto)
      then begin
         writeln(output,'position = ',position:1,' d=',d:1);
         writeln(output,'is out of the storage area');
halt
      end
      else begin
writeln(output);
writeln(output,'position = ',position:1,' d=',d:1);
writeln(output,' bins[position, d] = ', bins[position, d]:1);
writeln(output,'succ is: ',succ(bins[position, d]):1);
writeln(output,'filled ',position:1,' ',d:1);
}

      if (d >= wantminRi) and (d <= wantmaxRi)
         and
         (position >= wantfrom) and (position <= wantto)
      then begin
         bins[position, d] := succ(bins[position, d]);
      end;
   end;
end;

procedure fillbins;
(* fill the bins with data *)
begin
   writeln(output, 'filling bins ...');
   reset(data);
   if numbertoprocess < 0
   then while not eof(data) do fillonedatum
   else for c := 1 to numbertoprocess do fillonedatum;
   writeln(output, 'filled');
end;

procedure process;
(* process and output results *)
const
   wid = 7; (* character width of data output *)
   dec = 5; (* decimals for data output *)
var
   b: integer; (* bin value *)
   sum: integer; (* for determining the sum of data at a position *)
   r: real;  (* bin ratio = b/sum *)
   hue: real;  (* conversion of r to a color *)
   saturation: real;  (* conversion of r to a color *)
   brightness: real;  (* conversion of r to a color *)
{zzz}
begin
   writeln(output, 'processing...');
   writeln(den,'* columns:');
   writeln(den,'* 1: position');
   writeln(den,'* 2: bits');
   writeln(den,'* 3: total counts at this position');
   writeln(den,'* 4: (counts at this bit-position) /',
                   ' (total counts at this position)');
   writeln(den,'* 5: 1');
   writeln(den,'* 6: 1');
   writeln(den,'* 7: hue');
   writeln(den,'* 8: saturation');
   writeln(den,'* 9: brightness');
   for position := wantfrom to wantto do begin
      sum := 0;
      for d := wantminRi to wantmaxRi do sum := sum + bins[position,d];
      if sum > 0 then for d := wantminRi to wantmaxRi do begin
         b := bins[position,d];
{zzz}
         if b = 0 then begin
            brightness := 0.9;
            saturation := 0.0;
         end
         else begin
            r := b/sum;
            hue := 0.84*r + 0.16;
            saturation := 1.0;
            brightness := 1.0;
         end;
         writeln(den, position:4,              (* 1 *)
                      ' ',d:3,                 (* 2 *)
                      ' ',sum:wid,             (* 3 *)
                      ' ',b:wid,               (* 4 *)
                      ' ',1:wid,               (* 5 *)
                      ' ',1:wid,               (* 6 *)
                      ' ',hue:wid:dec,         (* 7 *)
                      ' ',saturation:wid:dec,  (* 8 *)
                      ' ',brightness:wid:dec); (* 9 *)
      end;
   end;
   writeln(output, 'processed');
end;

procedure mkxyplop(var xyplop: text);
(* make the xyplop file *)
var
  ix, iy: integer; (* intervals for x and y *)
begin
ix := wantto-wantfrom;
while ix > 20 do ix := trunc(ix/10);
iy := trunc((wantmaxRi-wantminRi)/10);
rewrite(xyplop);
writeln(xyplop,'3 6       zerox zeroy         ');
writeln(xyplop,'x ',wantfrom:1,' ',wantto:1,' zx min max          ');
writeln(xyplop,'y ',wantminRi:1,' ',wantmaxRi:1,'  zy min max          ');
writeln(xyplop,ix:1, ' ',iy:1,' 2 2 xinterval yinterval ');
writeln(xyplop,'4 4       xwidth    ywidth    ');
writeln(xyplop,'0 0       xdecimal  ydecimal  ');
writeln(xyplop,'15 15     xsize     ysize     ');

writeln(xyplop,'position, bases');
write(xyplop,'deltaRi, bits.');
writeln(xyplop,'  number of sequences at position 0: ',bins[0,0]:1);

writeln(xyplop,'a         zc                  ');
writeln(xyplop,'n 2       zxl base            ');
writeln(xyplop,'n 2       zyl base            ');
writeln(xyplop,'          ********************');
writeln(xyplop,'1 2       xcolumn   ycolumn   ');
writeln(xyplop,'0         symbol column       ');
writeln(xyplop,'5  6      xscolumn  yscolumn  ');
writeln(xyplop,'7 8 9     hue saturation brightnes');
writeln(xyplop,'          ********************');
writeln(xyplop,'R         symbol-to-plot      ');
writeln(xyplop,'i         symbol-flag         ');
writeln(xyplop,'-1.0      symbol sizex        ');
writeln(xyplop,'-1.0      symbol sizey        ');
writeln(xyplop,'n         connection size     ');
writeln(xyplop,'n  0.125  linetype  size      ');
writeln(xyplop,'          ********************');
writeln(xyplop,'.');
writeln(xyplop,'          **** version 8.50 of');
writeln(xyplop,'l 0 0 0.125  User defined line');
{
writeln(xyplop,'l -1 0 0.125  User defined line');
}
end;

begin
   writeln(output,'denri ',version:4:2);
   readparameters(denrip);
   rewrite(den);
   writeparameters(den);
   writeparameters(output);
   clearbins;
   fillbins;
   process;
   mkxyplop(xyplop);
   writeln(output, 'done');
end;
(* end module denri.themain *)

begin
   themain(data, denrip, den, xyplop);
1: end.