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

/*OLS______________________________________________________

  Pass 1 argument from Matlab to C : 
  	B : 2x1 or 1x2

  Return 1 scalar to Matlab: Sum (Y - B0 - B1*X)^2

  Syntax : n = OLS(B)
  Example:
	
	>> options = optimset('TolFun',1e-10) ;
	>> fminsearch('OLS',[1;1],options)

	ans =

		2.99999900523767
		2.00000010762575
*/

double Function(double *B); 

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

  //Check input and output
  if (nrhs != 1)
	  mexErrMsgTxt("One input required!");
  else if(nlhs > 1)
	  mexErrMsgTxt("Too many output arguments. One is enough.");

  int mrowsB = mxGetM(prhs[0]);
  int ncolsB = mxGetN(prhs[0]);
  if(!(mrowsB==1 && ncolsB==2 || mrowsB==2 && ncolsB==1))
	  mexErrMsgTxt("Input must by 2x1 or 1x2.");

  //Retrieving input
  B = mxGetPr(prhs[0]);
  
  //Setup output
  plhs[0] = mxCreateDoubleMatrix(1,1, mxREAL);
  N  = mxGetPr(plhs[0]);
  
  //Compute
  *N  = Function(B) ;
  
  return;
}


double Function(double *B) 
{
	double norm = 0.0 ;	// Final output

	// Create the data. This could be a complicated model.
	// Example: B are preferences and technology parameters.
	//			X is a path for prices
	//			2 X + 3 is model computations
	//          Y is actual data.
/*	for (int i=0;i<100;i++)			
	{							
		X[i] = i ;				
		Y[i] = 2 * X[i] + 3 ;	
	}	
*/	
	mxArray * Xarray = mexGetVariable("base","X");
	mxArray * Yarray = mexGetVariable("base","Y");

	int mrowsX = mxGetM(Xarray);
	int ncolsX = mxGetN(Xarray);
	int size;
	if (mrowsX == 1)
		size = ncolsX;
	else
		size = mrowsX;
	double * X = mxGetPr(Xarray);
	double * Y = mxGetPr(Yarray);
	
	for (int i=0;i<100;i++) 
	{
		// Distance between model and data.
		norm += pow ( Y[i] - X[i] * B[1] - B[0] , 2 ) ;  
	}
	mexPrintf("norm = %e\n", norm);	
	return norm ;
}