/*------------------------------------------------------------------------------*
 * 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			
 *------------------------------------------------------------------------------*/
 
#include <origin.h> // main Origin C header that is precompiled and already include most headers 

// add your additional include files here
// ***************************
// functions that accept scalars and return a scalar
// ***************************

//------- 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)
{
	return (da > db ? da : db);
}

// returns min of two doubles
double min(double da, double db)
{
	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
// ***************************

// 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;
}

// returns 25th percentile value of dataset
double y25(Dataset &aa)
{
	double perc, result;
	BOOL interp = false;
	perc=25;
	Data_percentiles(&aa, &perc, &result, 1,interp);
	return result;
}

// returns 50th percentile value of dataset
double y50(Dataset &aa)
{
	double perc, result;
	BOOL interp = false;
	perc=50;
	Data_percentiles(&aa, &perc, &result, 1,interp);
	return result;
}
		
// returns 75th percentile value of dataset
double y75(Dataset &aa)
{
	double perc, result;
	BOOL interp = false;
	perc=75;
	Data_percentiles(&aa, &perc, &result, 1,interp);
	return result;
}


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

// 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));
}		

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);	
}

//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;
}


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


// 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);	
		
	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));
		if(fabs(dy1)<fabs(dy2)) 
			return xatymax(cc);
		else
			return xatymin(cc);
	}      
}
//************** 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;
//	
//}