/* 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.
-----------------------------------------------------------------------------
*/

///////////////////////////////////////////////////////////////////////////////
// TriangulationOptimizer.hh
//
// This object computes the canonical (i.e. lex largest) class of a
// triangulation.
//
// The basic principle is the following: The GroupDecomposer decomposes the
// group in a way that allows one to find the canonical representative with a
// fast DFS. The TriangulationOptimizer object encapsulates this DFS.
///////////////////////////////////////////////////////////////////////////////

#ifndef TRIANGULATIONOPTIMIZER_HH
#define TRIANGULATIONOPTIMIZER_HH

#include "GroupDecomposer.hh"
#include "mptopcom_typedefs.hh"
#include <polymake/Bitset.h>
#include <utility>

class PackagedTriangulation {
   private:
      triang_vector_type tvec;
      const gkz_type& gkz;
   public:
      PackagedTriangulation(const TriangNode& t_in, const gkz_type& gkz_in):
         gkz(gkz_in), tvec(gkz_in.triangulationAsVector(t_in))
      {}
      
      PackagedTriangulation(const triang_vector_type& t, const gkz_type& gkz_in):
         gkz(gkz_in), tvec(t)
      {}

      const pm::Set<int> get_support(const int& fixed, const pm::Set<int>& in, bool& needsIdentity) const {
         pm::Set<int> intersection(tvec*in);
         if(intersection.empty()){
            needsIdentity = true;
            return in;
         } else {
            if(tvec.contains(fixed)){
               needsIdentity = true;
            }
            return intersection;
         }
      }

      bool operator >(const PackagedTriangulation& other) const {
         return other.tvec > tvec;
      }

      PackagedTriangulation& operator= (const PackagedTriangulation& other) {
         tvec = other.tvec;
         return *this;
      }

      PackagedTriangulation mutate(const pm::Array<int>* g) const {
         Symmetry tmp(array2symmetry(*g));
         triang_vector_type next;
         for(const auto& e: tvec){
            next += gkz.get_basis2int()[tmp(gkz.get_int2basis()[e])];
         }
         return PackagedTriangulation(next, gkz);
      }
};

class ReferenceCore {
   private:
      const pm::Array<int>* identity;
      std::map<int, std::map<int, const pm::Array<int>*>> switch_table;
      pm::Map<int, pm::Set<int>> supports;
      int bound;

   public:
      ReferenceCore(std::map<int, std::map<int, const pm::Array<int>*>> st, const pm::Array<int>* id):
         identity(id), switch_table(st), bound(0)
      {
         for(const auto& level : switch_table){
            for(const auto& sublevel : level.second){
               supports[level.first] += sublevel.first;
            }
            bound = level.first > bound ? level.first : bound;
         }
         bound++;
      }

      void extract_switches(int fixed, const pm::Set<int>& desired, std::list<const pm::Array<int>*>& switches) const {
         for(const auto& j : desired){
            switches.push_back(switch_table.at(fixed).at(j));
         }
      }

      int get_bound() const {
         return bound;
      }

      const pm::Array<int>& get_identity() const {
         return *identity;
      }

      bool support_exists(int fixed) const {
         return supports.exists(fixed);
      }

      const pm::Set<int>& get_support(int fixed) const {
         return supports[fixed];
      }

};

class TriangulationOptimizer {
   private:
      TriangNode optimal;
      gkz_vector_type optimal_gkz;
      const GroupDecomposer& gd;
      const gkz_type& gkz;
      Symmetry mult;


      // Find lex largest representative of class that also lex maximizes the
      // gkz vector.
      void optimize(){
         const polymake::group::SwitchTable& ST(gd.get_switch_table());
         std::pair<gkz_vector_type, pm::Array<int>> P(ST.lex_maximize_vector<Scalar>(optimal_gkz));
         mult = Symmetry(array2symmetry(P.second));
         optimal = mult(optimal);
         optimal_gkz = P.first;
         optimize_representative();
      }

      // Every set in result corresponds to a subset of indices of the gkz
      // vector with the same value. These need to be stabilized if we want to
      // change the representative while preserving the gkz vector.
      const signature_vector_type signature() const{
         signature_vector_type result(optimal_gkz.dim());
         pm::Map<Scalar, int> numbers, count;
         for(int i=0; i<result.dim(); i++){
            if(numbers.exists(optimal_gkz[i])){
               result[i] = numbers[optimal_gkz[i]];
               count[optimal_gkz[i]]++;
            } else {
               result[i] = i;
               numbers[optimal_gkz[i]] = i;
               count[optimal_gkz[i]] = 1;
            }
         }
         return result;
      }


      // Find the lex largest representative, while keeping the gkz-vector.
      void optimize_representative(){
         std::map<int, std::map<int, const pm::Array<int>*>> refST;
         signature_vector_type sig = signature();
         bool needsOpt = gd.populate_partition_stabilizer(sig, refST);
         if(needsOpt){
            ReferenceCore RC(std::move(refST), &(gd.get_identity()));
            PackagedTriangulation packaged(optimal, gkz);
            polymake::group::switchtable::Optimizer<ReferenceCore, PackagedTriangulation> opt(RC, packaged);
            opt.optimize();
            Symmetry g(array2symmetry(opt.get_optimal().second));
            optimal = g(optimal);
            mult = g * mult;
         }
      }



   public:
      TriangulationOptimizer(const TriangNode& t, 
         const gkz_vector_type& gkzt, 
         const GroupDecomposer& gd_in, 
         const gkz_type& gkz_in) :
         optimal(t), optimal_gkz(gkzt),
         gd(gd_in), gkz(gkz_in),
         mult(array2symmetry(gd.get_identity()))
   { }
      
      const TriangNode& get_optimal(){
         optimize();
         return optimal;
      }

      const gkz_vector_type& get_optimal_gkz(){
         return optimal_gkz;
      }

      const Symmetry& get_permutation() {
         return mult;
      }

};

#endif

