// $Id: WCMath.h,v 1.1 2011/01/31 16:06:52 samn Exp $ 
#pragma once
#include <cmath>
#include <vector>
#include "A2D.h"

using namespace std;

inline long Factorial(long n)
{
	if(n <= 1) return 1;
	return n * Factorial(n - 1);
}

inline long DFactorial(long n)
{
	if(n<=1) return 1;
	return n * DFactorial(n-2);
}

//NR version of logarithm of gamma function
inline float lngamma(float xx)
{
	double x,y,tmp,ser;
	static double cof[6]={76.18009172947146,-86.50532032941677,
		24.01409824083091,-1.231739572450155,
		0.1208650973866179e-2,-0.5395239384953e-5};
	int j;

	y=x=xx;
	tmp=x+5.5;
	tmp -= (x+0.5)*log(tmp);
	ser=1.000000000190015;
	for (j=0;j<=5;j++) ser += cof[j]/++y;
	return -tmp+log(2.5066282746310005*ser/x);
}

inline double Gamma(double R)
{
	int iR = (int) R;
	double d = iR;
	if(d == R)
	{
		return Factorial(iR-1);
	}
	else
	{
		const double PI=3.14159265358979323846;
		//R = n/2 + 1;
		//n = (R - 1)*2;		
		int n = static_cast<int>(( R - 1.0 ) * 2.0);
		return (sqrt(PI) * DFactorial(n)) / pow(2,(n+1)/2);
	}
}

//monotonically decreasing
template< class T >
inline bool MonoDec(std::vector<T>& v)
{
	int isz = v.size(), i, j;
	for(i=0;i<isz-1;i++)
	{
		for(j=i+1;j<isz;j++)
		{
			if(v[i]<v[j])
				return false;
		}
	}
	return true;
}

template< class T > 
T Sum(std::vector<T>& v)
{
	int i;
	T val = T(0);
	for(i=0;i<v.size();i++) val += v[i];
	return val;
}

template< class T > 
T Sum(std::vector<T>& v,int iStart,int iEnd)
{
	int i;
	T val = T(0);
	for(i=iStart;i<v.size() && i<iEnd;i++) val += v[i];
	return val;
}

template< class T > 
T Avg(std::vector<T>& v)
{
	if(!v.size()) return T(0);
	int i;
	T val = T(0);
	for(i=0;i<v.size();i++) val += v[i];
	return val / (T) v.size();
}

template< class T > 
T Avg(std::vector<T>& v,int iStart,int iEnd)
{
	if(!v.size()) return T(0);
	int i;
	T val = T(0);
	for(i=iStart;i<iEnd;i++) val += v[i];
	return val / (T) (iEnd-iStart);
}

template< class T >
T Variance(std::vector<T>& v,int iStart,int iEnd)
{	T avg = Avg(v,iStart,iEnd);
	int i;
	T var = T(0) , val = T(0);
	for(i=iStart;i<iEnd;i++)
	{	val = v[i]-avg;
		var += val*val;
	}
	var /= (T)(iEnd-iStart);
	return var;
}

//returns sum(1..N)
inline int IntegerSum(int N)
{
	return (N*N+N)/2;
}

using namespace std;

//#define SAMPLE_MEAN

// Variance-covariance matrix: creation
template<class REAL>
void CovarMat(const vector< vector<REAL> >& vDataIn,int iRows,int iCols,vector< vector<REAL> >& vCovarMat, vector<REAL>& vMean,bool bCorrelMat = false)
{	// Create iCols * iCols covariance matrix from given iRow * iCol data matrix. 
	int i, j, x2, x, y, j1 , j2;
	// Allocate storage for mean vector 
	vMean = vector<REAL>(iCols,0.0);

	vector< vector<REAL> > vData(vDataIn); // copy it!

#ifdef SAMPLE_MEAN
	REAL N = iRows - 1;
#else
	REAL N = iRows;
#endif

	// Determine mean of row vectors of input data matrix 
	for(x=0;x<iCols;x++)
	{	for(y=0;y<iRows;y++)
			vMean[x] += vData[y][x];
		vMean[x] /= N;
	}

	vector<REAL> vStdev(iCols);
	
	// Center the row vectors.
	for(y=0;y<iRows;y++)
    {	for(x=0;x<iCols;x++)
		{	vData[y][x] -= vMean[x];
		}
    }

	// compute std-dev
	//const bool bCorrelMat = false;
	if(bCorrelMat) for(x=0;x<iCols;x++)
	{	for(y=0;y<iRows;y++)
			vStdev[x]+=vData[y][x]*vData[y][x];
		vStdev[x] /= N;
		vStdev[x] = sqrt(vStdev[x]);
	}

	vCovarMat = vector<vector<REAL> >(iCols,vector<REAL>(iCols,0.0));

	// Calculate the iCols * iCols covariance matrix.
	for(j1=0;j1<iCols;j1++)
	{	for(j2=0;j2<=j1;j2++)
		{	vCovarMat[j1][j2]=0.0;
			for(i=0;i<iRows;i++)
				vCovarMat[j1][j2] += ( vData[i][j1] * vData[i][j2] );
			if(bCorrelMat)	//correlation matrix
				vCovarMat[j1][j2] /= (N * vStdev[j1] * vStdev[j2] );
			else	//covariance matrix
				vCovarMat[j1][j2] /= N;
			vCovarMat[j2][j1]=vCovarMat[j1][j2];
		}
	}
}

// Variance-covariance matrix: creation
template<class REALD , class REALC, class REALM>
void CovarMat(  A2D<REALD>& vData , int iRows, int iCols, A2D<REALC>& vCovarMat, vector<REALM>& vMean,bool bCorrelMat = false)
{	// Create iCols * iCols covariance matrix from given iRow * iCol data matrix. 
	int i, x, y, j1 , j2;
	// Allocate storage for mean vector 
	vMean = vector<REALM>(iCols);

#ifdef SAMPLE_MEAN
	REALM N = iRows - 1;
#else
	REALM N = iRows;
#endif

	// Determine mean of row vectors of input data matrix 
	for(x=0;x<iCols;x++)
	{	vMean[x] = 0.0;
		for(y=0;y<iRows;y++)
		{	vMean[x] += vData[y][x];
		}
		vMean[x] /= N;
	}

	vector<REALM> vStdev(iCols);
	
	// Center the row vectors.
	for(y=0;y<iRows;y++)
    {	for(x=0;x<iCols;x++)
		{	vData[y][x] -= vMean[x];
		}
    }

	// compute std-dev
	//const bool bCorrelMat = false;
	if(bCorrelMat) for(x=0;x<iCols;x++)
	{	for(y=0;y<iRows;y++)
			vStdev[x]+=vData[y][x]*vData[y][x];
		vStdev[x] /= N;
		vStdev[x] = sqrt(vStdev[x]);
	}

	vCovarMat.Init(iCols,iCols);

	if(bCorrelMat)	// Calculate the iCols * iCols correlation matrix.
	{	for(j1=0;j1<iCols;j1++)
		{	for(j2=0;j2<=j1;j2++)
			{	vCovarMat[j1][j2]=0.0;
				for(i=0;i<iRows;i++)
				{	vCovarMat[j1][j2] += (vData[i][j1] * vData[i][j2]);
				}
				vCovarMat[j1][j2] /= ((REALC) N * vStdev[j1] * vStdev[j2] );
				vCovarMat[j2][j1]=vCovarMat[j1][j2];
			}
		}
	}
	else	// Calculate the iCols * iCols covariance matrix.
	{	for(j1=0;j1<iCols;j1++)
		{	for(j2=0;j2<=j1;j2++)
			{	vCovarMat[j1][j2]=0.0;
				for(i=0;i<iRows;i++)
				{	vCovarMat[j1][j2] += (vData[i][j1] * vData[i][j2]);
				}
				vCovarMat[j1][j2] /= N;
				vCovarMat[j2][j1]=vCovarMat[j1][j2];
			}
		}
	}

	// Un-Center the row vectors.
	for(y=0;y<iRows;y++)
    {	for(x=0;x<iCols;x++)
		{	vData[y][x] += vMean[x];
		}
    }
}

// Variance-covariance matrix: creation
template<class REAL>
void CovarMat(const vector< REAL >& vDataIn,int iRows,int iCols,vector< vector<REAL> >& vCovarMat, vector<REAL>& vMean,bool bCorrelMat = false)
{	// Create iCols * iCols covariance matrix from given iRow * iCol data matrix. 
	int i, j, x2, x, y, j1 , j2;
	// Allocate storage for mean vector 
	vMean = vector<REAL>(iCols);

	vector< REAL > vData(vDataIn); // copy it!

#ifdef SAMPLE_MEAN
	REAL N = iRows - 1;
#else
	REAL N = iRows;
#endif

	// Determine mean of row vectors of input data matrix 
	for(x=0;x<iCols;x++)
	{	vMean[x] = 0.0;
		for(y=0;y<iRows;y++)
		{	vMean[x] += vData[y*iCols+x];
		}
		vMean[x] /= N;
	}

	vector<REAL> vStdev(iCols);
	
	// Center the row vectors.
	for(y=0;y<iRows;y++)
    {	for(x=0;x<iCols;x++)
		{	vData[y*iCols+x] -= vMean[x];
		}
    }

	// compute std-dev
	//const bool bCorrelMat = false;
	if(bCorrelMat) for(x=0;x<iCols;x++)
	{	for(y=0;y<iRows;y++)
			vStdev[x]+=vData[y*iCols+x]*vData[y*iCols+x];
		vStdev[x] /= N;
		vStdev[x] = sqrt(vStdev[x]);
	}

	vCovarMat = vector<vector<REAL> >(iCols,vector<REAL>(iCols,0.0));

	// Calculate the iCols * iCols covariance matrix.
	for(j1=0;j1<iCols;j1++)
	{	for(j2=0;j2<=j1;j2++)
		{	vCovarMat[j1][j2]=0.0;
			for(i=0;i<iRows;i++)
			{	vCovarMat[j1][j2] += (vData[i*iCols+j1] * vData[i*iCols+j2]);
			}
			if(bCorrelMat)	//correlation matrix
				vCovarMat[j1][j2] /= ( N * vStdev[j1] * vStdev[j2] );
			else	//covariance matrix
				vCovarMat[j1][j2] /= N;
			vCovarMat[j2][j1]=vCovarMat[j1][j2];
		}
	}
}

template< class T >
int GetMaxIndex(vector<T>& v)
{
	int iSz = v.size(), i = 0, iMaxIDX = 0;
	if(iSz<1) return -1;
	T maxVal = v[0];
	for(i=1;i<iSz;i++)
	{
		if(v[i]>maxVal)
		{
			maxVal=v[i];
			iMaxIDX=i;
		}
	}
	return iMaxIDX;
}

template< class T >
int GetMinIndex(vector<T>& v)
{
	int iSz = v.size(), i = 0, iMinIDX = 0;
	if(iSz<1) return -1;
	T minVal = v[0];
	for(i=1;i<iSz;i++)
	{
		if(v[i]<minVal)
		{
			minVal=v[i];
			iMinIDX=i;
		}
	}
	return iMinIDX;
}

template< class T >
T Energy(vector<T>& v)
{
	int iSz = v.size() , i = 0;
	T sum(0);
	if(iSz<1) return sum;
	for(;i<iSz;i++) sum += v[i]*v[i];
	sum = sqrt(sum);
	return sum / (double) iSz;
}

//chi-squared cumulative distribution function
//returns probability of chi-squared with df degrees of freedom has value <= x
double chisqCDF(double x, double df);

inline bool PowOf2(int N)
{
	if(N < 0) return false;
	while(N>1)
	{
		if(N%2==1) return false;
		N /= 2;
	}
	return true;
}

inline int CeilPowOf2(int N)
{
	int T = 1;
	while(T < N) T *= 2;
	return T;
}