#  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 HyperplaneArrangement {

file_suffix hyar

# @category Geometry
# A matrix containing the hyperplanes of the arrangement as rows.
#
# @example The same hyperplane with opposing directions.
# > $HA = new HyperplaneArrangement(HYPERPLANES=>[[1,-1],[-1,1],[1,1]]);
# > print $HA->HYPERPLANES;
# | 1 -1
# | -1 1
# | 1 1
#
# @example A hyperplane that does not cut through the [[SUPPORT]]
# > $HA = new HyperplaneArrangement(HYPERPLANES=>[[1,1]], "SUPPORT.INEQUALITIES"=>unit_matrix(2));
# > print $HA->HYPERPLANES;
# | 1 1
#
property HYPERPLANES = override VECTORS;

# @category Geometry
# A basis of the lineality space of the hyperplane arrangement.

property LINEALITY_SPACE : Matrix<Scalar>;

# @category Input property
# A cone being subdivided by the [[HYPERPLANES]]
# defaults to the whole space.
#
# @example Take the 2-dimensional positive orthant and slice it along the ray through (1,1)
# > $HA = new HyperplaneArrangement(HYPERPLANES=>[[-1,1]], "SUPPORT.INPUT_RAYS"=>[[1,0],[0,1]]);
# > $CD = $HA->CHAMBER_DECOMPOSITION;
# > print $CD->RAYS;
# | 0 1
# | 1 0
# | 1 1
# > print $CD->MAXIMAL_CONES;
# | {1 2}
# | {0 2}
#
# @example [require bundled:cdd] [prefer cdd] Subdivide the two-dimensional space along the axes
# > $HA = new HyperplaneArrangement(HYPERPLANES=>[[1,0],[0,1]]);
# > $CD = $HA->CHAMBER_DECOMPOSITION;
# > print $CD->RAYS;
# | -1 0
# | 0 -1
# | 0 1
# | 1 0
# > print $CD->MAXIMAL_CONES;
# | {2 3}
# | {1 3}
# | {0 2}
# | {0 1}
# > print $CD->COMPLETE;
# | true
property SUPPORT : Cone<Scalar>;

rule SUPPORT.INEQUALITIES : HYPERPLANES {
   $this->SUPPORT->INEQUALITIES = new Matrix<Scalar>(0, $this->HYPERPLANES->cols());
}
precondition : !exists(SUPPORT);

rule LINEALITY_SPACE : HYPERPLANES, SUPPORT.LINEAR_SPAN | SUPPORT.EQUATIONS, SUPPORT.FACETS | SUPPORT.INEQUALITIES {
   my $ineq = $this->SUPPORT->lookup("FACETS | INEQUALITIES");
   my $eq = $this->SUPPORT->lookup("LINEAR_SPAN | EQUATIONS");
   my $hyp = $this->HYPERPLANES;
   $this->LINEALITY_SPACE = null_space($ineq / ($eq / $hyp));
}


# @category Geometry
# Dimension of the space which contains the hyperplane arrangement.
property HYPERPLANE_AMBIENT_DIM = override VECTOR_AMBIENT_DIM;



# @category Geometry
# Slicing the [[SUPPORT]] along every hyperplane of [[HYPERPLANES]] one gets a
# polyhedral fan.
#
# @example Take the 2-dimensional positive orthant and slice it along the ray through (1,1)
# > $HA = new HyperplaneArrangement(HYPERPLANES=>[[-1,1]], "SUPPORT.INPUT_RAYS"=>[[1,0],[0,1]]);
# > $CD = $HA->CHAMBER_DECOMPOSITION;
# > print $CD->RAYS;
# | 0 1
# | 1 0
# | 1 1
# > print $CD->MAXIMAL_CONES;
# | {1 2}
# | {0 2}
property CHAMBER_DECOMPOSITION : PolyhedralFan<Scalar>;

# permuting the [[RAYS]]
permutation ConesPerm : PermBase;


rule ConesPerm.PERMUTATION : ConesPerm.CHAMBER_DECOMPOSITION.MAXIMAL_CONES, CHAMBER_DECOMPOSITION.MAXIMAL_CONES {
   $this->ConesPerm->PERMUTATION=find_permutation(rows($this->ConesPerm->CHAMBER_DECOMPOSITION->MAXIMAL_CONES), rows($this->CHAMBER_DECOMPOSITION->MAXIMAL_CONES));
}


rule CHAMBER_DECOMPOSITION.MAXIMAL_CONES : ConesPerm.CHAMBER_DECOMPOSITION.MAXIMAL_CONES, ConesPerm.PERMUTATION {
   $this->CHAMBER_DECOMPOSITION->MAXIMAL_CONES=permuted_rows($this->ConesPerm->CHAMBER_DECOMPOSITION->MAXIMAL_CONES, $this->ConesPerm->PERMUTATION);
}


# @category Combinatorics
# Incidences between [[CHAMBER_DECOMPOSITION.RAYS]] and [[HYPERPLANES]].
property RAYS_IN_HYPERPLANES : IncidenceMatrix;


# @category Combinatorics
# Number of [[HYPERPLANES]].
#
# @example Coordinate hyperplane arrangement in the plane.
# > $HA = new HyperplaneArrangement(HYPERPLANES=>[[1,0],[0,1]]);
# > print $HA->N_HYPERPLANES;
# | 2
property N_HYPERPLANES = override N_VECTORS;


# @category Combinatorics
# The i-th entry is the signature of the i-th maximal cone of the [[CHAMBER_DECOMPOSITION]] as Set<Int> # of indices of the [[HYPERPLANES]] that evaluate positively on this cone.
#
# @example Take the 2-dimensional positive orthant and slice it along the ray through (1,1)
# > $HA = new HyperplaneArrangement(HYPERPLANES=>[[-1,1]], "SUPPORT.INPUT_RAYS"=>[[1,0],[0,1]]);
# > $CD = $HA->CHAMBER_DECOMPOSITION;
# > print $CD->MAXIMAL_CONES;
# | {1 2}
# | {0 2}
# > print $HA->CHAMBER_SIGNATURES;
# | {}
# | {0}
# > print $HA->chamber_to_signature($CD->MAXIMAL_CONES->[0]);
# | {}
# > print $HA->chamber_to_signature($CD->MAXIMAL_CONES->[1]);
# | {0}
# > print $HA->signature_to_chamber($HA->CHAMBER_SIGNATURES->[0]);
# | 0
# > print $CD->MAXIMAL_CONES->[$HA->signature_to_chamber($HA->CHAMBER_SIGNATURES->[0])];
# | {1 2}
property CHAMBER_SIGNATURES : IncidenceMatrix;

rule CHAMBER_SIGNATURES : ConesPerm.CHAMBER_SIGNATURES, ConesPerm.PERMUTATION {
   $this->CHAMBER_SIGNATURES = permuted_rows($this->ConesPerm->CHAMBER_SIGNATURES,$this->ConesPerm->PERMUTATION);
}


# @category Combinatorics
# Given a maximal cone of [[CHAMBER_DECOMPOSITION]] as Set<Int> containing the
# indices of the rays spanning it, return the signature of the cone as Set<Int>
# of indices of the [[HYPERPLANES]] that evaluate negatively on this cone.
#
# @example Take the 2-dimensional positive orthant and slice it along the ray through (1,1)
# > $HA = new HyperplaneArrangement(HYPERPLANES=>[[-1,1]], "SUPPORT.INPUT_RAYS"=>[[1,0],[0,1]]);
# > $CD = $HA->CHAMBER_DECOMPOSITION;
# > print $CD->MAXIMAL_CONES;
# | {1 2}
# | {0 2}
# > print $HA->CHAMBER_SIGNATURES;
# | {}
# | {0}
# > print $HA->chamber_to_signature($CD->MAXIMAL_CONES->[0]);
# | {}
# > print $HA->chamber_to_signature($CD->MAXIMAL_CONES->[1]);
# | {0}
# > print $HA->signature_to_chamber($HA->CHAMBER_SIGNATURES->[0]);
# | 0
# > print $CD->MAXIMAL_CONES->[$HA->signature_to_chamber($HA->CHAMBER_SIGNATURES->[0])];
# | {1 2}
user_method chamber_to_signature : CHAMBER_DECOMPOSITION.MAXIMAL_CONES, CHAMBER_SIGNATURES {
   my($this, $cell) = @_;
   my $mc = $this->CHAMBER_DECOMPOSITION->MAXIMAL_CONES;
   my $cs = $this->CHAMBER_SIGNATURES;
   for(my $i=0; $i<$mc->rows(); $i++){
      if($mc->row($i) == $cell){
         return $cs->[$i];
      }
   }
   die "$cell could not be found among the maximal cones.";
}

# @category Combinatorics
# Given a signature as a Set<Int> of indices that indicate which
# [[HYPERPLANES]] should evaluate negatively (the remaining evaluate
# positively), return the maximal cone of [[CHAMBER_DECOMPOSITION]] associated to
# this signature. The result the index of the maximal cone in the maximal cones
# of [[CHAMBER_DECOMPOSITION]].
#
# @example Take the 2-dimensional positive orthant and slice it along the ray through (1,1)
# > $HA = new HyperplaneArrangement(HYPERPLANES=>[[-1,1]], "SUPPORT.INPUT_RAYS"=>[[1,0],[0,1]]);
# > $CD = $HA->CHAMBER_DECOMPOSITION;
# > print $CD->MAXIMAL_CONES;
# | {1 2}
# | {0 2}
# > print $HA->CHAMBER_SIGNATURES;
# | {}
# | {0}
# > print $HA->chamber_to_signature($CD->MAXIMAL_CONES->[0]);
# | {}
# > print $HA->chamber_to_signature($CD->MAXIMAL_CONES->[1]);
# | {0}
# > print $HA->signature_to_chamber($HA->CHAMBER_SIGNATURES->[0]);
# | 0
# > print $CD->MAXIMAL_CONES->[$HA->signature_to_chamber($HA->CHAMBER_SIGNATURES->[0])];
# | {1 2}
user_method signature_to_chamber : CHAMBER_DECOMPOSITION.MAXIMAL_CONES, CHAMBER_SIGNATURES {
   my($this, $signature) = @_;
   my $mc = $this->CHAMBER_DECOMPOSITION->MAXIMAL_CONES;
   my $cs = $this->CHAMBER_SIGNATURES;
   for(my $i=0; $i<$cs->rows(); $i++){
      if($cs->[$i] == $signature){
         return $i;
      }
   }
   die "$signature could not be found among the signatures.";
}

# Compute defining properties of the [[CHAMBER_DECOMPOSITION]] using reverse search.
rule CHAMBER_DECOMPOSITION.RAYS, CHAMBER_DECOMPOSITION.MAXIMAL_CONES, CHAMBER_DECOMPOSITION.LINEALITY_SPACE, \
      CHAMBER_SIGNATURES, RAYS_IN_HYPERPLANES : HYPERPLANES, SUPPORT.FACETS | SUPPORT.INEQUALITIES, \
      SUPPORT.LINEAR_SPAN | SUPPORT.EQUATIONS, SUPPORT.RAYS | SUPPORT.INPUT_RAYS, \
      SUPPORT.LINEALITY_SPACE | SUPPORT.INPUT_LINEALITY, SUPPORT.CONE_DIM, LINEALITY_SPACE {
      my($R, $MC, $S, $LS, $RH) = chamber_decomposition_rs($this);
      $this->CHAMBER_DECOMPOSITION->RAYS = $R;
      $this->CHAMBER_DECOMPOSITION->MAXIMAL_CONES = $MC;
      $this->CHAMBER_SIGNATURES = $S;
      $this->CHAMBER_DECOMPOSITION->LINEALITY_SPACE = $LS;
      $this->RAYS_IN_HYPERPLANES = $RH;
}
incurs ConesPerm;

rule RAYS_IN_HYPERPLANES : CHAMBER_DECOMPOSITION.RAYS, HYPERPLANES {
   my $r = $this->CHAMBER_DECOMPOSITION->RAYS;
   my $h = $this->HYPERPLANES;
   my @ia;
   foreach (@{$h*transpose($r)}) {
      my $zeros = new Set<Int>;
      for (my $i=0; $i<$r->rows; $i++) {
         $zeros->collect($i) if ($_->[$i]==0);
      }
      push @ia, $zeros; 
   } 
   $this->RAYS_IN_HYPERPLANES = new IncidenceMatrix<NonSymmetric>([@ia]);
} 

rule CHAMBER_DECOMPOSITION.LINEALITY_SPACE : LINEALITY_SPACE {
   $this->CHAMBER_DECOMPOSITION->LINEALITY_SPACE = $this->LINEALITY_SPACE;
}


# @category Visualization
# Unique names assigned to the [[HYPERPLANES]].

# For a polyhedral fan built from scratch, you should create this property by yourself,
# either manually in a text editor, or with a client program.

property HYPERPLANE_LABELS : Array<String> : mutable;

} #End Object HyperplaneArrangement

# @category Geometry
# This function computes the [[CHAMBER_DECOMPOSITION]] of a given hyperplane
# arrangement in a brute force way, by just considering every possible
# signature. Since not every signature gives a valid cell, it is much cheaper
# to traverse the cells of [[CHAMBER_DECOMPOSITION]] by flipping the walls. This
# method is here for verifying results of our other algorithms.
user_function chamber_decomposition_brute_force<Scalar>(HyperplaneArrangement<Scalar>) {
   my($HA) = @_;
   my $hyperplanes = $HA->HYPERPLANES;
   my $supportIneq = $HA->SUPPORT->lookup("INEQUALITIES | FACETS");
   my $supportEq = $HA->SUPPORT->lookup("EQUATIONS | LINEAR_SPAN");
   my $linSpace = $HA->LINEALITY_SPACE;
   my $joinedEq = $linSpace / $supportEq;
   my $expectedCellDim = $HA->SUPPORT->CONE_DIM - rank($linSpace);

   # my($hyperplanes, $supportIneq, $supportEq, $supportDim) = @_;
   my $numberedVectors = new Map<Vector<Scalar>, Int>();
   my $cells2signatures = new Map<Set<Int>, Set<Int>>();

   for(my $k = 0; $k<= $hyperplanes->rows; $k++){
      my $flippable = all_subsets_of_k(sequence(0, $hyperplanes->rows), $k);
      foreach my $signature (@$flippable) {
         my $cell = transform_to_chamber($hyperplanes, $supportIneq, $joinedEq, $expectedCellDim, $signature, $numberedVectors);
         if($cell->size() > 0){
            $cells2signatures->{$cell} = $signature;
         }
      }
   }

   my @result = canonicalize_ray_order($numberedVectors, $hyperplanes->cols(), $cells2signatures);
   return new HyperplaneArrangement(HYPERPLANES=>$hyperplanes, "SUPPORT.INEQUALITIES"=>$supportIneq, "SUPPORT.EQUATIONS"=>$supportEq, "CHAMBER_DECOMPOSITION.RAYS"=>$result[0], "CHAMBER_DECOMPOSITION.MAXIMAL_CONES"=>$result[1], CHAMBER_SIGNATURES=>$result[2], LINEALITY_SPACE=>$linSpace, "CHAMBER_DECOMPOSITION.LINEALITY_SPACE"=>$linSpace);
}




# Different algorithms might return the rays of the [[CHAMBER_DECOMPOSITION]] in
# different orders. This method sorts the rays and applies the resulting
# permutation to other variables from the algorithm. It then returns an array
# with three elements:
# 1. The sorted rays as a Matrix<Scalar>
# 2. The maximal cones of the [[CHAMBER_DECOMPOSITION]] as IncidenceMatrix
# 3. The signatures of the maximal cones as Array<Set<Int>>.
function canonicalize_ray_order<Scalar>(Map<Vector<Scalar>, Int>, Int, Map<Set<Int>, Set<Int>>) {
   my($numberedVectors, $dim, $cells2signatures) = @_;
   my $n = scalar keys %$numberedVectors;
   my $sortedRays = new Matrix<Scalar>(sort keys %$numberedVectors);
   my $permutation = new Array<Int>(map { $numberedVectors->{$_} } @$sortedRays);

   my $permutedCells2Signatures = group::action($permutation, $cells2signatures);
   my @keys = keys %$permutedCells2Signatures;
   my $I = new IncidenceMatrix(\@keys);
   my $signatures2cells = new Array<Set<Int>>(scalar keys %$cells2signatures);
   for(my $i=0; $i<$I->rows(); $i++){
      $signatures2cells->[$i] = $permutedCells2Signatures->{$I->[$i]};
   }
   return ($sortedRays, $I, $signatures2cells);
}




# Given a signature, assemble the associated cell as a Polytope.
sub assemble_chamber {
   my($hyperplanes, $supportIneq, $supportEq, $signature) = @_;
   my $inequalities = new Matrix($hyperplanes);
   $inequalities->minor(~$signature, All) *= -1;
   return new Cone(INEQUALITIES=>$inequalities/$supportIneq, EQUATIONS=>$supportEq);
}


# Given a cell as a polytope, extract the indices of the rays of the
# [[CHAMBER_DECOMPOSITION]] that span this cone. If some of the rays are new, add
# them and equip them with indices.
sub extract_vertex_indices {
   my($C, $expectedCellDim, $numberedVectors) = @_;
   my $result = new Set<Int>();
   if($C->CONE_DIM == $expectedCellDim){
      my $V = $C->RAYS;
      my $n = scalar keys %$numberedVectors;
      foreach my $vertex (@$V) {
         if(defined $numberedVectors->{$vertex}){
            $result += $numberedVectors->{$vertex};
         } else {
            $result += $n;
            $numberedVectors->{$vertex} = $n++;
         }
      }
   }
   return $result;
}


# Given a signature, get its representation as Set<Int> of indices of the rays
# of [[CHAMBER_DECOMPOSITION]]. This method is running while we have not found all
# rays, so it will remember new rays and equip them with indices.
sub transform_to_chamber {
   my($hyperplanes, $supportIneq, $supportEq, $expectedCellDim, $signature, $numberedVectors) = @_;
   my $C = assemble_chamber($hyperplanes, $supportIneq, $supportEq, $signature);
   return extract_vertex_indices($C, $expectedCellDim, $numberedVectors);
}


user_function make_hyperplanes_unique_in_support<Scalar>(Matrix<Scalar> , Cone<Scalar> ) {
   my($inputHyperplanes, $support) = @_;
   my $ineq = $support->lookup("FACETS | INEQUALITIES");
   my $eq = $support->lookup("LINEAR_SPAN | EQUATIONS");
   my $rays = $support->lookup("RAYS | INPUT_RAYS");
   my $lineality = $support->lookup("LINEALITY_SPACE | INPUT_LINEALITY");
   my $result = new Set<Int>();
   for my $hyp (0..($inputHyperplanes->rows() - 1)) {
      # Check whether we already saved an identical hyperplane (up to sign and multiples).
      my $id = scalar grep(rank($inputHyperplanes->minor(new Set<Int>($_, $hyp), All)) < 2, @$result);
      if($id == 0){
         # Check whether the given hyperplane is a facet of the support.
         $id = scalar grep(rank(new Matrix<Scalar>($inputHyperplanes->row($hyp),$_)) < 2, (@$ineq, @$eq));
         if($id == 0){
            # Check whether the given hyperplane intersects the support non-trivially.
            my($neg, $pos) = (0,0);
            foreach my $ray (@$rays){
               my $check = $inputHyperplanes->row($hyp) * $ray;
               if($check > 0) { $pos++;}
               if($check < 0) { $neg++;}
               if($neg > 0 && $pos > 0){
                  $result += $hyp;
                  last;
               }
            }
            if($neg == 0 || $pos == 0){
               foreach my $lin (@$lineality){
                  if($inputHyperplanes->row($hyp) * $lin != 0){
                     $result += $hyp;
                     last;
                  }
               }
            }
         }
      }
   }
   if($result->size() > 0){
      return $inputHyperplanes->minor($result, All);
   } else {
      return zero_matrix<Scalar>(0, $inputHyperplanes->cols());
   }
}

# @category Producing a hyperplane arrangement
user_function facet_arrangement<Scalar>(Polytope<Scalar>){
   my($P) = @_;
   return new HyperplaneArrangement(HYPERPLANES=>$P->FACETS,HYPERPLANE_LABELS=>$P->FACET_LABELS);
}

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