src/AlgebraicCore/ExternalLibs-Normaliz.C

//   Copyright (c)  2010-2015 Anna M. Bigatti, Christof Soeger
//   Authors:  2010-2015 Anna M. Bigatti, Christof Soeger

//   This file is part of the source of CoCoALib, the CoCoA Library.

//   CoCoALib 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 of the License, or
//   (at your option) any later version.

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

//   You should have received a copy of the GNU General Public License
//   along with CoCoALib.  If not, see <http://www.gnu.org/licenses/>.


// Source code for Normaliz integration

#include "CoCoA/ExternalLibs-Normaliz.H"

#ifdef CoCoA_WITH_NORMALIZ

#include "CoCoA/BigIntOps.H"
#include "CoCoA/BigRat.H"
#include "CoCoA/TmpPPVector.H"
#include "CoCoA/convert.H"
#include "CoCoA/error.H"
#include "CoCoA/matrix.H"
#include "CoCoA/SparsePolyRing.H"
#include "CoCoA/symbol.H"
#include "CoCoA/QuasiPoly.H"

#include "CoCoA/VectorOps.H"  // just for debugging

//#include "libnormaliz/libnormaliz.h"
using libnormaliz::Cone;
using libnormaliz::ConeProperties;


// #include <list>
// using std::list;
#include <map>
using std::map;
#include <vector>
using std::vector;
#include <ostream>
using std::ostream;
using std::endl;
#include <sstream>  // for ErrorMessage


namespace CoCoA
{

  namespace Normaliz
  {

    libnormaliz::InputType ToInputType(const std::string& TypeString)
    {
      return libnormaliz::to_type(TypeString);
    }

    class ConeImpl: protected IntrusiveReferenceCount
    {
      friend class SmartPtrIRC<const ConeImpl>; // Morally "friend Cone", so it can alter reference count.

      public:
        ConeImpl(libnormaliz::Type::InputType InputType, const std::vector<std::vector<BigInt> >& v);
        ConeImpl(const std::map< libnormaliz::InputType, std::vector<std::vector<BigInt> > >& m);
        ConeImpl(libnormaliz::Type::InputType InputType, const PPVector& v);
        friend void SetVerbosityLevel(cone& c, long n); // 0 off/ >1 on
        friend long VerbosityLevel(cone& c);
        friend std::vector<std::vector<BigInt> > HilbertBasis(const cone& c);
        friend std::vector<std::vector<BigInt> > ModuleGenerators(const cone& c);
        friend std::vector<std::vector<BigInt> > TriangulationGenerators(const cone& c);
        friend std::vector<std::vector<BigInt> > ExtremeRays(const cone& c);
        friend std::vector<std::vector<BigInt> > VerticesOfPolyhedron(const cone& c);
        friend std::vector<std::vector<BigInt> > Deg1Elements(const cone& c);
        friend std::vector<std::vector<BigInt> > GeneratorsOfToricRing(const cone& c);
        friend std::vector<std::vector<BigInt> > SupportHyperplanes(const cone& c);
        friend std::vector<std::vector<BigInt> > Equations(const cone& c);
        friend std::vector<std::vector<BigInt> > Congruences(const cone& c);
        friend std::vector<std::vector<BigInt> > ExcludedFaces(const cone& c);
        //friend std::vector<BigInt> HVector(const cone& c);
        friend HPSeries HilbertSeries(const cone& c);
        friend RingElem HilbertPoly(const cone& c);
        friend QuasiPoly HilbertQuasiPoly(const cone& c);
        friend BigRat multiplicity(const cone& c);
        friend std::vector<BigRat> grading(const cone& c);
        friend bool IsPointed(const cone& c);
        friend bool IsInhomogeneous(const cone& c);
        friend bool IsIntegrallyClosed(const cone& c);
        friend bool IsDeg1HilbertBasis(const cone& c);
        friend long EmbeddingDim(const cone& c);
        friend long rank(const cone& c);
        friend long RecessionRank(const cone& c);
        friend long AffineDim(const cone& c);
        friend long ModuleRank(const cone& c);
        friend std::vector<BigInt> dehomogenization(const cone& c);

      public:
        void myComputation(const libnormaliz::ConeProperties& CPs) const;
        void myComputation(libnormaliz::ConeProperty::Enum CP) const;
        void myComputation(libnormaliz::ConeProperty::Enum CP1, libnormaliz::ConeProperty::Enum CP2) const;
        void myComputation() const;   //default: compute everything possible
        bool isComputed(libnormaliz::ConeProperty::Enum CP) const;
        bool isComputed(libnormaliz::ConeProperty::Enum CP1, libnormaliz::ConeProperty::Enum CP2) const;

      private:
        mutable libnormaliz::Cone<mpz_class> myConeMPZ; // need "mutable" to cater for Normaliz interpretation of "const" (namely, they take it to mean that the repr is const, and not just the object being represented).
    };


    // implementation of cone

    cone::cone(libnormaliz::Type::InputType InputType, const std::vector<std::vector<BigInt> >& v) : mySmartPtr(new ConeImpl(InputType, v)) {}

    cone::cone(const std::map< libnormaliz::InputType, std::vector<std::vector<BigInt> > >& m) : mySmartPtr(new ConeImpl(m)) {}

    cone::cone(libnormaliz::Type::InputType InputType, ConstMatrixView M) : mySmartPtr(new ConeImpl(InputType, MatrixToVecVecBigInt(M))) {}

    cone::cone(const ConeImpl* ptr): mySmartPtr(ptr) {}
    cone::cone(const cone& c) : mySmartPtr(c.mySmartPtr) {}
    cone::~cone() {}


    // printing
    ostream& operator<< (ostream& out, const cone& C)
    {
      if (!out) return out;  // short-cut for bad ostreams
      using namespace libnormaliz;
      out << "cone(EmbeddingDim = " << EmbeddingDim(C);
      if (C.isComputed(ConeProperty::ExtremeRays)) {
        out << ", rank = " << rank(C);
        out << ", pointed = " << IsPointed(C);
      }
      if (C.isComputed(ConeProperty::Generators))
        out << ", NumGenerators = " << TriangulationGenerators(C).size();
      if (C.isComputed(ConeProperty::SupportHyperplanes))
        out << ", NumSupportHyperplanes = " << SupportHyperplanes(C).size();
      out << ")";
      return out;
    }

    // letVerbose(Default) sets the verbosity, and returns the old value
    void SetVerbosityLevel(long n) {
      libnormaliz::setVerboseDefault(n!=0);
    }

    long VerbosityLevel() {
      const bool b = libnormaliz::setVerboseDefault(0);
      libnormaliz::setVerboseDefault(b);
      return b;
    }

    void SetVerbosityLevel(cone& c, long n) {
      c->myConeMPZ.setVerbose(n!=0);
    }

    long VerbosityLevel(cone& c) {
      const bool b = c->myConeMPZ.setVerbose(0);
      c->myConeMPZ.setVerbose(b);
      return b;
    }

    // computations
    void cone::myComputation(const libnormaliz::ConeProperties& CPs) const { mySmartPtr->myComputation(CPs); }
    void cone::myComputation(libnormaliz::ConeProperty::Enum CP) const { mySmartPtr->myComputation(CP); }
    void cone::myComputation(libnormaliz::ConeProperty::Enum CP1, libnormaliz::ConeProperty::Enum CP2) const { mySmartPtr->myComputation(CP1, CP2); }
    void cone::myComputation() const { mySmartPtr->myComputation(); }
    bool cone::isComputed(libnormaliz::ConeProperty::Enum CP) const { return mySmartPtr->isComputed(CP); }
    bool cone::isComputed(libnormaliz::ConeProperty::Enum CP1, libnormaliz::ConeProperty::Enum CP2) const { return mySmartPtr->isComputed(CP1, CP2); }
    const ConeImpl* cone::operator->() const { return mySmartPtr.operator->(); }

    namespace  // conversion functions, implementation at the end of the file
    {
//      std::vector<BigInt> LongVToBigIntV(const std::vector<long>& VIn);
      std::vector<BigInt> MPZ_classVToBigIntV(const std::vector<mpz_class>& VIn);
//      std::vector<BigRat> LongVToBigRatV(const std::vector<long>& VIn, const BigInt& denom);
      std::vector<BigRat> MPZ_classVToBigRatV(const std::vector<mpz_class>& VIn, const BigInt& denom);
//      std::vector<long> BigIntVToLongV(const std::vector<BigInt>& VIn);
      std::vector<mpz_class> BigIntVToMPZ_classV(const std::vector<BigInt>& VIn);
//      void ConvertFromNormaliz(std::vector<std::vector<BigInt> >& v, const std::vector<std::vector<long> >& l);
      void ConvertFromNormaliz(std::vector<std::vector<BigInt> >& v, const std::vector<std::vector<mpz_class> >& l);
//      std::vector<std::vector<long> >  ReturnVecVecLong(const std::vector<std::vector<BigInt> >& v);
//      map< libnormaliz::InputType, vector<vector<long> > >  ReturnMapVecVecLong(const map< libnormaliz::InputType, vector<vector<BigInt> > >& m);
      std::vector<std::vector<mpz_class> >  ReturnVecVecMPZ_class(const std::vector<std::vector<BigInt> >& v);
      map< libnormaliz::InputType, vector<vector<mpz_class> > >  ReturnMapVecVecMPZ_class(const map< libnormaliz::InputType, vector<vector<BigInt> > >& m);
      std::vector<std::vector<BigInt> > PPVectorToVecVecBigInt(const PPVector& ppv);
      void VecVecBigIntToPPVector(PPVector& ppv, const std::vector<std::vector<BigInt> >& M);
    } // end of anonymous namespace


    // implementation of ConeImpl
    ConeImpl::ConeImpl(libnormaliz::Type::InputType InputType, const std::vector<std::vector<BigInt> >& v):
          myConeMPZ(InputType, ReturnVecVecMPZ_class(v))
    {
    }

    ConeImpl::ConeImpl(const std::map< libnormaliz::InputType, std::vector<std::vector<BigInt> > >& m):
      myConeMPZ(ReturnMapVecVecMPZ_class(m))
    {
    }

    void ConeImpl::myComputation(const libnormaliz::ConeProperties& CPs) const
    {
        libnormaliz::ConeProperties missing = myConeMPZ.compute(CPs);
        if (missing.any())
        {
          std::ostringstream os;
          os << missing;
          CoCoA_THROW_ERROR("Normaliz cannot compute "+ os.str(), "myComputation MPZ");
        }
    }

    void ConeImpl::myComputation(libnormaliz::ConeProperty::Enum CP) const
    {
      myComputation(ConeProperties(CP));
    }

    void ConeImpl::myComputation(libnormaliz::ConeProperty::Enum CP1,
                                 libnormaliz::ConeProperty::Enum CP2) const
    {
      myComputation(ConeProperties(CP1, CP2));
    }

    void ConeImpl::myComputation() const
    {
      myComputation(ConeProperties(libnormaliz::ConeProperty::DefaultMode));
    }

    bool ConeImpl::isComputed(libnormaliz::ConeProperty::Enum CP) const
    {
      return myConeMPZ.isComputed(CP);
    }

    bool ConeImpl::isComputed(libnormaliz::ConeProperty::Enum CP1,
                              libnormaliz::ConeProperty::Enum CP2) const
    {
        return isComputed(CP1) && isComputed(CP2);
    }

    // friend functions which are the interface to the user
    std::vector<std::vector<BigInt> > HilbertBasis(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::HilbertBasis);
      ConvertFromNormaliz(v, c->myConeMPZ.getHilbertBasis());
      return v;
    }

    std::vector<std::vector<BigInt> > ModuleGenerators(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::ModuleGenerators);
      ConvertFromNormaliz(v, c->myConeMPZ.getModuleGenerators());
      return v;
    }

    std::vector<std::vector<BigInt> > Deg1Elements(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::Deg1Elements);
      ConvertFromNormaliz(v, c->myConeMPZ.getDeg1Elements());
      return v;
    }

    std::vector<std::vector<BigInt> > TriangulationGenerators(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::Generators);
      ConvertFromNormaliz(v, c->myConeMPZ.getTriangulationGenerators());
      return v;
    }

    std::vector<std::vector<BigInt> > ExtremeRays(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::ExtremeRays);
      ConvertFromNormaliz(v, c->myConeMPZ.getExtremeRays());
      return v;
    }

    std::vector<std::vector<BigInt> > VerticesOfPolyhedron(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::VerticesOfPolyhedron);
      ConvertFromNormaliz(v, c->myConeMPZ.getVerticesOfPolyhedron());
      return v;
    }

    std::vector<std::vector<BigInt> > GeneratorsOfToricRing(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::OriginalMonoidGenerators);
      ConvertFromNormaliz(v, c->myConeMPZ.getOriginalMonoidGenerators());
      return v;
    }


    // The following constraints depend all on ConeProperty::SupportHyperplanes
    std::vector<std::vector<BigInt> > SupportHyperplanes(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::SupportHyperplanes);
      ConvertFromNormaliz(v, c->myConeMPZ.getSupportHyperplanes());
      return v;
    }

    std::vector<std::vector<BigInt> > Equations(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::Sublattice);
      ConvertFromNormaliz(v, c->myConeMPZ.getSublattice().getEquations());
      return v;
    }

    std::vector<std::vector<BigInt> > Congruences(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::Sublattice);
      ConvertFromNormaliz(v, c->myConeMPZ.getSublattice().getCongruences());
      return v;
    }

    // excluded faces are only available if they were explicitly given
    std::vector<std::vector<BigInt> > ExcludedFaces(const cone& c)
    {
      std::vector<std::vector<BigInt> > v;
      c->myComputation(libnormaliz::ConeProperty::ExcludedFaces);
      ConvertFromNormaliz(v, c->myConeMPZ.getExcludedFaces());
      return v;
    }

    HPSeries HilbertSeries(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::HilbertSeries);
      const libnormaliz::HilbertSeries& HS = c->myConeMPZ.getHilbertSeries();

      return HPSeries(MPZ_classVToBigIntV(HS.getNum()),
                      libnormaliz::to_vector(HS.getDenom()), HS.getShift());
    }

    RingElem HilbertPoly(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::HilbertSeries);
      const libnormaliz::HilbertSeries& HS = c->myConeMPZ.getHilbertSeries();
      if (HS.getPeriod() != 1)
      {
        CoCoA_THROW_ERROR("Hilbert function is a quasi-polynomial of period > 1.  This function works for regular polynomials only.", "HilbertPoly");
      }
      const vector<BigRat> coeffs = MPZ_classVToBigRatV( HS.getHilbertQuasiPolynomial()[0],
                                                         BigIntFromMPZ(HS.getHilbertQuasiPolynomialDenom().get_mpz_t()));
      PolyRing QQt = RingQQt(1);
      const RingElem t = indet(QQt,0);
      RingElem tpower = one(QQt);
      RingElem hpoly(QQt);
      for (long i=0; i<len(coeffs); ++i)
      {
          hpoly += coeffs[i] * tpower;   // NB tpower = t^i
          tpower *= t;
      }

      return hpoly;
    }

    QuasiPoly HilbertQuasiPoly(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::HilbertSeries);
      const libnormaliz::HilbertSeries& HS = c->myConeMPZ.getHilbertSeries();
      const long period = HS.getPeriod();
      if (period < 1)
      {
        CoCoA_THROW_ERROR("Hilbert function not computed.", "HilbertQuasiPoly");
      }

      const PolyRing QQt = RingQQt(1);
      const RingElem t = indet(QQt,0);
      RingElem tpower;
      vector<BigRat> coeffs;
      vector<RingElem> qp = vector<RingElem>(period, RingElem(QQt));
      for (long j=0; j<period; ++j)
      {
        coeffs = MPZ_classVToBigRatV( HS.getHilbertQuasiPolynomial()[j],
                                      BigIntFromMPZ(HS.getHilbertQuasiPolynomialDenom().get_mpz_t()));
        tpower = one(QQt);
        for (long i=0; i<len(coeffs); ++i)
        {
            qp[j]  += coeffs[i] * tpower;   // NB tpower = t^i
            tpower *= t;
        }
      }

      return QuasiPoly(qp);
    }

    BigRat multiplicity(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::Multiplicity);
      return BigRatFromMPQ(c->myConeMPZ.getMultiplicity().get_mpq_t());
    }


    bool IsPointed(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::ExtremeRays);
      return c->myConeMPZ.isPointed();
    }

    bool IsInhomogeneous(const cone& c)
    {
      // is always known
      return c->myConeMPZ.isInhomogeneous();
    }

    bool IsIntegrallyClosed(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::HilbertBasis);
      return c->myConeMPZ.isIntegrallyClosed();
    }

    bool IsDeg1HilbertBasis(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::HilbertBasis, libnormaliz::ConeProperty::Grading);
      return c->myConeMPZ.isDeg1HilbertBasis();
    }


    // dimension and rank invariants
    long EmbeddingDim(const cone& c)
    {
      // is always known
      return c->myConeMPZ.getEmbeddingDim();
    }

    long rank(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::ExtremeRays);
      return c->myConeMPZ.getRank();
    }

    std::vector<BigRat> grading(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::Grading);
      return MPZ_classVToBigRatV(c->myConeMPZ.getGrading(), BigIntFromMPZ(c->myConeMPZ.getGradingDenom().get_mpz_t()));
    }

    // only for inhomogeneous case:
    long RecessionRank(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::RecessionRank);
      return c->myConeMPZ.getRecessionRank();
    }

    long AffineDim(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::ExtremeRays);
      return c->myConeMPZ.getRank();
    }

    long ModuleRank(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::ModuleRank);
      return c->myConeMPZ.getModuleRank();
    }

    std::vector<BigInt> dehomogenization(const cone& c)
    {
      c->myComputation(libnormaliz::ConeProperty::Dehomogenization);
      return MPZ_classVToBigIntV(c->myConeMPZ.getDehomogenization());
    }

    // not so important in this library
    // size_t getTriangulationSize() const;
    // Integer getTriangulationDetSum() const;


/**************  applications to monomials  **************/

    namespace // anonymous for file local fns
    {
      PPVector HilbertBasis_PPVector(libnormaliz::InputType t, const PPVector& ppv)
      {
        cone c (t, PPVectorToVecVecBigInt(ppv));
        std::vector<std::vector<BigInt> > res = HilbertBasis(c);
        //create an emptyPPVector with the same PPMonoid
        PPVector ppv_res(PPM(ppv), DMR(ppv));
        VecVecBigIntToPPVector(ppv_res,res);
        return ppv_res;
      }

      //only for modes where normaliz adds a component, like polytope or rees_algebra
      PPVector HilbertBasis_PPVector(libnormaliz::Type::InputType t, const PPVector& ppv, long sym_pos)
      {
        const long dim = NumIndets(PPM(ppv));
        if (sym_pos < 0 || sym_pos >= dim)
          CoCoA_THROW_ERROR("Symbol position not valid", "HilbertBasis_PPVector");

        vector< vector<BigInt> > vv_input = PPVectorToVecVecBigInt(ppv);
        if (sym_pos != dim-1)
          CoCoA_THROW_ERROR("Symbol needs to be the last symbol of the PPMonoid (other implementations are missing)",
                      "HilbertBasis_PPVector (with extra symbol)");
        for (long i=0; i < len(vv_input); ++i)
        {
          if (vv_input[i][sym_pos] != 0)
            CoCoA_THROW_ERROR("Symbol must not be used in the power products", "HilbertBasis_PPVector (with extra symbol)");
          vv_input[i].erase(vv_input[i].begin()+sym_pos);
        }
        if (t != libnormaliz::Type::polytope && t != libnormaliz::Type::rees_algebra)
          CoCoA_THROW_ERROR("Invalid InputType for this method", "HilbertBasis_PPVector (with extra symbol)");


        cone c (t, vv_input);
        vector<vector<BigInt> > res = HilbertBasis(c);
        if (sym_pos != dim-1)
        {
          for (long i=0; i < dim; ++i)
          {
            swap(res[i][dim],res[i][sym_pos]);
          }
        }
        //create an emptyPPVector with the same PPMonoid
        PPVector ppv_res(PPM(ppv), DMR(ppv));
        VecVecBigIntToPPVector(ppv_res,res);

        return ppv_res;
      }

      PPVector HilbertBasis_PPVector(libnormaliz::InputType t, const std::vector<std::vector<BigInt> >& Mat, PPMonoid ppm)
      {
        cone c (t, Mat);
        std::vector<std::vector<BigInt> > res = HilbertBasis(c);
        //create an empty PPVector with the same PPMonoid
        PPVector ppv_res(ppm,NewDivMaskNull());
        VecVecBigIntToPPVector(ppv_res,res);
        return ppv_res;
      }

      PPVector HilbertBasis_PPVector(const cone& C, PPMonoid ppm)
      {
        std::vector<std::vector<BigInt> > res = HilbertBasis(C);
        //create an empty PPVector with the same PPMonoid
        PPVector ppv_res(ppm,NewDivMaskNull());
        VecVecBigIntToPPVector(ppv_res,res);
        return ppv_res;
      }
    } // end anonymous namespace

    PPVector NormalToricRing(const PPVector& ppv)
    {
      return HilbertBasis_PPVector(libnormaliz::Type::normalization, ppv);
    }

    PPVector IntClosureToricRing(const PPVector& ppv)
    {
      return HilbertBasis_PPVector(libnormaliz::Type::integral_closure, ppv);
    }

    PPVector IntClosureMonIdeal(const PPVector& ppv)
    {
      cone c (libnormaliz::Type::rees_algebra, PPVectorToVecVecBigInt(ppv));
//      const vector<vector<BigInt> > res = HilbertBasis(c);
      vector<vector<BigInt> > res1;
//      for (vector<vector<BigInt> >::const_iterator it=res.begin(); it!=res.end(); ++it)
      for (const auto& v: HilbertBasis(c))
      {
        //BigInt tmp = it->back();
        if (v.back() == 1 ) {
          res1.push_back(v);
          res1.back().pop_back();
        }
      }

      //create an empty PPVector with the same PPMonoid
      PPVector ppv_res1(PPM(ppv), DMR(ppv));
      VecVecBigIntToPPVector(ppv_res1,res1);
      return ppv_res1;
    }

    /* If you want also the whole normalization of the rees algebra then the
     * ring must have an unused symbol that we can use for homogenization.
     * You have to specify its position.
     */
    PPVector IntClosureMonIdeal(const PPVector& ppv, long sym_pos)
    {
      return HilbertBasis_PPVector(libnormaliz::Type::rees_algebra, ppv, sym_pos);
    }

    PPVector EhrhartRing(const PPVector& ppv, long sym_pos)
    {
      return HilbertBasis_PPVector(libnormaliz::Type::polytope, ppv, sym_pos);
    }

/**************  torus invariants and valuation rings  **************/
    PPVector TorusInvariants(const vector< vector<BigInt> >& T, const PPMonoid& ppm)
    {
      if (T.empty())
      {
        CoCoA_THROW_ERROR("Matrix should be non-empty", "TorusInvariants");
      }
      if (NumIndets(ppm) != len(T[0]))
      {
        CoCoA_THROW_ERROR("Number of columns in matrix does not match number of variables in ring", "TorusInvariants");
      }
      return HilbertBasis_PPVector(libnormaliz::Type::equations, T, ppm);
    }

    PPVector FiniteDiagInvariants(const vector< vector<BigInt> >& Cong, const PPMonoid& ppm)
    {
      if (Cong.empty())
      {
        CoCoA_THROW_ERROR("Matrix should be non-empty", "FiniteDiagInvariants");
      }
      if (NumIndets(ppm) != len(Cong[0])-1)
      {
        CoCoA_THROW_ERROR("Number of columns in matrix -1 does not match number of variables in ring", "FiniteDiagInvariants");
      }
      return HilbertBasis_PPVector(libnormaliz::Type::congruences, Cong, ppm);
    }

    PPVector DiagInvariants(const vector< vector<BigInt> >& T, const vector< vector<BigInt> >& Cong, const PPMonoid& ppm)
    {
      if (T.empty() && Cong.empty())
      {
        CoCoA_THROW_ERROR("At least one Matrix should be non-empty", "DiagInvariants");
      }
      if (T.empty())
        return FiniteDiagInvariants(Cong, ppm);
      if (Cong.empty())
        return TorusInvariants(T, ppm);

      if (NumIndets(ppm) != len(T[0]))
      {
        CoCoA_THROW_ERROR("Number of columns in matrix does not match number of variables in ring", "DiagInvariants");
      }
      if (NumIndets(ppm) != len(Cong[0])-1)
      {
        CoCoA_THROW_ERROR("Number of columns in matrix -1 does not match number of variables in ring", "DiagInvariants");
      }
      std::map< libnormaliz::InputType, std::vector<std::vector<BigInt> > > cone_input;
      cone_input[libnormaliz::Type::equations] = T;
      cone_input[libnormaliz::Type::congruences] = Cong;
      cone C(cone_input);
      return HilbertBasis_PPVector(C, ppm);
    }


    PPVector IntersectionValRings (const vector< vector<BigInt> >& V, const PPMonoid& ppm)
    {
      const long dim = len(V[0]);
      std::map< libnormaliz::InputType, std::vector<std::vector<BigInt> > > cone_input;
      cone_input[libnormaliz::Type::inequalities] = V;
      const vector< vector<BigInt> > positive_signs = vector< vector<BigInt> >(1,vector<BigInt>(dim,BigInt(1)));
      cone_input[libnormaliz::Type::signs] = positive_signs;
      cone C(cone_input);
      return HilbertBasis_PPVector(C, ppm);
    }


/********************************************************
 ***               conversion functions               ***
 ********************************************************/

    namespace  // conversion functions
    {

// Compiler says this is never used
      // std::vector<BigInt> LongVToBigIntV(const std::vector<long>& VIn)
      // {
      //   std::vector<BigInt> v;
      //   v.reserve(VIn.size());
      //   for (vector<long>::const_iterator it=VIn.begin(); it!=VIn.end(); ++it)
      //     v.push_back(BigInt(*it));
      //   return v;
      // }

      std::vector<BigInt> MPZ_classVToBigIntV(const std::vector<mpz_class>& VIn)
      {
        std::vector<BigInt> v; v.reserve(VIn.size());
//        for (vector<mpz_class>::const_iterator it=VIn.begin(); it!=VIn.end(); ++it)
        for (const mpz_class& N: VIn)
          v.push_back(BigIntFromMPZ(N.get_mpz_t()));
        return v;
      }

// Compiler says this is never used
      // std::vector<BigRat> LongVToBigRatV(const std::vector<long>& VIn, const BigInt& denom)
      // {
      //   std::vector<BigRat> v;
      //   v.reserve(VIn.size());
      //   for (vector<long>::const_iterator it=VIn.begin(); it!=VIn.end(); ++it)
      //     v.push_back(BigRat(BigInt(*it),denom));
      //   return v;
      // }

      std::vector<BigRat> MPZ_classVToBigRatV(const std::vector<mpz_class>& VIn, const BigInt& denom)
      {
        std::vector<BigRat> v;
        v.reserve(VIn.size());
//        for (vector<mpz_class>::const_iterator it=VIn.begin(); it!=VIn.end(); ++it)
        for(const mpz_class& N: VIn)
          v.push_back(BigRat(BigIntFromMPZ(N.get_mpz_t()), denom));
        return v;
      }


// Compiler says this is never used
      // std::vector<long> BigIntVToLongV(const std::vector<BigInt>& VIn)
      // {
      //   std::vector<long> v;
      //   long n;
      //   v.reserve(VIn.size());
      //   for (vector<BigInt>::const_iterator it=VIn.begin(); it!=VIn.end(); ++it)
      //     if (IsConvertible(n, *it))  v.push_back(n);
      //     else CoCoA_THROW_ERROR(ERR::BadConvert, "normaliz: BigIntVToLongV");
      //   return v;
      // }


      std::vector<mpz_class> BigIntVToMPZ_classV(const std::vector<BigInt>& VIn)
      {
        std::vector<mpz_class> v;
        v.reserve(VIn.size());
//        for (vector<BigInt>::const_iterator it=VIn.begin(); it!=VIn.end(); ++it)
        for (const BigInt& N: VIn)
          v.push_back(mpz_class(mpzref(N)));
        return v;
      }


// Compiler says this is never used
      // void ConvertFromNormaliz(std::vector<std::vector<BigInt> >& v, const std::vector<std::vector<long> >& l)
      // {
      //   v.clear();      // not exception safe
      //   v.reserve(l.size());
      //   //      transform(l.begin(), l.end(), v.begin(), LongVToBigIntV);
      //   for (vector<vector<long> >::const_iterator it=l.begin(); it!=l.end(); ++it)
      //     v.push_back(LongVToBigIntV(*it));
      // }


      void ConvertFromNormaliz(std::vector<std::vector<BigInt> >& v, const std::vector<std::vector<mpz_class> >& L)
      {
        v.clear();      // not exception safe
        v.reserve(L.size());
        //      transform(L.begin(), L.end(), v.begin(), LongVToBigIntV);
//        for (vector<vector<mpz_class> >::const_iterator it=L.begin(); it!=L.end(); ++it)
        for (const auto& vec: L)
          v.push_back(MPZ_classVToBigIntV(vec));
      }


// Compiler says this is never used
      // std::vector<std::vector<long> >  ReturnVecVecLong(const std::vector<std::vector<BigInt> >& v)
      // {
      //   std::vector<std::vector<long> > l;
      //   //      transform(v.begin(), v.end(), l.begin(), BigIntVToLongV);
      //   try
      //   {
      //       for (vector<vector<BigInt> >::const_iterator it=v.begin(); it!=v.end(); ++it)
      //         l.push_back(BigIntVToLongV(*it));
      //   } catch (...) {l.clear();}
      //   return l;
      // }


// Compiler says this is never used
      // map< libnormaliz::InputType, vector<vector<long> > >  ReturnMapVecVecLong(const map< libnormaliz::InputType, vector<vector<BigInt> > >& m)
      // {
      //   map< libnormaliz::InputType, vector<vector<long> > > ret;
      //   for (map<libnormaliz::InputType, vector<vector<BigInt> > >::const_iterator it=m.begin(); it!=m.end(); ++it)
      //     ret.insert(make_pair((*it).first, ReturnVecVecLong((*it).second)));
      //   return ret;
      // }


      std::vector<std::vector<mpz_class> >  ReturnVecVecMPZ_class(const std::vector<std::vector<BigInt> >& vlist)
      {
        std::vector<std::vector<mpz_class> > L; L.reserve(vlist.size());
        //      transform(v.begin(), v.end(), L.begin(), BigIntVToLongV);
//        for (vector<vector<BigInt> >::const_iterator it=v.begin(); it!=v.end(); ++it)
        for (const auto& vec: vlist)
          L.push_back(BigIntVToMPZ_classV(vec));
        return L;
      }


      map< libnormaliz::InputType, vector<vector<mpz_class> > >  ReturnMapVecVecMPZ_class(const map< libnormaliz::InputType, vector<vector<BigInt> > >& mlist)
      {
        map< libnormaliz::InputType, vector<vector<mpz_class> > > ret;
//        for (map<libnormaliz::InputType, vector<vector<BigInt> > >::const_iterator it=m.begin(); it!=m.end(); ++it)
        for (const auto& m: mlist)
          ret.insert(make_pair(m.first, ReturnVecVecMPZ_class(m.second)));
        return ret;
      }


      std::vector<std::vector<BigInt> > PPVectorToVecVecBigInt(const PPVector& ppv)
      {
        vector<BigInt> tmp;
        const long n =  len(ppv);
        vector<vector<BigInt> > v(n);
        for (long i=0; i<n; ++i)
        {
          BigExponents(tmp,PP(ppv[i]));
          v[i]=tmp;
        }
        return v;
      }

      void VecVecBigIntToPPVector(PPVector& ppv, const std::vector<std::vector<BigInt> >& M)
      {
        PPMonoid ppm = PPM(ppv);
        const long n =  len(M);
        for (long i=0; i < n; ++i)
        {
          ppv.myPushBack(PPMonoidElem(ppm, M[i]));
        }
      }

    } // end of anonymous namespace

    std::vector<std::vector<BigInt> > MatrixToVecVecBigInt(ConstMatrixView M)
    {
      const ErrorInfo ErrMesg("Matrix entries must be integer", "MatrixToVecVecBigInt (cone ctor)");

      vector<vector<BigInt> > v;
      for (long i=0; i<NumRows(M); ++i)
      {
        v.push_back(vector<BigInt>());
        for (long j=0; j<NumCols(M); ++j)
          v[i].push_back(ConvertTo<BigInt>(M(i,j), ErrMesg));
      }
      return v;
    }

    PPVector MonomialsToPPV(const std::vector<RingElem>& v)
    {
      if(!AreMonomials(v)) {
        CoCoA_THROW_ERROR("Expected list of monomials","MonomialsToPPV");
      }
      if(v.empty()) {
        CoCoA_THROW_ERROR("List of monomials has to be non-empty","MonomialsToPPV");
      }
      //convert it to a PPVector
      PPVector ppv(PPM(owner(v[0])), NewDivMaskNull());
      convert(ppv,v);
      return ppv;
    }


  } // namespace Normaliz
} // namespace CoCoA

// RCS header/log in the next few lines
// $Header: /Volumes/Home_1/cocoa/cvs-repository/CoCoALib-0.99/src/AlgebraicCore/ExternalLibs-Normaliz.C,v 1.46 2020/10/02 18:46:41 abbott Exp $
// $Log: ExternalLibs-Normaliz.C,v $
// Revision 1.46  2020/10/02 18:46:41  abbott
// Summary: Cleaned (actually removed) include directive; renamed Generators to TriangulationGenerators
//
// Revision 1.45  2020/06/17 15:49:22  abbott
// Summary: Changed CoCoA_ERROR into CoCoA_THROW_ERROR
//
// Revision 1.44  2020/02/18 11:27:04  abbott
// Summary: redmine 1346: new for loop syntax
//
// Revision 1.43  2019/11/14 17:56:07  abbott
// Summary: Removed commented out use of clog
//
// Revision 1.42  2018/05/22 14:16:39  abbott
// Summary: Split BigRat into BigRat (class defn + ctors) and BigRatOps
//
// Revision 1.41  2018/05/18 12:13:36  bigatti
// -- renamed IntOperations --> BigIntOps
//
// Revision 1.40  2018/05/17 15:44:36  bigatti
// -- renamed VectorOperations --> VectorOps
//
// Revision 1.39  2018/04/20 18:51:25  abbott
// Summary: Changed ctors for BigInt/BigRat from string or from MPZ/MPQ
//
// Revision 1.38  2018/04/18 14:36:21  abbott
// Summary: Commented out fns reported as unused by compiler
//
// Revision 1.37  2018/03/15 14:38:56  bigatti
// -- new file SparsePolyOps-ideal.H
//
// Revision 1.36  2018/01/17 10:29:05  abbott
// Summary: Several minor changes (e.g. adding const)
//
// Revision 1.35  2017/04/28 18:20:58  bigatti
// ++ verbosity for Normaliz: now as for cocoa, NmzSetVerbosityLevel
//
// Revision 1.34  2017/03/13 12:30:37  abbott
// Summary: Move #ifdef guard to after inclusion of header file (so that defns from PREPROCESSOR_DEFNS.H are visible)
//
// Revision 1.33  2016/11/11 14:15:32  abbott
// Summary: Added short-cut to operator<< when ostream is in bad state
//
// Revision 1.32  2015/09/10 16:04:48  abbott
// Summary: Improved a vague comment
//
// Revision 1.31  2015/09/03 10:18:41  bigatti
// -- changes by Christof Soeger (in Aarhus)
//
// Revision 1.30  2015/05/20 15:48:05  abbott
// Summary: Added commented out line about setting the verbose flag (just not to forget)
// Author: JAA
//
// Revision 1.29  2014/07/31 14:45:17  abbott
// Summary: Merged io.H and UtilsTemplate.H into new header VectorOperations.H
// Author: JAA
//
// Revision 1.28  2014/07/14 10:02:56  abbott
// Summary: Christof has added some fns which use quasi polys
// Author: JAA
//
// Revision 1.27  2014/07/11 10:12:51  abbott
// Summary: Christof added printing for cones
// Author: JAA
//
// Revision 1.26  2014/07/07 12:13:53  abbott
// Summary: Christof's updates (& removed AsSparsePolyRing)
// Author: JAA
//
// Revision 1.25  2014/07/01 15:31:21  bigatti
// -- fix by Christof Soeger
//
// Revision 1.24  2014/06/17 10:06:44  abbott
// Summary: Removed pointless AsPolyRing
// Author: JAA
//
// Revision 1.23  2014/05/12 14:33:51  bigatti
// -- updated from Christof Soeger
//
// Revision 1.22  2014/05/09 14:56:30  bigatti
// -- new fn by Christof Soeger
//
// Revision 1.21  2014/01/30 09:56:49  bigatti
// -- removed DOS newlines
//
// Revision 1.20  2014/01/29 17:38:35  bigatti
// -- added HilbertSeries (by Christof Soeger)
// -- added multiplicity (by Christof Soeger)
//
// Revision 1.19  2014/01/28 09:44:42  abbott
// Tidier impl of MatrixToVecVecBigInt using new ConvertTo syntax.
//
// Revision 1.18  2013/07/12 14:53:32  abbott
// Improved indentation.
//
// Revision 1.17  2013/03/15 17:49:17  abbott
// Minor cleaning; removed an out-of-date TODO comment.
//
// Revision 1.16  2012/10/08 13:52:02  bigatti
// -- more cleaning and updates by Christof Soeger
//
// Revision 1.15  2012/10/05 10:17:04  bigatti
// by Christof Soeger
// * Made the NewCone functions to constructors of cone.
// * Introduced some abbreviation in the method names:
//  IntegralClosure -> IntClosure
//  MonomialIdeal -> MonIdeal
//  Normaliz -> Nmz  (the prefix for CoCoA5 functions)
// * New function IntClosureMonIdeal
//
// Revision 1.14  2012/09/28 13:56:36  abbott
// Cleaned up code with Christof's help -- now cleaner & more symmetrical.
// Also fixed a subtle bug where mySmallConeIsGood was not updated as it should have been.
//
// Revision 1.13  2012/08/03 16:31:22  bigatti
// -- changed: procedural --> functional (by C.Soeger)
//
// Revision 1.12  2012/07/25 12:46:27  bigatti
// -- merged with C.Soeger changes for NormalizComputation (init from map, etc)
//
// Revision 1.11  2012/07/19 17:12:05  abbott
// Added NewCone -- unified pseudo-ctor so user does not have to choose between long and BigInt.
//
// Revision 1.10  2012/06/19 14:44:41  bigatti
// -- changed Ht1 --> Deg1, changed def of HVector (by C.Soeger)
//
// Revision 1.9  2011/11/07 11:09:32  bigatti
// -- new ctors taking ConstRefRingElem
// -- added PPVector NormalToricRing(const PPVector& ppv)
// -- removed void HilbertBasis(std::vector<std::vector<BigInt> >& v, const cone& c)
//
// Revision 1.8  2011/10/04 13:03:02  bigatti
// -- new logo for gui
//
// Revision 1.7  2011/09/30 12:55:44  bigatti
// -- introduced namespace "Normaliz" and removed Normaliz from function names
// -- input of Normaliz functions in CoCoA-5 is now a matrix instead of
//    vector<vector<BigInt>>
//
// Revision 1.6  2011/08/14 15:52:17  abbott
// Changed ZZ into BigInt (phase 1: just the library sources).
//
// Revision 1.5  2011/07/20 15:31:21  bigatti
// -- fixed InputType in (undefined) functions
//
// Revision 1.4  2011/07/20 13:49:37  bigatti
// -- added "Normaliz" postfix to Normaliz function calls
//
// Revision 1.3  2011/07/20 12:45:12  bigatti
// -- new normaliz interface (not yet public)
//
// Revision 1.2  2011/02/17 16:50:04  bigatti
// -- getting ready for new official veson on Normaliz: added HVector, removed Triangulation
//
// Revision 1.1  2010/11/05 14:27:18  bigatti
// -- was TmpNormaliz**
//
// Revision 1.6  2010/10/12 11:22:54  bigatti
// -- TmpNormaliz.H simplified:
//    now cone is a smart pointer to ConeBase
//    and the concrete classes are entirely in the .C file
// -- added Ht1Elements, SupportHyperplanes, Triangulation
// -- added some text in ex-Normaliz1 and 2
//
// Revision 1.5  2010/10/08 13:40:37  bigatti
// -- moved inclusion of libnormaliz.cpp and instantiation of types
//    into TmpNormalizTypes.C
//
// Revision 1.4  2010/10/08 10:39:32  bigatti
// -- extended interface for normaliz
//

#endif