program gwarpren;
{First written: September 30, 1994. Last revised, November 11,1994.
Authors: Isaias Hazarmabeth Salgado Ugarte, Makoto Shimizu and
Toru Taniuchi.
No results to screen, only to text file. This Pascal program calculates
the adjusted prediction error G(M) for warping regression Nadaraya-Watson
estimate, using modified algorithms and C programs modified 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;
  vectorint = array[0..maxsize2] of integer;

var
  n,selector,M,mstart,mend,mnumber,numbin,kerneltype,j,nl,binboundh,binboundl,nlweight,index2 :integer;
  start,delta,h,origin,maxval,minval,boundary,rm,fm,mm      : real;
  bin,counts,index,mv,indexweight                           : vectorint;
  x,y,ysum,ysquaresum,kwe,score,hv                          : 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;

function SelecFun (u,v :real): real;
begin
   case selector of
      1:   SelecFun:=1+2*u/v;
      2:   SelecFun:=1/Power((1-u/v),2);
      3:   SelecFun:=exp(2*u/v);
      4:   Selecfun:=(1.0+u/v)/(1.0-u/v);
      5:   Selecfun:=1.0/(1.0-2.0*u/v);
   end; {case}
end;

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

procedure InitialCalc;
var maxnoemp,ii: 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;
   binboundl:=trunc(((x[trunc(n*boundary/2)])-origin)/delta);
   binboundh:=trunc(((x[trunc(n-n*boundary/2)])-origin)/delta);
   if (n+1<numbin) then
      maxnoemp:=n+1
   else
      maxnoemp:=numbin;
   nl:=0;
   nlweight:=0;
   for ii:= 0 to numbin+1 do
      begin
        bin[ii]:=0;
        counts[ii]:=0;
        ysum[ii]:=0;
        ysquaresum[ii]:=0;
        kwe[ii]:=0;
        index[ii]:=0;
        score[ii]:=0;
      end;
end;

procedure BinningData;
var
ib : integer;
begin
   for ib:=1 to n do
    begin
      j:= trunc((x[ib]-origin)/delta);
      if (bin[j]=0) then
         begin
            nl:= nl+1;
            counts[nl]:=1;
            ysum[nl]:=y[ib];
            ysquaresum[nl]:=Power(y[ib],2);
            index[nl]:=j;
            bin[j]:=nl;
            if (j>=binboundl) and (j<=binboundh) then
               begin
                  nlweight:=nlweight+1;
                  indexweight[nlweight]:=nl;
               end
          end
      else
         begin
            counts[bin[j]]:=counts[bin[j]]+1;
            ysum[bin[j]]:=ysum[bin[j]]+y[ib];
            ysquaresum[bin[j]]:=ysquaresum[bin[j]]+Power(y[ib],2);
         end
    end;
end;

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

begin
  case kerneltype of
    1: begin
         for l:=0 to M-1 do
           kwe[l]:= 1;
       end;
    2: begin
         for l:=0 to M-1 do
           kwe[l]:=1 -(l/M);
        end;
    3: begin
         for l:=0 to M-1 do
           kwe[l]:= 1-Power((l/M),2);
       end;
    4: begin
         for l:=0 to M-1 do
           kwe[l]:=Power((1-Power((l/M),2)),2);
       end;
    5: begin
         for l:=0 to M-1 do
           kwe[l]:=Power((1-Power(l/M,2)),3);
       end;
     6: begin   {Gaussian kernel}
          for l:=0 to M-1 do
           kwe[l]:=exp(-8.0*(Power(l/M,2)));
        end;
    end;{case}
end;

procedure MainLoop;
var
jm,km,im :integer;
term1, term2 : real;

begin
   rm:=0;
   fm:=0;
   mm:=0;
   for M:=mstart to mend do
      begin
         CreateWeight;
         for jm:=1 to nlweight do
            begin
               km:=indexweight[jm];
               for im:=1-M to M-1 do
                   begin
                       index2:=bin[index[km]+im];
                       if index2>0 then
                           begin
                              rm:= rm+kwe[abs(im)]*ysum[index2];
                              fm:= fm+kwe[abs(im)]*counts[index2];
                           end; {endif}
                   end; {end im}
                mm:=rm/fm;
                term1:=ysquaresum[km]-2*mm*ysum[km]+counts[km]*Power(mm,2);
                term2:=SelecFun(kwe[0],fm);
                score[M-mstart+1]:=score[M-mstart+1]+term1*term2;
                rm:=0;
                fm:=0;
            end;  {endfor j}
         mv[M-mstart+1]:=M;
         hv[M-mstart+1]:=mv[M-mstart+1]*delta;
         score[M-mstart+1]:=score[M-mstart+1]/n;
   end;  {endfor M}
end; {Mainloop}

procedure ResulFile;
var
r1,rcount: integer;
begin


  assign(resul,'resfile');
  rewrite(resul);
  for rcount:=1 to mnumber do
    writeln(resul, mv[rcount]:5,score[rcount]:16:8,hv[rcount]:12:4);
  close(resul);
end;

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