/*------------------------------------------------------------------------------*
 * File Name: stats_utils.c 													*
 * Creation:																	*
 * Purpose: OriginC source file	containing LLOCs for statistics tools.			*
 * Copyright (c) OriginLab Corp.	2003, 2004, 2005, 2006, 2007, 2008, 2009	*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *	CPY 7/17/04 FIT_PARAM_ARRAY_NEED_SIZE										*
 *	Echo 7/21/04 MR_AVOID_ERROR_SETTING											*
 *	Max 9/09/04 MODIFY_FIXINTERCEPT												*
 *	ML 9/27/2004 USE_ONAG8_OC													*
 *	Max 10/13/04 REPLACE_KOLMOGOROV_SMIRNOV_TEST_WITH_NAG_FUNTION				*
 *  Echo 10/27/04 COMMENT_FOR_USING_NAG_ENUM									*
 *	Max 10/27/04 MODIFY_FIXINTERCEPTAT_CASE										*	
 *  Sonic 12/22/04	COMMENT_TO_FIX_BUG											*	
 *	Max 03/16/05 REWRITE_ONEWAYANOVA_MEAN_COMPARISON							*
 *	Max 03/16/05 REWRITE_TWOWAYANOVA_MEAN_COMPARISON							*
 *	Leo 08/02/05 TO_CENTRALIZE_POLYNOMIAL_FIT									*
 *	Max 08/24/05 MODIFY_DOF_OF_SHAPIRO_WILK_TEST								*
 *  Echo 9/12/05 MADE_MR_FASTER													*
 *	Leo 2006-4-7 FIX_NULL_REFERENCE_BUG											*
 *	Echo 9/12/06 NAG_FUNC_MOVED_TO_LLVC											*
 *	Arvin 01/25/08 QA70-10923 CALC_CONFIDENCE_AND_PREDICTION_BANDS_OF_POLYNOMIAL_FIT
 *	Fisher 12/08/08 QA80-12735  SUPPORT_USER_SUPPLIED_BIN_VECTOR_IN_FREQCOUNTS	*  
 *	Kyle 03/25/2009 QA80-13251-P4 SKIP_OR_EXTEND_THE_LAST_BIN_IF_IT_IS_NOT_A_FULL_BIN
 *	Hong 09/09/09 QA80-14287 FREQUENCY_COUNT_ADD_MIN_MAX_COUNT					*
 *------------------------------------------------------------------------------*/

#include <Origin.h>
//#include <OC_nag8.h>	//CPY 10/13/04 ///Echo 9/12/06 NAG_FUNC_MOVED_TO_LLVC
#include <stats_utils.h>		// Prototypes for LLOCs



///--- CPY 9/23/04
///Echo 10/27/04 COMMENT_FOR_USING_NAG_ENUM
//#define Nag_ClassBoundary int// need to replace with NAG enum
//#define Nag_ClassBoundaryUser 0// need to replace with NAG enum
//#define NE_NOERROR 0// need to replace with NAG enum
///end COMMENT_FOR_USING_NAG_ENUM
// begin ONAG calling old way
///////////////////////////////////////////////////////////////////////////////
///Echo 9/12/06 NAG_FUNC_MOVED_TO_LLVC
//// NAG functions that we still need to call in the old Origin7 way
//int onag_shapiro_wilk_test(
	//int n,	// the sample size, 3 <= n <= 2000
	//double x[],  // the ordered sample values
	//bool calc_wts, // calculate weight or not
	//double a[],	// weight
	//double* w,	// the value of the statistics, W
	//double* pw	// the significance level of W
//)	// Shapiro and Wilk's W test for Normality
//{
	//NagError		stNagErr;
	//nag_shapiro_wilk_test(n, x, calc_wts, a, w, pw, &stNagErr);
	//return stNagErr.code;
//}
//
//int onag_frequency_table(
	//int n, // the number of observations
	//const double x[], // the sample of observations, x[i]. i = 0,1,2,...,n 
	//int num_class,	// the number of classes desired in the frequency distribution
	//Nag_ClassBoundary class_val,
	//double cint[],
	//int ifreq[],
	//double *xmin,	// the smallest value in the sample
	//double *xmax	// the largest value in the sample
//)	// Frequency table from raw data
//{
	//NagError		stNagErr;
	//nag_frequency_table(n, x, num_class, class_val, cint, ifreq, xmin, xmax, &stNagErr);
	//return stNagErr.code;
//}

//// end ONAG calling old way

///END ///Echo 9/12/06 NAG_FUNC_MOVED_TO_LLVC


///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////


/////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//				Basic Functions: Descriptive Statistics, Percentile, Range, Quantile, Frenqucy Count									
//
/////////////////////////////////////////////////////////////////////////////////////////////////////////

/**
		Function to compute descriptive statistics (moments) with weighting passing results
*/

int stats_descriptive(const vector& vData, DescStatResults& dsr,
	 const vector* pvWeight, double dConfLevel, int iDenFlag) // dConfLevel = SU_DEFAULT_CONF_LEVEL, pvWeight = NULL
{
	//ocmathh_desc_stats will take care of missing values by making local copies of both data of weight vector
	if(pvWeight)		
	{
    	if ( pvWeight->GetSize() == vData.GetSize()) 		
			return ocmath_desc_stats(vData, vData.GetSize(), &dsr, *pvWeight, dConfLevel,iDenFlag);
    	else //if size of data is different from size of weight, return error.
    		return STATS_ERROR_WEIGHT_DIFF_SIZE;
	}
	else
	{
		return ocmath_desc_stats(vData, vData.GetSize(), &dsr, NULL, dConfLevel, iDenFlag);// when no weight is provided, each data value is assigned a weight of 1
	}
}

#define CHK_SET_DESC_STATS(_nn, _var) if (opt._nn) _var = &dsr._nn;

static void _init_desc_stat_res(DescStatResults& dsr)
{
	dsr.N = -1;
	dsr.Missing = -1;
	dsr.Sum = NANUM;
	dsr.Mode = NANUM;
	dsr.Mean = NANUM;
	dsr.GeoMean = NANUM;
	dsr.GeoSD = NANUM;
	dsr.LCL = NANUM;
	dsr.UCL = NANUM;
	dsr.SD = NANUM;
	dsr.SDx2 = NANUM;
	dsr.SDx3 = NANUM;
	dsr.SEM = NANUM;
	dsr.Variance = NANUM;
	dsr.COV = NANUM;
	dsr.Skewness = NANUM;
	dsr.Kurtosis = NANUM;
	dsr.MAD = NANUM;
	dsr.UncorSumSq = NANUM;
	dsr.CorSumSq = NANUM;
	dsr.WeightSum = NANUM;
	
}

int stats_desc_stats(const vector& vData, const DescStatOptions& opt, DescStatResults& dsr, 
					const vector* pvWeight, double dConfLevel, bool* pmultimodes) // dConfLevel = SU_DEFAULT_CONF_LEVEL, pvWeight = NULL, pmultimodes = NULL
{
	double *pSum, *pMean, *pMode, *pGeoMean, *pGSD, *pLCL, *pUCL, *pSD, *pSDx2, *pSDx3, *pSE, *pVar, *pCOV, *pSkew, *pKurt, *pMAD, *pUSS, *pCSS, *pWeightSum;
	int *pN, *pMissing;
	
	int nVarDivisor = opt.VarDivisor;
	int nSize = vData.GetSize();
	
	_init_desc_stat_res(dsr);
	
	CHK_SET_DESC_STATS(N, pN)
	CHK_SET_DESC_STATS(Missing, pMissing)
	CHK_SET_DESC_STATS(Sum, pSum)
	CHK_SET_DESC_STATS(Mode, pMode)
	CHK_SET_DESC_STATS(Mean, pMean)
	CHK_SET_DESC_STATS(GeoMean, pGeoMean)
	CHK_SET_DESC_STATS(GeoSD, pGSD)
	CHK_SET_DESC_STATS(LCL, pLCL)
	CHK_SET_DESC_STATS(UCL, pUCL)
	CHK_SET_DESC_STATS(SD, pSD)
	CHK_SET_DESC_STATS(SDx2, pSDx2)
	CHK_SET_DESC_STATS(SDx3, pSDx3)
	CHK_SET_DESC_STATS(SEM, pSE)
	CHK_SET_DESC_STATS(Variance, pVar)
	CHK_SET_DESC_STATS(COV, pCOV)
	CHK_SET_DESC_STATS(Skewness, pSkew)
	CHK_SET_DESC_STATS(Kurtosis, pKurt)
	CHK_SET_DESC_STATS(MAD, pMAD)
	CHK_SET_DESC_STATS(UncorSumSq, pUSS)
	CHK_SET_DESC_STATS(CorSumSq, pCSS)
	CHK_SET_DESC_STATS(WeightSum, pWeightSum)

	int nRet = ocmath_basic_summary_stats(nSize, vData, pN, pMean, pSD, pSE, pVar, pSum, pSkew, pKurt, pUSS, pCSS, pCOV, pMAD, pMissing, pWeightSum, pGeoMean, nVarDivisor, *pvWeight);
	if ( STATS_NO_ERROR != nRet)	
		return nRet;
	
	return nRet = ocmath_opt_summary_stats(nSize, vData, dConfLevel, pMode, pmultimodes, pGSD, pLCL, pUCL, pSDx2, pSDx3, nVarDivisor, *pvWeight);		
}


/**
		Function to compute quantiles including Min, iMin, Max, iMax, Range and IQR for a vector.
*/
int stats_quantiles(const vector& vData, const QuantileOptions& opt, QuantileResults& res,
	 const vector& vPercents, vector& vPercentiles, const vector* pvWeight) // pvPercents = NULL, pvPercentiles = NULL, pvWeight = NULL;
{	
	int nSize = vData.GetSize();
	//if (nSize <= 0)
		//return STATS_ERROR_SETTING;

	int nPercentSize = 0;
	if (vPercents)
		nPercentSize = vPercents.GetSize();

	double *pvPercentiles = NULL;
	double *pvPercents = NULL;
	if(nPercentSize)
	{
		vPercentiles.SetSize(nPercentSize);
		pvPercentiles = vPercentiles;
		pvPercents = vPercents;
	}
	return ocmath_quantiles(vData, nSize, &opt, &res, pvPercents, nPercentSize, pvPercentiles, *pvWeight);
}



/**
		Function to compute percentiles for a vector passing percents and percentiles in vectors.
*/
int stats_percentiles(const vector& vData, const vector& vPercents, vector& vPercentiles, int nInterpolate, const vector* pvWeight) // nInterpolate = INTERPOLATE_EDF
{
	int nSize = vData.GetSize();
	int nPercentSize = vPercents.GetSize();
	vPercentiles.SetSize(nPercentSize); //YuI 12/02/04
	if (pvWeight)
		return ocmath_percentiles(vData, nSize, vPercents, nPercentSize, vPercentiles, nInterpolate, *pvWeight);
	else 
		return ocmath_percentiles(vData, nSize, vPercents, nPercentSize, vPercentiles, nInterpolate);		
}



/**
		Function to compute range statistics for a vector passing results in a structure.
*/
int stats_range(const vector& vData, RangeResults& rr)
{
	return ocmath_range(vData, vData.GetSize(), &rr);

}



/**
		Function to perform frequency count
*/
/*
int stats_frequency_count(vector& vDataIn, FreqCountOptions& fcoOptions, vector& vBinCenters,
   vector& vAbsoluteCounts, vector* vCumulativeCounts )
{
	
	return ocmath_frequency_count(vDataIn, vDataIn.GetSize(), &fcoOptions, 
								vBinCenters, vBinCenters.GetSize(),
								vAbsoluteCounts, vAbsoluteCounts.GetSize(),
								*vCumulativeCounts, vCumulativeCounts->GetSize());
	
}
*/

/// YuI 11/02/04 V7.5158 QA70-7122 OC_DATA_PLOT_CREATION
// this function has too many mistakes
// not checking vCumulativeCounts pointer for NULL 
// assign double to int etc.
// I am rewriting it
/*
int stats_frequency_count(vector& vDataIn, FreqCountOptions& fcoOptions, vector& vBinCenters,
   vector& vAbsoluteCounts, vector* vCumulativeCounts, int nOption) 
{
	
	int ii;
	vector vBoundary;
	int nSourceSize = vDataIn.GetSize();
	int nBinCenterSize = vBinCenters.GetSize();
	int nAbsoluteCountSize = vAbsoluteCounts.GetSize();
	int nCumulativeCountSize = vCumulativeCounts->GetSize();
	
	for (ii = 0; ii < nBinCenterSize; ii++)
		vBinCenters[ii] = NANUM;
	for (ii = 0; ii < nAbsoluteCountSize; ii++)
		vAbsoluteCounts[ii] = NANUM;
	for (ii = 0; ii < nCumulativeCountSize; ii++)
		vCumulativeCounts[0][ii] = NANUM;		

	//Check size
	if(2 > nSourceSize)
		return STATS_ERROR_TOO_FEW_DATA_PTS;

	double	dMin, dMax, dInc;
	bool	bIncludeLTMin = true, bIncludeGEMax = true;
	int 	nIntervals;
	dMin = fcoOptions.FromMin;
	dMax = fcoOptions.ToMax;
	bIncludeLTMin = fcoOptions.IncludeLTMin;
	bIncludeGEMax = fcoOptions.IncludeGEMax;
	if (nOption == FC_STEPSIZE )
		dInc = fcoOptions.StepSize;
	else
		nIntervals = fcoOptions.StepSize;

	if (nOption == FC_STEPSIZE )
		nIntervals = (int)((dMax - dMin)/dInc);
	else
		dInc = (dMax - dMin)/nIntervals;
	int nBoundSize = nIntervals + 2;
	vBoundary.SetSize(nBoundSize-1);
	if(bIncludeLTMin) nIntervals++;
	if(bIncludeGEMax) nIntervals++;

	int		i;
	
	//Calculate Bin Center
	if( nBinCenterSize != nIntervals )
		return STATS_ERROR_INPUT_SIZE;

	double	dLowVal = dMin;
	for(i = 0; i < nBinCenterSize; i++)
	{
		if(bIncludeLTMin && i==0)  //if include LTMin, the center of first bin should be 0.
			vBinCenters[0] = dMin;  
		else
		{
			vBinCenters[i] = dLowVal + dInc/2;  
			dLowVal += dInc; 
		}
	}
	dLowVal = vBoundary[0] = dMin;
	for (i = 0; i < nBoundSize-2; i++)
	{
		vBoundary[i + 1] = dLowVal + dInc;
		dLowVal += dInc;
	}
	
	//Calculate absolute count
	if( nAbsoluteCountSize != nIntervals )
		return STATS_ERROR_INPUT_SIZE; 
	vector<int> vFreq;
	vFreq.SetSize(nBoundSize);
	double dXMin, dXMax;
	// From <NAG\OCN_g01.h>: g01aec nag_frequency_table: Conducts a frequency distribution of a variable
	int iErr = onag_frequency_table(nSourceSize, vDataIn, nBoundSize, Nag_ClassBoundaryUser, vBoundary, vFreq, &dXMin, &dXMax);
	if( iErr == NE_NOERROR )
	{
		//NAG Success so continue... 
		if (!bIncludeLTMin)
			vFreq.RemoveAt(0);
		int nTempSize = vFreq.GetSize();
		if (!bIncludeGEMax)
			vFreq.RemoveAt(nTempSize - 1);	
		vAbsoluteCounts = vFreq;
	}else
		return iErr;/// else stop and return error if NAG function returns error
	
	//Calculate Cumulative Counts	
	if( nCumulativeCountSize != nIntervals )
		return STATS_ERROR_INPUT_SIZE;

	vCumulativeCounts[0][0] = vAbsoluteCounts[0]; 
	for( i = 1; i< nCumulativeCountSize; i++ )
	{
		vCumulativeCounts[0][i] = vCumulativeCounts[0][i-1] + vAbsoluteCounts[i];
	}

	return STATS_NO_ERROR;
	
}
*/
///Echo 1/10/05 USE_LLVC_FOR_CENTRILAZING

//int stats_frequency_count(const vector& vDataIn, const FreqCountOptions& fcoOptions, vector& vBinCenters,
	   //vector& vAbsoluteCounts, vector& vCumulativeCounts, int nOption)
//{
	////Check size
	//if( vDataIn.GetSize() < 2 )
		//return STATS_ERROR_TOO_FEW_DATA_PTS;
	//
	//for( int ii = 0; ii < vBinCenters.GetSize(); ii++ )
		//vBinCenters[ii] = NANUM;
	//
	//for( ii = 0; ii < vAbsoluteCounts.GetSize(); ii++ )
		//vAbsoluteCounts[ii] = NANUM;
	//
	//if( NULL != vCumulativeCounts )
	//{
		//for( ii = 0; ii < vCumulativeCounts.GetSize(); ii++)
			//vCumulativeCounts[ii] = NANUM;		
	//}
	//
	//int nIntervals = (nOption == FC_STEPSIZE ) ? ((fcoOptions.ToMax - fcoOptions.FromMin) / fcoOptions.StepSize) : (fcoOptions.StepSize);
	//if (nIntervals < 1)
		//nIntervals = 1;
	//
	//double dInc = (nOption == FC_STEPSIZE ) ? (fcoOptions.StepSize) : ((fcoOptions.ToMax - fcoOptions.FromMin) / nIntervals);
	//
	//vector vBoundaries;
	//vBoundaries.SetSize(nIntervals + 1);
	//
	//vector<int> vFreq;
	//vFreq.SetSize(nIntervals + 2);
	//
	//if( fcoOptions.IncludeLTMin ) 
		//nIntervals++;
	//if( fcoOptions.IncludeGEMax )
		//nIntervals++;
//
	//vBinCenters.SetSize(nIntervals);
	//for( ii = 0; ii < nIntervals; ii++ )
	//{
		//if( ii == 0 )
		//{
			//vBinCenters[ii] = fcoOptions.IncludeLTMin ? fcoOptions.FromMin : fcoOptions.FromMin + dInc / 2;
		//}
		//else
		//{
			//vBinCenters[ii] = fcoOptions.FromMin + dInc * ((double)ii + 0.5);
		//}
	//}
//
	//for( ii = 0; ii < vBoundaries.GetSize(); ii++ )
		//vBoundaries[ii] = fcoOptions.FromMin + ii * dInc;
//
	//double dXMin, dXMax;
	//// From <NAG\OCN_g01.h>: g01aec nag_frequency_table: Conducts a frequency distribution of a variable
	//int iErr = onag_frequency_table(vDataIn.GetSize(), vDataIn, vFreq.GetSize(), Nag_ClassBoundaryUser, vBoundaries, vFreq, &dXMin, &dXMax);
	//if( iErr == NE_NOERROR )
	//{
		//if( !fcoOptions.IncludeLTMin )
			//vFreq.RemoveAt(0);
		//if( !fcoOptions.IncludeGEMax )
			//vFreq.RemoveAt(vFreq.GetSize() - 1);	
		//
		//vAbsoluteCounts = vFreq;
	//}
	//else
		//return iErr;/// else stop and return error if NAG function returns error
	//
	//if( vCumulativeCounts )
	//{
		//vCumulativeCounts.SetSize(vAbsoluteCounts.GetSize());
		//for( ii = 0; ii < vCumulativeCounts.GetSize(); ii++ )
		//{
			//if( 0 == ii )
				//vCumulativeCounts[ii] = vAbsoluteCounts[ii];
			//else
				//vCumulativeCounts[ii] = vCumulativeCounts[ii - 1] + vAbsoluteCounts[ii];
		//}
	//}
//
	//return STATS_NO_ERROR;
//}
/// end OC_DATA_PLOT_CREATION
/// Fisher 12/08/08         QA80-12735       SUPPORT_USER_SUPPLIED_BIN_VECTOR_IN_FREQCOUNTS 
int stats_frequency_count(vector& vDataIn, FreqCountOptions& fcoOptions, vector& vBinCenters,
   vector& vAbsoluteCounts, vector& vCumulativeCounts, int nOption, vector& vBins)
//int stats_frequency_count(vector& vDataIn, FreqCountOptions& fcoOptions, vector& vBinCenters,
// vector& vAbsoluteCounts, vector& vCumulativeCounts, int nOption)
/// End SUPPORT_USER_SUPPLIED_BIN_VECTOR_IN_FREQCOUNTS
{
	/// Fisher 12/08/08 QA80-12735 SUPPORT_USER_SUPPLIED_BIN_VECTOR_IN_FREQCOUNTS 
	int nIntervals;
	bool bUseUserBins;
	if(!vBins || vBins.GetSize() <= 1)
	{
		///Kyle 03/25/2009 QA80-13251-P4 SKIP_OR_EXTEND_THE_LAST_BIN_IF_IT_IS_NOT_A_FULL_BIN
		//nIntervals = (nOption == FC_STEPSIZE ) ? ceil((fcoOptions.ToMax - fcoOptions.FromMin) / fcoOptions.StepSize) : (fcoOptions.StepSize);
		nIntervals = (nOption == FC_STEPSIZE ) ? (int)((fcoOptions.ToMax - fcoOptions.FromMin) / fcoOptions.StepSize + 0.5) : (fcoOptions.StepSize);
		///End SKIP_OR_EXTEND_THE_LAST_BIN_IF_IT_IS_NOT_A_FULL_BIN
		if (nIntervals < 1)
			nIntervals = 1;
	
		if( fcoOptions.IncludeLTMin )	nIntervals++;
		if( fcoOptions.IncludeGEMax )	nIntervals++;
		/// Hong 09/09/09 QA80-14287 FREQUENCY_COUNT_ADD_MIN_MAX_COUNT
		if( fcoOptions.IncludeMin )	nIntervals++;
		if( fcoOptions.IncludeMax )	nIntervals++;
		/// end FREQUENCY_COUNT_ADD_MIN_MAX_COUNT
		
		bUseUserBins = false;
	}
	else
	{
		nIntervals = vBins.GetSize()-1;
		bUseUserBins = true;
	}
	
	vBinCenters.SetSize(nIntervals);	
	vAbsoluteCounts.SetSize(nIntervals);
	
	int nCumCounts; 
	
	double* pvCumulativeCounts = NULL;
	if (vCumulativeCounts)
	{
		nCumCounts = nIntervals; 
		vCumulativeCounts.SetSize(nCumCounts);
		pvCumulativeCounts = vCumulativeCounts;
	}
	
	if(bUseUserBins)
		return ocmath_frequency_count(vDataIn, vDataIn.GetSize(), &fcoOptions, 
								vBinCenters, vBinCenters.GetSize(),
								vAbsoluteCounts, vAbsoluteCounts.GetSize(),
								pvCumulativeCounts, nCumCounts, nOption, vBins, vBins.GetSize());	
	else
		return ocmath_frequency_count(vDataIn, vDataIn.GetSize(), &fcoOptions, 
								vBinCenters, vBinCenters.GetSize(),
								vAbsoluteCounts, vAbsoluteCounts.GetSize(),
								pvCumulativeCounts, nCumCounts, nOption);
	//int nIntervals = (nOption == FC_STEPSIZE ) ? ceil((fcoOptions.ToMax - fcoOptions.FromMin) / fcoOptions.StepSize) : (fcoOptions.StepSize);
	//if (nIntervals < 1)
		//nIntervals = 1;
	//
	//if( fcoOptions.IncludeLTMin )	nIntervals++;
	//if( fcoOptions.IncludeGEMax )	nIntervals++;
	//
	//vBinCenters.SetSize(nIntervals);	
	//vAbsoluteCounts.SetSize(nIntervals);
	//
	//int nCumCounts; 
	//
	//double* pvCumulativeCounts = NULL;
	//if (vCumulativeCounts)
	//{
		//nCumCounts = nIntervals; 
		//vCumulativeCounts.SetSize(nCumCounts);
		//pvCumulativeCounts = vCumulativeCounts;
	//}
	//
	//return ocmath_frequency_count(vDataIn, vDataIn.GetSize(), &fcoOptions, 
								//vBinCenters, vBinCenters.GetSize(),
								//vAbsoluteCounts, vAbsoluteCounts.GetSize(),
								//pvCumulativeCounts, nCumCounts, nOption);	
	/// End SUPPORT_USER_SUPPLIED_BIN_VECTOR_IN_FREQCOUNTS
}
/// END USE_LLVC_FOR_CENTRILAZING
/////////////////////////////////////////////////////////////////////////////////
//
//				Regression Functions									
//
/////////////////////////////////////////////////////////////////////////////////
/**
		General Linear Regression 
*/
static int stats_linear_regrssion_sum(const vectorbase& vX, const vectorbase& vY, const vectorbase& vW,
	double& dSumWeight, double& dxAve, double& dyAve, double& dSxx, double& dSxy, double& dSyy)
{
	vector vtemp;	

	vW.Sum(dSumWeight);

	// Average value of y
	vtemp = vY * vW;
	vtemp.Sum(dyAve);
	dyAve /= dSumWeight;

	// Syy
	vtemp *= vY;
	vtemp.Sum(dSyy);
	dSyy -= dyAve^2 * dSumWeight;

	// Average value of x
	vtemp = vX * vW;
	vtemp.Sum(dxAve);
	dxAve /= dSumWeight;

	// Sxx
	vtemp *= vX;
	vtemp.Sum(dSxx);
	dSxx -= dxAve^2 * dSumWeight;
	if( dSxx < 1e-15 )
		return STATS_ERROR_REGRESSION_X_RANGE;

	// Sxy
	vtemp = vW * vY * vX;
	vtemp.Sum(dSxy);
	dSxy -= dxAve * dSumWeight * dyAve;

	return STATS_NO_ERROR;
}
int stats_calc_internal(const vector& vX, const vector& vY, RegCalcInternal& sInternal, const vector& vWeight)
{
	vector vWTmp;
	vWTmp.SetSize(vX.GetSize());
	if (vWeight)
	{
		if (vWeight.GetSize() != vX.GetSize())
			return STATS_ERROR_WEIGHT_DIFF_SIZE;
		vWTmp = vWeight;
	}
	else
	{
		vWTmp = 1;
	}
	if (vX.GetSize() != vY.GetSize())
		return STATS_ERROR_DATA_DIFF_SIZE;
		
	int nRet = stats_linear_regrssion_sum(vX, vY, vWTmp, sInternal.dSumWeight, sInternal.dxAve, sInternal.dyAve, sInternal.dSxx, sInternal.dSxy, sInternal.dSyy);
	if( STATS_NO_ERROR != nRet )
		return nRet;
	
	return STATS_NO_ERROR;
		
}

/**
		Linear Regression Parameter Initialization
*/
///Echo 11/18/04
//For sake of centralize the code, I have rewriten stats_linear_fit.
//static void RegStatsNANUM(RegStats* psRegStats)
//{  
	//psRegStats->N = NANUM;
	//psRegStats->DOF = NANUM;
	//psRegStats->ReducedChiSq = NANUM;
	//psRegStats->SSR = NANUM;
	//psRegStats->Correlation = NANUM;
	//psRegStats->Rvalue = NANUM;
	//psRegStats->RSqCOD = NANUM;
	//psRegStats->AdjRSq = NANUM;
	//psRegStats->RMSESD = NANUM;
	//psRegStats->NormResiduals = NANUM;
//}
//static void RegANOVANANUM(RegANOVA* psRegANOVA)
//{
	//psRegANOVA->AnovaTable[ANOVA_ERR].DOF = NANUM;
	//psRegANOVA->AnovaTable[ANOVA_ERR].SSq = NANUM;
	//psRegANOVA->AnovaTable[ANOVA_ERR].MeanSq = NANUM;
	//
	//psRegANOVA->AnovaTable[ANOVA_TOTAL].DOF = NANUM;
	//psRegANOVA->AnovaTable[ANOVA_TOTAL].SSq = NANUM;
	//psRegANOVA->AnovaTable[ANOVA_TOTAL].MeanSq = NANUM;
	//
	//psRegANOVA->AnovaTable[ANOVA_MODEL].DOF = NANUM;
	//psRegANOVA->AnovaTable[ANOVA_MODEL].SSq = NANUM;
	//psRegANOVA->AnovaTable[ANOVA_MODEL].MeanSq = NANUM;
	//
	//psRegANOVA->Fvalue = NANUM;
	//psRegANOVA->Pvalue = NANUM;
//}

//void FitParameterNANUM(FitParameter& psFitParameter) //------ CPY 7/30/04 better use pointer, as ref in OC may not work when passing as psFitParameter[1]
//static void FitParameterNANUM(FitParameter* psFitParameter)
//{
	////psFitParameter.Value = NANUM;
	//psFitParameter->Value = NANUM;
	//psFitParameter->Error = NANUM; 
	//psFitParameter->tValue = NANUM;
	//psFitParameter->LCL = NANUM;
	//psFitParameter->UCL = NANUM;
	//psFitParameter->Prob = NANUM;
	//psFitParameter->Dependency = NANUM;
	//psFitParameter->ConfInterval = NANUM;
//}

/**
		General Linear Regression with vector inputs.
*/
int stats_linear_fit(const vector& vX, const vector& vY, const vector& vWeightIn,
const LROptions& sLROptions, 
FitParameter* psFitParameter, // [0] = intercept and [1] = slope
RegStats* psRegStats, RegANOVA* psRegANOVA, RegCalcInternal* psInternal, // = NULL,  = NULL = NULL
matrix& mCov, matrix& mCorr)
{
	if( sLROptions.FixIntercept && sLROptions.FixSlope )
		return STATS_ERROR_SETTING;
	
	if (0 == vWeightIn.GetSize())
	{
		vWeightIn.SetSize(vX.GetSize());
		for (int ii = 0; ii < vX.GetSize(); ii++)
			vWeightIn[ii] = 1;
	}
	
	if (vX.GetSize() != vWeightIn.GetSize())
			return STATS_ERROR_WEIGHT_DIFF_SIZE;
	
	if(vX.GetSize() != vY.GetSize()) 
			return STATS_ERROR_DATA_DIFF_SIZE;

	int nCovCorrRow, nCovCorrCol;

	if (mCov || mCorr)
		///Arvin 01/25/08 QA70-10923 CALC_CONFIDENCE_AND_PREDICTION_BANDS_OF_POLYNOMIAL_FIT
		//We should keep covariance and correlation matrix consistent with NLSF
		//If parameter is fixed, we should set corresponding row and col as missing values 
		//nCovCorrRow = nCovCorrCol = (sLROptions.FixIntercept || sLROptions.FixSlope ) ? 1 : 2;
		nCovCorrRow = nCovCorrCol = 2;
		///end CALC_CONFIDENCE_AND_PREDICTION_BANDS_OF_POLYNOMIAL_FIT

	//if( !sLROptions.FixIntercept && !sLROptions.FixSlope ) 
		//nCovCorrRow = nCovCorrCol = 2;
	//else if (sLROptions.FixIntercept && sLROptions.FixSlope)
			//nCovCorrRow = nCovCorrCol = 0;
	//else
		//nCovCorrRow = nCovCorrCol = 1;	
	double* pmCov = NULL;
	double* pmCorr = NULL;
	if (mCov)
	{
		mCov.SetSize(nCovCorrRow, nCovCorrCol);
		pmCov = mCov;
	}
	if(mCorr)
	{
		mCorr.SetSize(nCovCorrRow, nCovCorrCol);
		pmCorr = mCorr;
	}	
	
	vector vWeight;          // Weight function for calculation
	///Echo 11/23/04 
	//Weight would be process before entering the function
	// Check the weight function (if needed)
	//if (sLROptions.ErrBarWeight)
		//vWeight = 1 / vWeightIn^2;
	//else
		vWeight = vWeightIn;
		
	double da_Err, db_Err;		// The standard error of the intercept and slope
	double rsq, rss, dof;		// The R-square, residual sum of squares, and the degrees of freedom
	int mdof;					// The degrees of freedom of the Model

	double dSumWeight;			// The sum of weights
	double dSxx, dSxy, dSyy;
	double dxAve, dyAve;		// The average value of x and y
	
	///Echo 12/9/04 PSINTERNAL_COMPUTED_IN_LLVC
	//if( 0 != stats_linear_regrssion_sum(vX, vY, vWeight, dSumWeight, dxAve, dyAve, dSxx, dSxy, dSyy) )
		//return STATS_ERROR_REGRESSION_X_RANGE;			

	//if(psInternal)
	//{
		//psInternal->dSumWeight = dSumWeight;
		//psInternal->dSxx = dSxx;
		//psInternal->dSxy = dSxy;
		//psInternal->dSyy = dSyy;
		//psInternal->dxAve = dxAve;
		//psInternal->dyAve = dyAve;
	//}	
			
	return 	ocmath_linear_fit(vX, vY, vX.GetSize(),psFitParameter, vWeight, vWeight.GetSize(), &sLROptions, psRegStats, psRegANOVA, psInternal, pmCov, pmCorr, nCovCorrRow, nCovCorrCol);

}
///Echo 11/18/04
//For sake of centralize the code, I have rewriten stats_linear_fit.
/*
int stats_linear_fit(const vector& vX, const vector& vY, const vector& vWeightIn,
const LROptions& sLROptions,
FitParameter* psFitParameter, // [0] = intercept and [1] = slope
RegStats* psRegStats, RegANOVA* psRegANOVA, RegCalcInternal* psInternal)// = NULL,  = NULL = NULL
{
	if(psRegStats)
	{
		RegStatsNANUM(psRegStats);
	}
	if(psRegANOVA)
	{
		RegANOVANANUM(psRegANOVA);
	}		
	FitParameterNANUM(&psFitParameter[INTERCEPT]);
	FitParameterNANUM(&psFitParameter[SLOPE]);
	///Max 9/09/04 MODIFY_FIXINTERCEPT
	//if( sLROptions.ThroughZero && sLROptions.FixSlope )
	if( sLROptions.FixIntercept && sLROptions.FixSlope )
	///End MODIFY_FIXINTERCEPT
		return STATS_ERROR_SETTING;
	if( sLROptions.Confidence < 0 || sLROptions.Confidence > 1)
	{
		return STATS_ERROR_SETTING;
	}
	double da_Err, db_Err;		// The standard error of the intercept and slope
	double rsq, rss, dof;		// The R-square, residual sum of squares, and the degrees of freedom
	int mdof;					// The degrees of freedom of the Model

	double dSlope, dIntercept;
	double dSumWeight;			// The sum of weights
	double dSxx, dSxy, dSyy;
	double dxAve, dyAve;		// The average value of x and y

	double SD;
	
	vector vWeight;          // Weight function for calculation
	
	vector vANOVA(MULREG_Size);
	// Get the number of points of the curve
	int nPts = vY.GetSize();
	if( nPts < SU_MIN_NPTS || nPts != vX.GetSize() || nPts != vWeightIn.GetSize()) 
		if ( nPts < SU_MIN_NPTS )
			return STATS_ERROR_TOO_FEW_DATA_PTS;
		else if ( nPts != vX.GetSize() )
			return STATS_ERROR_DATA_DIFF_SIZE;
			else
				return STATS_ERROR_WEIGHT_DIFF_SIZE;
	// Check the weight function (if needed)
	if (sLROptions.ErrBarWeight)
		vWeight = 1 / vWeightIn^2;
	else
		vWeight = vWeightIn;

	
	// Linear Regression begin
	
	if( sLROptions.FixSlope ) // Fix slope
	{
		psFitParameter[SLOPE].Fix=1;
		psFitParameter[INTERCEPT].Fix=0;

		if( 0 != stats_linear_regrssion_sum(vX, vY, vWeight, dSumWeight, dxAve, dyAve, dSxx, dSxy, dSyy) )
			return STATS_ERROR_REGRESSION_X_RANGE;
			
			
		// Slope
		dSlope = sLROptions.FixSlopeAt;
		db_Err = 0;

		// Intercept		
		dIntercept = dyAve - dSlope * dxAve;

		// Residual sum of squares
		rss = dSyy - 2 * dSlope * dSxy + dSlope^2 * dSxx;

		// Degrees of freedom
		dof = nPts - 1;
		mdof = 0;

		// Standard error of intercept
		da_Err = sqrt(rss / dof / dSumWeight);

		// R-square
		rsq = dSxy^2 / dSxx / dSyy;
		
		//Standard deviation
		SD = sqrt(rss / dof);
		
	}
	else // Slope is not fixed
	{
		psFitParameter[SLOPE].Fix = 0 ;		
		vector vCoeff;
		matrix mCov;
		///Max 9/09/04 MODIFY_FIXINTERCEPT
		//if( sLROptions.ThroughZero == 0 ) // With constant
		if( sLROptions.FixIntercept == 0 ) // With changeable Intercept
		///End MODIFY_FIXINTERCEPT
		{	
			mCov.SetSize(2,2);
			vCoeff.SetSize(2);
			psFitParameter[INTERCEPT].Fix=0;
			
			///Echo 7/22/04 LR_EFFICIENCY_IMPROVE
			///Max 9/09/04 MODIFY_FIXINTERCEPT
			//int nRet = ocmath_non_sigular_ML_regression(vX, nPts, 1, vY, vWeight, vWeight.GetSize(), vCoeff, vCoeff.GetSize(), vANOVA, mCov, mCov.GetNumRows(), mCov.GetNumCols(), sLROptions.ThroughZero);			
			int nRet = ocmath_non_sigular_ML_regression(vX, nPts, 1, vY, vWeight, vWeight.GetSize(), vCoeff, vCoeff.GetSize(), vANOVA, mCov, mCov.GetNumRows(), mCov.GetNumCols(), sLROptions.FixIntercept);			
			///End MODIFY_FIXINTERCEPT
			if (nRet !=  STATS_NO_ERROR)
				return nRet;
			//ocmath_multiple_regression(vX, nPts, 1, vY, vWeight, vWeight.GetSize(), vCoeff, vCoeff.GetSize(), vANOVA, mCov, mCov.GetNumRows(), mCov.GetNumCols(), sLROptions.ThroughZero);			
			///end LR_EFFICIENCY_IMPROVE
			dIntercept = vCoeff[0];
			dSlope = vCoeff[1];
			da_Err = sqrt(mCov[0][0]);
			db_Err = sqrt(mCov[1][1]);
			
			// Get other statistics sums
			if( 0 != stats_linear_regrssion_sum(vX, vY, vWeight, dSumWeight, dxAve, dyAve, dSxx, dSxy, dSyy) )
				return STATS_ERROR_REGRESSION_X_RANGE;
			
			// Residual sum of squares
			rss = vANOVA[4];
	
			// Degrees of freedom
			dof = vANOVA[MULREG_EDOF];
			mdof = vANOVA[MULREG_DOF];
	
			// R-square
			rsq = vANOVA[MULREG_RSQ]/100.0;
			
			//Standard deviation
			SD = vANOVA[MULREG_SD];
		}
		else // With constant Intercept
		{
			/// Max 10/27/04 MODIFY_FIXINTERCEPTAT_CASE			
			//mCov.SetSize(1,1);
			//vCoeff.SetSize(1);
			//psFitParameter[INTERCEPT].Fix=1;
			/////Echo 7/22/04 LR_EFFICIENCY_IMPROVE
			/////Max 9/09/04 MODIFY_FIXINTERCEPT
			////int nRet = ocmath_non_sigular_ML_regression(vX, nPts, 1, vY, vWeight, vWeight.GetSize(), vCoeff, vCoeff.GetSize(), vANOVA, mCov, mCov.GetNumRows(), mCov.GetNumCols(), sLROptions.ThroughZero);	                          						                          
			//int nRet = ocmath_non_sigular_ML_regression(vX, nPts, 1, vY, vWeight, vWeight.GetSize(), vCoeff, vCoeff.GetSize(), vANOVA, mCov, mCov.GetNumRows(), mCov.GetNumCols(), sLROptions.FixIntercept);	                          						                          
			/////End MODIFY_FIXINTERCEPT
			//if (nRet !=  STATS_NO_ERROR)
				//return nRet;
			////ocmath_multiple_regression(vX, nPts, 1, vY, vWeight, vWeight.GetSize(), vCoeff, vCoeff.GetSize(), vANOVA, mCov, mCov.GetNumRows(), mCov.GetNumCols(), sLROptions.ThroughZero);	                          						                          
			/////end LR_EFFICIENCY_IMPROVE			
			/////Max 9/09/04 MODIFY_FIXINTERCEPT
			////dIntercept = 0;
			//dIntercept = sLROptions.FixInterceptAt;
			/////End MODIFY_FIXINTERCEPT
			////dSlope = dcoeff[0];
			//dSlope = vCoeff[0];
			///// End Jay 05/25/04
			//da_Err = 0;
			//db_Err = sqrt(mCov[0][0]);
			//
			psFitParameter[INTERCEPT].Fix=1;
			if( 0 != stats_linear_regrssion_sum(vX, vY, vWeight, dSumWeight, dxAve, dyAve, dSxx, dSxy, dSyy) )
				return STATS_ERROR_REGRESSION_X_RANGE;			
			
			// Intercept		
			dIntercept = sLROptions.FixInterceptAt;
			da_Err = 0;	// Standard error of intercept
							
			// Slope
			dSlope =(dSxy / dSumWeight + dxAve * dyAve - dIntercept * dxAve) / (dSxx / dSumWeight + dxAve^2) ;
				
			// Residual sum of squares
			rss = dSyy - 2 * dSlope * dSxy + dSlope^2 * dSxx;
	
			// Degrees of freedom
			dof = nPts - 1;
			mdof = 0;
	
			// Standard error of slope
			db_Err = sqrt(rss / dof / dSumWeight);	
			
			// R-square
			rsq = dSxy^2 / dSxx / dSyy;
			
			//Standard deviation
			SD = sqrt(rss / dof);
			/// End MODIFY_FIXINTERCEPTAT_CASE
		}
	}
	
	
	// t(1-alpha/2,dof) ------- alpha: significance level
	double tvalue = tTable(1 - (1.0 - sLROptions.Confidence) / 2, dof);//CPY 2/20/04, I think conf should come from LR option

	// Simple statistics results	
	if( sLROptions.ErrBarWeight && !sLROptions.UseReducedChiSq ) 
	{
		da_Err /= SD;
		db_Err /= SD;
	}
	
	
	// Intercept
	psFitParameter[INTERCEPT].Value = dIntercept;
	psFitParameter[INTERCEPT].Error = da_Err; 
	psFitParameter[INTERCEPT].tValue = dIntercept / da_Err;
	psFitParameter[INTERCEPT].LCI = dIntercept - da_Err * tvalue;
	psFitParameter[INTERCEPT].UCI = dIntercept + da_Err * tvalue;
	psFitParameter[INTERCEPT].Prob = 2 * (1 - invt(fabs(psFitParameter[0].tValue), dof));

	// Slope
	psFitParameter[SLOPE].Value = dSlope;
	psFitParameter[SLOPE].Error = db_Err;
	psFitParameter[SLOPE].tValue = dSlope / db_Err;
	psFitParameter[SLOPE].LCI = dSlope - db_Err * tvalue;
	psFitParameter[SLOPE].UCI = dSlope + db_Err * tvalue;
	psFitParameter[SLOPE].Prob = 2 * (1 - invt(fabs(psFitParameter[1].tValue), dof));

	// Stat - Correlation coefficient
	if(psRegStats)
	{
		psRegStats->Rvalue = sqrt(rsq);
		psRegStats->RSqCOD = rsq;
		psRegStats->AdjRSq = 1 - (1 - rsq) * (dof + 1) / dof;
		psRegStats->RMSESD = SD;
		psRegStats->N = nPts;
		psRegStats->DOF = dof;
		psRegStats->SSR = vANOVA[MULREG_RSS];
		psRegStats->ReducedChiSq = vANOVA[MULREG_EMS];
		psRegStats->NormResiduals = vANOVA[MULREG_ENORM];
	}
	
	// ANOVA
	if(psRegANOVA)
	{
		psRegANOVA->AnovaTable[ANOVA_ERR].DOF = dof;
		psRegANOVA->AnovaTable[ANOVA_ERR].SSq = vANOVA[MULREG_RSS];
		psRegANOVA->AnovaTable[ANOVA_ERR].MeanSq = vANOVA[MULREG_EMS];
		
		psRegANOVA->AnovaTable[ANOVA_TOTAL].DOF = vANOVA[MULREG_TDOF];
		psRegANOVA->AnovaTable[ANOVA_TOTAL].SSq = vANOVA[MULREG_TSS];
		psRegANOVA->AnovaTable[ANOVA_TOTAL].MeanSq = NANUM;
		
		psRegANOVA->AnovaTable[ANOVA_MODEL].DOF = vANOVA[MULREG_DOF];
		psRegANOVA->AnovaTable[ANOVA_MODEL].SSq = vANOVA[MULREG_ESS];
		psRegANOVA->AnovaTable[ANOVA_MODEL].MeanSq = vANOVA[MULREG_MMS];
		
		psRegANOVA->Fvalue = vANOVA[MULREG_F];
		psRegANOVA->Pvalue = vANOVA[MULREG_P];

	}
	if(psInternal)
	{
		psInternal->dSumWeight = dSumWeight;
		psInternal->dSxx = dSxx;
		psInternal->dSxy = dSxy;
		psInternal->dSyy = dSyy;
		psInternal->dxAve = dxAve;
		psInternal->dyAve = dyAve;
	}
		
	return STATS_NO_ERROR;
}
*/


/**
		Function to perform multiple linear regression 
		
*/
int stats_multiple_linear_regression(const matrix& mX, const vector& vY, vector& vWT, const LROptions& sLROptions,
FitParameter* psFitParameter, uint nSizeFitParams, RegStats* psRegStats, RegANOVA* psRegANOVA, matrix& mCov, matrix& mCorr, vector& vH, vector& vSE)
{
	UINT nPts = mX.GetNumRows();
	UINT nM = mX.GetNumCols();
	UINT nWT = 0;
	///Echo 7/21/04 MR_AVOID_ERROR_SETTING
	if (vY.GetSize() != nPts)
		return STATS_ERROR_DATA_DIFF_SIZE;
	///end MR_AVOID_ERROR_SETTING
		
	//------ CPY 7/17/04 FIT_PARAM_ARRAY_NEED_SIZE
	//UINT nSize = sizeof(psFitParameter)/sizeof(FitParameter);
	UINT  nSize = nSizeFitParams;
	//------

	int iErr, nCovCorrRow, nCovCorrCol;
	///Max 9/09/04 MODIFY_FIXINTERCEPT
	//if( sLROptions.ThroughZero == 0 ) // With constant
	///Arvin 01/25/08 QA70-10923 CALC_CONFIDENCE_AND_PREDICTION_BANDS
	//We should keep covariance and correlation matrix consistent with NLSF
	//If parameter is fixed, we should set corresponding row and col as missing values 
	/*
	if( sLROptions.FixIntercept == 0 ) // With constant
	///End MODIFY_FIXINTERCEPT
		nCovCorrRow = nCovCorrCol = nM + 1;
	else //without constant
		nCovCorrRow = nCovCorrCol = nM;
	*/
	nCovCorrRow = nCovCorrCol = nM + 1;
	///end CALC_CONFIDENCE_AND_PREDICTION_BANDS
	nWT = vWT.GetSize();
	//if (nWT == 0)
	//{
		//nWT = nPts;
		//vWT.SetSize(nWT);
		//vWT = 1;
	//}
//
	///Echo 10/22/04 USE_POINTER_TO_MATRIX_TO_CENTRALIZE_CALL_TO_OCMATH_FUNC
	//mCov.SetSize(nCovCorrRow, nCovCorrCol);
	//mCorr.SetSize(nCovCorrRow, nCovCorrCol);
	//
	//if (mCov && mCorr)
		//iErr = ocmath_multiple_linear_regression( mX, nPts, nM, vY, vWT, nWT, &sLROptions, psFitParameter, nSize, psRegStats, psRegANOVA, mCov, mCorr, nCovCorrRow, nCovCorrCol); 
	//else if (mCov)
			//iErr = ocmath_multiple_linear_regression( mX, nPts, nM, vY, vWT, nWT, &sLROptions, psFitParameter, nSize, psRegStats, psRegANOVA, mCov, NULL, nCovCorrRow, nCovCorrCol); 
		//else if (mCorr)
			//iErr = ocmath_multiple_linear_regression( mX, nPts, nM, vY, vWT, nWT, &sLROptions, psFitParameter, nSize, psRegStats, psRegANOVA, NULL, mCorr, nCovCorrRow, nCovCorrCol); 
			//else
				//iErr = ocmath_multiple_linear_regression( mX, nPts, nM, vY, vWT, nWT, &sLROptions, psFitParameter, nSize, psRegStats, psRegANOVA); 
	double* pmCov = NULL;
	double* pmCorr = NULL;
	if (mCov)
	{
		mCov.SetSize(nCovCorrRow, nCovCorrCol);
		pmCov = mCov;
	}
	if(mCorr)
	{
		mCorr.SetSize(nCovCorrRow, nCovCorrCol);
		pmCorr = mCorr;
	}
	
	double* pvH = NULL;
	if (vH)
	{
		vH.SetSize(nPts);
		pvH = vH;
	}
	
	double* pvSE = NULL;
	if (vSE)
	{
		vSE.SetSize(nPts);
		pvSE = vSE;
		
	}
	///Echo 9/12/05 MADE_MR_FASTER
	//DescStatResults dsr;
	int nNumPoints = mX.GetNumCols() * mX.GetNumRows();
	double dMean, dCOV;
	//if( iErr = ocmath_desc_stats(mX, nNumPoints, &dsr) )
	if( iErr = ocmath_basic_summary_stats(nNumPoints, mX, NULL, &dMean, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dCOV) )
		return iErr;
	
	if (abs(dMean) > 1e4 && dCOV < 1)
	{
		///Echo 1/2/05 X_NORMALIZE_MECHANISM
		//iErr = ocmath_multiple_linear_regression( mX, nPts, nM, vY, vWT, nWT, &sLROptions, psFitParameter, nSize, psRegStats, psRegANOVA, pmCov, pmCorr, nCovCorrRow, nCovCorrCol); ; 
		matrix mNormX(nPts, nM);
		vector vS(nCovCorrRow);
		///Kevin 11/29/05 CHANGE_OCMATH_FUNCTION_TO_LOWER_CASE
		//iErr = ocmath_Normalize_X(&sLROptions, mX, nPts, nM, mNormX, vS, vS.GetSize());
		iErr = ocmath_normalize_x(&sLROptions, mX, nPts, nM, mNormX, vS, vS.GetSize());
		///End CHANGE_OCMATH_FUNCTION_TO_LOWER_CASE
		if (iErr != STATS_NO_ERROR)
			return iErr;
		FitParameter* psFitParamTemp;
		psFitParamTemp = (FitParameter*)malloc(nSize*sizeof(FitParameter));
		matrix mCovTemp;
		mCovTemp.SetSize(nCovCorrRow, nCovCorrCol);
		iErr = ocmath_multiple_linear_regression( mNormX, nPts, nM, vY, vWT, nWT, &sLROptions, psFitParamTemp, nSize, psRegStats, psRegANOVA, mCovTemp, pmCorr, nCovCorrRow, nCovCorrCol, pvH, pvSE); 
		///Kevin 11/29/05 CHANGE_OCMATH_FUNCTION_TO_LOWER_CASE
		//ocmath_MR_Results_Corrected(&sLROptions, psRegStats, psFitParamTemp, nSize, mCovTemp, nCovCorrRow, nCovCorrCol, vS, vS.GetSize(), psFitParameter, pmCov);
		ocmath_MR_results_corrected(&sLROptions, psRegStats, psFitParamTemp, nSize, mCovTemp, nCovCorrRow, nCovCorrCol, vS, vS.GetSize(), psFitParameter, pmCov);
		///End CHANGE_OCMATH_FUNCTION_TO_LOWER_CASE
		///end X_NORMALIZE_MECHANISM	
		free(psFitParamTemp);///Forest 03/29/05		
	}else
		iErr = ocmath_multiple_linear_regression( mX, nPts, nM, vY, vWT, nWT, &sLROptions, psFitParameter, nSize, psRegStats, psRegANOVA, pmCov, pmCorr, nCovCorrRow, nCovCorrCol, pvH, pvSE);
		
	///end MADE_MR_FASTER
	
	//if (iErr != STATS_NO_ERROR)///Forest 03/29/05
	return iErr;				
}



/**
		Function to perform polynomial regression 
		
*/
int stats_polynomial_fit(const vector& vX, const vector& vY, const vector& vWeightIn,
		int nOrder, const LROptions& sLROptions, FitParameter* psFitParameter, uint nSizeFitParams, 
		RegStats* psRegStats, RegANOVA* psRegANOVA, matrix& mCov, matrix& mCorr, vector& vH)///Arvin 12/12/06 ADD_RESIDUAL_ANALYSIS_FOR_FITTING
{
///Leo 08/02/05 TO_CENTRALIZE_POLYNOMIAL_FIT
	double* pWT = vWeightIn;
	
	if (pWT)
		if (vX.GetSize() != vWeightIn.GetSize())
			return STATS_ERROR_WEIGHT_DIFF_SIZE;
	
	if (vX.GetSize() != vY.GetSize()) 
		return STATS_ERROR_DATA_DIFF_SIZE;
	
	if (mCov || mCorr)
	{
		int nCovCorrRow, nCovCorrCol;
		///Arvin 01/25/08 QA70-10923 CALC_CONFIDENCE_AND_PREDICTION_BANDS
		//We should keep covariance and correlation matrix consistent with NLSF
		//If parameter is fixed, we should set corresponding row and col as missing values 
		/*
		if( sLROptions.FixIntercept == 0 )
			nCovCorrRow = nCovCorrCol = nOrder + 1;
		else
			nCovCorrRow = nCovCorrCol = nOrder;
		*/
		nCovCorrRow = nCovCorrCol = nOrder + 1;
		///end CALC_CONFIDENCE_AND_PREDICTION_BANDS
		if (mCov)
			mCov.SetSize(nCovCorrRow, nCovCorrCol);
		if (mCorr)
			mCorr.SetSize(nCovCorrRow, nCovCorrCol);
	}
	/// Leo 2006-4-7 FIX_NULL_REFERENCE_BUG
	//return ocmath_polynomial_fit(vX.GetSize(), vX, vY, pWT, nOrder, &sLROptions, psFitParameter, nSizeFitParams, psRegStats, psRegANOVA, mCov, mCorr);	
	double* pmCov = NULL;
	double* pmCorr = NULL;
	if (mCov)
	{
		pmCov = mCov;
	}
	if (mCorr)
	{
		pmCorr = mCorr;
	}
	///Arvin 12/12/06 ADD_RESIDUAL_ANALYSIS_FOR_FITTING
	//return ocmath_polynomial_fit(vX.GetSize(), vX, vY, pWT, nOrder, &sLROptions, psFitParameter, nSizeFitParams, psRegStats, psRegANOVA, pmCov, pmCorr);
	double* pvH = NULL;
	if (vH)
	{
		vH.SetSize(vX.GetSize());
		pvH = vH;
	}
	
	return ocmath_polynomial_fit(vX.GetSize(), vX, vY, pWT, nOrder, &sLROptions, psFitParameter, nSizeFitParams, psRegStats, psRegANOVA, pmCov, pmCorr, pvH);
	///end ADD_RESIDUAL_ANALYSIS_FOR_FITTING
	/// end FIX_NULL_REFERENCE_BUG	
//	double dSumWeight;			// The sum of weights
//	double dSxx, dSxy, dSyy;
//	double dxAve, dyAve;		// The average value of x and y
//	
//	if (0 == vWeightIn.GetSize())
//	{
//		vWeightIn.SetSize(vX.GetSize());
//		for (int ii = 0; ii < vX.GetSize(); ii++)
//			vWeightIn[ii] = 1;
//	}
//	
//	if (vX.GetSize() != vWeightIn.GetSize())
//			return STATS_ERROR_WEIGHT_DIFF_SIZE;
//	
//	if(vX.GetSize() != vY.GetSize()) 
//			return STATS_ERROR_DATA_DIFF_SIZE;
//	
//	if( nOrder < 1 ) 
//		return  STATS_ERROR_SETTING;
//	
//	 //---	CPY 7/17/04 FIT_PARAM_ARRAY_NEED_SIZE
//	 if(nSizeFitParams < nOrder + 1)
//	 	return STATS_PARAMS_ARRAY_TOO_SMALL;
//	 //---
//
//	 ///Echo 12/9/04 ERROR_JUDGE_IN_LLVC
//	//if( 0 != stats_linear_regrssion_sum(vX, vY, vWeightIn, dSumWeight, dxAve, dyAve, dSxx, dSxy, dSyy) )
//		//return STATS_ERROR_REGRESSION_X_RANGE;
//
//	// Get the associated independent variable
//	matrix mIndepData(vX.GetSize(), nOrder);
//	int nOrder1 = nOrder + 1;
//
//	vector vtemp(vX.GetSize());
//	vtemp = 1;
//	for( int ii = 0; ii < nOrder; ii++)
//	{
//		vtemp *= vX;
//		mIndepData.SetColumn(vtemp, ii);
//	}
//
//	//int nCovCorrRow, nCovCorrCol;
//	///Max 9/09/04 MODIFY_FIXINTERCEPT
//	//if( sLROptions.ThroughZero == 0 ) // With constant
//	//if( sLROptions.FixIntercept == 0 ) // With constant
//	///End MODIFY_FIXINTERCEPT
//		//nCovCorrRow = nCovCorrCol = nOrder + 1;
//	//else //without constant
//		//nCovCorrRow = nCovCorrCol = nOrder;
//		
//	return stats_multiple_linear_regression(mIndepData,  vY, vWeightIn, sLROptions, psFitParameter, nSizeFitParams, psRegStats, psRegANOVA, mCov, mCorr);
//
///End TO_CENTRALIZE_POLYNOMIAL_FIT

	///Echo 11/18/04 Use stats_multiple_linear_regression instead. MR LLVC is rewriten and the speed is good enough. 
	//--- CPY 10/21/04 USE_POINTER_TO_MATRIX_TO_CENTRALIZE_CALL_TO_OCMATH_FUNC
	//double* pmCov = NULL;
	//double* pmCorr = NULL;
	//if (mCov!=NULL)
	//{
		//mCov.SetSize(nCovCorrRow, nCovCorrCol);
		//pmCov = mCov;
	//}
	//if(mCorr!=NULL)
	//{
		//mCorr.SetSize(nCovCorrRow, nCovCorrCol);
		//pmCorr = mCorr;
	//}
	//return ocmath_non_sigular_multiple_regression(mIndepData, nPts, nOrder, vY, vWeightIn, vWeightIn.GetSize(), 
			//&sLROptions, psFitParameter, nSizeFitParams, psRegStats, psRegANOVA, pmCov, pmCorr, nCovCorrRow, nCovCorrCol);
	//---
	/*
	if (mCov != NULL || mCorr != NULL)
	{
		mCov.SetSize(nCovCorrRow, nCovCorrCol);
		mCorr.SetSize(nCovCorrRow, nCovCorrCol);
	}
	
	if (mCov && mCorr)
		nRet = ocmath_non_sigular_multiple_regression(mIndepData, nPts, nOrder, vY, vWeightIn, vWeightIn.GetSize(), 
			&sLROptions, psFitParameter, nSizeFitParams, psRegStats, psRegANOVA, mCov, mCorr, nCovCorrRow, nCovCorrCol);
	else if (mCov)
			nRet = ocmath_non_sigular_multiple_regression(mIndepData, nPts, nOrder, vY, vWeightIn, vWeightIn.GetSize(), 
				&sLROptions, psFitParameter, nSizeFitParams, psRegStats, psRegANOVA, mCov, NULL, nCovCorrRow, nCovCorrCol);
		else if (mCorr)
			nRet = ocmath_non_sigular_multiple_regression(mIndepData, nPts, nOrder, vY, vWeightIn, vWeightIn.GetSize(), 
				&sLROptions, psFitParameter, nSizeFitParams, psRegStats, psRegANOVA, NULL, mCorr, nCovCorrRow, nCovCorrCol);
			else
				nRet = ocmath_non_sigular_multiple_regression(mIndepData, nPts, nOrder, vY, vWeightIn, vWeightIn.GetSize(), 
						&sLROptions, psFitParameter, nSizeFitParams, psRegStats, psRegANOVA);
	
	
		return nRet;
	*/
}

/**
   		 Function to find if exists duplicated x for given y.
*/
bool check_duplicate_x(const double dY, const Dataset dsSrcX, const Dataset dsSrcY)
{
	/// *** Terminate if error setting
	if( (dY < min(dsSrcY)) || (dY > max(dsSrcY)) ) return NANUM;
	
	// Now check if the dY for which user seeks x, exists duplicated exact point in the fit dataset
	int iSize = dsSrcY.GetSize();
	Dataset dsTemp = dsSrcY;
	int iIsPoint;
	bool bFind = false;
	do
	{
		///Sandy 2006-8-1 CHANGE_TO_USE_VECTOR_BASE_FUNCTION
		//iIsPoint = Data_list(dY, &dsTemp);
		/// If it is a point, just return the x corresponding to the index
		vector<uint> vecIndex;
		iIsPoint = dsTemp.Find(vecIndex, dY);
		//end 2006-8-1 CHANGE_TO_USE_VECTOR_BASE_FUNCTION
		
		if(iIsPoint != -1) 
		{
			if (bFind == true)
				return true;
			else
			{
				bFind = true;
				dsTemp[iIsPoint] = 0; //if x value found is not in the interval, find the next one. 
			}				
		}
	}while (iIsPoint != -1);
			
	// Now find if there may exists duplidate x values for y
	bFind = false;
	for(int ik = 0; ik < (iSize - 1); ik++)
	{
		double dSign = (dsSrcY[ik] - dY) * (dsSrcY[ik+1] - dY);//find exact point of curve of x interval in which yval belongs
		if (dSign < 0) 
		{
			if (bFind == true)
				return true;
			else
				bFind = true;
				
		}
	}
	return false;
	
}

/**
   		 Function to compute a Y value from a standard curve given X.
*/
double compute_x_from_y(int& ifault, double dY, double dXLo, double dXUp, const Dataset dsSrcX, const Dataset dsSrcY, int nMode, double dSmoothingFactor, vector* pvWeights)
{
	//CPY 10/26/04
	// this function is giving compile error, due to const related, and why was Dataset used? why pass in directly instead of by ref?
	// also, this type of function should be written in VC
	return 0;
	/*
	vector vSrcX, vSrcY;
	vSrcX = dsSrcX;
	vSrcY = dsSrcY;
	// If yval is outside of [ymin,ymax] of yfit dataset, return missing value
	if( (dY < min(dsSrcY)) || (dY > max(dsSrcY)) )	return NANUM;
	
	//if dXLo > dXUp, return missing value
	if (dXLo > dXUp)
	{
		ifault = STATS_ERROR_SETTING;
		return NANUM;		
	}
	
	// Now check if the dY for which user seeks x, is an exact point in the fit dataset
	int iSize = vSrcY.GetSize();
	Dataset dsTemp = dsSrcY;
	int iIsPoint;
	bool bFind;
	do
	{
		bFind =  false;
		iIsPoint = Data_list(dY, &dsTemp);
		// If it is a point, just return the x corresponding to the index
		if(iIsPoint != -1) 
		{
			if (dsSrcX[iIsPoint] > dXLo && dsSrcX[iIsPoint] < dXUp)//if x value found is in the interval, just return the x. 
			{
				bFind = true;
				break;	
			}else
				dsTemp[iIsPoint] = 0; //if x value found is not in the interval, find the next one. 
		}
	}while (iIsPoint != -1 );
	
	if (bFind == true)
		return dsSrcX[iIsPoint];
		
		
	// Now find (first) x interval in which yval belongs
	for(int ik = 0; ik < (iSize - 1); ik++)
	{
		double dSign = (dsSrcY[ik] - dY) * (dsSrcY[ik+1] - dY);//find exact point of curve of x interval in which yval belongs
		if (dSign < 0) 
		{
			if (dsSrcX[ik] > dXLo && dsSrcX[ik] < dXUp) //if lower interval found is in the given interval, let intersection be updated interval
			{
				dXLo = dsSrcX[ik];
				if (dsSrcX[ik+1] < dXUp) // if given upper bounding x value is larger than exact point of upper x interval, let the exact point be upper x interval
				{
					dXUp = dsSrcX[ik+1];						
					break;
				}else break;
			}else //if exact point of lower x interval is less, let given value be the lower x interval
			{
				if (dsSrcX[ik+1] > dXLo && dsSrcX[ik+1] < dXUp)// if upper interval found is in the given interval, let the intersection be updated interval.
				{
					dXUp = dsSrcX[ik+1];						
					break;					
				}//given interval does not includes this found x interval. Go on next search
			}//end judge whether the found x interval is in the given interval		
		}//end if the x interval is which yval belongs
	}//end loop
	if (ik == (iSize - 1)) return NANUM;
	
	// Now start iterative process in the interval [dXLo, dXUp]
	double dXMid, dYatXLo, dYatXHi, dYOld, dYNew, dYChange;
	dXMid = (dXLo + dXUp) / 2.0;	
	vector vX(1), vY(1);
	int nSize = 1;
	int	nSrcSize = dsSrcX.GetSize();
	vX[0] = dXLo;
	int nRet =  ocmath_interpolate(vX, vY, nSize, vSrcX, vSrcY, nSrcSize, nMode, dSmoothingFactor, *pvWeights);
	if (nRet != STATS_NO_ERROR)
	{
		ifault = nRet;
		return NANUM;
	}
	else
		dYatXLo = vY[0];	
	vX[0] = dXUp;
	nRet =  ocmath_interpolate(vX, vY, nSize, vSrcX, vSrcY, nSrcSize, nMode, dSmoothingFactor, *pvWeights);
	if (nRet != STATS_NO_ERROR)
	{
		ifault = nRet;
		return NANUM;
	}
	else
		dYatXHi = vY[0];
	vX[0] = dXMid;
	nRet =  ocmath_interpolate(vX, vY, nSize, vSrcX, vSrcY, nSrcSize, nMode, dSmoothingFactor, *pvWeights);
	if (nRet != STATS_NO_ERROR)
	{
		ifault = nRet;
		return NANUM;
	}
	else
		dYOld = vY[0];
			
	int ii = 0;
	do
	{
		ii++;
		if(dYatXLo > dYatXHi)
		{
			if(dYOld > dY) dXLo = dXMid;
			else dXUp = dXMid;
			
		}
		else
		{
			if(dYOld < dY) dXLo = dXMid;
			else dXUp = dXMid;
		}
		vX[0] = dXLo;
		nRet =  ocmath_interpolate(vX, vY, nSize, vSrcX, vSrcY, nSrcSize, nMode, dSmoothingFactor, *pvWeights);
		if (nRet != STATS_NO_ERROR)
		{
			ifault = nRet;
			return NANUM;
		}
		else
			dYatXLo = vY[0];	
		vX[0] = dXUp;
		nRet =  ocmath_interpolate(vX, vY, nSize, vSrcX, vSrcY, nSrcSize, nMode, dSmoothingFactor, *pvWeights);
		if (nRet != STATS_NO_ERROR)
		{
			ifault = nRet;
			return NANUM;
		}
		else
			dYatXHi = vY[0];		
		dXMid = (dXLo + dXUp) / 2.0;
		vX[0] = dXMid;
		nRet =  ocmath_interpolate(vX, vY, nSize, vSrcX, vSrcY, nSrcSize, nMode, dSmoothingFactor, *pvWeights);
		if (nRet != STATS_NO_ERROR)
		{
			ifault = nRet;
			return NANUM;
		}
		else
			dYNew = vY[0];
		
		// Compute percentage change in y
		dYChange = 100 * fabs( (dYNew - dY) / dY);
		dYOld = dYNew;
		// Loop until % change in y is less than 0.0001 and loop has not been executed 100 times
	} while ( (dYChange > 0.00001) && (ii < 200) ); 		
	
	// If number of iterations equals 200, did not converge and so return missing value
	if(ii == 200) return NANUM;
	return dXMid;	// Return the current dXMid value - this is the x corresponding to dY
	*/
	
}


/////////////////////////////////////////////////////////////////////////////////
//
//				t-Test Functions									
//
/////////////////////////////////////////////////////////////////////////////////
/**
		Function to perform a One Sample t-Test.
*/

int stats_one_sample_t_test(const vector& vData, const HypotTestOptions& opt,
	 tTestDescStats& DescRes, tTestResults& tTestRes)
{
	tTestMean1SampleResults res;
	int nRet = ocmath_one_sample_t_test( vData, vData.GetSize(), &opt, &res );
	//Sonic	12/22/04	COMMENT_TO_FIX_BUG	///Echo 4/1/05 UNDO_SONIC_COMMENT
	if (nRet != STATS_NO_ERROR)
		return nRet;
	else
	{
		DescRes.N = res.N;
		DescRes.Mean = res.Mean;
		DescRes.SD = res.SD;
		DescRes.SEM = res.SEM;
		tTestRes.DOF = res.DOF;
		tTestRes.tValue = res.tValue;
		tTestRes.Prob = res.Prob;
	}
//End	COMMENT_TO_FIX_BUG ///Echo 4/1/05 UNDO_SONIC_COMMENT
	return STATS_NO_ERROR;
}


/**
		Function to compute confidece levels for a mean.
*/
///Kevin 09/01/05 CHANGE_STATISTICS_NAME
//int stats_confidence_levels_for_mean(const vector vLevels, const tTestMean1SampleResults& tTestRes,	 vector& vLCIs, vector& vUCIs)
int stats_confidence_levels_for_mean(const vector vLevels, const tTestMean1SampleResults& tTestRes,
	 vector& vLCLs, vector& vUCLs)
 ///End CHANGE_STATISTICS_NAME
{
	int nLevels = vLevels.GetSize();
///Kevin 09/01/05 CHANGE_STATISTICS_NAME
	//vLCIs.SetSize(nLevels);
	//vUCIs.SetSize(nLevels);
	//return ocmath_confidence_levels_for_mean(vLevels, nLevels, &tTestRes, vLCIs, vUCIs);
	vLCLs.SetSize(nLevels);
	vUCLs.SetSize(nLevels);
	return ocmath_confidence_levels_for_mean(vLevels, nLevels, &tTestRes, vLCLs, vUCLs);
///End CHANGE_STATISTICS_NAME
}


/**
		Function to perform a Chi-Square Test for equal variance.
*/
///Echo 7/6/06 temp_comment
//int stats_chi_square_test(const HypotTestOptions& opt, const tTestMean1SampleResults& tTestRes, ChiSqVar1SampleResults& res)
//{
	//return ocmath_chi_square_test( &opt, &tTestRes, &res );
//}


/**
		Function to compute confidece levels for a variance.
*/
///Echo 7/6/06 temp_comment
/////Kevin 09/01/05 CHANGE_STATISTICS_NAME
////int stats_confidence_levels_for_variance(const vector vLevels, const tTestMean1SampleResults& tTestRes, vector& vLCIs, vector& vUCIs)
//int stats_confidence_levels_for_variance(const vector vLevels, const tTestMean1SampleResults& tTestRes, vector& vLCLs, vector& vUCLs)
/////End CHANGE_STATISTICS_NAME
//{
	//int nSize = vLevels.GetSize();
/////Kevin 09/01/05 CHANGE_STATISTICS_NAME
	////vLCIs.SetSize( nSize );
	////vUCIs.SetSize( nSize );
	////return ocmath_confidence_levels_for_variance( vLevels, nSize, &tTestRes, vLCIs, vUCIs );
	//vLCLs.SetSize( nSize );
	//vUCLs.SetSize( nSize );
	//return ocmath_confidence_levels_for_variance( vLevels, nSize, &tTestRes, vLCLs, vUCLs );
/////End CHANGE_STATISTICS_NAME
//}
//

/**
		Function to compute actual or hypothetical powers for a t-Test.
*/
int stats_tTest_powers(const double dSigLevel, const int iTailType, const vector& vSampleSizes, // Use vector double (vSampleSizes) for convenience
	 tTestMean1SampleResults& tTestRes, vector& vPowers)                                       //   as string::GetTokens requires doubles

{
	int nSize = vSampleSizes.GetSize();
	
	vPowers.SetSize( nSize );
	return ocmath_tTest_powers( dSigLevel, iTailType, vSampleSizes, nSize, &tTestRes, vPowers );
}


/**
		Function to perform a Two Sample t-Test for equal variance.
*/
int stats_two_sample_t_test(const vector& vData1, const vector& vData2, const HypotTestOptions& opt,
	 tTestDescStats& DescStats1, tTestDescStats& DescStats2, tTestDescStats& DescStatsDiff, tTestResults& tTestRes, tTestResults& tTestResUneqVar)
{
	tTestMean2SampleResults res;	
	uint nSize1 = vData1.GetSize();
	uint nSize2 = vData2.GetSize();
	
	tTestResults* ptTestResUneqVar;
	ptTestResUneqVar = (NULL != tTestResUneqVar)? &tTestResUneqVar : NULL;	
	int nRet = ocmath_two_sample_t_test(vData1, nSize1, vData2, nSize2, &opt, &res, ptTestResUneqVar);
	if (nRet != STATS_NO_ERROR)
		return nRet;
	else
		stats_get_ttest2_results(res, DescStats1, DescStats2, DescStatsDiff, tTestRes);	
		
	return STATS_NO_ERROR;
}
/**
		Function to preprocess Descriptive Statistics results in two sample ttest as required
*/
void stats_set_ttest2_results(const tTestDescStats& DescStats1, const tTestDescStats& DescStats2, 
const tTestDescStats& DescStatsDiff, const tTestResults& tTestRes, tTestMean2SampleResults& res)
{
		res.N1 = DescStats1.N;
		res.Mean1 = DescStats1.Mean;
		res.SD1 = DescStats1.SD;
		res.SEM1 = DescStats1.SEM;
		res.N2 = DescStats2.N ;
		res.Mean2 = DescStats2.Mean;
		res.SD2 = DescStats2.SD;
		res.SEM2 = DescStats2.SEM;
		res.DOF = tTestRes.DOF;
		res.tValue = tTestRes.tValue;
		res.Prob = tTestRes.Prob;
	
}

/**
		Function to preprocess two sample ttet in Descriptive Statistics results
*/
void stats_get_ttest2_results(const tTestMean2SampleResults& res, tTestDescStats& DescStats1, tTestDescStats& DescStats2, 
tTestDescStats& DescStatsDiff, tTestResults& tTestRes)
{
		DescStats1.N = res.N1;
		DescStats1.Mean = res.Mean1;
		DescStats1.SD = res.SD1;
		DescStats1.SEM = res.SEM1;
		DescStats2.N = res.N2;
		DescStats2.Mean = res.Mean2;
		DescStats2.SD = res.SD2;
		DescStats2.SEM = res.SEM2;
		DescStatsDiff.N = -1;
		DescStatsDiff.Mean = res.Mean1 - res.Mean2;
		DescStatsDiff.SD = NANUM;
		DescStatsDiff.SEM = res.SEM1 - res.SEM2;
		tTestRes.DOF = res.DOF;
		tTestRes.tValue = res.tValue;
		tTestRes.Prob = res.Prob;
	
}

/**
		Function to compute confidece levels for the mean of vData1 minus the mean of vData2 for equal variance.
*/
///Kevin 09/01/05 CHANGE_STATISTICS_NAME
//int stats_two_sample_confidence_levels_for_mean_equal_var(const vector vLevels, const tTestMean2SampleResults& tTestRes, vector& vLCIs, vector& vUCIs)
int stats_two_sample_confidence_levels_for_mean_equal_var(const vector vLevels, const tTestMean2SampleResults& tTestRes,
	 vector& vLCLs, vector& vUCLs)
 ///End CHANGE_STATISTICS_NAME
{
	int nLevels = vLevels.GetSize();
///Kevin 09/01/05 CHANGE_STATISTICS_NAME
	//vLCIs.SetSize( nLevels );
	//vUCIs.SetSize( nLevels );
	//return ocmath_two_sample_confidence_levels_for_mean_equal_var(vLevels, nLevels, &tTestRes, vLCIs, vUCIs);///Iris 4/9/04
	vLCLs.SetSize( nLevels );
	vUCLs.SetSize( nLevels );
	return ocmath_two_sample_confidence_levels_for_mean_equal_var(vLevels, nLevels, &tTestRes, vLCLs, vUCLs);
///End CHANGE_STATISTICS_NAME
}

int stats_two_sample_confidence_levels_for_mean_unequal_var(const vector vLevels, const tTestMean2SampleResults& tTestRes,
	 vector& vLCLs, vector& vUCLs)
{
	int nLevels = vLevels.GetSize();
	vLCLs.SetSize( nLevels );
	vUCLs.SetSize( nLevels );
	return ocmath_two_sample_confidence_levels_for_mean_unequal_var(vLevels, nLevels, &tTestRes, vLCLs, vUCLs);	
}

/**
		Function to perform a Chi-Square Test for equal variance.
*/
int stats_two_sample_chi_square_test(const HypotTestOptions& opt, const tTestMean2SampleResults& tTestRes, ChiSqVar2SampleResults& res)
{
	return ocmath_two_sample_chi_square_test(&opt, &tTestRes, &res);
}


/**
		Function to perform a F Test for equal variance.
*/
///Echo 7/6/06 temp_comment
//int stats_f_test(const HypotTestOptions& opt, const tTestMean2SampleResults& tTestRes, FVarResults& res)
//{
	//return ocmath_f_test(&opt, &tTestRes, &res);
//}


/**
		Function to compute confidece levels for a variance using f statistic.
*/	
///Echo 7/6/06 temp_comment
///Kevin 09/01/05 CHANGE_STATISTICS_NAME
//int stats_two_sample_confidence_levels_for_variance(const vector vLevels, const tTestMean2SampleResults& tTestRes, vector& vLCIs, vector& vUCIs)
//int stats_two_sample_confidence_levels_for_variance(const vector vLevels, const tTestMean2SampleResults& tTestRes,
	 //vector& vLCLs, vector& vUCLs)
 /////End CHANGE_STATISTICS_NAME
//{
	//int nSize = vLevels.GetSize();
/////Kevin 09/01/05 CHANGE_STATISTICS_NAME
	////vLCIs.SetSize(nSize);
	////vUCIs.SetSize(nSize);
	////return ocmath_two_sample_confidence_levels_for_variance(vLevels, nSize, &tTestRes, vLCIs, vUCIs);
	//vLCLs.SetSize(nSize);
	//vUCLs.SetSize(nSize);
	//return ocmath_two_sample_confidence_levels_for_variance(vLevels, nSize, &tTestRes, vLCLs, vUCLs);
/////End CHANGE_STATISTICS_NAME
//}	


/**
		Function to compute actual or hypothetical powers for a t-Test.
*/

int stats_two_sample_tTest_powers(const double dSigLevel, const HypotTestOptions& opt, const vector& vSampleSizes, // Use vector double (vSampleSizes) for convenience
	 tTestMean2SampleResults& tTestRes, vector& vPowers)                                        //   as string::GetTokens requires doubles
{
	int nSize = vSampleSizes.GetSize();
	vPowers.SetSize(nSize);
	return ocmath_two_sample_tTest_powers(dSigLevel, &opt, vSampleSizes, nSize, &tTestRes, vPowers);
}


/**
		Function to perform a Pair Two Sample t-Test.
*/
int stats_two_sample_t_test_pair(const vector& vData1, const vector& vData2, const HypotTestOptions& opt,
	 tTestMean2SampleResults& res)
{
	return ocmath_two_sample_t_test_pair(vData1, vData1.GetSize(), vData2, vData2.GetSize(), &opt, &res);
}


/**
		Function to compute confidece levels for mean in pair two sample t-Test.
*/
///Kevin 09/01/05 CHANGE_STATISTICS_NAME
//int stats_two_sample_confidence_levels_for_mean_pair(const vector vLevels, const tTestMean2SampleResults& tTestRes, const HypotTestOptions& opt,vector& vLCIs, vector& vUCIs)
int stats_two_sample_confidence_levels_for_mean_pair(const vector vLevels, const tTestMean2SampleResults& tTestRes,
	 const HypotTestOptions& opt,vector& vLCLs, vector& vUCLs)
 ///End CHANGE_STATISTICS_NAME
{
	int nSize = vLevels.GetSize();
///Kevin 09/01/05 CHANGE_STATISTICS_NAME
	//vLCIs.SetSize( nSize );
	//vUCIs.SetSize( nSize );	
	//return ocmath_two_sample_confidence_levels_for_mean_pair(vLevels, nSize, &tTestRes,&opt, vLCIs, vUCIs);
	vLCLs.SetSize( nSize );
	vUCLs.SetSize( nSize );	
	return	ocmath_two_sample_confidence_levels_for_mean_pair(vLevels, nSize, &tTestRes,&opt, vLCLs, vUCLs);
///End CHANGE_STATISTICS_NAME
}

/**
		Function to compute power and samplesize in pair two sample t-Test.
*/

int stats_two_sample_tTest_powers_pair(const double dSigLevel, const HypotTestOptions& opt, const vector& vSampleSizes, // Use vector double (vSampleSizes) for convenience
	 tTestMean2SampleResults& tTestRes, vector& vPowers)
{
    int nSize = vSampleSizes.GetSize();
	vPowers.SetSize( nSize );
	return ocmath_two_sample_tTest_powers_pair(dSigLevel, &opt, vSampleSizes, vSampleSizes.GetSize(), &tTestRes, vPowers);
}


/////////////////////////////////////////////////////////////////////////////////
//
//				Normality Test Functions									
//
/////////////////////////////////////////////////////////////////////////////////	 
/**
		Function to perform a Shapiro-Wilk Normality Test using inhouse code
*/
///Echo 9/12/06 NAG_FUNC_MOVED_TO_LLVC
int stats_shapiro_wilk_test(const vector& vData, NormTestResults& SWRes, int iOption)
{
	return ocmath_shapiro_wilk_test(vData.GetSize(), vData, &SWRes, iOption);
}

//int stats_shapiro_wilk_test(const vector& vData, NormTestResults& SWRes, int iOption)
//{
	//int nSize =vData.GetSize();
	////SWRes.DOF=nSize-1;	Max 08/24/05 MODIFY_DOF_OF_SHAPIRO_WILK_TEST
	//SWRes.DOF=nSize;
	//SWRes.TestStat = NANUM;
	//SWRes.Prob = NANUM;
	//
	//double shapiro_wilk_w;
	//double shapiro_wilk_pw;
	//int iFault;
	//double dMin, dMax;
	//vData.GetMinMax(dMin, dMax);
	//
	////when the number of observations <3, ocmath_shapiro_wilk_test_AS181 and ocmath_shapiro_wilk_test_R94 don't support. So return error
	//if(nSize<3)
		//return STATS_ERROR_TOO_FEW_DATA_PTS;
	//
	  //else if((nSize>5000)||((nSize>2000)&&(iOption==1)))
		//return STATS_ERROR_TOO_MANY_DATA_PTS;
	////All the values in vector must not be equal
	//else if (0 == (dMax - dMin))
		//return STATS_ERROR_ALL_ELEMENT_EQUAL;
	//else
	//{
		//if (iOption==0)			//Algorithm R94
		//{	
			//vector vDataCpy1;	//to devoid the changing of vData
			//vDataCpy1 = vData;
			//vDataCpy1.Sort();
			//ocmath_shapiro_wilk_test_R94(vDataCpy1,nSize,&shapiro_wilk_w,&shapiro_wilk_pw,&iFault);		
		//}
///*		else if (iOption==1)	//Algorithm AS181
			 //{
			 	//vector<float> vDataCpy2;	//to devoid the changing of vData
				//vDataCpy2 = vData;
				//vDataCpy2.Sort();
				//float f_shapiro_wilk_w;
				//float f_shapiro_wilk_pw;
			 	//ocmath_shapiro_wilk_test_AS181(vDataCpy2,nSize,&f_shapiro_wilk_w,&f_shapiro_wilk_pw,&iFault);
			 	//shapiro_wilk_w = f_shapiro_wilk_w;
			 	//shapiro_wilk_pw = f_shapiro_wilk_pw;
			 //}
//*/			//when success
		//else if (iOption == 1)
		//{
			//vector vDataCpy2;
			//vDataCpy2 = vData;
			//vDataCpy2.Sort();
			//vector vA;
			//vA.SetSize(nSize);
			//iFault = onag_shapiro_wilk_test( nSize, vDataCpy2, TRUE, vA, &shapiro_wilk_w, &shapiro_wilk_pw );
			//
		//}
		//if(0==iFault)
		//{
			//SWRes.TestStat=shapiro_wilk_w;
			//SWRes.Prob=shapiro_wilk_pw;
			//return STATS_NO_ERROR;
		//}
		//else 
			 //return STATS_ERROR_SHAPIRO_WILK;
	//}
//}
///end NAG_FUNC_MOVED_TO_LLVC

/**
		Function to perform a Kolmogorov-Smirnov Normality Test
*/

/// Max 10/13/04 REPLACE_KOLMOGOROV_SMIRNOV_TEST_WITH_NAG_FUNTION
/*
double fun(double dData)
{
	return (erf(dData/sqrt(2))+1)/2;
}

int stats_kolmogorov_smirnov_test(const vector& vData, NormTestResults& KSRes)
{

	double pd=0.0;
	double pprob;
	int n=vData.GetSize();
	
	if(ocmath_ksone(n, vData, fun, &pd, &pprob ))
	{
	
		KSRes.DOF=n;
		KSRes.TestStat=pd;
		KSRes.Prob=pprob;
		
		if (0 == n)
			KSRes.Prob = NANUM;

		return STATS_NO_ERROR;
	}
	else return STATS_FALSE;
}
*/
///Echo 9/12/06 NAG_FUNC_MOVED_TO_LLVC
int stats_kolmogorov_smirnov_test(const vector& vData, const int paratype, vector& par, NormTestResults& KSRes)
{
	return ocmath_kolmogorov_smirnov_test(vData.GetSize(), vData, &KSRes, paratype, 2, par);
	
}
//int stats_kolmogorov_smirnov_test(const vector& vData, const int paratype, vector& par, NormTestResults& KSRes)
//{
	//double pd, pz, pprob;
	//int nRet, n = vData.GetSize();
	//NagError fail;
	//
	//if (paratype == PARA_ESTIMATED)
	//{
		//nag_1_sample_ks_test(n, vData, Nag_Normal, par, Nag_ParaEstimated, Nag_TestStatisticsDAbs, &pd, &pz, &pprob, &fail);
		//KSRes.DOF = n;
		//KSRes.TestStat = pd;
		//KSRes.Prob = pprob;
		//nRet = fail.code;
	//}
	//else if (paratype == PARA_SUPPLIED)
		 //{
		 	//nag_1_sample_ks_test(n, vData, Nag_Normal, par, Nag_ParaSupplied, Nag_TestStatisticsDAbs, &pd, &pz, &pprob, &fail);
			//KSRes.DOF = n;
			//KSRes.TestStat = pd;
			//KSRes.Prob = pprob;
			//nRet = fail.code;	
		 //}
		 //else return STATS_FALSE;
	//return nRet;
//}
/// End REPLACE_KOLMOGOROV_SMIRNOV_TEST_WITH_NAG_FUNTION
///end NAG_FUNC_MOVED_TO_LLVC

int stats_lilliefors_test(const vector& vData, NormTestResults& LFRes)
{
	return ocmath_lilliefors_test(vData,vData.GetSize(),&LFRes);
}

/////////////////////////////////////////////////////////////////////////////////
//
//				One Way ANOVA Functions									
//
/////////////////////////////////////////////////////////////////////////////////
/**
	Function to get and pre-process all required input datasets
*/
int stats_anova_get_data_in_multiple_columns(const matrix mData, 
									vector& vData, int& iLevels, vector<int>& vIndex)
{
	int nColumns, nRows, ii;
	nColumns = mData.GetNumCols();
	nRows = mData.GetNumRows();
	
	/// Iris 10/05/04
	if(0 == nColumns || 0 == nRows)
		return STATS_ERROR_INPUT_SIZE;
	
	/// *** Get the number of group, which is equal to column of input matrix.
	iLevels = nColumns; 

	/// *** Get size of each group
	vIndex.SetSize(iLevels);
	for (ii = 0; ii < iLevels; ii++)
	{
		//If group are not equal size
		vector vTemp;
		vTemp.SetSize(nRows);
		mData.GetColumn(vTemp, ii);
		vTemp.Trim();
		int nSize = vTemp.GetSize();
		vIndex[ii] = nSize;
	}
	
	/// *** Get Data,put data in vector
	vData.SetSize(nColumns * nRows);
	mData.GetAsVector(vData, false);
	vData.Trim();
	
	return STATS_NO_ERROR;
	
}

/**
		Function to get and pre-process all required input datasets
*/
int stats_anova_get_data_in_group(const vector& vOrgData, const vector& cvGroup, vector& vData, int& iLevels, vector<int>& vIndex)
{
	//--- CPY 10/26/04 CONST_CHECK_ENABLE_NEED_FIXING_CODES
	vector vGroup;
	vGroup = cvGroup;
	//---
	/// *** Return error if error setting
	int nDataSize, nGroupSize;
	nDataSize = vOrgData.GetSize();
	nGroupSize = vGroup.GetSize();
	if (nDataSize != nGroupSize)
		return STATS_ERROR_DATA_DIFF_SIZE;

	/// *** Initial parameters.
	vData.SetSize(nDataSize);//Set size for vData
	vIndex.SetSize(nDataSize);//Set initial size for vIndex;
	vData = NANUM;
	iLevels = -1;
	vIndex = -1;

	/// *** Preprocess data as anova function required 
	
	//Sort vOrgData according to vGroup
	Curve csvTemp(vGroup, vOrgData);
	csvTemp.Sort();
	vData = csvTemp;
	Dataset dsTemp;
	csvTemp.AttachX(dsTemp);
	vGroup = dsTemp;
	
	//Compute iLevels and vIndex
	iLevels = 1;
	vIndex = 1;
	for (int ii = 1; ii < nDataSize; ii++)
	{
		if (vGroup[ii] != vGroup[ii-1])
		{						
			iLevels++;
		}else
			vIndex[iLevels-1]++;
	}
	vIndex.SetSize(iLevels);//resize vIndex as iLevels
	
	return STATS_NO_ERROR;
	
}


/**
		Function to perform a One-Way ANOVA.
*/
int stats_anova_one_way(const 	vector& vData,
						int 	iLevels,
						const 	vector<int>& vIndex,
						ANOVADescStats *DescStats,uint nSizeDescStat,
						ANOVAStats& TotalStats,
						RegANOVA	*ANOVATable)
{
	int TotalSize = 0;
	for (int i = 0; i< iLevels; i++)
	{	
		TotalSize += vIndex[i];
	}
	
	if (iLevels <= 0)
		return STATS_ERROR_BAD_PARAM;	
	
	if (TotalSize != vData.GetSize())
		return STATS_ERROR_DATA_DIFF_SIZE;
		
	int ret = ocmath_anova_one_way(vData, vData.GetSize(), iLevels, vIndex, DescStats, nSizeDescStat, &TotalStats, ANOVATable);

	return ret;
}


/**
		Function to test for equal variance of One-Way ANOVA. 
*/		
int stats_anova_equal_variance_test(const 	vector& vData,
									int 	iLevels,
									const 	vector<int>& vIndex,
									ANOVADescStats *DescStats,uint nSizeDescStat, 
									RegANOVA *ANOVATable,
									int 	iOption)
{
	int TotalSize = 0;
	for (int i = 0; i< iLevels; i++)
	{	
		TotalSize += vIndex[i];
	}
	
	if (iOption <= 0 && iLevels <= 0)
		return STATS_ERROR_BAD_PARAM;
	
	if (TotalSize != vData.GetSize())
		return STATS_ERROR_DATA_DIFF_SIZE;
	
	int ret = ocmath_anova_equal_variance_test(vData, vData.GetSize(),
												iLevels,
												vIndex,
												DescStats, nSizeDescStat,
												ANOVATable,
												iOption);		
	return ret;												
}

/**
		Function to compute One-Way ANOVA Mean Comparisons.
*/

/// Max 03/16/05 REWRITE_ONEWAYANOVA_MEAN_COMPARISON
/*
int ANOVA_Compute_StdErr_MeanDiffs( const int nGroupSize, const double dMSE, const int* pNPTS, const uint nSizeNPTs, double* pC, const uint nSizeC )
{
	int ii, jj; 
	
	// Loop through lower triangle of "matrix" indexing through vector as if it were a matrix
	//  For nGroupSize == 3, Non-zero numbers identify vector index for elements in Lower triangle of matrix:
	//     index = Col + Row * NumCols:                     O O O
	//														3 O O
	//														6 7 O
	
	for( ii = 1; ii < nGroupSize; ii++ )			// Loop on rows of "matrix" starting with second row (ii==1)
	{
		for( jj = 0; jj < ii; jj++ )			// Loop on columns of "matrix"
		{
			// Compute standard errors of the differences of the means 
			pC[jj + ii * nGroupSize] = sqrt( dMSE * ( 1.0 / pNPTS[ii] + 1.0 / pNPTS[jj] ) );
		}
	}

	return STATS_NO_ERROR;
}

int anova_mean_comparison(const	 double *pMean,	uint nSizeMean,
											const	 int *pNPTS, uint nSizeNPTs,
											const	 double dMSE, 
											const	 int nDOF,
											const	 double dAlpha,
											const	 int nGroupSize,
											double	 *pMeanDiffs, uint nSizeMeanDiff,
											double	 *pLCIs, uint nSizeLCI,
											double	 *pUCIs, uint nSizeUCI,
											double	 *pSe, uint nSizeSe,
											double	 *pSig, uint nSizeSig,
											int		 nIntervalType
											)
{
	/// *** Initialize variables
	int jj;
	for (jj =0; jj < nSizeMeanDiff; jj++)
	   pMeanDiffs[jj] = NANUM;

	for (jj =0; jj < nSizeLCI; jj++)
	   pLCIs[jj] = NANUM;

	for (jj =0; jj < nSizeUCI; jj++)
	   pUCIs[jj] = NANUM;

	for (jj =0; jj < nSizeSe; jj++)
	   pSe[jj] = NANUM;

	for (jj =0; jj < nSizeSig; jj++)
	   pSig[jj] = NANUM;

	//// *** Terminate if error setting
	if (dAlpha <= 0 || dAlpha >= 1)
		return STATS_ERROR_SETTING;

	if (nIntervalType != 1 && nIntervalType != 2 && nIntervalType != 3 && nIntervalType != 4)
		return STATS_ERROR_SETTING;

	if ((nSizeMean < nGroupSize) || (nSizeNPTs < nGroupSize))
		return STATS_PARAMS_ARRAY_TOO_SMALL;

	int nMeanDiffSize = nGroupSize*(nGroupSize - 1)/2;
	if ((nSizeMeanDiff < nMeanDiffSize) || (nSizeLCI < nMeanDiffSize) || (nSizeUCI < nMeanDiffSize) || (nSizeSe < nMeanDiffSize) || (nSizeSig < nMeanDiffSize))
		return STATS_PARAMS_ARRAY_TOO_SMALL;


	///////////Get Distribution
	
	///////////Get Mean Difference and SE
	int i,j, iErr;
	int nIndex = 0;
	for(i = 0; i< nGroupSize; i++)
	{
		for(j = 0; j< nGroupSize; j++)
		{
			if (i==j) continue;
			pMeanDiffs[nIndex] = pMean[i]-pMean[j];
			pSe[nIndex] = sqrt(dMSE*(1.0/pNPTS[i] + 1.0/pNPTS[j]));
			switch (nIntervalType)
			{
				case 1:
					//Bonferroni
					pSig[nIndex]	 = 2*nMeanDiffSize*(1 - invt(fabs(pMeanDiffs[nIndex]/pSe[nIndex]), nDOF));
					if (pSig[nIndex] > 1)
						pSig[nIndex] = 1;
					else if (pSig[nIndex] < 0)
						pSig[nIndex] = 0;
					break;
				case 2:
					//Turkey
					pSig[nIndex]	= 1 - ocmath_inverse_studentized_range(sqrt(2)*fabs(pMeanDiffs[nIndex]/pSe[nIndex]), nDOF, nGroupSize);
					if (pSig[nIndex] > 1)
						pSig[nIndex] = 1;
					else if (pSig[nIndex] < 0)
						pSig[nIndex] = 0;
					break;
				case 3:
					//Scheffe
					pSig[nIndex] =1 - invf(pow(pMeanDiffs[nIndex]/pSe[nIndex],2.0)/(nGroupSize-1), nGroupSize-1,nDOF);
					if (pSig[nIndex] > 1)
						pSig[nIndex] = 1;
					else if (pSig[nIndex] < 0)
						pSig[nIndex] = 0;
					break;
				default:
					break;
			}
			nIndex++;
		}
	}

#ifdef		NAG_H
	Nag_IntervalType NagType;
	switch (nIntervalType)
	{
		case 1:
			NagType = Nag_BonferroniInterval;
			break;
		case 2:
			NagType = Nag_TukeyInterval;
			break;
		case 3:
			NagType = Nag_ScheffeInterval;
			break;
		default:
			return STATS_ERROR_SETTING;	//input type error
			break;
	}

	//double *pC = 0;
	int nSizeC = nGroupSize * nGroupSize;
	vector pC;
	pC.SetSize(nSizeC);
	//pC=(double*)malloc(nSizeC*sizeof(double));

	iErr = ANOVA_Compute_StdErr_MeanDiffs(nGroupSize, dMSE, pNPTS, nSizeNPTs, pC, nSizeC );
	vector<long> piSig;
	piSig.SetSize(nMeanDiffSize);
	//piSig = (long*)malloc(nMeanDiffSize*sizeof(long));
	
	double dDOF = (double)nDOF;
	
	NagError NagFail;
	nag_anova_confid_interval( NagType, nGroupSize, pMean, dDOF, pC, nGroupSize, dAlpha, pLCIs,
			 pUCIs, piSig , &NagFail);

	//free(pC);
	//free(piSig);

	iErr = NagFail.code;
	if (iErr != NE_NOERROR)
		return iErr;

	return STATS_NO_ERROR;
#else		// !NAG_H
	return STATS_FALSE;
#endif		// !NAG_H
}
*/

int stats_anova_mean_comparison(const int nIntervalType, 
								const ANOVADescStats *DescStats, 
								const RegANOVA *ANOVATable,
								const double dAlpha,
								const int iLevels,
								MeanCompStats *pstMeanCompstats, uint nSizeMeanComp) 

{
	vector vMean;
	vector<int> vNPTs;
	vMean.SetSize(iLevels);
	vNPTs.SetSize(iLevels);	
	double dMeanSE;
	int iDOF;
	vector<int> vIndexComparison;
	vector<long> vSig;
///Kevin 09/01/05 CHANGE_STATISTICS_NAME
	//vector vMeanDiffs, vLCIs, vUCIs, vSE, vStat, vProb, vAlpha;
	vector vMeanDiffs, vLCLs, vUCLs, vSE, vStat, vProb, vAlpha;
///End CHANGE_STATISTICS_NAME
	int nSize = iLevels * (iLevels - 1)/2;
	vIndexComparison.SetSize(nSize);
	vMeanDiffs.SetSize(nSize);
///Kevin 09/01/05 CHANGE_STATISTICS_NAME
	//vLCIs.SetSize(nSize);
	//vUCIs.SetSize(nSize);
	vLCLs.SetSize(nSize);
	vUCLs.SetSize(nSize);
///End CHANGE_STATISTICS_NAME
	vSE.SetSize(nSize);
	vStat.SetSize(nSize);
	vProb.SetSize(nSize);
	vSig.SetSize(nSize);
	vAlpha.SetSize(nSize);
	
	if (nSizeMeanComp < nSize)
		return STATS_PARAMS_ARRAY_TOO_SMALL;
	
	for(int k = 0; k < iLevels; k++)
	{
		vNPTs[k] = DescStats[k].N;
		vMean[k] = DescStats[k].Mean;
	}	
	
	iDOF = ANOVATable->AnovaTable[ANOVA_ERR].DOF;
	dMeanSE = ANOVATable->AnovaTable[ANOVA_ERR].MeanSq;
	
	int nRet = ocmath_anova_mean_comparison(vMean, iLevels,
											vNPTs, iLevels,
											dMeanSE,
											iDOF,
											dAlpha,
											iLevels,
											vIndexComparison, nSize,
											vMeanDiffs, nSize,
										///Kevin 09/01/05 CHANGE_STATISTICS_NAME
											//vLCIs, nSize,
											//vUCIs, nSize,
											vLCLs, nSize,
											vUCLs, nSize,
										///End CHANGE_STATISTICS_NAME
											vSE, nSize,
											vStat, nSize,
											vProb, nSize,
											vSig, nSize,
											vAlpha, nSize,
											nIntervalType);
	for(int ii = 0; ii< nSize; ii++)
		{
			pstMeanCompstats[ii].IndexComp = vIndexComparison[ii];
			pstMeanCompstats[ii].MeanDiff  = vMeanDiffs[ii];
			pstMeanCompstats[ii].SEM 		= vSE[ii];
			pstMeanCompstats[ii].Stat 		= vStat[ii];
			pstMeanCompstats[ii].Prob 		= vProb[ii];
			pstMeanCompstats[ii].Sig 		= vSig[ii];
			pstMeanCompstats[ii].Alpha		= vAlpha[ii];
		///Kevin 09/01/05 CHANGE_STATISTICS_NAME
			//pstMeanCompstats[ii].LCI 		= vLCIs[ii];
			//pstMeanCompstats[ii].UCI 		= vUCIs[ii];			
			pstMeanCompstats[ii].LCL 		= vLCLs[ii];
			pstMeanCompstats[ii].UCL 		= vUCLs[ii];
		///End CHANGE_STATISTICS_NAME
		}									
	return nRet;											
}
/*
int stats_anova_mean_comparison(const int nIntervalType, 
								const ANOVADescStats *DescStats, 
								const RegANOVA *ANOVATable,
								const double dAlpha,
								const int iLevels,
								MeanCompStats *MeanCompstats, uint nSizeMeanComp) 
{
	vector vMean;
	vector<int> vNPTs;
	vMean.SetSize(iLevels);
	vNPTs.SetSize(iLevels);	
	double dMeanSE;
	int iDOF;
	vector vMeanDiffs, vLCIs, vUCIs, vSE, vSig;
	int nSize = iLevels * (iLevels - 1);
	vMeanDiffs.SetSize(nSize);
	vLCIs.SetSize(nSize);
	vUCIs.SetSize(nSize);
	vSE.SetSize(nSize);
	vSig.SetSize(nSize);
	
	if (nSizeMeanComp < nSize)
		return STATS_PARAMS_ARRAY_TOO_SMALL;
	
	for(int k = 0; k < iLevels; k++)
	{
		vNPTs[k] = DescStats[k].N;
		vMean[k] = DescStats[k].Mean;
	}	
	
	iDOF = ANOVATable->AnovaTable[ANOVA_ERR].DOF;
	dMeanSE = ANOVATable->AnovaTable[ANOVA_ERR].MeanSq;

	int nRet = anova_mean_comparison(vMean,iLevels,
											vNPTs,iLevels,
											dMeanSE,
											iDOF,
											dAlpha,
											iLevels,
											vMeanDiffs,nSize,
											vLCIs, nSize,
											vUCIs,nSize,
											vSE,nSize,
											vSig,nSize,
											nIntervalType);	
	for(int ii = 0; ii< nSize; ii++)
		{
			MeanCompstats[ii].MeanDiff  = vMeanDiffs[ii];
			MeanCompstats[ii].SEM 		= vSE[ii];
			MeanCompstats[ii].Prob 		= vSig[ii];
			MeanCompstats[ii].LCI 		= vLCIs[ii];
			MeanCompstats[ii].UCI 		= vUCIs[ii];
			
			
		}									
	return nRet;											
}
*/
/// END REWRITE_ONEWAYANOVA_MEAN_COMPARISON

/**
		Function to compute actual or hypothetical powers for one way ANOVA. 
*/
int stats_anova_powers(const double dSigLevel, const vector& vSampleSizes, 
	 RegANOVA* ANOVATable, vector& vPowers)
{
	UINT nSizes = vSampleSizes.GetSize();
	vPowers.SetSize(nSizes);
	return ocmath_anova_powers(dSigLevel, vSampleSizes,nSizes,ANOVATable,vPowers);
}



/////////////////////////////////////////////////////////////////////////////////
//
//				Independent Sample Sizes and Powers Functions									
//
/////////////////////////////////////////////////////////////////////////////////

/**
		Function to compute hypothetical powers for a one sample t-Test.
*/
int stats_tTest_powers1(const double dSigLevel, const int iTailType, const vector& vSampleSizes,
						PowerOpts& tPowerOpts, vector& vPowers)
{
    int nSizes = vSampleSizes.GetSize();
	vPowers.SetSize(nSizes);
	return ocmath_tTest_powers1(dSigLevel,iTailType,vSampleSizes,nSizes,&tPowerOpts, vPowers);	
}


/**
		Function to computer hypothetical sample sizes for a one sample t-Test.
*/
int stats_tTest_sample_size(const double dSigLevel, const int iTailType, const vector& vPowers, 
						   PowerOpts& tPowerOpts, vector& vSampleSizes, double& dLowestPow)
{
	dLowestPow = NANUM;
	int nSizes = vPowers.GetSize();
	vSampleSizes.SetSize(nSizes);
	int nRet = ocmath_tTest_samplesize(dSigLevel,iTailType,vPowers,nSizes,&tPowerOpts,vSampleSizes, &dLowestPow);
	if (nRet != STATS_NO_ERROR)
		return nRet;
	
	return STATS_NO_ERROR;
}


/**
		Function to compute hypothetical powers for a two sample t-Test.
*/
int stats_two_sample_tTest_powers1(const double dSigLevel, const HypotTestOptions& opt, const vector& vSampleSizes,
                                   PowerOpts& tPowerOpts, vector& vPowers)
{
	int nSizes = vSampleSizes.GetSize();
	vPowers.SetSize(nSizes);
	return ocmath_two_sample_tTest_powers1(dSigLevel,&opt,vSampleSizes,nSizes,&tPowerOpts,vPowers);
}

/**
		Function to compute hypothetical sample size for a two sample t-Test.
*/
int stats_two_sample_tTest_sample_size(const double dSigLevel, const HypotTestOptions& opt,
     									   const vector& vPowers, 
	 									   PowerOpts& tPowerOpts,
	 									   vector& vSampleSizes, double& dLowestPow)
{
	dLowestPow = NANUM;
	int nSizes = vPowers.GetSize();
	vSampleSizes.SetSize(nSizes);
	return ocmath_two_sample_tTest_samplesize(dSigLevel,&opt,vPowers,nSizes,&tPowerOpts,vSampleSizes,&dLowestPow);
}


/**
		Function to compute ypothetical powers for a pair t-Test.
*/
int stats_two_sample_tTest_powers_pair1(const double dSigLevel, const HypotTestOptions& opt, const vector& vSampleSizes,
	                                    PowerOpts& tPowerOpts, vector& vPowers)
{
	int nSizes = vSampleSizes.GetSize();
	vPowers.SetSize(nSizes);
	return ocmath_two_sample_tTest_powers_pair1(dSigLevel,&opt,vSampleSizes,nSizes,&tPowerOpts,vPowers);
}


/**
		Function to compute hypothetical Sample Size for a pair t-Test.
*/
int stats_two_sample_tTest_sample_size_pair(const double dSigLevel, const HypotTestOptions& opt,
     									   const vector& vPowers, 
	 									   PowerOpts& tPowerOpts,
	 									   vector& vSampleSizes, double& dLowestPow)
{
	dLowestPow = NANUM;
	int nSizes = vPowers.GetSize();
	vSampleSizes.SetSize(nSizes);
	return ocmath_two_sample_tTest_samplesize_pair(dSigLevel,&opt,vPowers,nSizes,&tPowerOpts,vSampleSizes,&dLowestPow);
}

/**
		Function to compute hypothetical powers for a one way ANOVA.
*/
int stats_anova_powers1(const double dSigLevel, const vector& vSampleSizes, 
	ANOVAPowerOpts& tPowerOpts, vector& vPowers)
{
	UINT nSizes = vSampleSizes.GetSize();
	vPowers.SetSize(nSizes);
	return ocmath_anova_powers1(dSigLevel, vSampleSizes,nSizes,&tPowerOpts,vPowers);
}


/**
		Function to compute hypothetical sample size for a one way ANOVA.
*/
int stats_anova_sample_size(const double dSigLevel, const vector& vPowers, 
	ANOVAPowerOpts& tPowerOpts, vector& vSampleSizes)
{
	UINT nSizes = vPowers.GetSize();
	vSampleSizes.SetSize(nSizes);	
    return ocmath_anova_samplesize(dSigLevel, vPowers,nSizes,&tPowerOpts,vSampleSizes);
}

/**
		Function to put Powers to a structure
*/	
void PowerResultsNANUM(PowerResults& PowerResult)
{
	PowerResult.Alpha = NANUM;
	PowerResult.SampleSize = NANUM;
	PowerResult.Power = NANUM;
}
int stats_power_output(const vector& vPowers, const vector& vSampleSize, const double dSigLevel, PowerResults* stPowers, const int nPowerResults)	
{
	/// *** Initialize the output
	for (int ii = 0; ii < nPowerResults; ii++)
	{
		PowerResultsNANUM(stPowers[ii])
	}
	
	/// *** Return error if error setting
	int nPowers = vPowers.GetSize();
	int nSampSizes = vSampleSize.GetSize();
	if (nPowers != nSampSizes)
		return STATS_ERROR_DATA_DIFF_SIZE;
	
	if (nPowerResults < nSampSizes)
		return STATS_PARAMS_ARRAY_TOO_SMALL;
	
	for (ii = 0; ii < nPowers; ii++)
	{
		stPowers[ii].Alpha = dSigLevel;
		stPowers[ii].SampleSize = vSampleSize[ii];
		stPowers[ii].Power = vPowers[ii];
		
	}
	return STATS_NO_ERROR;
}

/**
		Function to put Sample Sizes to structure
*/	
void SampleSizeResNANUM(SAMPResults& SampleSizeRes)
{
	SampleSizeRes.Alpha = NANUM;
	SampleSizeRes.SampleSize = NANUM;
	SampleSizeRes.Power = NANUM;
}
int stats_samplesize_output(const vector& vPowers, const vector& vSampleSize, const double dSigLevel, SAMPResults* psSampleSize, const int nSampleSizeRes)
{
	/// *** Initialize the output
	for (int ii = 0; ii < nSampleSizeRes; ii++)
	{
		SampleSizeResNANUM(psSampleSize[ii])
	}
	
	/// *** Return error if error setting
	int nPowers = vPowers.GetSize();
	int nSampSizes = vSampleSize.GetSize();
	if (nPowers != nSampSizes)
		return STATS_ERROR_DATA_DIFF_SIZE;
	
	if (nSampleSizeRes < nSampSizes)
		return STATS_PARAMS_ARRAY_TOO_SMALL;
	
	for (ii = 0; ii < nPowers; ii++)
	{
		psSampleSize[ii].Alpha = dSigLevel;
		psSampleSize[ii].SampleSize = vSampleSize[ii];
		psSampleSize[ii].Power = vPowers[ii];
		
	}
	return STATS_NO_ERROR;
}	


/////////////////////////////////////////////////////////////////////////////////
//
//				Two Way ANOVA Functions									
//
/////////////////////////////////////////////////////////////////////////////////

/**
		Two-Way ANOVA function to compute the Two-Way ANOVA F table.
*/
int stats_anova_two_way_compute_anova_table(int iInteractions, vector<int> &vFACTOR_A, int nLevelsA, vector<int> &vFACTOR_B,
	 int nLevelsB, vector<double> &vDEP_VAR, ANOVA2DescStats *ADescStats, ANOVA2DescStats *BDescStats, ANOVA2DescStats& ToDescStats, 
	 ANOVA2DescStats *ABDescStats, ANOVA2Table* arANOVA2_Table)
{
	int nPts = vDEP_VAR.GetSize();		
	int nRet = ocmath_stats_anova_two_way_compute_anova_table(  iInteractions, 
																vFACTOR_A, nLevelsA,
																vFACTOR_B, nLevelsB,
																vDEP_VAR, 
																nPts,
																ADescStats,
																BDescStats,
																&ToDescStats, 
																ABDescStats,
																arANOVA2_Table);
																
	return nRet;
}

/**
		Two-Way ANOVA function to perform a Means Comparison analysis
*/

/// Max 03/16/05 REWRITE_TWOWAYANOVA_MEAN_COMPARISON
int stats_anova_two_way_means_comparison(const int iTestType, const vector<int> &vFactor, 
										const int nLevels, const vector &vDEP_VAR,
										ANOVA2Table* arANOVA2_Table, double dSigLevel, 
	 									MeanCompStats *ANOVA2_MeanCompStats, const uint nSizeMeanComp)

{
	double dDFE, dMSE;
	dDFE = arANOVA2_Table[ANOVA2_ERR].DOF;
	dMSE = arANOVA2_Table[ANOVA2_ERR].MeanSq;
	vector<int> vIndexComparison;
	vector<long> vSig;
///Kevin 09/01/05 CHANGE_STATISTICS_NAME
	//vector	vMeanDiffs, vLCIs, vUCIs, vSe, vStat, vProb, vAlpha;	
	vector	vMeanDiffs, vLCLs, vUCLs, vSe, vStat, vProb, vAlpha;	
///End CHANGE_STATISTICS_NAME
	int nSize = nLevels * (nLevels - 1)/2;
	vIndexComparison.SetSize(nSize);
	vMeanDiffs.SetSize(nSize);
///Kevin 09/01/05 CHANGE_STATISTICS_NAME	
	//vLCIs.SetSize(nSize);
	//vUCIs.SetSize(nSize);
	vLCLs.SetSize(nSize);
	vUCLs.SetSize(nSize);
///End CHANGE_STATISTICS_NAME
	vSe.SetSize(nSize);
	vStat.SetSize(nSize);
	vProb.SetSize(nSize);
	vSig.SetSize(nSize);
	vAlpha.SetSize(nSize);

	int nPts = vDEP_VAR.GetSize();

	if (nSizeMeanComp < nSize)
		return STATS_PARAMS_ARRAY_TOO_SMALL;

	int nRet = osANOVA2Way_Means_Comparison( vFactor, 
											 nLevels,
											 vDEP_VAR,
											 nPts, 
											 dDFE, 
											 dMSE, 
											 dSigLevel,
											 vIndexComparison,nSize,
											 vMeanDiffs,nSize,
										///Kevin 09/01/05 CHANGE_STATISTICS_NAME
											 //vLCIs, nSize,
											 //vUCIs, nSize,
											 vLCLs, nSize,
											 vUCLs, nSize,
										 ///End CHANGE_STATISTICS_NAME
											 vSe, nSize,
											 vStat, nSize,
											 vProb, nSize,
											 vSig, nSize,
											 vAlpha, nSize,
											 iTestType);
	for(int ii = 0; ii< nSize; ii++)
	{
		ANOVA2_MeanCompStats[ii].IndexComp = vIndexComparison[ii];
		ANOVA2_MeanCompStats[ii].MeanDiff = vMeanDiffs[ii];
		ANOVA2_MeanCompStats[ii].SEM = vSe[ii];
		ANOVA2_MeanCompStats[ii].Stat = vStat[ii];
		ANOVA2_MeanCompStats[ii].Prob = vProb[ii];
		ANOVA2_MeanCompStats[ii].Sig = vSig[ii];
		ANOVA2_MeanCompStats[ii].Alpha = vAlpha[ii];		
	///Kevin 09/01/05 CHANGE_STATISTICS_NAME	
		//ANOVA2_MeanCompStats[ii].LCI = vLCIs[ii];
		//ANOVA2_MeanCompStats[ii].UCI = vUCIs[ii];			
		ANOVA2_MeanCompStats[ii].LCL = vLCLs[ii];
		ANOVA2_MeanCompStats[ii].UCL = vUCLs[ii];			
	///End CHANGE_STATISTICS_NAME
	}										
	return nRet;							
}
/*
int stats_anova_two_way_means_comparison(const int iTestType, const vector<int> &vFactor, 
										const int nLevels, const vector &vDEP_VAR,
										ANOVA2Table* arANOVA2_Table, double dSigLevel, 
	 									MeanCompStats *ANOVA2_MeanCompStats, const uint nSizeMeanComp)

{
	double dDFE, dMSE;
	dDFE = arANOVA2_Table[3].DOF;
	dMSE = arANOVA2_Table[3].MeanSq;
	vector	vMeanDiffs, vLCIs, vUCIs, vSe, vSig;	
	int nSize = nLevels * (nLevels - 1);
	vMeanDiffs.SetSize(nSize);
	vLCIs.SetSize(nSize);
	vUCIs.SetSize(nSize);
	vSe.SetSize(nSize);
	vSig.SetSize(nSize);
	int nPts = vDEP_VAR.GetSize();
	
 	///Echo 7/20/04 ADD_SIZE_FOR_POINTER
	if (nSizeMeanComp < nSize)
		return STATS_PARAMS_ARRAY_TOO_SMALL;
	///end ADD_SIZE_FOR_POINTER
	
	int nRet = osANOVA2Way_Means_Comparison( vFactor, 
											 nLevels,
											 vDEP_VAR,
											 nPts, 
											 dDFE, 
											 dMSE, 
											 dSigLevel,
											 vMeanDiffs,nSize,
											 vLCIs, nSize,
											 vUCIs, nSize,
											 vSe, nSize,
											 vSig, nSize,
											 iTestType);
	for(int ii = 0; ii< nLevels * (nLevels - 1); ii++)
	{
		ANOVA2_MeanCompStats[ii].MeanDiff = vMeanDiffs[ii];
		ANOVA2_MeanCompStats[ii].SEM = vSe[ii];
		ANOVA2_MeanCompStats[ii].Prob = vSig[ii];
		ANOVA2_MeanCompStats[ii].LCI = vLCIs[ii];
		ANOVA2_MeanCompStats[ii].UCI = vUCIs[ii];			
	}										
	return nRet;							
}	 						
*/
/// END REWRITE_TWOWAYANOVA_MEAN_COMPARISON

/**
		Two-Way ANOVA function used to compute and output Two-Way ANOVA Power Analysis.
*/
int stats_anova_two_way_power(int iInteractions, const double dSigLevel,
	  ANOVA2Table *arANOVA2_Table,vector &vSampleSizes, vector& vPowers, int iSource)
{
	int nSize = vSampleSizes.GetSize();
	/// Iris 05/10/2007 v8.0615 REMOVE_BRANCH_CHECK_BOX, for now support only out actual power
	/*
	if (iPower == 0)
		vPowers.SetSize(1);
	else if (iPower == 1)
		vPowers.SetSize(nSize + 1);
		//vPowers.SetSize(nSize);///Echo 12/10/04 ///Iris 12/28/04
	int nRet = ocmath_osANOVA2Way_Compute_Power( iPower,
												 iInteractions,
												 dSigLevel,
												 vSampleSizes,
												 nSize,
												 arANOVA2_Table,
												 vPowers,
												 iSource);	
	*/
	vPowers.SetSize(nSize);
	int nRet = ocmath_osANOVA2Way_Compute_Power( iInteractions,
												 dSigLevel,
												 vSampleSizes,
												 nSize,
												 arANOVA2_Table,
												 vPowers,
												 iSource);	
	///end REMOVE_BRANCH_CHECK_BOX
	return nRet;
}
	 


/////////////////////////////////////////////////////////////////////////////////
//
//				Survival Analysis Functions									
//
/////////////////////////////////////////////////////////////////////////////////
///Echo 11/30/2004 REWRITE_KAPLAN_USING_LLVC
/**
		Survival Analysis function used to compute the Kaplan-Meier Estimate of survival 
		probalities for a sample of failure times.
*/
int	stats_kaplan_meier_surviv_fn(const vector& cvTime, const vector& vCensor, const vector& vCensored, const double dConLevel,
	SAEstim *psSurvivEstim, uint nSizeSAEstm, SAQuantEst *psSAQuantEstim, uint nSizeQuanEstm, SASummary& stSASumStats)
{
	return ocmath_kaplan_meier_surviv_fn(cvTime, cvTime.GetSize(), vCensor, vCensor.GetSize(), vCensored, vCensored.GetSize(), dConLevel,
										psSurvivEstim, nSizeSAEstm, psSAQuantEstim, nSizeQuanEstm, &stSASumStats);
}
//static void SAEstimNANUM(SAEstim& stSAEstim)
//{
	//stSAEstim.Surv = NANUM;
	//stSAEstim.SD = NANUM;
	//stSAEstim.Events = -1;
	//stSAEstim.Remain = -1;	
//}
//
//static void SAQuantEstNANUM(SAQuantEst& stSAQuantEst)
//{
	//stSAQuantEst.Perc = NANUM;
	//stSAQuantEst.LCI = NANUM;
	//stSAQuantEst.UCI = NANUM;
	//
//}
//
//static  void SASummaryNANUM(SASummary& stSASummary)
//{
	//stSASummary.Total = -1;
	//stSASummary.Events = -1;
	//stSASummary.NumCen = -1;
	//stSASummary.PerCen = NANUM;
	//
//}
//
///**
		//List out the censored value in vector vCensor
//*/
//int stats_censored_value_list(const vector& vCensor, vector& vCensList)
//{
	//const int MaxCens = 20; //Max length of Censor List
	//int ii, jj;
	//int nSizeCens = vCensor.GetSize();
	//int nSizeList = vCensList.GetSize();
	//
	/////*** Terminate if error setting
	////Only list out the first 20 unduplicated censor values
	//if ((nSizeList > MaxCens) || (nSizeList <= 0))
		//return STATS_ERROR_INPUT_SIZE;
	//
	//if (nSizeCens <= 0)
		//return STATS_ERROR_TOO_FEW_DATA_PTS;
	//
	///// *** To devoid the changing of vData
	//vector vCensorCpy;
	//vCensorCpy = vCensor;
	//vCensorCpy.Sort();
	//
	///// *** Initialize the value of censor list
	//for (ii = 0; ii < MaxCens; ii++)
		//vCensList = NANUM;
	//
	///// *** Censor Value List	
	//vCensList[0] = vCensorCpy[0];
	//for (ii = 0; ii < MaxCens -1; ii++)
	//{
		//if ( NANUM == vCensList[ii] )
			//break;
		//for (jj = ii; jj < nSizeCens; jj++)
		//{
			//
			//if ((vCensorCpy[jj] > vCensList[ii])) 
			//{
				//vCensList[ii+1] = vCensorCpy[jj];
				//break;
			//}///end put value in vCensList
		//}//end loop in vCensor
	//}//end loop in vCensList
	//
	//return STATS_NO_ERROR;
	//
//}
//
///**
		//Survival Analysis function used to compute the Kaplan-Meier Estimate of survival 
		//probalities for a sample of failure times.
//*/
//
//int	stats_kaplan_meier_surviv_fn(const vector& cvTime, const vector& vCensor, const vector& vCensored, const double dConLevel,
	//SAEstim *psSurvivEstim, uint nSizeSAEstm, SAQuantEst *psSAQuantEstim, uint nSizeQuanEstm, SASummary& stSASumStats)
//{
	//int ii,jj,kk, nCensored;
	//int icen=1;
	//int nObs = cvTime.GetSize();
	//int nCol = 2;	//Input & Output are both of nObs*2.
	//vector vCensList(20);
	//nCensored = vCensored.GetSize();
//
	///// *** Terminate if error setting
	//if( nObs<2 )		//Must contain at least 2 data points.
		//return STATS_ERROR_TOO_FEW_DATA_PTS;
	//
	//if( nObs != vCensor.GetSize() )	//Size of vTime & vCensor should be equal
		//return STATS_ERROR_DATA_DIFF_SIZE;
	//
	//if( !( (dConLevel>0)&&(dConLevel<1) ) )	//Value for confidence level should >0 & <1
		//return STATS_ERROR_INVALID_CONF_LEVEL;
	//
	//if (nCensored > 20) //Number of censored value needed to be less than 20
		//return STATS_ERROR_INVALID_CENSOR;
	//
	//if (nSizeSAEstm < nObs)//Size of psSurvivEstim needs to be more than nObs
		//return STATS_PARAMS_ARRAY_TOO_SMALL;
	//
	//if (nSizeQuanEstm < 3)//Size of psSAQuantEstim needs to be more than nObs
		//return STATS_PARAMS_ARRAY_TOO_SMALL;
	//
	//// Values of vCnesored need to be in vCensor
	//int nRet = stats_censored_value_list(vCensor, vCensList);
	//if (nRet != STATS_NO_ERROR)
		//return nRet;
//
	////--- CPY 10/26/04 CONST_CHECK_ENABLE_NEED_FIXING_CODES
	//vector vTime;
	//vTime = cvTime;
	////----
	//
	//for (ii = 0; ii < nCensored; ii++)
	//{
		//nRet = STATS_ERROR_INVALID_CENSOR;
		//for (jj = 0; jj < 20; jj++)
		//{
			//if ( NANUM == vCensList[jj] )
				//break;
			//if (vCensored[ii] == vCensList[jj] )
			//{
				//nRet = STATS_NO_ERROR;
				//continue;
			//}
		//}
		//if (nRet != STATS_NO_ERROR)
			//return nRet;
	//}
//
	///// *** Initialize parameters
	//for (ii = 0;  ii < nSizeSAEstm; ii++)
		//SAEstimNANUM(psSurvivEstim[ii]);	
		//
	//for (ii = 0;  ii < nSizeQuanEstm; ii++)
		//SAQuantEstNANUM(psSAQuantEstim[ii]);
	//
	//SASummaryNANUM(stSASumStats);
	//
	///// *** Set Censor values
	//vector<int> vCensorVal;
	//vCensorVal.SetSize(nObs);
	//for(ii=0;ii<nObs;ii++)
	//{
		//for (jj = 0; jj < nCensored; jj++)
		//{
			//if( is_equal(vCensor[ii],vCensored[jj]) )
			//{
				//vCensorVal[ii] = 1;
				//break;
			//}
			//else
				//vCensorVal[ii] = 0;
		//}
	//}//end loop
			//
	///// *** Find the number of Censors.
	//int nCensorNum;
	//vCensorVal.Sum(nCensorNum);
	//int nFail = nObs - nCensorNum;
//
	//
	///// *** Kaplan-Meier Estimate
	//// Contains distinct ordered failure times on output from NAG (does not include duplicates or censored times)
	//vector<double> vTP;
	//vTP.SetSize( nObs );
//
	//// Contains Kaplan-Meier estimate of Survival Probability for each time in vTP on output from NAG 
	//vector<double> vP;
	//vP.SetSize( nObs );
//
	//// Contains Kaplan-Meier estimate of error (standard deviation) of each probability in vP on output from NAG
	//vector<double> vPSIG;
	//vPSIG.SetSize( nObs );
//
	//// From <NAG\OCN_g12.h>: g12aac nag_prod_limit_surviv_fn: Compute Kaplan-Meier estimate of
	//// Survivorship Function. Frequency is NULL (count each point once; nd is the number of distinct
	//// failure times output in vTP.	
	//int nd;
	//NagError fail;
	//nag_prod_limit_surviv_fn( nObs, vTime, vCensorVal, NULL, &nd, vTP, vP, vPSIG , &fail);
	//nRet = fail.code;
//
	//if( nRet != NE_NOERROR )
		//return nRet;
	//
	//nCensorNum = nObs - nFail;
	//
	///// *** Sort the columns (Time & Censor) first based on Censor value, then on Time value;
	//vector<uint> vIdx;
	//vIdx.SetSize(nObs);
	////vector vTemp;
	////vTemp.SetSize(nObs);
	//vTime.Sort(SORT_ASCENDING,true,vIdx);
	//vCensorVal.Reorder(vIdx);
	//vCensorVal.Sort(SORT_ASCENDING,true,vIdx);
	//vTime.Reorder(vIdx);
	////vTime.Sort(SORT_ASCENDING,true,vIdx);
	////for(ii=0,jj=0,kk=nFail;ii<nObs;ii++)
	////{
		////if(vCensorVal[vIdx[ii]]==0)
			////vTemp[jj++] = vTime[ii];
		////else
			////vTemp[kk++] = vTime[ii];
	////}
	////vTime = vTemp;
	////
	////if(nFail>0)
	////{
		////vCensorVal.SetSize(nFail);
		////vCensorVal = 0;
		////vTemp.SetSize(nCensorNum);
		////vTemp = 1;
		////vCensorVal.Append(vTemp);
	////}
	////else	//all are censor value
	////{
		////vCensorVal.SetSize(nCensorNum);
		////vTemp = 1;
	////}
	//
	//
	////Quartile Estimates
	//vector vPerc;
	//vPerc.SetSize(3);
	//vPerc = -1;
	//double dP, dPofii;
	//for(ii=0, jj=0; (ii<nFail) && (jj<3); ii++)
	//{
		//dP = 0.75 - 0.25 * jj;	// Compute integerized threshold for next percentile
		//dPofii = vP[ii];
		//if( round(dP,6) >= round(dPofii,6) )	// If survivorship function value is less than or equal to percentile threshold...
		//{
			//vPerc[jj] = vTP[ii];
			//ii--; // Decrement loop index to start looking for next percentile with current time value
			//jj++; // Got current percentile...now get next one
		//}
	//}
	//
	//vector vUpLimit, vLowLimit;
	//vUpLimit.SetSize(3);
	//vLowLimit.SetSize(3);
	//vUpLimit = -1;
	//vLowLimit = -1;
	//double dChi = ChisqDist(dConLevel, 1);
	//for(ii=0; ii<3; ii++)
	//{
		//// Compute non-integerized threshold for current percentile
		//dP = 0.75 - ii * 0.25;
		//// Loop through all distict failure times
		//for(jj= 0; jj<nFail; jj++)
		//{
			//// If condition is TRUE then time value may be upper or lower limit 
			//if( fabs( vP[jj] - dP ) <= sqrt( dChi * pow( vPSIG[jj],2 ) ) )
			//{
				//if( vLowLimit[ii] < 0 ) // Lower confidence limit is the first time the condition is TRUE
					//vLowLimit[ii] = vTP[jj];
				//if( jj + 1 < nFail ) // Upper confidence limit is the next time value the last time the condition is TRUE 
					//vUpLimit[ii] = vTP[jj + 1];				
			//}
		//}
	//}
	//
	//for(ii=0; ii<3; ii++)
	//{
		//if(vPerc[ii]>=0)
			//psSAQuantEstim[ii].Perc = vPerc[ii];
		//else
			//psSAQuantEstim[ii].Perc = NANUM;	//if vPerc[ii]=-1, it is NANUM.
		//if(vLowLimit[ii]>=0)
			//psSAQuantEstim[ii].LCI = vLowLimit[ii];
		//else
			//psSAQuantEstim[ii].LCI = NANUM;		//if vLowLimit[ii]=-1, it is NANUM.
		//if(vUpLimit[ii]>=0)
			//psSAQuantEstim[ii].UCI = vUpLimit[ii];
		//else
			//psSAQuantEstim[ii].UCI = NANUM;		//if vUpLimit[ii]=-1, it is NANUM.
	//}
	//
	////Survivorship Function
	//vTime.Sort(SORT_ASCENDING,true,vIdx);	//Resort the columns based on Time value.
	//
	////Check: if the Censor value is always behind the Fail value for the same Time point.
	//uint uiTemp;
	//for(ii=0;ii<nObs-1;ii++)
	//{
		//if(vTime[ii]==vTime[ii+1])
		//{
			//if(vIdx[ii]>vIdx[ii+1])
			//{
				//uiTemp = vIdx[ii];
				//vIdx[ii] = vIdx[ii+1];
				//vIdx[ii+1] = uiTemp;
			//}
		//}
	//}
	//
	//// Loop on all non-zero time values
	//int nEvents=0, nRemain=nObs;
	//int dPreUncensoredTime = -1;
	//jj = -1;
	//for(ii=0;ii<nObs;ii++)
	//{
		//psSurvivEstim[ii].Time = vTime[ii];
		//if(vCensorVal[vIdx[ii]]==1)// If censored value...
		//{
			////Output missing values
			//psSurvivEstim[ii].Surv = NANUM;
			//psSurvivEstim[ii].SD = NANUM;
			//psSurvivEstim[ii].Events = nEvents;
		//}
		//else// Else, not censored
		//{
			//// Output values of psSurvivEstim structure
			//// Must index into vP and vPSIG differently than vTIME because only non-repeated and uncensored
			//// survival/error values are in vP/vPSIG while all times are in vTIME 
			//if (jj < 0 || vTime[ii] != dPreUncensoredTime)
				//jj++;
			//dPreUncensoredTime = vTime[ii];
			//psSurvivEstim[ii].Surv = vP[jj];
			//if(vP[jj] == 0)
				//psSurvivEstim[ii].SD = 0;
			//else
				//psSurvivEstim[ii].SD = vPSIG[jj];
			//psSurvivEstim[ii].Events = ++nEvents;
		//}
		//psSurvivEstim[ii].Remain = --nRemain;
	//}
	//
	////Summary of Event and Censored Values
	//stSASumStats.Total = nObs;
	//stSASumStats.Events = nFail;
	//stSASumStats.NumCen = nCensorNum;
	//stSASumStats.PerCen = (double)nCensorNum/nObs;
	//
	//return STATS_NO_ERROR;
	//
//}



/** 
		Compute the Parameter estimates in Cox model and Weibull model
	Parameter:
		nCoeff: Input number of coefficients. 
		Coef: Input pointer to parameter estimates and associated statistics computed by nag
		SE: Input pointer to asymptotic standard error of the parameter estimates.
		psSAParaEstim:Output pointer to structure of parameter estimates. 
	Return:
		Returns STATS_NO_ERROR on successful exit 
*/
//static int SA_Compute_Param_Estim(int nCoeff, double* Coef, double* SE, SAParaEstim *psSAParaEstim)
//{
	//int ii;
	//for (ii = 0; ii < nCoeff; ii++)	
	//{
		//psSAParaEstim[ii].DOF = 1;
		//psSAParaEstim[ii].ParaEst = Coef[ii];
		//psSAParaEstim[ii].SE = SE[ii];
		//psSAParaEstim[ii].ChiSq = pow(psSAParaEstim[ii].ParaEst/psSAParaEstim[ii].SE, 2);//ChiSq = (Coeff/SE)^2
		//psSAParaEstim[ii].Prob = 1 - invchisq(psSAParaEstim[ii].ChiSq, psSAParaEstim[ii].DOF); //Prob = 1 - invchisq(ChiSq)
		//psSAParaEstim[ii].Hazard = exp(psSAParaEstim[ii].ParaEst);//Hazard = exp(Coeff)
	//}
	//return STATS_NO_ERROR;
	//
//}
//
///**
		//Survival Analysis function used to compute parameter estimates and other statistics associated with 
		//the Cox Proportional Hazards Model for fixed covariates.
//*/
//static void SAFitStatsNANUM(SAFitStats& stSAFitStats)
//{
	//stSAFitStats.Criter = NANUM;
	//stSAFitStats.ChiSq = NANUM;
	//stSAFitStats.DOF = -1;
	//stSAFitStats.Prob = NANUM;	
//}
//
//static void SAParaEstimNANUM(SAParaEstim& stSAParaEstim)
//{
	//stSAParaEstim.DOF = 0;
	//stSAParaEstim.ParaEst = 0;
	//stSAParaEstim.SE = NANUM;
	//stSAParaEstim.ChiSq = NANUM;
	//stSAParaEstim.Prob = NANUM;
	//stSAParaEstim.Hazard = NANUM;	
//}
//
int	stats_surviv_cox_model(const vector& vTime, const vector& vCensor, const vector& vCensored, const matrix& mCovars,
	SASummary& stSASummary, SAFitStats &stSAFitStats, SAParaEstim *psSAParaEstim, uint nSizeParaEstm, SAEstim *psSurvivEstim, uint nSizeSAEstm, matrix& mCov, matrix& mCorr)
{
	double* pmCov = NULL;
	double* pmCorr = NULL;
	if (mCov)
	{
		mCov.SetSize(mCovars.GetNumCols(), mCovars.GetNumCols());
		pmCov = mCov;
	}
	if(mCorr)
	{
		mCorr.SetSize(mCovars.GetNumCols(), mCovars.GetNumCols());
		pmCorr = mCorr;
	}
	return ocmath_surviv_cox_model(vTime, vTime.GetSize(), vCensor, vCensor.GetSize(), vCensored, vCensored.GetSize(), mCovars, mCovars.GetNumCols(), mCovars.GetNumRows(), 
	                               &stSASummary, &stSAFitStats, psSAParaEstim, nSizeParaEstm, psSurvivEstim, nSizeSAEstm, pmCov, pmCorr, mCovars.GetNumCols(), mCovars.GetNumCols());
}
//{
	//int ii,jj,kk;
	//int icen=1;
	//int nObs = vTime.GetSize();
	//int nEF = mCovars.GetNumCols(); 
	//vector vCoef(nEF);
	//int nCoef;
	//double dLogLik;
	//int nRet;
	//vector vCensList(20);
	//int nCensored = vCensored.GetSize();
		//
	///// *** Terminate if error setting
	////At least one covariate
	//if (nEF < 1) 
		//return STATS_ERROR_SETTING;
	//
	////Data points must be larger than 2
	//if (nObs < 2) 
		//return STATS_ERROR_TOO_FEW_DATA_PTS;
	//
 	////Number of censored value needed to be less than 20	
	//if (nCensored > 20)
		//return STATS_ERROR_INVALID_CENSOR;
	//
	////Size of vTime & vCensor should be equal
	//if( nObs != vCensor.GetSize() )
		//return STATS_ERROR_DATA_DIFF_SIZE;
	//
	////Size of Time, Censor and Covariates must be the same	
	//int nTime = vTime.GetSize();
	//int nCensor = vCensor.GetSize();
	//int nCovarRow = mCovars.GetNumRows();
	//if (nTime != nCensor || nCensor != nCovarRow || nTime != nCovarRow) 
		//return STATS_ERROR_DATA_DIFF_SIZE;
	//
	////size of psSAParaEstim needs to be larger than nEF	
	//if (nSizeParaEstm < nEF) 
		//return STATS_PARAMS_ARRAY_TOO_SMALL;
	//
	//// Values of pCnesor need to be in vCensor
	//nRet = stats_censored_value_list(vCensor, vCensList);
	//if (nRet != STATS_NO_ERROR)
		//return nRet;
	//
	//for (ii = 0; ii < nCensored; ii++)
	//{
		//nRet = STATS_ERROR_INVALID_CENSOR;
		//for (jj = 0; jj < 20; jj++)
		//{
			//if ( NANUM == vCensList[jj] )
				//break;
			//if (vCensored[ii] == vCensList[jj] )
			//{
				//nRet = STATS_NO_ERROR;
				//continue;
			//}
		//}
		//if (nRet != STATS_NO_ERROR)
			//return nRet;
	//}
	//
	///// *** Initialize parameters
	//SASummaryNANUM(stSASummary);
//
	//SAFitStatsNANUM(stSAFitStats);
		//
	//for (ii = 0;  ii < nSizeParaEstm; ii++)
		//SAParaEstimNANUM(psSAParaEstim[ii]);
		//
	//for (ii = 0;  ii < nSizeSAEstm; ii++)
		//SAEstimNANUM(psSurvivEstim[ii]);	
		//
	///// *** Set Censor values
	//vector<int> vCensorVal;
	//vCensorVal.SetSize(nObs);
	//for(ii=0;ii<nObs;ii++)
	//{
		//for (jj = 0; jj < nCensored; jj++)
		//{
			//if( is_equal(vCensor[ii],vCensored[jj]) )
			//{
				//vCensorVal[ii] = 1;
				//break;
			//}
			//else
				//vCensorVal[ii] = 0;
		//}
	//}
					//
	///// *** Cox Hazard Regression		
	//// vSZ[jj]==1 tells NAG that jjth Covariate is to be used in model...all are used so all == 1
	//vector<int> vSZ;
	//vSZ.SetSize( nEF );
	//vSZ = 1;
	//
	//// Output failure times...not used by CoxPHM but NAG requires
	//vector<double> vTP;
	//vTP.SetSize( nObs );
//
	//// Output Survivorship function...NAG expects matrix but second dimension is 1 so vector is used
	//vector<double> vSUR;
	//vSUR.SetSize( nObs );
//
	//// vSE[jj] is asymptotic standard error of jjth parameter estimate on exit from NAG 
	//vector<double> vSE;
	//vSE.SetSize( nEF );
	//vSE = 0;
	//
	//// Score Function not used by CoxPHM but NAG requires
	//vector<double> vSC;
	//vSC.SetSize( nEF );
	//
	//// Covariate matrix not used by CoxPHM but NAG requires
	//vector<double> vCOV;
	//vCOV.SetSize( nEF*( nEF + 1 )/2 );
//
	//// Residuals not used by CoxPHM but NAG requires
	//vector<double> vRES;
	//vRES.SetSize( nObs );
	//
	////int nFail;
	//// From <NAG\OCN_g12.h>: g12bac nag_surviv_cox_model: Compute estimates of parameters of Cox model
	////  and Survivorship function values. CoxPHM does not use strata (ns=0), Offsets and Stratum are not used
	////  and ==NULL, second dimension of vSUR is 1 (tdsur==1), do not print (iprint=0), and do not use Outfile
	////  which also ==NULL. 
	//NagError fail;
	//nag_surviv_cox_model( nObs, nEF, 0, mCovars, nEF, vSZ, nEF, vTime, vCensorVal, NULL, NULL, &dLogLik,
	//vCoef, vSE, vSC, vCOV, vRES, &nCoef, vTP, vSUR, 1, nObs, COX_MODEL_TOLERANCE, COX_MODEL_ITERATIONS, 0, NULL , &fail);
	//
	//nRet = fail.code;
	//
	//if( nRet != STATS_NO_ERROR && nRet != 58)
		//return nRet;
		//
	///// *** Compute the Parameter Estimates	
	//double nCensSum;
	//vCensorVal.Sum(nCensSum);
	//
	//if (nObs == nCensSum ) //if all censor values are the same
	//{
		//for (ii = 0;  ii < nSizeParaEstm; ii++)
			//SAParaEstimNANUM(psSAParaEstim[ii]);						
	//}else
		//nRet = SA_Compute_Param_Estim(nEF, vCoef, vSE, psSAParaEstim);
								//
	//// *** Compute the Fit Statistics
	//stSAFitStats.Criter = dLogLik; //-2 log likelihood
	//if (nObs == nCensSum || 0 == nCensSum) //if all censor values are the same
		//stSAFitStats.DOF = 0;
	//else
		//stSAFitStats.DOF = nEF;//Degree of freedom
//
	///// *** Compute the Chi-Square value of fit
	////The Chi-Square value of fit can be computed from the log-likelihood and 
	////reduced log-likelihood based upon no iterations. 	
	////Computer reduced -2 log likelihood
	//double dInLogLik;
	//int nInCoef;
	//vector vInCoef(nEF);
	//vInCoef = 0;
	//
	//nag_surviv_cox_model( nObs, nEF, 0, mCovars, nEF, vSZ, nEF, vTime, vCensorVal, NULL, NULL, &dInLogLik,
	//vInCoef, vSE, vSC, vCOV, vRES, &nInCoef, vTP, vSUR, 1, nObs, COX_MODEL_TOLERANCE, 0, 0, NULL , &fail);
		 	 	//
	//nRet = fail.code;
	//
	//stSAFitStats.ChiSq = dInLogLik - dLogLik; //Compute Chi-Square values of Fit
	//if (nObs == nCensSum || 0 == nCensSum) //if all censor values are the same
		//stSAFitStats.Prob = 0;
	//else
		//stSAFitStats.Prob = 1 - invchisq(stSAFitStats.ChiSq, stSAFitStats.DOF); //Compute P(Chi-Square) values
	//
	//
	///// *** Compute psSurvivEstim and stSASummary
	//double dConLevel = 0.95;
	//SAQuantEst psSAQuantEstim[3];
	//SASummary TempSASummary; 
	//SAEstim *TempSurvivEstim;
//
	//nRet = stats_kaplan_meier_surviv_fn(vTime, vCensor, vCensored, dConLevel,
	//psSurvivEstim, nSizeSAEstm, psSAQuantEstim, 3, TempSASummary);
	//
	////Do not return value until stats_kaplan_meier_surviv_fn run successful
	//if (nRet == STATS_NO_ERROR) 
	//{
		//stSASummary = TempSASummary;
	//}
	//else 
		//return nRet;
	//
	//return nRet;
//}
//
/// Iris 10/21/04 temporary function body
int stats_one_sample_wilcoxon(const vector& vData, const double dMedian,const int nTailType, NPHSummary& Wilco1Summary, NPHStats& Wilco1Stats)
{
	return ocmath_one_sample_wilcoxon(vData, vData.GetSize(), dMedian, nTailType, &Wilco1Summary, &Wilco1Stats);
}

int stats_mann_whitney(const vector& vData1, const vector& vData2, const int nTailType, NPHSummary *MannSummary, 
	NPHRanks *MannRanks, NPHStats &MannStats, const bool bExact)
{
	return ocmath_mann_whitney(vData1, vData1.GetSize(),vData2, vData2.GetSize() ,nTailType, MannSummary, MannRanks, &MannStats, bExact);
 
}

int stats_wilcoxon_sign_rank(const vector& vData1, const vector& vData2, const int nTailType, NPHSummary *WilSiSummary, 
	NPHRanks *WilSiRanks, NPHStats &WilSiStats)
{
	return ocmath_wilcoxon_sign_rank(vData1, vData1.GetSize(), vData2, vData2.GetSize(), nTailType, WilSiSummary, WilSiRanks, &WilSiStats);
	
}

int stats_kruskal_wallis_anova(const vector& vData, const int iLevels, const vector& vIndex,
	NPHSummary *KruSummary, NPHRanks *KruRanks, NPHStats &KruStats)
{
	return ocmath_kruskal_wallis_anova(vData, vData.GetSize(), iLevels, vIndex, KruSummary, KruRanks, &KruStats);

}

int stats_friedman_anova(const vector& vDEP_VAR, const vector vBlocks, const vector& vTreaments,
	NPHSummary *FrieSummary, NPHRanks *FrieRanks, NPHStats &FrieStats)
{
	return ocmath_friedman_test(vDEP_VAR, vDEP_VAR.GetSize(), vTreaments, vTreaments.GetSize(), vBlocks, vBlocks.GetSize(), FrieSummary, FrieRanks, &FrieStats);
	
}

/*
int stats_spearman_corr_coeff(const vector& vData, const vector<int>& vIndex, matrix& mCorr, matrix& mSig)
{
	return ocmath_spe_corr_coeff(vData, vData.GetSize(), vIndex.GetSize(), vIndex, mCorr, mSig);
	
}

int stats_pearson_corr_coeff(const vector& vData, const vector<int>& vIndex, matrix& mCorr, matrix& mSig)
{	
	return ocmath_pea_corr_coeff(vData, vData.GetSize(), vIndex.GetSize(), vIndex, mCorr, mSig);	
}

int stats_kendall_corr_coeff(const vector& vData, const vector<int>& vIndex, matrix& mCorr, matrix& mSig)
{	
	return ocmath_ken_corr_coeff(vData, vData.GetSize(), vIndex.GetSize(), vIndex, mCorr, mSig);	
}
*/


/////////////////////////////////////////////////////////////////////////////////
//
//				ROC Curve									
//
/////////////////////////////////////////////////////////////////////////////////
int stats_roc_curve(const vector vData, const vector vState, const vector vStaValList, ROCArea& sROCArea, ROCCuv *psROC, uint nROC, double dConfLevel, int nThdFlag, const vector* pvThres)
{
	if (vData.GetSize() <= 0)
		return STATS_ERROR_TOO_FEW_DATA_PTS;
	
	if (vState.GetSize() <= 0)
		return STATS_ERROR_SETTING;
	
	if (vStaValList.GetSize() <= 0)
		return STATS_ERROR_SETTING;
		
	if ( !pvThres )
		return ocmath_roc_curve(vData, vData.GetSize(), vState, vState.GetSize(), vStaValList, vStaValList.GetSize(), &sROCArea, psROC, nROC, dConfLevel, nThdFlag);
	else
		return ocmath_roc_curve(vData, vData.GetSize(), vState, vState.GetSize(), vStaValList, vStaValList.GetSize(), &sROCArea, psROC, nROC, dConfLevel, nThdFlag, *pvThres, pvThres->GetSize());

}

/////////////////////////////////////////////////////////////////////////////////
//
//				Others									
//
/////////////////////////////////////////////////////////////////////////////////
void InitQuantiles(QuantileOptions& opQuantiles)
{
	Tree trTemp;
	TreeNode trQuantiles = trTemp.AddNode("Quantiles");
	trQuantiles = opQuantiles;
	foreach(TreeNode trChild in trQuantiles.Children)
		trChild.nVal = 1;			
	opQuantiles = trQuantiles;
	
	//opQuantiles.Interpolate = 0;
	opQuantiles.Interpolate = INTERPOLATE_WEIGHT_AVER_RIGHT;
	
}

void InitDescStats(DescStatOptions& opDescStats)
{
	Tree trTemp;
	TreeNode trDescStats = trTemp.AddNode("DescStats");
	trDescStats = opDescStats;
	foreach(TreeNode trChild in trDescStats.Children)
		trChild.nVal = 1;			
	opDescStats = trDescStats;
	opDescStats.VarDivisor = 0;
		
}

void desc_stats_group(const vector& vData, const vector& vGroup, DescStatResults *pDescStats, QuantileResults *pQuantiles, int nOverall)
{
	vector<int> vIndex;
	vector vDataTmp;
	int iLevels;
	int nRet = stats_anova_get_data_in_group(vData, vGroup, vDataTmp, iLevels, vIndex);
	if (nRet != STATS_NO_ERROR)
		return;

	vector vSubData;
	vIndex.InsertAt(0,0);
	QuantileOptions stQuanOpt;
	InitQuantiles(stQuanOpt);
	DescStatOptions stDescStatsOpt;
	InitDescStats(stDescStatsOpt);
	stDescStatsOpt.Mode = false;
	double dEnd = vIndex[0], dBegin;
	for (int ii = 0; ii < iLevels; ii++)
	{
		dBegin += vIndex[ii];
		dEnd += vIndex[ii+1];
		vData.GetSubVector(vSubData, dBegin, dEnd-1);
		stats_desc_stats(vSubData, stDescStatsOpt, pDescStats[ii]);
		if (nRet != STATS_NO_ERROR)
			return;
		stats_quantiles(vSubData, stQuanOpt, pQuantiles[ii]);
		if (nRet != STATS_NO_ERROR)
			return;
	}
	
	if (nOverall)
	{
		stats_desc_stats(vData, stDescStatsOpt, pDescStats[iLevels]);
		if (nRet != STATS_NO_ERROR)
			return;
		stats_quantiles(vData, stQuanOpt, pQuantiles[iLevels]);
		if (nRet != STATS_NO_ERROR)
			return;
		
	}
		
}

///Echo 5/10/06 MOVE_TO_ANALYSIS_C
///move ananlysis related functions to analysis.c

//static void _remove(const vector<uint>& vnIndex, vector<string>& vstrData)
//{
	//vector<string> vstrTmp;
	//int jj = 0;
	//for (int ii = 0; ii < vstrData.GetSize(); ii++)
	//{
		//if (ii != vnIndex[jj])
			//vstrTmp.Add(vstrData[ii]);
		//else
			//jj = (jj < vnIndex.GetSize()-1) ? jj+1 : jj;
			//
	//}
	//vstrData = vstrTmp;
//}
//
//bool trim_nanum_1group(vector &vData, vector<string> &vstrGroup, int& nMissing)
//{
	//if (vData.GetSize() != vstrGroup.GetSize())
		//return false;
//
	//vector<uint> vnIndex;
	//nMissing = vData.Find(vnIndex, NANUM);
	//if (nMissing > 0)
	//{
		//_remove(vnIndex, vstrGroup);
		//vData.Trim();
	//}
	//
	//return true;
	//
//}
//
//bool trim_nanum_2group(vector &vData, vector<string> &vsGroup1, vector<string> &vsGroup2, int& nMissing)
//{
	//if (vData.GetSize() != vsGroup1.GetSize() || vsGroup1.GetSize() != vsGroup2.GetSize())
		//return false;
	//
	//vector<uint> vnIndex;
	//nMissing = vData.Find(vnIndex, NANUM);
	//if (nMissing > 0)
	//{
		//_remove(vnIndex, vsGroup1);
		//_remove(vnIndex, vsGroup2);
		//vData.Trim();
	//}
//
	//return true;
//}
//
/////Eric  06/3/22    
//int trim_matrix(matrix &mat, vector &vCensor, vector<string> &vGroup)
//{
	//vector vTemp;
	//int nTemp, jj = 0, nMissing = 0;
	//
	//for(int ii=0;ii<mat.GetNumRows();ii++)
	//{
	    //mat.GetRow(vTemp,ii);
	    //if(0 < (nTemp = ocmath_count(NANUM, vTemp.GetSize(), vTemp)))
	    //{
	    	//nMissing += nTemp;
	    	//if(vCensor != NULL)
	    	//{
	    		//vCensor.RemoveAt(ii - jj);
	    	//}
	    	//if(vGroup != NULL)
	    	//{
	    		//vGroup.RemoveAt(ii - jj);
	    	//}
	    	//jj++;
	    //}
	//}
	    //
	//mat.RemoveEmptyRows(false);
	//
	//return nMissing;
//}
//
//int  trim_nanum_sa(int* nMissing, vector &vTime, vector &vCensor, vector<string> &vGroup, matrix &mCovariat)
//{
	//if(vTime.GetSize() != vCensor.GetSize())
	//{
		//return OE_INVALID_SIZE;
	//}
	//
	//if(vGroup != NULL)
	//{
		//if(vTime.GetSize() != vGroup.GetSize())
		//{
			//return OE_INVALID_SIZE;
		//}
	//}
	//
	//if(mCovariat != NULL)
	//{
		//if(vTime.GetSize() != mCovariat.GetNumRows())
		//{
			//return OE_INVALID_SIZE;
		//}
	//}
	//
	//matrix mTemp;
	//
	//if(mCovariat != NULL)
	//{
		//mTemp.SetSize(vTime.GetSize(), 1 + mCovariat.GetNumCols());
		//mTemp.SetColumn(vTime, 0);
		//mTemp.SetSubMatrix(mCovariat, 1);
	//}	
	//else if(vTime != NULL)
	//{
		//mTemp.SetSize(vTime.GetSize(), 1);
		//mTemp.SetColumn(vTime, 0);
	//
	//}
	//
		//
	//int i;
	//*nMissing = 0;	
	//*nMissing += trim_matrix(mTemp, vCensor, vGroup);
	//
	//
    //if(mTemp.GetNumRows() == 0)
    //{
    	//return OE_RANGE_ZERO;
    //}
    //else
    //{
	  	//i = 0;
	  	//if(vTime != NULL)
		//{
	 		//mTemp.GetColumn(vTime, i++);
		//}
	//
		//if(mCovariat != NULL)
		//{
	    	//mTemp.GetSubMatrix(mCovariat, i);
		//}
	//
		//return OE_NOERROR;
    //}
//}
///END MOVE_TO_ANALYSIS_C

///Echo 5/10/06 MOVE_TO_REPORT_UTILS_C
///move report related functions to report_utils.c
//
/////Eric 2006/4/20
//static void _get_desc_stats(DescStatResults *pDescStats, QuantileResults *pQuantiles, int ii, vector &vStats)
//{
	//vStats.SetSize(6);
	//vStats[0] = pDescStats[ii].N;
	//vStats[1] = pQuantiles[ii].Median;
	//vStats[2] = pQuantiles[ii].Min;
	//vStats[3] = pQuantiles[ii].Q25;
	//vStats[4] = pQuantiles[ii].Q75;
	//vStats[5] = pQuantiles[ii].Max;	
//}
//
//static void _get_summary(SASummary stSASumStats, vector& vSummary)
//{
	//vSummary.SetSize(4);
	//vSummary[0] = stSASumStats.Total;
	//vSummary[1] = stSASumStats.Events;
	//vSummary[2] = stSASumStats.NumCen;
	//vSummary[3] = stSASumStats.PerCen;
//}
//
/////Eric 2006/4/17
//void report_tree_create(ReportTree& outRT, LPCSTR lpsczLabel)
//{
	//outRT.ID = 12345;
	//string strFunctionLabel = _L(lpsczLabel);
	//time_t ltime;
	//time(&ltime);
	//
	//SYSTEMTIME st;
	//TM tmLocal;
	//double dDate;
	//convert_time_to_local(&ltime , &tmLocal);
    //tm_to_systemtime(&tmLocal, &st);
	//SystemTimeToJulianDate(&dDate, &st);
		//
	//strFunctionLabel += " (" + get_date_str(dDate) + ")";
	//outRT.SetAttribute(STR_LABEL_ATTRIB, strFunctionLabel);
	//
	////////////////////////////////////////////////////////////////////
	//// Notes
	////////////////////////////////////////////////////////////////////
	//ReportTable rtNotes = outRT.CreateTable("Notes", _L("Notes"), IDST_REPORT_NOTES);	
	//rtNotes.SetAttribute(STR_ATTRIB_BRANCH, GETNBRANCH_KEEP_SIZE_ON_COLLAPSE);
	//rtNotes.AddRow("Description", _L("Description"), _L(lpsczLabel), IDE_NOTES_DESCRIPTION);
	//char szUserName[MAX_PATH];
	//DWORD dwSize = MAX_PATH;
	//if( GetUserName(szUserName, &dwSize) )
		//rtNotes.AddRow("User", _L("User Name"), szUserName, IDE_NOTES_USERNAME);
	//
	//rtNotes.AddRow("Time", _L("Time"), get_date_str(dDate), IDE_NOTES_CREATION_TIME);
//}
//
//
/////Eric 2006/4/20
//void rt_add_nph_descriptive(ReportTree& outRT, const DescStatResults *pDescStats, const QuantileResults *pQuantiles, const vector<string>& vstrGroup, const int nGroups)
//{
	//ReportTable rtStats = outRT.CreateTable("DescStats", _L("Descriptive Statistics"), IDST_DESC_STATS_RESULTS);
	//rtStats.SetAttribute(STR_ATTRIB_BRANCH, GETNBRANCH_KEEP_SIZE_ON_COLLAPSE);
	//
	//vector<string> vStatsLabels(6);
	//vStatsLabels[0] = _L("N");
	//vStatsLabels[1] = _L("Median");
	//vStatsLabels[2] = _L("Min");
	//vStatsLabels[3] = _L("Q1");
	//vStatsLabels[4] = _L("Q3");
	//vStatsLabels[5] = _L("Max");
//
	//vector vStats;
//
	//
	//for (int ii = 0; ii < nGroups; ii++)
	//{
		//_get_desc_stats(pDescStats, pQuantiles, ii, vStats);
	    //rtStats.AddRow(vStats, _L(vstrGroup[ii]), vStatsLabels, make_one_set_ID(IDST_DESC_STATS_RESULTS, ii));
	//}
	//
	//
//}
//
//void rt_add_sa_summary(ReportTree& outRT, SASummary stSASumStats, LPCSTR lpsczLab, const int nIndex)
//{
	//ReportTable rtSummary = outRT.CreateTable("Summary", _L("Summary of Event and Censored Values"), 0x00202015);  //IDST_DESC_SUMMARY_RESULTS);	
	//vector<string> vSummaryLabels(4);
	//vSummaryLabels[0] = _L("Total");
	//vSummaryLabels[1] = _L("Events");
	//vSummaryLabels[2] = _L("Censored");	
	//vSummaryLabels[3] = _L("Percent Censored");
		//
	//vector vSummary;
	//
//
	//_get_summary(stSASumStats, vSummary);
	//rtSummary.AddRow(vSummary, _L(lpsczLab), vSummaryLabels, 0x00202015+nIndex); //IDST_DESC_SUMMARY_RESULTS+nIndex);		
			//
//}
///END MOVE_TO_REPORT_UTILS_C

///Echo 5/10/06 MOVE_TO_ANALYSIS_C
///move ananlysis related functions to analysis.c
/////Echo 4/19/06
//int get_2_data_in_group(const Range& rInGroup, vector& vData1, vector& vData2, vector<string>& vstrLab, int& nMissing)
//{
	////get data from range
	//vector<string> vstrGroup;
	//vector vData, vGroup;
	//
	//rInGroup.GetData(vData, 0);
	//rInGroup.GetData(vstrGroup, 1);
	//
	//if (!trim_nanum_1group(vData, vstrGroup, nMissing))
		//throw(STATS_ERROR_DATA_DIFF_SIZE);	
		//
	//vector vIndex;
	//get_group(vstrGroup, vGroup, vstrLab, vIndex);
	//int nGroup = vIndex.GetSize();
	//if (nGroup != 2)
		//return STATS_ERROR_TOO_FEW_GROUPS;
	//
	////put data to 2 variables
	//vData.GetSubVector(vData1, 0, vIndex[0]-1);
	//vData.GetSubVector(vData2, vIndex[0], vData.GetSize()-1);
	//
	//return STATS_NO_ERROR;
	//
//}
//
//int get_2_data_in_var(const Range& rInVar, vector& vData1, vector& vData2, vector<string>& vstrLab, int& nMissing, const bool bPair)
//{
	//rInVar.GetData(vData1, 0);
	//rInVar.GetData(vData2, 1);
//
	//nMissing = ocmath_count(NANUM, vData1.GetSize(), vData1) + ocmath_count(NANUM, vData2.GetSize(), vData2);
	//if (!bPair)
	//{
		//if (0 < nMissing)
		//{
			//vData1.Trim();
			//vData2.Trim();
		//}
	//}else
	//{
		//if (vData1.GetSize() != vData2.GetSize())
			//return STATS_ERROR_DATA_DIFF_SIZE;
		//
		//matrix mData(vData1.GetSize(),2);
		//
		//mData.SetColumn( vData1, 0);
		//mData.SetColumn( vData2, 1);
	//
		//mData.RemoveEmptyRows(false);
		//
		//mData.GetColumn( vData1, 0);
		//mData.GetColumn( vData2, 1);
	//}
	//
	//
	//vstrLab.SetSize(2);
	//for (int ii = 0; ii < 2; ii++)
		//vstrLab[ii] = range_get_col_name(rInVar, ii);
	//
	//return STATS_NO_ERROR;
//}
///END MOVE_TO_ANALYSIS_C


///Echo 5/10/06 MOVE_TO_REPORT_UTILS_C
///move report related functions to report_utils.c

/////Eric  5/10/2006
//static int 	_get_report_data_table_id(int nIndex)
//{
	//return  REPORT_TABLE + nIndex;
//}
//
//void	report_tree_make_result_picture(ReportTable& rt, ReportData& rd, int nIndex, int imageIndex)
//{
	//string strImageName;
	//strImageName.Format("Image%d", imageIndex);	
	//TreeNode trImage = rt.GetNode(strImageName);
	//if( trImage )
	//{
		//rd.PlotIntoTreeNode(trImage, _get_report_data_table_id(nIndex), NULL, "RtSurvivalPlot", nIndex, TRUE);
	//}
//}
//
//
//BOOL	calculate_results_curve_update_report_data(ReportTree& outRT, ReportData& outData, ReportTable& rtFitCurves,  const TreeNode& plot, double conf, vector& vx, vector& vy, vector vErrors, int nPoints, int nIndex, int nPlot, bool bOneLayer)
//{
	//if( nPoints <= 1 )
		//return FALSE;
	//
	//vector vFitX, vFitY, vUpperConf, vLowerConf;
	//vFitX.SetSize(vx.GetSize());
	//vFitY.SetSize(vx.GetSize());
	//vUpperConf.SetSize(vx.GetSize());
	//vLowerConf.SetSize(vx.GetSize());
    //int index = nIndex*nPlot;
    //string strTableName;
    //int nOneLayer = 0;
    //if(bOneLayer) nOneLayer = plot.plot.nVal;
    	//
	//if(plot.SF.Use)
	//{
		//vFitX = vx;
	    //vFitY = vy;
	    //
	    //for(int ii=0; ii<vFitY.GetSize(); ii++)
	    //{
	       	//vUpperConf[ii] = vFitY[ii] + (norminv((conf+1)/2)) * vErrors[ii];
	       	//vLowerConf[ii] = vFitY[ii] - (norminv((conf+1)/2)) * vErrors[ii];
	    //}
	    //
	    //if(nOneLayer == 0) 
	    //{
	    	//index++;
	    	//strTableName.Format("SFPlot%d", nIndex);	
	    	//ReportTable rtSFCurves = outData.CreateTable(strTableName, "Survival Function"+nIndex, _get_report_data_table_id(index));
	    	//rtSFCurves.AddColumn(vFitX, "X", IDST_FIT_X, NULL, OKDATAOBJ_DESIGNATION_X);
	    	//if(plot.SF.SFCI.nVal == 0)
	    	//{
	    		//rtSFCurves.AddColumn(vFitY, "Y", IDST_FIT_Y, NULL, OKDATAOBJ_DESIGNATION_Y);
	    	//}
	    	//else
	    	//{
	    		//rtSFCurves.AddColumn(vFitY, "Y0", IDST_FIT_Y, NULL, OKDATAOBJ_DESIGNATION_Y);
	    		//rtSFCurves.AddColumn(vUpperConf, "Y1", IDST_FIT_Y+1, NULL, OKDATAOBJ_DESIGNATION_Y);
	    		//rtSFCurves.AddColumn(vLowerConf, "Y2", IDST_FIT_Y+2, NULL, OKDATAOBJ_DESIGNATION_Y);
	    	//}	
	    //}
	    //else
	    //{
	    	//ReportTable rtSFCurves = outData.CreateTable("SFPlot", "Survival Function", _get_report_data_table_id(1));
	    	//rtSFCurves.AddColumn(vFitX, "X"+nIndex, IDST_FIT_X+nIndex, NULL, OKDATAOBJ_DESIGNATION_X);
	    	//rtSFCurves.AddColumn(vFitY, "Y"+nIndex, IDST_FIT_Y+nIndex, NULL, OKDATAOBJ_DESIGNATION_Y);
	    //}
	//}
	//
	//if(plot.OMSF.nVal)
	//{
		//vFitX = vx;
		//vFitY = 1-vy;
		//
		//if(nOneLayer == 0)
	    //{
			//index++;
	    	//strTableName.Format("OMSFPlot%d", nIndex);
			//ReportTable rtOMSFCurves = outData.CreateTable(strTableName, "One Minus Survival Function"+nIndex, _get_report_data_table_id(index));
	    	//rtOMSFCurves.AddColumn(vFitX, "X", IDST_FIT_X, NULL, OKDATAOBJ_DESIGNATION_X);
	    	//rtOMSFCurves.AddColumn(vFitY, "Y", IDST_FIT_Y, NULL, OKDATAOBJ_DESIGNATION_Y);
	    //}
	    //else
	    //{
	    	//ReportTable rtSFCurves = outData.CreateTable("OMSFPlot", "One Minus Survival Function", _get_report_data_table_id(2));
	    	//rtSFCurves.AddColumn(vFitX, "X"+nIndex, IDST_FIT_X+nIndex, NULL, OKDATAOBJ_DESIGNATION_X);
	    	//rtSFCurves.AddColumn(vFitY, "Y"+nIndex, IDST_FIT_Y+nIndex, NULL, OKDATAOBJ_DESIGNATION_Y);
	    //}
	//}
	//
	//if(plot.LSF.nVal)
	//{
		//vFitX = vx;
		//vFitY[0] = 0;
		//for(int i=1;i<vy.GetSize();i++)
		//{
			//if(vy[i]>0) vFitY[i] = log(vy[i]);
		//}
//
		//if(nOneLayer == 0)
	    //{
			//index++;
	    	//strTableName.Format("LSFPlot%d", nIndex);
	    	//ReportTable rtLSFCurves = outData.CreateTable(strTableName, "Log Survival Function"+nIndex, _get_report_data_table_id(index));
	    	//rtLSFCurves.AddColumn(vFitX, "X", IDST_FIT_X, NULL, OKDATAOBJ_DESIGNATION_X);
	    	//rtLSFCurves.AddColumn(vFitY, "Y", IDST_FIT_Y, NULL, OKDATAOBJ_DESIGNATION_Y);
	    //}
	    //else
	    //{
	    	//ReportTable rtSFCurves = outData.CreateTable("LSFPlot", "Log Survival Function", _get_report_data_table_id(3));
	    	//rtSFCurves.AddColumn(vFitX, "X"+nIndex, IDST_FIT_X+nIndex, NULL, OKDATAOBJ_DESIGNATION_X);
	    	//rtSFCurves.AddColumn(vFitY, "Y"+nIndex, IDST_FIT_Y+nIndex, NULL, OKDATAOBJ_DESIGNATION_Y);
	    //}
	//}
	//
	//return TRUE;
//}
//
//
//void create_plot_curve_table(ReportTree& outRT, const TreeNode& plot, const int& nGroup, bool bOneLayer)
//{
	//ReportTable rtFitCurves = outRT.CreateTable("FitCurves", _L("Fitted Curves"), IDST_RESULT_GRAPHS);
	//
	//string strImageName, strGroupName;
	//int nOneLayer = 0;
    //if(bOneLayer) nOneLayer = plot.plot.nVal;
	//
	//if(plot.SF.Use)
	//{
		//ReportTable rtSFCurves = rtFitCurves.CreateTable("SF", _L("Survival Function"), IDST_RESULT_GRAPHS+1);
	    //if(nOneLayer == 0)
	    //{
	    	//for(int ii=0; ii<nGroup; ii++)
	    	//{
	            //strImageName.Format("Image%d", ii);
	            //strGroupName.Format("Group%d", ii);
	            //rtSFCurves.AddRow(strImageName, _L(strGroupName), "", IDE_IMAGE+ii);
	    	//}    
	    //}
	    //else
	    //{
		   //rtSFCurves.AddRow("Image0", _L("All Image"), "", IDE_IMAGE);
	    //}   
	    //rtSFCurves.SetAttribute(STR_TYPE_ATTRIB, TNVAL_TYPE_PICTURE);
	//}
	//
	//if(plot.OMSF.nVal)
	//{
		//ReportTable rtOMSFCurves = rtFitCurves.CreateTable("OMSF", _L("One Minus Survival Function"), IDST_RESULT_GRAPHS+2);
	    //if(nOneLayer == 0)
	    //{
	    	//for(int ii=0; ii<nGroup; ii++)
	    	//{
	            //strImageName.Format("Image%d", ii);
	            //strGroupName.Format("Group%d", ii);
	            //rtOMSFCurves.AddRow(strImageName, _L(strGroupName), "", IDE_IMAGE+ii);
	    	//}    
	    //}
	    //else
	    //{
		   //rtOMSFCurves.AddRow("Image0", _L("All Image"), "", IDE_IMAGE);
	    //}   
	    //rtOMSFCurves.SetAttribute(STR_TYPE_ATTRIB, TNVAL_TYPE_PICTURE);
	//}
	//
	//if(plot.LSF.nVal)
	//{
		//ReportTable rtLSFCurves = rtFitCurves.CreateTable("LSF", _L("Log Survival Function"), IDST_RESULT_GRAPHS+3);
	    //if(nOneLayer == 0)
	    //{
	    	//for(int ii=0; ii<nGroup; ii++)
	    	//{
	            //strImageName.Format("Image%d", ii);
	            //strGroupName.Format("Group%d", ii);
	            //rtLSFCurves.AddRow(strImageName, _L(strGroupName), "", IDE_IMAGE+ii);
	    	//}    
	    //}
	    //else
	    //{
		   //rtLSFCurves.AddRow("Image0", _L("All Image"), "", IDE_IMAGE);
	    //}   
	    //rtLSFCurves.SetAttribute(STR_TYPE_ATTRIB, TNVAL_TYPE_PICTURE);
	//}
//}
//
//
//void update_plot_curve(ReportTree& outRT, ReportData& outData, const TreeNode& plot, const int& nGroup, const int& nPlot, bool bOneLayer)
//{
	//int index, groups;
	//int nOneLayer = 0;
    //if(bOneLayer) nOneLayer = plot.plot.nVal;
    //
	//if(nOneLayer == 1)
	//{
		//groups = 1;
	//}
	//else groups = nGroup;
	//
	//for(int ii = 0; ii<groups; ii++)
	//{
		//index = ii*nPlot;
		//if(plot.SF.Use)
		//{
			//index++;
			//report_tree_make_result_picture(outRT.FitCurves.SF, outData, index, ii);
	    	//outData.PlotIntoLayer(NULL, _get_report_data_table_id(index), NULL, index, FALSE);
		//}
	//
		//if(plot.OMSF.nVal)
		//{
			//index++;
			//report_tree_make_result_picture(outRT.FitCurves.OMSF, outData, index, ii);
	    	//outData.PlotIntoLayer(NULL, _get_report_data_table_id(index), NULL, index, FALSE);
		//}
	//
		//if(plot.LSF.nVal)
		//{
			//index++;
			//report_tree_make_result_picture(outRT.FitCurves.LSF, outData, index, ii);
	    	//outData.PlotIntoLayer(NULL, _get_report_data_table_id(index), NULL, index, FALSE);
		//}
	//}
//}
//
///END MOVE_TO_REPORT_UTILS_C

////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////       Temporary Test Routines          /////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////