/* ----------------------------------------------------------------------------

uniden.cpp

Kernel density estimator.  Self documented, execute uniden with no arguments 
for documentation.

Copyright (C) 1991, 1992, 1993, 1994, 1999, 2002, 2006.

A. Ronald Gallant
Post Office Box 659 
Chapel Hill NC 27514-0659 
USA   

Permission to use, copy, modify, and distribute this software and its docu-
mentation for any purpose and without fee is hereby granted, provided that 
the above copyright notice appear in all copies and that both that copyright 
notice and this permission notice appear in supporting documentation.  This 
software is provided "as is" without any expressed or implied warranty.

-----------------------------------------------------------------------------*/

#include "libscl.h"
using namespace std;
using namespace scl;


INTEGER max_sample_size=1024;   // Can be changed with flag -m
INTEGER grid=100;               // Can be changed with flag -g
REAL    bandwidth;              // Set with flag -b, defaults otherwise.
REAL    xmin;                   // Set with flag -xmin, defaults otherwise.
REAL    xmax;                   // Set with flag -xmax, defaults otherwise.
int     predict = 0;            // The flag -p sets to true, estimate otw.
int     eformat = 0;            // The flag -e sets to true
int     cross_validate = 0;     // The flag -c sets to true, estimate otw.
int     normal_kernel = 1;      // The flag -q or -n sets to false.
int     quadratic_kernel = 0;   // The flag -q sets to true.
int     normal_plot = 0;        // The flag -n sets to true.
int     filter = 0;             // The flag -f sets to true.
int     Sheather_Jones = 1;     // The flag -o sets to false.
int     do_not_bin = 0;         // The flag -d sets to true.
int     compute_bandwidth = 1;
int     compute_xmax = 1;              
int     compute_xmin = 1;

INTEGER  n, nbin;
REAL     *Xsrt;
REAL     *Xsrt_top;
REAL     *Xbin;
REAL     *Xbin_top;
REAL     *Xwgt;
REAL     xmean,xsigma,xspread,xdist;

void usage_msg();
void process_file(istream& cfile);
void pred();
void estimate();
void bin(REAL b);
REAL crossval(int isw);
REAL normal(REAL z, int isw);
REAL quadratic(REAL z, int isw);
int  sjeqd(INTEGER n, INTEGER nbin, REAL *x, REAL *w, REAL xdist,
             REAL lambda, REAL & hsjd);

REAL const zero = (REAL) 0.0;

char    *progname;

void usage_msg()
{
  cerr << '\n'
    << "usage: " <<progname<< " [options] filename \n\n"
    << "options: -p             predict at data points only \n"
    << "         -e             output observations in e format with full precision\n"
    << "         -c             cross validate instead of estimate \n" 
    << "         -q             quadratic kernel, normal is default\n"
    << "         -o             override Sheather-Jones bandwidth selection\n"
    << "         -d             do not bin the data\n"
    << "         -f             acts as a filter, takes standard input\n"
    << "         -n             normal density points instead of kernel\n"
    << "         -m     value   maximum sample size \n" 
    << "         -b     value   bandwidth\n"
    << "         -g     value   number of plot points sent to cout\n"
    << "         -xmin  value   minimum plot point\n"
    << "         -xmax  value   maximum plot point\n\n"
    << "example: "<<progname<<" -b 6.0 -q -g 200 in.dat >out.dat\n\n"
    << "The input file contains x values, one per line.  The output\n"
    << "file contains y x pairs, one pair per line, y first.\n\n";
}

int main(int argc, char *argp[], char *envp[])
{
  int      argx, flagx, filex;
  char     ch, flag;
  REAL     cv, cv_min, bandwidth_min, decrement, tmp, A, silverman;
  int      sj_failure;

  ch = envp[0][0]; ch++; // This supresses a compiler warning.

  progname = argp[0];

  // begin: loop over flags and input files.

  argx = 0;  filex = 1;  

  while (argx < argc) {

    if (argx==0) argx++;

    // begin: read flags, if any.

    while (argx < argc && (argp[argx])[0]=='-') {

      // Read flags, if any.  Note that flags can occur singly as -a -b 
      // or they can be grouped as -ab.  Just insert code within the
      // appropriate case and don't worry about it.

      for (flagx=1; (flag=(argp[argx])[flagx])!=0; flagx++) {

        // If a flag can have more than one letter, deal with it within
        // the appropriate case by stripping off letters and incrementing
        // flagx.

        // If a flag indicates that argp[argx+1] is a file that should be read 
        // then, within the appropriate case, open arg[argx+1] or otherwise 
        // deal with the file name arg[argx+1], assign argp[argx] to 
        // argp[argx+1], and increment argx.  

        // If a flag indicates that argp[argx+1] is a string containing the  
        // value of a parameter, etc., then, within the appropriate case, deal 
        // with the string arg[argx+1], assign argp[argx] to argp[argx+1], and 
        // increment argx.  

        switch (flag) {

        case 'c' :
          cross_validate = 1;
          break;

        case 'o' :
          Sheather_Jones = 0;
          break;

        case 'd' :
          do_not_bin = 1;
          break;

        case 'p' :
          predict = 1;
          break;

        case 'e' :
          eformat = 1;
          break;

        case 'q' :
          normal_plot = 0;
          normal_kernel = 0;
          quadratic_kernel = 1;
          break;

        case 'f' :
          filter = 1;
          break;

        case 'n' :
          normal_plot = 1;
          normal_kernel = 0;
          quadratic_kernel = 0;
          break;

        case 'b' :
          bandwidth = (REAL) atof(argp[argx+1]); 
          compute_bandwidth = 0;
          argp[argx+1] = argp[argx];
          argx++;
          break;

        case 'm' :
          max_sample_size = (INTEGER) atoi(argp[argx+1]); 
          argp[argx+1] = argp[argx];
          argx++;
          break;

        case 'g' :
          grid = (INTEGER) atoi(argp[argx+1]); 
          argp[argx+1] = argp[argx];
          argx++;
          break;

        case 'x' : 
          if (memcmp(&(argp[argx][flagx]),"xmin",4)==0) {
            flagx += 3;
            xmin = atof(argp[argx+1]); 
            argp[argx+1] = argp[argx];
            compute_xmin = 0;
            argx++;
            break;
          }
          if (memcmp(&(argp[argx][flagx]),"xmax",4)==0) {
            flagx += 3;
            xmax = atof(argp[argx+1]); 
            argp[argx+1] = argp[argx];
            compute_xmax = 0;
            argx++;
            break;
          }
          goto badflag;

        badflag :
        default :
          cerr << "Flag was ";
          while ((flag=(argp[argx])[flagx++])!=0) 
            cerr << flag;
          cerr << "\nError: unknown flag.\n";
          exit(1);
          break;

        } // end: switch

      } // end: for 

    argx++;

    } // end: while

    // end: read flags, if any.

    if (filter) {
      process_file(cin);
    }
    else {
      if (argx >= argc) {
        usage_msg();
        exit(1);
      }

      ifstream infile(argp[argx]);
      if (!infile) {
        cerr << "Error, " << progname << ", Cannot open " << argp[argx] << "\n";
        usage_msg();
        exit(1);
      } 
      istream& from = infile;

      process_file(from);
      infile.close();
    }

    if (normal_plot) {
      if (predict) {
        pred();
      }
      else {
       estimate();
      }
      exit(0);
    }

    if (compute_bandwidth) {

      A = (xsigma < xspread) ? xsigma : xspread;
      silverman = bandwidth = 0.9*A/pow((double)n,0.2);

      decrement = bandwidth/10.0;
      bandwidth = silverman + 5.0*decrement;
      bandwidth_min = silverman - 7.001*decrement;

      bin(bandwidth_min/3.1);

      sj_failure = 1;

      if (Sheather_Jones && normal_kernel) {
         sj_failure = sjeqd(n, nbin, Xbin, Xwgt, xdist, 1.34*xspread, tmp);

         if (tmp < bandwidth_min/3.1) {
            sj_failure = 1;
            cerr << "\nSheather-Jones bandwidth ("
                 << tmp << ") is less than the bin width ("
                 << bandwidth_min/3.1 << ").\n";
            cerr << "\nUse -b option to get very small bandwidths.\n";
         }
         if (sj_failure) {
            cerr << "\nSheather-Jones failed, using cross validation.\n";
         }
         else {
            bandwidth = tmp;
            cerr <<progname<< ": Sheather-Jones bandwidth = "<<bandwidth<<'\n';
         }
      }

      if (quadratic_kernel || sj_failure) {
        tmp=bandwidth;
        cv_min = crossval(0);
        cerr << "\nTrials:\n";
        cerr << "cv " << cv_min << "  bw " << bandwidth << "\n";

        do{
          bandwidth -= decrement; 
          cv = crossval(0);
          if (fabs(bandwidth-silverman)<0.25*decrement) {  
            cerr << "cv " << cv <<"  bw " <<bandwidth
                << " = 0.9*min(sd,iqr/1.34)/pow(n,0.2)\n";
          }
          else {
            cerr << "cv " << cv << "  bw " << bandwidth << "\n";
          }
          if (cv < cv_min) { 
            cv_min = cv; 
            tmp = bandwidth;
          }
        } while (bandwidth > bandwidth_min);

        bandwidth = tmp;
        cerr << "\nSelection:\n";
        cerr << "bandwidth = " << bandwidth << "\n\n";
      }
    }
    else {
      bin(bandwidth/3.1);
    }

    if (cross_validate) {
      cv = crossval(1);
      exit(0);
    }

    if (predict) {
      pred();
      exit(0);
    }

    estimate();
    exit(0);

    filex++;
    argx++;
  } 

  // end: loop over flags and input files.

  exit(0);

}

void process_file(istream& from)
{ 
  REAL     *x;
  INTEGER  lo, hi;
  void usage_msg();

  REAL     *X;
  REAL     *X_top;

  X = new REAL[max_sample_size+1];

  n = 0;
  X_top=X;

  while( (from>>*X_top) && (n <= max_sample_size) ) { 
    X_top++; n++;
  }

  if ( n > max_sample_size ) {
    cerr << "\nError: The maximum sample size = " 
         << max_sample_size 
         << " is too small.\n"
         << "       Give a larger value with the -m flag.\n\n";
    usage_msg();
    exit(1);
  }

  Xsrt = new REAL[n];
  Xsrt_top = Xsrt + n;

  heapsort(X,n,Xsrt);

  lo = n/4;  hi = (3*n)/4;
  xspread = *(Xsrt+hi) - *(Xsrt+lo); 
  if (xspread<=0.0) xspread = *(Xsrt+n-1) - *Xsrt;
  xspread /= 1.34;

  if (compute_xmin) xmin=*Xsrt;
  if (compute_xmax) xmax=*(Xsrt+n-1);
  
  for (x=X, xmean=0.0; x < X_top; x++) {
    xmean += *x;
  }
  xmean /= (REAL)n;
  for (x=X, xsigma=0.0; x < X_top; x++) {
    xsigma += (*x-xmean)*(*x-xmean);
  }
  xsigma = sqrt(xsigma/(REAL)(n-1));  

  delete [] X; 

}

REAL crossval(int isw)
{
  REAL     *xs,   *xt;
  REAL     *ws,   *wt;
  REAL     u, sum, cv;

  u=sum=cv=0; //This stops a compiler warning.

  if (normal_kernel) {

    sum = zero;
    for (xt=Xbin, wt=Xwgt; xt < Xbin_top; xt++, wt++) {
    for (xs=xt,   ws=wt;   xs < Xbin_top; xs++, ws++) {
      if (xt != xs) {
        u = (*xt-*xs)/bandwidth;
        sum += 2.0*(*wt)*(*ws) * (normal(u,2) - 2.0*normal(u,1));
      }
      else {
        u = 0.0;
        sum += (*wt) * (normal(u,2) - 2.0*normal(u,1));
        u =(xdist)/bandwidth;
        if (*wt>1.0) sum += (*wt-1.0)*(*wt) * (normal(u,2) - 2.0*normal(u,1));
      }
    }
    }
    sum /= (REAL)n;
    sum += 2.0*normal(zero,1);
    sum /= bandwidth*(REAL)n;
    cv = sum;
  }

  if (quadratic_kernel) {
  // Fix this. It's inefficient.  There's no need to loop over 
  // the entire sample.

    sum = zero;
    for (xt=Xbin, wt=Xwgt; xt < Xbin_top; xt++, wt++) {
    for (xs=xt,   ws=wt;   xs < Xbin_top; xs++, ws++) {
      if (xt != xs) {
        u = (*xt-*xs)/bandwidth;
        sum += 2.0*(*wt)*(*ws) * (quadratic(u,2) - 2.0*quadratic(u,1));
      }
      else {
        u = 0.0;
        sum += (*wt) * (quadratic(u,2) - 2.0*quadratic(u,1));
        u =(xdist)/bandwidth;
        if (*wt>1.0) sum += (*wt-1.0)*(*wt)*(quadratic(u,2)-2.0*quadratic(u,1));
      }
    }
    }
    sum /= (REAL)n;
    sum += 2.0*quadratic(zero,1);
    sum /= bandwidth*(REAL)n;
    cv = sum;
  }

  if (isw) {
    cout << "Cross validation = " << cv 
         << " for bandwidth = " << bandwidth << "\n";
  }

  return cv;
}

void estimate()
{
  REAL     *x, *w;
  REAL     yplt, xplt;
  REAL     sum, weight, u;
  REAL     xinc;
  INTEGER  j;

  xinc = (xmax - xmin) / ( (REAL) (grid-1) );

  if (normal_kernel) {

    char str[80];

    for (j=0; j<grid; j++) {
      xplt = xmin + xinc*j;

      for (sum=zero, x=Xbin, w=Xwgt; x < Xbin_top; x++, w++) {
        u = (xplt-*x)/bandwidth;
        weight = (*w) * (normal(u,1)/bandwidth);
        sum += weight;
      }

      yplt = sum/(REAL)n;

      if (eformat) {
        sprintf(str,"%30.18e %30.18e\n",yplt,xplt); cout << str;
      }
      else {
        cout << yplt << ' ' << xplt << '\n';
      }
    }
  }
  
  if (normal_plot) {

    char str[80];

    for (j=0; j<grid; j++) {
      xplt = xmin + xinc*j;

      u = (xplt-xmean)/xsigma;

      yplt = normal(u,1)/xsigma;

      if (eformat) {
        sprintf(str,"%30.18e %30.18e\n",yplt,xplt); cout << str;
      }
      else {
        cout << yplt << ' ' << xplt << '\n';
      }
    }
  }

  if (quadratic_kernel) {
  // Fix this. It's inefficient.  There's no need to loop over 
  // the entire sample.

    char str[80];

    for (j=0; j<grid; j++) {
      xplt = xmin + xinc*j;

      for (sum=zero, x=Xbin, w=Xwgt; x < Xbin_top; x++, w++) {
        u = (xplt-*x)/bandwidth;
        weight = (*w) * (quadratic(u,1)/bandwidth);
        sum += weight;
      }

      yplt = sum/(REAL)n;

      if (eformat) {
        sprintf(str,"%30.18e %30.18e\n",yplt,xplt); cout << str;
      }
      else {
        cout << yplt << ' ' << xplt << '\n';
      }
    }
  }
}

void pred()
{
  REAL     *x, *xt, *w;
  REAL     yplt;
  REAL     sum, weight, u;

  if (normal_kernel) {

    char str[80];

    for (xt=Xsrt; xt < Xsrt_top; xt++) {

      for (sum=zero, x=Xbin, w=Xwgt; x < Xbin_top; x++, w++) {
        u = (*xt-*x)/bandwidth;
        weight = (*w) * (normal(u,1)/bandwidth);
        sum += weight;
      }

      yplt =sum/(REAL)n;

      if (eformat) {
        sprintf(str,"%30.18e %30.18e\n",yplt,*xt); cout << str;
      }
      else {
        cout << yplt << ' ' << *xt << '\n';
      }
    }
  }

  if (normal_plot) {

    char str[80];

    for (xt=Xsrt; xt < Xsrt_top; xt++) {
      u = (*xt-xmean)/xsigma;
      yplt = normal(u,1)/xsigma;

      if (eformat) {
        sprintf(str,"%30.18e %30.18e\n",yplt,*xt); cout << str;
      }
      else {
        cout << yplt << ' ' << *xt << '\n';
      }
    }
  }

  if (quadratic_kernel) {
  // Fix this. It's inefficient.  There's no need to loop over 
  // the entire sample.

    char str[80];

    for (xt=Xsrt; xt < Xsrt_top; xt++) {

      for (sum=zero, x=Xbin, w=Xwgt; x < Xbin_top; x++, w++) {
        u = (*xt-*x)/bandwidth;
        weight = (*w) * (quadratic(u,1)/bandwidth);
        sum += weight;
      }
    
      yplt = sum/(REAL)n;

      if (eformat) {
        sprintf(str,"%30.18e %30.18e\n",yplt,*xt); cout << str;
      }
      else {
        cout << yplt << ' ' << *xt << '\n';
      }
    }
  }
}

REAL normal(REAL z, int isw)
// Returns the value of the normal kernel at z if isw=1 and its 
// convolution with itself if isw=2.
{

  double const root_two_pi  = 2.5066282746310005024;
  double const root_four_pi = 3.5449077018110326520;
  REAL   const zero = 0.0;

  if ( (z < -6.0) || (z > 6.0) ) return zero;

  double       h;
  double       tmp;
  double       zz;

  zz = (double)z;

  switch (isw) {

    case 1 :
     tmp = -zz*zz/2.0;
     h = exp( tmp );
     h /= root_two_pi;
     return (REAL)h;

    case 2 :
     tmp = -zz*zz/4.0;
     h =  exp( tmp );
     h /= root_four_pi;
     return (REAL)h;

    default :
      cerr << "\nError: bad data to REAL normal(REAL z, int isw)\n";
      exit(1);
   }
   // These two statements can't be reached but stop a compiler warning.
   cerr << "\nError: bad data to REAL normal(REAL z, int isw)\n";
   return 1;
}    

REAL quadratic(REAL z, int isw)
// Returns the value of the quadratic kernel at z if isw=1 and its 
// convolution with itself if isw=2.
{
  REAL const three_fourths = 0.75;
  REAL const nine_sixteenths = 0.5625;
  REAL const one = 1.0;
  REAL       h;
  REAL       top,bot,sumtop,sumbot;

  switch (isw) {

    case 1 :
      if (-one <= z && z <= one) {
        h = three_fourths*(one - z*z);
      }
      else {
        h = zero;
      }
      return h;

    case 2 :

      if (z > 2.0 || z < -2.0 ) 
        return zero;;

      if (z > zero) {
        top = one;
        bot = z - one;
      }
      else {
        top = z + one;
        bot = -one;
      }

      sumtop =  zero;
      sumtop += one/5.0;
      sumtop *= top;
      sumtop += (-z/2.0);
      sumtop *= top;
      sumtop += (z*z-2.0)/3.0;
      sumtop *= top;
      sumtop += z;
      sumtop *= top;
      sumtop += one-z*z;
      sumtop *= top;

      sumbot =  zero;
      sumbot += one/5.0;
      sumbot *= bot;
      sumbot += (-z/2.0);
      sumbot *= bot;
      sumbot += (z*z-2.0)/3.0;
      sumbot *= bot;
      sumbot += z;
      sumbot *= bot;
      sumbot += one-z*z;
      sumbot *= bot;

      h = nine_sixteenths*(sumtop-sumbot);
      return h;

    default :
      cerr << "\nError: bad data to REAL quadratic(REAL z, int isw)\n";
      exit(1);
  }
  // These two statements can't be reached but stop a compiler warning.
  cerr << "\nError: bad data to REAL normal(REAL z, int isw)\n";
  return 1;
}

void bin(REAL b) 
{
  const REAL zero = 0.0;
  const REAL one = 1.0;

  INTEGER i,j;

  if ((n <= 50)||(do_not_bin)) {
    Xbin = new REAL[n];
    Xwgt = new REAL[n];
    Xbin_top = Xbin + n;
    for (i=0; i<n; i++) {
      Xbin[i]=Xsrt[i];
      Xwgt[i]=one;
    }
    nbin = n;
    xdist = 0;
    return;
  }

  double fn = (double)n;
  double sn = sqrt(fn);
  double qn = 0.1*fn;
  double d;

  d = (sn < qn) ? sn : qn;
  i = (INTEGER) d;

  INTEGER q10 = ( 0 > i ) ? 0 : i;
  INTEGER q90 = ( n-i < n-1 ) ? n-i : n-1;

  // bins:                 
  // left                         0,    ...,  q10
  // middle                 q10 + 1,    ...,  q10 + nbins_middle,
  // right   q10 + nbins_middle + 1,    ...,  q10 + nbins_middle + n - q90

  INTEGER nbins_middle = (INTEGER) ((Xsrt[q90]-Xsrt[q10])/b);
  INTEGER nbins_total = nbins_middle + q10 + 1 + n - q90;
  INTEGER bin_count = 0;

  Xbin = new REAL[nbins_total];
  Xwgt = new REAL[nbins_total];
  
  REAL tmp;
  tmp = nbins_middle;
  REAL xinc = (Xsrt[q90]-Xsrt[q10])/tmp;

  for (i=0; i<=q10; i++) {
    Xbin[i] = Xsrt[i];
    Xwgt[i] = one;
    bin_count++;
  }

  tmp = Xsrt[q10];

  for (i=q10+1; i<=q10 + nbins_middle; i++) {
    tmp += xinc;
    Xbin[i] = tmp;
    Xwgt[i] = zero;
    bin_count++;
  }

  for (i=q90; i<n; i++) {
    Xbin[q10 + nbins_middle + 1 + i - q90] = Xsrt[i];
    Xwgt[q10 + nbins_middle + 1 + i - q90] = one;
    bin_count++;
  }

  Xwgt[q10 + nbins_middle + 1] = zero;
  REAL * x = Xsrt;
  REAL low, high;
  for (i=q10+1; i<=q10 + nbins_middle + 1; i++) {
    low = Xbin[i-1];
    high = Xbin[i];
    while (*x<=high) {
      if( (low < *x) && (*x <= high) ) Xwgt[i]++;
      x++;
    }
  }

  if (nbins_total != bin_count) {
    cerr << "Error, "<<progname<< ", bin, this should never happen.";
    exit(1);
  }

//for (i=0; i<bin_count; i++) cout<< i <<' '<< Xbin[i] <<' '<< Xwgt[i] <<'\n';

  tmp = 0.5*xinc;
  for (i=q10+1; i<=q10 + nbins_middle; i++) {
    Xbin[i] -= tmp;
  }

  for (i=0, j=0; i<bin_count; i++) {
    if (Xwgt[i] > 0 ) {
      Xbin[j] = Xbin[i];
      Xwgt[j] = Xwgt[i];
      j++;
    }
  }

  bin_count = j;

  Xbin_top = Xbin + bin_count;

  nbin = bin_count;

  xdist = 0.5*xinc;

//for (i=0; i<bin_count; i++) cout << i <<' '<< Xbin[i] <<' '<< Xwgt[i] <<'\n';
}

