/*------------------------------------------------------------------------------*
 * File Name:matrix.c				 											*
 * Creation:July 22th, 2002														*
 * Purpose: OriginC Source C file												*
 * Copyright (c) OriginLab Corp.2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007	*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 * TCZ 07/29/02 QA70-2497	ADD_GLOBAL_STATS_SUMMARY_FUNCTION					*
 * ER 01/29/03 QA70_3802	ADD_VECTOR_FFT_IFFT
 * CPY 7/7/03 changed to include only origin.h to take advantage of PCH			
 *------------------------------------------------------------------------------*/
 
#include <origin.h> // main Origin C header that is precompiled and already include most headers 

// this file include most of the other header files except the NAG header, which takes longer to compile
// NAG routines
#include <NAG\OCN_c06.h> // this contains the FFT NAG headers, 
#include <NAG\OCN_f.h>//this contain the linear algebra function
#include <NAG\OCN_g01.h>//this file contain the simple statistics function
////////////////////////////////////////////////////////////////////////////////////
#include <matrix.h>
// start your functions here


//local functions for the internal use

 
// This function will do the 2d fft with the data form required by nag funtion.
static int fft2_for_complex_matrix(int nNumRows, int nNumCols,  matrix& matReal, matrix& matImage, matrix<complex>& matFFT2Result);


//  	This function will take a matrix and set the result matrix accoring to the parameter values.
//	in other words, it will do the neccessory padding with zero or truncting.
static int set_matrix_with_padding_truncting(matrix & matSource, matrix & matResult, int nRowProcess , int nColProcess);


//   This function will take care of the complex data matrix and sperate it into two double matrix.
static int seperate_complex_mat_to_real_and_imag_mat(matrix<complex>& matSource, matrix& matReal, matrix&  matImage, int nRowProcess, int nColProcess);

/**

	int FFT2(matrix& matSource, matrix<complex>& matFFT2Result, int nRowProcess  = -1, int nColProcess = -1);
	int FFT2(matrix<complex>& matSource, matrix<complex>& matFFT2Result, int nRowProcess  = -1, int nColProcess = -1);
			These two Functions will do the 2 dimension Fast Fourier Transform.
	Remarks:
			Right now, we only support the data matrix in double and complex form. If user have the data
			matrix in integer form, user should first transform it to double. It is a member function of matrix.
			Also, if user has a vector data, it should firstly be transformed into a matrix and then do the 
			FFT.
	Example:
	Parameters:
			matSource: The data source matrix
			matFFT2Result: The matrix containing the 2 dimension FFT results
			nRowProcess:
			nColProcess: The number of rows or columns user provided. If they are -1, there is no padding
	 					or truncating. If they are more than the sizes of the matrix, the function will pad zero to meet the size
						If less than, it will truncate the data matrix.
	Return:
			If succeed, it will return 0;	Otherwise it returns error indications.
			Error:
			-1: Provided column number or row number. Must be -1 or positive integer
			-2: The data source matrix is empty 
			-3: When calling member function, it fails.
			-4: When calling GetReal or GetImage member function of matrix class, it failed.
			Nag Error:
			 11: The column or the row size of matrix is less then 1
 			 73: System failed to allocate memory

*/
int FFT2(matrix& matSource, matrix<complex>& matFFT2Result, int nRowProcess, int nColProcess)
{
	matrix matReal;
	matrix matTemp(matSource);
	int nRet;
	if(nRet = set_matrix_with_padding_truncting(matTemp, matReal, nRowProcess, nColProcess))
		return nRet;
	matrix matImage;
		
	int nNumCols = matReal.GetNumCols();
	int nNumRows = matReal.GetNumRows();
	matImage.SetSize(nNumRows, nNumCols);
	matImage = 0;	
	if(nRet = fft2_for_complex_matrix(nNumRows, nNumCols, matReal, matImage, matFFT2Result))
		return nRet;
	return nRet;
}

static int fft2_for_complex_matrix(int nNumRows, int nNumCols,  matrix& matReal, matrix& matImage, matrix<complex>& matFFT2Result)
{
	vector vecTrigRows;
	vecTrigRows.SetSize(2 * nNumRows);
	//cal nage function to prepare the trig value
	int nRet;
	if( nRet = nag_fft_init_trig(nNumRows, vecTrigRows))
		return nRet;
	
	vector vecTrigCols;
	vecTrigCols.SetSize(2 * nNumCols);
	if( nRet = nag_fft_init_trig(nNumCols, vecTrigCols))
		return nRet;
	
	//call the nag function to do the fft2
	if(nRet = nag_fft_2d_complex(nNumRows, nNumCols, matReal, matImage, vecTrigRows, vecTrigCols))
		return nRet;
	//set the data back
	if( nRet = matFFT2Result.MakeComplex(matReal, matImage))
		return nRet;
	
}
		
int FFT2(matrix<complex>& matSource, matrix<complex>& matFFT2Result, int nRowProcess, int nColProcess )
{
	int nRet;
	matrix matReal, matImage;
	if(nRet = seperate_complex_mat_to_real_and_imag_mat(matSource, matReal, matImage, nRowProcess, nColProcess))
		return nRet;
	
	int nNumCols = matImage.GetNumCols();
	int nNumRows = matImage.GetNumRows();
	if(nRet = fft2_for_complex_matrix(nNumRows, nNumCols, matReal, matImage, matFFT2Result))
		return nRet;
	
}
		
static int set_matrix_with_padding_truncting(matrix & matSource, matrix & matResult, int nRowProcess , int nColProcess)
{
		
	
	if(nRowProcess < -1 || nRowProcess == 0 ||nColProcess < -1 ||nColProcess == 0)
		return FFT_ERROR_COULUMN_ROW_NUMBER_NEGATIVE;
	int nNumCols = matSource.GetNumCols();
	int nNumRows = matSource.GetNumRows();
	if(nNumCols == 0 || nNumRows == 0)
		return FFT_ERROR_SOURCE_MATRIX_EMPTY;
		
	int nIntendNumCols = nNumCols;
	int nIntendNumRows = nNumRows;
	if(nRowProcess != -1) 
		nIntendNumRows = nRowProcess;
	if(nColProcess != -1)
		nIntendNumCols = nColProcess;
	
	int nRet;
	matResult = matSource;
	
	if(nIntendNumRows > nNumRows)//so we should increase the dimension and do the padding
	{
		matrix matTemp;
		matTemp.SetSize(nIntendNumRows - nNumRows, nNumCols );
		matTemp = 0;
		if(nRet = matSource.Concatenate(1,matTemp, matResult))
			return nRet;
	}
	BOOL bOK = TRUE;
	if (nIntendNumRows < nNumRows)
    	bOK = matSource.GetSubMatrix(matResult, 0, -1, 0, nIntendNumRows - 1);
	if(!bOK)
		return FFT_ERROR_GETSUBMATRIX_FUNCTION_ERROR;
    		
	    
    if(nIntendNumCols > nNumCols)
	{
		matrix matTemp;
		matTemp.SetSize(nIntendNumRows, nIntendNumCols - nNumCols );
		matTemp = 0;
		if(nRet = matResult.Concatenate(2, matTemp, matResult))
			return nRet;
		
	}
	
	if(nIntendNumCols < nNumCols)
	{
		matrix matResultTemp(matResult); 
		bOK = matResultTemp.GetSubMatrix(matResult, 0, nIntendNumCols -1, 0, -1);
		if(!bOK)
			return FFT_ERROR_GETSUBMATRIX_FUNCTION_ERROR;
	}
	return 0;
}

/*
	int IFFT2(matrix<complex>& matSource, matrix<complex>& matIFFT2Result, int nRowProcess  = -1, int nColProcess = -1);
	Remarks:
			Right now, we only support the data matrix in complex form. If user have the data
			matrix in double form, user should add the imaginary part as zero to make a complex matrix.
			Matrix class has a member function to make the complex matrix. Also, if user has a vector data,
			it should firstly be transformed into a matrix and then one can do the inverse FFT.
	Example:
	Parameters:
			matSource: The data source matrix
			matIFFT2Result: The matrix containing the 2 dimension inverse FFT results
			nRowProcess:
			nColProcess:The number of rows or columns user provided. if they are -1, there are no padding
	 					or truncating. if they are larger than the size of the matrix, the function will pad zero to meet
						the size. If less than, it will truncate the data matrix.
	Return:
			If succeed, it will return 0;	Otherwise it returns error indications.
			Error:
			-1: Provided column number or row number. Must be -1 or positive integer
			-2: The data source matrix is empty 
			-3: When calling member function, it fails.
			-4: When calling GetReal or GetImage member function of matrix class, it failed.
			Nag Error:
			 11: The column or the row size of matrix is less then 1
 			 73: System failed to allocate memory

  
*/
int IFFT2(matrix<complex>& matSource, matrix<complex>& matFFT2Result, int nRowProcess, int nColProcess)
{

	int nRet;
	matrix<complex> matSourceCopy(matSource);
	if(nRet = matSourceCopy.Conjugate())
		return nRet;
	
	if(nRet = FFT2(matSourceCopy, matFFT2Result,nRowProcess, nColProcess ))
		return nRet;
	
	if(nRet = matFFT2Result.Conjugate())
		return nRet;
}		
/*
	int FFT(matrix& matSource, matrix<complex>& matFFTResult, int nColProcess = -1);
	int FFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess = -1);

				These two Functions will do the 1 dimension Fast Fourier Transform.
  	Remarks:
			Right now, we only support the data matrix in double and complex form. If user has the data
			matrix in integer form, user should firstly transform it to double. It is a member function of matrix.
			Also, if user has a vector data, he should firstly transform data into a matrix and then one can do the 
			FFT. Unlike the matlab, this function will do the FFT on each row (Matlab do FFT on each column).
	Example:
	Parameters:
			matSource: The data source matrix
			matFFTResult: The matrix containing the 1 dimension FFT results
			nColProcess: The number of columns user provided. If it is -1, function will not do any padding
	 					or truncating. But if it is larger than the size of the  matrix, it will pad zero to meet the size
						if less than, it will truncate the data matrix.
	Return:
			If succeed, it will return 0;	Otherwise it returns error indications.
			Error:
			-1: Provided column number or row number. Must be -1 or positive integer
			-2: The data source matrix is empty 
			-3: When calling member function, it fails.
			-4: When calling GetReal or GetImage member function of matrix class, it failed.
			Nag Error:
			 11: The column or the row size of matrix is less then 1
 			 73: System failed to allocate memory

*/

int FFT(matrix& matSource, matrix<complex>& matFFTResult, int nColProcess)
{
	matrix matTemp;
	int nRet;
	if(nRet = set_matrix_with_padding_truncting(matSource, matTemp, -1, nColProcess))
		return nRet;
	int nNumCols = matTemp.GetNumCols();
	int nNumRows = matTemp.GetNumRows();
	vector vecTrig;
	vecTrig.SetSize(2 * nNumCols); 
	if(nRet = nag_fft_init_trig(nNumCols, vecTrig))
		return nRet;
	if(nRet = nag_fft_multiple_real(nNumRows, nNumCols, matTemp, vecTrig))
		return nRet;
	matrix matReal, matImage;
	matReal.SetSize(nNumRows, nNumCols);
	matImage.SetSize(nNumRows, nNumCols);
	
	if(nRet = nag_multiple_hermitian_to_complex(nNumRows, nNumCols, matTemp, matReal, matImage))
		return nRet;
	if(nRet = matFFTResult.MakeComplex(matReal, matImage))
		return nRet;
	return nRet;
}
		
int FFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess)
{
	int nRet;
	matrix matReal, matImage;
	if(nRet = seperate_complex_mat_to_real_and_imag_mat(matSource, matReal, matImage, -1, nColProcess))
		return nRet;
	int nNumCols = matImage.GetNumCols();
	int nNumRows = matImage.GetNumRows();
	vector vecTrig;
	vecTrig.SetSize(2 * nNumCols); 
	if(nRet = nag_fft_init_trig(nNumCols, vecTrig))
		return nRet;
	if(nRet = nag_fft_multiple_complex(nNumRows, nNumCols, matReal, matImage, vecTrig ))
		return nRet;
	if(nRet = matFFTResult.MakeComplex(matReal, matImage))
		return nRet;
	return nRet;
}
/*	  
	int IFFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess = -1);
	
			This Function will do the 1 dimension inverse Fast Fourier Transform.
	Remarks:
			Right now, we only support the data matrix in complex form. If user has the data
			matrix in double form, user should add the imaginary part as zero to make a complex matrix.
			Matrix class has a member function to make the complex matrix. also, if user has a vector data,
			it should be firstly transformed into a matrix and then one can do the inverse FFT.

	Example:
	Parameters:
			matSource:The data source matrix
			matIFFT2Result: The matrix containing the 1 dimension inverse FFT results
			nColProcess:The number of columns user provided. if it is -1, there is no padding
	 					or truncting. if it is more than the size of the  matrix, it will pad zero to meet the size
						if less than, it will truncate the data matrix.
	Return:
			If succeed, it will return 0;	Otherwise it returns error indications.
			Error:
			-1: Provided column number or row number must be -1 or positve integer
			-2: The data source matrix is empty 
			-3: When calling member function, it fails.
			-4: When calling GetReal or GetImage member function of matrix class, it failed.
			Nag Error:
			 11: The coloumn or the row size of matrix is less then 1;
 			 73: System failed to allocate memory

  
*/


		
int IFFT(matrix<complex>& matSource, matrix<complex>& matFFTResult, int nColProcess)
{
	int nRet;
	matrix<complex> matSourceCopy(matSource);
	if(nRet = matSourceCopy.Conjugate())
		return nRet;

	if(nRet = FFT(matSourceCopy, matFFTResult, nColProcess))
		return nRet;

	if(nRet = matFFTResult.Conjugate())
		return nRet;	
		
}	
	
static int seperate_complex_mat_to_real_and_imag_mat(matrix<complex>& matSource, matrix& matReal, matrix&  matImage, int nRowProcess, int nColProcess)
{
	matrix matTemp;
	BOOL bOK;
	bOK = matSource.GetReal(matTemp);
	if(!bOK)
		return FFT_ERROR_GETREAL_FUNCTION_FAIL;
	int nRet;
	if(nRet = set_matrix_with_padding_truncting(matTemp, matReal, nRowProcess, nColProcess))
		return nRet;
	bOK = matSource.GetImaginary(matTemp);
	if(!bOK)
		return FFT_ERROR_GETREAL_FUNCTION_FAIL;
	if(nRet = set_matrix_with_padding_truncting(matTemp, matImage, nRowProcess, nColProcess))
		return nRet;
}

// ER 01/29/03 QA70_3802	ADD_VECTOR_FFT_IFFT
// The following function performs 1-D FFT on a complex vector.
// It uses the matrix-based 1-D FFT function.

int FFT(vector<complex>& vecSource, vector<complex>& vecFFTResult)
{
	// Get size of source vector
	int iSize = vecSource.GetSize();
	if(iSize < 1) return -5;
	
	// Create a complex matrix of size 1xn and copy source vector into matrix
	matrix<complex> matSource;
	matSource.SetSize(1,iSize);
	matSource.SetByVector(vecSource);
	
	// Create complex matrix to hold result of FFT operation
	matrix<complex> matFFTResult;
	
	// Call the matrix-based FFT function to perform FFT on this one row
	int iRet = FFT(matSource, matFFTResult);
	if(iRet) return iRet;

	// Extract FFT result from result matrix, into a complex matrix, and return
	matFFTResult.GetAsVector(vecFFTResult);
	
	// return 0 on success
	return 0;
}

// ER 01/29/03 QA70_3802	ADD_VECTOR_FFT_IFFT
// The following function performs 1-D inverse FFT on a complex vector.
// It uses the matrix-based 1-D inverse FFT function.

int IFFT(vector<complex>& vecSource, vector<complex>& vecFFTResult)
{
	// Get size of source vector
	int iSize = vecSource.GetSize();
	if(iSize < 1) return -5;
	
	// Create a complex matrix of size 1xn and copy source vector into matrix
	matrix<complex> matSource;
	matSource.SetSize(1,iSize);
	matSource.SetByVector(vecSource);
	
	// Create complex matrix to hold result of IFFT operation
	matrix<complex> matFFTResult;
	
	// Call the matrix-based IFFT function to perform IFFT on this one row
	int iRet = IFFT(matSource, matFFTResult);
	if(iRet) return iRet;

	// Extract IFFT result from result matrix, into a complex matrix, and return
	matFFTResult.GetAsVector(vecFFTResult);
	
	// return 0 on success
	return 0;
}

/**
	int SVD(matrix<complex>& matSource, matrix& matS, matrix<complex>& matU, matrix<complex>& matV);
	int SVD(matrix<complex>& matSource, matrix<complex>& matS); 
	int SVD(matrix& matSource, matrix& matS);
	int SVD(matrix& matSource, matrix& matS, matrix& matU, matrix& matV);

			These four functions are used to do the singular value decomposition. matSource = matU* matS* matV'
	Example:
	Parameters:
			matSource: The  m * n source data matrix
			matS: The min(m, n) diaganol sigular matrix 
			matU: The left side  m * m unitary matrix
			matV: The right side n* n unitary matrix
	return:
			if succeed, it will return 0

*/
	
int SVD(matrix& matSource, matrix & matS, matrix& matU, matrix& matV)
{
	//fist check the size of the matrix
	int nNumCols = matSource.GetNumCols();
	int nNumRows = matSource.GetNumRows();
	//no need to check size, because nag will take care for it
	//set the result matrix size
	matU.SetSize(nNumRows, nNumRows);
	matV.SetSize(nNumCols, nNumCols);
	vector vecE;
	int nMinSize = nNumCols > nNumRows ? nNumRows : nNumCols;
	vecE.SetSize(nMinSize -1);
	vector vecS;
	vecS.SetSize(nMinSize);
	//call the nag function
	int nRet;
	int nNumIter, nfailIfor;
	//we make a copy of source data here becasue we do not want to overwirte the data
	matrix matSourceCopy(matSource);
	if(nRet = nag_real_svd(nNumRows, nNumCols, matSourceCopy, nNumCols, 0, NULL, 0, TRUE, matU, nNumRows, 
							vecS, TRUE, matV, nNumCols, &nNumIter, vecE, &nfailIfor))
		return nRet;
	//set back the U V matrix
	if(nNumRows >= nNumCols)
		matU = matSourceCopy;
	else
		matV = matSourceCopy;
	matV.Transpose(); // Nag return MatV'
	//set back matS
	nRet = matS.SetDiagonal(vecS, 0);
	
	//here we do not consider nfailinfo, becasue it is rarely occurs.
	return nRet;				
				
}
	
int SVD(matrix<complex>& matSource, matrix& matS, matrix<complex>& matU, matrix<complex>& matV)
{
	//fist check the size of the matrix
	int nNumCols = matSource.GetNumCols();
	int nNumRows = matSource.GetNumRows();
	//no need to check size, because nag will take care for it
	//set the result matrix size
	matU.SetSize(nNumRows, nNumRows);
	matV.SetSize(nNumCols, nNumCols);
	vector vecE;
	int nMinSize = nNumCols > nNumRows ? nNumRows : nNumCols;
	vecE.SetSize(nMinSize -1 );
	vector vecS;
	vecS.SetSize(nMinSize);
	//call the nag function
	int nRet;
	int nNumIter, nfailIfor;
	//we make a copy of source data here becasue we do not want to overwirte the data
	matrix<complex> matSourceCopy(matSource); 
	if(nRet = nag_complex_svd(nNumRows, nNumCols, matSourceCopy, nNumCols, 0, NULL, 0, TRUE, matU, nNumRows, 
							vecS, TRUE, matV, nNumCols, &nNumIter, vecE, &nfailIfor))
		return nRet;
	//set back the U V matrix
	if(nNumRows >= nNumCols)
		matU = matSourceCopy;
	else
		matV = matSourceCopy;
	matV.Transpose(); // Nag return MatV'
	//set back matS
	nRet = matS.SetDiagonal(vecS, 0);

	//here we do not consider nfailinfo, becasue it is rarely occurs.
	return nRet;				
				
}
	
int SVD(matrix<complex>& matSource, matrix& matS) 
{
	matrix<complex> matU, matV;
	int nRet = 0;
	nRet = SVD(matSource, matS, matU, matV);
	return nRet;
}


int SVD(matrix& matSource, matrix& matS)
{
	matrix matU, matV;
	int nRet = 0;
	nRet = SVD(matSource, matS, matU, matV);
	return nRet;
}


/**
	int Trace(matrix<short> & matSource, int & nSum);
	int Trace(matrix<int> & matSource, int & nSum);
	int Trace(matrix<double> & matSource, double & dSum);
	int Trace(matrix<complex> & matSource, complex & cSum);

			These four functions will give the trace of a matrix
	Remarks:
			we do not support every matrix type.
 	Example:
	Parameters:
			matSource: The source matrix
			nSum:
			dSum:
			cSum:
				These are the variables to store the trace.
	Return:
			if succeed, it will return 0

*/

#define CHECK_MATRIX_FOR_TRACE {\
	int nNumCols = matSource.GetNumCols(); \
	int nNumRows = matSource.GetNumRows(); \
	if(!nNumCols || !nNumRows)  \
		return FFT_ERROR_SOURCE_MATRIX_EMPTY; \
}


int Trace(matrix<short> &matSource, int& nSum)
{
	CHECK_MATRIX_FOR_TRACE;

	vector<short> vecDiag;
	int nRet;
	if(nRet = matSource.GetDiagonal(vecDiag))
		return nRet;
	if(nRet = vecDiag.Sum(nSum))
		return nRet;

}

int Trace(matrix<int> &matSource, int& nSum)
{
	CHECK_MATRIX_FOR_TRACE;
	
	vector<int> vecDiag;
	int nRet;
	if(nRet = matSource.GetDiagonal(vecDiag))
		return nRet;
	if(nRet = vecDiag.Sum(nSum))
		return nRet;

}
int Trace(matrix<double> &matSource, double& dSum)
{
	CHECK_MATRIX_FOR_TRACE;
	
	vector<double> vecDiag;
	int nRet;
	if(nRet = matSource.GetDiagonal(vecDiag))
		return nRet;
	if(nRet = vecDiag.Sum(dSum))
		return nRet;

}	

int Trace(matrix<complex> & matSource, complex& cSum)
{
	CHECK_MATRIX_FOR_TRACE;

	vector<complex> vecDiag;
	int nRet;
	if(nRet = matSource.GetDiagonal(vecDiag))
		return nRet;
	if(nRet = vecDiag.Sum(cSum))
		return nRet;

}

///TCZ 07/29/02 QA70-2497	ADD_GLOBAL_STATS_SUMMARY_FUNCTION
//#define MAT_ERR_WEIGHTS_DATA_SIZE_NOT_MATCH (-5) //defined in the matrix.h
/**
			This function will give the summary of the statistics of a source data contained in a matrix	
	Remarks:
			If user does not have the weight, then just supply a matix without any initialization for weight 
	Example:
	Parameters:
			matData: The Data Source
			MatrixStatsSummary: The Results
			matWeights: The weight if any
	Return:
			Return 0 if succeed.
			MAT_ERR_WEIGHTS_DATA_SIZE_NOT_MATCH: The weight matrix size does not match that of the source
	int MatrixBasicStats(matrix & mData, MatrixStats & sMatStatsSum, matrix * mpWeights = NULL);

*/
int MatrixBasicStats(matrix & mData, MatrixStats & sMatStatsSum, matrix * mpWeights )
{
	
	//first check the size of the weight and data size 
	int nNumCols = mData.GetNumCols();
	int nNumRows = mData.GetNumRows();
	if(mpWeights != NULL) 
	{
		int nNumColsWeights = mpWeights->GetNumCols();
		int nNumRowsWeights = mpWeights->GetNumRows();
		if(nNumCols != nNumColsWeights || nNumRows != nNumRowsWeights)
			return MAT_ERR_WEIGHTS_DATA_SIZE_NOT_MATCH;
	}
	int nNumPoints = nNumCols * nNumRows;
	int nValid;
	double dMean, dSD, dMax, dMin, dWSum;
	int nRet;
	if(mpWeights == NULL)
	{
		if( nRet = nag_summary_stats_1var(nNumPoints, mData, NULL, &nValid, &dMean, &dSD, NULL, NULL, &dMin, &dMax, &dWSum))
			return nRet;
	}
	else
	{
		if( nRet = nag_summary_stats_1var(nNumPoints, mData, *mpWeights, &nValid, &dMean, &dSD, NULL, NULL, &dMin, &dMax, &dWSum))
			return nRet;
	
	}	
	//call the second nag function to get the median and hinge
	double dRes[5];
	if( nRet = nag_5pt_summary_stats(nNumPoints, mData, dRes))
		return nRet;
	//set the structure value
	sMatStatsSum.nMissingValue = nNumPoints - nValid;
	sMatStatsSum.dMean = dMean;
	sMatStatsSum.dSD = dSD;
	sMatStatsSum.dMax = dMax;
	sMatStatsSum.dMin = dMin;
	sMatStatsSum.dMedian = dRes[2];
	sMatStatsSum.dLowHinge = dRes[1];
	sMatStatsSum.dUpHinge = dRes[3];
	sMatStatsSum.dWSum = dWSum;
	
	return nRet;	
}
//END-------------------------------ADD_GLOBAL_STATS_SUMMARY_FUNCTION