#  Copyright (c) 1997-2024
#  Ewgenij Gawrilow, Michael Joswig, and the polymake team
#  Technische Universität Berlin, Germany
#  https://polymake.org
#
#  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 2, 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.
#-------------------------------------------------------------------------------


# A polytope all of whose vertex coordinates are integral.

declare object_specialization Lattice = Polytope<Rational> {

    precondition : LATTICE && BOUNDED;

}

# This alias is solely introduced for backward compatibility
# with some external extensions like polymake-gap interface.
# It will be removed after some transitional period.
# Please do not use it in polymake native code!
declare object LatticePolytope = Polytope<Rational>;

INCLUDE
  lattice_properties.rules
  lattice_visual.rules
  _4ti2.rules + require_libnormaliz.rules | lattice_missing_hilbert.rules
  _4ti2.rules | lattice_missing_hilbert_nonpointed.rules
  latte.rules + require_libnormaliz.rules | lattice_missing_ehrhart.rules


object Cone<Rational> {

    # @category Lattice points in cones
    # for a cone this method returns a Hilbert basis of the cone
    # for a polytope this method returns a Hilbert basis of the homogenization cone of the polytope
    # note: if the cone is not pointed (the polytope is not bounded) then the returned basis is not unique and usually not minimal
    # @return Matrix<Integer>
    user_method HILBERT_BASIS {
        my $c = shift;
        my $h = $c->HILBERT_BASIS_GENERATORS->[0];
        my $l = $c->HILBERT_BASIS_GENERATORS->[1];
        if ($l->rows) { return $h/$l/(-$l); }
        return $h;
    }


  rule N_HILBERT_BASIS : HILBERT_BASIS_GENERATORS {
    $this->N_HILBERT_BASIS = $this->HILBERT_BASIS_GENERATORS->[0]->rows;
  }
  precondition : POINTED;
  weight 0.1;

  rule SMOOTH_CONE : RAYS, SIMPLICIAL_CONE {
    if ( !$this->SIMPLICIAL_CONE ) {
      $this->SMOOTH_CONE = false;
    } else {
      my $rays = primitive($this->RAYS);
      my $g = 0;
      foreach (@{all_subsets_of_k(sequence(0, $rays->cols), $rays->rows)}) {
        $g = gcd($g,abs(det($rays->minor(All,$_))));
        last if ($g == 1);
      }
      $this->SMOOTH_CONE = ( $g == 1 );
    }
  }
  weight 1.10;

  rule Q_GORENSTEIN_CONE, Q_GORENSTEIN_CONE_INDEX : RAYS, CONE_DIM {
    my $g = q_gorenstein_cone($this->RAYS, $this->CONE_DIM);
    if ( $g->first ) {
      $this->Q_GORENSTEIN_CONE = 1;
      $this->Q_GORENSTEIN_CONE_INDEX = $g->second;
    } else {
      $this->Q_GORENSTEIN_CONE = 0;
    }
  }
  weight 1.10;

  rule GORENSTEIN_CONE : {
    $this->GORENSTEIN_CONE = 0;
  }
  precondition : !Q_GORENSTEIN_CONE;

  rule GORENSTEIN_CONE : Q_GORENSTEIN_CONE, Q_GORENSTEIN_CONE_INDEX {
      $this->GORENSTEIN_CONE = $this->Q_GORENSTEIN_CONE && ( $this->Q_GORENSTEIN_CONE_INDEX == 1 );
  }
  weight 0.1;

  rule GORENSTEIN_CONE, Q_GORENSTEIN_CONE, Q_GORENSTEIN_CONE_INDEX : {
    $this->GORENSTEIN_CONE = 1;
    $this->Q_GORENSTEIN_CONE = 1;
    $this->Q_GORENSTEIN_CONE_INDEX = 1;
  }
  precondition : SMOOTH_CONE;
  weight 0.1;


  rule HOMOGENEOUS : RAYS, LINEAR_SPAN {
      my $r = new Matrix<Rational>(primitive($this->RAYS));
      my $ls = $this->LINEAR_SPAN;
      $r = $r - repeat_row($r->[0],$r->rows);
      my $ns = null_space($r->minor(~[0],All));
      $this->HOMOGENEOUS = ( rank($ls/$ns) - $ls->rows > 0 );
  }
  precondition : POINTED;
  weight 1.10;

  rule cone_only : MONOID_GRADING : RAYS {
      my $r = primitive($this->RAYS);
      $r = $r - repeat_row($r->[0],$r->rows);
      my $ns = null_space($r->minor(~[0],All));
      $ns *= -1 if ($ns->rows > 0 and $ns->[0] * $this->RAYS->[0] < 0);
      $this->MONOID_GRADING = $ns->rows > 0 ? primitive($ns->[0]) : undef;
  }
  precondition : FULL_DIM ;
  precondition : POINTED ;
  precondition : HOMOGENEOUS ;
  weight 1.10;

  rule H_STAR_POLYNOMIAL : H_STAR_VECTOR {
      $this->H_STAR_POLYNOMIAL = new UniPolynomial($this->H_STAR_VECTOR, sequence(0, $this->H_STAR_VECTOR->dim));
  }
  weight 0.10;

  rule H_STAR_VECTOR : H_STAR_POLYNOMIAL , HILBERT_SERIES {
      my $len = denominator($this->HILBERT_SERIES)->deg;
      my $h_star_vector = new Vector( zero_vector<Integer>($len) );
      my $poly = $this->H_STAR_POLYNOMIAL;

      while($poly->deg>=0){
          my $coef = $poly->lc;
          my $deg = $poly->deg;
          $h_star_vector->[$deg] = $coef;
          $poly = $poly - new UniPolynomial([$coef],[$deg]);
      }
      $this->H_STAR_VECTOR = $h_star_vector;
  }
  weight 0.10;

  rule CONE_TORIC_IDEAL.BINOMIAL : {
      $this->CONE_TORIC_IDEAL->BINOMIAL = true;
  }
  weight 0.1;

}

object Polytope {

    rule N_LATTICE_POINTS : LATTICE_POINTS_GENERATORS {
        $this->N_LATTICE_POINTS = $this->LATTICE_POINTS_GENERATORS->[0]->rows;
    }
    precondition : BOUNDED;
    weight 0.1;

    rule BOUNDARY_LATTICE_POINTS, INTERIOR_LATTICE_POINTS : LATTICE_POINTS_GENERATORS , FACETS {
        my $lattice_points = $this->LATTICE_POINTS_GENERATORS->[0];
        my $facets = $this->FACETS;

        my (@blp, @ilp);
        foreach my $l (@$lattice_points) {
            my $interior=1;
            foreach my $f (@$facets) {
                if ($f * $l == 0) {
                    push @blp, new Vector<Integer>($l);
                    $interior=0;
                    last;
                }
            }
            if ($interior) {
                push @ilp, new Vector<Integer>($l);
            }
        }
        $this->BOUNDARY_LATTICE_POINTS=\@blp;
        $this->INTERIOR_LATTICE_POINTS=\@ilp;
    }
    precondition : BOUNDED;
    weight 2.10;

    rule N_BOUNDARY_LATTICE_POINTS : N_LATTICE_POINTS , N_INTERIOR_LATTICE_POINTS {
        $this->N_BOUNDARY_LATTICE_POINTS = $this->N_LATTICE_POINTS - $this->N_INTERIOR_LATTICE_POINTS;
    }
    weight 0.1;

    rule N_BOUNDARY_LATTICE_POINTS : BOUNDARY_LATTICE_POINTS {
        $this->N_BOUNDARY_LATTICE_POINTS = $this->BOUNDARY_LATTICE_POINTS->rows;
    }
    weight 0.1;

    rule N_INTERIOR_LATTICE_POINTS : INTERIOR_LATTICE_POINTS {
        $this->N_INTERIOR_LATTICE_POINTS = $this->INTERIOR_LATTICE_POINTS->rows;
    }
    weight 0.1;

    rule N_INTERIOR_LATTICE_POINTS : N_LATTICE_POINTS , N_BOUNDARY_LATTICE_POINTS {
        $this->N_INTERIOR_LATTICE_POINTS = $this->N_LATTICE_POINTS - $this->N_BOUNDARY_LATTICE_POINTS;
    }
    weight 0.1;

    rule LATTICE_WIDTH, LATTICE_WIDTH_DIRECTION : VERTICES | POINTS, LINEALITY_SPACE, CONE_AMBIENT_DIM {
       my $direction = unit_vector($this->CONE_AMBIENT_DIM,1);
       my $lp=$this->LP(LINEAR_OBJECTIVE=>$direction, temporary);
       my $sched = $lp->get_schedule("MAXIMAL_VALUE","MINIMAL_VALUE");
       $sched->apply($lp);
       my $width=$lp->MAXIMAL_VALUE - $lp->MINIMAL_VALUE;

       # we need LINEALITY_SPACE and VERTICES because of minkowski_sum.
       my $half_space = unit_vector<Scalar>($this->CONE_AMBIENT_DIM,1) + (0|(ones_vector<Scalar>($this->CONE_AMBIENT_DIM-1)/$width));
       my $S=new Polytope<Scalar>(INEQUALITIES=>(scale(polarize(minkowski_sum(1,$this, -1,$this)),$width)->FACETS)/$half_space);
       my $B=$S->LATTICE_POINTS;

       for (my $i=0; $i<$B->rows(); ++$i) {
           my $current_direction = 0 | primitive($B->row($i)->slice(range_from(1)));
           next if is_zero($current_direction); #because 0 is a stupid direction

           $lp=$this->LP(LINEAR_OBJECTIVE=>$current_direction, temporary);

           $sched->apply($lp);

           my $current_width=$lp->MAXIMAL_VALUE-$lp->MINIMAL_VALUE;

           if ($current_width<$width){
               $width=$current_width;
               $direction=$current_direction;
           }
       }
       $this->LATTICE_WIDTH=$width;
       $this->LATTICE_WIDTH_DIRECTION=$direction;
    }
    precondition : FULL_DIM;
    weight 5.50;

    # @category Lattice points in polytopes
    # Returns the lattice points in bounded Polytopes.
    # @return Matrix<Integer>
    user_method LATTICE_POINTS {
      my $p = shift;
      if (!($p->BOUNDED)) { die "input polytope must be bounded" }
      return $p->LATTICE_POINTS_GENERATORS->[0];
    }


label projection

    rule projection.integer_points : LATTICE_POINTS_GENERATORS : FACETS , AFFINE_HULL , VERTICES_IN_FACETS , VERTICES , CONE_AMBIENT_DIM , CONE_DIM {
        $this->LATTICE_POINTS_GENERATORS = [integer_points_projection($this,0), new Matrix(0,$this->CONE_AMBIENT_DIM),new Matrix(0,$this->CONE_AMBIENT_DIM)];
    }
    precondition : BOUNDED;
    weight 5.10;

label bbox

    rule bbox.integer_points : LATTICE_POINTS_GENERATORS : FACETS | INEQUALITIES , AFFINE_HULL | EQUATIONS, CONE_AMBIENT_DIM {
        $this->LATTICE_POINTS_GENERATORS = [integer_points_bbox($this), new Matrix(0,$this->CONE_AMBIENT_DIM), new Matrix(0,$this->CONE_AMBIENT_DIM)];
    }
    precondition : BOUNDED;
    weight 5.20;
}

object Polytope<Rational> {

    rule LATTICE_POINTS_GENERATORS : HILBERT_BASIS_GENERATORS, CONE_AMBIENT_DIM {
        $this->LATTICE_POINTS_GENERATORS = [ [ grep { $_->[0] == 1 } @{$this->HILBERT_BASIS_GENERATORS->[0]} ], new Matrix(0,$this->CONE_AMBIENT_DIM),new Matrix(0,$this->CONE_AMBIENT_DIM)];
    }
    precondition : BOUNDED;
    weight 1.10;

    rule N_LATTICE_POINTS : HILBERT_BASIS_GENERATORS {
        my $hilb = $this->HILBERT_BASIS_GENERATORS->[0];
        my $count = 0;
        foreach (@{$hilb->col(0)}) {
            if ($_ == 1) {
                $count++;
            }
        }
        $this->N_LATTICE_POINTS = $count;
    }
    precondition : BOUNDED;
    weight 1.10;

    rule N_INTERIOR_LATTICE_POINTS , N_LATTICE_POINTS : {
      $this->N_LATTICE_POINTS = 0;
      $this->N_INTERIOR_LATTICE_POINTS = 0;
    }
    weight 0.1;
    precondition : !FEASIBLE;

    rule LATTICE_POINTS_GENERATORS : DEGREE_ONE_GENERATORS, CONE_AMBIENT_DIM {
      $this->LATTICE_POINTS_GENERATORS = [$this->DEGREE_ONE_GENERATORS,new Matrix(0,$this->CONE_AMBIENT_DIM),new Matrix(0,$this->CONE_AMBIENT_DIM)];
    }
    precondition : BOUNDED;
    precondition : MONOID_GRADING { $this->MONOID_GRADING == unit_vector<Integer>($this->MONOID_GRADING->dim,0); };
    weight 0.10;


    rule MINKOWSKI_CONE : HASSE_DIAGRAM.ADJACENCY, HASSE_DIAGRAM.DECORATION, HASSE_DIAGRAM.INVERSE_RANK_MAP, HASSE_DIAGRAM.TOP_NODE, HASSE_DIAGRAM.BOTTOM_NODE, \
         GRAPH.ADJACENCY, GRAPH.EDGE_DIRECTIONS, FAR_FACE {
        $this->MINKOWSKI_CONE = minkowski_cone($this->HASSE_DIAGRAM, $this->GRAPH->ADJACENCY, $this->GRAPH->EDGE_DIRECTIONS, $this->FAR_FACE);
    }
    precondition : CONE_DIM { $this->CONE_DIM > 2; }

    # @category Geometry
    # returns the Minkowski summand of a polytope P given by
    # a coefficient vector to the rays of the [[MINKOWSKI_CONE]].
    # @param Vector<Rational> coeff coefficient vector to the rays of the Minkowski summand cone
    # @return Polytope<Rational>
    user_method MINKOWSKI_CONE_COEFF {
        my ($p, $coeff)  = @_;
        return minkowski_cone_coeff($coeff, $p->MINKOWSKI_CONE, $p->GRAPH, $p->FAR_FACE, $p->VERTICES->minor($p->FAR_FACE, All));
    }

    # @category Geometry
    # returns the Minkowski summand of a polytope P given by
    # a point in the [[MINKOWSKI_CONE]].
    # @param Vector<Rational> point point in the Minkowski summand cone
    # @return Polytope<Rational>
    user_method MINKOWSKI_CONE_POINT {
        my $p = shift;
        my $point = new Vector<Rational>(shift);
        my $c = $p->give("MINKOWSKI_CONE");
        foreach ( @{$point} ) {
            if ( $_ < 0 ) {  die "given point is not in Minkowski cone"; }
        }
        # foreach ( @{$c->LINEAR_SPAN * $point} ) {
        #     if ( $_ != 0 ) {  die "given point is not in Minkowski cone"; }
        # }
        my $vertices = $p->give("VERTICES");
        return minkowski_cone_point($point, $vertices->minor($p->FAR_FACE, All), $p->GRAPH, $p->FAR_FACE);
    }

    rule MONOID_GRADING : CONE_AMBIENT_DIM {
       $this->MONOID_GRADING = unit_vector<Integer>($this->CONE_AMBIENT_DIM,0);
    }
    precondition : BOUNDED;
    weight 0.10;
    override : SUPER::cone_only;

}


object_specialization Polytope::Lattice {

    rule TORIC_IDEAL.BINOMIAL : {
        $this->TORIC_IDEAL->BINOMIAL = true;
    }
    weight 0.10;

    rule TORIC_IDEAL = CONE_TORIC_IDEAL;
    precondition : NORMAL;

    rule CONE_TORIC_IDEAL = TORIC_IDEAL;
    precondition : NORMAL;

    rule LATTICE_BASIS : CONE_AMBIENT_DIM {
        $this->LATTICE_BASIS = unit_matrix<Rational>($this->CONE_AMBIENT_DIM-1);
    }
    weight 0.10;

    # @category Lattice points in cones
    # returns a polytope in the integer lattice basis if a [[LATTICE_BASIS]] is given
    # @param Polytope<Rational> P polytope
    # @return Polytope<Rational> Pnew polytope
    user_method POLYTOPE_IN_STD_BASIS {
        my $this=shift;
        my $p=new Polytope<Rational>;
        if (my ($v, $prop_name)=$this->lookup_with_name("VERTICES | POINTS")) {
            $p->take($prop_name, ones_vector<Rational>($v->rows) | $v->minor(All,~[0])*$this->LATTICE_BASIS);
        }

        if (my ($f, $prop_name)=$this->lookup_with_name("FACETS | INEQUALITIES")) {
            $p->take($prop_name, $f->col(0) | $f->minor(All,~[0])*inv($this->LATTICE_BASIS));
        }

        if (defined (my $f=$this->lookup("VERTICES_IN_FACETS"))) {
            $p->VERTICES_IN_FACETS=$f;
        }

        return $p;
    }

    rule N_LATTICE_POINTS : EHRHART_POLYNOMIAL {
        my $e = $this->EHRHART_POLYNOMIAL->coefficients_as_vector();
        my $dim = $e->dim();
        my $num = new Rational(0);
        for (my $i = 0; $i < $dim; $i++) {
            $num += $e->[$i];
        }
        $this->N_LATTICE_POINTS=$num;
    }
    weight 0.10;

    rule N_LATTICE_POINTS : H_STAR_VECTOR {
        $this->N_LATTICE_POINTS = $this->H_STAR_VECTOR->[1] + scalar(@{$this->H_STAR_VECTOR});
    }
    precondition : CONE_DIM { $this->CONE_DIM > 1; }
    weight 0.1;


    # @category Lattice points in cones
    # The number of [[LATTICE_POINTS]] in the //n//-th dilation of the polytope
    # @param Int n dilation factor
        # @return Int
    user_method N_LATTICE_POINTS_IN_DILATION {
        my ($this, $n)=@_;
        return abs($this->EHRHART_POLYNOMIAL->evaluate($n));
    }

    rule N_INTERIOR_LATTICE_POINTS : EHRHART_POLYNOMIAL {
        $this->N_INTERIOR_LATTICE_POINTS = abs($this->EHRHART_POLYNOMIAL->evaluate(-1));
    }
    weight 0.10;

    rule N_INTERIOR_LATTICE_POINTS : H_STAR_VECTOR {
        $this->N_INTERIOR_LATTICE_POINTS = $this->H_STAR_VECTOR->[-1];
    }
    weight 0.1;

        rule N_INTERIOR_LATTICE_POINTS : {
                $this->N_INTERIOR_LATTICE_POINTS = 1;
        }
        precondition : REFLEXIVE;
        weight 0.10;

    # Check, whether the polar dual is also a lattice polytope.
    rule REFLEXIVE : FACETS {
        my $n_int = $this->lookup("N_INTERIOR_LATTICE_POINTS");
        if (defined($n_int) && $n_int != 1) {
            $this->REFLEXIVE = 0;
            return;
        }
        my $fac = eliminate_denominators_in_rows($this->FACETS);
        foreach (@$fac) {
            my $a = $_->[0];
            for (my $i=1; $i<scalar(@$_); $i++) {
                if (($_->[$i])%$a != 0) {
                    $this->REFLEXIVE = 0;
                    return;
                }
            }
        }
        $this->REFLEXIVE = 1;
    }
    precondition : CENTERED;
    precondition : FULL_DIM;
    weight 1.10;

    rule REFLEXIVE : GORENSTEIN , GORENSTEIN_INDEX {
        $this->REFLEXIVE= $this->GORENSTEIN && ($this->GORENSTEIN_INDEX == 1);
    }
    precondition : FULL_DIM;
    weight 0.1;

    rule REFLEXIVE : H_STAR_VECTOR {
        my $h = $this->H_STAR_VECTOR;
        my $dim = $#$h;
        my $i = 0;
        if ($this->H_STAR_VECTOR->[$dim] == 0) {
            $this->REFLEXIVE = 0;
            return;
        }
        while ($i < $dim-$i) {
            if($this->H_STAR_VECTOR->[$i] != $this->H_STAR_VECTOR->[$dim-$i]) {
                $this->REFLEXIVE = 0;
                return;
            }
            $i++;
        }
        $this->REFLEXIVE = 1;
    }
    precondition : FEASIBLE;
    weight 0.10;

    rule GORENSTEIN , GORENSTEIN_INDEX : {
        $this->GORENSTEIN = 1;
        $this->GORENSTEIN_INDEX = 1;
    }
    precondition : REFLEXIVE;
    weight 0.1;

    rule GORENSTEIN : H_STAR_VECTOR {
        my $h = $this->H_STAR_VECTOR;
        my $max = $#$h;
        while ($this->H_STAR_VECTOR->[$max] == 0) {
            $max--;
        }
        my $i = 0;
        while ($i < $max-$i) {
            if($this->H_STAR_VECTOR->[$i] != $this->H_STAR_VECTOR->[$max-$i]) {
                $this->GORENSTEIN = 0;
                return;
            }
            $i++;
        }
        $this->GORENSTEIN = 1;
    }
    precondition : FEASIBLE;
    weight 0.10;

    #Read FACETS f1, . . . , fk (in lattice minimal form) and check whether there is a lattice minimal x ??? L such that hx, fii = 1 for all 1 . i . k. If so, then P is Gorenstein, and the first coordinate r of x of is the Gorenstein index of P. x is the unique inner lattice point of rP
    rule GORENSTEIN , GORENSTEIN_INDEX , GORENSTEIN_VECTOR : FACETS , CONE_AMBIENT_DIM {
        my $p = new Polytope<Rational>(POINTS => ones_vector<Integer>() | primitive($this->FACETS));
        if ($p->CONE_DIM <= $this->CONE_AMBIENT_DIM) {
            my $x = dehomogenize($p->AFFINE_HULL->row(0));
            if (is_integral($x)) {
                $this->GORENSTEIN = 1;
                my $lx = primitive($x->[0]>0 ? $x : -$x);
                $this->GORENSTEIN_INDEX = abs($lx->[0]);
                $lx->[0]=1;
                $this->GORENSTEIN_VECTOR = $lx;
                return;
            }
        }
        $this->GORENSTEIN = 0;
    }
    precondition : FULL_DIM;
    precondition : FEASIBLE;
    weight 3.10;

    rule GORENSTEIN_VECTOR : INTERIOR_LATTICE_POINTS {
        $this->GORENSTEIN_VECTOR = $this->INTERIOR_LATTICE_POINTS->row(0);
    }
    precondition : REFLEXIVE;
    weight 0.10;

    rule CANONICAL : N_INTERIOR_LATTICE_POINTS {
        $this->CANONICAL = ($this->N_INTERIOR_LATTICE_POINTS == 1);
    }
    weight 0.1;

    rule TERMINAL : CANONICAL , N_LATTICE_POINTS , N_VERTICES {
        $this->TERMINAL = ($this->CANONICAL && ($this->N_LATTICE_POINTS == $this->N_VERTICES + 1));
    }
    weight 0.1;

    rule LATTICE_EMPTY : N_LATTICE_POINTS , N_VERTICES {
        $this->LATTICE_EMPTY= ($this->N_LATTICE_POINTS == $this->N_VERTICES);
    }
    weight 0.10;

    user_method LATTICE_HOLLOW {
        my $p = shift;
        return $p->N_INTERIOR_LATTICE_POINTS == 0;
    }

    rule LATTICE_VOLUME : H_STAR_VECTOR {
        my $h = $this->H_STAR_VECTOR;
        my $dim = $#$h;
        my $vol = new Integer(0);
        $vol += $_ for @{$this->H_STAR_VECTOR};
        $this->LATTICE_VOLUME = $vol;
    }
    weight 0.10;

    rule LATTICE_VOLUME : VOLUME , CONE_DIM {
        $this->LATTICE_VOLUME = $this->VOLUME * fac($this->CONE_DIM-1);
    }
    precondition : FULL_DIM;
    weight 0.10;

    rule VOLUME : LATTICE_VOLUME , CONE_DIM {
        my $vol = new Rational($this->LATTICE_VOLUME);
        $this->VOLUME = $vol / fac($this->CONE_DIM-1);
    }
    precondition : FULL_DIM;
    weight 0.10;

    rule LATTICE_VOLUME : EHRHART_POLYNOMIAL , CONE_DIM {
        $this->LATTICE_VOLUME = $this->EHRHART_POLYNOMIAL->lc() * fac($this->CONE_DIM-1);
    }
    weight 0.10;

    rule LATTICE_DEGREE : LATTICE_CODEGREE , CONE_DIM {
        $this->LATTICE_DEGREE = $this->CONE_DIM - $this->LATTICE_CODEGREE;
    }
    weight 0.1;

    rule LATTICE_DEGREE : H_STAR_VECTOR {
        my $h = $this->H_STAR_VECTOR;
        my $i = $#$h;
        while ($h->[$i] == 0) {
            $i--;
        }
        $this->LATTICE_DEGREE = $i;
    }
    weight 0.10;

    rule LATTICE_CODEGREE : GORENSTEIN_INDEX {
        my $index = $this->GORENSTEIN_INDEX;
        if ($index > 0) {
            $this->LATTICE_CODEGREE = $index;
        }
    }
    precondition : GORENSTEIN;
    weight 0.1;

    rule LATTICE_CODEGREE : LATTICE_DEGREE , CONE_DIM {
        $this->LATTICE_CODEGREE = $this->CONE_DIM - $this->LATTICE_DEGREE;
    }
    weight 0.1;

    rule EHRHART_POLYNOMIAL : H_STAR_VECTOR , CONE_DIM {
        my $coeffs = new Vector<Rational>(binomial_to_power_basis($this->H_STAR_VECTOR->slice(range(0,$this->CONE_DIM-1))));
        my $exps = new Array<Int>(sequence(0,$coeffs->dim()));
        $this->EHRHART_POLYNOMIAL = new UniPolynomial($coeffs,$exps);
    }
    weight 0.10;


    sub get_ehrhart_coeff{
	my $poly = shift;
	my $monomials = new Vector<Int>( $poly->monomials_as_vector() );
	my $coeffs = new Vector<Rational>( $poly->coefficients_as_vector() );
	my $ord_coeffs = new Vector<Rational>( zero_vector<Rational>($poly->deg()+1));
        for (my $i =0; $i<$monomials->dim(); ++$i) {
	    $ord_coeffs->[$monomials->[$i]] = $coeffs->[$i];
	}
        return $ord_coeffs;
    }

    # @category Lattice points in cones
    # Vector containing the coefficients of the [[EHRHART_POLYNOMIAL]], ordered
    # by increasing degree of the corresponding term.
    # @return Vector<Rational>
    user_method EHRHART_POLYNOMIAL_COEFF() : EHRHART_POLYNOMIAL {
        my $this = shift;
        return get_ehrhart_coeff($this->EHRHART_POLYNOMIAL);
    }

    rule H_STAR_VECTOR : EHRHART_POLYNOMIAL , CONE_DIM {
               $this->H_STAR_VECTOR = power_to_binomial_basis(get_ehrhart_coeff($this->EHRHART_POLYNOMIAL)->slice(sequence(0, $this->CONE_DIM)));
    }
    weight 0.10;

    rule NORMAL : HILBERT_BASIS_GENERATORS {
        my $hilb = $this->HILBERT_BASIS_GENERATORS->[0];
        foreach (@{$hilb->col(0)}) {
            if ($_ != 1) {
                $this->NORMAL = 0;
                return;
            }
        }
        $this->NORMAL = 1;
    }
    precondition : BOUNDED;
    weight 1.10;

    rule NORMAL : N_LATTICE_POINTS , N_HILBERT_BASIS {
        $this->NORMAL = ($this->N_LATTICE_POINTS == $this->N_HILBERT_BASIS);
    }
    weight 0.1;

    rule SPANNING : LATTICE_POINTS_GENERATORS, VERTICES {
	my $v = new Vector<Integer>($this->VERTICES->[0]);
	my $A = new Matrix<Integer>($this->LATTICE_POINTS_GENERATORS->[0]);
	$A -= repeat_row($v, $A->rows);
	my $N = smith_normal_form($A)->[0];
	my $vol = 1;
	for my $j (0..$N->cols -2){ #ignoring zero row
		$vol *= $N->elem($j,$j);
	}
	$this->SPANNING = is_one($vol);
    }
    precondition : BOUNDED;

    rule FACET_WIDTH : FACET_WIDTHS {
        $this->FACET_WIDTH = maximum($this->FACET_WIDTHS);
    }
    weight 0.10;

    rule FACET_VERTEX_LATTICE_DISTANCES : VERTICES, FACETS {
        $this->FACET_VERTEX_LATTICE_DISTANCES = primitive($this->FACETS) * transpose(numerators($this->VERTICES));
    }
    precondition : FULL_DIM;


    # @category Lattice points in cones
    # Vector containing the distances of a given point //v// from all facets
    # @param Vector<Rational> v point in the ambient space of the polytope
    # @return Vector<Integer>
    user_method FACET_POINT_LATTICE_DISTANCES($) : FACETS {
        my ($p, $v_in)=@_;
        my $v = new Vector<Integer>([ 1, @$v_in ]);
        $v->dim == $p->CONE_DIM or die "polytope - point dimension mismatch\n";
        return primitive($p->FACETS) * $v;
    }
    precondition : FULL_DIM;

    rule FACET_WIDTHS : VERTICES , FACETS , CONE_AMBIENT_DIM {
        $this->FACET_WIDTHS = [
            map {
                my $f=$_;
                my $max;
                assign_max($max, $f * $_) for @{numerators($this->VERTICES)};
                $max;
            } @{ primitive($this->FACETS) }
        ];
    }
    precondition : FULL_DIM;
    weight 3.10;

    rule COMPRESSED : FACET_WIDTH {
        $this->COMPRESSED=($this->FACET_WIDTH==1);
    }
    weight 0.1;

    rule SMOOTH : N_VERTICES, GRAPH.ADJACENCY, GRAPH.EDGE_DIRECTIONS, SIMPLE {

        if (!$this->SIMPLE) { $this->SMOOTH = false; return; }

        for ( my $v=0; $v < $this->N_VERTICES; ++$v ) {
            my $m=new Matrix<Integer>;
            for (my $e = entire($this->GRAPH->ADJACENCY->out_edges($v)); $e; ++$e) {
                $m /= primitive($this->GRAPH->EDGE_DIRECTIONS->[$$e]);
            }

            my $g = 0;
            foreach (@{all_subsets_of_k(sequence(1, max(0, $m->cols-1)), $m->rows)}) {
                $g = gcd($g,abs(det($m->minor(All,$_))));
                last if ($g == 1);
            }

            if ($g != 1) { $this->SMOOTH = 0; return; }
        }

        $this->SMOOTH = true;
    }
    weight 3.10;

    rule POLAR_SMOOTH : VERTICES, VERTICES_IN_FACETS, INTERIOR_LATTICE_POINTS, SIMPLICIAL, CONE_DIM, CONE_AMBIENT_DIM {

        if (!$this->SIMPLICIAL) { $this->POLAR_SMOOTH = false; return; }

        my $ip = repeat_row($this->INTERIOR_LATTICE_POINTS->[0],$this->CONE_DIM-1);

        foreach my $v (@{$this->VERTICES_IN_FACETS}) {

          my $m=primitive(($this->VERTICES->minor($v,All)-$ip));

          my $g = 0;
          foreach (@{all_subsets_of_k(sequence(1, $this->CONE_AMBIENT_DIM-1), $this->CONE_DIM-1)}) {
            $g = gcd($g,abs(det( $m->minor(All,$_) )));
            last if ($g == 1);
          }

          if ($g != 1) { $this->POLAR_SMOOTH = false; return; }
        }

        $this->POLAR_SMOOTH = true;
    }
    precondition : REFLEXIVE;
    weight 3.10;

    rule VERY_AMPLE : GRAPH.ADJACENCY, GRAPH.EDGE_DIRECTIONS, N_VERTICES, VERTICES, FACETS | INEQUALITIES, AFFINE_HULL | EQUATIONS {
        my $sched;
        for (my $v = 0; $v < $this->N_VERTICES; ++$v) {
            my @dirs = ();
            push @dirs, unit_vector<Rational>($this->VERTICES->cols, 0);
            for (my $e = entire($this->GRAPH->ADJACENCY->out_edges($v)); $e; ++$e) {
                push @dirs, ($e->to_node > $e->from_node) ? -$this->GRAPH->EDGE_DIRECTIONS->[$$e] : $this->GRAPH->EDGE_DIRECTIONS->[$$e];
            }
            my $p = new Cone<Rational>(INPUT_RAYS=>\@dirs);
            $sched ||= $p->get_schedule("HILBERT_BASIS_GENERATORS")
                   || die "don't know how to compute HILBERT_BASIS: please install 4ti2 or activate bundled:libnormaliz\n";
            $sched->apply($p);
            foreach my $h (@{$p->HILBERT_BASIS_GENERATORS->[0]}) {
                next if ($h->[0] == 1);
                foreach (@{$this->give("FACETS | INEQUALITIES")}) {
                    if ($_ * ($h + $this->VERTICES->[$v]) < 0) {
                        $this->VERY_AMPLE = false;
                        return;
                    }
                }
                my $AH = $this->give("AFFINE_HULL | EQUATIONS");
                foreach (@$AH) {
                    if ($_ * ($h + $this->VERTICES->[$v]) != 0) {
                        $this->VERY_AMPLE = false;
                        return;
                    }
                }
            }
        }
        $this->VERY_AMPLE = true;
    }
    weight 5.50;

    rule VERY_AMPLE : { $this->VERY_AMPLE = true; }
    precondition : SMOOTH;
    weight 0.1;

    rule VERY_AMPLE : { $this->VERY_AMPLE = true; }
    precondition : NORMAL;
    weight 0.1;

    rule GRAPH.LATTICE_EDGE_LENGTHS : VERTICES, GRAPH.ADJACENCY {
        my $el = new EdgeMap<Undirected, Integer>($this->GRAPH->ADJACENCY);
        my $vert = numerators($this->VERTICES);

        for ( my $e=entire(edges($this->GRAPH->ADJACENCY)); $e; ++$e ) {
            my $w= $vert->[$e->from_node]-$vert->[$e->to_node];
            $el->[$$e] = gcd($w);
        }

        $this->GRAPH->LATTICE_EDGE_LENGTHS=$el;
    }

    rule GRAPH.LATTICE_ACCUMULATED_EDGE_LENGTHS : GRAPH.LATTICE_EDGE_LENGTHS {
        my $m = new Map<Integer, Int>;
        foreach my $e (@{$this->GRAPH->LATTICE_EDGE_LENGTHS}) {
            $m->{$e}++;
        }

        $this->GRAPH->LATTICE_ACCUMULATED_EDGE_LENGTHS=$m;
    }

    rule HOMOGENEOUS : {
       $this->HOMOGENEOUS = 1;
    }
    weight 0.1;


    rule FOLDABLE_COCIRCUIT_EQUATIONS : COMBINATORIAL_DIM, RAYS, RAYS_IN_FACETS, INTERIOR_RIDGE_SIMPLICES, MAX_INTERIOR_SIMPLICES {
       $this->FOLDABLE_COCIRCUIT_EQUATIONS = foldable_cocircuit_equations($this->COMBINATORIAL_DIM, $this->RAYS, $this->RAYS_IN_FACETS, $this->INTERIOR_RIDGE_SIMPLICES, $this->MAX_INTERIOR_SIMPLICES);
    }

}


# Local Variables:
# mode: perl
# cperl-indent-level:4
# indent-tabs-mode:nil
# End: