program smoothis(histog, smoothisp, smooth, output);
(* smoothis: smooth a histogram

  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.06; (* of smoothis.p 1999 December 4
origin 1999 Dec 3 *)
updateversion = 1.00; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.smoothis *)
(*
name
   smoothis: smooth a histogram

synopsis
   smoothis(histog: in, smoothisp: in, smooth: out, output: out)

files
   histog:  the output of the genhis program, consisting of comment lines
      that begin with "*" and data lines that have pairs of bin position and
      number of items at that bin position.

   smooth:  Same format as histog, but smoothed according to the parameters.

   smoothisp:  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 (integer): The range to smooth over.

   output: messages to the user

description

   The smoothis program smooths a histogram.  If the smoothing range is 1,
   then no smoothing is done.  If smooting is 2 then every bin is added to
   the previous bin and averaged.  The result is in the second bin position
   so that bin locations are on the high end.

examples

documentation

see also
   genhis.p,
   smoothisp

author
   Thomas Dana Schneider

bugs

technical notes

   memmax is a constant that defines the largest allowed smoothing range.

*)
(* end module describe.smoothis *)

(* begin module smoothis.const *)
   wid = 10; (* width of numbers to use for data *)
   dec = 5; (* number of decimal places to use for data *)
   memmax = 1000; (* largest allowed smoothing range *)
(* end module smoothis.const *)

var
   histog, (* file used by this program *)
   smoothisp, (* file used by this program *)
   smooth: 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 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 smoothis.themain *)
procedure themain(var histog, smoothisp, smooth: text);
(* the main procedure of the program *)
const
   debug = false; (* set to true to see debugging output *)
   debugwid = 5; (* width of numbers to use for debugging data *)
   debugdec = 2; (* number of decimal places to use for debugging data *)
   countmax = 10; (* how many data pieces to look at before halting while
                     debugging *)

var

   average: real; (* the average of values in the memory *)
   count: integer; (* counter for limiting debug output *)
   n: integer; (* index to memory for current number *)
   memory: array[1..memmax] of real; (* memory of past and current numbers *)
   memoryfull: boolean; (* have we filled enough of the memory yet? *)
   parameterversion: real; (* parameter version number *)
   position: real; (* the position in the histogram *)
   range: integer; (* range to average over *)
   total: real; (* current total in the memory *)

procedure cycle(var m: integer; range: integer);
(* step m through a cycle up to the range *)
begin
   m := succ(m);
   if m > range then m := 1;
end;

begin
   writeln(output,'smoothis ',version:4:2);
   reset(smoothisp);
   readln(smoothisp, parameterversion);
   if parameterversion < updateversion then begin
      writeln(output, 'You have an old parameter file!');
      halt
   end;

   readln(smoothisp, range);

   if range < 1 then begin
      writeln(output,'range must be positive');
      halt
   end;

   if range > memmax then begin
      writeln(output,'range must be less than memmax (',memmax:1,')');
      halt
   end;

   reset(histog);
   rewrite(smooth);

   writeln(output,'range is ',range:1);

   total := 0;
   memoryfull := false;
   for n := 1 to range do memory[n] := 0;;
   n := 1; (* start filling memory here *)

   if debug then count := 0;

   while not eof(histog) do begin
      if (eoln(histog)) or (histog^ = '*') then begin
         copyaline(histog, smooth);
      end
      else begin
         (* remove the most ancient value *) 
         total := total - memory[n];

         (* pull in the new value *) 
         readln(histog, position, memory[n]);

         (* add in the new value to the running total: *)
         total := total + memory[n];

if debug then write(output, position:(debugwid+5):debugdec,
         '   unsmoothed:  ',memory[n]:debugwid:debugdec);


	 (* when memory is filled, make a note of it *)
         if n = range then memoryfull := true;

	 if memoryfull then begin
            average := total / range;
            writeln(smooth, position:wid:dec, ' ',average:wid:dec);

if debug then write(output, '   average: ', average:debugwid:debugdec,
                            '   total: ',total:debugwid:debugdec);

{zzz}
         end;

         (* step cyclically through the memory *)
         cycle(n, range);

if debug then writeln(output);

{zzz}
if debug then count := succ(count);
if debug then if count >= countmax then halt;

      end
   end;

end;
(* end module smoothis.themain *)

begin
   themain(histog, smoothisp, smooth);
1: end.