/*------------------------------------------------------------------------------*
 * File Name: Internal.c														*
 * Creation: CPY 3/5/2001														*
 * Purpose: Support for using Origin C from inside Origin						*
 *			This is used primarily in places like NLSF, Set Matrix Values and	*
 *			function plotting and etc.											*
 * Copyright (c) OriginLab Corp.2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007	*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *   Bill Liang,1/16/02
 *   Bill Liang,1/17/02
 *	CPY v7.0365 QA70-2583 LABTALK_AND_OC_ABS_CONSISTENT		
 *  LAS 1/31/03 commented out the NLSFCntrl
 *              using NLSF = LabTalk.NLSF can be used instead
 * CPY 7/7/03 changed to include only origin.h to take advantage of PCH			
 * Danice 9/9/04 v8.0130 REPLACE_DATASET_WITH_VECTORBASE						*
 * Danice 9/11/04 v8.0132 MODIFY_REPLACE_DATASET_WITH_VECTORBASE				*
 * SDB 9/22/04 v7.5856 QA70-6185 REMOVE_FABS_FROM_CONDITIONAL					*
 * Leo 03/30/05 QA70-7449 REWRITE_FITPOLY										*
 * Leo 04/27/05 QA70-7449 SMOOTH_FFT_FILTER										*
 * Leo 06/10/05 USE_OCMATH_INTEGRATE 											*
 * Leo 06/11/05 REWRITE_XATYMIN_AND_XATYMAX										*
 * Cloud 03/24/07 QA80-9093 v8.0589 INITIAL_PARAMETERS                          *
 * Cloud 04/03/07 ADD_PEAK_HEIGHT_CALCULATION                                   *
 * Cloud 07/17/07 ADD_FFT_FOR_COMPATIBILITY										*
 * Arvin 08/10/07 QA70-10197 ADD_NONE_OPTION_FOR_ERROR_AS_WEIGHT_IN_LR_PR_AND_MR*
 * Arvin 08/30/07 SHOULD_TRIM_MISSING_VALUES_BEFORE_FIT							*
 * Echo 9/12/07	FITPOLY_ADD_BOUND_SETTING										*
 * Cloud 09/17/2007 FIT_WHOLE_CURVE_BY_DEFAULT									*
 * Cloud 12/06/07 FIX_BUG_OF_SECTION_LENGTH										*
 *	ML 12/19/2007 NEGATIVE_WIDTH_BUG											*
 * Cloud 12/25/07 UPDATE_METHOD_FOR_BASELINE									*
 * Cloud 01/14/08 REPRESENT_SCALE_BY_EXPONENT_FOR_HUGE_NUMBER					*
 * Cloud 02/21/08 PROVIDE_SEGMENT_CENTER										*
 * Cloud 02/21/08 REMOVE_0_FROM_SLOPES											*
 * Cloud 02/29/08 PX_NOT_RESTORE_BUG_BREAK_NLFIT_PARAM_INIT_FOR_EXP_DECAY		*
 * EJP 2008-04-08 v8.0840 QA80-11237 MAP_ROI_MATRIX_RECT_TO_IMAGE_RECT			*
 * Cloud 04/10/08 PASS_X_NOMALIZED_STATUS_BY_PARAMETER							*
 *	Folger 08/11/08 QA80-11995 COMPARE_LOCALIZED_GRAPH_OBJECT_STRING			*
 * Jack QA-12334 10/09/2008 ADD_ONE_PEAK_MATRIX_DATA_INITILIZE_CODE				*
 * CPY Iris 12/11/2008 v8.0985 QA80-12761 SUPPORT_LT_DATE_TIME_FUNCS			*
 *	Folger 12/12/08 QA80-12773 v8.0986 SUPPORT_SPECIFY_OC_FUNCTION_CATEGORY_IN_LABTALK_PRAGMA
 *	Jack 16/12/2008 QA80-12756 SUPPLY_OTHER_BASELINE_ESTIMATE_METHODS			*
 *	Fisher 05/27/09 QA80-13473 SOME_DATA_HAVE_LARGE_NOISE_LEVEL_WHICH_MAKE_SEGMENT_TESTING_UNABLE_TO_PASS_THREHOLD
 *	Sophy 6/26/2009 ADD_MORE_TRIANGLE_HYPERBOLIC_AND_INV_HYPERBOLIC_FUNCTIONS	*
 *	Sophy 8/26/2009 QA70-7298 COMPARE_BETWEEN_DOUBLE_TYPE_SHOULD_CHECK_MISSING_VALUES
 *	Sophy 9/15/2009 ADD_TOKEN_FUNCTION_ON_SET_COLUMNS_VALUE_DIALOG_MENU			*
 *  Iris 11/10/2009 QA80-13388 ZERO_CROSSING_DIRECTION							*
 *	Folger 11/28/09 BOLTZIV_FAILS_TO_INIT_PARAMS								*
 *	Folger 11/30/09 CENTRALIZE_CODE_ABOUT_FIND_ROOTS							*
 *	Sam 12/17/09 QA81-14596 MOVING_SLOPE_AND_RMS_FUNCTIONS						*
 *	CPY 2/19/10 OC_CALLING_SMOOTH_FUNC_WITH_TEXT_NUMERIC_COL_DATASET			*
 *	Folger 03/31/10 QA81-15258 XF_GETPTS_NEEDED									*
 *	Sim 04-08-2010 QA85-15284 SUPPORT_LT_FITTING_WITH_INTEGRAL					*
 *  Iris 6/09/2010 ORG-260 ADD_DERIVATIVE_AND_INTEGRAL_FUNC						*
 *------------------------------------------------------------------------------*/
 
#include <origin.h> // main Origin C header that is precompiled and already include most headers 

#pragma labtalk(0) //--- CPY 7/26/05 hide all functions in this file from LT access, as XF should be used instead
#include <labutils.h>  // CPY Iris 12/11/2008 v8.0985 QA80-12761 SUPPORT_LT_DATE_TIME_FUNCS

#include <..\OriginLab\fft.h> // Leo 04/27/05 SMOOTH_FFT_FILTER.

// Iris 12/29/2008 included for MDeterm function
#include	<NAG8\nag.h>
#include 	<NAG8\nagf03.h>
///---Sim 04-08-2010 QA85-15284 SUPPORT_LT_FITTING_WITH_INTEGRAL
#include <oc_nag8.h>
#include <okocUtils.h>
///---END QA85-15284 SUPPORT_LT_FITTING_WITH_INTEGRAL

///------ Folger 03/31/10 QA81-15258 XF_GETPTS_NEEDED
#include	<ocu.h>
///------ End XF_GETPTS_NEEDED

// *********************************************************************************
// OC functions for using in NLFit Parameter Init and NLFit Functions
// *********************************************************************************

///------ Folger 11/30/09 CENTRALIZE_CODE_ABOUT_FIND_ROOTS
// /// Fisher 4/1/2009 QA80-13388 ROOT_FINDING_FUNCTION
// // this function is used in the BoltsIV fitting function's parameter initialization code
// /// Iris 11/10/2009 QA80-13388 ZERO_CROSSING_DIRECTION
// //int		find_roots(const vector &vx, const vector &vy, vector &roots, double dBase)
// int		find_roots(const vector &vx, const vector &vy, vector &roots, double dBase, vector<int> &vnDirection)
// ///end ZERO_CROSSING_DIRECTION
// {
// 	int numRoots, nSize = vx.GetSize();
// 	if( 0 >= nSize && nSize != vy.GetSize() )
// 		return -1;
// 	
// 	vector vyy;
// 	vyy = vy - dBase;
// 	roots.SetSize(nSize);
// 	/// Iris 11/10/2009 QA80-13388 ZERO_CROSSING_DIRECTION, bad code, but no way to get int* from vector<int>.
// 	//ocmath_emd_get_zero_crossing(nSize, vx, vyy, &numRoots, roots);
// 	if( NULL == vnDirection )
// 		ocmath_emd_get_zero_crossing(nSize, vx, vyy, &numRoots, roots);
// 	else
// 	{
// 		vnDirection.SetSize(nSize);
// 		ocmath_emd_get_zero_crossing(nSize, vx, vyy, &numRoots, roots, vnDirection);
// 		///------ Folger 11/28/09 BOLTZIV_FAILS_TO_INIT_PARAMS
// 		vnDirection.SetSize(numRoots);
// 		///------ End BOLTZIV_FAILS_TO_INIT_PARAMS
// 	}
// 	///end ZERO_CROSSING_DIRECTION
// 	roots.SetSize(numRoots);
// 	///------ Folger 11/28/09 BOLTZIV_FAILS_TO_INIT_PARAMS
// 	//vnDirection.SetSize(numRoots);
// 	///------ End BOLTZIV_FAILS_TO_INIT_PARAMS
// 	return numRoots;
// }
// /// End ROOT_FINDING_FUNCTION
///------ End CENTRALIZE_CODE_ABOUT_FIND_ROOTS

//------- CPY v7.0365 QA70-2583 LABTALK_AND_OC_ABS_CONSISTENT	
// we need to provide this function such that the labtalk verison of abs
// will work the same as in Origin C. The C library function abs(n) would
// still be called if you are passing in an integer argument.
double abs(double x)
{
	return fabs(x);
}
//-------

// returns max of two doubles
double max(double da, double db)
{
	///Sophy 8/26/2009 QA70-7298 COMPARE_BETWEEN_DOUBLE_TYPE_SHOULD_CHECK_MISSING_VALUES
	if ( is_missing_value(da) || is_missing_value(db) )
		return get_missing_value();
	///end COMPARE_BETWEEN_DOUBLE_TYPE_SHOULD_CHECK_MISSING_VALUES
	return (da > db ? da : db);
}

// returns min of two doubles
double min(double da, double db)
{
	///Sophy 8/26/2009 QA70-7298 COMPARE_BETWEEN_DOUBLE_TYPE_SHOULD_CHECK_MISSING_VALUES
	if ( NANUM == da || NANUM == db )
		return NANUM;
	///end COMPARE_BETWEEN_DOUBLE_TYPE_SHOULD_CHECK_MISSING_VALUES
	return (da < db ? da : db);
}

// returns max of two floats
float max(float fa, float fb)
{
	return (fa > fb ? fa : fb);
}

// returns min of two floats
float min(float fa, float fb)
{
	return (fa < fb ? fa : fb);
}  

// returns max of two ints
int max(int ia, int ib)
{
	return (ia > ib ? ia : ib);
}

// returns min of two ints
int min(int ia, int ib)
{
	return (ia < ib ? ia : ib);
}  

// returns max of two uints
uint max(uint na, uint nb)
{
	return (na > nb ? na : nb);
}

// returns min of two uints
uint min(uint na, uint nb)
{
	return (na < nb ? na : nb);
}

// returns max of two shorts
short max(short na, short nb)
{
	return (na > nb ? na : nb);
}

// returns min of two shorts
short min(short na, short nb)
{
	return (na < nb ? na : nb);
}

// returns max of two ushorts
ushort max(ushort wa, ushort wb)
{
	return (wa > wb ? wa : wb);
}

// returns min of two ushorts
ushort min(ushort wa, ushort wb)
{
	return (wa < wb ? wa : wb);
}


// ***************************
// functions that accept a dataset and return values from its properties
// ***************************

///Danice REPLACE_DATASET_WITH_VECTORBASE
/*
// returns min value of dataset
double min(Dataset &aa)
{
	BasicStats	stat;
	stat.min = -1;
	Data_sum(&aa, &stat);
	return stat.min;
}

// returns max value of dataset
double max(Dataset &aa)
{
	BasicStats	stat;
	stat.max = -1;
	Data_sum(&aa, &stat);
	return stat.max;
}
*/
///Danice MODIFY_REPLACE_DATASET_WITH_VECTORBASE
// returns 25th percentile value of dataset
double y25(Dataset &aa)
{
	//Sandy 2006-8-2 REPLACE_WITH_OCMATH_FUNCTION
	//double perc, result;
	//BOOL interp = false;
	//perc=25;
	//Data_percentiles(&aa, &perc, &result, 1,interp);
	//return result;
	vector vData;
	vData = aa;
	vector vPerc(1);
	vPerc[0] = 25;
	vector vPercResult(1);
	vPercResult.SetSize(vPerc.GetSize());
	int nErr=ocmath_percentiles(vData, vData.GetSize(), vPerc, vPerc.GetSize(), vPercResult);
	if(nErr !=  STATS_NO_ERROR)
		return nErr;
	
	return vPercResult[0];
	
}

// returns 50th percentile value of dataset
double y50(Dataset &aa)
{
	//Sandy 2006-8-2 REPLACE_WITH_OCMATH_FUNCTION
	//double perc, result;
	//BOOL interp = false;
	//perc=50;
	//Data_percentiles(&aa, &perc, &result, 1,interp);
	//return result;
	
	vector vData;
	vData = aa;
	vector vPerc(1);
	vPerc[0] = 50;
	vector vPercResult(1);
	vPercResult.SetSize(vPerc.GetSize());
	int nErr=ocmath_percentiles(vData, vData.GetSize(), vPerc, vPerc.GetSize(), vPercResult);
	if(nErr !=  STATS_NO_ERROR)
		return nErr;
	
	return vPercResult[0];
}
		
// returns 75th percentile value of dataset
double y75(Dataset &aa)
{
	//Sandy 2006-8-2 REPLACE_WITH_OCMATH_FUNCTION
	//double perc, result;
	//BOOL interp = false;
	//perc=75;
	//Data_percentiles(&aa, &perc, &result, 1,interp);
	//return result;
	vector vData;
	vData = aa;
	vector vPerc(1);
	vPerc[0] = 75;
	vector vPercResult(1);
	vPercResult.SetSize(vPerc.GetSize());
	int nErr=ocmath_percentiles(vData, vData.GetSize(), vPerc, vPerc.GetSize(), vPercResult);
	if(nErr !=  STATS_NO_ERROR)
		return nErr;
	
	return vPercResult[0];
}
///end

// returns min value of vectorbase
double 	min(vectorbase &vec)
{
	double dMin, dMax;
	vec.GetMinMax(dMin, dMax);
	return dMin;
}

// returns max value of vectorbase
double	max(vectorbase &vec)
{
	double dMin, dMax;
	vec.GetMinMax(dMin, dMax);
	return dMax;
}

// returns 25th percentile value of vector
double	y25(vector &vec)
{
	QuantileOptions opt;
	opt.Q1=true;
	QuantileResults res;
	ocmath_quantiles(vec, vec.GetSize(), &opt, &res);
	return res.Q1;
}

// returns 50th percentile value of vector
double y50(vector &vec)
{
	QuantileOptions opt;
	opt.Median=true;
	QuantileResults res;
	ocmath_quantiles(vec, vec.GetSize(), &opt, &res);
	return res.Median;
}
		
// returns 75th percentile value of vector
double y75(vector &vec)
{
	QuantileOptions opt;
	opt.Q3=true;
	QuantileResults res;
	ocmath_quantiles(vec, vec.GetSize(), &opt, &res);
	return res.Q3;
}
///end

// ***************************
// functions that accept a curve and return values from its properties
// ***************************

///Leo 06/11/05 REWRITE_XATYMIN_AND_XATYMAX
// this function is used only for the two functions below, so must be set to static
// so that it is invisible outside this file
//static int minmax(Dataset &aa, BOOL bGetMin=TRUE)
//{
//	BasicStats	stat;
//	stat.iMin =0;
//	stat.iMax =0;
//	Data_sum(&aa, &stat);
//	return bGetMin? stat.iMin:stat.iMax;
//}
//
//// returns x value corresponding to minimum in y value
//double xatymin(Curve &cc)
//{
//	return Curve_x(&cc, minmax(cc));
//}
//
//// returns x value corresponding to minimum in y value
//double xatymax(Curve &cc)
//{
//	return Curve_x(&cc, minmax(cc,FALSE));
//}
// returns x value corresponding to minimum in y value
double xatymin(Curve &cc)
{
	Curve cc2(cc);
	Dataset dsX;
	cc2.AttachX(dsX);
	int iMin;
	ocmath_d_sum(cc2, cc2.GetSize(), NULL, NULL, &iMin);
	return dsX[iMin];
}
// returns x value corresponding to maximum in y value
double xatymax(Curve &cc)
{
	Curve cc2(cc);
	Dataset dsX;
	cc2.AttachX(dsX);
	int iMax;
	ocmath_d_sum(cc2, cc2.GetSize(), NULL, NULL, NULL, &iMax);
	return dsX[iMax];
}
///End REWRITE_XATYMIN_AND_XATYMAX

double yatxmax(Curve &cc)
{
	if (cc.HasX())
	{
		Dataset xdata();
		cc.AttachX(xdata);
		double xmax = max(xdata);
		return Curve_yfromX(&cc, xmax);
	}
	else
	{
		// GetUpperIndex returns the index immediately after the last point
		int n = cc.GetUpperBound();
		return cc[n-1];
	}	
}

double yatxmin(Curve &cc)
{
	if (cc.HasX())
	{
		Dataset xdata();
		cc.AttachX(xdata);
		double xmin = min(xdata);
		return Curve_yfromX(&cc, xmin);
	}
	else
	{
		return cc[0];
	}	
}

//returns x value correpsonding to mid y value
double xaty50(Curve &cc)
{
	double ymiddle = (max(cc) + min(cc))/2.0;
	return Curve_xfromY(&cc, ymiddle);	
}

///Leo 06/10/05 USE_OCMATH_INTEGRATE
//returns x width of dataset corresponding to half the max value in y
//double fwhm(Curve &cc, double yoffset = 0)
//double fwhm(Curve &cc, double yoffset)
//{
//	Curve ccLocal(cc);
//	ccLocal -= yoffset;
//	IntegrationResult result;
//	BOOL ipass = Curve_integrate(&ccLocal, &result);
//	return result.dxPeak;
//}
//
////returns area under the curve
////double area(Curve &cc, double yoffset = 0)
//double area(Curve &cc, double yoffset)
//{
//    Curve ccLocal(cc);
//	ccLocal -= yoffset;
//	IntegrationResult result;
//	BOOL ipass = Curve_integrate(&ccLocal, &result);
//	return result.Area;
//}
//Full width at half maximum of the curve
double fwhm(Curve &cc, double yoffset)
{
	Curve cc2(cc);
	cc2 -= yoffset;
	Dataset dsX;
	cc2.AttachX(dsX);
	///Arvin 03/19/07 IMPROVE_CALCULATION_CENTROID_METHOD
	//ocmath_IntegResult result;
	IntegrationResult result;
	///end IMPROVE_CALCULATION_CENTROID_METHOD
	ocmath_integrate(dsX, cc2, 0, dsX.GetSize()-1, &result);
	
	return result.dxPeak;
}
//Area under the curve
double area(Curve &cc, double yoffset)
{
	Curve cc2(cc);
	cc2 -= yoffset;	
	Dataset dsX;
	cc2.AttachX(dsX);
	///Arvin 03/19/07 IMPROVE_CALCULATION_CENTROID_METHOD
	//ocmath_IntegResult result;
	IntegrationResult result;
	///end IMPROVE_CALCULATION_CENTROID_METHOD
	ocmath_integrate(dsX, cc2, 0, dsX.GetSize()-1, &result);
	return result.Area;
}
///End USE_OCMATH_INTEGRATE

// sorts the curve
BOOL sort(Curve &cc)
{
	cc.TrimLeft(TRUE);
	cc.Sort();// missing values are sorted to the end
	cc.TrimRight();
	return true;
}	


//--- CPY QA70-6893 9/10/04 v8.0132 OCMATH_SMOOTH_IF_DOUBLE_CURVE
// smooth the vector
enum {SMOOTHING_ADJAVE, SMOOTHING_FFTFILTER, SMOOTHING_SG, SMOOTHING_MEDIAN_FILTER};
static bool _smooth(vector &vY, int method, int leftpts , int rightpts, int polydeg)
{
	vector vysrc;
	vysrc = vY;
	int iRet = -1;
	switch(method)
	{
	case SMOOTHING_ADJAVE:
		//return ocmath_adjave_smooth(vysrc, vY, vY.GetSize(), leftpts);	//Leo 08/29/05 QA70-7847 TO_BE_CONSISTENT_WITH_75
		///Mouqx/Leo 2005-9-26 QA70-8110 TO_CENTRALIZE_ERROR_CODES_IN_OCMATH
		//return ocmath_adjave_smooth(vysrc, vY, vY.GetSize(), leftpts/2);
		iRet = ocmath_adjave_smooth(vysrc, vY, vY.GetSize(), leftpts/2);
		return (iRet == OE_NOERROR);
	case SMOOTHING_FFTFILTER:
		///Leo 2006-02-06 BRING_BACK_PARABOLIC_LOW_PASS_FILTER
		///Leo 04/27/05 QA70-7449 SMOOTH_FFT_FILTER
		//DescStatResults dsr;
		//ocmath_desc_stats(vysrc, vysrc.GetSize(), &dsr);
		//double sigma = dsr.SD; // sigma is the noise std
		//return fft_wiener_filter(vY.GetSize(), vysrc, vY, sigma);
		if(leftpts < 0)
		{
			return false;
		}
		if (0 != ocmath_count(NANUM, vY.GetSize(), vY))
		{
			return false;
		}
		vector vWindow(vY.GetSize());
		if(0 != fft_filter_window(vY.GetSize(), LOW_PASS_PARABOLIC_FILTERING, vWindow, 0, 1./leftpts, 1./vY.GetSize(), true))
		{
			return false;
		}
		if(0 != fft_fft_filter(vY.GetSize(), vY, vWindow))
		{
			return false;
		}
		return true;
	case SMOOTHING_SG:
		//---CPY 2/19/10 OC_CALLING_SMOOTH_FUNC_WITH_TEXT_NUMERIC_COL_DATASET
		//return ocmath_savitsky_golay(vysrc, vY, vY.GetSize(), leftpts, rightpts, polydeg);
		iRet = ocmath_savitsky_golay(vysrc, vY, vY.GetSize(), leftpts, rightpts, polydeg); 
		return (iRet == OE_NOERROR);
		//---
		/// Leo 07/26/05 QA70-7088 ADD_OCMATH_MEDIAN_FILTER
	case SMOOTHING_MEDIAN_FILTER:
		///Cheney 2006-8-15 ADD_PADDING
		//int iRet = ocmath_median_filter(vY.GetSize(), vysrc, vY, leftpts, rightpts, polydeg);
		iRet = ocmath_median_filter(vY.GetSize(), vysrc, vY, leftpts, rightpts, EDGEPAD_ZERO, polydeg);
		return (iRet == OE_NOERROR);
		///end ADD_PADDING
	}
	return false;
}
///Sophy 5/7/2009 KEEP_BACKWARD_COMPATIBLE
bool smooth(Curve &crv, int imethod/*=0*/, int ileftpts/*= 3*/, int irightpts/*= 3*/, int ipolydeg/*= 2*/)
{
	vector vX, vTmpY;
	crv.CopyData(vX, vTmpY);
	int nSize = crv.GetSize();
	if ( _smooth(vTmpY, imethod, ileftpts, irightpts, ipolydeg) )
	{
		crv = vTmpY;
	 	return true;
	}
	return false;
}
///end KEEP_BACKWARD_COMPATIBLE
//------ CPY 2/19/10 OC_CALLING_SMOOTH_FUNC_WITH_TEXT_NUMERIC_COL_DATASET
// this function will need to support passing in Dataset, see http://beta.originlab.com/forum/topic.asp?TOPIC_ID=3100
bool smooth(vector &vY, int method, int leftpts , int rightpts, int polydeg)
{
	vector vTmpY;
	vTmpY = vY;
	if(_smooth(vTmpY, method, leftpts, rightpts, polydeg))
	{
		vY = vTmpY;
		return true;
	}
	return false;
}
//------

////--- end

////--- CPY QA70-6893 9/10/04 v8.0132 OCMATH_SMOOTH_IF_DOUBLE_CURVE
//enum {SMOOTH_ADJAVE, SMOOTH_FFTFILTER, SMOOTH_SG};
//bool vector_smooth(vector& vy, int method = SMOOTH_ADJAVE, int leftpts = 3, int rightpts = 3, int polydeg = 2)
//{
//	vector vysrc;
//	vysrc = vy;
//	switch(method)
//	{
//	case SMOOTH_ADJAVE:
//		return ocmath_adjave_smooth(vysrc, vy, vy.GetSize(), leftpts);
//	case SMOOTH_SG:
//		return ocmath_savitsky_golay(vysrc, vy, vy.GetSize(), leftpts, rightpts, polydeg);
//	case SMOOTH_FFTFILTER:
//		out_str("SMOOTH_FFTFILTER not implemented yet");
//		break;
//	}
//	return false;
//		
//}
////--- end
//
//// smooth the curve
////BOOL smooth(Curve &cc, int method=0, int leftpts = 3, int rightpts = 3, int polydeg = 2)
//BOOL smooth(Curve &cc, int method, int leftpts , int rightpts, int polydeg)
//{
//	string ydataname;
//	cc.GetName(ydataname);	
//	//--- CPY QA70-6893 9/10/04 v8.0132 OCMATH_SMOOTH_IF_DOUBLE_CURVE
//	DWORD dwInfo = Project.GetDatasetInfo(ydataname);
//	if(HIWORD(dwInfo) == FSI_DOUBLE && method != SMOOTH_FFTFILTER) // later will need to add SMOOTH_FFTFILTER
//	{
//		return vector_smooth(cc, method, leftpts, rightpts, polydeg);
//	}
//	//---
//	
//	if(method == 0)								// perform adjacent averaging
//	{
//		_LT_Obj
//		{
//			curve.reset();
//			curve.data$ = ydataname;
//			curve.result$ = ydataname;
//			curve.smoothPts = leftpts;
//			curve.adjave();
//		}
//		return true;
//	}
//	else if (method == 1)						// perform fft filtering
//	{
//		_LT_Obj
//		{
//			curve.reset();
//			curve.data$ = ydataname;
//			curve.result$ = ydataname;
//			curve.smoothPts = leftpts;
//			curve.FFTSmooth();
//		}
//		return true;
//	}
//	else										// perform Savitzky-Golay filtering
//	{
//		_LT_Obj
//		{
//			curve.reset();
//			curve.data$ = ydataname;
//			curve.result$ = ydataname;
//			curve.smoothLeftPts = leftpts;
//			curve.smoothRightPts = rightpts;
//			curve.smoothPts = leftpts + rightpts + 1;
//			curve.polyDeg = polydeg;
//			curve.SGSmooth();
//		}
//		return true;
//	}
//}
//


//-------- Bill Liang,1/16/02
// return the asymtoptic y value as follow: 
// double-branch: (yatxmax+yatxmin)/2
// single-branch: yatxmax or yatxmin with the smaller slope
double yatasymt(Curve &cc)
{
	Dataset xd;
	cc.AttachX(xd);
	if(fabs(xatymax(cc)-xatymin(cc))<(max(xd)-min(xd))*2.0/3)
	{ // double-branch
		return (yatxmin(cc)+yatxmax(cc))/2.0;
	}		
	else 
	{ // single-branch
		double dy1=Curve_yfromX(&cc,0.9*min(xd)+0.1*max(xd))-yatxmin(cc);
		double dy2=yatxmax(cc)-Curve_yfromX(&cc,0.1*min(xd)+0.9*max(xd));
		if(fabs(dy1)<fabs(dy2)) 
			return yatxmin(cc);
		else
			return yatxmax(cc);
	}
      
}

// return the asymtoptic x value as follow: 
// double-branch: (xatymax+xatymin)/2
// single-branch: xatymax or xatymin with the smaller 1/slope
//************** Bill Liang,1/17/02
double xatasymt(Curve &cc)
{
	Dataset xd;
	cc.AttachX(xd);
	if(fabs(yatxmax(cc)-yatxmin(cc))<(max(cc)-min(cc))*2.0/3)
	{ // double-branch
		if(fabs(xatymax(cc)-xatymin(cc))>(max(xd)-min(xd))*1.0/3) 
		{ // y-symemtric double-branch
			if(fabs(min(cc)-yatasymt(cc))<fabs(max(cc)-yatasymt(cc)))
				return xatymax(cc);
			else
				return xatymin(cc);
		}
		else	
			return (xatymin(cc)+xatymax(cc))/2.0;
	}		
	else 
	{ // single-branch
		double dy1=Curve_yfromX(&cc,0.9*min(xd)+0.1*max(xd))-yatxmin(cc);
		double dy2=yatxmax(cc)-Curve_yfromX(&cc,0.1*min(xd)+0.9*max(xd));
		//SDB v7.5856 QA70-6185 REMOVE_FABS_FROM_CONDITIONAL
		//if(fabs(dy1)<fabs(dy2)) 
		if(dy1<dy2) 
			return xatymax(cc);
		else
			return xatymin(cc);
	}      
}

///Leo 30/03/05 QA70-7449 REWRITE_FITPOLY CPY 5/17/05
static void _init_fit_options(LROptions& sLROptions)
{
	// Fill in LROptions with default values.
	sLROptions.FixSlope = false;		// Fixed value of slope
	sLROptions.FixSlopeAt = 0;			// Fixed value of slope
	///Arvin 08/10/07 QA70-10197 ADD_NONE_OPTION_FOR_ERROR_AS_WEIGHT_IN_LR_PR_AND_MR
	//sLROptions.ErrBarWeight = false;	// Errors for weight   
	sLROptions.ErrBarWeight = ERRBARWEIGHT_DIRECT_WEIGHTING;  // Errors for weight   
	///end ADD_NONE_OPTION_FOR_ERROR_AS_WEIGHT_IN_LR_PR_AND_MR
	sLROptions.UseReducedChiSq = false;	// Reduced Chi-square                    
	sLROptions.Confidence = 0.95;		// Confidence value for parameter U/L        
	sLROptions.ApparentFit = false;		// Apparent scale                        
	sLROptions.FixIntercept = false;	// Force intercept to be fixed               
	sLROptions.FixInterceptAt = 0;		// Fixed value of intercept                  
}

bool fitlinear(const vector& vx, const vector& vy, vector& vCoeff, const vector& vWt, RegStats* psRegStats, bool bThroughZero)
{
	if(vx.GetSize() != vy.GetSize()) 
			return false;
	
	LROptions sLROptions;
	_init_fit_options(sLROptions);
	sLROptions.FixIntercept = bThroughZero;
	
	vector vWeight;
	if(vWt)
		vWeight = vWt;
	
	if (vWeight.GetSize() < vy.GetSize())
	{
		int nOld = vWeight.GetSize();
		vWeight.SetSize(vy.GetSize());
		for (int ii = nOld; ii < vy.GetSize(); ii++)
			vWeight[ii] = 1;
	}

	FitParameter sFitParameter[2];
			
	int nErr = ocmath_linear_fit(vx, vy, vx.GetSize(), sFitParameter, vWeight, vWeight.GetSize(), &sLROptions, psRegStats);

	if(STATS_NO_ERROR == nErr)
	{
		vCoeff.SetSize(4);
		vCoeff[0] = sFitParameter[0].Value;
		vCoeff[1] = sFitParameter[1].Value;
		vCoeff[2] = sFitParameter[0].Error;
		vCoeff[3] = sFitParameter[1].Error;
		return true;
	}
	out_int("fitlinear err = ", nErr);
	return false;
}
//----------- CPY 1/26/2008 FITPR_XF
/*
// Polynomial fit base
static bool _fitpoly(int nSize, double* pX, double* pY, int nOrder, vector& vCoeff, vector& vError, RegStats* psRegStats)
{
	if (vCoeff && vCoeff.GetSize() <= nOrder)
		vCoeff.SetSize(nOrder + 1);
		
	if (vError && vError.GetSize() <= nOrder)
		vError.SetSize(nOrder + 1);

	LROptions sLROptions;	//[input]
	FitParameter* psFitParameter;	//[output]
	psFitParameter = (FitParameter*) malloc( (nOrder+1)*sizeof(FitParameter));
	
	_init_fit_options(sLROptions);

//	matrix mX(nSize, nOrder);
//	vector vTemp(nSize);
//
//	vTemp = 1;
//	for (int i = 0; i < nOrder; i++)
//	{
//		for (int j = 0; j < nSize; j++)
//			vTemp[j] *= pX[j];
//		mX.SetColumn(vTemp, i);
//	}
//	int nErr = ocmath_multiple_linear_regression(mX, nSize, nOrder, pY, NULL, 0, &sLROptions, psFitParameter, nOrder+1, 
//				psRegStats, NULL, NULL, NULL, 0, 0);

	int nErr = ocmath_polynomial_fit(nSize, pX, pY, NULL, nOrder, &sLROptions, psFitParameter, nOrder+1, 
					psRegStats, NULL, NULL, NULL);

	if (nErr != STATS_NO_ERROR)
		return false;

	for (int j = 0; j <= nOrder; j++)
	{
		if (vCoeff) vCoeff[j] = psFitParameter[j].Value;
		if (vError) vError[j] = psFitParameter[j].Error;
	}
	
	free(psFitParameter);

	return true;
}*/
int fitpoly(int nOrder, double* pX, double* pY, int nSize, vector& vCoeff, vector& vError, RegStats* psRegStats, bool* pbXNorm)
{
	if (vCoeff && vCoeff.GetSize() <= nOrder)
		vCoeff.SetSize(nOrder + 1);
		
	if (vError && vError.GetSize() <= nOrder)
		vError.SetSize(nOrder + 1);

	LROptions sLROptions;	//[input]
	FitParameter* psFitParameter;	//[output]
	psFitParameter = (FitParameter*) malloc( (nOrder+1)*sizeof(FitParameter));
	
	_init_fit_options(sLROptions);
	int nErr = ocmath_polynomial_fit(nSize, pX, pY, NULL, nOrder, &sLROptions, psFitParameter, nOrder+1, 
					psRegStats, NULL, NULL, NULL);

	/// Cloud 04/10/08 PASS_X_NOMALIZED_STATUS_BY_PARAMETER
	//if (nErr == STATS_NO_ERROR)
	//{
	if (nErr >= STATS_NO_ERROR)
	{
		if (pbXNorm!=NULL)
		{
			*pbXNorm = nErr>STATS_NO_ERROR ? true : false;
		}
		nErr = STATS_NO_ERROR;
	/// End Cloud 04/10/08 PASS_X_NOMALIZED_STATUS_BY_PARAMETER
		for (int j = 0; j <= nOrder; j++)
		{
			if (vCoeff) vCoeff[j] = psFitParameter[j].Value;
			if (vError) vError[j] = psFitParameter[j].Error;
		}
	}
	free(psFitParameter);
	return nErr;
}
// Polynomial fit base
static bool _fitpoly(int nSize, double* pX, double* pY, int nOrder, vector& vCoeff, vector& vError, RegStats* psRegStats)
{
	return fitpoly(nOrder, pX, pY, nSize, vCoeff, vError, psRegStats) == STATS_NO_ERROR? true:false;
}
//----- end FITPR_XF


// Fit polynomial to a curve bounded from iBeg to iEnd and return coefficients and erros
bool fitpoly(Curve& cv, int nOrder, vector& vCoeff, vector& vError, RegStats* psRegStats, int iBeg, int iEnd, double* pxmin, double *pxmax)
{
	//Curve cv2(cv); CPY 5/13/05, need to use more reliable constructor that takes iBeg and iEnd
	int nMissing, nSrcOffset;
	Curve cv2(cv, nMissing, nSrcOffset, CURVECOPY_SKIP_MISSING_INSIDE, iBeg, iEnd);
	
	//Leo 8/26/05 Curve is invalid if iBeg > iEnd.
	if (!cv2.IsValid())
		return false;
	
	Dataset dsX;
	cv2.AttachX(dsX);
	
	// Here pass cv2 to double* since cv2 is a copy.
	if(_fitpoly(cv2.GetSize(), dsX, cv2, nOrder, vCoeff, vError, psRegStats))
	{
		if(pxmin || pxmax)
		{
			double min, max;
			dsX.GetMinMax(min, max);
			if(pxmin) *pxmin = min;
			if(pxmax) *pxmax = max;
		}
		return true;
	}
	return false;
}

// Fit polynomial to X Y vectors and return coefficients and erros
bool fitpoly(vector& vX, vector& vY, int nOrder, vector& vCoeff, vector& vError, RegStats* psRegStats)
{
	if (vX.GetSize() != vY.GetSize())
		return false;
	
	///Arvin 08/30/07 SHOULD_TRIM_MISSING_VALUES_BEFORE_FIT
	//return 	_fitpoly(vX.GetSize(), vX, vY, nOrder, vCoeff, vError, psRegStats);
	Curve cuvTemp(vX, vY);
	int nMissing, nSrcOffset;
	Curve cuv(cuvTemp, nMissing, nSrcOffset, CURVECOPY_SKIP_MISSING_INSIDE);
	vector vxTmp, vyTmp;
	cuv.CopyData(vxTmp, vyTmp);
	return 	_fitpoly(vxTmp.GetSize(), vxTmp, vyTmp, nOrder, vCoeff, vError, psRegStats);
	///end SHOULD_TRIM_MISSING_VALUES_BEFORE_FIT
}
///END REWRITE_FITPOLY

///LEO (30/03/05) this function is obstacle, see REWRITE_FITPOLY
BOOL fitpoly_range(Curve &cc,int ibeg,int iend,int ipolyorder, double *coeff)
{
	vector vCoeff(ipolyorder+1);
	bool bRet = fitpoly(cc, ipolyorder, vCoeff, NULL, NULL, ibeg, iend);
	if (bRet)
		for (int i = 0; i <= ipolyorder; i++)
			coeff[i] = vCoeff[i];
		
	return bRet;
}

///LEO (30/03/05) this function is obstacle, see REWRITE_FITPOLY
BOOL fitpoly(Curve &crv, int ipolyorder, double *coeff, int inumer, int idenom)
{
	///Echo 9/12/07	FITPOLY_ADD_BOUND_SETTING
	int ilower, iupper, ilen, iseg, ibeg, iend;
	if  ( (idenom != 0) )
	{
		if (inumer > idenom) return false;
		ilower = crv.GetLowerBound();
		iupper = crv.GetUpperBound();
		ilen = iupper - ilower;
		if (ilen <= idenom) return false;
		iseg = ilen / idenom;
		ibeg = iseg * (inumer - 1);
		iend = ibeg + iseg;
	}
	/// Cloud 09/17/2007 FIT_WHOLE_CURVE_BY_DEFAULT
	else
	{
		ibeg = 0;
		iend = -1;
	}
	/// End FIT_WHOLE_CURVE_BY_DEFAULT
	//return fitpoly_range(crv, 0, -1, ipolyorder, coeff);
	return fitpoly_range(crv, ibeg, iend, ipolyorder, coeff);
	///END FITPOLY_ADD_BOUND_SETTING
}

////************** Bill Liang,1/17/02
// fit polynomial to curve on data points indexed from ibeg to iend
//BOOL fitpoly_range(Curve &cc,int ibeg,int iend,int ipolyorder, double *coeff)
//{
//	int ilower=cc.GetLowerBound();
//	int iupper=cc.GetUpperBound();
//	if(ibeg<ilower || iend>iupper || ibeg>=iend) 
//		return false;		
//	cc.SetLowerIndex(ibeg);
//	cc.SetUpperIndex(iend);
//	fitpoly(cc,ipolyorder,coeff);	
//	// restore original indics
//	cc.SetLowerIndex(ilower);
//	cc.SetUpperIndex(iupper);
//	return true;
//}
//
////-----------------	end Bill Liang,1/16/02

// fit polynomial to curve and return coefficients
//BOOL fitpoly(Curve &cc, int ipolyorder, double *coeff, int inumer, int idenom)
//{
//	string strYdata;
//	cc.GetName(strYdata);
//
//	int ilower, iupper, ilen, iseg, ibeg, iend;
//
//	if  ( (idenom != 0) )
//	{
//		if (inumer > idenom) return false;
//		ilower = cc.GetLowerBound();
//		iupper = cc.GetUpperBound();
//		ilen = iupper - ilower;
//		if (ilen <= idenom) return false;
//		iseg = ilen / idenom;
//		ibeg = iseg * (inumer - 1);
//		iend = ibeg + iseg;
//		cc.SetLowerIndex(ibeg);
//		cc.SetUpperIndex(iend);
//		}
//
//	_LT_Obj
//	{
//		stat.reset();
//		stat.apparent=0;
//		stat.pr.order = ipolyorder;
//		stat.data$ = strYdata;
//		stat.pr();
//		coeff[0] = stat.pr.a;
//	}
//
//	for( int i = 1; i <= ipolyorder; i++ )
//	{
//		string str;
//		str.Format("stat.pr.b%d", i);
//		LT_get_var(str, &coeff[i]);
//	}
//
//	if  ( (idenom != 0) )
//	{
//		cc.SetLowerIndex(ilower);
//		cc.SetUpperIndex(iupper);
//	}
//
//	return true;
//}

// the following is no longer needed and users should use 
// using NLSF = LabTalk.NLSF;  // Point to the NLSF object

// perform initialization of NLSF and set all control strucutre variables to sensible values
//BOOL initNLSF(NLSFCntrl &control)
//{
//	_LT_Obj
//	{
//		nlsf.init();
//		lstrcpy(control.szYdataName, nlsf.y$);
//		lstrcpy(control.szFuncName, nlsf.func$);
//		lstrcpy(control.szWtDataName, nlsf.w$);
//		control.Tolerance = nlsf.tolerance;
//		control.Confidence = nlsf.conf;
//		control.Prediction = nlsf.pred;
//		control.imaxIter = nlsf.maxIter;		
//		control.idataBegin = nlsf.dataBegin;
//		control.idataEnd = nlsf.dataEnd;
//		control.idataStep = nlsf.dataStep;
//		control.iwType = nlsf.wType;
//
//		for(int ii=0; ii <200; ii++)
//		{
//			control.Par[ii] = 0;
//			control.Err[ii] = 0;
//			control.Dep[ii] = 0;
//			control.ConfIntv[ii]  = 0;
//		}	
//		
//		control.openGraph = true;
//		control.createCurve = true;
//		control.confBands = false;
//	 	control.predBands = false;
//	 	control.createWks = false;
//	 	control.pasteToPlot = true; 
//	 	control.pasteToResults = true; //control.sameXasData = true;
//	 	
//		
//	}
//	return true;
//}
//
// perform NLSF fitting using the NLSFControl structure
//BOOL fitNLSF(NLSFCntrl &control)
//{
//	_LT_Obj
//	{
//		nlsf.msgPrompt = 0;												// suppress error messages!
//		nlsf.func$ = control.szFuncName;								// set fitting function
//		nlsf.fitdata$ = control.szYdataName;							// set dependent variable
//		
//		
//		Dataset dWeight(control.szWtDataName);							// create dWeight dataset variable
//					
//		if( (control.iwType > 0) && (control.iwType < 5) )
//		{ 
//			nlsf.wType = control.iwType; 								// set weight method
//			if( (dWeight.IsValid()) && (nlsf.wType != 0 && nlsf.wType !=2) )
//			{
//				nlsf.w$ = control.szWtDataName;							// set weighting dataset
//			}
//		}
//		
//
//		nlsf.execute("parainit");										// initialize parameters
//		
//		if (control.imaxIter > 0) nlsf.iterateEX(control.imaxIter);
//		else nlsf.iterateEX(100);										// iterate
//		
//		if(control.openGraph)
//		{		
//			if(control.confBands)
//			{
//				nlsf.createCurves("C");									// create confidence bands
//			}
//			if(control.predBands)
//			{
//				nlsf.createCurves("P");									// create prediction bands
//			}
//			if(control.pasteToPlot)
//			{
//				nlsf.pasteParams("Plot");								// create fit label on graph
//			}
//			if(control.createCurve)
//			{
//				nlsf.end(12);											// create fit curve on graph
//			}	
//		}
//		if(control.pasteToResults)
//		{
//			nlsf.pasteParams("Results");								// create results in Results Log
//		}
//		if(control.createWks)
//		{
//			nlsf.paramWKS("Parameters");								// create Parameters worksheet
//		}		
//		
//		// assign nlsf object properties to NLSFCntrl structure
//		control.ChisqDoF = nlsf.chiSqr;
//		control.SSR = nlsf.ssr;
//		control.Correlation = nlsf.cor; 
//		control.COD = nlsf.cod;
//		control.MuFinal = nlsf.muMin;
//		control.Mu = nlsf.mu;
//		control.DerivStep = nlsf.derivStep;
//		control.Tolerance = nlsf.tolerance;
//		control.Confidence = nlsf.conf;
//		control.Prediction = nlsf.pred;
//		control.inPara = nlsf.nPara;
//		control.inParaVary = nlsf.nParaVary;
//		control.imaxIter = nlsf.maxIter;
//		control.idataBegin = nlsf.dataBegin;
//		control.idataEnd = nlsf.dataEnd;
//		control.idataStep = nlsf.dataStep;
//		control.inPoints = nlsf.nPoints;
//		/*
//		control.ixBegin = nlsf.xBegin;
//		control.ixEnd = nlsf.xEnd;
//		control.ixPoints = nlsf.xPoints;
//		*/
//		control.iDOF = nlsf.dof;
//		control.inConstr = nlsf.nConstr;
//		control.inConstrEff = nlsf.nConstrEff;
//		control.iwType = nlsf.wType;
//		
//		// set values of results from fit
//		for(int ii = 1; ii <= nlsf.nPara; ii++)
//		{
//			control.Par[ii-1] = nlsf.p$(ii);
//			control.Err[ii-1] = nlsf.e$(ii);
//			control.Dep[ii-1] = nlsf.d$(ii);
//			control.ConfIntv[ii-1] = nlsf.c$(ii);
//		}
//	}	
//	return true;		
//}


// overloaded fitNLSF function to allow for quick fitting with weighting
//double fitNLSF(string strYData, string strFitFunc, int iwType, string strWtData)
//{	
//	NLSFCntrl control;													// instantiate structure variable
//	initNLSF(control);													// initialize fit 
//	
//	lstrcpy(control.szYdataName, strYData);								// initialize dependent variable
//	lstrcpy(control.szFuncName, strFitFunc);							// initialize fitting function
//	
//	control.iwType = iwType;											// initialize weight type
//	lstrcpy(control.szWtDataName, strWtData); 							// initialize weighting dataset name
//	
//	fitNLSF(control);													// fit
//	return control.ChisqDoF;
//	
//}	

/// Cloud 03/26/07 QA80-9093 v8.0589 INITIAL_PARAMETERS
/// Cloud 02/21/08 PROVIDE_SEGMENT_CENTER
//bool fit_polyline(int iSize, double *px, double *py, int n, vector &slope, vector &ave, vector &cutoff)
bool fit_polyline(int iSize, double *px, double *py, int n, vector &slope, vector &ave, vector &cutoff, vector &xave)
/// End PROVIDE_SEGMENT_CENTER
{
	/// Cloud 12/06/07 FIX_BUG_OF_SECTION_LENGTH
	//if (NULL==px || NULL==py || iSize<n*2)
	if (NULL==px || NULL==py || iSize-1<n)
	/// End FIX_BUG_OF_SECTION_LENGTH
		return FALSE;

	double *pTmpx = px; /// Cloud 02/29/08 PX_NOT_RESTORE_BUG_BREAK_NLFIT_PARAM_INIT_FOR_EXP_DECAY
	double *pTmp = py;
	/// Cloud 12/06/07 FIX_BUG_OF_SECTION_LENGTH
	//int iLength = iSize / n;
	//int iLast = iSize - iLength*(n-1);
	int iLength = (iSize-1) / n;
	int iLast = iSize-1 - iLength*(n-1);
	/// End FIX_BUG_OF_SECTION_LENGTH
	LROptions sLROptions;
	_init_fit_options(sLROptions);
	FitParameter sFitParam[2];
	vector vSlope;
	vector vCut;
	
	for (int ii=0; ii<n-1; ii++)
	{
		/// Cloud 12/06/07 FIX_BUG_OF_SECTION_LENGTH
		//ocmath_polynomial_fit(iLength, px, py, NULL, 1, &sLROptions, sFitParam, 2);
		ocmath_polynomial_fit(iLength+1, px, py, NULL, 1, &sLROptions, sFitParam, 2);
		/// End FIX_BUG_OF_SECTION_LENGTH
		vSlope.Add(sFitParam[1].Value);
		vCut.Add(sFitParam[0].Value);
		px += iLength;
		py += iLength;
	}
	/// Cloud 12/06/07 FIX_BUG_OF_SECTION_LENGTH
	//ocmath_polynomial_fit(iLast, px, py, NULL, 1, &sLROptions, sFitParam, 2);
	ocmath_polynomial_fit(iLast+1, px, py, NULL, 1, &sLROptions, sFitParam, 2);
	/// End FIX_BUG_OF_SECTION_LENGTH
	vSlope.Add(sFitParam[1].Value);
	vCut.Add(sFitParam[0].Value);

	if(slope!=NULL)
		slope = vSlope;

	px = pTmpx; /// Cloud 02/29/08 PX_NOT_RESTORE_BUG_BREAK_NLFIT_PARAM_INIT_FOR_EXP_DECAY
	py = pTmp;
	if(ave!=NULL)
	{
		double dAve;
		for (ii=0; ii<n-1; ii++)
		{
			/// Cloud 12/06/07 FIX_BUG_OF_SECTION_LENGTH
			//dAve = ocmath_d_sum(py, iLength) / iLength;
			dAve = ocmath_d_sum(py, iLength+1) / (iLength+1);
			/// End FIX_BUG_OF_SECTION_LENGTH
			ave.Add(dAve);
			py += iLength;
		}
		/// Cloud 12/06/07 FIX_BUG_OF_SECTION_LENGTH
		//dAve = ocmath_d_sum(py, iLast) / iLast;
		dAve = ocmath_d_sum(py, iLast+1) / (iLast+1);
		/// End FIX_BUG_OF_SECTION_LENGTH
		ave.Add(dAve);
	}

	/// Cloud 02/21/08 PROVIDE_SEGMENT_CENTER
	if(xave!=NULL)
	{
		double dAve;
		for (ii=0; ii<n-1; ii++)
		{
			dAve = ocmath_d_sum(px, iLength+1) / (iLength+1);
			xave.Add(dAve);
			px += iLength;
		}
		dAve = ocmath_d_sum(px, iLast+1) / (iLast+1);
		xave.Add(dAve);
	}
	/// End PROVIDE_SEGMENT_CENTER

	if(cutoff!=NULL)
		cutoff = vCut;

	return TRUE;
}

/// Cloud 04/03/07 ADD_PEAK_HEIGHT_CALCULATION
//double peak_pos(Curve &data, double *pWidth, double *pBaseline, double *pArea, bool bSortDataFirst, int *pErr)
/// Fisher 2008-10-8 BUG_IN_SURFACE_FITTING_INIT_CODES
//double peak_pos(Curve &data, double *pWidth, double *pBaseline, double *pArea, double *pHeight, double *pCentroid, bool bHaveBase, bool bSortDataFirst, int *pErr)
double peak_pos(Curve &data, double *pWidth, double *pBaseline, double *pArea, double *pHeight, double *pCentroid, bool bHaveBase, bool bSortDataFirst, int *pErr, int *pnDir)
/// End BUG_IN_SURFACE_FITTING_INIT_CODES
{
	if (bSortDataFirst)
		sort(data);

	Curve cc(data);
	Dataset dsX;
	int iSize = cc.GetSize();
	cc.AttachX(dsX);
	
	vector vSlope;
	/// Cloud 12/25/07 UPDATE_METHOD_FOR_BASELINE
	//vector vAve;
	/// Cloud 12/06/07 FIX_BUG_OF_SECTION_LENGTH
	//int iSec = 5;
	//if (iSize<10)
		//iSec = iSize / 2;
	//int iSec;
	double dBaseline;
	//if (iSize>5)
		//iSec = 5;
	//else
		//iSec = iSize - 1;
	/// End FIX_BUG_OF_SECTION_LENGTH
	//if (bHaveBase)
		//dBaseline = 0.5 * (cc[0] + cc[iSize-1]);
	//else
		//dBaseline = 0;

	/********************************************************************
		Compare the slopes at start and end of the curve and
		consider the smaller one is the position of the baseline
		(Need improvement)
	 ********************************************************************/
	//if (fit_polyline(iSize, dsX, cc, iSec, vSlope, vAve))
	//{
		////vector<uint> vIndex;
		//vSlope = abs(vSlope);
		////vSlope.Sort(SORT_ASCENDING, TRUE, vIndex);
		//dBaseline = (vSlope[0]<vSlope[iSec-1]) ? vAve[0] : vAve[iSec-1];
		////dBaseline = vAve[vIndex[0]];
	//}
	int nSegs;
	vector<int> vLength(iSize);
	vector vAve(iSize);
	/// Cloud 01/02/08 USE_NOISE_AS_THRESHOLD
	//nSegs = ocmath_cut_to_segments_by_height(iSize, dsX, cc, 0.05, iSize, vLength, vAve);
	double dThr = ocmath_get_noise_level(iSize, cc, dsX);
	
	/// Fisher 05/27/09 QA80-13473 SOME_DATA_HAVE_LARGE_NOISE_LEVEL_WHICH_MAKE_SEGMENT_TESTING_UNABLE_TO_PASS_THREHOLD
	// Default is 0.7 which depends on how we defined that a noised data fits the gaussian model, the larger, the stricter.
	double dd;
	LT_get_var("@NG", &dd);
	dThr *= dd;
	/// End SOME_DATA_HAVE_LARGE_NOISE_LEVEL_WHICH_MAKE_SEGMENT_TESTING_UNABLE_TO_PASS_THREHOLD
	
	nSegs = ocmath_cut_to_segments_by_height(iSize, dsX, cc, dThr, iSize, vLength, vAve);
	/// End USE_NOISE_AS_THRESHOLD
	if (nSegs<0)
		return NANUM;
	vLength.SetSize(nSegs);
	vAve.SetSize(nSegs);
	dBaseline = 0.5 * (vAve[0]+vAve[nSegs-1]);
	/// End UPDATE_METHOD_FOR_BASELINE
	
	/// Jack 16/12/2008 QA80-12756 SUPPLY_OTHER_BASELINE_ESTIMATE_METHODS
	// if the original estimate of baseline is not proper
	// i.e. stay in the middle around mean
	// then try to reestimate
	double dMin, dMax, dMean;
	cc.GetMinMax(dMin, dMax);
	dMean = (dMin+dMax)/2.0; //initilize dmean
	cc.Sum(dMean); 
	if(cc.GetSize()>0)
		dMean = dMean/cc.GetSize(); // calculate mean
	
	double dPercentile = 0.05;
	if( fabs(dBaseline-dMean)<= 0.05*abs(dMean) )// the baseline stay around mean
		dBaseline = (cc[0]+cc[cc.GetSize()-1])/2.0;  // assume using the last and first point as the baseline point 
	if( fabs(dBaseline-dMean)<= 0.05*abs(dMean) )  // another estimate method  use the most distant from mean as base
		dBaseline =  fabs(cc[0]-dMean)>fabs(cc[cc.GetSize()-1]-dMean) ? cc[0] : cc[cc.GetSize()-1];
	///End SUPPLY_OTHER_BASELINE_ESTIMATE_METHODS
	
	cc -= dBaseline;
	
	/// Fisher 2008-10-8 BUG_IN_SURFACE_FITTING_INIT_CODES
	if(pnDir && *pnDir != 0)
	{
		double dMin, dMax;
		cc.GetMinMax(dMin, dMax);
		if(*pnDir > 0)
			cc -= dMin;
		else
			cc -= dMax;
	}
	/// End BUG_IN_SURFACE_FITTING_INIT_CODES
	
	IntegrationResult result;
	int err = ocmath_integrate(dsX, cc, 0, dsX.GetSize()-1, &result);
	
	if (err)
	{
		if (pErr!=NULL)
			*pErr = err;
		return NANUM;
	}
	
	/// Fisher 2008-10-8 BUG_IN_SURFACE_FITTING_INIT_CODES
	if(pnDir && *pnDir == 0)
	{
		*pnDir = result.nPeakDir;
	}
	/// End BUG_IN_SURFACE_FITTING_INIT_CODES	
	
	double dArea = result.Area;
	
	if (pWidth!=NULL)
	{
		/// Cloud 12/18/07 USE_AREA_AND_HEIGHT_MORE_ACCURATE
		//*pWidth = result.dxPeak;
		//if (NANUM==*pWidth)
		//{
		/// ML 12/19/2007 NEGATIVE_WIDTH_BUG
		//*pWidth = dArea / (1.2533*result.yPeak);
		*pWidth = fabs(dArea / (1.2533*result.yPeak));
		/// end NEGATIVE_WIDTH_BUG
		//dArea *= 2;
		//}
		/// End USE_AREA_AND_HEIGHT_MORE_ACCURATE
	}
	
	if (pHeight!=NULL)
		*pHeight = result.yPeak;
	
	if (pBaseline!=NULL)
		*pBaseline = dBaseline;
	
	if (pArea!=NULL)
		*pArea = dArea;

	if (pCentroid!=NULL)
		*pCentroid = result.xCentroid;

	return result.xPeak;
}

//int init_sigmoidal(Curve &data, vector vX, vector vY, double &bottom, double &top, double &center, double &slope)
//{
	//sort(data);
	//bottom = min(vY);
	//top = max(vY);
	//center = xaty50(data);
	//double xmin;
	//double xmax;
	//vX.GetMinMax(xmin, xmax);
	//int iSize = vX.GetSize();
	//int iCenter = (center-xmin) / (xmax-xmin) * iSize;
	//int iInterval = iSize / 20;
	//slope = (vY[iCenter+iInterval] - vY[iCenter-iInterval]) / (vX[iCenter+iInterval] - vX[iCenter-iInterval]);
	//
	//return 0;
//}
double sig_inflex(Curve &data, double &slope, double *pLeft, double *pRight, bool bSortDataFirst)
{
	if (bSortDataFirst)
		sort(data);
	
	Curve cc(data);
	int iSize = data.GetSize();
	Dataset dsX;
	cc.AttachX(dsX);
	
	int iSec = 5;
//	deal with the curve with few points
	if (iSize/iSec < 2)
		iSec = iSize / 2;
	
	vector vSlope;
	vector vAve;
	if (!fit_polyline(iSize, dsX, cc, iSec, vSlope, vAve))
		return NANUM;
	
	vector<uint> vIndex;
	vector vTmp;
	///Sandy 2009-1-8 SPEED_UP_VECTOR_ABS
	//vTmp = abs(vSlope);
	vTmp = vSlope;
	vTmp.Abs();
	//end SPEED_UP_VECTOR_ABS
	
	vTmp.Sort(SORT_DESCENDING, TRUE, vIndex);
	slope = vSlope[vIndex[0]];
	
//  Estimate the index of the inflexion using the center of the line with the largest slope
	int iInflex = (vIndex[0] + 0.5) / iSec * iSize;
	double dLeft;
	double dRight;
	if (iInflex < iSize/2)
	{
		dRight = vAve[iSec-1];
		dLeft = 2 * cc[iInflex] - dRight;
	}
	else
	{
		dLeft = vAve[0];
		dRight = 2 * cc[iInflex] - dLeft;
	}
	
	if (pLeft!=NULL)
		*pLeft = dLeft;
	
	if (pRight!=NULL)
		*pRight = dRight;

	return dsX[iInflex];
}

/// Cloud 01/14/08 REPRESENT_SCALE_BY_EXPONENT_FOR_HUGE_NUMBER
//double get_exponent(vector &x_data, vector &y_data, double *pBase, double *pScale, bool bSortDataFirst)
double get_exponent(vector &x_data, vector &y_data, double *pBase, double *pScaleExp, int *pScaleSign, bool bSortDataFirst)
/// End REPRESENT_SCALE_BY_EXPONENT_FOR_HUGE_NUMBER
{
	if (bSortDataFirst)
	{
		vector<uint> vi;
		x_data.Sort(SORT_ASCENDING, TRUE, vi);
		y_data.Reorder(vi);
	}

	int iNum = x_data.GetSize();
	/// Cloud 12/06/07 FIX_BUG_OF_SECTION_LENGTH
	//int iSec = 5;
	//if (iNum/iSec < 2)
		//iSec = iNum / 2;
	int iSec;
	if (iNum>5)
		iSec = 5;
	else
		iSec = iNum - 1;
	/// End FIX_BUG_OF_SECTION_LENGTH
	
	vector vk;
	/// Cloud 02/21/08 REMOVE_0_FROM_SLOPES
	/*
	fit_polyline(iNum, x_data, y_data, iSec, vk);
	if (vk[0]==0)
		vk[0] = 1e-5 * vk[iSec-1];

	double dExp = log(abs(vk[iSec-1]/vk[0])) / (x_data[iNum-1]-x_data[0]);
	*/
	vector vxave;
	fit_polyline(iNum, x_data, y_data, iSec, vk, NULL, NULL, vxave);
	vector vkabs;
	///Sandy 2009-1-8 SPEED_UP_VECTOR_ABS
	//vkabs = abs(vk);
	vkabs = vk;
	vkabs.Abs();
	//end SPEED_UP_VECTOR_ABS
	ocmath_replace_data_in_vector(vkabs.GetSize(), vkabs, 0, NANUM);
	matrix mm(iSec, 2);
	mm.SetColumn(vkabs, 0);
	mm.SetColumn(vxave, 1);
	mm.RemoveEmptyRows(false);
	mm.GetColumn(vkabs, 0);
	mm.GetColumn(vxave, 1);
	double dMin, dMax;
	int iMin, iMax;
	vkabs.GetMinMax(dMin, dMax, &iMin, &iMax);
	double dExp = log(dMax/dMin) / (vxave[iMax]-vxave[iMin]);
	/// End REMOVE_0_FROM_SLOPES
	/// Cloud 01/14/08 REPRESENT_SCALE_BY_EXPONENT_FOR_HUGE_NUMBER
	//double y2 = exp(dExp * x_data[iNum-1]);
	//double y1 = exp(dExp * x_data[0]);
	double y2 = exp(dExp * (x_data[iNum-1]-x_data[0]));
	double y1 = 1;
	double dBase = (y_data[0]*y2 - y_data[iNum-1]*y1) / (y2 - y1);
	double dScale = (y_data[iNum-1] - y_data[0]) / (y2 - y1);
	int nSign = (dScale>0) ? 1 : -1;
	/// End REPRESENT_SCALE_BY_EXPONENT_FOR_HUGE_NUMBER

	if (pBase!=NULL)
		*pBase = dBase;
	/// Cloud 01/14/08 REPRESENT_SCALE_BY_EXPONENT_FOR_HUGE_NUMBER
	//if (pScale!=NULL)
		//*pScale = dScale;
	if (pScaleExp!=NULL)
		*pScaleExp = log(abs(dScale)) - dExp * x_data[0];
	if (pScaleSign!=NULL)
		*pScaleSign = nSign;
	/// End REPRESENT_SCALE_BY_EXPONENT_FOR_HUGE_NUMBER
	return dExp;
}


/// Jack  QA-12334 10/09/2008 ADD_ONE_PEAK_MATRIX_DATA_INITILIZE_CODE
// for usage in one peak Matrix data initial code
double peak_pos_without_baseline(Curve &data, double *pWidth, double *pArea, double *pHeight)
{
	Curve cc(data);
	Dataset dsX;
	int iSize = cc.GetSize();
	cc.AttachX(dsX);
	
	IntegrationResult result;
	int err = ocmath_integrate(dsX, cc, 0, dsX.GetSize()-1, &result);
	if(err)
	{
		return NANUM;
	}
	
		double dArea = result.Area;
	
	if (pWidth!=NULL)
	{
		*pWidth = fabs(dArea / (1.2533*result.yPeak));

	}
	
	if (pHeight!=NULL)
		*pHeight = result.yPeak;
	

	if (pArea!=NULL)
		*pArea = dArea;
	
	return result.xPeak;
	
}
/// End ADD_ONE_PEAK_MATRIX_DATA_INITILIZE_CODE


/// Cloud 07/17/07 ADD_FFT_FOR_COMPATIBILITY
/// Fisher 08/19/08 QA80-12043
/*
int FFT(const vector<complex>& vc, vector<complex>& vr)
{
	vr = vc;
	return fftw_fft_complex(vr.GetSize(), vr, FFT_FORWARD);
}
*/
/// End 08/19/08 QA80-12043
/// End ADD_FFT_FOR_COMPATIBILITY
/// End INITIAL_PARAMETERS

/// EJP 2008-04-08 v8.0840 QA80-11237 MAP_ROI_MATRIX_RECT_TO_IMAGE_RECT
//---- CPY 12/6/06 QA70-12736 LT_NEED_MORE_CNTRL_ON_ALLOWING_OC_FUNCS
//#pragma labtalk(1) // Let the following functions be called from LabTalk
#pragma labtalk(2)
//----

#define SWAP_DOUBLE(_d1,_d2) {double _dTmp = _d1; _d1 = _d2; _d2 = _dTmp;}
int get_active_ROI_rect(int& top, int& left, int& bottom, int& right)
{
	MatrixLayer ml = Project.ActiveLayer();
	if( !ml.IsValid() )
		return 1; // matrix not found

	GraphObject grObj;
	grObj = ml.GraphObjects("ROI");
	if( !grObj.IsValid() )
		return 2; // ROI object not found

	//------ Folger 08/11/08 QA80-11995 COMPARE_LOCALIZED_GRAPH_OBJECT_STRING
	//string strObjectType = grObj.GetObjectType();
	//if( strObjectType.Compare("Rectangle") )
	int nGrObjType = GROBJ_TN_UNKNOWN;
	grObj.GetObjectType(&nGrObjType);
	if( nGrObjType != GROBJ_TN_RECT )
	//------
		return 3; // ROI object is not a rectangle

	MatrixObject matObj(ml, -1);
	if( !matObj.IsValid() )
		return 4; // Matrix not found

	// Get matrix dimensions
	int nCols = matObj.GetNumCols();
	int nRows = matObj.GetNumRows();

	// Get matrix X min,max and Y min,max
	double dXMin, dXMax, dYMin, dYMax;
	matObj.GetXY(dXMin, dYMin, dXMax, dYMax);

	// Check the X min,max and calculate the X step
	if( dXMin > dXMax )
		SWAP_DOUBLE(dXMin, dXMax);
	double dXStep = (dXMax - dXMin) / (double)nCols;
	dXMin -= (dXStep / 2.0);
	dXMax += (dXStep / 2.0);

	// Calculate the left cell index
	double dLeft = grObj.X - grObj.DX / 2; // object center minus half it's width
	if( dLeft < dXMin )
		dLeft = dXMin;
	left = (dLeft - dXMin) / dXStep;

	// Calculate the right cell index
	double dRight = grObj.X + grObj.DX / 2; // object center plus half it's width  
	if( dRight > dXMax )
		dRight = dXMax;
	right = (dRight - dXMin) / dXStep;

	// Check the Y min,max and calculate the Y step
	if( dYMin > dYMax )
		SWAP_DOUBLE(dYMin, dYMax);
	double dYStep = (dYMax - dYMin) / (double)nRows;
	dYMin -= (dYStep / 2.0);
	dYMax += (dYStep / 2.0);

	// Calculate the top cell index
	double dTop = grObj.Y - grObj.DY / 2; // object center minus half it's height
	if( dTop < dYMin )
		dTop = dYMin;
	top = (dTop - dYMin) / dYStep;

	// Calculate the bottom cell index
	double dBottom = grObj.Y + grObj.DY / 2; // object center minus half it's height
	if( dBottom > dYMax )
		dBottom = dYMax;
	bottom = (dBottom - dYMin) / dYStep;

	return 0;
}
/// end MAP_ROI_MATRIX_RECT_TO_IMAGE_RECT
/////////////////////////////////////////////////////////////////////////////////
//CPY 11/16/08 useful functions for Set Col Value to process strings in columns
//--------CPY 10/29/09 QA70-14563 SUPPORT_OC_FUNC_NO_NEED_TO_BE_USED_AS_LT_CMD
//	replace all "#pragma labtalk(3," from here to "#pragma labtalk(3,"
//--------
/////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////
// Begin LabTalk Set Values dialog functions
/////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////
#pragma labtalk(3, Statistics) //3=can be used as function only, cannot be used as LT command

//using a vector copy since LT might pass in dataset if by ref and it might be text&numeric, or float
//so cannot use pass by ref, this applies to all functions that can receive a dataset as argument
double Min(vector vv)
{
	double y1 = NANUM,y2=NANUM;
	vv.GetMinMax(y1, y2);
	return y1;
}
double Max(vector vv)
{
	double y1 = NANUM,y2=NANUM;
	vv.GetMinMax(y1, y2);
	return y2;
}
double Mean(vector vv)
{
	double dSum, dMean = NANUM, dSD;
	int N;
	ocmath_basic_summary_stats(vv.GetSize(), vv, &N, &dMean, &dSD, NULL, NULL, &dSum);
	return dMean;
}
double Total(vector vv)
{
	double dd = NANUM;
	vv.Sum(dd);
	return dd;
}

/*
|-
! 0 (default)
| Counts the number of elements in the vector, size of the vector
|-
! 1
| Counts the number of Numeric Values
|-
! 2
| Counts the number of Missing Values
|-
*/
int Count(vector vv, int option = 0)
{
	int nSize = vv.GetSize();
	if(0==option)
		return nSize;
    int numMissing=ocmath_count(NANUM, vv.GetSize(), vv);
    if(2 == option)
    	return numMissing;
    if(1 == option)
    	return nSize - numMissing;
    
    //not supported yet
    return -1;
}

static double _std_dev(vector& vv, int option)
{
	int nn = 0;
	double mean = NANUM;
	double SD = NANUM;
	int MomentDenFlag = 0 == option? DS_SAS1_DOF : DS_SAS2_N;
	
	ocmath_basic_summary_stats(vv.GetSize(), vv, &nn, &mean, &SD, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, MomentDenFlag);
	
	return SD;	
}

double StdDev(vector vv)
{
	return _std_dev(vv, 0);	
}

double StdDevP(vector vv)
{
	return _std_dev(vv, 1);
}

// nInterpolation here same as the option Computation Control -> Interpolation of Quantiles on Statistics on Column dialog.
double Median(vector<double> vdData, int nInterpolation = 0)
{
	QuantileOptions opt;
	opt.Interpolate = nInterpolation;	
	opt.Median = true;
	
	QuantileResults res;
	if(STATS_NO_ERROR == ocmath_quantiles(vdData, vdData.GetSize(), &opt, &res))
		return res.Median;
	return NANUM;
}

/// Sam 12/17/09 QA81-14596 v8.1.11.89  MOVING_SLOPE_AND_RMS_FUNCTIONS
// Return the slope of the linear regression line through each npt data points in two vectors. 
vector movslope(vector vX, vector vY, int nPt = -1/*, int& nErr = 0*/)
{	
	vector vS;
	
	if(vX.GetSize() != vY.GetSize() || vX.GetSize() < nPt)
	{	
		//nErr = 1; // Vectors X and Y should have a same size and be no less than nPt.
		return vS;
	}
	
	if (nPt == -1 || nPt == vX.GetSize())
	{
		vS.SetSize(1);
		vS[0] = slope(vX, vY);
		return vS;
	}
	
	if(mod(nPt, 2) == 0)
		nPt = nPt + 1;
		
	if(nPt < 3)
	{
		//nErr = 2; // nPt should be an odd number and greater than 2.
		return vS;
	}
	
	vS.SetSize(vX.GetSize());
	vS=NANUM;
	
	vector vXtemp, vYtemp;
	for(int nI = 0; nI <= vX.GetSize() - nPt; nI++)
	{
		vX.GetSubVector(vXtemp, nI, nI+nPt-1);
		vY.GetSubVector(vYtemp, nI, nI+nPt-1);
		
		int nIndex = nI + floor(nPt/2);

		vS[nIndex] = slope(vXtemp, vYtemp);
	}
	
	//nErr = 0; // No error
	return vS;
}

// Return the slope of the linear regression line through data points in two vectors. 
static double slope(const vector& vX, const vector& vY, int& nErr = 0)
{
	double dSp = NANUM;
	if(vX.GetSize() != vY.GetSize())
	{
		nErr = 1; //Vectors X and Y should have a same size.
		return dSp;
	}
	
	double dMx, dMy;
	ocmath_basic_summary_stats(vX.GetSize(), vX, NULL, &dMx);
	ocmath_basic_summary_stats(vY.GetSize(), vY, NULL, &dMy);
	
	double 	dSx = 0, dSy = 0;
	int nPt = vX.GetSize();
	for (int nI = 0; nI < nPt; nI++)
	{
		dSx = dSx + (vX[nI] - dMx)*(vX[nI] - dMx);
		dSy = dSy + (vX[nI] - dMx)*(vY[nI] - dMy);
	}
	
	double 	eps = 2.2e-16;
	if (dSx < eps)
		dSp = NANUM;
	else
		dSp = dSy/dSx;
	nErr = 0; //No error
	return dSp;
}

//Return the root mean square of a vector.
double rms(vector vX)
{
	double dRms;
	int nCount;  //Number of non missing values
	
	nCount = vX.GetSize() - ocmath_count(NANUM, vX.GetSize(), vX);
	double dSumsq;
	ocmath_basic_summary_stats(vX.GetSize(), vX, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dSumsq);
	dRms = sqrt(dSumsq/nCount);
	
	return dRms;
}
/// End MOVING_SLOPE_AND_RMS_FUNCTIONS
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
#pragma labtalk(3, String)
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////

int Len(string str)
{
	return str.GetLength();
}

string Upper(string str)
{
	str.MakeUpper();
	return str;
}
string Lower(string str)
{
	str.MakeLower();
	return str;
}


//option = 
		//0 (default)  Trim leading and trailing spaces  
		//1  Trim all spaces, including those inside the string 
string Trim(string str, int option = 0)
{
	if(0 == option)
	{
		str.TrimLeft();
		str.TrimRight();
	}
	else
	{
		str.Remove(' ');
	}
	return str;
}

//Returns TRUE if both strings are an exact match (case and length) 
BOOL Exact(string str1, string str2)
{
	return 0 == str1.Compare(str2);
}

//Replace n2 characters in string1 starting at n1th position with string2. String2 may be different length than n2. 
string Replace(string str1, int n1, int n2, string str2)
{
	n1--; // 1 offset
	if(n1 < 0 || n2 <=0)
		return "";
	
	str1.Delete(n1, n2);
	str1.Insert(n1, str2);
	return str1;
}

//Substitute string3 with string2 when found in string1. n=0 Substitute all, otherwise only nth instance.
string Substitute(string str1, string str2, string str3, int n = 0)
{
	//int nn = str1.Replace(str3, str2);
	int start = 0;
	int instance = 1;
	while( start < str1.GetLength() )
	{
		int ii = str1.Find(str3, start);		
		if( ii < 0 ) // not found
			break;
		
		if(0 == n || 0 != n && n == instance)
		{		
			str1 = Replace(str1, ii+1, str3.GetLength(), str2); // should use ii+1 here since index of Replace offset is 1
			start = ii + str2.GetLength(); 	
		}
		else
		{
			start = ii + str3.GetLength();
		}
		instance++;
	}
	
	return str1;
}

///Sophy 9/15/2009 ADD_TOKEN_FUNCTION_ON_SET_COLUMNS_VALUE_DIALOG_MENU
//Get the Nth token using specified delimiter from a string, iDelimiter = 0 for any space charaters, 124 for '|', 32 for ' '(one white space)
string Token(string str, int iToken, int iDelimiter = 0)
{
	return str.GetToken(iToken - 1, iDelimiter);
}
///end ADD_TOKEN_FUNCTION_ON_SET_COLUMNS_VALUE_DIALOG_MENU
string Right(string str, int n)
{
	return str.Right(n);
}

string Left(string str, int n)
{
	return str.Left(n);
}

//Returns n2 characters of string starting from n1th position (numbered from 1) 	
string Mid(string str, int n1, int n2)
{
	// offset is 1
	n1--; 
	
	if(n1 < 0 || n2 <= 0)
		return "";
	
	if(n1 + n2 >= str.GetLength())
		n2 = str.GetLength() - n1; //to end
	
	return str.Mid(n1, n2);	
}

//return location of a substring in a string. The search is case-sensitive. To do case insensitive search, use the Search function. 
int Find(string str1, string str2, int StartPos = 1)
{
	StartPos--; // StartPos is LT offset from 1
	
	int nFind = str1.Find(str2, StartPos); // return -1 if str2 is not found. 
	if(nFind < 0)
		return -1;
	
	return nFind+1; // change to 1 offset   
}

//return the location of a substring in a string. the search is NOT case-sensitive. To do case-sensitive search, use the Find function.
int Search(string str1, string str2, int StartPos = 1)
{
	str1.MakeUpper();
	str2.MakeUpper();
	StartPos--; // StartPos is LT offset from 1
	
	int nFind = str1.Find(str2, StartPos); // return -1 if str2 is not found. 
	if(nFind < 0)
		return -1;
	
	return nFind+1; // change to 1 offset
}

//return the beginning location of a substring in a string by using a string pattern containing wildcard characters like * or ?. 
int MatchBegin(string str1, string str2, int StartPos = 1, bool bCaseSensitive = false)
{
	StartPos--;
	
	if(StartPos < 0)
		return -1;
	
	if(StartPos > 0)
		str1 = str1.Right(strlen(str1)-StartPos);
	
	str2 = "*" + str2 + "*"; // intput str2 here is partial pattern 

	int nBegin;
	if(pattern_match(str2, str1, bCaseSensitive, &nBegin))
		nBegin += StartPos;
	else
		return -1;
	
	return nBegin + 1; // LT offset is 1
}

//return the ending location of a substring in a string by using a string pattern containing wildcard characters like * or ?. 
int MatchEnd(string str1, string str2, int StartPos = 1, bool bCaseSensitive = false)
{
	StartPos--;
	
	if(StartPos < 0)
		return -1;
	
	if(StartPos > 0)
		str1 = str1.Right(strlen(str1)-StartPos);

	str2 = "*" + str2 + "*"; // intput str2 here is partial pattern 
	
	int nEnd;
	if(pattern_match(str2, str1, bCaseSensitive, NULL, &nEnd))	
		nEnd += StartPos;
	else
		return -1;
	
	return nEnd + 1;// LT offset is 1
}

bool Compare(string str1, string str2, bool bCaseSensetive = true)
{
	if(bCaseSensetive)
		return 0 == str1.Compare(str2);
	else
		return 0 == str1.CompareNoCase(str2);
}

// number should in 0-255 else return "--"
string Char(int number)
{
	string str = NANUM;// will display -- in column for error
	
	if(number >=0 && number <= 255)
	{
		str.Format("%c", number);
	}
	return str;
}

int Code(string str)
{
	if(str.IsEmpty())
		return 0;

	char ch = str.GetAt(0);
	return (int)ch;
}

//LabTalk formatting string, "" will use @SD significant digits. Use "*" for Origin's global setting. 
string Format(double data, string fmt = "*")
{
	return ftoa(data, fmt);
}

///////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
#pragma labtalk(3, Date and Time)
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////

int Year(double date)
{
	SYSTEMTIME st;
	if(!JulianDateToSystemTime(&date, &st))
		return -1;
	
	return st.wYear;
}
int Month(double date)
{
	SYSTEMTIME st;
	if(!JulianDateToSystemTime(&date, &st))
		return -1;
	
	return st.wMonth;
}

static bool _get_Jan1_julian_date(int nYear, double& date)
{
	SYSTEMTIME stJan1;
	stJan1.wYear = nYear;
	stJan1.wMonth = 1;
	stJan1.wDay = 1;
	
	if(!SystemTimeToJulianDate(&date, &stJan1))
		return false;	
	return true;
}

int Day(double date, int option = 1)
{
	SYSTEMTIME st;
	if(!JulianDateToSystemTime(&date, &st))
		return -1;
	
	if(1 == option)
	{
	return st.wDay;
	}
	else if(2 == option)
	{
		double dateJan1;
		if(!_get_Jan1_julian_date(st.wYear, dateJan1))
			return -2;
		return (int)date - (int)dateJan1 + 1;
	}
	
	return -3; // invalid option value	
}

static double _add_day_info_for_time_value(double time)
{
	SYSTEMTIME st;
    GetSystemTime(&st); 
    
    double dDate;
    SystemTimeToJulianDate(&dDate, &st);
    dDate = (int)dDate + time; 	    
	return dDate;
}

int Hour(double value)
{
	double date = value;
	if( value < 1 ) // only time value not included year, month and day info
	{
	    date = _add_day_info_for_time_value(value);
	}
	
	SYSTEMTIME st;
	if(!JulianDateToSystemTime(&date, &st))
		return -1;
	
	return st.wHour;
}


int Minute(double value)
{
	double date = value;
	if( value < 1 ) // only time value not included year, month and day info
	{
	    date = _add_day_info_for_time_value(value);
	}

	SYSTEMTIME st;
	if(!JulianDateToSystemTime(&date, &st))
		return -1;
	
	return st.wMinute;
}

double Second(double value)
{
	double date = value;
	if( value < 1 ) // only time value not included year, month and day info
	{
	    date = _add_day_info_for_time_value(value);
	}

	SYSTEMTIME st;
	if(!JulianDateToSystemTime(&date, &st))
		return -1;
	
	double dSecond = st.wSecond + st.wMilliseconds / 1000.0 ;
	return dSecond;
}

double Now()
{
	time_t ltime;
	time( &ltime );
	
	TM* tmNow = localtime( &ltime );

	SYSTEMTIME SysTime; 
	tm_to_systemtime(tmNow, &SysTime);	
    
    double date;
    SystemTimeToJulianDate(&date, &SysTime); 
    return date;
}

double Today()
{
	return (int)Now();
}

#define JULIAN_DAYS_1900	(2415019)

string MonthName(double dateOrIndex, int nAbbriv = 3)
{
	int nMonth;
	if(dateOrIndex >= JULIAN_DAYS_1900)
		nMonth = Month(dateOrIndex);
	else
		nMonth = nint(dateOrIndex);
		
	char szTemp[MAXLINE];
	if(get_month_name(nMonth, nAbbriv, szTemp, MAXLINE))
		return szTemp;
	
	return "";
}

///Army 08/06/09 QA81-14051 V8.1076 ADD_WEEKDAYNAME_FUNCTION
//nOption = 0, return 0 to 6 (Sunday to Saturday)
		//= 1, return 1 to 7 (Sunday to Saturday)
		//= 2, return 0 to 6 (Monday to Sunday)
		//= 3, return 1 to 7 (Monday to Sunday)
string WeekDayName(double dateOrIndex, int nAbbriv = 3, int option = 3)
{
	int nDayOfWeek;
	if(dateOrIndex >= JULIAN_DAYS_1900)
		nDayOfWeek = WeekDay(dateOrIndex);
	else
	{
		nDayOfWeek = nint(dateOrIndex);	
		
		switch(option)
		{
			case 0:
			default:
				// do nothing
				break;
				
			case 1:
				nDayOfWeek--;
				break;	
				
			case 2:
				nDayOfWeek++;
				if(7 == nDayOfWeek)
					nDayOfWeek = 0;
				break;
				
			case 3:
				if(7 == nDayOfWeek)
					nDayOfWeek = 0;
				break;
		}
	}
	char szTemp[MAXLINE];
	if(get_weekday_name(nDayOfWeek, nAbbriv, szTemp, MAXLINE))
		return szTemp;
	
	return "";
}
///End ADD_WEEKDAYNAME_FUNCTION

//option:
// 0 = 2 digits like 99, 08
// 1 = '99, '08 form
// 2 = 1999, 2008
// dateOrYear must be either Date or 4 digit year 
string YearName(double dateOrYear, int option = 1)
{
	int nn;
	if(dateOrYear >= JULIAN_DAYS_1900)
		nn = Year(dateOrYear);
	else
		nn = nint(dateOrYear);
	
	if(2 == option)
		return (string) nn;
	
	int nYearSince1900 = nn - 1900;
	int n2digit = nYearSince1900 % 100;
	char szTemp[20];
	int_to_fixed_str(n2digit, 2, szTemp);	 
	if(0 == option)
		return szTemp;
	
	
	return "'" + szTemp;
}


//nOption = 0, return 0 to 6 (Sunday to Saturday)
		//= 1, return 1 to 7 (Sunday to Saturday)
		//= 2, return 0 to 6 (Monday to Sunday)
		//= 3, return 1 to 7 (Monday to Sunday)
int WeekDay(double date, int option = 0)
{
	SYSTEMTIME st;
	if(!JulianDateToSystemTime(&date, &st))
		return -1;
	
	int nDayOfWeek = st.wDayOfWeek; // this day is from 0 to 6 (Sunday to Saturday)
	if( nDayOfWeek < 0 || nDayOfWeek > 6)
		return -1;
	
	switch(option)
	{
	case 0:
	default:
		// do nothing
		break;
		
	case 1:
		nDayOfWeek++;
		break;
		
	case 2:
		nDayOfWeek--;
		if(nDayOfWeek < 0)
			nDayOfWeek = 6;
		break;
		
	case 3:
		if(0 == nDayOfWeek)
			nDayOfWeek = 7;
		break;
	}
	
	return nDayOfWeek;
}

//Returns 1 to 53, indicating the calendar week number of the year, 
//where option = 1 means week starts on Sunday and option = 2 means week starts on Monday 
int WeekNum(double date, int option = 1)
{
	SYSTEMTIME st;
	if(!JulianDateToSystemTime(&date, &st))
		return -1;	
	
	// get Julian date from Jan1 (d0) of that year of date(d1)
	double dateJan1;
	if(!_get_Jan1_julian_date(st.wYear, dateJan1))
		return -2;		
	
	// take integer part of (d1-d0) and divide it by 7
	int nn = (int)date - (int)dateJan1;
	int nWeekNum = nn / 7 + 1;
	
	// option = 1 means week starts on Sunday and option = 2 means week starts on Monday
	int nOptionWeekDay = 1 == option? 1: 3;
	int nWeekDay = WeekDay(date, nOptionWeekDay);
	int nWeekDayOfJan1 = WeekDay(dateJan1, nOptionWeekDay);
	if( nWeekDay < nWeekDayOfJan1 )// if the week day of input date large than Jan 1, then should be next week
		nWeekNum++;
	
	return nWeekNum;
	
}

double Quarter(double date)
{
	int month = Month(date);
	if( -1 == month )
		return -1;
	
	if( month >= 1 && month <= 3)
		return 1;
	
	if( month >= 4 && month <= 6)
		return 2;
	
	if( month >= 7 && month <= 9)
		return 3;
	
	if( month >= 10 && month <= 12)
		return 4;
	
	return -1;
}
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
#pragma labtalk(3, Math)
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////

double asinh(double x)
{	
	return log (x + sqrt(1 + pow(x,2)));
}

double acosh(double x)
{
	if ( fabs(x) < 1 )
		return NANUM;
	return 2 * log (sqrt((x+1)/2.0) + sqrt((x-1)/2.0));
}

double atanh(double x)
{
	if ( fabs(x) >= 1 )
		return NANUM;
	return (log (1+x) - log (1-x))/2.0;
}

///Sophy 6/26/2009 ADD_MORE_TRIANGLE_HYPERBOLIC_AND_INV_HYPERBOLIC_FUNCTIONS
double cot(double x)
{
	if ( is_equal(tan(x), 0) )
		return NANUM;
	return 1/tan(x);
}

double acot(double x)
{
	return PI/2 - atan(x);
}

//coth(x) = cosh(x)/sinh(x) = (exp(x) + exp(-x))/(exp(x) - exp(-x))
double coth(double x)
{
	if ( is_equal(sinh(x), 0) )
		return NANUM;
	return cosh(x)/sinh(x);
}

double acoth(double x)
{
	if ( fabs(x) <= 1 || is_equal(x, 0) )
		return NANUM;
	return atanh(1/x);
}

//if named as sec, will conflict with labtalk timer command sec -i.
double secant(double x)
{
	if ( is_equal(cos(x), 0) )
		return NANUM;
	return 1/cos(x);
}

double csc(double x)
{
	if ( is_equal(sin(x), 0) )
		return NANUM;
	return 1/sin(x);
}

double asec(double x)
{
	if ( fabs(x) < 1 )
		return NANUM;
	return acos(1/x);
}

double acsc(double x)
{
	if ( fabs(x) < 1 )
		return NANUM;
	return asin(1/x);
}

//sech(x) = 1/cosh(x);
double sech(double x)
{
	return 1/cosh(x);
}

//csch(x) = 1/sinh(x)csch
double csch(double x)
{
	if ( is_equal(sinh(x), 0) )
		return NANUM;
	return 1/sinh(x);
}

//asech(x) = ln(sqrt(1/x^2 - 1) + 1/x)
double asech(double x)
{
	if ( x > 1 || x <= 0 )
		return NANUM;
	return log(sqrt(1/pow(x, 2) - 1) + 1/x);
}

//acsch(x) = ln(sqrt(1 + 1/x^2) + 1/x);
double acsch(double x)
{
	if ( is_equal(x, 0) )
		return NANUM;
	return log(sqrt(1/pow(x, 2) + 1) + 1/x);
}

///end ADD_MORE_TRIANGLE_HYPERBOLIC_AND_INV_HYPERBOLIC_FUNCTIONS

/// Iris Matrix& not support now, so commented out the following functions 
/*
double MDeterm(Matrix& mat)
{
	vector<int> vPivot;
	vPivot.SetSize(mat.GetNumRows());
	
	double derf;
	int dete;
	NagError fail;
	nag_real_lu(mat.GetNumRows(), mat, mat.GetNumCols(), vPivot, &derf, &dete, &fail);
	if(fail.code == NE_NOERROR)
		return (derf * pow(2, dete));   
	
	return NANUM;
}

Matrix MInverse(Matrix& mat)
{
	matrix matResult = mat;
	if(0 == matResult.Inverse()) // returns 0 on success
		return matResult;
	
	return mat; // if failed, return input matrix
}

Matrix MMult(Matrix& mat1, Matrix& mat2)
{
	Matrix mResult;
	mResult.SetSize(mat1.GetNumRows(), mat1.GetNumCols());
	mResult = mat1;
	
	if(mResult.DotMultiply(mat2))
		return mResult;	
	
	mResult = NANUM;
	return mResult;
}
*/

int sign(double x)
{
	if( x > 0 )
		return 1;
	else if( x < 0 )
		return -1;
	else
		return 0;	
}

double Degrees(double angle)
{
	return (180 * angle / PI);
}

double Radians(double angle)
{
	return angle * pi / 180;
}

double Distance(double x1, double y1, double x2, double y2)
{
	return sqrt( pow((x1-x2), 2) + pow((y1-y2), 2) );
}

double Distance3D(double x1, double y1, double z1, double x2, double y2, double z2)
{
	return sqrt( pow((x1-x2), 2) + pow((y1-y2), 2) + pow((z1-z2), 2));	
}

//in degrees if unit=1 or radians if unit=0, default is radians
double Angleint1(double x1, double y1, double x2, double y2, int unit = 0, int direction = 0)
{
	double radians = NANUM;
	
	if( is_equal(x1, x2) && is_equal(y1, y2) )
		return NANUM;
	
	/// Iris 4/13/2009 FIX_ANGLE_FUNC_BUG_REPORTED_BY_GREG
	//if( is_equal(x1, x2) )
			//radians = PI/2;
	//else
	///end FIX_ANGLE_FUNC_BUG_REPORTED_BY_GREG	
	{	
		if(0 == direction)
		{
			/// Iris 4/13/2009 FIX_ANGLE_FUNC_BUG_REPORTED_BY_GREG
			if( is_equal(x1, x2) )
				radians = PI/2;	
			else
			///end FIX_ANGLE_FUNC_BUG_REPORTED_BY_GREG
			{		
				double k = (y1 - y2) / (x1 - x2);		
		        radians = atan(k);	
			}
		}
		else// vector angle
		{				
			/// Iris 4/13/2009 FIX_ANGLE_FUNC_BUG_REPORTED_BY_GREG
			if( is_equal(x1, x2) )
				radians = PI/2;	
			else
			///end FIX_ANGLE_FUNC_BUG_REPORTED_BY_GREG
			{
				radians = acos( (x2 - x1) / sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) ) );
			}
	 		
	 		if( y2 - y1 < 0 )
	 			radians = 2*PI - radians;
		}
	}
	
	if(1 == unit)
		return Degrees(radians);
	return radians;	
}

double Angleint2 (double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, int unit = 0, int direction = 0)
{
	if( 1 == direction ) // angle for two vectors
	{
		double ang1 = Angleint1(x1, y1, x2, y2, unit, 1);
		double ang2 = Angleint1(x3, y3, x4, y4, unit, 1);
		double ang = ang1 - ang2; 
		
		if( ang <= 0 )
		{
			if( 0 == unit )
				ang = 2*PI + ang; 
			else	
				ang = 360 + ang;
		}		
		return ang;
	}
	
	double radians = NANUM;
	
	if( is_equal(x1, x2) && is_equal(x3, x4) )
	{
		return 0;
	}
	else if( is_equal(x1, x2) )
	{		
		radians = PI/2 - abs(Angleint1(x4, y4, x3, y3, 0));
		
		double k2 = (y3 - y4) / (x3 - x4);
		if( k2 < 0 )
			radians = - radians;
	}
	else if( is_equal(x3, x4) )
	{
		radians = PI/2 - abs(Angleint1(x1, y1, x2, y2, 0) );	  
		
		double k1 = (y1 - y2) / (x1 - x2);
		/// Iris 4/13/2009 FIX_ANGLE_FUNC_BUG_REPORTED_BY_GREG
		//if( k1 > 0 )
		//	radians = - radians;	
		///end FIX_ANGLE_FUNC_BUG_REPORTED_BY_GREG
	}
	else
	{	
		double k1 = (y1 - y2) / (x1 - x2);
		double k2 = (y3 - y4) / (x3 - x4);
		
		if( k1 * k2 == -1 ) //right angle 
		{
			radians = PI / 2;
			if( k1 > k2 )
				radians = - radians;
		}
		else
		{
			radians = atan( (k2 - k1) / (1 + k1 * k2) );		
		}
	}
	
	if(1 == unit)
		return Degrees(radians);
	return radians;	
}

/// Iris 6/09/2010 ORG-260 ADD_DERIVATIVE_AND_INTEGRAL_FUNC
vector Derivative(vector vec, int order = 1) 
{
	vector vOutput;
	vOutput = vec;
	
	vector vX;
	vX.Data(1, vec.GetSize(), 1);	
	
	for(int nn = 0; nn < order; nn++)
	{
		if( OE_NOERROR != ocmath_derivative(vX, vOutput, vec.GetSize() ) )
		{
			vOutput = NANUM;
			break;
		}
	}
	return vOutput;
}

vector Integral(vector vec)
{
	vector vOutput(vec.GetSize());	
	vOutput = NANUM;
	
	vector vX;
	vX.Data(1, vec.GetSize(), 1);	
	
	IntegrationResult result;
	if( OE_NOERROR != ocmath_integrate(vX, vec, 0, vOutput.GetSize()-1, &result, vOutput) )
		vOutput = NANUM;
		
	return vOutput;	
}
///End ADD_DERIVATIVE_AND_INTEGRAL_FUNC

////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
#pragma labtalk(3, Miscellaneous)
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////

bool ISNA(double dd)
{
	return is_missing_value(dd);
}

double NA()
{
	return NANUM;
}

/*
bool BitAND(int a, int b)
{
	return a & b;
}

bool BitOR(int a, int b)
{
	return a | b;
}

bool BitXOR(int a, int b)
{
	return a ^ b;
}

*/

enum
{
	BIT_OP_AND,
	BIT_OP_OR,
	BIT_OP_XOR
};

#pragma xor(push, FALSE)

static int _bit_op(int a, int b, int option)
{
	string strHexA, strHexB;
	strHexA.Format("%p", a);
	strHexB.Format("%p", b);
	
	vector<byte> va, vb;
	BitsHex bithexConvert;
	bithexConvert.HexStrToBits(strHexA, va);
	bithexConvert.HexStrToBits(strHexB, vb);
	ASSERT(va.GetSize() == vb.GetSize());
	
	string strResult;
	for(int ii=0; ii <va.GetSize(); ii++)
	{
		byte bb;
		switch(option)
		{
		case BIT_OP_AND: //AND
			bb = va[ii] & vb[ii];
			break;
		case BIT_OP_OR: //OR
			bb = va[ii] | vb[ii];
			break;
		case BIT_OP_XOR:
			bb = va[ii] ^ vb[ii];
			break;
		default:
			ASSERT(false);
			break;
		}
		
		strResult += (string)bb;
	}
	
	char* pch = strResult.GetBuffer(ii);
	int result = Bin2Dec(pch);
	return result;		
}

int BitAND(int a, int b)
{	
	return _bit_op(a, b, BIT_OP_AND);	
}

int BitOR(int a, int b)
{
	return _bit_op(a, b, BIT_OP_OR);	
}

int BitXOR(int a, int b)
{
	return _bit_op(a, b, BIT_OP_XOR);	
}


#pragma labtalk(3, Dataset Information)

enum
{
	INDEX_NEAREST_ANY = 0,
	INDEX_NEAREST_LEFT,
	INDEX_NEAREST_RIGHT
};

int Index(double x, vector Data, int ctrl = INDEX_NEAREST_ANY)  
{
	vector 	vSrc;  
	vSrc = Data;
	int 	nRet = ocmath_is_monotonic(vSrc, vSrc.GetSize(), true);
	
	if( MONO_INCREASE == nRet || MONO_DECREASE == nRet )
	{			
		vector 	vIndex;
		vIndex.Data(1, vSrc.GetSize(), 1);
		
		double 	dDataIndex;   
		int nRet = ocmath_interpolate(&x, &dDataIndex, 1, vSrc, vIndex, vSrc.GetSize(), INTERP_TYPE_LINEAR, 1, NULL, NULL, OCMATH_SRCDATA_X_MONOTONIC);
		if(0 != nRet) 
			return -1; // error, fail to do interpolate
		
		int 	nn = -1;
		switch( ctrl )
		{
		case INDEX_NEAREST_ANY:
			nn = nint(dDataIndex);
			break;
		case INDEX_NEAREST_LEFT:
			nn = (int)dDataIndex;
			break;
		case INDEX_NEAREST_RIGHT:
			nn = (dDataIndex - (int)dDataIndex) == 0 ? dDataIndex : (int)dDataIndex + 1;
			break;
		}
		
		if( nn > vSrc.GetSize() ) // over the size of Data
			return vSrc.GetSize(); 
		
		return nn;
	}
		
	return -2; // Data is not monotonic
}



/////////////////////////////////////////////////////////////////////////////////
// End LabTalk Set Values dialog functions
/////////////////////////////////////////////////////////////////////////////////

//---- CPY 10/29/09 QA70-14563 SUPPORT_OC_FUNC_NO_NEED_TO_BE_USED_AS_LT_CMD
// function from here on are for LT usage that should not appear in SCV dialog and they can be used both as command and as functions
#pragma labtalk(1)
//----
class GetGraphPointsEx : public GetGraphPoints
{
public:
	GetGraphPointsEx()
	{
		#ifdef		_DEBUG
		out_str("GetGraphPointsEx");
		#endif		/// _DEBUG
	}
	
	virtual void OnNewPoint(int nPoint)
	{
		//out_str("OnNewPoint");
		
		GraphLayer gl = Project.ActiveLayer();
		if(gl)
		{
			DataPlot dp = gl.DataPlots(-1); // active data plot
			if(dp)
			{				
				DataPlotStrings stPlotString;
				dp.GetPlotType(&stPlotString);
				m_vsPlots.Add(stPlotString.szDepRange);
				
				#ifdef		_DEBUG
				printf("%s - %d(%s)\n", gl.GetPage().GetName(), dp.GetIndex(), stPlotString.szDepRange);
				#endif		/// _DEBUG

				///------ Folger 03/31/10 QA81-15258 XF_GETPTS_NEEDED
				string	str;
				int		nSize = m_vsPlots.GetSize();
				ocu_load_msg_str(XFINFO_PICK_POINTS_FROM_GRAPH, &str, NULL, &nSize);
				SetStatusBarText(str);
				///------ End XF_GETPTS_NEEDED
			}
		}
	}	
	
	int GetPlotStrings(vector<string>& vsPlots)
	{
		vsPlots = m_vsPlots;
		return vsPlots.GetSize();
	}
	
private:
	vector<string> m_vsPlots;
};

/**$ 
	Parameters:
		strPage = [input] graph page to get points
		nPoints = [modified] number of points to be read
		nReaderType = [input] Data Reader(0) or Screen Reader(1)
		vsPlots =  [output] the name of the plots
		vx = [output] x coordinate value of the points
		vy = [output] y coordinate value of the points
		vIndices = [output] index of the points

	Return: 
		Returns the number of selcted points if successfully, else return
		-1, if nPoints <= 0
		-2, not graph page specified by strPage in current Project
		-3, not find plot specified by strPlot in graph
		
*/
int GetPoints(string strPage, int& nPoints, int nReaderType, vector<string>& vsPlots, vector& vx, vector& vy, vector& vIndices)
{
	if( nPoints <=0 )
		return -1;
	
	GraphPage gp(strPage);
	if( !gp )
		return -2;	
		
	GraphLayer gl = gp.Layers(-1); // active layer
	
	GetGraphPointsEx sp;  
	bool bSpanToData = 0 == nReaderType? true : false;
	sp.SetFollowData(bSpanToData); 
	sp.GetPoints(nPoints, gl);
	
	vector vxTemp, vyTemp;
	vector<int> vnIndicesTemp;
	int nRet = sp.GetData(vxTemp, vyTemp, vnIndicesTemp);
	
	///------ Folger 03/31/10 QA81-15258 XF_GETPTS_NEEDED
	if ( vx )
	///------ End XF_GETPTS_NEEDED
		vx = vxTemp;
	///------ Folger 03/31/10 QA81-15258 XF_GETPTS_NEEDED
	if ( vy )
	///------ End XF_GETPTS_NEEDED
		vy = vyTemp;

	///------ Folger 03/31/10 QA81-15258 XF_GETPTS_NEEDED
	if ( vIndices )
	///------ End XF_GETPTS_NEEDED
	{
		if( bSpanToData )
		{
			vIndices = vnIndicesTemp;	
		}
		else
		{
			vIndices.SetSize(vnIndicesTemp.GetSize());
			vIndices = NANUM;
		}
	}
	///------ Folger 03/31/10 QA81-15258 XF_GETPTS_NEEDED
	if ( vsPlots )
	///------ End XF_GETPTS_NEEDED
	{
		sp.GetPlotStrings(vsPlots);
		///------ Folger 03/31/10 QA81-15258 XF_GETPTS_NEEDED
		//ASSERT(vsPlots.GetSize() == vx.GetSize());
		ASSERT(vsPlots.GetSize() == vxTemp.GetSize());
		///------ End XF_GETPTS_NEEDED
	}
	
	gp.Refresh(); // refresh graph to remove redundant cursors
	
	///------ Folger 03/31/10 QA81-15258 XF_GETPTS_NEEDED
	//nPoints = vx.GetSize();
	nPoints = vxTemp.GetSize();
	///------ End XF_GETPTS_NEEDED
	return nPoints;
}


/// Fisher qa80-14828 12/10/09 PARAM_INITIAL_VAL_TOO_BAD_FOR_SOME_DATA
void get_exponent_ex(vector &x_data, vector &y_data, double *pBase, double *pScaleExp, double *pScaleX)
{
	vector vy;
	double x0, y0;
	y0 = min( y_data );
	vy = ln( y_data - y0 );
	vector coeff;
	coeff.SetSize(2);
	fitpoly( x_data, vy, 1, coeff);  
	
	if(pBase)
		*pBase = y0;
	if(pScaleExp)
		*pScaleExp = exp(coeff[0]);
	if(pScaleX)
		*pScaleX = coeff[1];
}
/// END PARAM_INITIAL_VAL_TOO_BAD_FOR_SOME_DATA

///---Sim 04-08-2010 QA85-15284 SUPPORT_LT_FITTING_WITH_INTEGRAL
struct LTIntegralFunction
{
	//string strLTFnName;
	LPCSTR lpcszLTFnName;
	//vector<double> vParams;
	double *pParams;
	int nNumParams;
};
static double NAG_CALL lt_integral_callback(double x, Nag_User *comm)  // x is the independent variable of the integrand
{
	/*
	struct LTIntegralFunction *sp = (struct LTIntegralFunction *)(comm->p);
	
	string strParams;
	if ( sp->pParams )
	{
		double dVal;
		for ( int ii = 0; ii < sp->nNumParams; ii++ )
		{
			dVal = sp->pParams[ii];
			strParams += ",";
			strParams += dVal;
		}
	}
	
	string strLT;
	strLT += "double dIntegFuncVal;";
	strLT += "dIntegFuncVal=" + (string)(sp->lpcszLTFnName) + "(" + x + strParams + ");";
	if ( !LT_execute(strLT) )
	{
		error_report("failed execute labtalk function " + (string)(sp->lpcszLTFnName) + '()!');
		return NANUM;
	}
	
	double v = NANUM;
	LT_get_var("dIntegFuncVal", &v);
	*/
	struct LTIntegralFunction *sp = (struct LTIntegralFunction *)(comm->p);
	
	double v = NANUM;
	okoc_lt_integral_func(sp->lpcszLTFnName, &v, x, sp->pParams, sp->nNumParams);	
	
	return v;
}
 
double integralFF(string strLTFnName, double dX1, double dX2, vector<double> vParams)
{
	double epsabs = 0.00001, epsrel = 0.0000001, result = NANUM, abserr;
	Integer max_num_subint = 500;  
    // you may use epsabs and epsrel and this quantity to enhance your desired precision 
    // when not enough precision encountered
 
	Nag_QuadProgress qp;
	static NagError fail;
 
	// the parameters parameterize the integrand can be input to the call_back function
    // through the Nag_User communication struct 
    Nag_User comm;	
	struct LTIntegralFunction s;
	s.lpcszLTFnName = strLTFnName;
	//s.vParams = vParams;
	s.nNumParams = vParams.GetSize();
	s.pParams = (double*)malloc(s.nNumParams * sizeof(double));
	if ( !s.pParams )
	{
		error_report("failed alloc memory in integral function!");
		return NANUM;
	}
	for ( int nParam = 0; nParam < s.nNumParams; nParam++ )
		s.pParams[nParam] = vParams[nParam];
	
    comm.p = (Pointer)&s;
 
    okoc_lt_integral_func(strLTFnName, NULL, NANUM, NULL, 0, LT_INTEG_FUNC_INIT_CACHE);
    
    if ( !is_missing_value(dX1) && !is_missing_value(dX2))
    {
		d01sjc(lt_integral_callback, dX1, dX2, epsabs, epsrel, max_num_subint, &result, &abserr, &qp, &comm, &fail);
    }
    else
    {
    	Nag_BoundInterval bound = Nag_LowerSemiInfinite;
    	double x = dX2;
    	if ( is_missing_value(dX1) && is_missing_value(dX2) )
    	{
    		bound = Nag_Infinite;
    		x = NANUM;
    	}
    	if ( is_missing_value(dX1) )
    	{
    		bound = Nag_LowerSemiInfinite;
    		x = dX2;
    	}
    	if ( is_missing_value(dX2) )
    	{
    		bound = Nag_UpperSemiInfinite;
    		x = dX1;
    	}
    		
		d01smc(lt_integral_callback, bound, x, epsabs, epsrel, max_num_subint, &result, &abserr, &qp, &comm, &fail);
    } 
 
    okoc_lt_integral_func(strLTFnName, NULL, NANUM, NULL, 0, LT_INTEG_FUNC_DISPOSE_CACHE);
    // you may want to exam the error by printing out error message, just uncomment the following lines
	// if (fail.code != NE_NOERROR)
        // printf("%s\n", fail.message);
 
 
	// For the error other than the following three errors which are due to bad input parameters 
	// or allocation failure  NE_INT_ARG_LT  NE_BAD_PARAM   NE_ALLOC_FAIL
	// You will need to free the memory allocation before calling the integration routine again to 
        // avoid memory leakage
	if (fail.code != NE_INT_ARG_LT && fail.code != NE_BAD_PARAM && fail.code != NE_ALLOC_FAIL)
	{
		NAG_FREE(qp.sub_int_beg_pts);
		NAG_FREE(qp.sub_int_end_pts);
		NAG_FREE(qp.sub_int_result);
		NAG_FREE(qp.sub_int_error);
	}
 
	if ( s.pParams )
		free(s.pParams);
	
	return result;
}
///---END QA85-15284 SUPPORT_LT_FITTING_WITH_INTEGRAL
