/*------------------------------------------------------------------------------*
 * File Name:math_utils.c 														*
 * Creation: SDB 5/14/03														*
 * Purpose: Origin's basic math related routines								*
 * Copyright (c) Originlab Corp.	2003										*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *	Leo 08/29/05 ALL_FUNCTIONS_HAVE_LLVC_COUNTERPARTS							*
 *------------------------------------------------------------------------------*/
 
///Leo 08/29/05 ALL_FUNCTIONS_HAVE_LLVC_COUNTERPARTS, Thus commented out

//////////////////////////////////////////////////////////////////////////////////////
//// you can include just this typical header file for most Origin built-in functions and classes
//// and it takes a reasonable amount of time to compile, 
//#include <origin.h> // The Project and Application class 
//#include <tree.h>
//#include <tree_utils.h>
//#include <analysis_utils.h>
//#include <matrix.h>
//#include <NAG\OCN_e02.h> // this contains all the NAG headers, 
//
//
//////////////////////////////////////////////////////////////////////////////////////
//#define _DBINT(_STR, _INT)	out_int(_STR, _INT);
//#define _DBMGS(_STR)		out_str(_STR);
//#define _DB_LIST(_STR)		nlsf_list_params(_STR);	
//	
///*
//	This function perform Adjacent Averaging Smoothing on the curve.
//	The parameter nOption specifies the methods to handle the ends of the data.
//Parameter:
//	crvSignal: Input&Output, the original data need to smooth and the result after smoothing
//	nPts: Input, the degree of the smoothing which specified the number of points of 
//	          the smoothing window
//	nOption: Input, the methods to handle the ends of the data
//		0: Missing Values Method. 
//	  	1: Reflected Method. Uses w[i] in place of the missing w[-i] and w[n-i-1] in place 
//	  	     of the missing w[n+i-1]. 
//		2: Wrap Method. Uses w[n-i] in place of the missing w[-i] and vice-versa.
//		3: End Weighted Method. Uses w[0] in place of the missing w[-i] and w[n-1] in place of 
//		     the missing w[n+i-1].
//Return:
//	0: Successfully called. 
//	-1: Invalid method number.
//	-2: Invalid curve or empty curve.
//*/
//int smooth_adjave(Curve& crvSignal, int nPts = 5, int nOption = 0)
//{
//	// Check the input parameters
//	if(!crvSignal.IsValid()||crvSignal.GetSize()<1) return -2;
//	
//	if(nPts<=1) return 0;
//	
//	//the number of points on the left of the target one which are used in the calculation
//	int nleft=nPts/2;
//	
//	// Cut down the invalid size of the smoothing window for method 0, 1, 2.
//	int nSize=crvSignal.GetSize();
//	if(nOption<3)
//	{
//		while(nleft>=nSize)
//			nleft-=1;
//	}
//	
//	// Smoothing
//	switch(nOption)
//	{
//	case 0: smooth_missvalue(crvSignal, nleft);
//			break;
//	case 1: smooth_reflected(crvSignal, nleft);
//			break;
//	case 2: smooth_wrap(crvSignal, nleft);
//			break;
//	case 3: smooth_endweighted(crvSignal, nleft);
//			break;
//	default: return -1;
//	}
//	
//	return 0;
//}
//
//// Adjacent Averaging Smoothing with missing value end-points-method.
//void smooth_missvalue(vectorbase& vbSignal, int nleft)
//{
//	int nSize=vbSignal.GetSize();
//	vector vectemp(nSize);
//	
//	double sum=0.;
//	for(int i=0; i<=nleft; i++)
//	{
//		sum+=vbSignal[i];
//	}
//	
//	int k=nleft+1;
//	vectemp[0]=sum/k;
//	
//	for(i=1; i<nSize; i++)
//	{
//		if(i-nleft>0)
//		{
//			sum-=vbSignal[i-nleft-1];
//			k--;
//		}
//		if(i+nleft<nSize)
//		{
//			sum+=vbSignal[i+nleft];
//			k++;
//		}
//		vectemp[i]=sum/k;
//	}
//	
//	vbSignal=vectemp;
//}
//
////Adjacent Averaging Smoothing with reflected end-points-method.
//void smooth_reflected(vectorbase& vbSignal, int nleft)
//{
//	int nSize=vbSignal.GetSize();   
//	vector vectemp(nSize);
//	
//	double sum=0.0;
//	for(int i=1; i<=nleft; i++)
//	{
//		sum+=2*vbSignal[i];
//	}
//	
//	int k=2*nleft+1;
//	sum+=vbSignal[0];
//	vectemp[0]=sum/k;
//	
//	int nSize1=nSize-1;
//	int nleft1=nleft+1;
//	
//	for(i=1; i<nSize; i++)
//	{
//		sum-=vbSignal[abs(i-nleft1)];
//		sum+=vbSignal[nSize1-abs(i+nleft-nSize1)];
//		vectemp[i]=sum/k;
//	}
//	
//	vbSignal=vectemp;
//}	
//
//// Adjacent Averaging Smoothing with wrap end-points-method.
//int smooth_wrap(vectorbase& vbSignal, int nleft)
//{
//	int nSize=vbSignal.GetSize();   
//	vector vectemp(nSize);
//
//	double sum=0.0;
//	for(int i=1; i<=nleft; i++)
//	{
//		sum+=vbSignal[i]+vbSignal[nSize-i];
//	}
//	
//	int k=2*nleft+1;
//	sum+=vbSignal[0];
//	vectemp[0]=sum/k;	
//	
//	int nSize1=nSize-1;
//	
//	for(i=1; i<nSize; i++)
//	{
//		sum-=vbSignal[mod(i-nleft+nSize1, nSize)];
//		sum+=vbSignal[mod(i+nleft, nSize)];
//		vectemp[i]=sum/k;
//	}
//	
//	vbSignal=vectemp;
//	return nleft;
//}
//	
//// Adjacent Averaging Smoothing with end weighted method.
//void smooth_endweighted(vectorbase& vbSignal, int nleft)
//{
//	int nSize=vbSignal.GetSize();  
//	vector vectemp(nSize);
//	
//	vectemp=vbSignal;
//	vectemp.InsertAt(nSize, vbSignal[nSize-1], nleft);
//	vectemp.InsertAt(0, vbSignal[0], nleft);
//
//	double sum=0.0;
//	for(int i=0; i<=2*nleft; i++)
//	{
//		sum+=vectemp[i];
//	}
//	
//	int k=2*nleft+1;	
//	vbSignal[0]=sum/k;
//	
//	for(i=1; i<nSize; i++)
//	{
//		sum=sum-vectemp[i-1]+vectemp[i+k-1];
//		vbSignal[i]=sum/k;
//	}
//}
//
///*
//	This function smooth the curve using the Savitzky-Golay filter method.
//Parameter:
//	crvSignal: Input&Output, the curve need to smooth and the smoothing results.
//	nPolyOrder: Input, the degree of the polynomial filters.
//	nPtsLeft: Input, the number of left points to the target one which are included in
//	        the calculation.
//	nPtsRight: Input, the number of points on the right of the target one.
//Return:
//	0: Successfully called
//	-1: The degree nPolyOrder must less than nPtsLeft+nPtsRight.
//	-2: nPtsLeft+nPtsRight must less than the size of the curve.
//	-3: Invalid curve or empty curve.
//	-4: Positive needed for nPolyOrder and nPtsLeft and nPtsRight.
//	positive: NAG error codes 
//*/
//int smooth_sg(Curve& crvSignal, int nPolyOrder = 2, int nPtsLeft = 2, int nPtsRight = 2)
//{
//	// Check the input parameter
//	int nSize=crvSignal.GetSize();
//	if(!crvSignal.IsValid()||nSize<1)
//		return -3;
//	crvSignal.Sort();
//	
//	int nPs=nPtsLeft+nPtsRight+1;
//	if(nPs>nSize) return -2;
//	if(nPolyOrder>=nPs-1)
//		return -1;
//	
//	if(nPolyOrder<0||nPtsLeft<0||nPtsRight<0)
//		return -4;
//	
//	// Define variable for nag function
//	int kplus=nPolyOrder+1;             //order+1
//	
//	vector vectempx(nPs),vectempy(nPs);
//	vector vecW(nPs);      				//weights
//	vecW=1;
//	
//	matrix mcoeff(kplus,kplus);			//the coefficients of polynomial
//	vector vecCoeff(kplus);
//	
//	vector vecS(kplus);					//the root mean square residual 
//	
//	vector vecResult(nSize);            //the smoothing results
//	
//	// Obtain the associated x column of the curve
//	Dataset dsX;
//	if(crvSignal.HasX())
//	{
//		crvSignal.AttachX(dsX);
//	}
//	else
//	{
//		dsX.SetSize(nSize);
//		dsX.Data(1, nSize);
//	}
//	
//	int nRet;
//	int iMax, iMin;  // the two ends of the smooth window
//	
//	// Deal with the points in the begin
//	iMin=0;
//	iMax=nPs-1;
//	
//	for(int i=0;i<=iMax; i++)
//	{
//		vectempy[i]=crvSignal[i];
//		vectempx[i]=dsX[i];
//	}
//	
//	// 1D least-squares polynomial fitting
//	nRet=nag_1d_cheb_fit(nPs,kplus,kplus,vectempx,vectempy,vecW,mcoeff,vecS);  
//	if(nRet!=0)
//		return nRet;
//	
//	for(i=0; i<kplus; i++)
//		vecCoeff[i]=mcoeff[nPolyOrder][i];
//	
//	for(i=0;i<=nPtsLeft;i++)
//	{
//		double xcap=(2*dsX[i]-dsX[iMin]-dsX[iMax]) / (dsX[iMax]-dsX[iMin]);
//		nag_1d_cheb_eval(kplus,vecCoeff,xcap,&vecResult[i]);
//	}
//	
//	// Deal with the points in the middle
//	for(i=nPtsLeft+1;i<nSize-nPtsRight;i++)
//	{
//		iMax=i+nPtsRight;
//		iMin=i-nPtsLeft;
//		
//		vectempy.RemoveAt(0);
//		vectempy.Add(crvSignal[iMax]);
//		vectempx.RemoveAt(0);
//		vectempx.Add(dsX[iMax]);
//		
//		nRet=nag_1d_cheb_fit(nPs,kplus,kplus,vectempx,vectempy,vecW,mcoeff,vecS);
//		if(nRet!=0)
//			return nRet;
//			
//		for(int j=0; j<kplus; j++)
//			vecCoeff[j]=mcoeff[nPolyOrder][j];
//
//		double xcap=(2*dsX[i]-dsX[iMin]-dsX[iMax]) / (dsX[iMax]-dsX[iMin]);
//		nag_1d_cheb_eval(kplus,vecCoeff,xcap,&vecResult[i]);
//	}
//
//	// Deal with the end points 
//	for(i=nSize-nPtsRight;i<nSize;i++)
//	{
//		double xcap=(2*dsX[i]-dsX[iMin]-dsX[iMax]) / (dsX[iMax]-dsX[iMin]);
//		nag_1d_cheb_eval(kplus,vecCoeff,xcap,&vecResult[i]);
//	}
//	
//	crvSignal=vecResult;
//	return 0;
//}
//
//////////////////////////////////////////////////////////////////////////////////////
//// start your functions here
///*
//	This function smooth the curve by FFT filter
//Parameter:
//	crvSignal: Input&Output, the curve need to smooth and the results after 
//smoothing
//	nPts: Input, the number of data points considered at a time
//Return:
//	0: Successfully called
//	-1: Invalid curve.
//	-2: nPts should be positive
//	-3: Error occurs when perform FFT
//	-4: Error occurs when perform IFFT
//*/
//int smooth_fft(Curve& crvSignal, int nPts = 5)
//{
//	if(!crvSignal.IsValid())
//		return -1;
//
//	if(nPts<0) return -2;
//	else if(nPts==0) return 0;
//
//	int nSize=crvSignal.GetSize();
//
//	vector<complex> vecFFT(nSize);
//	vecFFT=crvSignal;
//
//	int nRet=FFT(vecFFT, vecFFT);
//	if(nRet!=0) return -3;
//
//	// Remove the high frequency components
//	int nFc1=nSize/nPts;
//	int nFc=nFc1+1;
//
//	double dWindow;
//	for(int i=1;i<nFc;i++)
//	{
//		dWindow=-1.0*i^2/nFc1^2+1;
//		vecFFT[i]*=dWindow;
//		vecFFT[nSize-i]*=dWindow;
//	}
//	for(i=nFc; i<=nSize/2; i++)
//	{
//		vecFFT[i]=0;
//		vecFFT[nSize-i]=0;
//	}
//
//	// Compute inverse FFT, extract real part, and place in signal curve
//	nRet = IFFT(vecFFT, vecFFT);
//	if( nRet !=0 )
//		return -4;
//
//	vector vecReal(nSize);
//	vecFFT.GetReal(vecReal);
//	crvSignal = vecReal;
//
//	return 0;
//}
//

///End ALL_FUNCTIONS_HAVE_LLVC_COUNTERPARTS