#include "libscl.h"
using namespace std;

/* 
Fortran documentation:

        subroutine sjeqd(n,x,lambda,hsjd)
C
C       From sish@agsm.unsw.edu.au Wed Mar 23 17:34 EST 1994
C       To: arg@eshprg.econ.unc.edu (Ron Gallant)
C       Date: Thu, 24 Mar 1994 08:25:55 +1000 (EET)
C
C       This subroutine calculates the Sheather-Jones
C       "solve-the-equation" bandwidth for a Gaussian kernel estimate.
C       (Journal of the Royal Statistical Society, Series B,
C       1991, Vol. 53, pp 683-690). 
C       Double precision is used throughout.
C       IMPORTANT INPUT VARIABLES
C       n = sample size
C       x = vector containing the data
C       lambda = interquartile range (iqr)
C       IMPORTANT OUTPUT VARIABLES
C       hsjd = value of Sheather-Jones bandwidth
C       OTHER IMPORTANT VARIABLES
C       hnml = bandwidth based on iqr assuming data is normal
C              (used as a starting value in this subroutine)
C       chat = estimate of constant c in bandwidth a
C       (i.e. a = c. h**(5/7))
C 
        real*8 x(1600),lambda,l2
        real*8 k1,k2,k3
        implicit real*8(a-h,o-z)
        integer n,loop
        parameter(tola=1.0d-5,tolg=1.0d-4,loop=15)
        xn = dble(n) 
        k3 = (xn - 1.d0)/dsqrt(2.d0)
        threen = 3.d0*xn
        pi = 3.141592654d0
        rt2pi = dsqrt(2.d0*pi)
        hnml = 0.79d0*lambda/(xn**0.2)
        l2=lambda*lambda
C       *** ESTIMATE VALUE OF c ***
        ac = (0.920d0*lambda)/(xn**(1.d0/7.d0))
        bc = (0.912d0*lambda)/(xn**(1.d0/9.d0))
        sa = 0.d0
        sb = 0.d0
        do 2 i = 1, n-1
          do 3 j = i+1, n
            da = ((x(i) - x(j))/ac)**2
             db = ((x(i) - x(j))/bc)**2
             ea = dexp(-0.5d0*da)
             eb = dexp(-0.5d0*db)
             sa = sa + (da**2 - 6.d0*da + 3.d0)*ea
             sb = sb + (db**3 - 15.d0*db**2 + 45.d0*db - 15.d0)*eb
3         continue
2       continue
        rhat2 = (2.d0*sa)/(xn*(xn-1.d0)*(ac**5)*rt2pi)
        rhat2 = rhat2 + 3.d0/(rt2pi*(xn-1.d0)*(ac**5))
        rhat3 = (-2.d0*sb)/(xn*(xn-1.d0)*(bc**7)*rt2pi)
        rhat3 = rhat3 + 15.d0/(rt2pi*(xn-1.d0)*(bc**7))
        chat = 1.357d0*((dabs(rhat2/rhat3))**(1.d0/7.d0))
        chat7 = chat**7.d0
C       *** USE NEWTON-RAPHSON METHOD TO CALCULATE h ***
C       INITIAL VALUES FOR h AND a
        firsth = hnml 
        ainit = chat*(firsth**(5.d0/7.d0))
99      a = 1.5d0*ainit
        do 4 m=1,loop
          s=0.d0
          t=0.d0
          do 5 i=1,n-1
            do 5 j=i+1,n
              y = (x(i)-x(j))/a
              y2 = y*y
              e = dexp(-0.5d0*y2)
              s = s + (y2*y2 - 6.d0*y2 + 3.d0)*e
              t = t + y2*(y2*y2-10.d0*y2+15.d0)*e
5           continue
          t = 2.0d0 * t
          s = 2.0d0 * s + threen
          k1 = a*a/chat7
          k1 = dsign(dabs(k1)**0.2d0,k1)
          k2 = k3/s
          k2 = dsign(dabs(k2)**0.2d0,k2)
          gprime = -0.2d0*(k2*(t/s-5.0d0) + 7.0d0*k1)
          g = a*(k2-k1)
          anew = a - g/gprime
C         CALCULATE h FROM a
          k1 = anew*anew/chat7
          k1 = dsign(dabs(k1)**0.2d0,k1)
          h=k1*anew
          if(anew.le.0.0d0)goto 98
          if(dabs(anew-a).lt.tola)goto100
          if(dabs(g).lt.tolg)goto100
          a=anew
4       continue
98      write(6,910)
        write(7,910)
        ainit = 1.5d0*ainit
        goto 99
        stop
100     continue
        hsjd = h
910     format(' NEWTON`S METHOD UNABLE TO FIND SOLUTION TO sjeqd')
        return
        end

Translated to C++ 5/4/94 by Ron Gallant.  The extra arguments nbin and
*w are to acccomdate binned input.  Returns 0 for success, returns 1 on 
failure.
*/

int sjeqd(INTEGER n, INTEGER nbin, REAL * x, REAL * w, REAL xdist,
            REAL lambda, REAL & hsjd)
{
  REAL const tola=1.0e-5;
  REAL const tolg=1.0e-4;

  INTEGER const loop=15;
  INTEGER const fail_limit = 25;

  INTEGER fail_count = 0;

  REAL const pi = 3.141592654;

  INTEGER i,j,m;

  REAL xn = n;

  REAL sswts = 0.0;
  for (i=0; i<nbin; i++) if (w[i]>1.0) sswts += (w[i]-1.0)*w[i]/2.0;

  REAL k1, k2;

  REAL k3 = (xn - 1.0)/sqrt(2.0);

  REAL threen = 3.0*xn;

  REAL rt2pi = sqrt(2.0*pi);

  REAL hnml = 0.79*lambda/pow(xn,0.2);

  //REAL l2 =lambda*lambda;

//*** estimate value of C ***

  REAL ac = (0.920*lambda)/(pow(xn,1.0/7.0));

  REAL bc = (0.912*lambda)/(pow(xn,1.0/9.0));

  REAL sa = 0.0;
  REAL sb = 0.0;
  REAL da, db, ea, eb;

  REAL xij, wij;

  for (i=1; i <= nbin-1; i++) {
    for (j=i+1; j <= nbin; j++) {
      xij = x[i-1] - x[j-1];
      wij = w[i-1]*w[j-1];
      da = (xij/ac)*(xij/ac);
      db = (xij/bc)*(xij/bc);
      ea = exp(-0.5*da);
      eb = exp(-0.5*db);
      sa += (da*da - 6.0*da + 3.0)*ea*wij;
      sb += (db*db*db - 15.0*db*db + 45.0*db - 15.0)*eb*wij;
    }
  }

  if (n != nbin) {
    da = (xdist/ac)*(xdist/ac);
    db = (xdist/bc)*(xdist/bc);
    ea = exp(-0.5*da);
    eb = exp(-0.5*db);
    sa += (da*da - 6.0*da + 3.0)*ea*sswts;
    sb += (db*db*db - 15.0*db*db + 45.0*db - 15.0)*eb*sswts;
  }

  REAL rhat2, rhat3, chat, chat7;

  rhat2 = (2.0*sa)/(xn*(xn-1.0)*pow(ac,5.0)*rt2pi);
  rhat2 += 3.0/(rt2pi*(xn-1.0)*pow(ac,5.0));
  rhat3 = (-2.0*sb)/(xn*(xn-1.0)*pow(bc,7.0)*rt2pi);
  rhat3 += 15.0/(rt2pi*(xn-1.0)*pow(bc,7.0));
  chat = 1.3570*pow(fabs(rhat2/rhat3),1.0/7.0);
  chat7 = pow(chat,7.0);


//*** USE NEWTON-RAPHSON METHOD TO CALCULATE h ***
//INITIAL VALUES FOR h AND a

  REAL firsth = hnml;
  REAL ainit = chat*pow(firsth,5.0/7.0);

  REAL a;

  REAL s, t, y, y2, e, gprime, g, anew, h;

L99:
   
  a = 1.5*ainit;

  for (m=1; m<=loop; m++) {
    s=0.0;
    t=0.0;
    for (i=1; i<=nbin-1; i++) {
      for (j=i+1; j<=nbin; j++) {
        xij = x[i-1] - x[j-1];
        wij = w[i-1]*w[j-1];
        y = xij/a;
        y2 = y*y;
        e = exp(-0.5*y2);
        s += wij*(y2*y2 - 6.0*y2 + 3.0)*e;
        t += wij*y2*(y2*y2 - 10.0*y2 + 15.0)*e;
      }
    }

    if (n != nbin) {
      y = xdist/a;
      y2 = y*y;
      e = exp(-0.5*y2);
      s += sswts*(y2*y2 - 6.0*y2 + 3.0)*e;
      t += sswts*y2*(y2*y2 - 10.0*y2 + 15.0)*e;
    }

    t *= 2.0;
    s *= 2.0;
    s += threen;
    k1 = a*a/chat7;
    k1 = (k1 >= 0.0 ) ? pow(fabs(k1),0.2) : -pow(fabs(k1),0.2);
    k2 = k3/s;
    k2 = (k2 >= 0.0 ) ? pow(fabs(k2),0.2) : -pow(fabs(k2),0.2);
    gprime = -0.2*(k2*(t/s-5.0) + 7.0*k1);
    g = a*(k2-k1);
    anew = a - g/gprime;

//  calculate h from a

    k1 = anew*anew/chat7;
    k1 = (k1 >= 0.0 ) ? pow(fabs(k1),0.2) : -pow(fabs(k1),0.2);
    h = k1*anew;
    if (anew <= 0.0) goto L98;
    if (fabs(anew-a) <= tola) goto L100;
    if (fabs(g) < tolg) goto L100;
    a=anew;
  }

L98:
    fail_count ++;
    if (fail_count > fail_limit) return 1;

    ainit *= 1.5;
    goto L99;

L100: 
    hsjd = h;
    return 0;
}







