#ifndef CANTOR_FAULT_FACTORY_HH
#define CANTOR_FAULT_FACTORY_HH

#include <math.h>

#include <dune/faultnetworks/utils/debugutils.hh>
#include <dune/faultnetworks/faultfactories/oscunitcube.hh>
#include <dune/faultnetworks/levelinterfacenetwork.hh>
#include <dune/faultnetworks/interfacenetwork.hh>

#include <dune/grid/common/mcmgmapper.hh>


class CantorIndexHierarchy {
	public:
		typedef std::map<std::pair<int,int>, bool> LevelCantorIndices;
		
	private:
		const size_t maxLevel_;
		std::vector<size_t> levelN_;
		std::vector<LevelCantorIndices> cantorIndexHierarchy_;

	    void prolongIndices(size_t currentLevel, size_t newLevel) {
			const LevelCantorIndices& indices = cantorIndexHierarchy_[currentLevel];
			LevelCantorIndices& newIndices = cantorIndexHierarchy_[newLevel];
			
            size_t offset = levelN_[currentLevel];
			
			std::map<std::pair<int,int>, bool>::const_iterator it = indices.begin();
			std::map<std::pair<int,int>, bool>::const_iterator endIt = indices.end();
			for (; it!=endIt; ++it) {
				int xID = it->first.first;
				int yID = it->first.second;
				newIndices[std::make_pair(xID, yID)] = true;
				newIndices[std::make_pair(xID+offset, yID)] = true;
				newIndices[std::make_pair(xID, yID+offset)] = true;
			}

            size_t doubleOffset = 2*offset;
            for (size_t i=offset+1; i<=doubleOffset; i=i+2) {
                newIndices[std::make_pair(doubleOffset, i)] = true;
                newIndices[std::make_pair(i, doubleOffset)] = true;
            }
		}

        void print(const LevelCantorIndices& indices) const {
            std::cout << "LevelCantorIndices: " << std::endl;

            std::map<std::pair<int,int>, bool>::const_iterator it = indices.begin();
            std::map<std::pair<int,int>, bool>::const_iterator endIt = indices.end();
            for (; it!=endIt; ++it) {
                std::cout << "(" << it->first.first << ", " << it->first.second << ")"<< std::endl;
            }
        }

	public:
		CantorIndexHierarchy(size_t maxLevel) :
			maxLevel_(maxLevel)
		{
			levelN_.resize(maxLevel_+1);
			cantorIndexHierarchy_.resize(maxLevel_+1);
			
			// init levelCantorIndices on level 0
			LevelCantorIndices& initCantorIndices = cantorIndexHierarchy_[0];
			initCantorIndices[std::make_pair(0, 1)] = true;
			initCantorIndices[std::make_pair(1, 0)] = true;
			initCantorIndices[std::make_pair(1, 2)] = true;
			initCantorIndices[std::make_pair(2, 1)] = true;

			for (size_t i=0; i<maxLevel_; i++) {
				levelN_[i] = std::pow(2, i+1);
				prolongIndices(i, i+1);
			}

            levelN_[maxLevel_] = std::pow(2, maxLevel_+1);
		}
		
		const LevelCantorIndices& levelCantorIndices(size_t i) const {
				return cantorIndexHierarchy_[i];
		}
		
        size_t levelN(const size_t i) const {
			return levelN_[i];
		}
};


template <class GridType>
class CantorFaultFactory
{
    //! Parameter for mapper class
    template<int dim>
    struct FaceMapperLayout
    {
        bool contains (Dune::GeometryType gt)
        {
            return gt.dim() == dim-1;
        }
    };

    protected:
        typedef OscUnitCube<GridType, 2> GridOb;
        static const int dimworld = GridType::dimensionworld ;
        static const int dim = GridType::dimension;
        typedef typename GridType::ctype ctype;
        typedef typename GridType::LevelGridView GV;

        typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;

    private:
        const std::map<size_t, size_t>& levelToCantorLevel_;
        const int coarseResolution_;
        const size_t maxLevel_;
        const size_t maxCantorLevel_;
        const CantorIndexHierarchy cantorIndexHierarchy_;
        const int coarseGridN_;

        std::shared_ptr<GridType> grid_;
        std::vector<double> levelResolutions_;
        
        std::shared_ptr<InterfaceNetwork<GridType>> interfaceNetwork_;   

private:
    typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type almost_equal(ctype x, ctype y, int ulp) const {
        return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp || std::abs(x-y) < std::numeric_limits<ctype>::min();
    }

    bool computeID(Dune::FieldVector<ctype, dimworld> vertex, const size_t gridN, std::pair<size_t, size_t>& IDs) const {
        ctype xID = vertex[0]*gridN;
        ctype yID = vertex[1]*gridN;

        ctype x = 0;
        ctype y = 0;

        bool xIsInt = almost_equal(std::modf(xID, &x), 0.0, 2);
        bool yIsInt = almost_equal(std::modf(yID, &y), 0.0, 2);

        if (xIsInt && yIsInt) {
            IDs = std::make_pair((size_t) x, (size_t) y);
            return true;
        }

        if (xIsInt) {
            y = std::ceil(yID);
            if ((size_t) y % 2==0)
                y = y-1;
            IDs = std::make_pair((size_t) x, (size_t) y);
            return true;
        }

        if (yIsInt) {
            x = std::ceil(xID);
            if ((size_t) x % 2==0)
                x = x-1;
            IDs = std::make_pair((size_t) x, (size_t) y);
            return true;
        }

        return false;
    }

public:
        //setup 
    CantorFaultFactory(const std::map<size_t, size_t>& levelToCantorLevel, const int coarseResolution, const size_t maxLevel, const size_t maxCantorLevel) :
        levelToCantorLevel_(levelToCantorLevel),
        coarseResolution_(coarseResolution),
        maxLevel_(maxLevel),
        maxCantorLevel_(maxCantorLevel),
        cantorIndexHierarchy_(maxCantorLevel_),
        coarseGridN_(std::pow(2, coarseResolution_)),
        interfaceNetwork_(nullptr)
        {   
            Dune::UGGrid<dim>::setDefaultHeapSize(4000);
            GridOb unitCube(coarseGridN_);
            grid_ = unitCube.grid();

            grid_->globalRefine(maxLevel_);

            levelResolutions_.resize(maxLevel_+1);

            // init interface network
            interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_);

            for (size_t i=0; i<=maxLevel_; i++) {
                if (i>0) {
                    interfaceNetwork_->prolongLevel(i-1, i);
                }

                if (levelToCantorLevel_.count(i)) {

                    levelResolutions_[i] = std::pow(2, coarseResolution_+i);
                    const size_t levelResolution = levelResolutions_[i];
                    const size_t cantorResolution = cantorIndexHierarchy_.levelN(levelToCantorLevel_.at(i));
				
                    if (2*levelResolution<cantorResolution)
                        DUNE_THROW(Dune::Exception, "Level " + std::to_string(i) + " does not support cantor set with resolution " + std::to_string(cantorResolution));

                    const typename CantorIndexHierarchy::LevelCantorIndices& levelCantorIndices = cantorIndexHierarchy_.levelCantorIndices(levelToCantorLevel_.at(i));
                    std::shared_ptr<LevelInterfaceNetwork<GV>> levelInterfaceNetwork = interfaceNetwork_->levelInterfaceNetworkPtr(i);
					
                    const GV& gridView = levelInterfaceNetwork->levelGridView();
                    FaceMapper intersectionMapper(gridView);
                    std::vector<bool> intersectionHandled(intersectionMapper.size(),false);

                    for (const auto& elem:elements(gridView)) {
                        for (const auto& isect:intersections(gridView, elem)) {

                            if (intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)])
                                continue;

                            intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)]=true;

                            if (isect.boundary())
                                continue;

                            std::pair<size_t, size_t> intersectionID;
                            bool isAdmissibleID = computeID(isect.geometry().center(), cantorResolution, intersectionID);

                            if (isAdmissibleID && levelCantorIndices.count(intersectionID)) {
                                levelInterfaceNetwork->addIntersection(isect, i);
                            }
                        }
                    }
                }
            }
        }

    /*
    void prolongToAll() {
        // prolong all faults to all subsequent levels
        for (int i=maxLevel_-1; i>=0; i--) {
            if (interfaceNetwork_->size(i)>0) {
                std::set<int> toLevels;
                for (size_t j=i+1; j<=maxLevel_; j++) {
                    toLevels.insert(j);
                }

                interfaceNetwork_->prolongLevelInterfaces(i, toLevels);
            }
        }
        interfaceNetwork_->build();
    }*/

    /*void prolongToAll() {
        // prolong all faults to all subsequent levels
        for (int i=interfaceNetwork_->size()-1; i>=0; i--) {
            interfaceNetwork_->prolongLevelInterfaces(i, maxLevel_);
        }
        interfaceNetwork_->build();
    }*/

    const GridType& grid() const {
        return *grid_;
	}

    InterfaceNetwork<GridType>& interfaceNetwork() {
        return *interfaceNetwork_;
    }
};
#endif