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

# The [[POINTS]] of an object of type PointConfiguration encode a not necessarily convex finite point set.
# The difference to a parent [[VectorConfiguration]] is that the points have homogeneous coordinates, i.e.
# they will be normalized to have first coordinate 1 without warning.
# @tparam Scalar default: [[Rational]]
declare object PointConfiguration<Scalar=Rational> : VectorConfiguration<Scalar>;

INCLUDE
  point_configuration_properties.rules

object PointConfiguration {


method construct(polytope::Polytope) {
  my $polytope=$_[1];
  if($polytope->LINEALITY_DIM!=0){ die "The polytope has non trivial lineality." }
  return new PointConfiguration(POINTS=>$polytope->VERTICES);
}

rule CONVEX_HULL.POINTS = POINTS;

rule AFFINE_HULL = CONVEX_HULL.AFFINE_HULL;

rule BOUNDED : FAR_POINTS {
   $this->BOUNDED= !@{$this->FAR_POINTS};
}
weight 0.1;


rule BOUNDED : POINTS {
   foreach my $v (@{$this->POINTS}) {
      if (! $v->[0]) {
         $this->BOUNDED=0;
         return;
      }
   }
   $this->BOUNDED=1;
}
weight 1.10;

rule BOUNDED = CONVEX_HULL.BOUNDED;

rule CONVEX_HULL.BOUNDED = BOUNDED;

rule CONVEX : N_POINTS, CONVEX_HULL.N_VERTICES {
   $this->CONVEX= $this->N_POINTS==$this->CONVEX_HULL->N_VERTICES;
}
weight 0.1;


rule CONVEX_HULL.VERTICES = POINTS;
precondition : CONVEX;

rule CONVEX_HULL.VERTEX_POINT_MAP : POINTS, CONVEX_HULL.VERTICES{
    $this->CONVEX_HULL->VERTEX_POINT_MAP=vertex_point_map($this->CONVEX_HULL->VERTICES,$this->POINTS);
}
weight 2.10;

rule VERTEX_POINT_MAP=CONVEX_HULL.VERTEX_POINT_MAP;

rule NON_VERTICES : POINTS, CONVEX_HULL.VERTICES {
   $this->NON_VERTICES = non_vertices($this->POINTS,$this->CONVEX_HULL->VERTICES);
}
weight 2.10;

rule GRAPH.ADJACENCY : CONVEX_HULL.POINTS_IN_FACETS, CONVEX_HULL.FACETS, POINTS, VECTOR_DIM {
   $this->GRAPH->ADJACENCY=points_graph_from_incidence($this->POINTS,$this->CONVEX_HULL->POINTS_IN_FACETS,$this->CONVEX_HULL->FACETS,$this->VECTOR_DIM-1);
}
precondition : VECTOR_DIM { $this->VECTOR_DIM>2 }
weight 3.10;

rule GRAPH.ADJACENCY : N_POINTS, VERTEX_POINT_MAP {
   my $g = new GraphAdjacency($this->N_POINTS);
   $g->edge($this->VERTEX_POINT_MAP->[0],$this->VERTEX_POINT_MAP->[1]);
   $this->GRAPH->ADJACENCY = $g;
}
precondition : VECTOR_DIM { $this->VECTOR_DIM==2 }
weight 1.10;

rule GRAPH.ADJACENCY : N_POINTS {
   my $g = new GraphAdjacency($this->N_POINTS);
   $this->GRAPH->ADJACENCY = $g;
}
precondition : VECTOR_DIM { $this->VECTOR_DIM==1 }
weight 1.10;

rule GRAPH.ADJACENCY = CONVEX_HULL.GRAPH.ADJACENCY;
precondition : CONVEX;

rule FAR_POINTS : POINTS {
   $this->FAR_POINTS=far_points($this->POINTS);
}
weight 1.10;


rule BARYCENTER : POINTS {
   $this->BARYCENTER=barycenter($this->POINTS);
}
precondition : BOUNDED;
weight 1.10;

rule CENTERED = CONVEX_HULL.CENTERED;

rule TRIANGULATION.MASSIVE_GKZ_VECTOR : POINTS, TRIANGULATION.FACETS, CONVEX_HULL.POINTS_IN_FACETS, VECTOR_DIM {
   $this->TRIANGULATION->MASSIVE_GKZ_VECTOR=massive_gkz_vector($this, $this->TRIANGULATION, $this->VECTOR_DIM-1);
}

rule TRIANGULATION.GKZ_VECTOR : POINTS, TRIANGULATION.FACETS {
   $this->TRIANGULATION->GKZ_VECTOR=gkz_vector($this->POINTS,$this->TRIANGULATION->FACETS);
}

rule SPLITS : POINTS, GRAPH.ADJACENCY, CONVEX_HULL.FACETS, VECTOR_DIM {
   $this->SPLITS = splits($this->POINTS,$this->GRAPH->ADJACENCY,$this->CONVEX_HULL->FACETS,$this->VECTOR_DIM-1);
}

rule SPLIT_COMPATIBILITY_GRAPH.ADJACENCY : CONVEX_HULL.FACETS, CONVEX_HULL.AFFINE_HULL, SPLITS {
   $this->SPLIT_COMPATIBILITY_GRAPH->ADJACENCY = split_compatibility_graph($this->SPLITS,$this->CONVEX_HULL);
}

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

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

rule CONVEX_HULL.TRIANGULATION_INT = TRIANGULATION.FACETS;

rule LABELS : N_POINTS {
   my @labels = (0..$this->N_POINTS-1);
   $this->LABELS="@labels";
}
weight 0.10;


sub translate_to_points {
    my ($arr, $vpm) = @_;
    foreach (@$$arr) {
        foreach (@$_) {
            $_ = $vpm->[$_];
        }
    }
} 


#Visualization

rule PIF_CYCLIC_NORMAL : CONVEX_HULL.VIF_CYCLIC_NORMAL, CONVEX_HULL.VERTEX_POINT_MAP {
    my $pif=new Array<Array<Int>>($this->CONVEX_HULL->VIF_CYCLIC_NORMAL);
    translate_to_points(\$pif, $this->CONVEX_HULL->VERTEX_POINT_MAP);
    $this->PIF_CYCLIC_NORMAL(temporary)=$pif;
}

# @category Combinatorics
# Output the faces of a given dimension
# @param PointConfiguration p the input point configuration
# @return Array<Set<Int>>
user_method faces_of_dim($) {
    my ($this, $d) = @_;
    return new Array<Set<Int>>(map { [map {$this->CONVEX_HULL->VERTEX_POINT_MAP->[$_]} @$_] } @{$this->CONVEX_HULL->faces_of_dim($d)})
}


# Triangulation

rule TRIANGULATION(any).VERTEX_LABELS = LABELS;

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

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

rule TRIANGULATION(any).VOLUME = CONVEX_HULL.VOLUME;

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

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

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

#rule TRIANGULATION(any).BOUNDARY.G_DIM = AMBIENT_DIM;

# @category Combinatorics
# Tells the full-dimensional simplices that contain no points except for the vertices.
property MAX_INTERIOR_SIMPLICES : Array<Set>;

# @category Combinatorics
# Tells the number of MAX_INTERIOR_SIMPLICES
property N_MAX_INTERIOR_SIMPLICES : Int;

# @category Combinatorics
# Tells the full-dimensional simplices on the boundary that contain no points except for the vertices.
property MAX_BOUNDARY_SIMPLICES : Array<Set>;

# @category Combinatorics
# Tells the number of MAX_BOUNDARY_SIMPLICES
property N_MAX_BOUNDARY_SIMPLICES : Int;

# @category Combinatorics
# Tells the number of codimension 1 simplices that are not on the boundary
property INTERIOR_RIDGE_SIMPLICES : Array<Set>;

# @category Combinatorics
# Tells the cocircuit equations that hold for the configuration, one for each interior ridge
property COCIRCUIT_EQUATIONS : SparseMatrix;

# @category Combinatorics
# A lower bound for the minimal number of simplices in a triangulation
property SIMPLEXITY_LOWER_BOUND : Int;

# @category Combinatorics
# The full-dimensional simplices that contain no points except for the vertices.
rule MAX_INTERIOR_SIMPLICES : VECTOR_DIM, POINTS, CONVEX_HULL.POINTS_IN_FACETS {
    $this->MAX_INTERIOR_SIMPLICES = max_interior_simplices($this);
}

# @category Combinatorics
# The number of MAX_INTERIOR_SIMPLICES
rule N_MAX_INTERIOR_SIMPLICES : MAX_INTERIOR_SIMPLICES {
    $this->N_MAX_INTERIOR_SIMPLICES = $this->MAX_INTERIOR_SIMPLICES->size();
}


# @category Combinatorics
# The number of MAX_BOUNDARY_SIMPLICES
rule N_MAX_BOUNDARY_SIMPLICES : MAX_BOUNDARY_SIMPLICES {
    $this->N_MAX_BOUNDARY_SIMPLICES = $this->MAX_BOUNDARY_SIMPLICES->size();
}

# @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 : VECTOR_DIM, POINTS, CONVEX_HULL.POINTS_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 C
rule COCIRCUIT_EQUATIONS : VECTOR_DIM, POINTS, CONVEX_HULL.VERTICES_IN_FACETS, INTERIOR_RIDGE_SIMPLICES, MAX_INTERIOR_SIMPLICES {
    $this->COCIRCUIT_EQUATIONS = cocircuit_equations($this->VECTOR_DIM-1, $this->POINTS, $this->CONVEX_HULL->VERTICES_IN_FACETS, $this->INTERIOR_RIDGE_SIMPLICES, $this->MAX_INTERIOR_SIMPLICES);
}

# @category Combinatorics
# A lower bound for the number of simplices needed to triangulate the convex hull
# Symmetries are NOT taken into account
rule SIMPLEXITY_LOWER_BOUND : CONVEX_HULL.VOLUME, MAX_INTERIOR_SIMPLICES, POINTS, VECTOR_DIM, COCIRCUIT_EQUATIONS {
    $this->SIMPLEXITY_LOWER_BOUND = simplexity_lower_bound($this->VECTOR_DIM-1, $this->POINTS, $this->MAX_INTERIOR_SIMPLICES, $this->CONVEX_HULL->VOLUME, $this->COCIRCUIT_EQUATIONS);
}

} # end object PointConfiguration


# A point configuration with an exact coordinate type, like Rational.
declare object_specialization ExactCoord<Scalar> = PointConfiguration<Scalar> [is_ordered_field_with_unlimited_precision(Scalar)] {

rule default.triangulation.pc: TRIANGULATION(new).FACETS : POINTS {
   $this->TRIANGULATION->FACETS=placing_triangulation($this->POINTS);
}
weight 4.10;

rule TRIANGULATION.FACETS : POINTS, TRIANGULATION.WEIGHTS {
   $this->TRIANGULATION->FACETS=regular_subdivision($this->POINTS, $this->TRIANGULATION->WEIGHTS);
}
weight 3.10;
incurs TRIANGULATION.FacetPerm;

rule TRIANGULATION.REGULAR, TRIANGULATION.WEIGHTS : POINTS, TRIANGULATION.FACETS {
    my $pair = is_regular($this->POINTS,$this->TRIANGULATION->FACETS);
    if ($this->TRIANGULATION->REGULAR = $pair->first) {
        $this->TRIANGULATION->WEIGHTS = $pair->second;
    }
}


} # end object_specialization ExactCoord

# @category Producing a point configuration
# Produces the Minkowski sum of //P1// and //P2//.
# @param PointConfiguration P1
# @param PointConfiguration P2
# @return PointConfiguration
# @example
# > $P1 = new PointConfiguration(POINTS=>simplex(2)->VERTICES);
# > $P2 = new PointConfiguration(POINTS=>[[1,1,1],[1,-1,1],[1,1,-1],[1,-1,-1],[1,0,0]]);
# > $m = minkowski_sum($P1,$P2);
# > print $m->POINTS;
# | 1 1 1
# | 1 -1 1
# | 1 1 -1
# | 1 -1 -1
# | 1 0 0
# | 1 2 1
# | 1 0 1
# | 1 2 -1
# | 1 0 -1
# | 1 1 0
# | 1 1 2
# | 1 -1 2
# | 1 1 0
# | 1 -1 0
# | 1 0 1

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

# @category Producing a point configuration
# Produces the polytope //lambda//*//P1//+//mu//*//P2//, where * and + are scalar multiplication
# and Minkowski addition, respectively.
# @param Scalar lambda
# @param PointConfiguration P1
# @param Scalar mu
# @param PointConfiguration P2
# @return PointConfiguration
# @example
# > $P1 = new PointConfiguration(POINTS=>simplex(2)->VERTICES);
# > $P2 = new PointConfiguration(POINTS=>[[1,1,1],[1,-1,1],[1,1,-1],[1,-1,-1],[1,0,0]]);
# > $m = minkowski_sum(1,$P1,3,$P2);
# > print $m->POINTS;
# | 1 3 3
# | 1 -3 3
# | 1 3 -3
# | 1 -3 -3
# | 1 0 0
# | 1 4 3
# | 1 -2 3
# | 1 4 -3
# | 1 -2 -3
# | 1 1 0
# | 1 3 4
# | 1 -3 4
# | 1 3 -2
# | 1 -3 -2
# | 1 0 1

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


# @category Producing a cone
# Computes the normal cone of //p// at a face //F// (or vertex //v//).
# By default this is the inner normal cone.
# @param PointConfiguration p
# @param Set<Int> F (or Int v) point 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 [prefer cdd] [require bundled:cdd] To compute the outer normal cone of a doubled 2-cube, do this:
# > $v = cube(2)->VERTICES;
# > $p = new PointConfiguration(POINTS=>($v/$v));
# > print normal_cone($p, 4, outer=>1)->RAYS;
# | 0 -1
# | -1 0
user_function normal_cone<Scalar>(PointConfiguration<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, "CONVEX_HULL.FACETS_THRU_POINTS", "CONVEX_HULL.POINTS", "CONVEX_HULL.FACETS", $options);
}

# @category Triangulations, subdivisions and volume
# find the maximal interior simplices of a point configuration 
# Symmetries of the configuration are NOT taken into account.
# @param PointConfiguration P the input point configuration
# @return Array<Set>
# @example [prefer cdd] [require bundled:cdd] To calculate the maximal interior simplices of a point configuration, type
# > $p=new PointConfiguration(POINTS=>[[1,0,0],[1,1,0],[1,0,1],[1,1,1],[1,1/2,1/2]]);
# > print max_interior_simplices($p);
# | {0 1 2}
# | {0 1 3}
# | {0 1 4}
# | {0 2 3}
# | {0 2 4}
# | {1 2 3}
# | {1 3 4}
# | {2 3 4}
user_function max_interior_simplices<Scalar>(PointConfiguration<Scalar>) {
    max_interior_simplices_impl<Scalar>(@_);
}


# @category Triangulations, subdivisions and volume
# Calculate the universal polytope of a point configuration A. It is a 0/1 polytope with one vertex for every triangulation of A.
# Each coordinate of the ambient space corresponds to a simplex in the configuration.
# @tparam Scalar the underlying number type
# @param PointConfiguration<Scalar> PC the point configuration
# @return Polytope
# @example [prefer cdd] [require bundled:cdd] To calculate the universal polytope of a point configuration, type
# > $p=new PointConfiguration(POINTS=>[[1,0,0],[1,1,0],[1,0,1],[1,1,1],[1,1/2,1/2]]);
# > print universal_polytope($p)->VERTICES;
# | 1 0 1 0 1 0 0 0 0
# | 1 1 0 0 0 0 1 0 0
# | 1 0 0 1 0 1 0 1 1
# Notice how the first vertex corresponds to the triangulation using all points, and the other ones to the triangulations that don't use the inner point.

user_function universal_polytope<Scalar>(PointConfiguration<Scalar>) {
    my $p = shift;
    return universal_polytope_impl($p->VECTOR_DIM-1, $p->POINTS, $p->MAX_INTERIOR_SIMPLICES, $p->CONVEX_HULL->VOLUME, $p->COCIRCUIT_EQUATIONS);
}

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