program sphere(spherep, sigma, xyin, output);
(* sphere: plot density of shannon spheres

  Dr. Thomas D. Schneider
  National Cancer Institute
  Laboratory of Experimental and Computational Biology
  Frederick, Maryland  21702-1201
  toms@ncifcrf.gov
  permanent email: toms@alum.mit.edu
  http://www.lecb.ncifcrf.gov/~toms/

   *)

label 1; (* end of program *)

const
(* begin module version *)
version = 1.47; (* of sphere.p 2001 July 27
2001 Jul 25, 1.43: document spherep.fractal
2001 Jul 25, 1.43: noninteger D allowed
2001 Jul 25, 1.42: upgrade documentation.
2001 Jul 25, 1.41: merge documentation with published ccmm version
1998 Jun  1, 1.40: add sphere.xyplop pointer
1994 Sep  6, 1.39: last previous change (used for ccmm paper)
1988 Jan 24, 1.00: origin
*)
(* end module version *)

(* begin module describe.sphere *)
(*
name
   sphere: plot density of shannon spheres

synopsis
   sphere(spherep: in, sigma: out, xyin: out, output:out)

files
   spherep: parameters.
     The first line is the step size interval (0.01 works well).

     the second line is the maximum radius to calculate out to (= maxr,
         3.1 works well).

     Each following line is a dimension to plot.

     If the dimension number is negative, it must be followed on the same line
     by the coordinates of the position to place the dimension numeral.
     The absolute value of the dimension is used in the calculation.
     
     If the dimension is negative AND not an integer, the coordinates of the
     position must be followed by the number of decimal places to display
     the dimension.

   sigma: lists the estimates for Rmaximum +/- sigma,
      taken as the radius when the curve passes through exp(-1/2).

   xyin: input to xylop, the plot

   output: messages to the user

description
   Create a graph of radius versus density of Shannon spheres
   at various given dimensions.  The output is run through xyplo.

   The function is:

   pd(R) = R^(D-1) * exp( sqr(R)/ (2* sqr(sigma)))

   where '^' means to exponentiate and
   where sqr(sigma) * (D-1) - sqr(Rmaximum)
   so setting Rmaximum = 1 relates sigma and D.

   The graph is in the range (0,0) to (r=maxr,1)).
   The curve is normalized so that its maximum is at (1,1).
   (except when dimension = 1, where it is at (1,0).

   Since xyplo can't plot several separate curves, without being
   told each symbol, this program simply starts at (0,pd(r)), draws
   the curve to (maxr,pd(maxr)), then circles back by drawing lines
   to the x axis (2*maxr,0) and then the origin (0,0).  By setting
   the region that xyplo plots below maxr, one gets nice, fully
   correct curves that do not appear to be connected.

documentation
   [1988 jan 23,5]

@article{Schneider.ccmm,
author = "T. D. Schneider",
title = "Theory of Molecular Machines.
{I. Channel} Capacity of Molecular Machines",
journal = "J. Theor. Biol.",
volume = "148",
thenumber = "1",
pages = "83-123",
comment = "{(Note: The figures were printed out of order!
Fig. 1 is on p. 97)}",
note = "http://www-lecb.ncifcrf.gov/$\sim$toms/paper/ccmm/",
year = 1991}

see also

   {Schneider.ccmm paper:}
   http://www.lecb.ncifcrf.gov/~toms/paper/ccmm

   {example parameters:} spherep
   {example plotting file for xyplo:} sphere.xyplop
   {plotting program:} xyplo.p

   {resulting graph, postscript ... :} sphereinteger.ps
   {resulting graph, jpg .......... :}
   sphereinteger.jpg

   {example of fractal graph parameters:} spherep.fractal
   {resulting graph, postscript ... :} spherefractal.ps
   {resulting graph, jpg .......... :}
   spherefractal.jpg

   {Related programs:}
   {compress D dimensional sphere to 2D: ...} ring.p
   {plot output of ring: ...................} riden.p
   {Match fdr curve: .......................} fdr.p

author
   Thomas Dana Schneider

bugs
   none known

*)
(* end module describe.sphere *)

(* begin module sphere.const *)
   wid = 10; (* width of numbers *)
   dec = 8; (* number of decimal places of numbers *)
   Rmaximum = 1.0; (* where the curves are to have their maximum *)
(* end module sphere.const *)

var
   spherep: text; (* parameters *)
   sigma: text; (* list of sigma values *)
   xyin: text; (* input to xyplo *)

function near(n: real; spot: real): boolean;
(* return true if n is near spot *)
const
   tolerance = 5e-2; (* how close we have to be to be "near" *)
begin
   near:=(n > spot - tolerance) and (n < spot + tolerance)
end;

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

function lnpd(d: real; r: real; twosigma2: real): real;
(* natural log f the pd function, to make calculation of high
dimensions possible *)
(* probability at a given radius in a given dimension *)
begin
   lnpd := (-sqr(r)/twosigma2) + ((d-1)*ln(r))
end;

(* begin module sphere.themain *)
procedure themain(var spherep, sigma, xyin: text);
(* the main procedure of the program *)
var
   D: real; (* dimension, fractal dimension allowed! *)
   Dx: real; (* x coordinate to plot the value of D *)
   Dy: real; (* y coordinate to plot the value of D *)
   decimals: integer; (* number of decimal places to display the dimension *)
   exphalf: real; (* exponent of -1/2 *)
   interval: integer; (* what segment of the curve we are currently plotting.
      1: before first intercept with exphalf,
      2: before rmax has been reached,
      3: before the second intercept with exphalf has been reached,
      4: tail of the curve *)
   lastvalue: real; (* the previous value of the curve *)
   lnpdmax: real; (* ln of maximum value of pd curve, for normalizing *)
   lnvalue: real; (* intermediate value: ln of the final result *)
   maxr: real; (* maximum radius to calculate to *)
   pdmax: real; (* maximum value of pd curve, for normalizing *)
   r: real; (* current radius *)
   sigmaplus, sigmaminus: real; (* estimates of sigma above and below rmax *)
   step: real; (* the step between calculations *)
   twosigma2: real; (* 2 times sigma squared, where
                       sigma is the standard deviation of the gaussian *)
   value: real; (* the final result *)
begin
   writeln(output,'sphere ',version:4:2);
   rewrite(sigma);
   writeln(sigma,'* sphere ',version:4:2);
   writeln(sigma,'* D: dimension');
   writeln(sigma,'* sigma-: deviation to inside of sphere');
   writeln(sigma,'* sigma+: deviation to outside of sphere');
   writeln(sigma,'* average: average of absolute values of the sigmas');
   writeln(sigma,'* estimate: of sigma (1/sqrt(2*(D-1))');
   writeln(sigma,'*');

   reset(spherep);
   readln(spherep, step);
   if step <= 0 then begin
      step := 0.2;
      writeln(output, 'step must be positive, using ',step:wid:dec);
   end;
   readln(spherep,maxr);

   rewrite(xyin);
   writeln(xyin,'* sphere ',version:4:2);
   writeln(xyin,'* step: ',step:wid:dec);
   writeln(xyin,'* maximum r: ',maxr:wid:dec);
   writeln(xyin,'* kind (d=dimension numeral, c = curve) r, pd(r)');
   exphalf := exp(-1/2);

   while not eof(spherep) do begin
      read(spherep,D);
      if D < 0 
      then begin (* plot the dimension number at the coordinates given *)
         read(spherep, Dx, Dy);
         D := abs(D);
         if D = round(D)
         then decimals := 0
         else read(spherep, decimals);
         readln(spherep);

         if D = round(D)
         then begin
            writeln(xyin,'D=',round(D):1, ' ', Dx:wid:dec, ' ', Dy:wid:dec);
            writeln(output,'D = ',round(D):1,
                        ' plotted at ', Dx:wid:dec, ' ', Dy:wid:dec);
         end
         else begin
            writeln(xyin,'D=',D:1:decimals, ' ', Dx:wid:dec, ' ', Dy:wid:dec);
            writeln(output,'D = ',D:1:decimals,
                        ' plotted at ', Dx:wid:dec, ' ', Dy:wid:dec);
         end;

{zzz}
{
writeln(output,'decimals =',decimals:1);
writeln(output,'round(D) =',round(D):1);
writeln(output,'      D  =',D:1:dec);
if D = round(D)
then writeln(output,'D  = round(D) SAME')
else writeln(output,'D <> round(D) DIFFERENT: fractal!');
writeln(output);
}

      end
      else begin
         readln(spherep);
         if D = round(D)
         then writeln(output,'D = ',D:wid)
         else writeln(output,'D = ',D:wid:dec);
      end;

      write(sigma,'D:',D:5,' ');
      if D > 1
      then begin
         twosigma2 := 2*sqr(Rmaximum)/(D-1);
         (* pdmax := pd(D,1.0, twosigma2); *)
         lnpdmax := lnpd(D,1.0, twosigma2);
      end
      else begin (* make gaussian *)
         twosigma2 := 2*sqr(Rmaximum);
         pdmax := 1.0;
         lnpdmax := ln(pdmax);
      end;

      sigmaplus := 0;
      sigmaminus := 0;
      interval := 1;
      lastvalue := -maxint; (* preset for the way up *)
      r := step;
      while r <= maxr do begin
         (* writeln(xyin,r:wid:dec,' ',
                        (pd(D,r, twosigma2)/pdmax):wid:dec) *);
         lnvalue := (lnpd(D,r, twosigma2)-lnpdmax);

         (* by testing for the log of this value, we can avoid
            exponentiation of very small values.  the result is not
            be affected, for these values come up only when the final
            result is nearly zero anyway.
            Just don't even plot them! *)
         if lnvalue >= -25
         then begin
            value := exp(lnvalue);
            writeln(xyin,'c ',r:wid:dec,' ',value:wid:dec);

            case interval of
               1: begin (* before first intercept *)
                     if value > exphalf then begin
                        sigmaminus :=
                             ((r-step
                               +step*(exphalf-lastvalue)
                                    /(value-lastvalue))
                                    -Rmaximum);
                        lastvalue := maxint; (* reset for the way down *)
                        interval := interval + 1;
                     end
                     else lastvalue := value
                  end;
               2: begin (* before peak *)
                     (* wait for the maximum of the curve *)
                     if r >= Rmaximum then interval := interval + 1;
                  end;
               3: begin (* before second intercept *)
                     if value < exphalf then begin
                        sigmaplus :=
                             ((r+step
                               -step*(exphalf-value)
                                    /(lastvalue-value))
                                    -Rmaximum);
                        write(sigma,
                           '  sigma-: ',sigmaminus:5:3,
                           '  sigma+: ',sigmaplus:5:3,
                           '  average: ',((sigmaplus-sigmaminus)/2.0):5:3);
                        if D > 1
                        then writeln(sigma,
                           '  estimate: ',(1/sqrt(2*(D-1))):5:3)
                        else writeln(sigma,
                           '  exactly!: ',exphalf:5:3);
                        interval := interval + 1;
                     end
                     else lastvalue := value
                  end;
               4: ; (* after second intercept *)
            end;
         end;
         r := r + step
      end;

      (* make a line to the axis, off the edge of the plot: *)
      write(xyin,'c ',(2*maxr):wid:dec,' ',0.0:wid:dec);
      writeln(xyin, ' move to axis');

      (* move back to the origin along the axis: *)
      write(xyin,'c ',0.0:wid:dec,' ',0.0:wid:dec);
      writeln(xyin, ' move to origin');
   end;
end;
(* end module sphere.themain *)

begin
   themain(spherep, sigma, xyin);
1: end.