/*
interpolate_rule.c
Include description here later...
*/

#include "mex.h"

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

/* 
  Check for proper number and formatting of inputs and outputs
  Input 1: ordered sample x
  Input 2: sample y
  Input 3: ordered desired xd
  Input 4: check whether ordered or not
  Output 1: yd corresponding to desired xd
*/

  if (nrhs != 4 || nlhs != 1) {
    mexErrMsgIdAndTxt("MATLAB:amatmult:invalidNumInputOrOutput", "Four input arguments and one output argument required");
  }
  
  if (mxGetN(prhs[0]) > 1 || mxGetN(prhs[1]) > 1 || mxGetN(prhs[2]) > 1){
  	mexErrMsgIdAndTxt("MATLAB:amatmult:invalidInputFormat", "only vector inputs desired");
  }
  
  mwSize n = mxGetM(prhs[0]);
  mwSize nd = mxGetM(prhs[2]);
  
  if (mxGetM(prhs[1]) != n){
  	mexErrMsgIdAndTxt("MATLAB:amatmult:invalidInputFormat", "sample x and y should be same length");
  }
  
  if (!mxIsLogical(prhs[3])){mexErrMsgIdAndTxt("MATLAB:amatmult:invalidInputFormat", "fourth input should be logical");}
  if(!(mxIsDouble(prhs[0]) && mxIsDouble(prhs[1]) && mxIsDouble(prhs[2]))){mexErrMsgIdAndTxt("MATLAB:amatmult:invalidInputFormat", "first three inputs should be numeric");}
  
  double * xgrid = mxGetPr(prhs[0]);
  double * ygrid = mxGetPr(prhs[1]);
  double * xdgrid = mxGetPr(prhs[2]);
  
  plhs[0] = mxCreateDoubleMatrix(nd,1,mxREAL);
  double * ydgrid = mxGetPr(plhs[0]);
    
  mwIndex i;
    
  /*
  Test for order if flagged
  */
  int * flagpointer = (int*) mxGetLogicals(prhs[2]);
  int flag = flagpointer[0];
    
  if(flag==1){
  	for(i=0; i<n-1; i++){
  		if(xgrid[i] >= xgrid[i+1]){
  			mexErrMsgIdAndTxt("MATLAB:amatmult:invalidInputFormat", "sample x should be strictly increasing");
  		}
  	}
  	for(i=0; i<nd-1; i++){
  		if(xdgrid[i] >= xdgrid[i+1]){
  			mexErrMsgIdAndTxt("MATLAB:amatmult:invalidInputFormat", "desired x should be strictly increasing");
  		}
  	}
  }
  
  /*
  Now interpolate for every point
  */
  mwIndex cur_xi=0;
  mwIndex act_xi;
  double cur_xd, share_xlow;
  
  for(i=0; i<nd; i++){
  	cur_xd = xdgrid[i];
  	while(cur_xi<n-1){
  			if(xgrid[cur_xi] >= cur_xd){break;}
  			cur_xi++;
  	}
  	act_xi = cur_xi-1;
  	if(act_xi == -1){act_xi = 0;}
  	share_xlow = (xgrid[act_xi+1]-cur_xd)/(xgrid[act_xi+1]-xgrid[act_xi]);
  	ydgrid[i] = share_xlow*ygrid[act_xi] + (1-share_xlow)*ygrid[act_xi+1]; 
  }
}