program warpingc(input,output);
{(snp9: STB-27)
First written on November, 11, 1994; last correction November 18, 1994.
Authors: Isaias Hazarmabeth Salgado-Ugarte, Makoto Shimizu and
Toru Taniuchi. Reference: Stata Technical Bulletin XX, snp6.1

This version does not display results to screen, instead write them to
a text file in order to be used by Stata's command "infile" as a part
of three Stata's ado files (warpstep.ado, warpoly.ado and warping.ado).
This Pascal program performs Haerdle's WARPing density estimation
using modified algorithms from Haerdle (1991) "Smoothing Techniques
with Implementations in S. Springer-Verlag Series in Statistics
New York, and Scott (1992) "Multivariate Density Estimation: Theory,
Practice, and Visualization. John Wiley and Sons, New York.}

const
  maxsize = 2000;
  maxsize2 = 2000; {it is a really big limit!}

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

var
  binumber,n,i,M,numbin,kerneltype,index   : integer;
  start, delta,h,origin,maxval,minval,cm   : real;
  binmesh,counts,J                         : vectorint;
  x,wM,fh,midval                           : vector;
  inpval,data,resul                        : text;


procedure InputValues;
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,h,M,kerneltype);
   close(inpval);
end;

procedure InitialCalc;
var maxnoemp,ib: integer;

begin
   maxval:=x[n];
   minval:=x[1];
   if (kerneltype=6) then
      h:=4*h;
   delta:= h/M;
   numbin:= trunc((maxval-minval)/delta)+2*(M+1+round((M/10)+0.5));
   start:=(minval-h)-delta*0.1;
   if (start<0) then
      origin:=(round((start/delta)-0.5)-0.5)*delta
   else
      origin:=(trunc(start/delta)-0.5)*delta;
   if (n+1<numbin) then
      maxnoemp:=n+1
   else
      maxnoemp:=numbin;
   binumber:=0;
   index:=0;
   for ib:= 0 to numbin+1 do
      begin
        binmesh[ib]:=0;
        counts[ib]:=0;
        J[ib]:=0;
      end;
end;

procedure BinmeshCalc;
var
   ia : integer;

begin
   for i:=1 to n do
    begin
      index:= trunc((x[i]-origin)/delta);
      if (binmesh[index]=0) then
         begin
            binumber:= binumber+1;
            J[binumber]:=index;
            counts[binumber]:=1;
            binmesh[index]:=binumber;
          end
      else
        counts[binmesh[index]]:=counts[binmesh[index]]+1;
    end;

    for ia:=0 to numbin-1 do
         midval[ia]:=((0.5+ia)*delta)+origin;
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      {uniform kernel}
         cm:=M/((2*M-1)*(n*h));
         for l:=0 to M-1 do
           wM[l]:= cm;
       end;
    2: begin {Triangle kernel (Original Average-Shifted-Histogram kernel)}
         cm:=1/(n*h);
         for l:=0 to M-1 do
           wM[l]:= cm*(1 -l/M);
        end;
    3: begin     {Epanechnikov kernel}
         cm:=3*CleverPower(M,2)/((4*CleverPower(M,2)-1)*(n*h));
         for l:=0 to M-1 do
           wM[l]:= cm*(1-CleverPower((l/M),2));
       end;
    4: begin   {Quartic Kernel}
         cm:=0.9375/((1-0.0625/(CleverPower(M,2)))*(n*h));
         for l:=0 to M-1 do
           wM[l]:=cm*CleverPower(1-CleverPower(l/M,2),2);
       end;
    5: begin    {Triweight kernel}
         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
           wM[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
           wM[l]:=cm*exp(-8.0*(sqr(l/M)));
        end;
    end;{case}
end;

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

begin
  for iw:=0 to numbin-1 do
     fh[iw]:=0;
  for kw:=1 to binumber do
      for lw:=1-M to M-1 do
        fh[J[kw]+lw]:= fh[J[kw]+lw]+wM[abs(lw)]*counts[kw];
end;

procedure ResulFile;
var rcount: integer;
begin
  assign(resul,'resfile');
  rewrite(resul);
  for rcount:=0 to numbin-1 do
    writeln(resul,fh[rcount]:12:8,midval[rcount]:12:6);
  close(resul);
end;

begin {main}
    InputValues;
    InitialCalc;
    BinmeshCalc;
    CreateWeight;
    WeightBins;
    ResulFile;
end.