mls.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "linalg.h"
#include "combinadic.h"
#include "trees.h"
#include "math_config.h"

typedef struct
{
	/**
	 * The number of points used in the interpolation;
	 * this is also the number of interpolating basis
	 * functions used in the reconstruction
	 */
	int np;

	/**
	 * The dimension of the points given as input
	 */
	int dim;

	/**
	 * The degree of interpolating polynomial to seek as
	 * the weighted reconstruction
	 */
	int deg;

	/**
	 * The list of points used in the interpolation
	 */
	double *pts;

	/**
	 * The dilation parameters for each point in pts
	 */
	double *dlts;

	/**
	 * The list of function values at their corresponding
	 * point in pts; these can be changed arbitrarily during
	 * the process of evaluating the interpolant; this is a
	 * single real value corresponding to each dim-dimensional
	 * point in pts
	 */
	double *vals;

	/**
	 * Function pointer to the weight function; this function
	 * is shifted to properly center in about each interpolating
	 * point stored in pts; the first argument is the dimension;
	 * the second is the actual point, which must be centered on
	 * the appropriate point by the user, at which to evaluate
	 */
	double (*wfs)(int,double,double*); /* Add another double to the arg list for setting the local radius of the basis function */
} mls_t;

typedef struct
{
	/**
	 * The number of interpolation points, i.e. number of
	 * basis functions used to approximate a function
	 */
	int np;

	/**
	 * The dimension of the domain
	 */
	int dim;

	/**
	 * The degree of the local basis to generate
	 */
	int deg;

	/**
	 * Pointer to the list of vertices
	 */
	double *pts;

	/**
	 * List of indexes into pts corresponding to boundary
	 * particles; not yet used, but implement somehow soon
	 */
	int *bpts;

	/**
	 * The dilation parameters for each point in pts
	 */
	double *dlts;

	/**
	 * Values of the interpolant at the desired function
	 * points stored in pts
	 */
	double *vals;

	/**
	 * Global quadrature weights used for evaluating the
	 * interpolant; there are as many values in gqw as
	 * there are points in pts
	 */
	double *gqw;

	/**
	 * Pointer to the weight functions prototype to use
	 */
	double (*wfs)(int,double,double*);

	/**
	 * Pointer to function which returns the gradient of
	 * the weight functions; just for testing purposes
	 */
	void (*wfsg)(int,double,double*,double*);

	/**
	 * Radius of the function wfs when a_in is unity
	 */
	double wrad;

	/**
	 * The Gram matrix defining the reproducing kernel of
	 * the given desired order using the window function
	 * given as wfs_in
	 */
	double *grm;

	/**
	 * Use this space to stored the LU-decomposed grm matrix
	 * internally; using this space so as not to disturb the
	 * original undilated Gram matrix
	 */
	double *grmlu;

	/**
	 * Store the coefficients here for the correction function;
	 * this is currently used for only one correction function,
	 * but must store np * pdim coefficients; assumes all particle
	 * basis functions have the same polynomial order
	 */
	double *coeffs; /* Adjust this to be used as the storage for the INTERNAL basis polynomial coefficients */

	/**
	 * Wavelet coefficients generated by wavelet generator
	 */
	double *wcoeffs; /* Do something similar here */

	/**
	 * A multi-use mode specifier
	 */
	int mode;
} rkp_t;

/**
 * Set up the data set and parameters for use in building the interpolant;
 * right now this just copies the input parameters into the mls_t struct;
 * maybe we can do without the mls_t struct all together and just use arrays
 * directly
 */
int mls_init( mls_t *obj_in, int np_in, int dim_in, int deg_in, double *pts_in, double *dlts_in, double *vals_in, double (*wfs_in)(int,double,double*) )
{
	obj_in->np = np_in;
	obj_in->dim = dim_in;
	obj_in->deg = deg_in;
	obj_in->pts = pts_in; /* Passing by reference here; do not free pts_in until finished */
	obj_in->dlts = dlts_in; /* Passing the dilation parameters */
	obj_in->vals = vals_in; /* Same here */
	obj_in->wfs = wfs_in;

	return 0;
}

/**
 * This function initializes the system by doing the integrals using Monte Carlo
 * averaging; this is not the best way to do this, but it will work as 
 * an optional method; this is why this function ends in mc, indicating the
 * initialization using Monte Carlo instead of a systematic quadrature
 * @param obj_in The reproducing kernel particle object
 * @param np_in The number of points about which basis particles are placed
 * @param dim_in The dimension of the vertices
 * @param deg_in The order of the polynomial space to use to generate the moment matrix
 * @param pts_in The set of interpolation points at which to place the RKPs
 * @param vals_in The values of the interpolant at each entry in pts_in
 * @param wfs_in The pointer to the window function to use for the system
 * @param rad_in This is the radius of the window function defining the support of the local function
 * @return Returns 0 if all is well or -1 if some error
 */
int rkp_init( rkp_t *obj_in, int np_in, int dim_in, int deg_in, double *pts_in, double *dlts_in, double *gqw_in, double *vals_in, double (*wfs_in)(int,double,double*), double wrad_in )
{
	int i;

	obj_in->np = np_in;
	obj_in->dim = dim_in;
	obj_in->deg = deg_in;
	obj_in->pts = (double*) malloc( dim_in * np_in * sizeof(double) );
	if( obj_in->pts == NULL )
		return -1;
	for(i=0;i<dim_in*np_in;i++)
		obj_in->pts[i] = pts_in[i];
	obj_in->dlts = (double*) malloc( np_in * sizeof(double) );
        if( obj_in->dlts == NULL )
                return -1;
        for(i=0;i<np_in;i++)
                obj_in->dlts[i] = dlts_in[i];
	obj_in->vals = (double*) malloc( np_in * sizeof(double) );
        if( obj_in->vals == NULL )
                return -1;
        for(i=0;i<np_in;i++)
                obj_in->vals[i] = vals_in[i];
	obj_in->gqw = (double*) malloc( np_in * sizeof(double) );
	if( obj_in->gqw == NULL )
		return -1;
	for(i=0;i<np_in;i++)
		obj_in->gqw[i] = gqw_in[i];
	obj_in->wfs = wfs_in;
	obj_in->wrad = wrad_in;
	obj_in->mode = 0;

	int n = binomial( dim_in + deg_in, dim_in );
	obj_in->grm = (double*) malloc( n * n * sizeof(double) ); /* Stores a matrix */
	obj_in->grmlu = (double*) malloc( n * n * sizeof(double) ); /* Stores a matrix */
	obj_in->coeffs = (double*) malloc( np_in * n * sizeof(double) ); /* Stores a vector */
	obj_in->wcoeffs = (double*) malloc( np_in * n * ( n - 1 ) * sizeof(double) ); /* Stores higher order wavelet coefficients */

	if( obj_in->grm == NULL || obj_in->grmlu == NULL || obj_in->coeffs == NULL || obj_in->wcoeffs == NULL )
		return -1;

	return 0;
}

int mls_free( mls_t *obj_in )
{
	return 0;
}

int rkp_free( rkp_t *obj_in )
{
	if( obj_in->vals != NULL )
		free( obj_in->vals );
	if( obj_in->pts != NULL )
		free( obj_in->pts );
	if( obj_in->dlts != NULL )
		free( obj_in->dlts );
	if( obj_in->gqw != NULL )
		free( obj_in->gqw );
	if( obj_in->grm != NULL )
		free( obj_in->grm );
	if( obj_in->grmlu != NULL )
		free( obj_in->grmlu );
	if( obj_in->coeffs != NULL )
		free( obj_in->coeffs );
	if( obj_in->wcoeffs != NULL )
		free( obj_in->wcoeffs );

	return 0;
}

/**
 * This function evaluate the polynomial basis vector at the specified
 * point given in x_in and outputs the elements to vec_out, which is
 * presumed to have enough space already allocated.
 * @param obj_in The MLS object
 * @param x_in The point at which to evaluate the basis; should have same dimension as mls_t
 * @param vec_out The space to which the p-vector is output; presumed allocated; will be assumed to have C(deg+dim,dim) entries
 * @return Returns 0 if all is well otherwise < 0
 */
int mls_basis_evaluate( mls_t *obj_in, double *x_in, double *vec_out ) /* TODO: Pre-allocate cmb and occ so we don't have to do it repeatedly */
{
	int i,j,k,m,n,*cmb,*occ;
	int pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim ); /* Make sure vec_out points to at least this many doubles */
	double prd;

	/* Allocate space for iterating combinations */
	cmb = (int*) malloc( obj_in->dim * sizeof(int) ); /* Only really need obj_in->dim - 1 but padding won't hurt */
	occ = (int*) malloc( obj_in->dim * sizeof(int) );
	for(i=0,m=0;i<=obj_in->deg;i++)
	{
		n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
		rcombinadic_init( i, obj_in->dim, cmb );
		for(j=0;j<n;j++)
		{
			rcombinadic_occupancy( i, obj_in->dim, cmb, occ );
			prd = 1.0;
			for(k=0;k<obj_in->dim;k++)
				prd *= pow( x_in[k], (double) occ[k] );
			vec_out[m++] = prd; /* Place each entry following those from its previous degrees continuously */
			rcombinadic_next( i, obj_in->dim, cmb );
		}
	}
	free( cmb );
	free( occ );

	return 0;
}

/**
 * This function evaluate the important A(x) matrix which must be inverted
 * in order to generate the global interpolant via the nodal basis functions
 * @param obj_in The MLS object
 * @param x_in The point at which to evaluate the matrix
 * @param mat_out The matrix as output; make sure this is allocated with enough space
 * @return Returns 0 if all went well
 */
int mls_matrix_evaluate( mls_t *obj_in, double dlt_in, double *x_in, double *mat_out )
{
	int i,j,k,n,p,pdim;
	double *x,*vec;

	/* Allocate the right size of vector for output from mls_basis_evaluate */
	pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );
	vec = (double*) malloc( pdim * sizeof(double) );

	/* Allocate a temporary dim-vector */
	x = (double*) malloc( obj_in->dim * sizeof(double) );

	/* Zero the matrix before forming the sum of outer products */
	for(i=0;i<pdim*pdim;i++)
		mat_out[i] = 0.0;

	/* Sum all outer products pj pj^T */
	for(i=0;i<obj_in->np;i++)
	{
		/* Evaluate the basis vector at each value of i */
		mls_basis_evaluate( obj_in, obj_in->pts + i * obj_in->dim, vec );

		/* Function obj_in->wfs at point i is obtained by translating it to be centered at obj_in->pts[i*obj_in->dim+...] */
		for(j=0;j<obj_in->dim;j++)
			x[j] = x_in[j] - obj_in->pts[i*obj_in->dim+j];

		/* Fill in only the upper triangular first; then copy to the lower */
		for(j=0;j<pdim;j++)
			for(k=j;k<pdim;k++)
				mat_out[j*pdim+k] += vec[j] * vec[k] * obj_in->wfs( obj_in->dim, dlt_in, x );
	}

	/* Copy the upper triangular part of mat_out into the lower triangular */
	for(j=0;j<pdim;j++)
		for(k=0;k<j;k++)
			mat_out[j*pdim+k] = mat_out[k*pdim+j];

	free( x );
	free( vec );

	return 0;
}

/**
 * This is where the major difference between RKP and MLS becomes
 * evident; this function will use Monte Carlo integration if the
 * option is chosen during initialization, or it will use a basic
 * stencil method based on polynomials if it is so initialized;
 * this matrix is also called the moment matrix in literature
 * @param obj_in The RKP object
 * @param x_in The point at which to evaluate the matrix
 * @return Returns 0 if all went well
 */
int rkp_matrix_evaluate( rkp_t *obj_in, double *x_in )
{
	int i,j,k,m,n,p,nq,ret,pdim,*cmb,*occ,*ord;
	double prd,*x,*vec;

	/* Set up iteration through monomials */
	cmb = (int*) malloc( obj_in->dim * sizeof(int) ); /* Only really need obj_in->dim - 1 but padding won't hurt */
        occ = (int*) malloc( obj_in->dim * sizeof(int) );

	/* Allocate the right size of vector for output from mls_basis_evaluate */
        pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

	/* Allocate a temporary dim-vector */
        x = (double*) malloc( obj_in->dim * sizeof(double) );
	vec = (double*) malloc( pdim * sizeof(double) );

	/* Zero the matrix before forming the sum of outer products */
        for(i=0;i<pdim*pdim;i++)
                obj_in->grm[i] = 0.0;

	/* Just iterate through all nodes */
	for(i=0;i<obj_in->np;i++)
	{
		/* Generate the point */
		for(j=0;j<obj_in->dim;j++)
			x[j] = obj_in->pts[i*obj_in->dim+j] - x_in[j];

		/* Generate the polynomial vector */
		for(p=0,j=0;j<=obj_in->deg;j++)
		{
			rcombinadic_init( j, obj_in->dim, cmb );
			n = binomial( j + obj_in->dim - 1, obj_in->dim - 1 );
			for(k=0;k<n;k++)
			{
				rcombinadic_occupancy( j, obj_in->dim, cmb, occ );
				prd = 1.0;
				for(m=0;m<obj_in->dim;m++)
					prd *= pow( x[m], (double) occ[m] );
				vec[p++] = prd;
				rcombinadic_next( j, obj_in->dim, cmb );
			}
		}

		/* Now form the outer product */
		for(j=0;j<pdim;j++)
		{
			for(k=j;k<pdim;k++)
			{
				prd = vec[j] * vec[k] * obj_in->wfs( obj_in->dim, obj_in->dlts[i], x );
				obj_in->grm[j*pdim+k] += obj_in->gqw[i] * prd;
			}
		}
	}

	/* Copy into lower triangular area */
	for(i=0;i<pdim;i++)
		for(j=0;j<i;j++)
			obj_in->grm[i*pdim+j] = obj_in->grm[j*pdim+i];

	free( cmb );
	free( occ );
	free( x );
	free( vec );

	return 0;
}

/**
 * Does the same thing as rkp_matrix_evaluate with the
 * exception that the x_in values are scaled by the dilation
 * parameters of each local particle
 * @param obj_in RKP object
 * @param x_in Point at which to evaluate moments
 */
int rkp_matrix_evaluate_scaled( rkp_t *obj_in, double *x_in )
{
    int i,j,k,m,n,p,nq,ret,pdim,*cmb,*occ,*ord;
    double prd,*x,*vec;

    /* Set up iteration through monomials */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) ); /* Only really need obj_in->dim - 1 but padding won't hurt */
    occ = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Allocate the right size of vector for output from mls_basis_evaluate */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate a temporary dim-vector */
    x = (double*) malloc( obj_in->dim * sizeof(double) );
    vec = (double*) malloc( pdim * sizeof(double) );

    /* Zero the matrix before forming the sum of outer products */
    for(i=0;i<pdim*pdim;i++)
        obj_in->grm[i] = 0.0;

    /* Just iterate through all nodes */
    for(i=0;i<obj_in->np;i++)
    {
        /* Generate the point */
        for(j=0;j<obj_in->dim;j++)
            x[j] = obj_in->pts[i*obj_in->dim+j] - x_in[j];

        /* Generate the polynomial vector */
        for(p=0,j=0;j<=obj_in->deg;j++)
        {
            rcombinadic_init( j, obj_in->dim, cmb );
            n = binomial( j + obj_in->dim - 1, obj_in->dim - 1 );
            for(k=0;k<n;k++)
            {
                rcombinadic_occupancy( j, obj_in->dim, cmb, occ );
                prd = 1.0;
                for(m=0;m<obj_in->dim;m++)
                    prd *= pow( x[m] / obj_in->dlts[i], (double) occ[m] );
                vec[p++] = prd;
                rcombinadic_next( j, obj_in->dim, cmb );
            }
        }

        /* Now form the outer product */
        for(j=0;j<pdim;j++)
        {
            for(k=j;k<pdim;k++)
            {
                prd = vec[j] * vec[k] * obj_in->wfs( obj_in->dim, obj_in->dlts[i], x );
                obj_in->grm[j*pdim+k] += obj_in->gqw[i] * prd;
            }
        }
    }

    /* Copy into lower triangular area */
    for(i=0;i<pdim;i++)
        for(j=0;j<i;j++)
            obj_in->grm[i*pdim+j] = obj_in->grm[j*pdim+i];

    free( cmb );
    free( occ );
    free( x );
    free( vec );

    return 0;
}

/**
 * Use this function to deal with widely varying particle sizes
 * when generating wavelets for use in evaluating derivatives
 * @param obj_in RKP object to use
 * @param x_in Point at which to evaluate the moments
 * @param rho_in Scale on which to evaluate the moments
 * @return Returns 0 if all is well
 */
int rkp_matrix_evaluate_scaled_const( rkp_t *obj_in, double *x_in, double rho_in )
{
    int i,j,k,m,n,p,nq,ret,pdim,*cmb,*occ,*ord;
    double prd,*x,*vec;

    /* Set up iteration through monomials */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) ); /* Only really need obj_in->dim - 1 but padding won't hurt */
    occ = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Allocate the right size of vector for output from mls_basis_evaluate */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate a temporary dim-vector */
    x = (double*) malloc( obj_in->dim * sizeof(double) );
    vec = (double*) malloc( pdim * sizeof(double) );

    /* Zero the matrix before forming the sum of outer products */
    for(i=0;i<pdim*pdim;i++)
        obj_in->grm[i] = 0.0;

    /* Just iterate through all nodes */
    for(i=0;i<obj_in->np;i++)
    {
        /* Generate the point */
        for(j=0;j<obj_in->dim;j++)
            x[j] = obj_in->pts[i*obj_in->dim+j] - x_in[j];

        /* Generate the polynomial vector */
        for(p=0,j=0;j<=obj_in->deg;j++)
        {
            rcombinadic_init( j, obj_in->dim, cmb );
            n = binomial( j + obj_in->dim - 1, obj_in->dim - 1 );
            for(k=0;k<n;k++)
            {
                rcombinadic_occupancy( j, obj_in->dim, cmb, occ );
                prd = 1.0;
                for(m=0;m<obj_in->dim;m++)
                    prd *= pow( x[m] / rho_in, (double) occ[m] );
                vec[p++] = prd;
                rcombinadic_next( j, obj_in->dim, cmb );
            }
        }

        /* Now form the outer product */
        for(j=0;j<pdim;j++)
        {
            for(k=j;k<pdim;k++)
            {
                /* Need to scale the input to the window function here */
                prd = vec[j] * vec[k] * obj_in->wfs( obj_in->dim, obj_in->dlts[i], x );
                obj_in->grm[j*pdim+k] += obj_in->gqw[i] * prd;
            }
        }
    }

    /* Copy into lower triangular area */
    for(i=0;i<pdim;i++)
        for(j=0;j<i;j++)
            obj_in->grm[i*pdim+j] = obj_in->grm[j*pdim+i];

    free( cmb );
    free( occ );
    free( x );
    free( vec );

    return 0;
}

/**
 * This function generates the matrix at the given point and
 * inverts it to give the coefficients for the correction function
 * by calling rkp_matrix_evaluate to get the moments
 * @param obj_in The RKP object
 * @param x_in The point at which to generate
 * @param dlt_in The dilation factor to use
 * @param coeffs_out Where to put the generated coefficients for the correction function
 * @return Returns 0 if all is well, < 0 otherwise
 */
int rkp_matrix_generate( rkp_t *obj_in, double *x_in, double *coeffs_out )
{
    int i,j,ret,pdim,*ipiv;

    /* Build the matrix */
    rkp_matrix_evaluate( obj_in, x_in );

    /* Set the dimension and allocate pivot space */
    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    ipiv = (int*) malloc( pdim * sizeof(int) );

    /* Copy into grmlu */
    for(i=0;i<pdim*pdim;i++)
        obj_in->grmlu[i] = obj_in->grm[i];

    /* Invert the matrix and put the LU-decomposed matrix in grmlu */
    j = 1;
    coeffs_out[0] = 1.0;
    for(i=1;i<pdim;i++)
        coeffs_out[i] = 0.0;
    dgesv_( &pdim, &j, obj_in->grmlu, &pdim, ipiv, coeffs_out, &pdim, &ret );

    free( ipiv );

    return ret;
}

/**
 * This function generates the matrix at the given point and
 * inverts it to give the coefficients for the correction function
 * by calling rkp_matrix_evaluate to get the moments
 * @param obj_in The RKP object
 * @param x_in The point at which to generate
 * @param dlt_in The dilation factor to use
 * @param coeffs_out Where to put the generated coefficients for the correction function
 * @return Returns 0 if all is well, < 0 otherwise
 */
int rkp_matrix_generate_scaled( rkp_t *obj_in, double *x_in, double *coeffs_out )
{
    int i,j,ret,pdim,*ipiv;

    /* Build the matrix */
    rkp_matrix_evaluate_scaled( obj_in, x_in );

    /* Set the dimension and allocate pivot space */
    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    ipiv = (int*) malloc( pdim * sizeof(int) );

    /* Copy into grmlu */
    for(i=0;i<pdim*pdim;i++)
        obj_in->grmlu[i] = obj_in->grm[i];

    /* Invert the matrix and put the LU-decomposed matrix in grmlu */
    j = 1;
    coeffs_out[0] = 1.0;
    for(i=1;i<pdim;i++)
        coeffs_out[i] = 0.0;
    dgesv_( &pdim, &j, obj_in->grmlu, &pdim, ipiv, coeffs_out, &pdim, &ret );

    free( ipiv );

    return ret;
}

/**
 * This function generates the matrix at the given point and
 * inverts it to give the coefficients for the correction function
 * by calling rkp_matrix_evaluate_scaled_const to get the moments;
 * This function only generates the correction function for the partition
 * of unity; it does not generate the wavelets
 * @param obj_in The RKP object
 * @param x_in The point at which to generate
 * @param coeffs_out Where to put the generated coefficients for the correction function
 * @param rho_in The position-dependent dilation factor to use to generate the Gram matrix of moments
 * @return Returns 0 if all is well, < 0 otherwise
 */
int rkp_matrix_generate_scaled_const( rkp_t *obj_in, double *x_in, double *coeffs_out, double rho_in )
{
    int i,j,ret,pdim,*ipiv;

    /* Build the matrix */
    rkp_matrix_evaluate_scaled_const( obj_in, x_in, rho_in );

    /* Set the dimension and allocate pivot space */
    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    ipiv = (int*) malloc( pdim * sizeof(int) );

    /* Copy into grmlu */
    for(i=0;i<pdim*pdim;i++)
        obj_in->grmlu[i] = obj_in->grm[i];

    /* Invert the matrix and put the LU-decomposed matrix in grmlu */
    j = 1;
    coeffs_out[0] = 1.0;
    for(i=1;i<pdim;i++)
        coeffs_out[i] = 0.0;
    dgesv_( &pdim, &j, obj_in->grmlu, &pdim, ipiv, coeffs_out, &pdim, &ret );

    free( ipiv );

    return ret;
}

/**
 * Generate wavelet wdx_in about point x_in; it is important here to note
 * that if wdx_in is equal to 0 (zero), then this will just generate the
 * particle functions, in the same way as rkp_matrix_generate.
 */
int rkp_wavelet_generate_order( rkp_t *obj_in, int wdx_in, double *x_in, double *coeffs_out )
{
    int i,j,k,m,ret,pdim,*ipiv,*exp;

    /* Build the matrix centered at x_in */
    rkp_matrix_evaluate( obj_in, x_in );

    /* Set the dimensionality */
    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    ipiv = (int*) malloc( pdim * sizeof(int) );
    exp = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Copy into obj_in->grmlu for calculation */
    for(i=0;i<pdim*pdim;i++)
        obj_in->grmlu[i] = obj_in->grm[i];

    global_polynomial_vector( wdx_in, obj_in->dim, exp );

    /* Solve the right system */
    j = 1;
    for(i=0;i<pdim;i++)
        if( i == wdx_in )
        {
            coeffs_out[i] = 1.0;
            m = 1;
            for(k=0;k<obj_in->dim;k++)
                m *= factorial( exp[k] );
            coeffs_out[i] /= (double) m;
        }
        else
            coeffs_out[i] = 0.0;
    dgesv_( &pdim, &j, obj_in->grmlu, &pdim, ipiv, coeffs_out, &pdim, &ret );

    free( ipiv );
    free( exp );

    return 0;
}

/**
 * Generate wavelet wdx_in about point x_in; it is important here to note
 * that if wdx_in is equal to 0 (zero), then this will just generate the
 * particle functions, in the same way as rkp_matrix_generate. This is the
 * scaled version of the function.
 */
int rkp_wavelet_generate_order_scaled( rkp_t *obj_in, int wdx_in, double *x_in, double *coeffs_out )
{
    int i,j,k,m,ret,pdim,*ipiv,*exp;

    /* Build the matrix centered at x_in */
    rkp_matrix_evaluate_scaled( obj_in, x_in );

    /* Set the dimensionality */
    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    ipiv = (int*) malloc( pdim * sizeof(int) );
    exp = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Copy into obj_in->grmlu for calculation */
    for(i=0;i<pdim*pdim;i++)
        obj_in->grmlu[i] = obj_in->grm[i];

    global_polynomial_vector( wdx_in, obj_in->dim, exp );

    /* Solve the right system */
    j = 1;
    for(i=0;i<pdim;i++)
        if( i == wdx_in )
        {
            coeffs_out[i] = 1.0;
            m = 1;
            for(k=0;k<obj_in->dim;k++)
                m *= factorial( exp[k] );
            coeffs_out[i] /= (double) m;
        }
        else
            coeffs_out[i] = 0.0;
    dgesv_( &pdim, &j, obj_in->grmlu, &pdim, ipiv, coeffs_out, &pdim, &ret );

    free( ipiv );
    free( exp );

    return 0;
}

/**
 * Constantly scaled version of rkp_wavelet_generate_order
 */
int rkp_wavelet_generate_order_scaled_const( rkp_t *obj_in, int wdx_in, double *x_in, double rho_in, double *coeffs_out )
{
    int i,j,k,m,ret,pdim,*ipiv,*exp;

    /* Build the matrix centered at x_in */
    rkp_matrix_evaluate_scaled_const( obj_in, x_in, rho_in );

    /* Set the dimensionality */
    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    ipiv = (int*) malloc( pdim * sizeof(int) );
    exp = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Copy into obj_in->grmlu for calculation */
    for(i=0;i<pdim*pdim;i++)
        obj_in->grmlu[i] = obj_in->grm[i];

    global_polynomial_vector( wdx_in, obj_in->dim, exp );

    /* Solve the right system */
    j = 1;
    for(i=0;i<pdim;i++)
        if( i == wdx_in )
        {
            coeffs_out[i] = 1.0;
            m = 1;
            for(k=0;k<obj_in->dim;k++)
                m *= factorial( exp[k] );
            coeffs_out[i] /= (double) m;
        }
        else
            coeffs_out[i] = 0.0;
    dgesv_( &pdim, &j, obj_in->grmlu, &pdim, ipiv, coeffs_out, &pdim, &ret );

    free( ipiv );
    free( exp );

    return 0;
}

/**
 * Wavelet generator; solves the same problem as the above, but keeps
 * different columns of the resulting matrix inverse
 * @param obj_in The RKP object
 * @param x_in The point at which to center the matrix
 * @param dlt_in The dilation parameter to use about this center
 * @param coeffs_out The output storage of the coefficients generated
 * @return Returns 0 if all went well, < 0 otherwise
 */
int rkp_wavelet_generate( rkp_t *obj_in, double *x_in, double *coeffs_out )
{
	int i,j,k,ret,pdim,*ipiv;
	double *tmp;

	/* Build the matrix */
	rkp_matrix_evaluate( obj_in, x_in );

	/* Set the dimension and allocate pivot space */
	pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
	ipiv = (int*) malloc( pdim * sizeof(int) );

	/* Copy into grmlu */
	for(i=0;i<pdim*pdim;i++)
		obj_in->grmlu[i] = obj_in->grm[i];

	/* Allocate quick storage then release it */
	tmp = (double*) malloc( pdim * pdim * sizeof(double) );

	/* Save the system matrix in tmp temporarily */
	for(i=0;i<pdim*pdim;i++)
		tmp[i] = obj_in->grmlu[i];

	/* Solve the system and apply it to the correct righthand side pdim - 1 times */
	for(k=1;k<pdim;k++)
	{
		/* Copy the system into grmlu */
		for(i=0;i<pdim*pdim;i++)
			obj_in->grmlu[i] = tmp[i];

		/* Invert the matrix and put the LU-decomposed matrix in grmlu */
		j = 1;
		for(i=0;i<pdim;i++)
			if( i == k )
				coeffs_out[(k-1)*pdim+i] = 1.0;
			else
				coeffs_out[(k-1)*pdim+i] = 0.0;
		dgesv_( &pdim, &j, obj_in->grmlu, &pdim, ipiv, coeffs_out + ( k - 1 ) * pdim, &pdim, &ret );
	}

	/* Clean everything up */
	free( ipiv );
	free( tmp );

	return 0;
}

/**
 * Scaled version of rkp_wavelet_generate
 */
int rkp_wavelet_generate_scaled( rkp_t *obj_in, double *x_in, double *coeffs_out )
{
    int i,j,k,ret,pdim,*ipiv;
    double *tmp;

    /* Build the matrix */
    rkp_matrix_evaluate_scaled( obj_in, x_in );

    /* Set the dimension and allocate pivot space */
    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    ipiv = (int*) malloc( pdim * sizeof(int) );

    /* Copy into grmlu */
    for(i=0;i<pdim*pdim;i++)
        obj_in->grmlu[i] = obj_in->grm[i];

    /* Allocate quick storage then release it */
    tmp = (double*) malloc( pdim * pdim * sizeof(double) );

    /* Save the system matrix in tmp temporarily */
    for(i=0;i<pdim*pdim;i++)
        tmp[i] = obj_in->grmlu[i];

    /* Solve the system and apply it to the correct righthand side pdim - 1 times */
    for(k=1;k<pdim;k++)
    {
        /* Copy the system into grmlu */
        for(i=0;i<pdim*pdim;i++)
            obj_in->grmlu[i] = tmp[i];

        /* Invert the matrix and put the LU-decomposed matrix in grmlu */
        j = 1;
        for(i=0;i<pdim;i++)
            if( i == k )
                coeffs_out[(k-1)*pdim+i] = 1.0;
            else
                coeffs_out[(k-1)*pdim+i] = 0.0;
        dgesv_( &pdim, &j, obj_in->grmlu, &pdim, ipiv, coeffs_out + ( k - 1 ) * pdim, &pdim, &ret );
    }

    /* Clean everything up */
    free( ipiv );
    free( tmp );

    return 0;
}

/**
 * Constantly scaled version of rkp_wavelet_generate
 */
int rkp_wavelet_generate_scaled_const( rkp_t *obj_in, double *x_in, double rho_in, double *coeffs_out )
{
    int i,j,k,ret,pdim,*ipiv;
    double *tmp;

    /* Build the matrix */
    rkp_matrix_evaluate_scaled_const( obj_in, x_in, rho_in );

    /* Set the dimension and allocate pivot space */
    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    ipiv = (int*) malloc( pdim * sizeof(int) );

    /* Copy into grmlu */
    for(i=0;i<pdim*pdim;i++)
        obj_in->grmlu[i] = obj_in->grm[i];

    /* Allocate quick storage then release it */
    tmp = (double*) malloc( pdim * pdim * sizeof(double) );

    /* Save the system matrix in tmp temporarily */
    for(i=0;i<pdim*pdim;i++)
        tmp[i] = obj_in->grmlu[i];

    /* Solve the system and apply it to the correct righthand side pdim - 1 times */
    for(k=1;k<pdim;k++)
    {
        /* Copy the system into grmlu */
        for(i=0;i<pdim*pdim;i++)
            obj_in->grmlu[i] = tmp[i];

        /* Invert the matrix and put the LU-decomposed matrix in grmlu */
        j = 1;
        for(i=0;i<pdim;i++)
            if( i == k )
                coeffs_out[(k-1)*pdim+i] = 1.0;
            else
                coeffs_out[(k-1)*pdim+i] = 0.0;
        dgesv_( &pdim, &j, obj_in->grmlu, &pdim, ipiv, coeffs_out + ( k - 1 ) * pdim, &pdim, &ret );
    }

    /* Clean everything up */
    free( ipiv );
    free( tmp );

    return 0;
}

/**
 * This function fills in obj_in->coeffs for all of the
 * particle basis functions. Most of the near the middle
 * of the domain will be identical in their coefficients.
 * WARNING: This function needs serious work. It is pointless
 * to evaluate the correction polynomial with a fixed dilation
 * factor since dlt_in depends on the particle being evaluated
 * and the particle evaluated and the point at which the particle
 * is being evaluated are different, meaning that calling
 * rkp_matrix_generate with obj_in->dlts[i] and with
 * obj_in->pts[i*obj_in->dim] is both wrong and useless.
 * @param obj_in The RKP object to generate
 * @return Returns 0 if all went well
 */
int rkp_basis_generate( rkp_t *obj_in )
{
    int i,pdim;

    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    for(i=0;i<obj_in->np;i++)
        rkp_matrix_generate( obj_in, obj_in->pts + i * obj_in->dim, obj_in->coeffs + i * pdim );

    return 0;
}

/**
 * Scaled version of rkp_basis_generate
 */
int rkp_basis_generate_scaled( rkp_t *obj_in )
{
    int i,pdim;

    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    for(i=0;i<obj_in->np;i++)
        rkp_matrix_generate_scaled( obj_in, obj_in->pts + i * obj_in->dim, obj_in->coeffs + i * pdim );

    return 0;
}

/**
 * Constantly scaled version of rkp_basis_generate
 */
int rkp_basis_generate_scaled_const( rkp_t *obj_in )
{
    int i,pdim;

    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    for(i=0;i<obj_in->np;i++)
        rkp_matrix_generate_scaled_const( obj_in, obj_in->pts + i * obj_in->dim, obj_in->coeffs + i * pdim, obj_in->dlts[i] );

    return 0;
}

/**
 * This function does the same as rkp_basis_generate except it
 * generates the wavelet coefficients for all orders up to and
 * including order obj_in->deg-1.
 * WARNING: This function needs to be rethought for the same
 * reason shown above for rkp_basis_generate. Another means
 * of storing correction function coefficients needs to be
 * found in order to generalize this method
 * @param obj_in The RKP object to generate
 * @return Returns 0 if all went well
 */
int rkp_wavelet_basis_generate( rkp_t *obj_in )
{
    int i,pdim;

    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    for(i=0;i<obj_in->np;i++)
        rkp_wavelet_generate( obj_in, obj_in->pts + i * obj_in->dim, obj_in->wcoeffs + i * pdim * ( pdim - 1 ) );

    return 0;
}

/**
 * Scaled version of rkp_wavelet_basis_generate
 */
int rkp_wavelet_basis_generate_scaled( rkp_t *obj_in )
{
    int i,pdim;

    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    for(i=0;i<obj_in->np;i++)
        rkp_wavelet_generate_scaled( obj_in, obj_in->pts + i * obj_in->dim, obj_in->wcoeffs + i * pdim * ( pdim - 1 ) );

    return 0;
}

/**
 * Constantly scaled version of rkp_wavelet_basis_generate
 */
int rkp_wavelet_basis_generate_scaled_const( rkp_t *obj_in )
{
    int i,pdim;

    pdim = binomial( obj_in->dim + obj_in->deg, obj_in->dim );
    for(i=0;i<obj_in->np;i++)
        rkp_wavelet_generate_scaled_const( obj_in, obj_in->pts + i * obj_in->dim, obj_in->dlts[i], obj_in->wcoeffs + i * pdim * ( pdim - 1 ) );

    return 0;
}

/**
 * Perform steps to evaluate the MLS interpolant at the given point;
 * First evaluate A(x), then evaluate p(xj), W(x-xj)and finally p(x); then
 * invert A(x) and calculate p^T(x) * A^-1(x) p(xj) W(x-xj) = phij(x)
 * @param obj_in The MLS object
 * @param idx_in The index of the basis function to evaluate at x_in
 * @param x_in The point at which to evaluate the basis function
 * @return Returns the value of the basis function (term) at x_in
 */
double mls_term_evaluate( mls_t *obj_in, int idx_in, double *x_in )
{
    double val,*tmp,*x,*px,*pxi,*mat;
    int i,j,pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    tmp = (double*) malloc( pdim * sizeof(double) );
    x = (double*) malloc( obj_in->dim * sizeof(double) );
    px = (double*) malloc( pdim * sizeof(double) );
    pxi = (double*) malloc( pdim * sizeof(double) );
    mat = (double*) malloc( pdim * pdim * sizeof(double) );

    /* Calculate the x-dependent vectors and matrix */
    mls_basis_evaluate( obj_in, x_in, px );
    mls_matrix_evaluate( obj_in, obj_in->dlts[idx_in], x_in, mat );
    mls_basis_evaluate( obj_in, obj_in->pts + idx_in * obj_in->dim, pxi );

    /* Now invert A(x) */
    int ret;
    int nrhs = 1;
    int *ipiv = (int*) malloc( pdim * sizeof(int) );

    /* Invert the matrix and put it in apply it to p(x_idx_in) */
    dgesv_( &pdim, &nrhs, mat, &pdim, ipiv, pxi, &pdim, &ret );

    /* Now take the inner product with p(x) */
    val = 0.0;
    for(i=0;i<pdim;i++)
        val += px[i] * pxi[i];

    /* Localize the x variable around index idx_in in obj_in->pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = x_in[i] - obj_in->pts[idx_in*obj_in->dim+i];

    val = val * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x );

    free( x );
    free( tmp );
    free( px );
    free( pxi );
    free( mat );
    free( ipiv );

    return val;
}

/**
 * Evaluate a basis function centered on an interpolation point
 * at a given global position; this function assumes that the
 * matrix has already been generated and inverted and that the
 * resulting reproducing polynomial coefficients are stored in
 * obj_in->grmlu; always call rkp_matrix_generate before calling
 * this function if a matrix with a different dilation parameter
 * is needed; this function assumes that obj_in->coeffs is filled
 * in with the polynomial coefficients of the associated correction
 * functions
 * @param obj_in The RKP object
 * @param idx_in The index of the basis function to evaluate
 * @param x_in The point at which to evaluate the basis function
 * @return Returns the value of the function in double
 */
double rkp_term_evaluate( rkp_t *obj_in, int idx_in, double *x_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,*cf,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* If x_in is more than obj_in->dlts[idx_in] away from obj_in->pts[idx_in], then don't calculate this */
    sum = 0.0;
    for(i=0;i<obj_in->dim;i++)
        sum += pow( x_in[i] - obj_in->pts[idx_in*obj_in->dim+i], 2.0 );
    if( sum > obj_in->wrad * obj_in->dlts[idx_in] * obj_in->wrad * obj_in->dlts[idx_in] ) /* WARNING: Assumes that wfs has radius 2 * dlts[idx_in], not dlts[idx_in] */
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );
    cf = (double*) malloc( pdim * sizeof(double) );

    /* Generate the coefficients for the correction function */
    rkp_matrix_generate( obj_in, x_in, cf );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = obj_in->pts[idx_in*obj_in->dim+i] - x_in[i];

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= cf[p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );
    free( cf );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x ); /* This is 1.0 because already divided by dilation parameter!!! */
}

/**
 * Scaled version of rkp_term_evaluate
 */
double rkp_term_evaluate_scaled( rkp_t *obj_in, int idx_in, double *x_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,*cf,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* If x_in is more than obj_in->dlts[idx_in] away from obj_in->pts[idx_in], then don't calculate this */
    sum = 0.0;
    for(i=0;i<obj_in->dim;i++)
        sum += pow( x_in[i] - obj_in->pts[idx_in*obj_in->dim+i], 2.0 );
    if( sum > obj_in->wrad * obj_in->dlts[idx_in] * obj_in->wrad * obj_in->dlts[idx_in] ) /* WARNING: Assumes that wfs has radius 2 * dlts[idx_in], not dlts[idx_in] */
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );
    cf = (double*) malloc( pdim * sizeof(double) );

    /* Generate the coefficients for the correction function */
    rkp_matrix_generate_scaled( obj_in, x_in, cf );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = ( obj_in->pts[idx_in*obj_in->dim+i] - x_in[i] ) / obj_in->dlts[idx_in]; /* Important: Scale the input */

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= cf[p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );
    free( cf );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    for(i=0;i<obj_in->dim;i++)
        x[i] *= obj_in->dlts[idx_in]; /* Must undo dilation to give input to obj_in->wfs */
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x ); /* This is 1.0 because already divided by dilation parameter!!! */
}

/**
 * Constantly scaled version
 */
double rkp_term_evaluate_scaled_const( rkp_t *obj_in, int idx_in, double *x_in, double rho_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,*cf,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* If x_in is more than obj_in->dlts[idx_in] away from obj_in->pts[idx_in], then don't calculate this */
    sum = 0.0;
    for(i=0;i<obj_in->dim;i++)
        sum += pow( x_in[i] - obj_in->pts[idx_in*obj_in->dim+i], 2.0 );
    if( sum > obj_in->wrad * obj_in->dlts[idx_in] * obj_in->wrad * obj_in->dlts[idx_in] ) /* WARNING: Assumes that wfs has radius 2 * dlts[idx_in], not dlts[idx_in] */
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );
    cf = (double*) malloc( pdim * sizeof(double) );

    /* Generate the coefficients for the correction function */
    rkp_matrix_generate_scaled_const( obj_in, x_in, cf, rho_in );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = ( obj_in->pts[idx_in*obj_in->dim+i] - x_in[i] ) / rho_in; /* Important: Scale the input */

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= cf[p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );
    free( cf );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    for(i=0;i<obj_in->dim;i++)
        x[i] *= rho_in; /* Must undo dilation to give input to obj_in->wfs */
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x ); /* This is 1.0 because already divided by dilation parameter!!! */
}

/**
 * Evaluate a basis function centered on an interpolation point
 * at a given interpolation node; this function assumes that the
 * matrix has already been generated and inverted and that the
 * resulting reproducing polynomial coefficients are stored in
 * obj_in->grmlu; always call rkp_matrix_generate before calling
 * this function if a matrix with a different dilation parameter
 * is needed; this function assumes that obj_in->coeffs is filled
 * in with the polynomial coefficients of the associated correction
 * functions
 * @param obj_in The RKP object
 * @param idx_in The index of the basis function to evaluate
 * @param jdx_in The index of the interpolation node at which to evaluate the basis function
 * @return Returns the value of the function in double
 */
double rkp_term_evaluate_node( rkp_t *obj_in, int idx_in, int jdx_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Check quickly */
    sum = 0.0;
    for(i=0;i<obj_in->dim;i++)
        sum += pow( obj_in->pts[jdx_in*obj_in->dim+i] - obj_in->pts[idx_in*obj_in->dim+i], 2.0 );
    if( sum > obj_in->wrad * obj_in->dlts[idx_in] * obj_in->wrad * obj_in->dlts[idx_in] ) /* WARNING: Assumes that wfs has radius 2 * dlts[idx_in], not dlts[idx_in] */
        return 0.0;

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = obj_in->pts[idx_in*obj_in->dim+i] - obj_in->pts[jdx_in*obj_in->dim+i];

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= obj_in->coeffs[jdx_in*pdim+p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x ); /* This is 1.0 because already divided by dilation parameter!!! */
}

/**
 * Scaled version of rkp_term_evaluate_node
 */
double rkp_term_evaluate_node_scaled( rkp_t *obj_in, int idx_in, int jdx_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Check quickly */
    sum = 0.0;
    for(i=0;i<obj_in->dim;i++)
        sum += pow( obj_in->pts[jdx_in*obj_in->dim+i] - obj_in->pts[idx_in*obj_in->dim+i], 2.0 );
    if( sum > obj_in->wrad * obj_in->dlts[idx_in] * obj_in->wrad * obj_in->dlts[idx_in] ) /* WARNING: Assumes that wfs has radius 2 * dlts[idx_in], not dlts[idx_in] */
        return 0.0;

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = ( obj_in->pts[idx_in*obj_in->dim+i] - obj_in->pts[jdx_in*obj_in->dim+i] ) / obj_in->dlts[idx_in];

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= obj_in->coeffs[jdx_in*pdim+p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    for(i=0;i<obj_in->dim;i++)
        x[i] *= obj_in->dlts[idx_in];
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x ); /* This is 1.0 because already divided by dilation parameter!!! */
}

double rkp_term_evaluate_node_scaled_const( rkp_t *obj_in, int idx_in, int jdx_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Check quickly */
    sum = 0.0;
    for(i=0;i<obj_in->dim;i++)
        sum += pow( obj_in->pts[jdx_in*obj_in->dim+i] - obj_in->pts[idx_in*obj_in->dim+i], 2.0 );
    if( sum > obj_in->wrad * obj_in->dlts[idx_in] * obj_in->wrad * obj_in->dlts[idx_in] ) /* WARNING: Assumes that wfs has radius 2 * dlts[idx_in], not dlts[idx_in] */
        return 0.0;

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = ( obj_in->pts[idx_in*obj_in->dim+i] - obj_in->pts[jdx_in*obj_in->dim+i] ) / obj_in->dlts[jdx_in];

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= obj_in->coeffs[jdx_in*pdim+p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    for(i=0;i<obj_in->dim;i++)
        x[i] *= obj_in->dlts[jdx_in];
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x ); /* This is 1.0 because already divided by dilation parameter!!! */
}

/**
 * Evaluate a specified wavelet at the given point
 */
double rkp_wavelet_term_evaluate( rkp_t *obj_in, int idx_in, int wdx_in, double *x_in, double rho_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,*cf,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* If x_in is more than obj_in->dlts[idx_in] away from obj_in->pts[idx_in], then don't calculate this */
    sum = 0.0;
    for(i=0;i<obj_in->dim;i++)
        sum += pow( x_in[i] - obj_in->pts[idx_in*obj_in->dim+i], 2.0 );
    if( sum > obj_in->wrad * obj_in->dlts[idx_in] * obj_in->wrad * obj_in->dlts[idx_in] ) /* WARNING: Assumes that wfs has radius 2 * dlts[idx_in], not dlts[idx_in] */
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );
    cf = (double*) malloc( pdim * sizeof(double) );

    /* Generate the coefficients for wavelet wdx_in at point x_in */
    rkp_wavelet_generate_order( obj_in, wdx_in, x_in, cf );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = obj_in->pts[idx_in*obj_in->dim+i] - x_in[i];

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= cf[p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );
    free( cf );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x );
}

/**
 * Scaled version of rkp_wavelet_term_evaluate
 */
double rkp_wavelet_term_evaluate_scaled( rkp_t *obj_in, int idx_in, int wdx_in, double *x_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,*cf,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* If x_in is more than obj_in->dlts[idx_in] away from obj_in->pts[idx_in], then don't calculate this */
    sum = 0.0;
    for(i=0;i<obj_in->dim;i++)
        sum += pow( x_in[i] - obj_in->pts[idx_in*obj_in->dim+i], 2.0 );
    if( sum > obj_in->wrad * obj_in->dlts[idx_in] * obj_in->wrad * obj_in->dlts[idx_in] ) /* WARNING: Assumes that wfs has radius 2 * dlts[idx_in], not dlts[idx_in] */
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );
    cf = (double*) malloc( pdim * sizeof(double) );

    /* Generate the coefficients for wavelet wdx_in at point x_in */
    rkp_wavelet_generate_order_scaled( obj_in, wdx_in, x_in, cf );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = ( obj_in->pts[idx_in*obj_in->dim+i] - x_in[i] ) / obj_in->dlts[idx_in];

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= cf[p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );
    free( cf );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    for(i=0;i<obj_in->dim;i++)
        x[i] *= obj_in->dlts[idx_in];
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x );
}

/**
 * Constant scaled version of rkp_wavelet_term_evaluate
 */
double rkp_wavelet_term_evaluate_scaled_const( rkp_t *obj_in, int idx_in, int wdx_in, double *x_in, double rho_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,*cf,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* If x_in is more than obj_in->dlts[idx_in] away from obj_in->pts[idx_in], then don't calculate this */
    sum = 0.0;
    for(i=0;i<obj_in->dim;i++)
        sum += pow( x_in[i] - obj_in->pts[idx_in*obj_in->dim+i], 2.0 );
    if( sum > obj_in->wrad * obj_in->dlts[idx_in] * obj_in->wrad * obj_in->dlts[idx_in] )
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );
    cf = (double*) malloc( pdim * sizeof(double) );

    /* Generate the coefficients for wavelet wdx_in at point x_in */
    rkp_wavelet_generate_order_scaled_const( obj_in, wdx_in, x_in, rho_in, cf );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = ( obj_in->pts[idx_in*obj_in->dim+i] - x_in[i] ) / rho_in;

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= cf[p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );
    free( cf );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    for(i=0;i<obj_in->dim;i++)
        x[i] *= rho_in;
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x );
}

/**
 * Evaluate a specified wavelet at the given point
 */
double rkp_wavelet_term_evaluate_node( rkp_t *obj_in, int idx_in, int wdx_in, int jdx_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = obj_in->pts[idx_in*obj_in->dim+i] - obj_in->pts[jdx_in*obj_in->dim+i];

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= obj_in->wcoeffs[jdx_in*pdim*(pdim-1)+wdx_in*pdim+p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x );
}

/**
 * Scaled version of rkp_wavelet_term_evaluate_node
 */
double rkp_wavelet_term_evaluate_node_scaled( rkp_t *obj_in, int idx_in, int wdx_in, int jdx_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = ( obj_in->pts[idx_in*obj_in->dim+i] - obj_in->pts[jdx_in*obj_in->dim+i] ) / obj_in->dlts[idx_in];

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= obj_in->wcoeffs[jdx_in*pdim*(pdim-1)+wdx_in*pdim+p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    for(i=0;i<obj_in->dim;i++)
        x[i] *= obj_in->dlts[idx_in];
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x );
}

/**
 * Constantly scaled version of rkp_wavelet_term_evaluate_node
 */
double rkp_wavelet_term_evaluate_node_scaled_const( rkp_t *obj_in, int idx_in, int wdx_in, int jdx_in )
{
    int i,j,k,n,p,pdim,*cmb,*occ;
    double sum,prd,x[obj_in->dim];

    if( idx_in < 0 || idx_in > obj_in->np - 1 )
        return 0.0;

    /* Space dimension */
    pdim = binomial( obj_in->deg + obj_in->dim, obj_in->dim );

    /* Allocate space for polynomial counting */
    cmb = (int*) malloc( obj_in->dim * sizeof(int) );
    occ = (int*) malloc( obj_in->dim * sizeof(int) );

    /* Generate the x value taken as input relative to the particle position in pts */
    for(i=0;i<obj_in->dim;i++)
        x[i] = ( obj_in->pts[idx_in*obj_in->dim+i] - obj_in->pts[jdx_in*obj_in->dim+i] ) / obj_in->dlts[jdx_in]; /* NOTE: This is an important change! */

    /* Now evaluate the particle function */
    sum = 0.0;
    for(p=0,i=0;i<=obj_in->deg;i++)
    {
        /* Initialize the polynomial iterator */
        rcombinadic_init( i, obj_in->dim, cmb );
        n = binomial( i + obj_in->dim - 1, obj_in->dim - 1 );
        for(j=0;j<n;j++)
        {
            /* Evaluate the vector of exponents for each dimension */
            rcombinadic_occupancy( i, obj_in->dim, cmb, occ );

            /* Calculate the monomial value at x */
            prd = 1.0;
            for(k=0;k<obj_in->dim;k++)
                prd *= pow( x[k], (double) occ[k] );

            /* Multiply by the coefficient on this term from the matrix inverse */
            prd *= obj_in->wcoeffs[jdx_in*pdim*(pdim-1)+wdx_in*pdim+p++];

            /* Make sure you call rkp_matrix_generate with the right dilation parameter */
            sum += prd;

            /* Step next */
            rcombinadic_next( i, obj_in->dim, cmb );
        }
    }

    /* Clean up */
    free( cmb );
    free( occ );

    /* Call the window function to multiply the integral using the correct dilation parameter for idx_in */
    for(i=0;i<obj_in->dim;i++)
        x[i] *= obj_in->dlts[jdx_in];
    return obj_in->gqw[idx_in] * sum * obj_in->wfs( obj_in->dim, obj_in->dlts[idx_in], x );
}

/**
 * Evaluate the global interpolant at any point
 * @param obj_in The MLS object
 * @param x_in The point at which to evaluate the global interpolant
 */
double mls_evaluate( mls_t *obj_in, double *x_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * mls_term_evaluate( obj_in, i, x_in );

    return val;
}

/**
 * Evaluates the entire interpolant based on the values given
 * at the function points stored in vals
 * @param obj_in The RKP object
 * @param x_in The position at which to evaluate
 * @return Returns the value of the function
 */
double rkp_evaluate( rkp_t *obj_in, double *x_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * rkp_term_evaluate( obj_in, i, x_in );

    return val;	
}

/**
 * Scaled version of rkp_evaluate
 */
double rkp_evaluate_scaled( rkp_t *obj_in, double *x_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * rkp_term_evaluate_scaled( obj_in, i, x_in );

    return val;
}

double rkp_evaluate_scaled_const( rkp_t *obj_in, double *x_in, double rho_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * rkp_term_evaluate_scaled_const( obj_in, i, x_in, rho_in );

    return val;
}

/**
 * Evaluates the entire interpolant based on the values given
 * at the function points stored in vals; the point at which
 * to evaluate must be an interpolation node
 * @param obj_in The RKP object
 * @param jdx_in The index of the interpolation node at which to evaluate
 * @return Returns the value of the function
 */
double rkp_evaluate_node( rkp_t *obj_in, int jdx_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * rkp_term_evaluate_node( obj_in, i, jdx_in );

    return val;	
}

/**
 * Scaled version of rkp_evaluate_node
 */
double rkp_evaluate_node_scaled( rkp_t *obj_in, int jdx_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * rkp_term_evaluate_node_scaled( obj_in, i, jdx_in );

    return val;
}

/**
 * Constant scaled version of rkp_evaluate_node
 */
double rkp_evaluate_node_scaled_const( rkp_t *obj_in, int jdx_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * rkp_term_evaluate_node_scaled_const( obj_in, i, jdx_in );

    return val;
}

double rkp_wavelet_evaluate( rkp_t *obj_in, int wdx_in, double *x_in, double rho_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * rkp_wavelet_term_evaluate( obj_in, i, wdx_in, x_in, rho_in );

    return val;
}

double rkp_wavelet_evaluate_scaled( rkp_t *obj_in, int wdx_in, double *x_in, double rho_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * rkp_wavelet_term_evaluate_scaled( obj_in, i, wdx_in, x_in );

    return val;
}

/**
 * Evaluates the entire refinement interpolant based on the values given
 * at the function points stored in vals; the point at which
 * to evaluate must be an interpolation node
 * @param obj_in The RKP object
 * @param wdx_in The order of the wavelet
 * @param jdx_in The index of the interpolation node at which to evaluate
 * @return Returns the value of the function
 */
double rkp_wavelet_evaluate_node( rkp_t *obj_in, int wdx_in, int jdx_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * rkp_wavelet_term_evaluate_node( obj_in, i, wdx_in, jdx_in );

    return val;	
}

/**
 * Scaled version of rkp_wavelet_evaluate_node
 */
double rkp_wavelet_evaluate_node_scaled( rkp_t *obj_in, int wdx_in, int jdx_in )
{
    int i;
    double val;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        val += obj_in->vals[i] * rkp_wavelet_term_evaluate_node_scaled( obj_in, i, wdx_in, jdx_in );

    return val;
}

/**
 * Evaluate the inner products of pairs of functions. When evaluating
 * these inner products, remember you are only using vertices which are
 * used for generating the RKP basis! This means that you only need to
 * know reproducing polynomial coefficients for a fixed set of vertices
 * in the system. These happen to already have been saved! Don't forget!
 */
double rkp_evaluate_inner_prod( rkp_t *obj_in, int idx_in, int jdx_in )
{
    int i;
    double val,tmp;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        tmp = obj_in->gqw[i] * rkp_term_evaluate_node( obj_in, idx_in, i ) * rkp_term_evaluate_node( obj_in, jdx_in, i ), val += tmp;

    return val;
}

/**
 * Evaluate the Poisson inner products
 */
double rkp_evaluate_poisson_inner_prod( rkp_t *obj_in, int idx_in, int jdx_in )
{
    int i,j;
    double val,tmp,ttmp;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
    {
        tmp = 0.0;
        for(j=1;j<=obj_in->dim;j++)
            ttmp = obj_in->gqw[i] * rkp_wavelet_term_evaluate_node( obj_in, idx_in, j - 1, i ) * rkp_wavelet_term_evaluate_node( obj_in, jdx_in, j - 1, i ), tmp += ttmp;
        val += tmp;
    }

    return val;
}

double rkp_evaluate_functional( rkp_t *obj_in, int idx_in, double (*func_in)(int,double*) )
{
    int i;
    double val,tmp;

    val = 0.0;
    for(i=0;i<obj_in->np;i++)
        tmp = obj_in->gqw[i] * rkp_term_evaluate_node( obj_in, idx_in, i ) * func_in( obj_in->dim, obj_in->pts + i * obj_in->dim ), val += tmp;

    return val;	
}

void rkp_evaluate_gradient( rkp_t *obj_in, int idx_in, double *grad_out )
{

}

// vim: ts=4:sts=4:sw=4:et:sta