/*------------------------------------------------------------------------------*
 * File Name: NAGTutorial.C														*
 * Creation: GJL 7/20/2002														*
 * Purpose: Origin C file containing NAGTutorial								*
 * Copyright (c) OriginLab Corp.	2002-2007									*
 * All Rights Reserved															*
 *------------------------------------------------------------------------------*/

// Include needed header files 
#include <OC_const.h>         // OC Constants
#include <Common.h>           // Basic types
#include <stdio.h>            // Printf
#include <Data.h>             // Vectors and matrices
#include <NAG\OCN_g02.h>      // NAG function

/**
		NAGTutorial demonstrates how to call a NAG function.
	Example:
		int iErr;
		iErr = NAGTutorial();
	Return: 
		Returns 0 on successful exit or an error code on failure.
*/
int NAGTutorial()
{
	// Declare local variable holding error code
	int iErr;

	// *** Declare and initialize input parameters to function nt_multiple_linear_regression ***
	// *** most of which are also passed into the NAG function nag_regsn_mult_linear.        ***      
	int nPts;        // Input number of rows (or data points) in the input matrix
	int nTdx;        // Input number of columns (or independent variables) in the input matrix
	matrix mX;       // Input matrix containing data points of the independent variables
	vector vY;       // Input vector containing data points of dependent variable
	vector vWT;      // Input vector containing weights of data points, NULL if weighting not used

	// There are 10 rows of data in sample data
	nPts = 10;
	
	// There are 3 independent variables in data
	nTdx = 3;
	
	// Get independent variable data sets
	Worksheet wks("Data1");           // Attach OC Worksheet object to Data1 worksheet
	mX.CopyFromWks( wks, 1, 3, 0, 9); // Copy independent variables from columns 2 to 4 and rows
	                                  // 1 to 10 of worksheet Data1 into matrix
	
	// Get dependent data set
	Dataset ds( "Data1_A" );          // Attach to Data1_A data set         
	vY = ds;                          // Copy dependent variable from Dataset to vector
	
	// Set weights - this example uses all weights of 1
	vWT.SetSize( nPts );
	vWT = 1;
	
	// *** Declare output parameters to nt_multiple_linear_regression most of which are also ***
	// *** returned from the NAG function nag_regsn_mult_linear.                             ***	
	double dRss;     // Output residual sum of squares for the regression 
	double dDf;      // Output degrees of freedom associated with dRss
	vector vB;       // Output least-squares estimates of the parameters of the regression model
	vector vSE;      // Output standard errors of the parameter estimates given in vB
	vector vCOV;     // Output variance-covariance matrix of estimated parameters in vB
	vector vRES;     // Output weighted residuals
	vector vH;       // Output diagonal elements of H
	matrix mQ;       // Output results of the QR decomposition
	BOOL bSvd;       // Output flag set TRUE if singular value decomposition was performed, otherwise FALSE
	int iRank;       // Output rank of the independent variables
	vector vP;       // Output details of QR decomposition if SVD is used
	vector vCOMAR;   // Output information which is needed by nag_regsn_mult_linear_newyvar_if bSvd is TRUE 

	// *** Call generalized Origin C function that calls NAG multiple linear regression function ***
	iErr = nt_multiple_linear_regression( nPts, nTdx, mX, vY, vWT, dRss, dDf, vB, vSE, vCOV, vRES, vH, mQ, bSvd, iRank, vP, vCOMAR );

	// *** Print out equation with parameter values ***
	printf( "Equation: Y = %g + %g*X1 + %g*X2 + %g*X3", vB[0], vB[1], vB[2], vB[3] );

	return iErr;
}

/**
		Function to perform multiple linear regression using the NAG function nag_regsn_mult_linear.
	Parameters:
		nPts=Input number of rows (or data points) in the input matrix
		nTdx=Input number of columns (or independent variables) in the input matrix
		mX=Input matrix containing data points of the independent variables
		vY=Input vector containing data points of dependent variable
		vWT=Input vector containing weights of data points, NULL if weighting not used
		dRss=Output residual sum of squares for the regression 
		dDf=Output degrees of freedom associated with dRss
		vB=Output least-squares estimates of the parameters of the regression model
		vSE=Output standard errors of the parameter estimates given in vB
		vCOV=Output variance-covariance matrix of estimated parameters in vB
		vRES=Output weighted residuals
		vH=Output diagonal elements of H
		mQ=Output results of the QR decomposition
		bSvd=Output flag set TRUE if singular value decomposition was performed, otherwise FALSE
		iRank=Output rank of the independent variables
		vP=Output details of QR decomposition if SVD is used
		vCOMAR=Output information which is needed by nag_regsn_mult_linear_newyvar_if bSvd is TRUE 
		dTol=Input tolerance used to decide rank of independent variables, default is 0.000001 
		bThroughZero=Input flag forcing regression line through origin, default is FALSE
	Return:
		Returns the results of the NAG nag_regsn_mult_linear function.
*/
int nt_multiple_linear_regression( int nPts, int nTdx, matrix& mX, vector& vY, vector& vWT,
	double& dRss, double& dDf, vector& vB, vector& vSE, vector& vCOV, vector& vRES,
	vector& vH, matrix& mQ, BOOL& bSvd, int& iRank, vector& vP,	vector& vCOMAR,
	double dTol = 0.000001, BOOL bThroughZero = FALSE )
{
	// Declare local variable holding error code
	int iErr;

	// *** Declare, size, prepare, and check all NAG Linear Regression arguments (as needed) ***
	// Force regression line through zero or not...
	int iP;
	Nag_IncludeMean iMean;                     // NAG variable type from header file \OriginC\system\NAG\nag_types.h
	if( bThroughZero )                         // If input flag says to force through zero...
	{
		iMean = Nag_MeanZero;                  // Set NAG constant to force regression line through zero
		iP = nTdx;                             // Set number P of independent variables in the model not including the mean (or intercept)		
	}
	else
	{
		iMean = Nag_MeanInclude;               // Set NAG constant to NOT force regression line through zero
		iP = nTdx + 1;                         // Set number P of independent variables in the model including the mean
	}

	// NAG nag_regsn_mult_linear function requires nPts > 1
	if( nPts < 2 )                             // If less than 2 points return error
		return 1;

	// The number of independent variables equals the second dimension of the matrix
	int nM = nTdx;
	
	// Indicates which of the potential independent variables are to be included in the model...use all so set to 1
	vector<int> vSX;
	vSX.SetSize( nM );
	vSX = 1;
	
	// The least-squares estimates of the parameters of the regression model
	vB.SetSize( iP );
	
	// The standard errors of the iP parameter estimates given in B
	vSE.SetSize( iP );
	
	// The variance-covariance matrix of estimated parameters in B
	vCOV.SetSize( iP * ( iP + 1 ) / 2);
	
	// The (weighted) residuals
	vRES.SetSize( nPts );
	
	// The diagonal elements of H
	vH.SetSize( nPts );
	
	// The second dimension of the array Q
	int nTdq;
	nTdq = iP + 1;
	
	// The results of the QR decomposition
	mQ.SetSize( nPts, nTdq );
	
	// Details of the QR decomposition and SVD if used
	vP.SetSize( 2 * iP + iP * iP );
	
	// Information which is needed by nag_regsn_mult_linear_newyvar_if bSvd is TRUE
	vCOMAR.SetSize( 5 * ( iP - 1 ) + iP * iP );
	
	// *** Call NAG function nag_regsn_mult_linear to perform linear regression ***
	// From <NAG\OCN_g02.h>: g02dac nag_regsn_mult_linear: Perform multiple linear regression
	iErr = nag_regsn_mult_linear( iMean, nPts, mX, nTdx, nM, vSX, iP, vY, vWT, &dRss, &dDf,
		vB, vSE, vCOV, vRES, vH, mQ, nTdq, &bSvd, &iRank, vP, dTol, vCOMAR );

	return iErr;
}