/*------------------------------------------------------------------------------*
 * File Name: CallingNAGFunctionsEx												*
 * Creation: GJL 02/01/05														*
 * Purpose: OriginC Source C file containing NAG examples.						*
 * Copyright (c) OriginLab Corp.	2005, 2006, 2007, 2008, 2009, 2010			*
 * All Rights Reserved															*
 * 																				*
 *------------------------------------------------------------------------------*/

#include <Origin.h>
#include <OC_nag8.h>

// Example program calling simple NAG function (Cosh)
double SimpleNAGEx(double dx)
{
	NagError neNagError;
	double dCosh = nag_cosh(dx, &neNagError);
	if( NE_NOERROR == neNagError.code )
	{
		printf("Using NAG Cosh(%g) = %g\n", dx, dCosh);
		return dCosh;
	}
	printf("Error computing Cosh(%g) using NAG\n", dx);
	return NANUM;
}

// Example program calling NAG function having one-dimensional array arguments
int NAGOneDimesionArrayEx(string strDatasetName)
{
	int nPts, nValid;
	double dMean, dSD, dSkewness, dKurtosis, dMin, dMax, dWeightSum;
	NagError neFail;

	DataRange dr;
	Worksheet wks = Project.ActiveLayer(); //Get the active layer
	dr.Add("X", wks, 0, 0, -1, 0);		//Get the first column as DataRange "X"

	if(dr.IsValid())
	{
		// Calling NAG function weighting data with 1 and Origin missing values with 0
		vector vData1, vWeights;
		dr.GetData(&vData1, 0);		//get the datarange data
		dr.GetData(&vWeights, 0);
		
		vWeights.Replace(NANUM, 1); // Replace non-NANUMs (valid data) with 1
		vWeights.Replace(NANUM, 0, MATREPL_TEST_EQUAL); // Replace NANUMs with 0
		nPts = vData1.GetSize();
		nag_summary_stats_1var(nPts, vData1, vWeights, &nValid, &dMean, &dSD,
			&dSkewness, &dKurtosis, &dMin, &dMax, &dWeightSum, &neFail);
		if( !neFail.code ) // If 0 then success so print results 
			printf("\nnPts = %d\nnValid = %d\ndMean = %g\ndSD = %g\ndSkewness = %g\ndKurtosis = %g\ndMin = %g\ndMax = %g\ndWeightSum = %g\n",
				nPts, nValid, dMean, dSD, dSkewness, dKurtosis, dMin, dMax, dWeightSum);
		else
			printf("Error computing summary stats using NAG\n");

		 //Calling NAG function first removing Origin missing values from data
		vector vData2;         
		dr.GetData(&vData2, 0);
		vData2.Trim();             //remove the missing values
		nPts = vData2.GetSize();
		nag_summary_stats_1var(nPts, vData2, NULL, &nValid, &dMean, &dSD,
			&dSkewness, &dKurtosis, &dMin, &dMax, &dWeightSum, &neFail);
		if( !neFail.code ) // If 0 then success so print results 
			printf("\nnPts = %d\nnValid = %d\ndMean = %g\ndSD = %g\ndSkewness = %g\ndKurtosis = %g\ndMin = %g\ndMax = %g\ndWeightSum = %g\n",
				nPts, nValid, dMean, dSD, dSkewness, dKurtosis, dMin, dMax, dWeightSum);
		else
			printf("Error computing summary stats using NAG\n");

		return 0;
	}
	printf("Error attaching to dataset %s\n", strDatasetName);
	return 1;
}

// Example program calling NAG function having two-dimensional array arguments
// Assumes that active worksheet contains dependent variable Y in column 1 and independent variables
//   in columns 2 to N
int NAGTwoDimesionArrayEx()
{
	// *** Declare, size, prepare, and check all NAG Linear Regression arguments ***
	// Input flag to force regression line through zero (or not)
	Nag_IncludeMean iMean; // NAG variable type from header file \OriginC\system\NAG8\nag_types.h
	iMean = Nag_MeanZero;  // Set NAG constant to force regression line through zero

	// Attach Worksheet object to active layer (worksheet)
	Worksheet wks = Project.ActiveLayer();
	if( !wks.IsValid() ) // If wks not valid return with error
	{
		printf("Error: Active page must be a worksheet.\n");
		return 1;
	}

	// Input matrix (dynamic 2D array) holding independent variables
	matrix mX;
	mX.CopyFromWks(wks, 1, wks.GetNumCols()); // Copy independent variables from columns 2 to N of wks
	if( mX.RemoveEmptyRows(FALSE) ) // If any missing values in independent variables then error
	{
		printf("Error: Input independent variables contain missing values.\n");
		return 1;
	}

	// Input number of rows (or data vales) in the input matrix
	int nPts; // NAG nag_regsn_mult_linear function requires nPts > 1
	nPts = mX.GetNumRows();
	if( nPts < 2 ) // If less than 2 data vales return error
	{
		printf("Error: Input independent variables must have 2 or more data values.\n");
		return 1;
	}

	// Input number of columns (or independent variables) in the input matrix
	int nTdx;
	nTdx = mX.GetNumCols();	

	// Input number of independent variables equals the second dimension of the matrix
	int nM = nTdx;
	if( nM < 1 ) // If less than 1 independent variable return error
	{
		printf("Error: There must be 1 or more input independent variables.\n");
		return 1;
	}

	// Input vector indicating which of the potential independent variables are to be included in the model
	vector<int> vSX;
	vSX.SetSize(nM);
	vSX = 1; // Use all so set entire vector to 1

	// Input number of independent variables
	int iP;
	iP = nTdx; // Set number P of independent variables in the model not including the mean or intercept

	// Input vector containing data points of dependent variable
	matrix mY; // Use 1 column matrix for vector...provides access to CopyFromWks and RemoveEmptyRows methods 
	mY.CopyFromWks(wks, 0, 0); // Copy dependent variable from column 1 of wks
	if( mY.RemoveEmptyRows(FALSE) ) // If any missing values in independent variables then error
	{
		printf("Error: Input dependent variable contains missing values.\n");
		return 1;
	}
	if( mY.GetNumRows() != nPts ) // Dependent and independent variable must have same number of data values
	{
		printf("Error: Dependent and independent variables must have same number of data values.\n");
		return 1;
	}

	// Output residual sum of squares for the regression 
	double dRss;

	// Output degrees of freedom associated with dRss
	double dDf;

	// Output least-squares estimates of the parameters of the regression model
	vector vB;
	vB.SetSize(iP);

	// Output standard errors of the parameter estimates given in vB
	vector vSE;
	vSE.SetSize(iP);

	// Output variance-covariance matrix of estimated parameters in vB
	vector vCOV;
	vCOV.SetSize(iP * ( iP + 1 ) / 2);

	// Output weighted residuals
	vector vRES;
	vRES.SetSize(nPts);

	// Output diagonal elements of H
	vector vH; 
	vH.SetSize(nPts);

	// Input second dimension of the array Q
	int nTdq;
	nTdq = iP + 1;

	// Output results of the QR decomposition
	matrix mQ; 	
	mQ.SetSize(nPts, nTdq);

	// Output flag set TRUE if singular value decomposition (SVD) was performed, otherwise FALSE
	BOOL bSvd;

	// Output rank of the independent variables
	int iRank;

	// Output details of QR decomposition if SVD is used
	vector vP;
	vP.SetSize(2 * iP + iP * iP);

	// Input Tolerance used to determine the rank of the independent variables
	double dTol = 0.000001;

	// Output information which is needed by nag_regsn_mult_linear_newyvar if SVD is used
	vector vCOMAR;
	vCOMAR.SetSize(5 * ( iP - 1 ) + iP * iP);

	// Output NAG error structure
	NagError neFail;

	// *** Call NAG function nag_regsn_mult_linear to perform linear regression ***
	// From <NAG8\nagg02.h>: g02dac nag_regsn_mult_linear: Perform multiple linear regression
	nag_regsn_mult_linear(iMean, nPts, mX, nTdx, nM, vSX, iP, mY, NULL, &dRss, &dDf,
		vB, vSE, vCOV, vRES, vH, mQ, nTdq, &bSvd, &iRank, vP, dTol, vCOMAR, &neFail);

	// *** Output results of linear regression ***
	if( neFail.code )
		printf("Error: nag_regsn_mult_linear error =  %d.\n", neFail.code);
	else
	{
		printf("\nY = %g*X1", vB[0]); // Assume at least 1 independent variable
		for( int ii = 1; ii < nTdx; ii++ ) // Print remaining independent variables
			printf(" + %g*X%d", vB[ii], ii + 1);
		printf("\n");
	}

	return neFail.code;
}
