///////////////////////////////////////////////////////////////////////////////
// File: GKZ.hh
// This file contains a class that can compute the GKZ vector of a
// triangulation. Given n points in RR^k, it stores all n choose k determinants
// one can form. These are then used to sum up the GKZ-vector for a
// triangulation.
//
// TODO:
// 1. Better algorithm for initialization, analogously to the algorithm
// building the Chirotope (avoid zeros, use minors already computed)
// 2. Switch to polymake vectors, since these are sparse
///////////////////////////////////////////////////////////////////////////////

#ifndef GKZ_HH
#define GKZ_HH
#include "SimplicialComplex.hh"

#include "CommandlineOptions.hh"

#include "Chirotope.hh"
#include "Flip.hh"
#include "MarkedFlips.hh"
#include "TriangNode.hh"
#include "Flip.hh"

#include <polymake/client.h>
#include <polymake/Map.h>
#include <polymake/Matrix.h>
#include <polymake/Vector.h>
#include <polymake/Set.h>
#include <polymake/linalg.h>
#include <polymake/SparseVector.h>

#include "mptopcom_typedefs.hh"
#include "CerrWrappers.hh"
#include "RegularityCheck.hh"

#include <stdexcept>


// This class computes our GKZ-vectors. At construction it computes the volumes
// of the simplices that one can build (i.e. the determinants of d+1 minors of
// the points). The method computing an actual GKZ-vector of a triangulation is
// called evaluate. Furthermore one can also compute the GKZ-vector of a flip.
template<typename Scalar, typename Vectortype, typename TriangVectortype>
class GKZ {
   private:
      const PointConfiguration& pts;
      const Chirotope& chiro;
      // HashMap mapping simplices to their volumes
      pm::Map<basis_type, Scalar> abs_determinants;
      // Number all possible simplices
      int total_simplices;
      std::set<basis_type> all_bases;
      pm::Map<basis_type, int> basis2int;
      pm::Map<int, basis_type> int2basis;
      // Save the rank = number of points, and dim = ambient dim of pts.
      parameter_type rk, dim;
      // Matrix to populate abs_determinants
      pm::Matrix<Scalar> pmMat;
      pm::Vector<Scalar> pointNorms;
      pm::Set<int> indices;
      Scalar det_of_submat(const basis_type& basis){
         int i=0;
         indices.clear();
         for(basis_type::const_iterator iter = basis.begin(); iter != basis.end(); ++iter){
            indices.insert(*iter);
            i++;
         }
         return pm::det(pmMat.minor(indices, pm::All));
      }

      void read_points(const PointConfiguration& pts){
         pm::Rational x;
         for(int i=0; i<rk; i++){
            for(int j=0; j<dim; j++){
               x.copy_from(pts[i][j].get_mpq_t());
               pmMat(i,j) = pm::convert_to<Scalar>(x);
            }
         }
      }

      void populate_pointNorms(){
         for(int i=0; i<pointNorms.dim(); i++){
            pointNorms[i] = pmMat[i] * pmMat[i] - 1;
         }
      }

      // int abs value det_of_submat
      Scalar abs_det_of_submat(const basis_type& basis){
         // FIXME: We should be much more careful here, but for now this is ok.
         Scalar val = det_of_submat(basis);
         if(val < 0){
            return -val;
         } else {
            return val;
         }
      }
      // Walk over all card dim subsets of pts and save the absolute values of
      // the determinants.
      // FIXME: This method computes many subdeterminants multiple times. Can
      // we use the same optimized algorithm as for the Chirotope in TOPCOM?
      void populate_abs_determinants(int pts_index, int mat_index, basis_type basis){
         if(mat_index == dim){
            Scalar val = abs_det_of_submat(basis);
            if(val != 0){
               all_bases.insert(basis);
               total_simplices++;
               abs_determinants[basis] = val;
            }
         } else if (pts_index >=rk){
            std::cout << "Fell in here!" << std::endl;
            return;
         } else {
            int corank = dim - basis.card();
            for(int i=pts_index; i<rk-corank+1; i++){
               basis += i;
               populate_abs_determinants(i+1, mat_index+1, basis);
               basis -= i;
            }
         }
      }

      void randomly_number_simplices(){
         int n = all_bases.size();
         int a = 2321443;
         a %= n;
         if(a<0 || a>=n){
            std::cout << "Something went wrong." << std::endl;
         }
         int c = a;
         for(const auto& basis : all_bases){
            basis2int[basis] = c;
            int2basis[c] = basis;
            c += a;
            c %= n;
            if(c < 0){
               c += n;
            }
            if(c >= n){
               c -= n;
            }
         }
      }

   public:
      GKZ(const PointConfiguration& points, const Chirotope& c): pmMat(points.no(), points.rank()), indices(0), pts(points), chiro(c), pointNorms(points.no()){
         rk = points.no();
         dim = points.rank();
         read_points(points);
         populate_pointNorms();
         total_simplices = 0;
         populate_abs_determinants(0,0,basis_type());
         randomly_number_simplices();
         if(CommandlineOptions::verbose()){
            cerrFirstWorker << "GKZ done." << std::endl;
         }
      }

      bool isRegular(const TriangNode& t) const{
         RegularityCheck regcheck(pts, chiro, t);
         return regcheck.is_regular();
      }

      const pm::Map<basis_type, int>& get_basis2int() const{
         return basis2int;
      }
      
      const pm::Map<int, basis_type>& get_int2basis() const{
         return int2basis;
      }

      int size() const{
         return abs_determinants.size();
      }

      Scalar fine_coordinate(const gkz_vector_type& gkzvec) const{
         return (-1) * (pointNorms * gkzvec);
      }
      
      TriangVectortype triangulationAsVector(const TriangNode& triangulation) const{
         TriangVectortype result(total_simplices);
         for(TriangNode::const_iterator triangle = triangulation.begin();
               triangle != triangulation.end();
               ++triangle){
            if(!basis2int.exists(*triangle)){
               std::cout << "ARGH"<<std::endl;
            }
            result += basis2int[*triangle];
         }
         return result;
      }
      
      // Check whether we might have overflow
      void init_check(const TriangNode& triangulation) const{
         Scalar val(0);
         cerrFirstWorker << "Initial triangulation: ";
         for(TriangNode::const_iterator triangle = triangulation.begin();
               triangle != triangulation.end();
               ++triangle){
            cerrFirstWorker << *triangle << " ";
            val += abs_determinants[*triangle];
            if(val < 0){
               throw std::runtime_error("Normalized volume makes Scalar type overflow.");
            }
         }
         cerrFirstWorker << std::endl;
         cerrFirstWorker << "Normalized volume is: " << val << std::endl;
      }

      // Compute GKZ-vector of a triangulation
      Vectortype evaluate(const TriangNode& triangulation) const{
         Vectortype result(rk);
         Scalar val;
         for(TriangNode::const_iterator triangle = triangulation.begin();
               triangle != triangulation.end();
               ++triangle){
            val = abs_determinants[*triangle];
            for(basis_type::const_iterator entry = (*triangle).begin();
                  entry != (*triangle).end();
                  ++entry){
               result[*entry] += val;
            }
         }
         return result;
      }

      Vectortype evaluate_fine(const TriangNode& triangulation) const {
         Vectortype result = evaluate(triangulation);
         return fine_coordinate(result) | result;
      }

      std::string stringify(const Vectortype& v) const{
         std::string result = "[";
         for(int i=0; i<v.dim(); i++){
            result += std::to_string(v[i]);
            if(i<rk-1){
               result += ",";
            }
         }
         result += "]";
         return result;
      }

      // Return the GKZ vector of a triangulation as string for printing.
      std::string gkz_as_string(const TriangNode& triang) const{
         return stringify(evaluate(triang));
      }
      std::string gkz_as_string(const Flip& flip) const{
         return stringify(evaluate(flip));
      }

      // Compute GKZ-vector of a simplicial complex, multiply it by a factor.
      // The simplicial complex does not have to be a triangulation and most of
      // the time it will not be one. This method is needed to compute the
      // GKZ-vector of a flip.
      void evaluate(Vectortype& result, const FastSimplicialComplex& fsc, int factor) const{
         Scalar val;
         for(FastSimplicialComplex::const_iterator triangle = fsc.begin();
               triangle != fsc.end();
               ++triangle){
            val = abs_determinants[*triangle];
            for(basis_type::const_iterator entry = (*triangle).begin();
                  entry != (*triangle).end();
                  ++entry){
               result[*entry] += factor*val;
            }
         }
      }

      // Compute GKZ-vector of a flip, i.e. the difference of the GKZ-vectors
      // of the corresponding two triangulations.
      Vectortype evaluate(const Flip& flip) const{
         Vectortype result(rk);
         evaluate(result, flip.first, -1);
         evaluate(result, flip.second, 1);
         return result;
      }
      
      Vectortype evaluate_fine(const Flip& flip) const {
         Vectortype result = evaluate(flip);
         return fine_coordinate(result) | result;
      }

      int is_lex_positive(const gkz_vector_type& v) const{
         auto index = pm::indices(pm::ensure(v, (pm::cons<pm::end_sensitive, pm::pure_sparse>*)0));
         if(!index.begin().at_end()){
            Scalar check = v[index.front()];
            if(check > 0){
               return 1;
            } else if (check < 0){
               return -1;
            } else {
               cerrFirstWorker << "First entry was zero, this is a bug." << std::endl;
               return 0;
            }
         } else {
            return 0;
         }
      }
      
};


#endif
