#include "spline.h"


extern"C" {
	void sgram_(double *sg0, double *sg1, double *sg2,
				double *sg3, double* knot, int* nk);
	void stxwx_(double *xs, double *ys, double *ws, int *n,
				double *knot, int *nk, double *xwy, double *hs0, double *hs1, double *hs2, double *hs3);
	void sslvrg_(double *penalt, double *dofoff, double *xs, double *ys, double *ws, double *ssw, int *n,
				 double *knot, int *nk,
				 double *coef, double *sz, double *lev, double *crit, int *icrit, double *lspar, double *xwy,
				 double *hs0, double *hs1, double *hs2, double *hs3,
				 double *sg0, double *sg1, double *sg2, double *sg3, double *abd,
				 double *p1ip, double *p2ip, int *ld4, int *ldnk, int *ier);
}
/***********************************************************/
// used as reference - http://www.koders.com/fortran/fid8AA63B49CF22F0138E9B3DBDC405696F4A62C1CF.aspx
//  http://www.koders.com/c/fidD995301A8A5549CE0361F4E7FFDFD3CDC4B4E4A3.aspx
/* A Cubic B-spline Smoothing routine.
 
 // sbart.f -- translated by f2c (version 20010821).
 // ------- and f2c-clean,v 1.9 2000/01/13
 //
 // According to the GAMFIT sources, this was derived from code by
 // Finbarr O'Sullivan.
 //
 
 
 The algorithm minimises:
 
 (1/n) * sum ws(i)^2 * (ys(i)-sz(i))^2 + lambda* int ( s"(x) )^2 dx
 
 lambda is a function of the spar which is assumed to be between 0 and 1
 
 INPUT
 -----
 penalt			A penalty > 1 to be used in the gcv criterion
 dofoff			either `df.offset' for GCV or `df' (to be matched).
 n				number of data points
 ys(n)			vector of length n containing the observations
 ws(n)			vector containing the weights given to each data point
 NB: the code alters the values here.
 xs(n)			vector containing the ordinates of the observations
 ssw			`centered weighted sum of y^2'
 nk				number of b-spline coefficients to be estimated
 nk <= n+2
 knot(nk+4)		vector of knot points defining the cubic b-spline basis.
 To obtain full cubic smoothing splines one might
 have (provided the xs-values are strictly increasing)
 spar			penalised likelihood smoothing parameter
 ispar			indicating if spar is supplied (ispar=1) or to be estimated
 lspar, uspar	lower and upper values for spar search;  0.,1. are good values
 tol, eps		used in Golden Search routine
 isetup			setup indicator [initially 0
 icrit			indicator saying which cross validation score is to be computed
 0: none ;  1: GCV ;  2: CV ;  3: 'df matching'
 ld4			the leading dimension of abd (ie ld4=4)
 ldnk			the leading dimension of p2ip (not referenced)
 
 OUTPUT
 ------
 coef(nk)		vector of spline coefficients
 sz(n)			vector of smoothed z-values
 lev(n)			vector of leverages
 crit			either ordinary or generalized CV score
 spar			if ispar != 1
 lspar			== lambda (a function of spar and the design)
 iter			number of iterations needed for spar search (if ispar != 1)
 ier			error indicator
 ier = 0 ___  everything fine
 ier = 1 ___  spar too small or too big
 problem in cholesky decomposition
 
 Working arrays/matrix
 xwy				X'Wy
 hs0,hs1,hs2,hs3	the diagonals of the X'WX matrix
 sg0,sg1,sg2,sg3	the diagonals of the Gram matrix SIGMA
 abd (ld4,nk)		[ X'WX + lambda*SIGMA ] in diagonal form
 p1ip(ld4,nk)		inner products between columns of L inverse
 p2ip(ldnk,nk)		all inner products between columns of L inverse
 where  L'L = [X'WX + lambda*SIGMA]  NOT REFERENCED
 */

int Spline::sbart
(double *penalt, double *dofoff,
 double *xs, double *ys, double *ws, double *ssw,
 int *n, double *knot, int *nk, double *coef,
 double *sz, double *lev, double *crit, int *icrit,
 double *spar, int *ispar, int *iter, double *lspar,
 double *uspar, double *tol, double *eps, int *isetup,
 int *ld4, int *ldnk, int *ier)
{
	try{
	/* A Cubic B-spline Smoothing routine.
	 
	 The algorithm minimises:
	 
	 (1/n) * sum ws(i)^2 * (ys(i)-sz(i))^2 + lambda* int ( s(x) )^2 dx
	 
	 lambda is a function of the spar which is assumed to be between 0 and 1
	 
	 INPUT
	 -----
	 penalt	A penalty > 1 to be used in the gcv criterion
	 dofoff	either `df.offset' for GCV or `df' (to be matched).
	 n		number of data points
	 ys(n)	vector of length n containing the observations
	 ws(n)	vector containing the weights given to each data point
	 NB: the code alters the values here.
	 xs(n)	vector containing the ordinates of the observations
	 ssw          `centered weighted sum of y^2
	 nk		number of b-spline coefficients to be estimated
	 nk <= n+2
	 knot(nk+4)	vector of knot points defining the cubic b-spline basis.
	 To obtain full cubic smoothing splines one might
	 have (provided the xs-values are strictly increasing)
	 spar		penalised likelihood smoothing parameter
	 ispar	indicating if spar is supplied (ispar=1) or to be estimated
	 lspar, uspar lower and upper values for spar search;  0.,1. are good values
	 tol, eps	used in Golden Search routine
	 isetup	setup indicator [initially 0
	 icrit	indicator saying which cross validation score is to be computed
	 0: none ;  1: GCV ;  2: CV ;  3: 'df matching'
	 ld4		the leading dimension of abd (ie ld4=4)
	 ldnk		the leading dimension of p2ip (not referenced)
	 
	 OUTPUT
	 ------
	 coef(nk)	vector of spline coefficients
	 sz(n)	vector of smoothed z-values
	 lev(n)	vector of leverages
	 crit		either ordinary or generalized CV score
	 spar         if ispar != 1
	 lspar         == lambda (a function of spar and the design)
	 iter		number of iterations needed for spar search (if ispar != 1)
	 ier		error indicator
	 ier = 0 ___  everything fine
	 ier = 1 ___  spar too small or too big
	 problem in cholesky decomposition
	 
	 Working arrays/matrix
	 xwy			XWy
	 hs0,hs1,hs2,hs3	the diagonals of the XWX matrix
	 sg0,sg1,sg2,sg3	the diagonals of the Gram matrix SIGMA
	 abd (ld4,nk)		[ XWX + lambda*SIGMA ] in diagonal form
	 p1ip(ld4,nk)		inner products between columns of L inverse
	 p2ip(ldnk,nk)	all inner products between columns of L inverse
	 where  LL = [XWX + lambda*SIGMA]  NOT REFERENCED
	 */
	
#define CRIT(FX) (*icrit == 3 ? FX - 3. : FX)
		/* cancellation in (3 + eps) - 3, but still...informative */
		
#define BIG_f (1e100)
		
		/* c_Gold is the squared inverse of the golden ratio */
		static const double c_Gold = 0.381966011250105151795413165634;
		/* == (3. - sqrt(5.)) / 2. */
		
		/* Local variables */
		static double ratio;/* must be static (not needed in R) */
		
		double a, b, d, e, p, q, r, u, v, w, x;
		double ax, fu, fv, fw, fx, bx, xm;
		double t1, t2, tol1, tol2;
		double* xwy = new double[*nk];
		double* hs0 = new double[*nk];
		double* hs1 = new double[*nk];
		double* hs2 = new double[*nk];
		double* hs3 = new double[*nk];
		double* sg0 = new double[*nk];
		double* sg1 = new double[*nk];
		double* sg2 = new double[*nk];
		double* sg3 = new double[*nk];
		double* abd = new double[*nk*(*ld4)]; 
		double* p1ip = new double[*nk*(*ld4)]; 
		double* p2ip = new double[*nk];	
		int i, maxit;

		/* unnecessary initializations to keep  -Wall happy */
		d = 0.; fu = 0.; u = 0.;
		ratio = 1.;
		
		/*  Compute SIGMA, X' W X, X' W z, trace ratio, s0, s1.
		 
		 SIGMA	-> sg0,sg1,sg2,sg3
		 X' W X	-> hs0,hs1,hs2,hs3
		 X' W Z	-> xwy
		 */
		
		/* trevor fixed this 4/19/88
		 * Note: sbart, i.e. stxwx() and sslvrg() {mostly, not always!}, use
		 *	 the square of the weights; the following rectifies that */
		for (i = 0; i < *n; ++i)
			if (ws[i] > 0.)
				ws[i] = sqrt(ws[i]);
		
		if (*isetup == 0) {
			/* SIGMA[i,j] := Int  B''(i,t) B''(j,t) dt  {B(k,.) = k-th B-spline} */
			sgram_(sg0, sg1, sg2, sg3, knot, nk);
			stxwx_(xs, ys, ws, n,
							knot, nk,
							xwy,
							hs0, hs1, hs2, hs3);
			/* Compute ratio :=  tr(X' W X) / tr(SIGMA) */
			t1 = t2 = 0.;
			for (i = 3 - 1; i < (*nk - 3); ++i) {
				t1 += hs0[i];
				t2 += sg0[i];
			}
			ratio = t1 / t2;
			*isetup = 1;
		}
		/*     Compute estimate */
		
		if (*ispar == 1) { /* Value of spar supplied */
			*lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0));
			sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n,
							 knot, nk,
							 coef, sz, lev, crit, icrit, lspar, xwy,
							 hs0, hs1, hs2, hs3,
							 sg0, sg1, sg2, sg3, abd,
							 p1ip, p2ip, ld4, ldnk, ier);
			/* got through check 2 */
			return 0;
		}
		
		/* ELSE ---- spar not supplied --> compute it ! ---------------------------
		 
		 Use Forsythe Malcom and Moler routine to MINIMIZE criterion
		 f denotes the value of the criterion
		 
		 an approximation	x  to the point where	f  attains a minimum  on
		 the interval  (ax,bx)  is determined.
		 */
		ax = *lspar;
		bx = *uspar;
		
		/* INPUT
		 
		 ax	 left endpoint of initial interval
		 bx	 right endpoint of initial interval
		 f	 function subprogram which evaluates  f(x)  for any  x
         in the interval  (ax,bx)
		 tol	 desired length of the interval of uncertainty of the final
         result ( >= 0 )
		 
		 OUTPUT
		 
		 fmin	 abcissa approximating the point where	f  attains a minimum
		 */
		
		/*
		 The method used is a combination of  golden  section  search  and
		 successive parabolic interpolation.	convergence is never much slower
		 than	 that  for  a  fibonacci search.  if  f	 has a continuous second
		 derivative which is positive at the minimum (which is not  at  ax  or
		 bx),	 then  convergence  is	superlinear, and usually of the order of
		 about  1.324....
		 the function  f  is never evaluated at two points closer together
		 than	 eps*abs(fmin) + (tol/3), where eps is	approximately the square
		 root	 of  the  relative  machine  precision.	  if   f   is a unimodal
		 function and the computed values of	 f   are  always  unimodal  when
		 separated by at least  eps*abs(x) + (tol/3), then  fmin  approximates
		 the abcissa of the global minimum of	 f  on the interval  ax,bx  with
		 an error less than  3*eps*abs(fmin) + tol.  if   f	is not unimodal,
		 then fmin may approximate a local, but perhaps non-global, minimum to
		 the same accuracy.
		 this function subprogram is a slightly modified	version	 of  the
		 algol  60 procedure	localmin  given in richard brent, algorithms for
		 minimization without derivatives, prentice - hall, inc. (1973).
		 
		 Double	 a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w
		 Double	 fu,fv,fw,fx,x
		 */
		
		/*  eps is approximately the square root of the relative machine
		 precision.
		 
		 -	 eps = 1e0
		 - 10	 eps = eps/2e0
		 -	 tol1 = 1e0 + eps
		 -	 if (tol1 > 1e0) go to 10
		 -	 eps = sqrt(eps)
		 R Version <= 1.3.x had
		 eps = .000244     ( = sqrt(5.954 e-8) )
		 -- now eps is passed as argument
		 */
		
		/* initialization */
		
		maxit = *iter;
		*iter = 0;
		a = ax;
		b = bx;
		v = a + c_Gold * (b - a);
		w = v;
		x = v;
		e = 0.;
		*spar = x;
		*lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0));
		sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n,
						 knot, nk,
						 coef, sz, lev, crit, icrit, lspar, xwy,
				hs0, hs1, hs2, hs3,
				sg0, sg1, sg2, sg3, abd,
				p1ip, p2ip, ld4, ldnk, ier);
		fx = *crit;
		fv = fx;
		fw = fx;
		
		/* main loop
		 --------- */
		while(*ier == 0) { /* L20: */
			
			if (m->control_pressed) { return 0; }
			
			xm = (a + b) * .5;
			tol1 = *eps * fabs(x) + *tol / 3.;
			tol2 = tol1 * 2.;
			++(*iter);
			
			
			/* Check the (somewhat peculiar) stopping criterion: note that
			 the RHS is negative as long as the interval [a,b] is not small:*/
			if (fabs(x - xm) <= tol2 - (b - a) * .5 || *iter > maxit)
				goto L_End;
			
			
			/* is golden-section necessary */
			
			if (fabs(e) <= tol1 ||
				/*  if had Inf then go to golden-section */
				fx >= BIG_f || fv >= BIG_f || fw >= BIG_f) goto L_GoldenSect;
			
			/* Fit Parabola */
			
			r = (x - w) * (fx - fv);
			q = (x - v) * (fx - fw);
			p = (x - v) * q - (x - w) * r;
			q = (q - r) * 2.;
			if (q > 0.)
				p = -p;
			q = fabs(q);
			r = e;
			e = d;
			
			/* is parabola acceptable?  Otherwise do golden-section */
			
			if (fabs(p) >= fabs(.5 * q * r) ||
				q == 0.)
			/* above line added by BDR;
			 * [the abs(.) >= abs() = 0 should have branched..]
			 * in FTN: COMMON above ensures q is NOT a register variable */
				
				goto L_GoldenSect;
			
			if (p <= q * (a - x) ||
				p >= q * (b - x))			goto L_GoldenSect;
			
			
			
			/* Parabolic Interpolation step */
			d = p / q;
			u = x + d;
			
			/* f must not be evaluated too close to ax or bx */
			if ((u - a < tol2 || b - u < tol2)) {
				d = abs(tol1) * sgn(xm - x);
			}	
			
			goto L50;
			/*------*/
			
		L_GoldenSect: /* a golden-section step */
			
			if (x >= xm)    e = a - x;
			else/* x < xm*/ e = b - x;
			d = c_Gold * e;
			
			
		L50:
			u = x + ((fabs(d) >= tol1) ? d : (abs(tol1)*sgn(d)));
			/*  tol1 check : f must not be evaluated too close to x */
			
			*spar = u;
			*lspar = ratio * pow(16.0, (*spar * 6.0 - 2.0));
			sslvrg_(penalt, dofoff, xs, ys, ws, ssw, n,
							 knot, nk,
							 coef, sz, lev, crit, icrit, lspar, xwy,
					hs0, hs1, hs2, hs3,
					sg0, sg1, sg2, sg3, abd,
					p1ip, p2ip, ld4, ldnk, ier);
			fu = *crit;
			
			if(isnan(fu)) {
				fu = 2. * BIG_f;
			}
			
			/*  update  a, b, v, w, and x */
			
			if (fu <= fx) {
				if (u >= x) a = x; else b = x;
				
				v = w; fv = fw;
				w = x; fw = fx;
				x = u; fx = fu;
			}
			else {
				if (u < x)  a = u; else b = u;
				
				if (fu <= fw || w == x) {		        /* L70: */
					v = w; fv = fw;
					w = u; fw = fu;
				} else if (fu <= fv || v == x || v == w) {	/* L80: */
					v = u; fv = fu;
				}
			}
		}/* end main loop -- goto L20; */
		
	L_End:
		
		*spar = x;
		*crit = fx;
				
		//free memory
		delete [] xwy;
		delete [] hs0;
		delete [] hs1;
		delete [] hs2;
		delete [] hs3;
		delete [] sg0;
		delete [] sg1;
		delete [] sg2;
		delete [] sg3;
		delete [] abd;
		delete [] p1ip;
		delete [] p2ip;
		
		return 0;
		/* sbart */

	}catch(exception& e) {
		m->errorOut(e, "Spline", "sbart");
		exit(1);
	}
}	

/***********************************************************/