program riden(color, ridenp, xyin, output);
(* riden: ring density graph

  Tom Schneider
  National Cancer Institute
  Laboratory of Mathematical Biology
  Frederick, Maryland  21701-1013
  toms@ncifcrf.gov

  1989
 *)

label 1; (* end of program *)

const
(* begin module version *)
version = 1.29; (* of riden.p 1996 November 18
origin 1989 November 19 *)
(* end module version *)

(* begin module describe.riden *)
(*
name
   riden: ring density graph

synopsis
   riden(color: in, ridenp: in, xyin: out, output: out)

files
   color: output of the ring program
   ridenp: parameter file for this program.  Two lines:
           First line: largest radial distance recorded
           Second line: number of bins to store the data in
   xyin:  histogram of the density
   output: messages to the user

description
   This program converts the graph generated by the ring program into
   a form that allows one to see if the results are as expected.

examples

documentation

see also
   ring.p

author
   Thomas Dana Schneider

bugs
   Program only works for D=2.  The curves don't match for D=4, but
   do for higher dimensions.  It is not obvious why.

technical notes

*)
(* end module describe.riden *)

(* more constants *)
   binsmax = 1000; (* maximum number of storage locations *)
   Rmaximum = 1.0; (* where the curves are to have their maximum *)

type
   (* an array for storing integers *)
   storage = array[0..binsmax] of integer;

var
   color: text; (* output of the ring program *)
   ridenp: text; (* output, input to xyplo *)
   xyin: text; (* parameter file of 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 = 2.46; (@ of ring.p 1989 Nov 19 *)

(* 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 = 2.46; (@ of ring.p 1989 Nov 19 *)

function pd(d: integer; r: real; twosigma2: real): real;
(* probability at a given radius in a given dimension *)
begin
   if r = 0.0 then pd := 0.0
              else pd := exp(-r*r/twosigma2) * exp((d-1)*ln(r))
end;

(* begin module riden.themain *)
procedure themain(var color, ridenp, xyin: text);
(* the main procedure of the program *)
var
   biggest: integer; (* the largest point in the histogram *)
   bins: integer; (* number of divisions of rmax to make *)
   D: integer; (* the dimensionality *)
   pdmax: real; (* the largest expected density value *)
   r: real;    (* radial distance from the origin coordinate of the point *)
   rmax: real; (* maximum radial distance to record *)
   s: integer; (* index for finding the dimensionality,
                   and indexing the store *)
   skipped: integer; (* number of points lost *)
   stored: integer; (* number of points captured *)
   store: storage; (* the place to build the histogram *)
   twosigma2: real; (* 2 times sigma squared, where
                       sigma is the standard deviation of the gaussian *)
   x: real; (* x coordinate of the point *)
   y: real; (* y coordinate of the point *)
begin
   writeln(output,'riden ',version:4:2);

   reset(ridenp);
   readln(ridenp,rmax);
   readln(ridenp,bins);
   if bins > binsmax then begin
      writeln(output,'too many bins, increase binsmax');
      halt
   end;

   (* determine the dimensonality *)
   reset(color);
   for s := 1 to 9 do readln(color);
   get(color);
   readln(color,D);
   if D < 2 then begin
      writeln(output,'D must be > 1');
      halt
   end;
   writeln(output, 'dimensionality: D = ',D:1);

   rewrite(xyin);
   writeln(xyin,'* riden ',version:4:2);
   reset(color);
   while color^ = '*' do copyaline(color,xyin);

   (* clear the store *)
   for s := 0 to bins do store[s] := 0;
   skipped := 0;

   (* build the histogram *)
   while not eof(color) do begin
      if color^ = 's' then begin (* pick up the simulated points only *)
         get(color); (* skip the s *)
         readln(color,x,y); (* grab the coordinates *)
         r := sqrt(x*x + y*y);
         s := round( (bins/rmax) * r);

         if s <= bins then store[s] := store[s] + 1
         else begin
            skipped := skipped + 1;
            writeln(output,'skipping r=',r:10:8)
         end
      end
      else readln(color)
   end;

   (* determine the largest stored value, and the number stored *)
   stored := 0;
   biggest := 0;
   for s := 0 to bins do begin
      if store[s] > biggest then biggest := store[s];
      stored := store[s] + stored;
   end;

{
         (* #1 *)
   twosigma2 := 2*sqr(Rmaximum)/sqrt(D-1); (* this worked with #1! *)
   (* ie, the fdr curve (4d) matched the data, but both were shifted high *)
}
        (* #2 *)
   twosigma2 := 2*sqr(Rmaximum)/(D-1);
   (* this forces the fdr curve to pass through 1 *)

   pdmax := pd(D,1.0, twosigma2);

   (* write it out *)
   writeln(xyin,'* symbol,',
                ' radius,',
                ' normalized histogram,',
                ' raw histogram,',
                ' expected fD(r)');

   (* identification tag *)
   writeln(xyin,'D=',D:1,
                   ' ', 2.0:10:8,
                   ' ', 0.9:10:8,
                   ' ', 0:5,
                   ' ', 0:10);

   for s := 0 to bins do begin
      r := rmax * s/bins;
      writeln(xyin,'s',
                   ' ', r:10:8,
                   ' ',(store[s]/biggest):10:8,
                   ' ', store[s]:5,
                   ' ', (pd(D, r, twosigma2)/pdmax):10:8)
   end;

   (* now give the expected data in the normalized column to allow plotting *)
   for s := 0 to bins do begin
      r := rmax * s/bins;
      writeln(xyin,'e',
                   ' ', r:10:8,
                   ' ', (pd(D, r, twosigma2)/pdmax):10:8,
                   ' ', store[s]:5,
                   ' ',(store[s]/biggest):10:8)
   end;

   writeln(output,'number of points       skipped: ',skipped:5);
   writeln(output,'number of points  histogrammed: ',stored:5);
   writeln(output,'number of points in input file: ',(skipped+stored):5);

end;
(* end module riden.themain *)

begin
   themain(color, ridenp, xyin);
1: end.