/*

   gcvspl.C	(spline smooth)

   A Matlab mex interface for GCVSPL.
   Use 'cmex splsmoth.c gcvspl.f' to compile it   
 
   Returns the values at XX of the cubic smoothing spline for the
    given data (X,Y) by calling Waltring's GCVSPL (general cross vailidation
    spline). 
 
  The calling syntax is:
  [YY [, PP [,ERR ]]]=gcvspl(Y)
    equal interval,  assume X=[0:lenght(Y)-1] with default smooth parameters
  [YY [,PP, [ERR]]]=gcvspl(Y, smooth_para)
    equal interval,  assume X=[0:lenght(Y)]
  [YY[, PP[, ERR]]]=gcvspl (Y, X, smooth_para[])
    smooth Y at X,  return YY, PP as the same X
  [YY[, PP[, ERR]]]=gcvspl (Y, X, XX, smooth_para[])
    smooth Y at X,  return YY, PP at XX (resampling)

   smooth_para[0]=spline_mode  ;M in gcvspl default: 3, quintic
   smooth_para[1]=opt_mode     ; MD in gcvspl   default: 2, cross validation
   smooth_para[2]=VAL          ; mod value	default: 1,

   ERR=[error, WK(0:5)] of gcvspl
                       WK(1) = Generalized Cross Validation value
                       WK(2) = Mean Squared Residual.
                       WK(3) = Estimate of the number of degrees of
                               freedom of the residual sum of squares
                               per dataset, with 0.lt.WK(3).lt.N-M.
                       WK(4) = Smoothing parameter p, multiplicative
                               with the splines' derivative constraint.
                       WK(5) = Estimate of the true mean squared error
                               (different formula for |MD| = 3).
                       WK(6) = Gauss-Markov error variance.

                       If WK(4) -->  0 , WK(3) -->  0 , and an inter-
                       polating spline is fitted to the data (p --> 0).
                       A very small value > 0 is used for p, in order
                       to avoid division by zero in the GCV function.

                       If WK(4) --> inf, WK(3) --> N-M, and a least-
                       squares polynomial of order M (degree M-1) is
                       fitted to the data (p --> inf). For numerical
                       reasons, a very high value is used for p.

                       Upon return, the contents of WK can be used for
                       covariance propagation in terms of the matrices
                       B and WE: see the source listings. The variance
                       estimate for dataset J follows as WK(6)/WY(J).

   
   opt_mode     ( I )   Optimization mode switch:
	       | opt_mode| = 1: Prior given value for p in VAL
			 (VAL.ge.ZERO). This is the fastest
			 use of GCVSPL, since no iteration
			is performed in p.
	       | opt_mode| = 2: Generalized cross validation.
	       | opt_mode| = 3: True predicted mean-squared error,
			 with prior given variance in VAL.
	       | opt_mode| = 4: Prior given number of degrees of
			 freedom in VAL (ZERO.le.VAL.le.N-M).
			 X, W, M, N, and WK have not been
			 modified since the previous invoca-
			 tion of GCVSPL. 
			 If MD < -1, WK(4)
			 is used as an initial estimate for
			 the smoothing parameter p.  At the
			 first call to GCVSPL, MD must be > 0.
	       Other values for |MD|, and inappropriate values
	       for VAL will result in an error condition, or
	       cause a default value for VAL to be selected.
	       After return from MD.ne.1, the same number of
	       degrees of freedom can be obtained, for identical
	       weight factors and knot positions, by selecting
	       |MD|=1, and by copying the value of p from WK(4)
	       into VAL. In this way, no iterative optimization
	       is required when processing other data in Y. 
	       
   by C.X. Tian,  Last updated: Wed Jun 12 09:02:43 MST 1996
   
   tian@motorcortex.eas.asu.edu,  
   Computational Motor Control Lab, Arizona State University
	
  On PC platforms, must be compiled with /Alfw and /Gs option to generate
  proper code.

*/
#include <math.h>
#include "mex.h"

/* Input Arguments */

#define	IN1	prhs[0]
#define	IN2	prhs[1]
#define	IN3	prhs[2]
#define	IN4	prhs[3]

/* Output Arguments */

#define	OUT1	plhs[0]
#define	OUT2	plhs[1]
#define	OUT3	plhs[2]


extern void gcvspl_();
extern double splder_();

void mexFunction(
	int		nlhs,
	Matrix	*plhs[],
	int		nrhs,
	Matrix	*prhs[]
	)
{
	double	*yy, *pp, *err;
	double	*x, *y, *xx, *smooth_para;
	double  default_smooth_para[]= {2.0, /*cubic spline*/
		         	        2.0, /*general cross validation*/
			                1.0 };
			           
/* parameters for gcvspl_ and spl*/
        int	i, j, k, l;
	int     num_raw_data, num_output_data,  num_data_sets, deriv, iwork, index;
	int	spline_mode, opt_mode,error;
	double	 opt_val;
	double	time_instant;
	double  *weight_x, *weight_y, *work, *fortran_data;

	/* set different initial conditions accroding to nrhs */
	
	if (nrhs>0){
	    y=mxGetPr(IN1);
	    num_raw_data=num_output_data=mxGetM(IN1);
	    num_data_sets=mxGetN(IN1);
	} else 
	    mexErrMsgTxt("Too few input parameters");
	    
	switch (nrhs)  {
	    case 1:
		x=(double *)mxMalloc(num_raw_data, sizeof(double));
		smooth_para=default_smooth_para;
		for (i=0; i<num_raw_data; i++) x[i]=(double) i;
		xx=x;   /*output interval same as input*/
		break;
	    case 2:
		smooth_para=mxGetPr(IN2);
		x=(double *)mxMalloc(num_raw_data, sizeof(double));
		xx=x;     /*output interval same as input*/
		for (i=0; i<num_raw_data; i++) x[i]=(double) i;
		break;
	    case 3:
		x=mxGetPr(IN2);
		smooth_para=mxGetPr(IN3);
		xx=x;
		break;
	    case 4:
		x=mxGetPr(IN2);
		xx=mxGetPr(IN3);
		num_output_data=mxGetM(IN3);    /*input xx should be a col vector*/
		smooth_para=mxGetPr(IN4);
		break;
	    default:
		mexErrMsgTxt("calling syntax invalid.");
		break;
	} 

  	/* Create matries for the return argument */
	OUT1 = mxCreateFull(num_output_data, num_data_sets, REAL); /*yy*/
	OUT2 = mxCreateFull(num_output_data, num_data_sets, REAL); /*pp*/
	OUT3 = mxCreateFull(7, 1, REAL); /*ERR*/ 

	/* Assign pointers to the various parameters */
	yy = mxGetPr(OUT1);
	pp = mxGetPr(OUT2);
	err = mxGetPr(OUT3);
	
	/* Do the actual computations in a subroutine */
	/*prepare input for gcvspl*/
	spline_mode=(int)smooth_para[0];
	opt_mode=(int)smooth_para[1];
	opt_val=smooth_para[2];
        deriv = 0;   
	/*allocate memory*/
	weight_x=(double *)mxCalloc(num_raw_data, sizeof(double));
	weight_y=(double *)mxCalloc(num_raw_data, sizeof(double));
	work=(double *)mxCalloc(6*(num_raw_data*spline_mode+1)+num_raw_data, sizeof(double));
	
	for(i=0; i<num_raw_data; i++)   {
	    weight_x[i] = 1.0;
	    weight_y[i] = 1.0;
	}   

   	gcvspl_(x,
	    y,
	    &num_raw_data,
	    weight_x,
	    weight_y,
	    &spline_mode,
	    &num_raw_data,
	    &num_data_sets,
	    &opt_mode,
	    &opt_val,
	    pp,
	    &num_raw_data,
	    work,
	    &error );
	
	err[0]=error;
	for (i=1;i<7;i++) err[i]=work[i]; 
		
/*	mexPrintf("W[1]= %f  W[2]= %f  W[3]=%f  W[4]=%f\n", work[1], work[2], work[3], 
		work[4]);  
*/	/*
	** Use the program splder() to reconstruct various representations of marker
	** trajectories given the b-spline coefficients calculated using the gcvspl
	*/
	
    
	for (j = 0; j < num_data_sets; j++) {
	    iwork = 0; 
	    for (l = 0; l < num_output_data; l++){
		time_instant=xx[l];
		*yy++ =splder_(&deriv,
		    &spline_mode,
		    &num_raw_data,
		    &time_instant,
		    xx,
		    &pp[num_raw_data*j],
		    &iwork,
		    work);
	    }
	}
	
	/*clean up*/
	return;
}