program bcvwarpf;
{(snp9: STB-27) 
First written: September 27, 1994. Last revised: March 7, 1995.
Authors: Isaias Hazarmabeth Salgado Ugarte, Makoto Shimizu and 
Toru Taniuchi. Univ. of Tokyo, Fac. of Agriculture, Dept. of 
Fisheries. Yayoi 1-1-1, Bunkyo-ku, Tokyo 113, Japan. 
No results to screen, only to text file. This Pascal program performs
least squares cross validation for WARPing density estimation, using
modified algorithms from Haerdle (1991) "Smoothing Techniques with
Implementations in S. Springer-Verlag Series in Statistics, New York }

const
  maxsize = 1000;
  maxsize2 = 1000;

type
  vector = array[0..maxsize] of real;
  striname = string[30];
  vectorint = array[0..maxsize2] of integer;

var
  n,i,M,mstart,mend,mnumber,numbin,kerneltype,indexi,nl        :integer;
  start,delta,h,origin,maxval,minval,cm,score,tau,kappa,zeta   : real;
  bin,frequ,index,mv                                           : vectorint;
  x,kwe,fm,bcv,hv,ma                                           : vector;
  inpval,data,resul                                            : text;

function Power(x:real; n:integer): real;
  begin
  if n=1 then
    Power:=x
  else if odd(n) then
    Power:=x*sqr(Power(x, n div 2))
  else
    Power:=sqr(Power(x, n div 2))
  end;

procedure InputValues;
var i: integer;
begin
   assign(data,'_data2.raw');
   reset(data);
   i:=1;
   while (not eof(data)) do
    begin
       read(data,x[i]);
       i:=i+1;
    end;
   close(data);
   n:=i-2;
   assign(inpval,'_inpval.raw');
   reset(inpval);
   read(inpval,delta,kerneltype,mstart,mend);
   close(inpval);
end;

procedure InitialCalc;
var maxnoemp,ib: integer;

begin
   mnumber:=mend-mstart+1;
   maxval:=x[n];
   minval:=x[1];
   start:=minval-delta*(Mend+0.1);
    if (start<0) then
      origin:=(round((start/delta)-0.5)-0.5)*delta
   else
      origin:=(trunc(start/delta)-0.5)*delta;
   numbin:=round(((maxval+delta*(Mend+0.1)-origin)/delta)+0.5)+1;
   if (n+1<numbin) then
      maxnoemp:=n+1
   else
      maxnoemp:=numbin;
   nl:=0;
   indexi:=0;
   for ib:= 0 to numbin+1 do
      begin
        bin[ib]:=0;
        frequ[ib]:=0;
        index[ib]:=0;
      end;
end;

procedure BinningData;
var
   ia,bindx,bi : integer;
   sum1, sum2,interval : real;
begin
   for i:=1 to n do
    begin
      indexi:= trunc((x[i]-origin)/delta);
      if (bin[indexi]=0) then
         begin
            nl:= nl+1;
            frequ[nl]:=1;
            index[nl]:=indexi;
            bin[indexi]:=nl;
          end
      else
        frequ[bin[indexi]]:=frequ[bin[indexi]]+1;
    end;
end;

procedure CalcFunctional;
var
   divisor : real;

begin
   case kerneltype of
      1: begin {Quartic}
         tau:=(16*ma[4]+32*ma[2]-13)/(112*ma[2]+28);
         divisor:=(16*ma[4]-1)*(4*ma[2]+1);
         kappa:=15*ma[1]/7*(1-(32*ma[4]-28*ma[2]+1)/divisor);
         zeta:=30/7*(336*ma[3]-420*ma[2]+85*ma[1]+104)/divisor;
         end;
      2: begin {Triweight}
         tau:=(32*ma[6]+80*ma[4]+62*ma[2]-69)/(288*ma[4]+72*ma[2]+60);
         divisor:=(4*ma[2]-1)*Power((24*ma[4]+6*ma[2]+5),2);
         kappa:=350*ma[1]/143*(1-(4492.8*ma[8]-979.2*ma[6]+2124*ma[4]-866.4*ma[2]-2755.2)/(3*divisor));
         zeta:=630*(128*ma[7]-448*ma[5]+560*ma[4]-280*ma[3]-280*ma[2]+320*ma[1]+35)/divisor;
         end;
      end; {case}
end;

procedure CreateWeight2;
var
l: integer;
part1,part2 : real;

begin
  case kerneltype of
    1: begin
         cm:=0.9375/((1-0.0625/(Power(M,4)))*(n*h));
         for l:=0 to M-1 do
           kwe[l]:=cm*Power((1-Power((l/M),2)),2);
       end;
    2: begin
         part1:=1+0.14583333/Power(M,4);
         part2:=0.052083333/Power(M,6);
         cm:= 1.09375/((part1-part2)*(n*h));
         for l:=0 to M-1 do
           kwe[l]:=cm*Power(1-Power((l/M),2),3);
       end;
    end;{case}
end;

procedure WeightBins;
var kw,lw,iw: integer;

begin
  for kw:=1 to nl do
      for lw:=1-M to M-1 do
        fm[index[kw]+lw]:= fm[index[kw]+lw]+kwe[abs(lw)]*frequ[kw];
end;

procedure MainLoop;
var
im,kl,z :integer;

begin
   for M:= mstart to mend do
      begin
         ma[1]:=M;
         for im:=2 to 8 do
              ma[im]:=Power(ma[1],im);
         h:=M*delta;
         score:=0;
         for z:=0 to numbin do
            fm[z]:=0;
         CalcFunctional;
         CreateWeight2;
         WeightBins;
         for z:=1 to numbin-1 do
            score:=score+Power((fm[z+1]-2*fm[z]+fm[z-1]),2);
         score:=score+Power((fm[numbin-1]-fm[numbin-2]),2);
         score:=score+Power((2*fm[0]-fm[1]),2);
         score:=score+Power((fm[0]),2)+Power(fm[numbin-1],2);
         score:=score*delta*0.25;
         score:=score-zeta/(4*n*ma[2]*delta);
         score:=score*Power(tau,2)+1/180;
         bcv[M-mstart+1]:=score+kappa/(3*n*ma[2]*delta);
         mv[M-mstart+1]:=M;
         hv[M-mstart+1]:=mv[M-mstart+1]*delta;
      end;
end;

procedure ResulFile;
var rcount: integer;
begin
  assign(resul,'resfile');
  rewrite(resul);
  for rcount:=1 to mnumber do
    writeln(resul, mv[rcount]:12,bcv[rcount]:12:8,hv[rcount]:12:4);
  close(resul);
end;

begin {main}
  InputValues;
  InitialCalc;
  BinningData;
  MainLoop;
  ResulFile;
end.