/* Copyright (c) 2018
   Charles Jordan (Hokkaido University, Sapporo, Japan),
   Michael Joswig, Lars Kastner (Technische Universitaet Berlin, Germany)

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 3, or (at your option) any
   later version: http://www.gnu.org/licenses/gpl.txt.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
-----------------------------------------------------------------------------
*/

#ifndef AUGMENTED_POINT_CONFIGURATION_CYCLIC_HH
#define AUGMENTED_POINT_CONFIGURATION_CYCLIC_HH

#include "CommandlineOptions.hh"
#include "Chirotope.hh"
#include "Symmetry.hh"
#include "PointConfiguration.hh"
#include "mptopcom_typedefs.hh"

#include "Flip.hh"
#include "TriangNode.hh"
#include "mptopcom_typedefs.hh"


#include <polymake/Main.h>
#include <polymake/Matrix.h>
#include <polymake/Set.h>
#include <polymake/Bitset.h>
#include <polymake/SparseVector.h>
#include <polymake/IncidenceMatrix.h>
#include <polymake/Integer.h>
#include <polymake/Map.h>
#include <polymake/Vector.h>
#include <polymake/PowerSet.h>
#include <bitset>


template<typename checkFine, typename smallGroup>
class AugmentedPointConfigurationCyclic {

   private:
      const PointConfiguration&           points;
      const Chirotope&                    chiro;
      const SymmetryGroup&                symmetries;
      // const gkz_type                      gkz;
      
      // Strategy:
      // circuits holds all d+2 subsets. These correspond to the circuits.
      // Then we have two maps:
      // Map<Circuit, UpperTriang>
      // Map<Circuit, LowerTriang>
      // Given a triangulation, we check for every circuit whether the upper or
      // lower triang is contained. If it is, we found a flip.
      const int d;
      const int n;
      const int faceThreshold;
      pm::Map<pmCircuit, pmSimpComp> upperHulls;
      pm::Map<pmCircuit, pmSimpComp> lowerHulls;
      pm::Map<pmCircuit, Flip> topcomDownflips;
      pm::Map<pmCircuit, Flip> topcomUpflips;
      const pmSimplex empty;
      
      void init_polymake(){
         static polymake::Main pmm("none");
         pmm.set_application("polytope");
      }

      void compute_hulls(const pmCircuit& circuit){
         pmSimpComp upper, lower;
         pmSimplex all(circuit);
         int j = 0;
         for(const auto& i : circuit){
            if(j%2 != 0){
               upper += all-i;
            } else {
               lower += all-i;
            }
            j++;
         }
         upperHulls.insert(circuit, upper);
         lowerHulls.insert(circuit, lower);
         topcomDownflips.insert(circuit, pmFlip2Topcom(upper, lower));
         topcomUpflips.insert(circuit, pmFlip2Topcom(lower, upper));
      }
      
      pmSimpComp triangNode2SetSetInt(const TriangNode& t) const {
         pmSimpComp result;
         for(const auto& triangle : t){
            result += static_cast<const pmSimplex&>(triangle);
            // pmSimplex tmp;
            // for(const auto& i: triangle){
            //    tmp += i;
            // }
            // result += tmp;
         }
         return result;
      }

      Flip pmFlip2Topcom(const pmSimpComp& before, const pmSimpComp& after) const {
         SimplicialComplex bsc, asc;
         for(const auto& s: before){
            Simplex tmp;
            for(const auto& i: s){
               tmp += i;
            }
            bsc += tmp;
         }
         for(const auto& s: after){
            Simplex tmp;
            for(const auto& i: s){
               tmp += i;
            }
            asc += tmp;
         }
         Pair<SimplicialComplex, SimplicialComplex> pss(bsc, asc);
         Flip flip;
         flip.first = bsc;
         flip.second = asc;
         return flip;
      }


   public:
      typedef smallGroup smallGroupBool;
      typedef checkFine checkFineBool;
      typedef std::true_type trivialSymmetry;
      typedef pmFlipHash<pmCircuit, Flip> pmFlipContainer;

      AugmentedPointConfigurationCyclic(
            const PointConfiguration& points_in,
            const Chirotope&          chiro_in,
            const SymmetryGroup&      symmetries_in) :
         points(points_in), chiro(chiro_in),
         symmetries(symmetries_in),
         empty(),
         n(points_in.no()), d(points_in.rank()-1),
         faceThreshold((d+2)/2-1)
         {
            init_polymake();
            for(const auto& circuit: pm::all_subsets_of_k(pm::range(0, n-1), d+2)){
              compute_hulls(pmCircuit(pm::Set<int>(circuit)));
            }
         }
      
      bool check_triang_equality(const TriangNode& t1, const TriangNode& t2) const{
         return t1==t2;
      }

      const PointConfiguration& get_points() const{
         return points;
      }

      const Chirotope& get_chiro() const{
         return chiro;
      }

      const SymmetryGroup& get_symmetries() const{
         return symmetries;
      }


      std::string get_extra_info(const TriangNode& triang) const {
         std::string result;
         return result;
      }
      
      std::string get_triangulation_data(const TriangNode& triang) const {
         std::ostringstream bos;
         bos << triang;
         return bos.str();
      }

      void printInitialTriang(CerrWrapper<false>& ost, const TriangNode& triang) const {
         ost << "Initial triang: " << triang << std::endl;
      }

      void mutate(const Flip& flip, const TriangNode& target, pmFlipContainer& upflips, pmFlipContainer& downflips) const {
         pmSimplex support;
         for(const auto& simplex: flip.first){
            for(const auto& i: simplex){
               support += i;
            }
         }
         remove_destroyed_flips(support, upflips);
         remove_destroyed_flips(support, downflips);
         populate_flips_inner<true>(target, upflips, downflips, support);
      }

      void remove_destroyed_flips(const pmSimplex& support, pmFlipContainer& flips) const{
         std::list<pmCircuit> toRemove;
         for(const auto& entry: flips){
            if(( support * pmSimplex(entry.first)).size() > d){
               toRemove.push_back(pmCircuit(entry.first));
            }
         }
         for(const auto& circuit: toRemove){
            flips.erase(circuit);
         }
      }
      
      void populate_flips(const TriangNode& source, pmFlipContainer& upflips, pmFlipContainer& downflips) const {
         populate_flips_inner<false>(source, upflips, downflips, empty);
      }

      template<bool mutationBool>
      void populate_flips_inner(const TriangNode& source, pmFlipContainer& upflips, pmFlipContainer& downflips, const pmSimplex& support) const {
         pmSimpComp pmTriang = triangNode2SetSetInt(source);
         for(const auto& simplex: pmTriang){
            std::list<pmCircuit> candidatesUp, candidatesDown;
            // This applies Gale's eveness criterion.
            if(d%2 == 0){
               get_candidates<mutationBool>(simplex, candidatesDown, candidatesUp, support);
            } else {
               get_candidates<mutationBool>(simplex, candidatesUp, candidatesDown, support);
            }
            for(const auto& circuit: candidatesUp){
               if(isApplicableFlip(lowerHulls[circuit], pmTriang)){
                  upflips.insert(circuit, topcomUpflips[circuit]);
               }
            }
            for(const auto& circuit: candidatesDown){
               if(isApplicableFlip(upperHulls[circuit], pmTriang)){
                  downflips.insert(circuit, topcomDownflips[circuit]);
               }
            }
         }
      }
      
      // mutationBool indicates whether we come from a flip and just want to
      // get the new flips that are generated by applying the flip.
      template<bool mutationBool>
      void get_candidates(const pmSimplex& simplex, std::list<pmCircuit>& type1, std::list<pmCircuit>& type2, const pmSimplex& support) const {
         if(!mutationBool){
            get_candidates_from_scratch(simplex, type1, type2);
         } else {
            get_candidates_from_support(simplex, type1, type2, support);
         }
      }

      void get_candidates_from_support(pmSimplex simplex, std::list<pmCircuit>& type1, std::list<pmCircuit>& type2, const pmSimplex& support) const {
         int intersectionSize = (simplex * support).size();
         if(intersectionSize < d){ return;}
         // For small intersections we only have the choice to add a point from
         // the support.
         bool smallIntersection = intersectionSize == d;
         int last = simplex.back();
         if(last < n-1){
            pmSimplex newElements = pmSimplex(pm::Set<int>(pm::range(last+1, n-1)));
            if(smallIntersection){ newElements *= support; }
            for(int i : newElements){
               simplex.push_back(i);
               type1.push_back(simplex);
               simplex.pop_back();
            }
         }
         simplex.pop_back();
         int secondlast = simplex.back();
         if(secondlast < last-1){
            pmSimplex newElements = pmSimplex(pm::Set<int>(pm::range(secondlast+1, last-1)));
            if(smallIntersection){ newElements *= support; }
            for(int i : newElements){
               simplex.push_back(i);
               simplex.push_back(last);
               type2.push_back(simplex);
               simplex.pop_back();
               simplex.pop_back();
            }
         }
      }


      void get_candidates_from_scratch(pmSimplex simplex, std::list<pmCircuit>& type1, std::list<pmCircuit>& type2) const {
         int last = simplex.back();
         if(last < n-1){
            for(const auto i : pm::range(last+1, n-1)){
               simplex.push_back(i);
               type1.push_back(simplex);
               simplex.pop_back();
            }
         }
         simplex.pop_back();
         int secondlast = simplex.back();
         if(secondlast < last-1){
            for(const auto i : pm::range(secondlast+1, last-1)){
               simplex.push_back(i);
               simplex.push_back(last);
               type2.push_back(simplex);
               simplex.pop_back();
               simplex.pop_back();
            }
         }
      }

      bool isApplicableFlip(const pmSimpComp& sourceSimplices, const pmSimpComp triang) const {
         for(const auto& s : sourceSimplices){
            if(!triang.contains(s)){ return false;}
         }
         return true;
      }
      
      int size() const {
         return lowerHulls.size();
      }

      void init_check(const TriangNode& t) const {
      }

      TriangNode canonicalize(const TriangNode& triang) const {
         return triang;
      }

};

#endif

