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

# @category Geometry
object Cone {

rule N_INPUT_RAYS : INPUT_RAYS {
    $this->N_INPUT_RAYS=$this->INPUT_RAYS->rows;
}
weight 0.1;

rule N_INPUT_LINEALITY : INPUT_LINEALITY {
    $this->N_INPUT_LINEALITY=$this->INPUT_LINEALITY->rows;
}
weight 0.1;

rule N_RAYS : RAYS {
    $this->N_RAYS=$this->RAYS->rows;
}
weight 0.1;

rule N_RAYS : F_VECTOR {
    $this->N_RAYS=$this->F_VECTOR->[0];
}
precondition : COMBINATORIAL_DIM { $this->COMBINATORIAL_DIM >= 1 }
weight 0.1;

rule N_FACETS : FACETS {
    $this->N_FACETS=$this->FACETS->rows;
}
weight 0.1;

rule N_FACETS : F_VECTOR, COMBINATORIAL_DIM {
   if ($this->COMBINATORIAL_DIM>0) {
      $this->N_FACETS=$this->F_VECTOR->[$this->COMBINATORIAL_DIM-1];
   } else {
      if ($this->COMBINATORIAL_DIM==0) {
         $this->N_FACETS=1; # a vertex or just a ray
      } else {
         $this->N_FACETS=0; # empty polytope or trivial cone
      }
   }
}
weight 0.1;

rule RAY_LABELS : FACET_LABELS, FACETS_THRU_RAYS {
   $this->RAY_LABELS=induced_labels($this->FACET_LABELS, $this->FACETS_THRU_RAYS);
}

# the following rule exists to force the existence of default labels, e.g., to allow meaningful FACET_LABELS

rule RAY_LABELS : N_RAYS {
   my @labels = (0..$this->N_RAYS-1);
   $this->RAY_LABELS="@labels";
}
weight 4.10;

rule FACET_LABELS : RAY_LABELS, RAYS_IN_FACETS {
   $this->FACET_LABELS=induced_labels($this->RAY_LABELS, $this->RAYS_IN_FACETS);
}

# The purpose of the following rule is to have simpler preconditions in subsequent rules.
# The scheduler efficiently handles preconditions which amount to checking a boolean.
rule FULL_DIM : CONE_AMBIENT_DIM, CONE_DIM {
  $this->FULL_DIM=($this->CONE_AMBIENT_DIM==$this->CONE_DIM);
}
weight 0.1;

rule CONE_DIM : CONE_AMBIENT_DIM, LINEAR_SPAN {
   $this->CONE_DIM= $this->CONE_AMBIENT_DIM - $this->LINEAR_SPAN->rows;
}
weight 0.1;

rule LINEAR_SPAN : CONE_AMBIENT_DIM {
   $this->LINEAR_SPAN=new Matrix<Scalar>(0, $this->CONE_AMBIENT_DIM);
}
precondition : FULL_DIM;
weight 0.1;

rule CONE_DIM : RAYS | INPUT_RAYS {
   my $cmp=$this->prepare_computations;
   my $r;
   if (defined (my $l=$this->lookup("LINEALITY_SPACE | INPUT_LINEALITY"))) {
       $r = $this->give("RAYS | INPUT_RAYS") / $l;

   } else {
       $r = $this->give("RAYS | INPUT_RAYS");
   }
   $this->CONE_DIM=rank($r);
}
weight 1.10;

rule COMBINATORIAL_DIM : CONE_DIM, LINEALITY_DIM {
   $this->COMBINATORIAL_DIM = $this->CONE_DIM - $this->LINEALITY_DIM-1;
}
weight 0.10;

rule COMBINATORIAL_DIM : F_VECTOR {
   $this->COMBINATORIAL_DIM = $this->F_VECTOR->dim;
}
weight 0.1;

rule CONE_DIM : COMBINATORIAL_DIM, LINEALITY_DIM {
   $this->CONE_DIM = $this->COMBINATORIAL_DIM + $this->LINEALITY_DIM+1;
}
weight 0.1;

rule LINEALITY_DIM : COMBINATORIAL_DIM, CONE_DIM {
   $this->LINEALITY_DIM = $this->CONE_DIM - $this->COMBINATORIAL_DIM-1;
}
weight 0.1;

rule LINEALITY_DIM : LINEALITY_SPACE {
   $this->LINEALITY_DIM=$this->LINEALITY_SPACE->rows;
}
weight 0.1;

rule POINTED : LINEALITY_DIM {
   $this->POINTED= $this->LINEALITY_DIM == 0;
}
weight 0.1;

rule POINTED : RAYS {
  if (defined (my $l=$this->lookup("INPUT_LINEALITY | LINEALITY_SPACE"))) {
    $this->POINTED = $l->rows==0;
  } else {
    $this->POINTED = true;
  }
}
weight 0.1;


rule POINTED : INPUT_RAYS {
   if (defined (my $l=$this->lookup("INPUT_LINEALITY | LINEALITY_SPACE"))) {
      if ($l->rows>0) { $this->POINTED=false; return; }
   }

   my $v = zero_vector<Scalar>() | (ones_vector<Scalar>() | $this->INPUT_RAYS);
   my $w = new Vector<Scalar>($v->cols);
   $w->[0] = 1;
   $w->[1] = -1;
   my $lp = new LinearProgram<Scalar>(LINEAR_OBJECTIVE => unit_vector<Scalar>($v->cols, 1));
   my $p = new Polytope<Scalar>(INEQUALITIES => -$v/$w, LP => $lp);
   $this->POINTED = $p->LP->MAXIMAL_VALUE > 0;
}
weight 3.10;

rule POINTED : FACETS | INEQUALITIES, CONE_AMBIENT_DIM {
   my $f = $this->give("FACETS | INEQUALITIES");
   my $l = $this->lookup("LINEAR_SPAN | EQUATIONS");
   $this->POINTED = rank(defined($l) && $l->rows > 0 ? $f/$l : $f) == $this->CONE_AMBIENT_DIM;
}

rule cone_only : LINEALITY_SPACE : FACETS | INEQUALITIES {
   my $F;
   my $AH = $this->lookup("LINEAR_SPAN | EQUATIONS");
   if (defined($AH) && $AH->rows > 0) {
      $F = $this->give("FACETS | INEQUALITIES") / $AH;
   } else {
      $F = $this->give("FACETS | INEQUALITIES");
   }
   my $cmp = $this->prepare_computations;
   $this->LINEALITY_SPACE = common::null_space($F);
}
precondition : CONE_AMBIENT_DIM { $this->CONE_AMBIENT_DIM >= 1 }
weight 1.10;

rule LINEALITY_DIM, LINEALITY_SPACE : CONE_AMBIENT_DIM {
   $this->LINEALITY_DIM = 0;
   $this->LINEALITY_SPACE = new Matrix<Scalar>(0, $this->CONE_AMBIENT_DIM);
}
precondition : POINTED;
weight 0.10;

rule ONE_RAY : RAYS {
  $this->ONE_RAY = $this->RAYS->[0];
}
precondition : N_RAYS;
weight 0.10;

rule cone_only : REL_INT_POINT : RAYS {
   $this->REL_INT_POINT=barycenter($this->RAYS);
}
precondition : N_RAYS;
weight 1.10;

rule cone_only : REL_INT_POINT : LINEALITY_DIM, CONE_AMBIENT_DIM {
   $this->REL_INT_POINT=$this->LINEALITY_DIM ? zero_vector<Scalar>($this->CONE_AMBIENT_DIM) : undef;
}
precondition : !N_RAYS;
weight 0.1;

rule N_RAYS : RAYS_IN_FACETS {
   $this->N_RAYS = $this->RAYS_IN_FACETS->cols;
}
weight 0.1;

rule N_FACETS : RAYS_IN_FACETS {
   $this->N_FACETS=$this->RAYS_IN_FACETS->rows;
}
weight 0.1;

rule N_RIDGES : RAYS_IN_RIDGES {
   $this->N_RIDGES=$this->RAYS_IN_RIDGES->rows;
}
weight 0.1;


rule ALTSHULER_DET : RAYS_IN_FACETS {
   $this->ALTSHULER_DET=altshuler_det($this->VERTICES_IN_FACETS);
}


rule GRAPH.NODE_LABELS = RAY_LABELS;

rule DUAL_GRAPH.NODE_LABELS = FACET_LABELS;

rule N_RAY_FACET_INC : RAYS_IN_FACETS {
   my $n=0;
   foreach (@{$this->RAYS_IN_FACETS}) { $n+=@$_ }
   $this->N_RAY_FACET_INC=$n;
}
weight 1.10;

rule SIMPLICIAL : { $this->SIMPLICIAL=1 }
precondition : ESSENTIALLY_GENERIC;
weight 0.1;

rule SIMPLICIAL_CONE : N_RAYS, LINEALITY_DIM, COMBINATORIAL_DIM {
   if ( $this->LINEALITY_DIM > 0 ) {
      $this->SIMPLICIAL_CONE = 0;
   } else {
      $this->SIMPLICIAL_CONE = ( $this->N_RAYS == $this->COMBINATORIAL_DIM+1 );
   }
}
weight 0.1;

rule SIMPLICIAL : COMBINATORIAL_DIM, RAYS_IN_FACETS {
   foreach (@{$this->RAYS_IN_FACETS}) {
      if (@$_ != $this->COMBINATORIAL_DIM) {
         $this->SIMPLICIAL=0;
         return;
      }
   }
   $this->SIMPLICIAL=1;
}
weight 1.10;

rule SIMPLE : COMBINATORIAL_DIM, RAYS_IN_FACETS {
   foreach (@{transpose($this->RAYS_IN_FACETS)}) {
      if (@$_ != $this->COMBINATORIAL_DIM) {
         $this->SIMPLE=0;
         return;
      }
   }
   $this->SIMPLE=1;
}
weight 1.10;

# @category Combinatorics
# Ray degrees of the cone
# @return Vector<Int> - in the same order as [[RAYS]]
user_method VERTEX_DEGREES = GRAPH.NODE_DEGREES;

# @category Combinatorics
# Facet degrees of the polytope.
# The __degree__ of a facet is the number of adjacent facets.
# @return Vector<Int> - in the same order as [[FACETS]]
user_method FACET_DEGREES = DUAL_GRAPH.NODE_DEGREES;

# @category Backward compatibility
# The diameter of the [[GRAPH]] of the cone
# @return Int
user_method DIAMETER = GRAPH.DIAMETER;

# @category Backward compatibility
# The diameter of the [[DUAL_GRAPH]]
# @return Int
user_method DUAL_DIAMETER = DUAL_GRAPH.DIAMETER;

# @category Backward compatibility
# True if the [[GRAPH]] contains no triangle
# @return Bool
user_method TRIANGLE_FREE = GRAPH.TRIANGLE_FREE;

# @category Backward compatibility
# True if the [[DUAL_GRAPH]] contains no triangle
# @return Bool
user_method DUAL_TRIANGLE_FREE = DUAL_GRAPH.TRIANGLE_FREE;

# @category Combinatorics
# True if the [[GRAPH]] is bipartite
# @return Bool
user_method EVEN = GRAPH.BIPARTITE;

# @category Combinatorics
# True if the [[DUAL_GRAPH]] is bipartite
# @return Bool
user_method DUAL_EVEN = DUAL_GRAPH.BIPARTITE;

# @category Topology
# Difference of the black and white nodes if the [[GRAPH]] is [[BIPARTITE]].
# Otherwise -1.
# @return Int
user_method GRAPH_SIGNATURE = GRAPH.SIGNATURE;

# @category Topology
# Difference of the black and white nodes if the [[DUAL_GRAPH]] is [[BIPARTITE]].
# Otherwise -1.
# @return Int
user_method DUAL_GRAPH_SIGNATURE = DUAL_GRAPH.SIGNATURE;

# @category Combinatorics
# Connectivity of the [[GRAPH]]
# this is the minimum number of nodes that have to be removed from the [[GRAPH]] to make it disconnected
# @return Int
user_method CONNECTIVITY = GRAPH.CONNECTIVITY;

# @category Combinatorics
# Connectivity of the [[DUAL_GRAPH]]
# this is the minimum number of nodes that have to be removed from the [[DUAL_GRAPH]] to make it disconnected
# @return Int
user_method DUAL_CONNECTIVITY = DUAL_GRAPH.CONNECTIVITY;


rule GRAPH.CONNECTED : {
  $this->GRAPH->CONNECTED = 1;
}
weight 0.1;

rule DUAL_GRAPH.CONNECTED : {
  $this->DUAL_GRAPH->CONNECTED = 1;
}
weight 0.1;

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

rule F2_VECTOR : HASSE_DIAGRAM.ADJACENCY, HASSE_DIAGRAM.DECORATION, HASSE_DIAGRAM.INVERSE_RANK_MAP, HASSE_DIAGRAM.TOP_NODE, HASSE_DIAGRAM.BOTTOM_NODE, COMBINATORIAL_DIM {
  my $d=$this->COMBINATORIAL_DIM;
  if ( $d < 0 ) {  # polytope empty or point
     $this->F2_VECTOR=undef;
  } else {
   $this->F2_VECTOR=f2_vector($this->HASSE_DIAGRAM);
  }
}
weight 3.10;

# only for COMBINATORIAL_DIM >= 1
# as diagonal requires at least one row/column in F2_VECTOR
rule F_VECTOR : F2_VECTOR {
   $this->F_VECTOR=$this->F2_VECTOR->diagonal;
}
precondition : COMBINATORIAL_DIM { $this->COMBINATORIAL_DIM >= 1 }
weight 0.10;

rule F_VECTOR : HASSE_DIAGRAM.N_NODES, HASSE_DIAGRAM.DECORATION, HASSE_DIAGRAM.INVERSE_RANK_MAP, HASSE_DIAGRAM.TOP_NODE, HASSE_DIAGRAM.BOTTOM_NODE {
  my $n=$this->HASSE_DIAGRAM->N_NODES;
  if ( $n >= 2 ) {
     my $rk = $this->HASSE_DIAGRAM->rank();
     my @f= map { scalar(@{$this->HASSE_DIAGRAM->nodes_of_rank($_)}) } (1 .. $rk-1);
     $this->F_VECTOR= \@f;
  }
}
weight 0.10;

rule F_VECTOR : N_FACETS, N_RAYS, COMBINATORIAL_DIM {
   my $dim = $this->COMBINATORIAL_DIM;
   if ($dim>=0) {
      my $vec = new Vector<Integer>($dim);
      if ($dim>0) {
         $vec->[0] = $this->N_RAYS;
         if ($dim>1) {
            $vec->[$dim-1] = $this->N_FACETS;
         }
         if ($dim==3) {
            $vec->[1] = $vec->[0]+$vec->[2]-2; # euler
         }
      }
      $this->F_VECTOR=$vec;
   }
}
precondition : COMBINATORIAL_DIM { $this->COMBINATORIAL_DIM < 4 }
weight 0.10;

rule N_RIDGES = DUAL_GRAPH.N_EDGES;

rule N_EDGES = GRAPH.N_EDGES;

rule F_VECTOR : N_FACETS, N_RAYS, GRAPH.N_EDGES, DUAL_GRAPH.N_EDGES, COMBINATORIAL_DIM {
    my $dim =  $this->COMBINATORIAL_DIM;
    my $vec = new Vector<Integer>($dim);
    $vec->[0] = $this->N_RAYS;
    $vec->[$dim-1] = $this->N_FACETS;
    $vec->[1] = $this->GRAPH->N_EDGES;
    $vec->[$dim-2] = $this->DUAL_GRAPH->N_EDGES;
    if ($dim==5) {
       $vec->[2] = -$vec->[0]+$vec->[1]+$vec->[3]-$vec->[4]+2 ; # euler
    }
    $this->F_VECTOR=$vec;
}
precondition : COMBINATORIAL_DIM { $this->COMBINATORIAL_DIM == 4 || $this->COMBINATORIAL_DIM == 5 }
weight 0.10;

rule GRAPH.ADJACENCY : RAYS_IN_FACETS {
   $this->GRAPH->ADJACENCY=graph_from_incidence($this->RAYS_IN_FACETS);
}
weight 3.10;

rule DUAL_GRAPH.ADJACENCY : RAYS_IN_FACETS {
   $this->DUAL_GRAPH->ADJACENCY=dual_graph_from_incidence($this->RAYS_IN_FACETS);
}
weight 3.10;

rule HASSE_DIAGRAM.ADJACENCY, HASSE_DIAGRAM.DECORATION, HASSE_DIAGRAM.INVERSE_RANK_MAP, HASSE_DIAGRAM.TOP_NODE, HASSE_DIAGRAM.BOTTOM_NODE : RAYS_IN_FACETS, COMBINATORIAL_DIM {
   $this->HASSE_DIAGRAM=hasse_diagram($this->RAYS_IN_FACETS, $this->COMBINATORIAL_DIM+1);
}
weight 6.20;

rule COMBINATORIAL_DIM : RAYS_IN_FACETS {
   $this->COMBINATORIAL_DIM=dim_from_incidence($this->RAYS_IN_FACETS);
}
weight 3.10;

rule COMBINATORIAL_DIM : HASSE_DIAGRAM.ADJACENCY, HASSE_DIAGRAM.DECORATION, HASSE_DIAGRAM.INVERSE_RANK_MAP, HASSE_DIAGRAM.TOP_NODE, HASSE_DIAGRAM.BOTTOM_NODE {
   $this->COMBINATORIAL_DIM=$this->HASSE_DIAGRAM->dim;
}
weight 0.1;

rule cone_only : 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);
}

rule RAY_SIZES : RAYS_IN_FACETS {
   $this->RAY_SIZES(temporary)=[ map { scalar(@$_) } @{transpose($this->RAYS_IN_FACETS)} ];
}
weight 1.10;

rule FACET_SIZES : RAYS_IN_FACETS {
   $this->FACET_SIZES(temporary)=[ map { scalar(@$_) } @{$this->RAYS_IN_FACETS} ];
}
weight 1.10;

rule POSITIVE : RAYS | INPUT_RAYS {
   foreach my $v (@{$this->give("RAYS | INPUT_RAYS")}) {
      foreach my $x (@$v) {
         if ($x<0) {
            $this->POSITIVE=0;
            return;
         }
      }
   }
   $this->POSITIVE=1;
}
precondition : POINTED;
weight 1.10;

rule LINEAR_SPAN : RAYS | INPUT_RAYS {
   my $cmp=$this->prepare_computations;
   my $r =new Matrix<Scalar>($this->give("RAYS | INPUT_RAYS"));
   if (defined (my $l = $this->lookup("LINEALITY_SPACE | INPUT_LINEALITY"))) {
     $r /= $l;
   }
   $this->LINEAR_SPAN=null_space($r);
}
weight 1.10;

rule RAYS_IN_FACETS : RAYS, FACETS {
   my $cmp=$this->prepare_computations;
   $this->RAYS_IN_FACETS=incidence_matrix($this->FACETS, $this->RAYS);
}

rule INPUT_RAYS_IN_FACETS : INPUT_RAYS, FACETS {
   my $cmp=$this->prepare_computations;
   $this->INPUT_RAYS_IN_FACETS=incidence_matrix($this->FACETS, $this->INPUT_RAYS);
}

rule RAYS_IN_INEQUALITIES : INEQUALITIES, RAYS {
   my $cmp=$this->prepare_computations;
   $this->RAYS_IN_INEQUALITIES=incidence_matrix($this->INEQUALITIES, $this->RAYS);
}

rule RAYS_IN_FACETS, RAYS, LINEALITY_SPACE : INPUT_RAYS_IN_FACETS, INPUT_RAYS {
   my $cmp=$this->prepare_computations;
   compress_incidence_primal($this);
   $this->remove("INPUT_RAYS_IN_FACETS");
}
incurs VertexPerm;

rule RAYS_IN_FACETS, FACETS, LINEAR_SPAN : RAYS_IN_INEQUALITIES, INEQUALITIES {
   my $cmp=$this->prepare_computations;
   compress_incidence_dual($this);
   $this->remove("RAYS_IN_INEQUALITIES");
}
incurs FacetPerm;

rule RAYS_IN_RIDGES : RAYS_IN_FACETS, COMBINATORIAL_DIM {
   if ($this->COMBINATORIAL_DIM >= 2) {
      my $uhd = upper_hasse_diagram($this->RAYS_IN_FACETS, $this->COMBINATORIAL_DIM, $this->COMBINATORIAL_DIM - 2);
      my $nodes = $uhd->nodes_of_rank($this->COMBINATORIAL_DIM - 2);
      my @rir = map {$uhd->FACES->[$_]} (@{$nodes});
      $this->RAYS_IN_RIDGES = new IncidenceMatrix(\@rir);
   } else {
      $this->RAYS_IN_RIDGES = new IncidenceMatrix(0, $this->RAYS_IN_FACETS->cols);
   }
}
precondition : !exists(HASSE_DIAGRAM);

rule RAYS_IN_RIDGES : HASSE_DIAGRAM {
   if ($this->HASSE_DIAGRAM->dim >= 2) {
      my $nodes = $this->HASSE_DIAGRAM->nodes_of_dim($this->HASSE_DIAGRAM->dim - 2);
      my @rir = map {$this->HASSE_DIAGRAM->FACES->[$_]} (@{$nodes});
      $this->RAYS_IN_RIDGES = new IncidenceMatrix(\@rir);
   } else {
      $this->RAYS_IN_RIDGES = new IncidenceMatrix(0, $this->HASSE_DIAGRAM->dim + 1);
   }
}

rule FACETS, LINEAR_SPAN : RAYS, LINEALITY_SPACE, RAYS_IN_FACETS {
   my $cmp=$this->prepare_computations;
   facets_from_incidence($this);
}
precondition : N_RAYS;

rule RAYS, LINEALITY_SPACE : FACETS, LINEAR_SPAN, RAYS_IN_FACETS {
   my $cmp=$this->prepare_computations;
   vertices_from_incidence($this);
}

# helper clients for various visualization tasks

# FIXME check dependence on dim
rule RIF_CYCLIC_NORMAL, NEIGHBOR_FACETS_CYCLIC_NORMAL : CONE_DIM, RAYS, LINEAR_SPAN, RAYS_IN_FACETS, DUAL_GRAPH.ADJACENCY {
   neighbors_cyclic_normal_primal($this);
}
precondition : POINTED;
precondition : CONE_DIM { $this->CONE_DIM>=3 && $this->CONE_DIM<=4 }

# FIXME check dependence on dim
rule RIF_CYCLIC_NORMAL : RAYS_IN_FACETS {
   # types are Array<Array<Int>> and IncidenceMatrix
   $this->RIF_CYCLIC_NORMAL = rows($this->RAYS_IN_FACETS);
}
precondition : COMBINATORIAL_DIM { $this->COMBINATORIAL_DIM<=1 }
weight 1.10;

rule FTR_CYCLIC_NORMAL : FACETS_THRU_RAYS {
   # types are Array<Array<Int>> and IncidenceMatrix
   $this->FTR_CYCLIC_NORMAL = rows($this->FACETS_THRU_RAYS);
}
precondition : COMBINATORIAL_DIM { $this->COMBINATORIAL_DIM<=1 }
weight 1.10;

# FIXME check dependence on dim
rule FTR_CYCLIC_NORMAL, NEIGHBOR_RAYS_CYCLIC_NORMAL : CONE_DIM, FACETS, RAYS_IN_FACETS, GRAPH.ADJACENCY {
   neighbors_cyclic_normal_dual($this);
}
precondition : POINTED;
precondition : CONE_DIM { $this->CONE_DIM>=3 && $this->CONE_DIM<=4 }

# @category Combinatorics
# The interior //d//-dimensional simplices of a cone of combinatorial dimension //d//
# symmetries of the cone are NOT taken into account
rule MAX_INTERIOR_SIMPLICES : COMBINATORIAL_DIM, RAYS, RAYS_IN_FACETS {
    $this->MAX_INTERIOR_SIMPLICES = max_interior_simplices($this);
}

# @category Combinatorics
# The interior //d-1//-dimensional simplices of a cone of combinatorial dimension //d//
# symmetries of the cone are NOT taken into account
rule INTERIOR_RIDGE_SIMPLICES, MAX_BOUNDARY_SIMPLICES : COMBINATORIAL_DIM, RAYS, RAYS_IN_FACETS {
    my $pair=interior_and_boundary_ridges($this);
    $this->INTERIOR_RIDGE_SIMPLICES = $pair->first;
    $this->MAX_BOUNDARY_SIMPLICES = $pair->second;
}

# @category Combinatorics
# A matrix whose rows contain the cocircuit equations of a cone C
# symmetries of the cone are NOT taken into account
rule COCIRCUIT_EQUATIONS : COMBINATORIAL_DIM, RAYS, RAYS_IN_FACETS, INTERIOR_RIDGE_SIMPLICES, MAX_INTERIOR_SIMPLICES {
    $this->COCIRCUIT_EQUATIONS = cocircuit_equations($this->COMBINATORIAL_DIM, $this->RAYS, $this->RAYS_IN_FACETS, $this->INTERIOR_RIDGE_SIMPLICES, $this->MAX_INTERIOR_SIMPLICES);
}

# @category Combinatorics
# A matrix whose rows contain the cocircuit equations of a cone C
# symmetries of the cone are NOT taken into account
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);
}

# @category Combinatorics
# The vertex degrees add up to twice the number of edges.
rule EXCESS_RAY_DEGREE : N_RAYS, GRAPH.N_EDGES, CONE_DIM {
  $this->EXCESS_RAY_DEGREE = 2 * $this->GRAPH->N_EDGES - $this->N_VERTICES * ($this->CONE_DIM-1);
}

# @category Combinatorics
# The vertex degrees add up to twice the number of edges.
rule EXCESS_FACET_DEGREE : N_FACETS, DUAL_GRAPH.N_EDGES, CONE_DIM {
  $this->EXCESS_FACET_DEGREE = 2 * $this->DUAL_GRAPH->N_EDGES - $this->N_FACETS * ($this->CONE_DIM-1);
}

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

# @category Geometry
# returns the geometric dimension of the cone (including the lineality space)
# for the dimension of the pointed part ask for [[COMBINATORIAL_DIM]]
# @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;
}

# @category Geometry
# checks whether a given cone is containeed in another
# @param Cone<Scalar> C_in
# @return Bool

user_method contains(Cone){
  my ($self, $C_in)=@_;
  return cone_contains($C_in, $self);
}

# @category Geometry
# checks whether a given point is contained in a cone
# @param Vector<Scalar> v point
# @return Bool
user_method contains (Vector) {
  my ($self,$v)=@_;
  my $zero = 0;

  # check if $v is zero 
  foreach (@{$v}){
    if ( $_ != $zero ) {
       return cone_contains(new Cone<Scalar>(INPUT_RAYS=>[$v]), $self);
    }
  }

  # 0 is in any cone
  return true;
}


# @category Geometry
# checks whether a given point is contained in the strict interior of a cone
# @param Vector<Scalar> v point
# @return Bool
user_method contains_in_interior {
  my ($self,$v)=@_;
  if (defined (my $f=$self->lookup("FACETS | INEQUALITIES"))) {
    if (defined (my $ah=$self->lookup("LINEAR_SPAN | INPUT_LINEALITY")) ) {
      if ( $ah->rows ) {
        my $b=$ah*$v;
        foreach (@{$b}) { if ( $_ != 0 ) { return false; } }
      }
    }
    my $b=$f*$v;
    foreach (@{$b}) { if ( $_ <= 0 ) { return false; } }
    return true;
  }else {
    return cone_contains_point($self,$v,in_interior=>1);
  }
}

# @category Combinatorics
# The interior //d//-dimensional simplices of a cone of combinatorial dimension //d//
# symmetries of the cone are NOT taken into account
rule MAX_INTERIOR_SIMPLICES : COMBINATORIAL_DIM, RAYS, RAYS_IN_FACETS {
    $this->MAX_INTERIOR_SIMPLICES = max_interior_simplices($this);
}

# @category Combinatorics
# The interior //d-1//-dimensional simplices of a cone of combinatorial dimension //d//
# symmetries of the cone are NOT taken into account
rule INTERIOR_RIDGE_SIMPLICES, MAX_BOUNDARY_SIMPLICES : COMBINATORIAL_DIM, RAYS, RAYS_IN_FACETS {
    my $pair=interior_and_boundary_ridges($this);
    $this->INTERIOR_RIDGE_SIMPLICES = $pair->first;
    $this->MAX_BOUNDARY_SIMPLICES = $pair->second;
}

# @category Combinatorics
# A matrix whose rows contain the cocircuit equations of a cone C
# symmetries of the cone are NOT taken into account
rule COCIRCUIT_EQUATIONS : COMBINATORIAL_DIM, RAYS, RAYS_IN_FACETS, INTERIOR_RIDGE_SIMPLICES, MAX_INTERIOR_SIMPLICES {
    $this->COCIRCUIT_EQUATIONS = cocircuit_equations($this->COMBINATORIAL_DIM, $this->RAYS, $this->RAYS_IN_FACETS, $this->INTERIOR_RIDGE_SIMPLICES, $this->MAX_INTERIOR_SIMPLICES);
}


} # end object Cone

# @category Triangulations, subdivisions and volume
# Create a simplicial complex as a barycentric subdivision of a given cone or polytope.
# Each vertex in the new complex corresponds to a face in the old complex.
# @param Cone c input cone or polytope
# @option Bool no_labels Do not generate [[VERTEX_LABELS]] from the faces of the original cone. default: 0
# @option Bool geometric_realization create a [[topaz::GeometricSimplicialComplex]]; default is true
# @return topaz::SimplicialComplex
user_function barycentric_subdivision<Scalar>(Cone<Scalar> { no_labels=>0, pin_hasse_section=>"HASSE_DIAGRAM", label_section=>"RAY_LABELS", geometric_realization=>1, coord_section=>"RAYS", ignore_top_node=>0 }) {
  topaz::barycentric_subdivision_impl<BasicDecoration,Sequential,Scalar>(@_);
}

# @category Triangulations, subdivisions and volume
# Create a simplicial complex as an iterated barycentric subdivision of a given cone or polytope.
# @param Cone c input cone or polytope
# @param Int n how many times to subdivide
# @option Bool no_labels Do not generate [[VERTEX_LABELS]] from the faces of the original cone. default: 0
# @option Bool geometric_realization create a [[topaz::GeometricSimplicialComplex]]; default is false
# @return topaz::SimplicialComplex
user_function iterated_barycentric_subdivision<Scalar>(Cone<Scalar> $ { no_labels=>0, ignore_top_node=>0, pin_hasse_section=>"HASSE_DIAGRAM", label_section=>"VERTEX_LABELS", pout_section=>"TRIANGULATION.FACETS", geometric_realization=>0, coord_section=>"RAYS" }) {
    topaz::iterated_barycentric_subdivision_impl<BasicDecoration,Sequential,Scalar>(@_);
}

# @category Other
# @param Matrix M
# Create the Lawrence matrix $ Lambda(M) $ corresponding to M.
# If //M// has //n// rows and //d// columns, then Lambda(M) equals
# ( //M//       //I_n// )
# ( //0_{n,d}// //I_n// ).
# @return Matrix
user_function lawrence_matrix<Scalar>(Matrix<Scalar>) {
    my $m = shift;
    my $n = $m->rows();
    my $d = $m->cols();
    my $m1 = $m | unit_matrix<Scalar>($n);
    my $m2 = ones_vector<Scalar>($n) | zero_matrix<Scalar>($n,$d-1) | unit_matrix<Scalar>($n);
    return $m1/$m2;
}


# @category Producing a polytope from polytopes
# Create the Lawrence polytope $ Lambda(P) $ corresponding to P.
# $ Lambda(P) $ has the property that
# $ Gale( Lambda(P) ) = Gale(P) union -Gale(P) $.
# @param Cone P an input cone or polytope
# @return Cone the Lawrence cone or polytope to P
user_function lawrence<Scalar>(Cone<Scalar>) {
    my $p = shift;
    return new Polytope<Scalar>(VERTICES=>lawrence_matrix($p->VERTICES));
}

# @category Triangulations, subdivisions and volume
# Find the maximal interior simplices of a polytope P.
# Symmetries of P are NOT taken into account.
# @param Polytope P the input polytope
# @return Array<Set>
# @example
# > print max_interior_simplices(cube(2));
# | {0 1 2}
# | {0 1 3}
# | {0 2 3}
# | {1 2 3}
user_function max_interior_simplices<Scalar>(Cone<Scalar>) {
    max_interior_simplices_impl<Scalar>(@_);
}


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