#ifndef __GEOWORLDVMPDG_H__
#define __GEOWORLDVMPDG_H__

#include <mars_grids.h>
#include <mars_ptr.h>
#include <mars_util.h>

#include <vector>
#include <set>

namespace geoworld
{
	template< typename T >
	class MPDGVecType
	{
	public:
		typedef T Precision;

		T altitude;
		unsigned short mindex;

		MPDGVecType ()
			: altitude(0), mindex(0) {}

		MPDGVecType (const MPDGVecType & copy)
			: altitude(copy.altitude), mindex(copy.mindex) {}
		MPDGVecType (const T altitude, const unsigned short nIdx)
			: altitude(altitude), mindex(nIdx) {}

		inline MPDGVecType & operator += (const MPDGVecType & v)
			{ return operator += (v.altitude); }
		inline MPDGVecType & operator -= (const MPDGVecType & v)
			{ return operator -= (v.altitude); }
		inline MPDGVecType operator + (const MPDGVecType & v) const
			{ return operator + (v.altitude); }
		inline MPDGVecType operator - (const MPDGVecType & v) const
			{ return operator - (v.altitude); }
		inline MPDGVecType & operator += (const T & n)
		{
			altitude += n;
			return *this;
		}
		inline MPDGVecType & operator -= (const T & n)
		{
			altitude -= n;
			return *this;
		}
		inline MPDGVecType operator + (const T & n) const
		{
			return MPDGVecType(altitude + n, mindex);
		}
		inline MPDGVecType operator - (const T & n) const
		{
			return MPDGVecType(altitude - n, mindex);
		}
		inline MPDGVecType operator << (const unsigned short n) const
		{
			return MPDGVecType(altitude << n, mindex);
		}
		inline MPDGVecType operator >> (const unsigned short n) const
		{
			return MPDGVecType(altitude >> n, mindex);
		}

		inline bool operator < (const MPDGVecType & b) const
			{ return mindex < b.mindex || (mindex == b.mindex && altitude < b.altitude); }
		inline bool operator > (const MPDGVecType & b) const
			{ return mindex > b.mindex || (mindex == b.mindex && altitude > b.altitude); }

		inline bool operator <= (const MPDGVecType & b) const
			{ return mindex <= b.mindex || (mindex == b.mindex && altitude <= b.altitude); }
		inline bool operator >= (const MPDGVecType & b) const
			{ return mindex >= b.mindex || (mindex == b.mindex && altitude >= b.altitude); }
	};

	template< typename T >
	struct MPDGVecTypeNumerics
	{ 
		typedef long double UP; 
		typedef double POS;

		static inline geoworld::MPDGVecType< T > zero () { return geoworld::MPDGVecType< T > (); }
		static inline geoworld::MPDGVecType< T > step() { return geoworld::MPDGVecType< T > (mars::Numerics< T >::step(), 0); }
		static inline geoworld::MPDGVecType< T > min () { return geoworld::MPDGVecType< T >(std::numeric_limits< T >::lowest(), std::numeric_limits< unsigned short >::lowest()); }
		static inline geoworld::MPDGVecType< T > max () { return geoworld::MPDGVecType< T >(std::numeric_limits< T >::max(), std::numeric_limits< unsigned short >::max()); }
	};

	template< typename E, typename T >
	class ElementalFieldView : public mars::ScalarFieldView< MPDGVecType< T >, MPDGVecTypeNumerics< T > >
	{
	public:
		typedef std::vector< E > ElementList;
        
		using ScalarFieldView = typename mars::ScalarFieldView< MPDGVecType< T >, MPDGVecTypeNumerics< T > >;        
        using Offset = typename ScalarFieldView::Offset;

	private:
		ElementList * _pElements;

	protected:
		ElementalFieldView (mars::MemoryBlock< MPDGVecType< T > > * pBlock, const mars::BBox< unsigned short > & bbox, ElementList * pElemList, const mars::GridType engt = mars::Normal)
			: _pElements(pElemList), ScalarFieldView (pBlock, bbox, engt) {}
		ElementalFieldView (mars::MemoryBlock< MPDGVecType< T > > * pBlock, const unsigned short nLeft, const unsigned short nTop, const unsigned short nWidth, const unsigned short nHeight, ElementList * pElemList, const mars::GridType engt = mars::Normal)
			: _pElements(pElemList), ScalarFieldView (pBlock, nLeft, nTop, nWidth, nHeight, engt) {}

		// Reserved for use by the ElementalField and subclasses, used to duplicate the memory block and element list
		ElementalFieldView (const ElementalFieldView & copy, mars::MemoryBlock< E > * const pBlock, ElementList * const pElements)
			: ScalarFieldView(copy, pBlock), _pElements(pElements) {}

		// Copy constructor
		ElementalFieldView (const ElementalFieldView & copy)
			: ScalarFieldView(copy), _pElements(copy._pElements) {}

		ElementList * elements() { return _pElements; }

		static ElementalFieldView< E, T > createHosted (const unsigned short nWidth, const unsigned short nHeight, const mars::GridType engt) 
			{ return ElementalFieldView< E, T > (new mars::MemoryBlock< MPDGVecType< T > > (nWidth, nHeight), 0, 0, nWidth, nHeight, new ElementList(), engt); }

		mars::Range< MPDGVecType< T >, MPDGVecTypeNumerics< T > > rangeImpl (const int x0, const int y0, const int xN, const int yN) const
		{
			T 
				fAltMax = std::numeric_limits <T>::min(),
				fAltMin = std::numeric_limits <T>::max();
			unsigned short
				nIdxMax = std::numeric_limits <unsigned short>::lowest(),
				nIdxMin = std::numeric_limits <unsigned short>::max();

			int x, y;
			int iti;
			
#ifdef _OPENMP
#define RANGEIMPL_THREADCOUNT 4
#else
#define RANGEIMPL_THREADCOUNT 1
#endif
			unsigned short vnIdxMin4[RANGEIMPL_THREADCOUNT], vnIdxMax4[RANGEIMPL_THREADCOUNT];
			T vfAltMin4[RANGEIMPL_THREADCOUNT], vfAltMax4[RANGEIMPL_THREADCOUNT];

			for (iti = 0; iti < RANGEIMPL_THREADCOUNT; ++iti)
			{
				vfAltMin4[iti] = fAltMin;
				vfAltMax4[iti] = fAltMax;
				vnIdxMin4[iti] = nIdxMin;
				vnIdxMax4[iti] = nIdxMax;
			}

#ifdef _OPENMP
	#pragma omp parallel for private (x, y) num_threads(RANGEIMPL_THREADCOUNT)
#endif
			for (x = x0; x <= xN; ++x)
			{
				const int ti = 
#ifdef _OPENMP
					omp_get_thread_num();
#else
					0;
#endif
				for (y = y0; y <= yN; ++y)
				{
					const MPDGVecType< T > & elem = ScalarFieldView::get(x, y);

					if (elem.altitude > vfAltMax4[ti])
						vfAltMax4[ti] = elem.altitude;
					if (elem.altitude < vfAltMin4[ti])
						vfAltMin4[ti] = elem.altitude;

					if (elem.mindex > vnIdxMax4[ti])
						vnIdxMax4[ti] = elem.mindex;
					if (elem.mindex < vnIdxMin4[ti])
						vnIdxMin4[ti] = elem.mindex;
				}
			}

			for (iti = 0; iti < RANGEIMPL_THREADCOUNT; ++iti)
			{
				if (vfAltMax4[iti] > fAltMax)
					fAltMax = vfAltMax4[iti];
				if (vfAltMin4[iti] < fAltMin)
					fAltMin = vfAltMin4[iti];

				if (vnIdxMax4[iti] > nIdxMax)
					nIdxMax = vnIdxMax4[iti];
				if (vnIdxMin4[iti] < nIdxMin)
					nIdxMin = vnIdxMin4[iti];
			}

			return mars::Range< MPDGVecType< T >, MPDGVecTypeNumerics < T > > (MPDGVecType< T >(fAltMin, nIdxMin), MPDGVecType< T >(fAltMax, nIdxMax));
		}

		friend class mars::ScalarFieldView< MPDGVecType< T >, MPDGVecTypeNumerics< T > >;

	public:
		typedef E Element;
		typedef ElementalFieldView< E, T > View;
		typedef T AltitudePrecision;
        
        ElementalFieldView(ElementalFieldView && move) : ScalarFieldView(move), _pElements(move._pElements) {}

		ElementalFieldView < E, T > * createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom);
		const ElementalFieldView < E, T > * createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom) const;

		inline const E & getElement (const signed int x, const signed int y) const
		{ 
			return _pElements->at(
				ScalarFieldView::get(
					x + static_cast< const ScalarFieldView * const > (this) ->left, 
					y + static_cast< const ScalarFieldView * const > (this) ->top
				).mindex
			); 
		}
		inline size_t elementCount () const 
			{ return _pElements->size(); }
		inline const E & getElement (const size_t nIdx) const
			{ return _pElements->at(nIdx); }
		inline E & getElement (const size_t nIdx)
			{ return _pElements->at(nIdx); }

		signed int findElementIndex (const E & elem) const
		{
			unsigned short nIdx = 0;
			for (typename ElementList::const_iterator i = _pElements->begin(); i != _pElements->end(); ++i, ++nIdx)
				if (*i == elem)
					return nIdx;

			return -1;
		}

		size_t countViewElements () const
		{
			std::set< unsigned short >	setMIdx;

			for (size_t j = 0; j < ScalarFieldView::height; ++j)
				for (size_t i = 0; i < ScalarFieldView::width; ++i)
					setMIdx.insert(ScalarFieldView::get(i, j).mindex);

			return setMIdx.size();
		}

		void fetchViewElements (std::set< E > & st) const
		{
			for (size_t j = 0; j < ScalarFieldView::height; ++j)
				for (size_t i = 0; i < ScalarFieldView::width; ++i)
					st.insert(getElement(i, j));
		}
		void fetchViewIndices (std::set< unsigned short > & st) const
		{
			for (size_t j = 0; j < ScalarFieldView::height; ++j)
				for (size_t i = 0; i < ScalarFieldView::width; ++i)
					st.insert(ScalarFieldView::operator () (i, j) .mindex);
		}

		inline typename ElementalFieldView< E, T >::ElementList::const_iterator beginElements() const 
			{ return _pElements->begin(); }
		inline typename ElementalFieldView< E, T >::ElementList::const_iterator endElements() const 
			{ return _pElements->end(); }

		template< typename L >
		inline void reg (const typename L::const_iterator iBegin, const typename L::const_iterator iEnd, const L * = NULL)
			{ _pElements->insert(_pElements->end(), iBegin, iEnd); }

		inline void reg (const E & e)
			{ _pElements->push_back(e); }

		inline void clearRegs ()
			{ _pElements->clear(); }

		inline const T & altitude (const signed int x, const signed int y) const 
			{ return ScalarFieldView::get(x, y).altitude; }
		inline T & altitude (const signed int x, const signed int y)
			{ return ScalarFieldView::get(x, y).altitude; }

		inline void operator += (const T & n)
			{ ScalarFieldView::operator += (MPDGVecType< T >(n, 0)); }
		inline void operator -= (const T & n)
			{ ScalarFieldView::operator -= (MPDGVecType< T >(n, 0)); }

		inline virtual mars::Range< MPDGVecType< T >, MPDGVecTypeNumerics< T > > range () const
			{ return rangeImpl (0, 0, ScalarFieldView::right - ScalarFieldView::left, ScalarFieldView::bottom - ScalarFieldView::top); }

		inline void operator << (const ElementalFieldView< E, T > & copy)
		{
			_pElements->clear();
			_pElements->insert(_pElements->end(), copy.beginElements(), copy.endElements());
			ScalarFieldView::operator << (copy);
		}
		inline void operator >> (ElementalFieldView< E, T > & dest) const
			{ dest << *this; }

		inline void operator >> (Offset & dest) const
			{ ScalarFieldView::operator >> (dest); }
		inline void operator << (const Offset & src)
			{ ScalarFieldView::operator << (src); }
	};
	
	template< typename E, typename T >
	ElementalFieldView < E, T > * ElementalFieldView< E, T >::createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom)
		{ return ScalarFieldView::template createTypedView< ElementalFieldView< E, T > > (nLeft, nTop, nRight, nBottom, _pElements); }

	template< typename E, typename T >
	const ElementalFieldView < E, T > * ElementalFieldView< E, T >::createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom) const
		{ return ScalarFieldView::template createTypedView< ElementalFieldView< E, T > > (nLeft, nTop, nRight, nBottom, _pElements); }

	template <typename E, typename T> class VarriedMPDG;
	
	template <typename E, typename T>
	class ElementalField : public ElementalFieldView< E, T >
	{
	public:
        using MemoryBlock = mars::MemoryBlock< MPDGVecType< T > >;
        using ScalarFieldView = mars::ScalarFieldView< MPDGVecType< T >, MPDGVecTypeNumerics< T > >;
        using VarriedMPDGType = VarriedMPDG< E, T >;
        using ElementList = typename ElementalFieldView< E, T >::ElementList;
        
		ElementalField (const unsigned short width, const unsigned short height, const mars::ScalarField< MPDGVecType< T > > & copy, const unsigned short x, const unsigned short y, const mars::GridType engt = mars::Normal)
			: ElementalFieldView< E, T >(VarriedMPDGType::createHosted(width, height, engt)) 
			{ blitFrom (copy, x, y); }

		ElementalField (const unsigned short width, const unsigned short height, const mars::GridType engt = mars::Normal)
			: ElementalFieldView< E, T >(ElementalFieldView< E, T >::createHosted(width, height, engt)) 
			{ memset(this->data(), 0, this->block() -> bytesize); }
		ElementalField (const ElementalFieldView< E, T > & copy)
			: ElementalFieldView< E, T >(copy, new MemoryBlock (copy.getMemoryBlock()), new ElementList (*copy._pElements)) {}

		virtual ~ElementalField()
			{ delete this->elements(); }

		inline std::istream & operator << (std::istream & ins)
			{ return *this->block() << ins; }
		inline std::ostream & operator >> (std::ostream & outs) const
			{ return *this->block() >> outs; }
	
		inline void operator >> (typename ScalarFieldView::Offset & dest) const
			{ ElementalFieldView< E, T >::operator >> (dest); }
		inline void operator << (const typename ScalarFieldView::Offset & src)
			{ ElementalFieldView< E, T >::operator << (src); }
	};

	template <typename E, typename T>
	class VarriedMPDG : public mars::MidpointDisplacementGrid< MPDGVecType< T >, MPDGVecTypeNumerics< T > >, public ElementalFieldView< E, T >
	{
	private:
		const enum StorageMode
		{
			Concrete,
			Virtual
		} _ensmStorageMode;

	protected:
        using MPDG = mars::MidpointDisplacementGrid< MPDGVecType< T >, MPDGVecTypeNumerics< T > >;
        using Elems = ElementalFieldView< E, T >;

		class SeedFunctor : public MPDG::PinkFn
		{
		private:
			const unsigned short _maxidx;
			float _base;

		public:
			inline SeedFunctor (unsigned short maxidx, const unsigned short nDim, const float fCoarseness)
				: MPDG::PinkFn(nDim, fCoarseness), _maxidx(maxidx) {}

			inline void iteration (const unsigned short nStride, const unsigned short c)
			{
				MPDG::PinkFn::iteration(nStride, c);
				_base = this->computeRandomAmount(0) / 2;
			}

			inline MPDGVecType< T > operator () (int & state, const int x, const int y, const MPDGVecType< T > & a, const MPDGVecType< T > & b, const MPDGVecType< T > & c)
			{
				MPDGVecType< T > result;

				result.mindex = mars::RAND(_maxidx);
				result.altitude = static_cast< T > (this->applyFlux (_base));
				return result;
			}
			inline MPDGVecType< T > operator () (int & state, const int x, const int y, const MPDGVecType< T > & a, const MPDGVecType< T > & b, const MPDGVecType< T > & c, const MPDGVecType< T > & d)
			{
				return operator() (state, x, y, a, b, c);
			}
		};

		struct LocalFunctorState
		{
			T sumalt[3];
			unsigned char ctalt[3];
			MPDGVecType< T > result;
		};

		class LocalFunctor : public MPDG::template AbstractFunctor<LocalFunctorState>
		{
		public:
			inline LocalFunctor (const unsigned short nDim, const float fCoarseness)
				: MPDG::template AbstractFunctor<LocalFunctorState>(nDim, fCoarseness) {}

			MPDGVecType< T > operator () (LocalFunctorState & state, const int x, const int y, const MPDGVecType< T > & a, const MPDGVecType< T > & b, const MPDGVecType< T > & c)
			{
				// Count cells matching 'a'
				state.ctalt[0] = 0;
				state.sumalt[0] = a.altitude;
				if (a.mindex == b.mindex)
				{
					state.sumalt[0] += b.altitude;
					++state.ctalt[0];
				}
				if (a.mindex == c.mindex)
				{
					state.sumalt[0] += c.altitude;
					++state.ctalt[0];
				}

				// Count cells matching 'b'
				state.ctalt[1] = 0;
				state.sumalt[1] = b.altitude;
				if (b.mindex == c.mindex)
				{
					state.sumalt[1] += c.altitude;
					++state.ctalt[1];
				}

				// Flip a coin to pick between a dominant alpha or auxiliary beta cell
				const unsigned int pick = mars::RAND(3);

				if (pick < 2)
					state.result.altitude = this->applyFlux(state.sumalt[pick] / (state.ctalt[pick] + 1));
				else
					state.result.altitude = this->applyFlux(c.altitude);

				switch (pick)
				{
				case 0: 
					state.result.mindex = a.mindex;
					break;
				case 1:
					state.result.mindex = b.mindex;
					break;
				case 2:
					state.result.mindex = c.mindex;
					break;
				}

				return state.result;
			}
			MPDGVecType< T > operator () (LocalFunctorState & state, const int x, const int y, const MPDGVecType< T > & a, const MPDGVecType< T > & b, const MPDGVecType< T > & c, const MPDGVecType< T > & d)
			{
				// Count cells matching 'a'
				state.ctalt[0] = 0;
				state.sumalt[0] = a.altitude;
				if (a.mindex == b.mindex)
				{
					state.sumalt[0] += b.altitude;
					++state.ctalt[0];
				}
				if (a.mindex == c.mindex)
				{
					state.sumalt[0] += c.altitude;
					++state.ctalt[0];
				}
				if (a.mindex == d.mindex)
				{
					state.sumalt[0] += d.altitude;
					++state.ctalt[0];
				}

				// Count cells matching 'b'
				state.ctalt[1] = 0;
				state.sumalt[1] = b.altitude;
				if (b.mindex == c.mindex)
				{
					state.sumalt[1] += c.altitude;
					++state.ctalt[1];
				}
				if (b.mindex == d.mindex)
				{
					state.sumalt[1] += d.altitude;
					++state.ctalt[1];
				}

				// Count cells matching 'c'
				state.ctalt[2] = 0;
				state.sumalt[2] = c.altitude;
				if (c.mindex == d.mindex)
				{
					state.sumalt[2] += d.altitude;
					++state.ctalt[2];
				}

				// Flip a coin to pick between a dominant alpha or auxiliary beta cell

				const unsigned int pick = mars::RAND(4);

				if (pick < 3)
					state.result.altitude = this->applyFlux(state.sumalt[pick] / (state.ctalt[pick] + 1));
				else
					state.result.altitude = this->applyFlux(d.altitude);

				switch (pick)
				{
				case 0: 
					state.result.mindex = a.mindex;
					break;
				case 1:
					state.result.mindex = b.mindex;
					break;
				case 2:
					state.result.mindex = c.mindex;
					break;
				case 3:
					state.result.mindex = d.mindex;
					break;
				}

				return state.result;
			}
		};

		VarriedMPDG (const short nIterCount, ElementalFieldView< E, T > v, const StorageMode ensm)
			: MPDG(nIterCount, v), ElementalFieldView<E, T>(v), _ensmStorageMode(ensm) {}
		
		static VarriedMPDG * createHosted (const short nIterCount, const unsigned short nWidth, const unsigned short nHeight, const mars::GridType engt = mars::Normal)
			{ return new VarriedMPDG(nIterCount, ElementalFieldView< E, T >::createHosted(nWidth, nHeight, engt), Concrete); }

	public:
   		using ScalarFieldView = typename mars::ScalarFieldView< MPDGVecType< T >, MPDGVecTypeNumerics< T > >;    

		virtual ~VarriedMPDG()
		{
			switch (_ensmStorageMode)
			{
			case Concrete: delete this->elements(); break;
			}
		}

		inline static VarriedMPDG * createInstance (const unsigned short nMinWidth, const unsigned short nMinHeight, const mars::GridType engt)
		{
			assert(nMinHeight > 0 && nMinWidth > 0);
			return createHosted (
				MPDG::calcIterCount(std::min(nMinHeight, nMinWidth)), 
				MPDG::dimensionFor(nMinWidth, engt),
				MPDG::dimensionFor(nMinHeight, engt),
				engt
			);
		}
		inline static VarriedMPDG * createInstance (Elems & view)
		{
			return new VarriedMPDG (
				calcIterCount(std::min(view.width, view.height) - (view.gridtype == mars::Normal ? 1 : 0)), 
				view,
				Virtual
			);
		}

		inline virtual mars::Range< MPDGVecType< T >, MPDGVecTypeNumerics< T > > range () const
			{ return Elems::range(); }
		inline virtual void clear()
			{ Elems::clear(); }
		inline void stitch(const int nStitch)
			{ Elems::stitch(nStitch); }
		void stitch(const int nStitch, const ScalarFieldView * pOther)
			{ Elems::stitch(nStitch, pOther); }
		void stitchAll(const int nStitch)
			{ Elems::stitchAll(nStitch); }
		inline void stitchAll(const int nStitch, const ScalarFieldView * pLeft, const ScalarFieldView * pTop, const ScalarFieldView * pRight, const ScalarFieldView * pBottom)
			{ Elems::stitchAll(nStitch, pLeft, pTop, pRight, pBottom); }
			
        // TODO: Type hierarchy seems broken (bogus diamond pattern)
		inline void operator -= (const T & n)
		{
            Elems::operator -= (n);
		}
	
		void seed (unsigned short nIterations, const float fCoarseNess, const signed int nFlux = -1)
		{
			SeedFunctor f (Elems::elementCount(), nFlux < 0 ? this->mindim : static_cast< unsigned short > (nFlux), fCoarseNess);

			this->iterations(nIterations, f);
		}

		void run (const float fCoarseNess, const int nStitch = 0)
		{
			LocalFunctor f (this->mindim, fCoarseNess);

			this->iterations (this->itercount - this->index(), f, nStitch);
		}
	};
}

#endif