// MBEstNLPart.c	
// Using the current estimated parameters to estimate the nonlinear 
// part of the semiparametric regression, i.e., B(x)'*beta.
//
//   The algorithm for MicroBayes is described in the following paper
//
//   D. Zhang, M.T. Wells, C.D. Smart and W.E. Fry (2003). Bayesian
//       Normalization and Identification for Differential Gene
//       Expression Data.
//
//   PLASE CITE THIS PAPER AFTER YOU USE THIS SOFTWARE.
//
//	Dabao Zhang		February 24, 2003
//	Revised by Dabao Zhang
//       March 9, 2003 -- Consider genes' location effect
//
//	Copyright (c) 2003 by Dabao Zhang.

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

void MBEstNLPart(int degree,int nobs,int nknots,double knots[],
                 double x[],double beta[],double eNL[])
{
    int n, k;

    for(n=0; n<nobs; n++)
    {
        eNL[n] = beta[0];
        for(k=0; k<degree; k++)
        {
            eNL[n] = eNL[n] + beta[k+1]*pow(x[n],k+1);
        }
        
        for(k=0; k<nknots; k++)
        {
            if( x[n]>knots[k] )
            {
                eNL[n] = eNL[n] + beta[k+degree+1]*pow(x[n]-knots[k],degree);
            }
        }
    }
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *eNL,*knots,*beta,*x;
    int degree,nobs,nknots;

    // Check for proper number of arguments.
    if( nrhs!=3 )
    {
        mexErrMsgTxt("Three inputs required.");
    }
    else if( nlhs>1 )
    {
        mexErrMsgTxt("Too many output arguments.");
    }
    
    // RHS: config, iota=x, beta
    degree = (int)mxGetScalar(mxGetField(prhs[0],0,"degree"));
    nobs = (int)mxGetScalar(mxGetField(prhs[0],0,"nobs"));
    nknots = (int)mxGetScalar(mxGetField(prhs[0],0,"nknots"));
    knots = mxGetPr(mxGetField(prhs[0],0,"knots"));
    
    x = mxGetPr(prhs[1]);

    beta = mxGetPr(prhs[2]);
    
    plhs[0] = mxCreateDoubleMatrix(nobs,1,mxREAL);
    eNL = mxGetPr(plhs[0]);

    MBEstNLPart(degree,nobs,nknots,knots,x,beta,eNL);
}