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

///////////////////////////////////////////////////////////////////////////////
// RegularityCheckCache.hh
// 
// This file contains a cache for the regularity checks. When we walk down a
// flip and then check all the flips of the target triangulations for
// regularity, many of them will already have been previously checked, so a
// cache should be useful.
///////////////////////////////////////////////////////////////////////////////
#ifndef REGULARITYCHECKCACHE_HH
#define REGULARITYCHECKCACHE_HH
#include "mptopcom_typedefs.hh"
#include "TriangNode.hh"
#include "LRUCache.hh"
#include "GKZ.hh"


#include <polymake/Main.h>
#include <polymake/Matrix.h>
#include <polymake/Set.h>
#include <polymake/SparseVector.h>
#include <polymake/polytope/is_regular.h>

// Old polymake does not offer to choose LP solver
#ifdef NO_LP_SOLVE
#include <polymake/polytope/soplex_interface.h>
#else
#include <polymake/polytope/solve_LP.h>
#endif

namespace std {

template <>
struct iterator_traits< __sc_const_iterator<IntegerSet> > {
  typedef forward_iterator_tag iterator_category;
  typedef IntegerSet value_type;
  typedef void difference_type;
  typedef void pointer;
  typedef void reference;
};

// template <>
// struct iterator_traits< __is_const_iterator > {
//   typedef forward_iterator_tag iterator_category;
//   typedef size_type value_type;
//   typedef void difference_type;
//   typedef void pointer;
//   typedef void reference;
// };

}

namespace pm {

   template <>
      struct iterator_cross_const_helper<__sc_const_iterator<IntegerSet>, true> {
         typedef __sc_const_iterator<IntegerSet> iterator;
         typedef __sc_const_iterator<IntegerSet> const_iterator;
      };


   template <>
      struct can_initialize<IntegerSet, pm::Set<int>>
      : std::true_type {};


   template<>
      class conv<IntegerSet,pm::Set<int>> {
         public:
            typedef IntegerSet argument_type;
            typedef pm::Set<int> result_type;

            result_type operator() (const argument_type& x) const { 
               result_type res;
               for(auto i : x) 
                  res += i;
               return res;
            }
      };

}


class RegularityCheckCache {
   private:
      mutable LRUCache<TriangNode, bool> cache;
      const PointConfiguration& pts;
      const Chirotope& chiro;
      pm::Matrix<Scalar> pmPts;

      void init(const std::string& lp_client);

   public: 
      RegularityCheckCache(int size, const PointConfiguration& pts_in, const Chirotope& chiro_in, const std::string& lp_client):
         pts(pts_in), chiro(chiro_in), cache(size), pmPts(pts.no(), pts.rank())
         {
            init(lp_client);
            // Convert topcom pts to polymake pts.
            pm::Rational x;
            for(int i=0; i<pts.no(); i++){
               for(int j=0; j<pts.rank(); j++){
                  x.copy_from(pts[i][j].get_mpq_t());
                  pmPts(i,j) = pm::convert_to<Scalar>(x);
               }
            }
         };

      bool isRegular(const TriangNode& t) const{
         bool result;
         if(cache.exist(t)){
            result = cache.get(t);
         } else {
            result = isRegular_inner(t);
            cache.put(t, result);
         }
         return result;
      }

#ifndef NO_LP_SOLVE
      bool isRegular_inner(const TriangNode& t) const{

         // Convert triangulation to polymake datatype
         pm::Array<pm::Set<int>> pmMC(t.card());
         int i = 0;
         for(auto s : t) {
            pmMC[i] = pm::convert_to<pm::Set<int>>(s);
            i++;
         }
         // Make secondary cone ineqs
         std::pair<const pm::Matrix<Scalar>, const pm::Matrix<Scalar>> SC = polymake::polytope::secondary_cone_ineq(pmPts, pmMC, pm::perl::OptionSet("lift_face_to_zero", 0));
         bool result = is_regular_lp(SC);
         return result;
      }

      bool is_regular_lp(std::pair<const pm::Matrix<Scalar>, const pm::Matrix<Scalar>> SC) const {
         const int n_vertices=pts.no();
         pm::Matrix<pm::Rational> H(SC.first);
         pm::Matrix<pm::Rational> E(SC.second);
         const pm::Vector<pm::Rational> Obj(Scalar(0) | pm::ones_vector<Scalar>(n_vertices));
         bool maximize = false;
         E = pm::zero_matrix<pm::Rational>(E.rows(),1) | E;
         H = ((pm::Rational(-1))*(pm::ones_matrix<pm::Rational>(H.rows(),1))) | H;

         const auto S = polymake::polytope::solve_LP(H, E, Obj, maximize);
         switch (S.status) {
            case polymake::polytope::LP_status::valid:
            case polymake::polytope::LP_status::unbounded:
               return true;
         }
         return false;
      }

#else
# ifndef WITH_SOPLEX
      // Old method without soplex
      bool isRegular_inner(const TriangNode& t) const{
         RegularityCheck regcheck(pts, chiro, t);
         return regcheck.is_regular();
      }
# else
      bool isRegular_inner(const TriangNode& t) const{

         // Convert triangulation to polymake datatype
         pm::Array<pm::Set<int>> pmMC(t.card());
         int i = 0;
         for(auto s : t) {
            pmMC[i] = pm::convert_to<pm::Set<int>>(s);
            i++;
         }
         // Make secondary cone ineqs
         std::pair<const pm::Matrix<Scalar>, const pm::Matrix<Scalar>> SC = polymake::polytope::secondary_cone_ineq(pmPts, pmMC, pm::perl::OptionSet("lift_face_to_zero", 0));
         bool result = is_regular_soplex(SC);
         return result;
      }

      bool is_regular_soplex(std::pair<const pm::Matrix<Scalar>, const pm::Matrix<Scalar>> SC) const {
         const int n_vertices=pts.no();
         typedef polymake::polytope::soplex_interface::solver Solver;
         pm::Matrix<pm::Rational> H(SC.first);
         pm::Matrix<pm::Rational> E(SC.second);
         const pm::Vector<pm::Rational> Obj(Scalar(0) | pm::ones_vector<Scalar>(n_vertices));
         bool maximize = false;
         E = pm::zero_matrix<pm::Rational>(E.rows(),1) | E;
         H = ((pm::Rational(-1))*(pm::ones_matrix<pm::Rational>(H.rows(),1))) | H;

         try
         {
            Solver solver;
            typename Solver::lp_solution S = solver.solve_lp(H, E, Obj, maximize);
            return true;
         }
         catch (polymake::infeasible)
         {
            // std::cerr << "Infeasible." << std::endl;
            return false;
         }
         catch (polymake::polytope::unbounded)
         {
            // std::cerr << "Unbounded." << std::endl;
            return true;
         }
         return false;
      }
# endif
#endif

};


#endif

