#ifndef LAPLACIAN_MATRIX_1D_INCLUDED
#define LAPLACIAN_MATRIX_1D_INCLUDED
#include "LinearAlgebra/SparseMatrix.h"
#include "FunctionBasis/PPolynomial.h"

#define MISHA_CODE 0

#define NEW_LAPLACIAN_CODE 1
enum
{
	ZERO_VALUE,
	ZERO_DERIVATIVE,
	PERIODIC_BOUNDARY,
	TOTAL_BOUNDARY_TYPES
};
#define DERIVATIVE(Type) ( ((~(Type&1))&1) | (Type&2) )

static const char* BoundaryTypes[]=
{
	"Zero Boundary Values",
	"Zero Boundary Derivatives",
	"Periodic Boundary"
};

//////////////////////////
// BinomialCoefficients //
//////////////////////////
template<int Degree>
class BinomialCoefficients
{
public:
	int coeffs[Degree+1];
	BinomialCoefficients(void);
};

/////////////////////////
// FiniteDifferences1D //
/////////////////////////
template<class Real,int Type>
class FiniteDifferences1D
{
	static void SetValues(int dim,int i,MatrixEntry<Real>* values,bool makePositive);
public:
	static void GetMatrix		(const int& dim,SparseMatrix<Real>& m,bool makePositive);
	static void StripDiagonal	(const SparseMatrix<Real>& M,Vector<MatrixEntry<Real> >& D,bool clearDiagonal);
};


//////////////////////
// FiniteElements1D //
//////////////////////
template< class Real , int Type , int Degree >
class FiniteElements1D
{
	template<class Data>
	static bool UpSample(Vector<Data>& in,Vector<Data>& out,const Real* weights,int weightCount);
	template<class Data>
	static bool DownSample(Vector<Data>& in,Vector<Data>& out,const Real* weights,int weightCount);
	static bool UpSampleMatrix(SparseMatrix<Real>& M,int lowD,const Real* weights,int weightCount,int& highD);
	static bool DownSampleMatrix(SparseMatrix<Real>& M,int highD,const Real* weights,int weightCount,int& lowD);
public:
	static int DimensionOffset(void);
	static int DomainSize(int dim);
	static int Dimension(int domainSize);

	static int DerivativeMatrix(int dim,SparseMatrix<Real>& m);

	template<class Data>
	static bool UpSample(Vector<Data>& in,Vector<Data>& out);
	template<class Data>
	static bool DownSample(Vector<Data>& in,Vector<Data>& out);

	static bool IsDownSamplable(int highD,int& lowD);
	static bool IsUpSamplable(int lowD,int& highD);
	static bool UpSampleMatrix(SparseMatrix<Real>& M,int lowD,int& highD);
	static bool DownSampleMatrix(SparseMatrix<Real>& M,int highD,int& lowD);

	static bool IdentityMatrix(int dim,SparseMatrix<Real>& m,bool rowMajor=true);
	static bool DotProductMatrix(int dim,SparseMatrix<Real>& m);
	static bool DotProductMatrix(int dim,SparseMatrix<Real>& m,Real min,Real max);
	static bool LaplacianMatrix(int dim,SparseMatrix<Real>& m,bool weakForm,bool negate=false);
	static bool LaplacianMatrix(int dim,SparseMatrix<Real>& m,Real min,Real max,bool weakForm,bool negate=false);
	static void StripDiagonal(const SparseMatrix<Real>& M,Vector<MatrixEntry<Real> >& D,bool clearDiagonal);


	template<int DotType=Type,int DotDegree=Degree>
	class DotProduct
	{
	public:
		class Helper
		{
			double shift,dotShift;
			static const int Start,Stop;
			PPolynomial< Degree    , double > F1;
			PPolynomial< DotDegree , double > F2;
			PPolynomial< Degree+DotDegree , double > F1F2[Degree+DotDegree+1];
			double min[Degree+DotDegree+1],max[Degree+DotDegree+1];
			double fullValues[Degree+DotDegree+1];
		public:
			Helper(void);
			void setDerivatives(int d1,int d2);
			double GetValue(int i,int j,double min,double max) const;
			void SetValues(int dim,int i,MatrixEntry<Real>* values,Real min,Real max);

			static int StartOffset(void);
			static int StopOffset(void);
		};
		static bool DerivativeMatrix(int dim,SparseMatrix<Real>& m,int d1,int d2,bool negate=false);
		static bool DerivativeMatrix(int dim,SparseMatrix<Real>& m,Real min,Real max,int d1,int d2,bool negate=false);
		class FullDotProductStencil
		{
		public:
			class DotProductStencil
			{
			public:
				Real values[Degree+DotDegree+1];
			};
			DotProductStencil caseTable[2*DotDegree+1];
		};
		class FullPaddedDotProductStencil
		{
		public:
			class PaddedDotProductStencil
			{
			public:
				Real values[Degree+DotDegree+1];
			};
			PaddedDotProductStencil caseTable[2*DotDegree+3];
		};
		static bool DotProductStencil(int dim,FullDotProductStencil& s,int d1,int d2,bool negate=false);
		static bool PaddedDotProductStencil( int dim , FullPaddedDotProductStencil& s , int d1 , int d2 , bool negate=false );
	};
	class FullProlongationStencil
	{
	public:
		class ProlongationStencil
		{
		public:
			static const int Size=Degree+2;
			static int Start(int i) { return (i<<1)-((Degree+1)>>1); }
			Real values[Size];
		};
		ProlongationStencil caseTable[2*Degree+1];
		void Normalize(int lowD);
	};
	static bool ProlongationStencil(int lowD,FullProlongationStencil &s,int& highD);

#if NEW_LAPLACIAN_CODE
	// BADNESS!!! This is only going to work for Type==ZERO_DERIVATIVE indexing
	class FullRestrictionStencil
	{
	public:
		class RestrictionStencil
		{
		public:
			static const int Size=((Degree+1)>>1)<<1;
//			static const int Size=(Degree+3)>>1;
			static int Start(int i) { return (i-Degree+1)>>1; }
			Real values[Size];
		};
		RestrictionStencil caseTable[2*Degree+2];	// Need +2 to handle even/odd parity in the interior
	};
	static bool RestrictionStencil(int highD,FullRestrictionStencil &s,int& lowD);
#endif // NEW_LAPLACIAN_CODE
};

template<class Real,int Type,int Degree>
void CombineStencils(const typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullDotProductStencil& dStencil,
					 const typename FiniteElements1D<Real,Type,Degree>::FullProlongationStencil& pStencil,int dim,
					 typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullDotProductStencil& newDStencil);
template<class Real,int Type,int Degree>
void CombineStencils(const typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullDotProductStencil& dStencil,
					 const typename FiniteElements1D<Real,Type,Degree>::FullProlongationStencil& pStencil,int lowDim,int highDim,
					 typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullDotProductStencil& newDStencil);
template<class Real,int Type,int Degree>
void CombineStencils(const typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullDotProductStencil& dStencil,
					 const typename FiniteElements1D<Real,Type,Degree-1>::FullProlongationStencil& pStencil,int dim,
					 typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullDotProductStencil& newDStencil);
template<class Real,int Type,int Degree>
void CombineStencils(const typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullDotProductStencil& dStencil,
					 const typename FiniteElements1D<Real,Type,Degree+1>::FullProlongationStencil& pStencil,int dim,
					 typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullDotProductStencil& newDStencil);

template<class Real,int Type,int Degree>
void CombineStencils(const typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullPaddedDotProductStencil& dStencil,
					 const typename FiniteElements1D<Real,Type,Degree>::FullProlongationStencil& pStencil,int dim,
					 typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullPaddedDotProductStencil& newDStencil);
template<class Real,int Type,int Degree>
void CombineStencils(const typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullPaddedDotProductStencil& dStencil,
					 const typename FiniteElements1D<Real,Type,Degree>::FullProlongationStencil& pStencil,int lowDim,int highDim,
					 typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullPaddedDotProductStencil& newDStencil);
template<class Real,int Type,int Degree>
void CombineStencils(const typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullPaddedDotProductStencil& dStencil,
					 const typename FiniteElements1D<Real,Type,Degree-1>::FullProlongationStencil& pStencil,int dim,
					 typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullPaddedDotProductStencil& newDStencil);
template<class Real,int Type,int Degree>
void CombineStencils(const typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullPaddedDotProductStencil& dStencil,
					 const typename FiniteElements1D<Real,Type,Degree+1>::FullProlongationStencil& pStencil,int dim,
					 typename FiniteElements1D<Real,Type,Degree>::template DotProduct<Type,Degree>::FullPaddedDotProductStencil& newDStencil);
#include "LaplacianMatrix/LaplacianMatrix1D.inl"
#endif // LAPLACIAN_MATRIX_INCLUDED