#  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.
#-------------------------------------------------------------------------------

function far_points(Matrix) : c++ (include => "polymake/linalg.h");


object Polytope {

# Here we have to compute the intersection of the cone with one linear halfspace.
# This is done implicitly in the initial rules.  Beware: the combinatorics may change.

# with the check for already_homogeneous we provide for the following special case:
# we might have a cone that already *is* a homogenized polytope, 
# e.g. when we extract one polytope from a polyhedral complex (which is done via the corresponding function for the fan).
# for this, we just run over the first column of  RAYS and LINEALITY_SPACE and check whether we only have zeros and ones
# we also need to make sure that the polytope has at least one bounded vertex
# if true, we make a true copy, otherwise we use FACETS to construct the polytope
method construct(Cone) {
  my ($proto, $cone)=@_;
  my $already_homogeneous=0;
  my $have_point=0;
  if (defined(my $verts=$cone->lookup("RAYS"))) {
      $already_homogeneous=1;
      foreach (@{$verts->col(0)}) {
         if ($_ == 1) {
            $have_point = 1;
         } elsif ($_ != 0) {
            $already_homogeneous = 0;
            last;
         }
      }
      my $lin_space;
      if ( $already_homogeneous && defined( $lin_space=$cone->lookup("LINEALITY_SPACE") ) && $lin_space->rows ) {
          foreach (@{$lin_space->col(0)}) { if ( $_!=0 && $_!=1 ) { $already_homogeneous=0; last; } }
      }
  }
  if ($already_homogeneous && $have_point) {
      Polymake::Core::BigObject::new_copy($proto, $cone);
  } else {
      Polymake::Core::BigObject::new_filled($proto, $cone->name, CONE_AMBIENT_DIM=>$cone->CONE_AMBIENT_DIM, INEQUALITIES=>$cone->FACETS, EQUATIONS=>$cone->LINEAR_SPAN);
  }
}

rule LINEALITY_SPACE : FACETS | INEQUALITIES {
   my $F;
   if (defined (my $AH=$this->lookup("AFFINE_HULL | EQUATIONS"))) {
      $F=$this->give("FACETS | INEQUALITIES") / $AH;
   } else {
      $F=$this->give("FACETS | INEQUALITIES");
   }
   my $cmp=$this->prepare_computations;
   $this->LINEALITY_SPACE=common::lineality_space($F);
}
precondition : CONE_AMBIENT_DIM { $this->CONE_AMBIENT_DIM >= 1 }
precondition : FEASIBLE;
weight 1.10;
override : SUPER::cone_only;


rule VALID_POINT : VERTICES | POINTS {
   my $p=$this->give("VERTICES | POINTS");
   if ($p->rows > 0) {
      $this->VALID_POINT(temporary) = $p->row(0);
   }
}
weight 0.10;


# IMPORTANT
# For the combinatorial rules below bear in mind the following semantics:
# (1) The combinatorics of an unbounded and pointed polyhedron is defined to be the combinatorics
#     of the projective closure.
# (2) The combinatorics of an unbounded polyhedron which is not pointed is defined to be the
#     combinatorics of the quotient modulo the lineality space.
# Therefore: Each property which is grouped under "Combinatorics" always refers to a polytope.

# simplex case
rule SIMPLICITY, SIMPLICIALITY : COMBINATORIAL_DIM {
   $this->SIMPLICIALITY=$this->SIMPLICITY=$this->COMBINATORIAL_DIM;
}
precondition : COMBINATORIAL_DIM, N_VERTICES {
   $this->COMBINATORIAL_DIM >= 0 && $this->COMBINATORIAL_DIM+1 == $this->N_VERTICES
}
weight 0.1;

rule SIMPLICIALITY : F_VECTOR, F2_VECTOR {
   # a polytope is k-simplicial if each k-face is a simplex,
   # i.e. f_{0,k} = (k+1) * f_k
   my $k=0;
   my $d=$#{$this->F_VECTOR};
   ++$k while ($k<=$d && $this->F2_VECTOR->elem(0,$k) == ($k+1)*$this->F_VECTOR->[$k]);
   $this->SIMPLICIALITY=$k-1;
}
precondition : COMBINATORIAL_DIM, N_VERTICES {
   $this->COMBINATORIAL_DIM+1 != $this->N_VERTICES
}

rule SIMPLICITY : F_VECTOR, F2_VECTOR {
   # dual notion to simpliciality,
   # i.e. f_{k,d-1} = (k+1) * f_k
   my $k=0;
   my $d=$#{$this->F_VECTOR};
   ++$k while ($k<=$d && $this->F2_VECTOR->elem($d-$k,$d) == ($k+1)*$this->F_VECTOR->[$d-$k]);
   $this->SIMPLICITY=$k-1;
}
precondition : COMBINATORIAL_DIM, N_VERTICES {
   $this->COMBINATORIAL_DIM+1 != $this->N_VERTICES
}

rule FACE_SIMPLICITY : COMBINATORIAL_DIM {
   $this->FACE_SIMPLICITY=$this->COMBINATORIAL_DIM;
}
precondition : SIMPLE;

rule FACE_SIMPLICITY : F2_VECTOR {
   # a polytope is k-face-simple if all k-faces are simple
   # i.e. 2 f_{1,k} = f_{0,1,k} = k * f_{0,k}
   my $k=1;
   my $dim=$this->F2_VECTOR->rows-1;
   ++$k while ($k<=$dim && 2 * $this->F2_VECTOR->elem($k,1) == $k * $this->F2_VECTOR->elem($k,0));
   $this->FACE_SIMPLICITY=$k-1;
}
precondition : !SIMPLE;

rule CUBICALITY : F2_VECTOR {
   # a polytope is k-cubical if its graph is bipartite and all j-faces have 2^j vertices for all j<=k
   # i.e. f_{j,0} = 2^j * f_j
   my $k=1;
   my $k2=2; # 2^k
   my $dim=$this->F2_VECTOR->rows-1;
   my $n_vertices=$this->F2_VECTOR->elem(0,0);
   while ($k<$dim) {
      $k2*=2;
      last if $this->F2_VECTOR->elem($k+1,0) != $k2 * $this->F2_VECTOR->elem($k+1,$k+1);
      ++$k;
   }
   ++$k if ($k==$dim && $n_vertices==2*$k2);    # it's a cube
   $this->CUBICALITY=$k;
}
precondition : CONE_DIM { $this->CONE_DIM >= 2 };
precondition : GRAPH.BIPARTITE;

rule N_VERTEX_FACET_INC : F2_VECTOR {
   $this->N_VERTEX_FACET_INC=$this->F2_VECTOR->[0]->[-1];
}
weight 0.1;

rule CUBICAL : FACET_SIZES {
   # this is a special case of the previous rule which is useful for 4-polytopes with a very large face lattice
   my $cubical=1;
   foreach my $fs (@{$this->FACET_SIZES}) {
      $cubical=0, last if $fs != 8;
   }
   $this->CUBICAL=$cubical;
}
precondition : GRAPH.BIPARTITE, COMBINATORIAL_DIM {
  $this->GRAPH->BIPARTITE && $this->COMBINATORIAL_DIM==4
}
weight 1.10;

rule GRAPH.ADJACENCY : HASSE_DIAGRAM.ADJACENCY, HASSE_DIAGRAM.DECORATION, HASSE_DIAGRAM.INVERSE_RANK_MAP, HASSE_DIAGRAM.TOP_NODE, HASSE_DIAGRAM.BOTTOM_NODE {
   $this->GRAPH->ADJACENCY=vertex_graph($this->HASSE_DIAGRAM);
}
precondition : COMBINATORIAL_DIM { $this->COMBINATORIAL_DIM>=2 };
weight 1.20;

rule DUAL_GRAPH.ADJACENCY : HASSE_DIAGRAM.ADJACENCY, HASSE_DIAGRAM.DECORATION, HASSE_DIAGRAM.INVERSE_RANK_MAP, HASSE_DIAGRAM.TOP_NODE, HASSE_DIAGRAM.BOTTOM_NODE {
   $this->DUAL_GRAPH->ADJACENCY=facet_graph($this->HASSE_DIAGRAM);
}
precondition : COMBINATORIAL_DIM { $this->COMBINATORIAL_DIM>=2 };
weight 1.20;

rule CUBICALITY : COMBINATORIAL_DIM {
   $this->CUBICALITY= $this->COMBINATORIAL_DIM>0;
}
precondition : !GRAPH.BIPARTITE;
weight 0.1;

rule CUBICAL : CUBICALITY, COMBINATORIAL_DIM {
   $this->CUBICAL= $this->CUBICALITY >= $this->COMBINATORIAL_DIM-1;
}
weight 0.1;

rule COCUBICALITY : F2_VECTOR {
   # a polytope is k-cocubical if its dual is k-cubical
   # f_{d-1-j,d-1} = 2^j * f_j for all j<=k
   my $k=1;
   my $k2=2; # 2^k
   my $dim=$this->F2_VECTOR->rows-1;
   my $n_facets=$this->F2_VECTOR->elem($dim,$dim);
   while ($k<$dim) {
      $k2*=2;
      last if $this->F2_VECTOR->elem($dim-($k+1),$dim) != $k2 * $this->F2_VECTOR->elem($dim-($k+1),$dim-($k+1));
      ++$k;
   }
   ++$k if ($k==$dim && $n_facets==2*$k2);      # it's a cross polytope
   $this->COCUBICALITY=$k;
}
precondition : DUAL_GRAPH.BIPARTITE;

rule COCUBICALITY : COMBINATORIAL_DIM {
   $this->COCUBICALITY= $this->COMBINATORIAL_DIM>0;
}
precondition : !DUAL_GRAPH.BIPARTITE;
weight 0.1;

rule COCUBICAL : COCUBICALITY, COMBINATORIAL_DIM {
   $this->COCUBICAL= $this->COCUBICALITY >= $this->COMBINATORIAL_DIM-1;
}
weight 0.1;

rule NEIGHBORLINESS : F_VECTOR {
   # a polytope is k-neighborly if any k vertices form a face (which is necessarily simplicial)
   # i.e. f_{k-1} = binomial(n,k), where n = N_VERTICES
   my $d=$#{$this->F_VECTOR};
   if ( $d >= 0 ) {
     my $k=1;
     my $n_vertices=$this->F_VECTOR->[0];
     ++$k while ($k<=$d && $this->F_VECTOR->[$k] == binomial($n_vertices, $k+1));
     $this->NEIGHBORLINESS=$k;
   } else {
     $this->NEIGHBORLINESS=0;
   }
}

rule BALANCE : F_VECTOR {
   # dual notion to neighborliness
   # i.e. f_{d-k} = binomial(n,k), where n = N_FACETS
   my $d=$#{$this->F_VECTOR};
   if ( $d >= 0 ) {
     my $k=0;
     my $n_facets=$this->F_VECTOR->[$d];
     ++$k while ($k<=$d && $this->F_VECTOR->[$d-$k] == binomial($n_facets, $k+1));
     $this->BALANCE=$k;
   } else {
     $this->BALANCE=0;
   }
}

rule SIMPLICIAL : SIMPLICIALITY, COMBINATORIAL_DIM {
   $this->SIMPLICIAL= $this->SIMPLICIALITY >= $this->COMBINATORIAL_DIM-1;
}
weight 0.1;

rule SIMPLE : SIMPLICITY, COMBINATORIAL_DIM {
   $this->SIMPLE= $this->SIMPLICITY >= $this->COMBINATORIAL_DIM-1;
}
weight 0.1;

rule NEIGHBORLY : NEIGHBORLINESS, COMBINATORIAL_DIM {
   $this->NEIGHBORLY= $this->NEIGHBORLINESS >= floor($this->COMBINATORIAL_DIM/2);
}
weight 0.1;

rule BALANCED : BALANCE, COMBINATORIAL_DIM {
   $this->BALANCED= $this->BALANCE >= floor($this->COMBINATORIAL_DIM/2);
}
weight 0.1;

rule FATNESS : F_VECTOR {
   $this->FATNESS = (new Float($this->F_VECTOR->[1]+$this->F_VECTOR->[2]-20))/(new Float($this->F_VECTOR->[0]+$this->F_VECTOR->[3]-10));
}
precondition : COMBINATORIAL_DIM { $this->COMBINATORIAL_DIM==4 }
precondition : N_VERTICES { $this->N_VERTICES>=6 }

rule COMPLEXITY : F_VECTOR, F2_VECTOR {
   $this->COMPLEXITY = (new Float($this->F2_VECTOR->elem(0,3)-20))/(new Float($this->F_VECTOR->[0]+$this->F_VECTOR->[3]-10));
}
precondition : COMBINATORIAL_DIM { $this->COMBINATORIAL_DIM==4 }
precondition : N_VERTICES { $this->N_VERTICES>=6 }

rule TWO_FACE_SIZES : HASSE_DIAGRAM.ADJACENCY, HASSE_DIAGRAM.DECORATION, HASSE_DIAGRAM.INVERSE_RANK_MAP, HASSE_DIAGRAM.TOP_NODE, HASSE_DIAGRAM.BOTTOM_NODE {
   $this->TWO_FACE_SIZES=two_face_sizes($this->HASSE_DIAGRAM);
}
precondition : COMBINATORIAL_DIM {
$this->COMBINATORIAL_DIM>=2
}

rule TWO_FACE_SIZES : GRAPH.ADJACENCY, VERTICES_IN_FACETS {
   two_face_sizes_simple($this);
}
precondition : SIMPLE;
precondition : COMBINATORIAL_DIM {
   $this->COMBINATORIAL_DIM>=2
}

rule SUBRIDGE_SIZES : HASSE_DIAGRAM.ADJACENCY, HASSE_DIAGRAM.DECORATION, HASSE_DIAGRAM.INVERSE_RANK_MAP, HASSE_DIAGRAM.TOP_NODE, HASSE_DIAGRAM.BOTTOM_NODE {
   $this->SUBRIDGE_SIZES=subridge_sizes($this->HASSE_DIAGRAM);
}
precondition : COMBINATORIAL_DIM {
   $this->COMBINATORIAL_DIM>=2
}

rule SUBRIDGE_SIZES : DUAL_GRAPH.ADJACENCY, VERTICES_IN_FACETS {
   subridge_sizes_simple($this);
}
precondition : SIMPLICIAL;
precondition : COMBINATORIAL_DIM {
   $this->COMBINATORIAL_DIM>=2
}

rule GRAPH.EDGE_DIRECTIONS : GRAPH.ADJACENCY, VERTICES, FAR_FACE {
   $this->GRAPH->EDGE_DIRECTIONS = edge_directions($this->GRAPH,$this->VERTICES,$this->FAR_FACE);
}

rule GRAPH.SQUARED_EDGE_LENGTHS : VERTICES, GRAPH.ADJACENCY, FAR_FACE {
   my $el = new EdgeMap<Undirected, Scalar>($this->GRAPH->ADJACENCY);
   my $vert = dehomogenize($this->VERTICES);

   for ( my $e=entire(edges($this->GRAPH->ADJACENCY)); $e; ++$e ) { 
      my $w= $vert->[$e->from_node]-$vert->[$e->to_node];
      $el->[$$e] = $w * $w;
   }
   $this->GRAPH->SQUARED_EDGE_LENGTHS=$el;
}
precondition : FAR_FACE { $this->FAR_FACE->size()==0 }


rule DUAL_GRAPH.DIHEDRAL_ANGLES : FACETS, DUAL_GRAPH.ADJACENCY {
   my $da = new EdgeMap<Undirected, Float>($this->DUAL_GRAPH->ADJACENCY);
   my $facets = $this->FACETS;
   for ( my $e=entire(edges($this->DUAL_GRAPH->ADJACENCY)); $e; ++$e ) {
      $da->[$$e] = dihedral_angle($facets->[$e->from_node],$facets->[$e->to_node]);
   }
   $this->DUAL_GRAPH->DIHEDRAL_ANGLES = $da;
}
precondition : FULL_DIM;
precondition : FAR_FACE { $this->FAR_FACE->size()==0 }

rule BOUNDED : VERTICES | POINTS {
   my $lin_space;
   if (defined($lin_space = $this->lookup("LINEALITY_SPACE | INPUT_LINEALITY")) && $lin_space->rows > 0) {
      $this->BOUNDED = 0;
   } else {
      foreach my $v (@{$this->give("VERTICES | POINTS")}) {
         if ($v->[0]==0) {
            $this->BOUNDED=0;
            return;
         }
      }
      $this->BOUNDED=1;
   }
}
weight 1.10;


rule BOUNDED : FACETS | INEQUALITIES, LINEALITY_SPACE {
  $this->BOUNDED = H_input_bounded($this);
}
weight 3.10;


rule POINTED : {
   $this->POINTED=1;
}
precondition : BOUNDED;
weight 0.1;

rule BOUNDED: {
   $this->BOUNDED=0;
}
precondition : !POINTED;
weight 0.1;


rule SPLITS : VERTICES, FACETS, GRAPH.ADJACENCY, CONE_DIM {
   $this->SPLITS = splits($this->VERTICES,$this->GRAPH->ADJACENCY,$this->FACETS,$this->CONE_DIM-1);
}
precondition : CONE_DIM {
   $this->CONE_DIM > 1;
}
precondition : POINTED;
 
rule SPLIT_COMPATIBILITY_GRAPH.ADJACENCY : FACETS, AFFINE_HULL, SPLITS {
   $this->SPLIT_COMPATIBILITY_GRAPH->ADJACENCY = split_compatibility_graph($this->SPLITS,$this);
}

rule TRIANGULATION.REFINED_SPLITS : SPLITS, VERTICES, TRIANGULATION.FACETS {
   $this->TRIANGULATION->REFINED_SPLITS=splits_in_subdivision($this->VERTICES,$this->TRIANGULATION->FACETS,$this->SPLITS);
}

rule WEAKLY_CENTERED : FACETS, AFFINE_HULL {
   my $zero=0;
   foreach my $f (@{$this->FACETS}) {
      if ($f->[0] < 0) {
         $this->WEAKLY_CENTERED=0;
         return;
      }
   }
   foreach my $f (@{$this->AFFINE_HULL}) {
      if ($f->[0] != $zero) {
         $this->WEAKLY_CENTERED=0;
         return;
      }
   }
   $this->WEAKLY_CENTERED=1;
}
weight 1.10;

rule WEAKLY_CENTERED : VERTICES | POINTS, LINEALITY_SPACE | INPUT_LINEALITY, CONE_AMBIENT_DIM{
   my $pt = new Polytope<Scalar>(POINTS=>[unit_vector<Scalar>($this->CONE_AMBIENT_DIM, 0)]);
    $this->WEAKLY_CENTERED = contains_V_V_via_LP($pt, $this);
}
weight 4.10;

rule CENTERED : FACETS, AFFINE_HULL {
   foreach my $f (@{$this->FACETS}) {
      if ($f->[0] <= 0) {
         $this->CENTERED=0;
         return;
      }
   }
   foreach my $f (@{$this->AFFINE_HULL}) {
      if ($f->[0] != 0) {
         $this->CENTERED=0;
         return;
      }
   }
   $this->CENTERED=1;
}
weight 1.10;

rule CENTERED : VERTICES | POINTS, LINEALITY_SPACE | INPUT_LINEALITY, CONE_AMBIENT_DIM{
    $this->CENTERED = $this->contains_in_interior(unit_vector($this->CONE_AMBIENT_DIM, 0));
}
weight 4.10;

rule WEAKLY_CENTERED : {
    $this->WEAKLY_CENTERED=1;
}
precondition : CENTERED;
weight 0.1;

rule CENTRALLY_SYMMETRIC, CS_PERMUTATION : VERTICES {
    cs_permutation($this);
}
weight 1.10;

rule TRIANGULATION.BOUNDARY.FACETS, TRIANGULATION.BOUNDARY.VERTEX_MAP, TRIANGULATION.BOUNDARY.FACET_TRIANGULATIONS : TRIANGULATION.FACETS, VERTICES_IN_FACETS {
   ($this->TRIANGULATION->BOUNDARY->FACETS, $this->TRIANGULATION->BOUNDARY->VERTEX_MAP, $this->TRIANGULATION->BOUNDARY->FACET_TRIANGULATIONS)=
     triang_boundary($this->TRIANGULATION->FACETS, $this->VERTICES_IN_FACETS);
}

rule TRIANGULATION(any).BOUNDARY.FACETS, TRIANGULATION(any).BOUNDARY.VERTEX_MAP, TRIANGULATION(any).BOUNDARY.FACET_TRIANGULATIONS : VERTICES_IN_FACETS, N_FACETS, N_VERTICES {
   $this->TRIANGULATION->BOUNDARY->FACETS=[ @{$this->VERTICES_IN_FACETS} ];
   $this->TRIANGULATION->BOUNDARY->VERTEX_MAP = sequence(0,$this->N_VERTICES);
   $this->TRIANGULATION->BOUNDARY->FACET_TRIANGULATIONS=[ map { scalar2set($_) } 0 .. $this->N_FACETS-1 ];
}
precondition : SIMPLICIAL;
weight 1.10;

rule SPECIAL_FACETS : VERTEX_BARYCENTER, VERTICES_IN_FACETS, VERTICES {
   my $v=dehomogenize($this->VERTEX_BARYCENTER);
   my $vif=$this->VERTICES_IN_FACETS;
   my @special_facets=();
   for (my $f=0; $f<$vif->rows(); ++$f) {
      my $vertices_in_f=$this->VERTICES->minor($vif->[$f],~[0]);
      my $cone_f=new Cone<Scalar>(RAYS=>$vertices_in_f);
      my $facets_of_cone=$cone_f->FACETS;
      my $is_special=1;
      for (my $ff=0; $ff<$facets_of_cone->rows(); ++$ff) {
         if ($v*$facets_of_cone->[$ff]<0) {
            $is_special=0;
            last;
         }
      }
      push @special_facets, $f if $is_special;
   }
   $this->SPECIAL_FACETS=\@special_facets;
}
precondition : BOUNDED;
precondition : CENTERED;

# @category Triangulation and volume
# the orientation of the simplices of [[TRIANGULATION_INT]] in the given order of the [[POINTS]]
# @return Array<Int> - +1/-1 array specifying the sign of the determinant of each simplex
user_method TRIANGULATION_INT_SIGNS : TRIANGULATION_INT, POINTS {
   my $this=shift;
   triang_sign($this->TRIANGULATION_INT, $this->POINTS);
}


# clients capable of handling different coordinate types


rule VERTEX_BARYCENTER : VERTICES, N_VERTICES {
   $this->VERTEX_BARYCENTER=$this->N_VERTICES ? barycenter($this->VERTICES) : undef;
}
precondition : BOUNDED;
weight 1.10;

rule default.volume: CENTROID, VOLUME : VERTICES, TRIANGULATION.FACETS {
   centroid_volume($this, $this->VERTICES, $this->TRIANGULATION->FACETS);
}
precondition : BOUNDED;
precondition : FULL_DIM;
weight 2.30;

rule CENTROID, VOLUME : POINTS, TRIANGULATION_INT {
   centroid_volume($this, $this->POINTS, $this->TRIANGULATION_INT);
}
precondition : BOUNDED;
precondition : FULL_DIM;
weight 2.40;

rule VOLUME : {
   $this->VOLUME = 0;
}
precondition : !FULL_DIM;
weight 0.1;

rule VOLUME : {
   $this->VOLUME = 0;
}
precondition : CONE_AMBIENT_DIM {
   $this->CONE_AMBIENT_DIM == 0;
}
weight 0.1;

rule VOLUME : {
   $this->VOLUME = 1;
}
precondition : FULL_DIM;
precondition : CONE_AMBIENT_DIM {
   $this->CONE_AMBIENT_DIM == 1;
}
weight 0.1;

rule default.squared_relative_volumes: SQUARED_RELATIVE_VOLUMES : VERTICES, TRIANGULATION.FACETS {
   $this->SQUARED_RELATIVE_VOLUMES=squared_relative_volumes($this->VERTICES, $this->TRIANGULATION->FACETS);
}
precondition : BOUNDED;
precondition : CONE_DIM {
   $this->CONE_DIM>1;
}
weight 2.20;

rule default.squared_relative_volumes: SQUARED_RELATIVE_VOLUMES : POINTS, TRIANGULATION_INT {
   $this->SQUARED_RELATIVE_VOLUMES=squared_relative_volumes($this->POINTS, $this->TRIANGULATION_INT);
}
precondition : BOUNDED;
precondition : CONE_DIM {
   $this->CONE_DIM>1;
}
weight 2.30;

rule MAHLER_VOLUME : VOLUME {
  $this->MAHLER_VOLUME = $this->VOLUME * polarize($this)->VOLUME;
  # calling "polarize" makes it necessary to add WEAKLY_CENTERED to the preconditions below
  # even if it is implied by CENTERED
}
precondition : BOUNDED;
precondition : WEAKLY_CENTERED;
precondition : CENTERED;
precondition : FULL_DIM;
    
rule MINIMAL_VERTEX_ANGLE : VERTICES, VERTEX_BARYCENTER {
   $this->MINIMAL_VERTEX_ANGLE = minimal_vertex_angle($this);
}

# Triangulation

# quick fix for #302
# compute minimal non-faces of a simplicial polytope 
# using the corresponding functions from application topaz
rule MINIMAL_NON_FACES : HASSE_DIAGRAM.ADJACENCY, HASSE_DIAGRAM.DECORATION, HASSE_DIAGRAM.INVERSE_RANK_MAP, HASSE_DIAGRAM.TOP_NODE, HASSE_DIAGRAM.BOTTOM_NODE {    
    $this->MINIMAL_NON_FACES=topaz::minimal_non_faces($this->HASSE_DIAGRAM);
}
weight 4.10;
precondition : SIMPLICIAL;
precondition : CONE_DIM { $this->CONE_DIM <= 12; }
# this is not really a strict precondition for the algorithm
# it used to be implemented as a dynamic weight but those don't work as expected (#7)
# and this causes the lawler client to run for some medium sized (5-dim) examples where it
# is a lot slower than the one via the hasse-diagram

rule MINIMAL_NON_FACES : VERTICES_IN_FACETS, CONE_DIM, N_VERTICES {
    my $f = new Array<Set<Int>>(@{$this->VERTICES_IN_FACETS});
    $this->MINIMAL_NON_FACES=topaz::lawler_minimal_non_faces($f,$this->N_VERTICES);
}
weight 4.10;
precondition : SIMPLICIAL;
precondition : CONE_DIM { $this->CONE_DIM > 12; }
# this is not really a strict precondition for the algorithm
# it used to be implemented as a dynamic weight but those don't work as expected (#7)
# and this causes the lawler client to run for some medium sized (5-dim) examples where it
# is a lot slower than the one via the hasse-diagram


rule TRIANGULATION(any).VERTEX_LABELS = VERTEX_LABELS;


rule TRIANGULATION(any).COORDINATES : VERTICES {
   $this->TRIANGULATION->COORDINATES(temporary)=$this->VERTICES->minor(All,~[0]);
}
weight 0.10;

rule TRIANGULATION(any).BALL : {
   $this->TRIANGULATION->BALL=1;
}
weight 0.1;

rule TRIANGULATION(any).VOLUME = VOLUME;

rule TRIANGULATION(any).G_DIM : CONE_AMBIENT_DIM {
   $this->TRIANGULATION->G_DIM=$this->CONE_AMBIENT_DIM-1;
}

rule TRIANGULATION(any).BOUNDARY.SPHERE : {
   $this->TRIANGULATION->BOUNDARY->SPHERE=1;
}
weight 0.1;

# @category Geometry
# returns the dimension of the ambient space of the polytope
# @return Int
user_method AMBIENT_DIM() : CONE_AMBIENT_DIM {
  my ($self)=@_;
  return ($self->CONE_AMBIENT_DIM-1);
}

# @category Geometry
# returns the dimension of the polytope
# @return Int
user_method DIM {
  my ($self)=@_;
  if (!defined ($self->lookup("LINEALITY_SPACE | INPUT_LINEALITY | INPUT_RAYS | RAYS | INEQUALITIES | EQUATIONS | FACETS | LINEAR_SPAN"))) {
    return $self->COMBINATORIAL_DIM;
  }
  return $self->CONE_DIM-1;
}

}


# @category Producing a cone
# retuns the recession cone (tail cone, characteristic cone) of a polytope
# @param Polytope<Scalar> P polytope
# @return Cone<Scalar>

user_function recession_cone<Scalar>(Polytope<Scalar>) {
  my $self=shift;
  my @r= map { $_->[0]==0 ? $_ : (); } @{$self->VERTICES};
  my $rays=new Matrix<Scalar>(@r);
  my $c=new Cone<Scalar>(RAYS=>$rays->minor(All,~[0]), LINEALITY_SPACE=>$self->LINEALITY_SPACE->minor(All,~[0]));
  return $c;
}


# @category Geometry
# Compute the dihedral angle between two (oriented) affine or linear hyperplanes.
# @param Vector<Scalar> H1 : first hyperplane
# @param Vector<Scalar> H2 : second hyperplane
# @option Bool deg output in degrees rather than radians, default is false
# @option Bool cone hyperplanes seen as linear hyperplanes, default is false
# @return Float
# @example
# > $H1 = new Vector(1,5,5);
# > $H2 = new Vector(1,-5,5);
# > print dihedral_angle($H1,$H2,deg=>1);
# | 90

user_function dihedral_angle<Scalar> (Vector<Scalar>, Vector<Scalar>; { deg=>0, cone=>0 }) {
  # this code is slightly more elaborate than necessary in order to compute exactly as long as possible
  my ($H1,$H2,$options) = @_;
  if (!$options->{cone}) {
    $H1=$H1->slice(range_from(1)); $H2=$H2->slice(range_from(1));
  }
  my ($norm1, $norm2) = (sqrt(convert_to<Float>($H1*$H1)), sqrt(convert_to<Float>($H2*$H2)));
  $H1=(1/$norm1)*convert_to<Float>($H1); $H2=(1/$norm2)*convert_to<Float>($H2);
  my $pi=2*Math::Trig::asin(1);
  my $angle=$pi-Math::Trig::acos($H1*$H2);
  if ($options->{deg}) {
    $angle=$angle/$pi*180;
  }
  return $angle;
}

# @category Geometry
# Compute the mean or median distance of the [[VERTICES]] to the [[VERTEX_BARYCENTER]].
# @param Polytope p
# @option Bool median use median instead of arithmetic mean
# @return Float
user_function center_distance (Polytope; {median=>0}) {
    my ($p,$options) = @_;
    die "center_distance only for bounded polytopes!" unless $p->BOUNDED;
    my @dists = map { sqrt(convert_to<Float>(sqr($_))) } @{$p->VERTICES - repeat_row($p->VERTEX_BARYCENTER,$p->N_VERTICES)};
    if ($options->{median}) {
        return median(\@dists);
    } else {
        return average(\@dists);
    }
}

# @category Producing a polytope from polytopes
# Produces the Minkowski sum of //P1// and //P2//.
# @param Polytope P1
# @param Polytope P2
# @return Polytope
# @example [prefer cdd] The following stores the minkowski sum of a square and a triangle in the variable $p
# and then prints its vertices.
# > $p = minkowski_sum(cube(2),simplex(2));
# > print $p->VERTICES;
# | 1 -1 -1
# | 1 2 -1
# | 1 -1 2
# | 1 2 1
# | 1 1 2

user_function minkowski_sum<Scalar> (Polytope<type_upgrade<Scalar>>, Polytope<type_upgrade<Scalar>> ) {
    my ($p1,$p2) = @_;
    return minkowski_sum<Scalar>(1,$p1,1,$p2);
}

# @category Producing a polytope from polytopes
# Produces the polytope //lambda//*//P1//+//mu//*//P2//, where * and + are scalar multiplication
# and Minkowski addition, respectively.
# @param Scalar lambda
# @param Polytope P1
# @param Scalar mu
# @param Polytope P2
# @return Polytope
# @example [prefer cdd] The following stores the minkowski sum of a scaled square and a triangle in the variable $p
# and then prints its vertices.
# > $p = minkowski_sum(2,cube(2),1,simplex(2));
# > print $p->VERTICES;
# | 1 -2 -2
# | 1 3 -2
# | 1 -2 3
# | 1 3 2
# | 1 2 3

user_function minkowski_sum<Scalar> (type_upgrade<Scalar>, Polytope<type_upgrade<Scalar>>, type_upgrade<Scalar>, Polytope<type_upgrade<Scalar>> ) {
    my ($l1,$p1,$l2,$p2) = @_;
    my $v1=$p1->give("VERTICES | POINTS");
    my $v2=$p2->give("VERTICES | POINTS");
    my $p=minkowski_sum_client<Scalar>($l1,$v1,$l2,$v2);
    my $ls1=$p1->give("LINEALITY_SPACE");
    my $ls2=$p2->give("LINEALITY_SPACE");
    my $p_out = new Polytope<Scalar>(POINTS=>$p, INPUT_LINEALITY=>$ls1/$ls2);
    $p_out->description = "Minkowski sum of ".$l1."*".$p1->name." and ".$l2."*".$p2->name;
    return $p_out;
}


# @category Producing a polytope from polytopes
# Computes the ([[Polytope::VERTICES]] of the) __Minkowski sum__ of a list of polytopes using the algorithm by Fukuda described in
#	   Komei Fukuda, From the zonotope construction to the Minkowski addition of convex polytopes, J. Symbolic Comput., 38(4):1261-1272, 2004.
# @param Array<Polytope> summands
# @return Polytope
# @example [nocompare] > $p = minkowski_sum_fukuda([cube(2),simplex(2),cross(2)]);
# > print $p->VERTICES;
# | 1 3 -1
# | 1 3 1
# | 1 -1 -2
# | 1 1 3
# | 1 -1 3
# | 1 2 -2
# | 1 -2 2
# | 1 -2 -1
user_function minkowski_sum_fukuda<E>(Polytope<type_upgrade<E>> +) : c++ (include => "polymake/polytope/minkowski_sum_fukuda.h");

# @category Producing a polytope from scratch
# Create the vertices of a zonotope from a matrix whose rows are input points or vectors.
# @param Matrix M
# @option Bool centered_zonotope default 1
# @return Matrix
# @example [nocompare]
# The following stores the vertices of a parallelogram with the origin as its vertex barycenter and prints them:
# > $M = new Matrix([[1,1,0],[1,1,1]]);
# > print zonotope_vertices_fukuda($M);
# | 1 0 -1/2
# | 1 0 1/2
# | 1 -1 -1/2
# | 1 1 1/2
user_function zonotope_vertices_fukuda<E>(Matrix<E> { centered_zonotope => 1 }) : c++ (include => "polymake/polytope/minkowski_sum_fukuda.h");


# @category Geometry
# Given two matrices L (n x d) and R (m x d) such that (L/R) has rank r, select all (r+1-n)-subsets C of rows of R such that (L,S) or (S,L) is a circuit.
# Optionally, if d > r, a basis H for the orthogonal span of the affine hull of (L/R) may be given.
# @param Matrix L
# @param Matrix R
# @option Matrix H
# @return Array<Set>
# @example Divide the vertex set of the 3-cube into a body diagonal L and six remaining vertices R.
# To find the subsets of R that complete L to a circuit, type
# > $c = cube(3);
# > $L = $c->VERTICES->minor([0,7],All);
# > $R = $c->VERTICES->minor([1,2,3,4,5,6],All);
# > print circuit_completions($L,$R);
# | {0 1 3}
# | {2 4 5}
user_function circuit_completions<Scalar>(Matrix<Scalar>,Matrix<Scalar>; Matrix<Scalar> = new Matrix<Scalar>()) {
    my ($L,$R,$H) = @_;
    if ($H->rows() == 0) {
        $H = new Matrix<Scalar>(0, $R->cols());
    }
    my $a= circuit_completions_impl($L, $R, $H);
    return $a;
}

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