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

object SubdivisionOfVectors {


rule LINEAR_SPAN : VECTORS {
   $this->LINEAR_SPAN=null_space($this->VECTORS);
}
weight 1.10;

rule VECTOR_DIM : VECTORS {
   $this->VECTOR_DIM=rank($this->VECTORS);
}
weight 1.10;

rule VECTOR_AMBIENT_DIM : VECTORS {
   $this->VECTOR_AMBIENT_DIM=$this->VECTORS->cols;
}
precondition : VECTORS { $this->VECTORS->rows > 0 }
weight 0.10;

rule VECTOR_DIM : VECTOR_AMBIENT_DIM, LINEAR_SPAN {
   $this->VECTOR_DIM= $this->VECTOR_AMBIENT_DIM - $this->LINEAR_SPAN->rows;
}
weight 0.10;

rule LINEAR_SPAN : {
   $this->LINEAR_SPAN=new Matrix<Scalar>;
}
precondition : FULL_DIM;
weight 0.10;

rule N_VECTORS : VECTORS {
    $this->N_VECTORS=$this->VECTORS->rows;
}
weight 0.10;

rule N_MAXIMAL_CELLS : MAXIMAL_CELLS {
   $this->N_MAXIMAL_CELLS = $this->MAXIMAL_CELLS->rows();
}
weight 0.10;

# @category Geometry
# Parameters for user method [[secondary_cone]].
options %secondary_cone_options=(
    # Matrix system of linear equation the cone is cut with
    equations => undef,
    # Set<Int> restrict lifting function to zero at points designated
    lift_to_zero => undef,
    # Bool restrict lifting functions to zero at the entire face spanned by points designated
    lift_face_to_zero => undef,
    # Bool throws an exception if subdivision is not regular
    test_regularity => undef
    );
    
# @category Geometry
# The secondary cone is the polyhedral cone of all lifting functions on the [[VECTORS]] which induce the subdivision given by the [[MAXIMAL_CELLS]].
# If the subdivision is not regular, the cone will be the secondary cone of the finest regular coarsening.
# @return Cone<Scalar>
user_method secondary_cone(%secondary_cone_options) : VECTORS, MAXIMAL_CELLS {
    my ($self, $options) = @_;
    my $vectors=full_dim_projection($self->VECTORS);
    my $cells=new Array<Set>(rows($self->MAXIMAL_CELLS));

    my $n=$vectors->rows();
    my $rank=rank($vectors);
    if ($cells->size()==1 && $cells->[0]->size()==$n && $rank==$n) {
        return new Cone<Scalar>(RAYS => new Matrix<Scalar>(0,$n),
                                CONE_AMBIENT_DIM => $n,
                                LINEALITY_SPACE => unit_matrix<Scalar>($n));
    }

    my $sc_ineq=secondary_cone_ineq($vectors,$cells,$options);
    my $sc=new Cone<Scalar>(INEQUALITIES=>$sc_ineq->first, EQUATIONS=>$sc_ineq->second);
    if ($options->{test_regularity}) {
        my $w=$sc->REL_INT_POINT;
        my $slack=$sc_ineq->first*$w;
        for (my $i=0; $i<$slack->dim; ++$i) {
            die "subdivision not regular" if $slack->[$i]==0;
        }
    }
    
    return $sc;
}

rule MIN_WEIGHTS : VECTORS, MAXIMAL_CELLS {
    my $n_vectors = $this->VECTORS->rows();
    my $secT = $this->secondary_cone(); # no normalization etc
    my $secTineqs = new Matrix<Rational>(primitive($secT->give("FACETS|INEQUALITIES")));
    my $minus_ones = new Vector<Rational>(-ones_vector<Rational>($secTineqs->rows()));
    my $positive_orthant = zero_vector<Rational>|unit_matrix<Rational>($n_vectors);
    my $ineqs = ($minus_ones|$secTineqs) / $positive_orthant; # shifted ineqs intersected with positive orthant
    my $eqs = zero_vector<Rational> | $secT->give("LINEAR_SPAN|EQUATIONS");
    my $C = new Cone(INEQUALITIES=>$ineqs, EQUATIONS=>$eqs);
    my $inner = primitive($C->REL_INT_POINT);
    my $objective = 0|ones_vector<Rational>($n_vectors);
    my $val = $objective * $inner;
    $ineqs = new Matrix($ineqs / new Vector($val | -ones_vector<Rational>($n_vectors)));
    my $ilp = new polytope::MixedIntegerLinearProgram<Rational>(LINEAR_OBJECTIVE=>$objective);
    my $secTilp = new polytope::Polytope(INEQUALITIES=>$ineqs, EQUATIONS=>$eqs, MILP=>$ilp);
    $this->MIN_WEIGHTS = $secTilp->MILP->MINIMAL_SOLUTION->slice(~[0]);
}
weight 6.10;

} # object SubdivisionOfVectors


object SubdivisionOfPoints {

rule POLYHEDRAL_COMPLEX.VERTICES, POLYHEDRAL_COMPLEX.MAXIMAL_POLYTOPES : POINTS, MAXIMAL_CELLS {
   my $points=$this->POINTS;
   my $n_points=$points->rows();
   my $max_cells=$this->MAXIMAL_CELLS;
   my $vertices=new Set<Vector<Scalar>>;
   # VERTICES are the union of all vertices of MAXIMAL_CELLS
   foreach my $cell (@{$max_cells}) {
      my $v=$points->minor($cell,All);
      my $p=new polytope::Polytope<Scalar>(POINTS=> $v);
      $vertices+=$_ for @{rows($p->VERTICES)};
   }
   # find the non-vertex points
   my $vertex_ord=0;
   my @vertex_indices;
   my $i=-1;
   my @point_map=map {
      ++$i;
      if ($vertices->contains($_)) {
         push @vertex_indices, $i;
         $vertex_ord++;
      } else {
         -1;
      }
   } @{$this->POINTS};
   $this->POLYHEDRAL_COMPLEX->VERTICES=$points->minor(\@vertex_indices,All);
   # re-index MAXIMAL_CELLS
   # requires that MAXIMAL_CELLS are really maximal
   $this->POLYHEDRAL_COMPLEX->MAXIMAL_POLYTOPES=[ map { new Set(grep { $_>=0 } @point_map[@$_]) } @$max_cells ];
}
weight 2.10;

rule POLYHEDRAL_COMPLEX.MAXIMAL_POLYTOPES = MAXIMAL_CELLS;
precondition : CONVEX;

rule POLYHEDRAL_COMPLEX.VERTICES = POINTS;
precondition : CONVEX;

rule REGULAR : {
   $this->REGULAR = 1;
}
precondition : defined(WEIGHTS);

rule MAXIMAL_CELLS : POINTS, WEIGHTS {
   $this->MAXIMAL_CELLS = polytope::regular_subdivision($this->POINTS, $this->WEIGHTS);
}
weight 3.10;
incurs CellPerm;

rule REGULAR, WEIGHTS : POINTS, MAXIMAL_CELLS {
   my $pair = polytope::is_regular($this->POINTS, rows($this->MAXIMAL_CELLS));
   if ($this->REGULAR = $pair->first) {
      $this->WEIGHTS = $pair->second;
   }
}
weight 3.10;

rule UNIMODULAR : POINTS, MAXIMAL_CELLS {
    $this->UNIMODULAR=unimodular( $this->POINTS, rows($this->MAXIMAL_CELLS), true );
}

rule TIGHT_SPAN.HASSE_DIAGRAM.ADJACENCY, TIGHT_SPAN.HASSE_DIAGRAM.DECORATION, TIGHT_SPAN.HASSE_DIAGRAM.INVERSE_RANK_MAP, \
     TIGHT_SPAN.HASSE_DIAGRAM.TOP_NODE, TIGHT_SPAN.HASSE_DIAGRAM.BOTTOM_NODE : \
     POLYHEDRAL_COMPLEX.MAXIMAL_POLYTOPES, POLYHEDRAL_COMPLEX.MAXIMAL_POLYTOPES_INCIDENCES, \
     POLYHEDRAL_COMPLEX.MAXIMAL_POLYTOPES_COMBINATORIAL_DIMS, POLYHEDRAL_COMPLEX.COMBINATORIAL_DIM {
   $this->TIGHT_SPAN->HASSE_DIAGRAM = tight_span_lattice_for_subdivision(
      $this->POLYHEDRAL_COMPLEX->MAXIMAL_POLYTOPES,
      $this->POLYHEDRAL_COMPLEX->MAXIMAL_POLYTOPES_INCIDENCES,
      $this->POLYHEDRAL_COMPLEX->COMBINATORIAL_DIM);
}
weight 6.10;

rule TIGHT_SPAN.VERTICES : MAXIMAL_CELLS, POINTS, WEIGHTS {
   $this->TIGHT_SPAN->VERTICES = tight_span_vertices($this->POINTS, $this->MAXIMAL_CELLS, $this->WEIGHTS);
}
precondition : VECTOR_DIM, VECTOR_AMBIENT_DIM {$this->VECTOR_DIM+1 >= $this->VECTOR_AMBIENT_DIM}
weight 3.10;

rule TIGHT_SPAN.VERTEX_LABELS : TIGHT_SPAN.VERTICES {
   $this->TIGHT_SPAN->VERTEX_LABELS = [0.. $this->TIGHT_SPAN->VERTICES->rows-1];
}

} # object SubdivisionOfPoints

object polytope::Polytope {

property POLYTOPAL_SUBDIVISION {

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

rule POLYTOPAL_SUBDIVISION(any).POINTS = VERTICES;

rule POLYTOPAL_SUBDIVISION.REFINED_SPLITS : SPLITS, VERTICES, POLYTOPAL_SUBDIVISION.MAXIMAL_CELLS {
   $this->POLYTOPAL_SUBDIVISION->REFINED_SPLITS=polytope::splits_in_subdivision($this->VERTICES,$this->POLYTOPAL_SUBDIVISION->MAXIMAL_CELLS,$this->SPLITS);
}

} # object polytope::Polytope

object polytope::PointConfiguration {

rule POLYTOPAL_SUBDIVISION(any).POINTS = POINTS;

rule POLYTOPAL_SUBDIVISION(any).CONVEX = CONVEX;

rule POLYTOPAL_SUBDIVISION.REFINED_SPLITS : SPLITS, POINTS, POLYTOPAL_SUBDIVISION.MAXIMAL_CELLS {
   $this->POLYTOPAL_SUBDIVISION->REFINED_SPLITS=polytope::splits_in_subdivision($this->POINTS,$this->POLYTOPAL_SUBDIVISION->MAXIMAL_CELLS,$this->SPLITS);
}

} # object polytope::PointConfiguration

# @category Triangulations, subdivisions and volume
# Calculate the secondary fan of a point or vector configuration, or polytope.
# @param polytope::VectorConfiguration V (or polytope) the input configuration
# @option Array<Set> initial_subdivision a seed subdivision of //V//
# @option Matrix restrict_to the equations defining a subspace that the secondary fan should be restricted to
# @option Int seed controls the outcome of the random number generator for generating a randomized initial subdivision
# @return PolyhedralFan<Scalar>
user_function secondary_fan<Scalar>(VectorConfiguration<Scalar> { initial_subdivision=>undef, restrict_to=>undef, seed=>undef }) {
    my ($V, $options) = @_;
    return secondary_fan_impl($V->VECTORS, $options);
}

# @category Triangulations, subdivisions and volume
user_function secondary_fan<Scalar>(Cone<Scalar> { initial_subdivision=>undef, restrict_to=>undef, seed=>undef }) {
    my ($V, $options) = @_;
    return secondary_fan_impl($V->RAYS, $options);
}


# @category Geometry
# Calculate the subdivision induced on a point configuration by a height function h.
# The height function is specified as the sum of a set of rows of a matrix.
# Using the RAYS of the secondary_fan of the configuration works well.
# @tparam Scalar the underlying number type
# @param VectorConfiguration<Scalar> pc (or polytope/cone) the input configuration
# @param Matrix<Scalar> R a matrix such that R->cols() == pc->N_VECTORS
# @param Set I (or ARRAY) a set of indices that select rows from R
# @option Bool verbose print the final height function used=? Default 0
# @return Set<Set> the subdivision induced on the configuration by the final height function
user_function induced_subdivision<Scalar>(VectorConfiguration<Scalar>, Matrix<Scalar>, $; { verbose=>0 }) {
    my ($pc, $rays, $cone, $options) = @_;

    if ($cone =~ /^\d+/ || ref($cone) == "ARRAY") {
        $cone = new Set($cone);
    }

    my $d = $pc->VECTORS->cols();
    my $heights = ones_vector<Scalar>($cone->size()) * $rays->minor($cone,All);
    if ($options->{verbose}) {
        print "heights: $heights\n";
    }
    my $p = new Polytope<Scalar>(POINTS=>($pc->VECTORS | $heights));
    my $F = $p->FACETS;
    my $pif = $p->POINTS_IN_FACETS;

    my @subdiv=();
    for (0..$F->rows-1) {
        if ($F->elem($_,$d) > 0) {
            push @subdiv, $pif->[$_];
        }
    }
    return new Set<Set>(\@subdiv);
}

# @category Geometry
# Calculate the subdivision induced on a polytope by a height function h.
user_function induced_subdivision<Scalar>(Cone<Scalar>, Matrix<Scalar>, $; { verbose=>0 }) {
    my ($c, $rays, $cone, $options) = @_;
    return induced_subdivision(new VectorConfiguration<Scalar>(VECTORS=>$c->RAYS), $rays, $cone, $options);
}

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