program l2cvwarf;
{(snp9: STB-27)
First written: September 25, 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. Read data from _data2 and
input values from _inpval (ASCII). 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,intsquare,estexp : real;
  bin,frequ,index,mv                                     : vectorint;
  x,kwe,fm,cv,hv                                         : vector;
  inpval,data,resul                                      : text;


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];
   if kerneltype=6 then
      delta:=delta*4;
   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:= 1 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 CreateWeight;
var
l: integer;
part1,part2 : real;

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

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

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

begin
  for kw:=0 to nl-1 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
kl,z :integer;

begin
   for M:=mstart to mend do
      begin
         h:=M*delta;
         intsquare:=0;
         estexp:=0;
         for z:=1 to numbin do
            fm[z]:=0;
         CreateWeight;
         WeightBins;
         for kl:=1 to nl do
            estexp:=estexp+frequ[kl]*fm[index[kl]];
         estexp:=estexp-kwe[0]*n;
         estexp:=estexp/(n-1);
         for z:=1 to numbin do
            intsquare:=intsquare+sqr(fm[z]);
         intsquare:=intsquare*delta;
         cv[M-mstart+1]:=intsquare-2*estexp;
         mv[M-mstart+1]:=M;
         if kerneltype=6 then
           hv[M-mstart+1]:=mv[M-mstart+1]*delta/4
         else
           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,cv[rcount]:12:8,hv[rcount]:12:4);
  close(resul);
end;

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