/* INTERP_SHELL.C: shell/driver program for interpolation routines,
 *                 interpolate.c, interpolate_wts.c 
 *                 interp.c, interp_wts.c

                   John Rust, University of Maryland, June, 2005 */

#include "mex.h"
#include <math.h>
#include "numrec.h"

double *vf,*p,*iv,*sv,*grid,*gridnumd;
int *gridnum;
int n,d,intmeth;

#define minc(x,y) ((x <= y)? x:y)
#define maxc(x,y) ((x >= y)? x:y)
#define sign(a) ((a>0.0)? 1.0:-1.0)

#include "mysort.c"   /* sorts a list of numbers */
#include "binsucc.c"  /* binary successor operation */
#include "power.c"    /* computes binary exponential, 2^d */
#include "kseq.c"     /* converts an array of indices into an integer
			 giving the position number of the array element
			 specified by the vector if indices as a single
			 integer in the vectorized array (i.e. a linear
			 sequence of numbers) */
#include "vertwt.c"   /* computes the probability weights for each vertex
			 of the grid hypercube containing a point to be
			 interpolated by interpolate.c or interp.c such that
			 the expectation of these weights with the function
			 to be interpolated (evaluated at the grid points)
			 equals the interpolated value. This is used to represent
			 the interpolation operator as an expectation operator
			 with respect to the function to be interpolated at the
			 grid points */
#include "interp.c"   
     /* routine that actually performs the multidimensional
        interpolation (calls kseq, binsucc and mysort) */
#include "interp_wts.c"  
     /* routine that actually performs the multidimensional
        interpolation (calls kseq, binsucc, vertwt and mysort).
        Differs from interp.c by computing a vector p that 
        represents the interpolation operation, p'vf=v(s) */
#include "interp_extrap.c" 
     /* routine that actually performs the multidimensional
        interpolation (calls kseq, binsucc and mysort) 
        This is same as interp.c except it also extrapolates
        and computes interpolated gradient */
#include "interp_extrap_wts.c" 
     /* routine that actually performs the multidimensional
        interpolation (calls kseq, binsucc, vertwt, and mysort) 
        This is same as interp_extrap.c except it also extrapolates
        and computes interpolated gradient and a vector p that
        represents the interpolation, i.e. p'vf=v(s), when all
        components of s are within the upper and lower grid points */

void interp_shell(void)
{
   int i,j;
   double *grad;

   grad=(double *) malloc(d*sizeof(double));

   for (i=0; i<n; i++)
   {
/* iv[i]=interp_wts(sv+d*i,p,vf,grid,gridnum,d,intmeth); */
/*  iv[i]=interp(sv+d*i,vf,grid,gridnum,d,intmeth); */
/* iv[i]=interp_extrap(sv+d*i,grad,vf,grid,gridnum,d,intmeth); */
 iv[i]=interp_extrap_wts(sv+d*i,grad,p,vf,grid,gridnum,d,intmeth); 
   }

   free(grad);
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
     int i;

     if (nrhs != 9)
     {
        mexErrMsgTxt("Error in interp_shell.c: 9 inputs required.");
     }

     n=(int) mxGetScalar(prhs[0]);
     d=(int) mxGetScalar(prhs[1]);
     grid=mxGetPr(prhs[2]);
     gridnumd=mxGetPr(prhs[3]);
     sv=mxGetPr(prhs[4]);
     vf=mxGetPr(prhs[5]);
     p=mxGetPr(prhs[6]);
     iv=mxGetPr(prhs[7]);
     intmeth=(int) mxGetScalar(prhs[8]);

     gridnum=(int *) malloc(d*sizeof(int));

     for (i=0; i<d; i++)
     {
        gridnum[i]=(int) gridnumd[i];	
     }

     interp_shell();

     free(gridnum);

}