program fdr(histog,fdrp,outfile,output);
(* computes f_D(r) *)

const
(* begin module version *)
  version = 1.04; (* of fdr.p 1995 Jan 14 *)
(*
  Revision history: v. 1.00 started   = 10 Jan 1995
  Written at National Cancer Institute, Frederick, Maryland
*)
(* end module version *)
 
(* begin module describe.fdr *)
(* 
name
  fdr:  computes best values of sigma/D for f_D(r) fit of real data.

synopsis
  fdr(histog: in, fdrp:in, outfile: out, output: out)

files
  histog:  output file from genhis.
  fdrp:  parameter file for fdr.
    sigma,sigmastep (reals): initial value and incremental increase in sigma.
    dim,dimstep (reals): initial value and incremental increase in D.
    temperature (real): Kelvin temperature of system.
  outfile:  writes best value of sigma/D and error in f_D(r) fit.
  output:  sends messages to the user.

description
              D-1      2       2
             r   exp(-r /2sigma ) 
  f (r) =  -----------------------
   D                       D D/2-1
           Gamma(D/2) sigma 2

  We compute the values of D, sigma that minimze the least-squares error
  between f_D(r) and the actual data read in from histog.

examples 
  See fdr.eg (not yet written ) 

documentation
  See Appendix 3 to TDS's CCMM paper.

see also
  genhis.p

author
  V. Jejjala (vishnu@ncifcrf.gov) 

bugs
  Unknown but existent (existential bugs, as it were... :-) ).  If you find a 
  bug, please forward your discovery to the author.

disclaimer
  The author of this code is in no way responsible for any havoc that this
  program may cause to your computer system, your files, or your reputation.  
  Use this program and the data it generates at your own risk.  

technical notes
  Use large data sets in histog when possible.  

*) 
(* end module describe.fdr *)

(* begin module globals *)

(* more constants *) 
  k=1.380662e-23; (* Boltzmann's constant *)
  maxval=50; (* maximum number of points to consider on f_D(r) *)
  maxpass=100; (* maximum number of passes of program *)
  maxpassp=101; (* maxpass + 1 *)
  big=1000000.0; (* very big real number *)
  machep=-500.0; (* exp(machep) ~ machine epsilon *)

type
  singlept=record
    start,histval:real
  end;
  recdata=record
    allpts:array[1..maxval] of singlept; 
    sig,D,err:real
  end;
  rdata=array[0..maxpassp] of real;

var
  histog, (* output of genhis from which data are read *)
  fdrp, (* parameter file *)
  outfile:text; (* file to which data are written *)
(* end module globals *)

(* begin module expo *)
function expo(a,b:real):real;
(* yields a^b *)
begin
  if a=0.0 then expo:=0
  else if b*ln(a)<machep then expo:=0
  else expo:=exp(b*ln(a))
end;
(* end module expo *)

(* begin module gamma *)
function gamma(xx:real):real;
(* yields gamma(xx) using expansion formula/recurrance relation;
   code modified from _Numerical Recipes_ *)
const
  stp=2.50662827465; (* constant to compute gamma *)
var
  x,tmp,ser, (* temporary variables *)
  gammaln:real; (* ln(gamma(x)) *)
begin
  if xx<1 then x:=xx
  else x:=xx-1.0;
  tmp:=x+5.5;
  tmp:=(x+0.5)*ln(tmp)-tmp;
  ser:=1.0+76.18009173/(x+1.0)-86.50532033/(x+2.0)+24.01409822/(x+3.0)-
       1.231739516/(x+4.0)+0.120858003e-2/(x+5.0)-0.536382e-5/(x+6.0);
  gammaln:=tmp+ln(stp*ser);
  if xx<1 then gamma:=exp(gammaln)/xx
  else gamma:=exp(gammaln)
end;
(* end module gamma *)

(* begin module f *)
function f(r,sigma,dim:real):real;
(* computes f_D(r) *)
var
  g,h:real; (* temporary variables *)
begin 
  if -1*expo(r,2)/(2*expo(sigma,2))<machep then g:=0.0
  else g:=exp(-1*expo(r,2)/(2*expo(sigma,2)));
  h:=gamma(dim/2)*expo(sigma,dim)*expo(2,(dim/2-1));
  f:=expo(r,(dim-1))*g/h
end;
(* end module f *)

(* begin module readdata *)
procedure readdata(var histog:text;var datareal:recdata;var Rseq,Rstart,
                   Rstep:real);
(* reads in data from histog *)
var
  total:real; (* total number of sequences *)
  dummy, (* dummy variable *)
  count1,count2:integer; (* for loops *)
begin
  reset(histog);
  count1:=0;
  total:=0.0;
  while not eof(histog) do begin
    if histog^='*' then readln(histog)
    else begin
      count1:=count1+1;
      readln(histog,datareal.allpts[count1].start,dummy);
      datareal.allpts[count1].histval:=dummy;
      total:=dummy+total
    end;
  end;
  for count2:=1 to count1 do begin 
    Rseq:=datareal.allpts[count2].histval*datareal.allpts[count2].start+Rseq;
    datareal.allpts[count2].histval:=datareal.allpts[count2].histval/total
  end;
  Rseq:=Rseq/total; 
  Rstart:=datareal.allpts[1].start;
  Rstep:=datareal.allpts[2].start-datareal.allpts[1].start;
  for count2:=count1+1 to maxval do begin
    datareal.allpts[count2].start:=Rstart+(count2-1)*Rstep;
    datareal.allpts[count2].histval:=0
  end
end;
(* end module readdata *)

(* begin module readpar *)
procedure readpar(var fdrp:text;var sigma,sigmastep,dim,dimstep,temp:real);
(* reads in parameters from fdrp *)
begin
  reset(fdrp);
  readln(fdrp,sigma,sigmastep);
  readln(fdrp,dim,dimstep);
  readln(fdrp,temp)  
end;
(* end module readpar *)

(* begin module computeerror *)
procedure computeerror(var databest:recdata;var error:rdata;datareal:recdata;
                       sigma,dim,kT,Rstart,Rstep,sigmastep,dimstep:real;pass:
                       integer);
(* computes least-squares error--(sigma_(f_D(r)))^2 *)
var
  datameas:recdata; (* f_D(r) data for this particular value of D,sig,r *)
  r, (* current value of r *)
  diff, (* difference between f_D(r) and actual curve *) 
  err:real; (* least-squares error in f_D(r)--sum of diff over all r *)
  count:integer; (* for loops *)
begin
  datameas.sig:=sigma;
  datameas.D:=dim;
  for count:=1 to maxval do begin
    r:=Rstart+(count-1)*Rstep;
    datameas.allpts[count].start:=r;
    datameas.allpts[count].histval:=f(r,sigma,dim)/kT;
    diff:=datareal.allpts[count].histval-datameas.allpts[count].histval;
    if diff<0 then diff:=-1*diff;
    err:=expo(diff,2)+err
  end;
  datameas.err:=err/maxval;
{temp below}
  writeln(output,sigma,dim,datameas.err/225*0.84+0.16,1,1,sigmastep,dimstep);
{temp above}
  if datameas.err<error[pass] then error[pass]:=datameas.err;
  if datameas.err<databest.err then begin
    for count:=1 to maxval do begin 
      databest.allpts[count].start:=datameas.allpts[count].start;
      databest.allpts[count].histval:=datameas.allpts[count].histval
    end;
    databest.sig:=datameas.sig;
    databest.D:=datameas.D;
    databest.err:=datameas.err
  end
end;
(* end module computeerror *)

(* begin module writefile *)
procedure writefile(var outfile:text;databest:recdata;sigmastep,dimstep:real;
                    pass:integer);
(* writes out best sigma/D to outfile *)
var
  count:integer; (* for loops *)
begin
  writeln(outfile,'* Sigma:  ',databest.sig,' +/-',sigmastep);
  writeln(outfile,'* Dimension:  ',databest.D,' +/-',dimstep);
  writeln(outfile,'* Std. Deviation:  ',expo(databest.err,0.5));
  writeln(outfile,'* in ',pass,' passes');
  for count:=1 to maxval do begin
    write(outfile,databest.allpts[count].start);
    writeln(outfile,databest.allpts[count].histval)
  end
end;
(* end module writefile *)

(* begin module themain *)
procedure themain(var histog,fdrp,outfile:text);
(* main procedure of fdr.p *)
var
  datareal,databest:recdata; (* actual and best-fit data *) 
  error:rdata; (* array of data at various values o dim/sigma *)
  Rseq, (* rough estimate of R_seq *)
  Rstart,Rstep, (* initial value and steps in gradation of r *)
  sigma,sigmastep, (* initial value and steps in gradation of sigma *)
  dim,dimstep, (* initial value and steps in gradation of D *)
  dime, (* actual value of D used at any given time *)
  temp,kT:real; (* temperature and k*temperature *)
  pass, (* number of tries done *)
  count:integer; (* for loops *)
  result:boolean; (* decides when we're done *)
begin
  rewrite(outfile);
  readdata(histog,datareal,Rseq,Rstart,Rstep);
  pass:=1;
  result:=false;
  databest.err:=big;
  error[0]:=big;
  repeat
    if pass=1 then begin
      readpar(fdrp,sigma,sigmastep,dim,dimstep,temp);
      sigma:=sigma;
      kT:=1
    end;
    dime:=dim;
    sigma:=sigma+sigmastep;
    for count:=0 to 3*round(rseq) do begin
      dime:=dime+dimstep;
      error[pass]:=big;
      computeerror(databest,error,datareal,sigma,dime,kT,Rstart,Rstep,
                   sigmastep,dimstep,pass)
    end;
    if ((error[pass]<error[pass-1]) and (pass>1)) or (pass>maxpass) then 
    begin
      writefile(outfile,databest,sigmastep,dimstep,pass);
      result:=true
    end;
    pass:=pass+1
  until result=true
end;
(* end module themain *)

(* begin module mainroutine *)
begin
  themain(histog,fdrp,outfile)
end.
(* end module mainroutine *)