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

# An affine cone with an exact coordinate type, like Rational.
declare object_specialization ExactCoord<Scalar> = Cone<Scalar> [is_ordered_field_with_unlimited_precision(Scalar)] {

rule TRIVIAL : LINEAR_SPAN, CONE_AMBIENT_DIM {
    $this->TRIVIAL = rank($this->LINEAR_SPAN) == $this->CONE_AMBIENT_DIM;
}	

# A Cone defined with an empty [[RAYS]] or [[INPUT_RAYS]] matrix gets trivial [[FACETS]] assigned
# Note that all convex hull clients need at least one ray.
rule FACETS, LINEAR_SPAN : CONE_AMBIENT_DIM {
   my $ambientdim = $this->CONE_AMBIENT_DIM;
   $this->FACETS = new Matrix<Scalar>(0, $ambientdim);
   my $ls = $this->lookup("LINEALITY_SPACE | INPUT_LINEALITY");
   if (defined($ls) && $ls->rows > 0) {
      $this->LINEAR_SPAN = null_space($ls);
   } else {
      $this->LINEAR_SPAN = unit_matrix<Scalar>($ambientdim);
   }
}
precondition : defined(INPUT_RAYS | RAYS);
precondition : INPUT_RAYS | RAYS {
   $this->give("INPUT_RAYS | RAYS")->rows == 0;
}
weight 1.10;

# @category Convex hull computation
# Use the sequential (beneath-beyond) convex hull algorithm. It performs well at lower dimensions
# and produces a triangulation of the polytope as a byproduct.
label beneath_beyond

rule beneath_beyond.convex_hull.primal, default.triangulation.poly, beneath_beyond.convex_hull: \
     FACETS, LINEAR_SPAN, RAYS_IN_FACETS, DUAL_GRAPH.ADJACENCY, TRIANGULATION(new).FACETS, ESSENTIALLY_GENERIC : RAYS {
   beneath_beyond_find_facets($this, non_redundant => true);
}
weight 4.10;
incurs FacetPerm;

rule beneath_beyond.convex_hull.primal, beneath_beyond.convex_hull: \
     FACETS, RAYS, LINEAR_SPAN, LINEALITY_SPACE, RAYS_IN_FACETS, DUAL_GRAPH.ADJACENCY, TRIANGULATION_INT : INPUT_RAYS {
   beneath_beyond_find_facets($this);
}
weight 4.10;
incurs FacetPerm;

rule beneath_beyond.convex_hull.dual, beneath_beyond.convex_hull: \
     RAYS, LINEALITY_SPACE, RAYS_IN_FACETS, GRAPH.ADJACENCY : FACETS {
   beneath_beyond_find_vertices($this, non_redundant => true);
}
weight 4.10;
incurs VertexPerm;

rule beneath_beyond.convex_hull.dual, beneath_beyond.convex_hull: \
     FACETS, RAYS, LINEAR_SPAN, LINEALITY_SPACE, RAYS_IN_FACETS, GRAPH.ADJACENCY : INEQUALITIES {
   beneath_beyond_find_vertices($this);
}
weight 4.10;
incurs VertexPerm;

rule default.triangulation.poly: TRIANGULATION(new).FACETS : RAYS {
   $this->TRIANGULATION->FACETS=placing_triangulation($this->RAYS, non_redundant => true);
}
weight 4.10;

rule LINEALITY_SPACE : INPUT_RAYS, CONE_AMBIENT_DIM {
   $this->LINEALITY_SPACE = lineality_via_lp<Scalar>($this->INPUT_RAYS, $this->lookup("INPUT_LINEALITY") // new Matrix<Scalar>(0,$this->CONE_AMBIENT_DIM));
}
weight 4.5;

rule LINEAR_SPAN : INEQUALITIES, CONE_AMBIENT_DIM {
   $this->LINEAR_SPAN = lineality_via_lp<Scalar>($this->INEQUALITIES, $this->lookup("EQUATIONS") // new Matrix<Scalar>(0,$this->CONE_AMBIENT_DIM));
}
weight 4.5;

}

object Polytope {

label jarvis

rule jarvis.convex_hull.primal, jarvis.convex_hull: VERTICES : POINTS {
  $this->VERTICES=jarvis($this->POINTS);
}
precondition : BOUNDED, CONE_AMBIENT_DIM { $this->BOUNDED and $this->CONE_AMBIENT_DIM==3 }
weight 1.50;
incurs VertexPerm;


rule VERTICES, AFFINE_HULL : ZONOTOPE_INPUT_POINTS, CENTERED_ZONOTOPE {
   my $m = zonotope_vertices_fukuda($this->ZONOTOPE_INPUT_POINTS, centered_zonotope => $this->CENTERED_ZONOTOPE);
   $this->VERTICES = $m;
   $this->AFFINE_HULL = null_space($m);
}
weight 2.50;
incurs VertexPerm;

rule FEASIBLE, CONE_AMBIENT_DIM : ZONOTOPE_INPUT_POINTS {
   $this->FEASIBLE=$this->ZONOTOPE_INPUT_POINTS->rows > 0;
   $this->CONE_AMBIENT_DIM=$this->ZONOTOPE_INPUT_POINTS->cols;
}
weight 0.10;

rule TILING_LATTICE : VERTICES, VERTICES_IN_FACETS, VERTEX_BARYCENTER, ZONOTOPE_INPUT_POINTS {
   $this->TILING_LATTICE = zonotope_tiling_lattice($this);
}
weight 2.50;

rule SIMPLEXITY_LOWER_BOUND : COMBINATORIAL_DIM, VERTICES, MAX_INTERIOR_SIMPLICES, VOLUME, COCIRCUIT_EQUATIONS {
   $this->SIMPLEXITY_LOWER_BOUND = simplexity_lower_bound($this->COMBINATORIAL_DIM, $this->VERTICES, $this->MAX_INTERIOR_SIMPLICES, $this->VOLUME, $this->COCIRCUIT_EQUATIONS);
}


}  # /Polytope



# @category Producing a polytope from scratch
# Create a zonotope from a matrix whose rows are input points or vectors.
#
# This method merely defines a Polytope object with the property
# [[ZONOTOPE_INPUT_POINTS]].
# @param Matrix<Scalar> M input points or vectors
# @option Bool rows_are_points true if M are points instead of vectors; default true
# @option Bool centered true if output should be centered; default true
# @return Polytope<Scalar> the zonotope generated by the input points or vectors
# @example [nocompare]
# The following produces a parallelogram with the origin as its vertex barycenter:
# > $M = new Matrix([[1,1,0],[1,1,1]]);
# > $p = zonotope($M);
# > print $p->VERTICES;
# | 1 0 -1/2
# | 1 0 1/2
# | 1 -1 -1/2
# | 1 1 1/2
# @example [nocompare]
# The following produces a parallelogram with the origin being a vertex (not centered case):
# > $M = new Matrix([[1,1,0],[1,1,1]]);
# > $p = zonotope($M,centered=>0);
# > print $p->VERTICES;
# | 1 1 0
# | 1 0 0
# | 1 1 1
# | 1 2 1
user_function zonotope<Scalar> (Matrix<Scalar>, { rows_are_points => 1, centered => 1 }) {
   my ($M, $options) = @_;
   my $z = new Polytope<Scalar>(ZONOTOPE_INPUT_POINTS => ($options->{"rows_are_points"} ? $M : ones_vector<Scalar>() | $M),
				CENTERED_ZONOTOPE => $options->{"centered"});
   $z->description = "Zonotope generated by input " . ($options->{"rows_are_points"} ? "points" : "vectors");
   return $z;
}


# @category Producing a polytope from polytopes
# Orthogonally project a pointed polyhedron to a coordinate subspace.
#
# The subspace the polyhedron //P// is projected on is given by indices in the set //indices//.
# The option //revert// inverts the coordinate list.
# The client scans for all coordinate sections and produces proper output from each.
# If a description in terms of inequalities is found, the client performs Fourier-Motzkin elimination
# unless the //nofm// option is set.  Setting the //nofm// option is useful if the corank of the projection
# is large; in this case the number of inequalities produced grows quickly.
# @param Cone P
# @param Array<Int> indices
# @option Bool revert inverts the coordinate list
# @option Bool nofm suppresses Fourier-Motzkin elimination
# @return Cone
# @example [prefer cdd] [require bundled:cdd] project the 3-cube along the first coordinate, i.e. to the subspace
# spanned by the second and third coordinate:
# > $p = projection(cube(3),[1],revert=>1);
# > print $p->VERTICES;
# | 1 1 -1
# | 1 1 1
# | 1 -1 1
# | 1 -1 -1

user_function projection<Scalar>(Cone<Scalar>; $=[ ], { revert=>0, nofm=>0 }) {
   my ($P, $indices, $options) = @_;
   projection_cone_impl<Scalar>($P, $indices, $options);
}


# @category Producing a polytope from polytopes
# Orthogonally project a polyhedron to a coordinate subspace such that redundant columns are omitted,
# i.e., the projection becomes full-dimensional without changing the combinatorial type.
# The client scans for all coordinate sections and produces proper output from each.
# If a description in terms of inequalities is found, the client performs Fourier-Motzkin elimination
# unless the //nofm// option is set.  Setting the //nofm// option is useful if the corank of the projection
# is large; in this case the number of inequalities produced grows quickly.
# @param Cone P
# @option Bool nofm suppresses Fourier-Motzkin elimination
# @option Bool no_labels Do not copy [[VERTEX_LABELS]] to the projection. default: 0
# @return Cone
user_function project_full<Scalar>(Cone<Scalar>; {nofm => 0, no_labels=>0}) {
   projection(@_);
}

# @category Producing a polytope from polytopes
# Construct a new polyhedron that projects to a given array of polyhedra.
# If the n polyhedra are d_1, d_2, ..., d_n-dimensional and all have m vertices,
# the resulting polyhedron is (d_1+...+d_n)-dimensional, has m vertices, and
# the projection to the i-th d_i coordinates gives the i-th input polyhedron.
# @param Array<Cone> P_Array
# @return Cone
# @example
# > $p = projection_preimage(cube(2),cube(2));
# > print $p->VERTICES;
# | 1 -1 -1 -1 -1
# | 1 1 -1 1 -1
# | 1 -1 1 -1 1
# | 1 1 1 1 1
user_function projection_preimage<Scalar>(Cone<Scalar> +) {
    my $a = new Array<Cone<Scalar>>(@_);
    projection_preimage_impl($a);
}


# @category Producing a polytope from polytopes
# Construct a new polyhedron as the free sum of two given bounded ones.
# @param Polytope P1
# @param Polytope P2
# @option Bool force_centered if the input polytopes must be centered. Defaults to true.
# @option Bool no_coordinates produces a pure combinatorial description. Defaults to false.
# @return Polytope
# @example
# > $p = free_sum(cube(2),cube(2));
# > print $p->VERTICES;
# | 1 -1 -1 0 0
# | 1 1 -1 0 0
# | 1 -1 1 0 0
# | 1 1 1 0 0
# | 1 0 0 -1 -1
# | 1 0 0 1 -1
# | 1 0 0 -1 1
# | 1 0 0 1 1
user_function free_sum<Scalar>(Cone<Scalar> Cone<Scalar>; { force_centered=>1, no_coordinates=>0 }) {
   my ($P1, $P2, $options) = @_;
   if (!$P1->isa("Polytope") && $P2->isa("Polytope") ||
       !$P2->isa("Polytope") && $P1->isa("Polytope")) {
       die "free_sum: cannot mix cones and polytopes";
   }
   my $first_coord = ($P1->isa("Polytope") ? 1 : 0);
   free_sum_impl($P1, $P2, "CONE", "LINEAR_SPAN", $first_coord, $options);
}

# @category Producing a polytope from polytopes
# Decompose a given polytope into the free sum of smaller ones
# @param Polytope P
# @return Array<Polytope>
user_function free_sum_decomposition<Scalar>(Polytope<Scalar>) {
    my ($p) = @_;
    my $indices = free_sum_decomposition_indices($p);
    my $summands = new Array<Polytope<Scalar>>($indices->size());
    foreach my $i (0..$indices->size()-1) {
	my $q = new Polytope<Scalar>(VERTICES=>$p->VERTICES->minor($indices->[$i], All), N_VERTICES=>$indices->[$i]->size(), CENTERED=>1);
	$summands->[$i] = $q;
    }
    return $summands;
}


# @category Producing a cone
# Computes the normal cone of //p// at a face //F// (or a vertex //v//).
# By default this is the inner normal cone.
# @param Cone p
# @param Set<Int> F (or Int v) vertex indices which are not contained in the far face
# @option Bool outer Calculate outer normal cone?  Default value is 0 (= inner)
# @option Bool attach Attach the cone to //F//? Default 0 (ie, return the cone inside the hyperplane at infinity)
# @return Cone
# @example To compute the outer normal cone at a vertex of the 3-cube, do this:
# > $c = normal_cone(cube(3), 0, outer=>1);
# > print $c->RAYS;
# | -1 0 0
# | 0 -1 0
# | 0 0 -1
# @example To compute the outer normal cone along an edge of the 3-cube, do this:
# > print normal_cone(cube(3), [0,1], outer=>1)->RAYS;
# | 0 -1 0
# | 0 0 -1
# @example If you want to attach the cone to the polytope, specify the corresponding option:
# > print normal_cone(cube(3), [0,1], outer=>1, attach=>1)->RAYS;
# | 1 -1 -1 -1
# | 1 1 -1 -1
# | 0 0 -1 0
# | 0 0 0 -1
user_function normal_cone<Scalar>(Cone<Scalar> $; { outer => 0, attach => 0 }) {
    my ($c, $F, $options) = @_;
    if ($F =~ /^\d+/ || ref($F) == "ARRAY") {
	$F = new Set($F);
    }
    return normal_cone_impl<Scalar>($c, $F, "FACETS_THRU_RAYS", "RAYS", "FACETS", $options);
}

# @category Producing a cone
# Computes the inner cone of //p// at a face //F// (or a vertex //v//).
# @param Cone p
# @param Set<Int> F (or Int v) vertex indices which are not contained in the far face
# @option Bool outer Make it point outside the polytope? Default value is 0 (= point inside)
# @option Bool attach Attach the cone to //F//? Default 0 (ie, return the cone inside the hyperplane at infinity)
# @return Cone
# @example To compute the inner cone at a vertex of the 3-cube, do this:
# > $c = inner_cone(cube(3), 1);
# > print $c->RAYS;
# | -1 0 0
# | 0 1 0
# | 0 0 1
# @example [nocompare] To compute the inner cone along an edge of the 3-cube, and make it point outside the polytope, do this:
# > print inner_cone(cube(3), [0,1], outer=>1)->RAYS;
# | 0 0 -1
# | 0 -1 0
# @example If you want to attach the cone to the polytope, specify the corresponding option:
# > print normal_cone(cube(3), [0,1], attach=>1)->RAYS;
# | 1 -1 -1 -1
# | 1 1 -1 -1
# | 0 0 1 0
# | 0 0 0 1
user_function inner_cone<Scalar>(Cone<Scalar> $; { outer => 0, attach => 0 }) {
    my ($c, $F, $options) = @_;
    if ($F =~ /^\d+/ || ref($F) == "ARRAY") {
	$F = new Set($F);
    }
    return inner_cone_impl<Scalar>($c, $F, $options);
}


# @category Geometry
# For a face //F// of a cone or polytope //P//, return the polyhedral cone //C// such that
# taking the convex hull of //P// and any point in //C// destroys the face //F//
# @param Cone P
# @param Set F
# @return Cone
# @example To find the occluding cone of an edge of the 3-cube, type
# > $c=occluding_cone(cube(3), [0,1]);
# > print $c->FACETS;
# | -1 0 -1 0
# | -1 0 0 -1
user_function occluding_cone<Scalar>(Cone<Scalar> $) {
    my ($c, $face, $options) = @_;
    if ($face =~ /^\d+/ || ref($face) == "ARRAY") {
	$face = new Set($face);
    }
    my $F   = $c->FACETS;
    my $ls  = $c->LINEAR_SPAN;
    my $vif = $c->RAYS_IN_FACETS;

    my @occluding_facets;
    for (0..$F->rows-1) {
	if (incl($face, $vif->[$_]) <= 0) {
	    push @occluding_facets, $_;
	}
    }
    if (scalar @occluding_facets == 0) {
	croak("The set $face does not index a face of the input polytope or cone");
    }

    my $ineqs = -$F->minor(\@occluding_facets, All);
    return new Cone<Scalar>(INEQUALITIES=>$ineqs, EQUATIONS=>$ls, CONE_AMBIENT_DIM=>$F->cols);
}

# @category Combinatorics
# Calculate the codegree of a cone or polytope P.
# This is the maximal positive integer c such that every subset of size < c lies in a common facet of conv P.
# Moreover, the relation degree(P) + codegree(P) = dim(P) + 1 holds.
# @param Cone P
# @tparam Scalar the underlying number type,
# @example To find the codegree of the 3-cube, type
# > print codegree(cube(3));
# | 1
user_function codegree<Scalar>(Cone<Scalar>) {
    my $p = shift;
    return codegree_impl($p->COMBINATORIAL_DIM, $p->RAYS_IN_FACETS);
}

# @category Combinatorics
# Calculate the codegree of a point configuration P.
# This is the maximal positive integer c such that every subset of size < c lies in a common facet of conv P.
# Moreover, the relation degree(P) + codegree(P) = dim(P) + 1 holds.
# @param PointConfiguration P
# @tparam Scalar the underlying number type,
user_function codegree<Scalar>(PointConfiguration<Scalar>) {
    my $p = shift;
    return codegree_impl($p->CONVEX_HULL->COMBINATORIAL_DIM, $p->CONVEX_HULL->POINTS_IN_FACETS);
}

# @category Combinatorics
# Calculate the degree of a cone, polytope or point configuration P.
# This is the maximal dimension of an interior face of P,
# where an interior face is a subset of the points of P whose convex hull does not lie on the boundary of P.
# Moreover, the relation degree(P) + codegree(P) = dim(P) + 1 holds.
# @param PointConfiguration P (or Cone or Polytope)
# @tparam Scalar the underlying number type,
# @example To find the degree of the 3-cube, type
# > print degree(cube(3));
# | 3
user_function degree($) {
    my $p = shift;
    my $d   = $p->isa("PointConfiguration") ? $p->CONVEX_HULL->COMBINATORIAL_DIM : $p->COMBINATORIAL_DIM;
    my $pif = $p->isa("PointConfiguration") ? $p->CONVEX_HULL->POINTS_IN_FACETS : $p->RAYS_IN_FACETS;
    return $d + 1 - codegree_impl($d, $pif);
}


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