#  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
# An object of type VectorConfiguration deals with properties of row vectors,
# assembled into an n x d matrix called [[VECTORS]]. 
# The entries of these row vectors are interpreted as non-homogeneous coordinates. 
# In particular, the coordinates of a VECTOR will *NOT* be normalized to have a leading 1.
# @tparam Scalar default: [[Rational]]
declare object VectorConfiguration<Scalar=Rational>;

INCLUDE
  vector_configuration_properties.rules

object VectorConfiguration {

rule CHIROTOPE : VECTORS {
   $this->CHIROTOPE = chirotope($this->VECTORS);
}
precondition : FULL_DIM;
weight 6.11;

method construct(polytope::Cone) {
  my $cone=$_[1];
  return new VectorConfiguration(VECTORS=>$cone->RAYS/$cone->LINEALITY_SPACE/-$cone->LINEALITY_SPACE);
}

rule MULTIPLE_VECTORS : VECTORS {
   $this->MULTIPLE_VECTORS=detect_multiple($this->VECTORS); 
}
weight 1.10;

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

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

# 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 : VECTOR_AMBIENT_DIM, VECTOR_DIM {
  $this->FULL_DIM=($this->VECTOR_AMBIENT_DIM==$this->VECTOR_DIM);
}
weight 0.1;

rule VECTOR_AMBIENT_DIM : VECTORS {
   $this->VECTOR_AMBIENT_DIM=$this->VECTORS->cols;
}
weight 0.1;

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

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

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


rule POSITIVE : VECTORS {
   foreach my $v (@{$this->VECTORS}) {
      foreach my $x (@$v) {
	 if ($x<0) {
	    $this->POSITIVE=0;
	    return;
	 }
      }
   }
   $this->POSITIVE=1;
}
weight 1.10;

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


}

# @category Combinatorics
# Convert [[CIRCUITS]] or [[COCIRCUITS]] to a 0/+1/-1 matrix, with one row for each circuit/cocircuit, 
# and as many columns as there are VECTORs/POINTS.
# @param Set<Pair<Set<Int>,Set<Int>>> co/circuits a set of circuits or cocircuits
# @return SparseMatrix<Rational>
user_function circuits2matrix(Set<Pair<Set<Int>,Set<Int>>>) {
   my $circuits = shift;
   my $m = $circuits->size();
   my $n = 0;
   my $mat = new SparseMatrix($m, $n+1);
   my $i=0;
   for (my $e=entire($circuits); $e; ++$e, $i++) {
      if ($$e->first->size()) {
	 $n = max($n, $$e->first->back());
      }
      if ($$e->second->size()) {
	 $n = max($n, $$e->second->back());
      }
      $mat->resize($m, $n+1);
      for (my $pos = entire($$e->first); $pos; ++$pos) {
	 $mat->[$i]->[$$pos] = 1;
      }
      for (my $neg = entire($$e->second); $neg; ++$neg) {
	 $mat->[$i]->[$$neg] = -1;
      }
   }
   return $mat;
}


# @category Combinatorics
# Contract a vector configuration //C// along a specified vector //v//.
# @param VectorConfiguration C
# @param Int v index of the vector to contract
user_function contraction(VectorConfiguration $) {
   my ($vc, $index) = @_;
   my $vcc = new Matrix($vc->VECTORS);
   my $vcm = new Matrix($vc->VECTORS->minor([$index], All));
   project_to_orthogonal_complement($vcc, $vcm);
   my @labels = @{$vc->LABELS};
   splice @labels, $index, 1;
   return new VectorConfiguration(VECTORS => $vcc->minor(~[$index],All), LABELS => "@labels");
}

# @category Combinatorics
# Delete a specified vector //v// from a vector configuration //C//.
# @param VectorConfiguration C
# @param Int v index of the vector to delete
user_function deletion(VectorConfiguration $) {
   my ($vc, $index) = @_;
   my @labels = @{$vc->LABELS};
   splice @labels, $index, 1;
   return new VectorConfiguration(VECTORS => $vc->VECTORS->minor(~[$index],All), LABELS => "@labels");
}

# @category Producing a vector configuration
# Orthogonally project a vector configuration to a coordinate subspace.
# 
# The subspace the VectorConfiguration //P// is projected on is given by indices in the set //indices//.
# The option //revert// inverts the coordinate list.
# @tparam Scalar coordinate type
# @param VectorConfiguration P
# @param Array<Int> indices
# @option Bool revert inverts the coordinate list
# @return VectorConfiguration
user_function projection<Scalar>(VectorConfiguration<Scalar>; $=[ ], { revert=>0 }) {
   my ($P, $indices, $options) = @_;
   projection_vectorconfiguration_impl($P, $indices, $options);
}


# @category Producing a vector configuration
# Orthogonally project a vector configuration to a coordinate subspace such that redundant columns are omitted,
# i.e., the projection becomes full-dimensional without changing the combinatorial type.
# @tparam Scalar coordinate type
# @param VectorConfiguration P
# @option Bool no_labels Do not copy [[VERTEX_LABELS]] to the projection. default: 0
# @return VectorConfiguration
user_function project_full<Scalar>(VectorConfiguration<Scalar>; {no_labels=>0}) { 
   projection(@_); 
}

# @category Producing a vector configuration
# Project a vector configuration V to the subspace with a given basis B and express the result in that basis.
# A boolean flag make_affine (by default undef) determines whether the resulting coordinates acquire an extra leading '1'.
# The return type is a VectorConfiguration, unless 
#   (i) P is a PointConfiguration, 
#  (ii) the first column of B is zero,
# (iii) make_affine is not 0,
# in which case it is a PointConfiguration. 
# The return type depends on the input:
# (1) If V is a VectorConfiguration, the result is also a VectorConfiguration.
# (2a) If V is a PointConfiguration and all rows in B start with a 0, the result is a PointConfiguration.
#      If, furthermore, make_affine is undef, it is set to 1.
# (2b) If V is a PointConfiguration and some row of B has a non-zero first entry, the result is a VectorConfiguration.
# The reasoning here is that B having a zero first column or not influences the hyperplane at infinity.
# @tparam Scalar coordinate type
# @param VectorConfiguration V
# @param Matrix B a matrix whose rows contain the basis of the image space
# @option Bool make_affine. If undef (default), apply the above reasoning. If 1 (0), unconditionally (don't) add leading 1's.
# @return VectorConfiguration or PointConfiguration
user_function project_to<Scalar>(VectorConfiguration<Scalar> Matrix<Scalar>; { make_affine=>undef }) { 
    my ($v,$B,$options) = @_;
    my $proj = $v->VECTORS * transpose($B);
    my $make_affine = $options->{make_affine};
    my $first_col_zero = ($B->col(0) == zero_vector<Scalar>($B->rows));
    if ($v->isa("PointConfiguration") && $first_col_zero && $make_affine eq undef) {
	$make_affine = 1;
    }
    if ($make_affine eq 1) {
	$proj = ones_vector<Scalar>() | $proj;
    }
    my $ret;
    if ($v->isa("PointConfiguration") && $first_col_zero) {
	$ret = new PointConfiguration<Scalar>(POINTS=>$proj);
    } else {
	$ret = new VectorConfiguration<Scalar>(VECTORS=>$proj);
    }
    if (defined (my $labels = $v->lookup("LABELS"))) {
	$ret->LABELS = $labels;
    }
    return $ret;
}

# @category Producing a vector configuration
# Project a Polytope or Cone to the subspace with a given basis, and express the result in that basis		       
# A boolean flag make_affine (by default undef) determines whether the resulting coordinates acquire an extra leading '1'.
# The return type is a Cone, unless 
#   (i) P is a Polytope, 
#  (ii) the first column of B is zero,
# (iii) make_affine is not 0,
# in which case it is a Polytope. 
# If make_affine is undef and (ii) is true, make_affine is set to 1.
# The reasoning here is that B having a zero first column or not influences the hyperplane at infinity.
# @tparam Scalar coordinate type
# @param Cone C
# @param Matrix B a matrix whose rows contain the basis of the image space
# @return Cone or Polytope
user_function project_to<Scalar>(Cone<Scalar> Matrix<Scalar>; { make_affine=>undef }) { 
    my ($v,$B, $options) = @_;
    my $proj = $v->give("RAYS | INPUT_RAYS") * transpose($B);
    my $make_affine = $options->{make_affine};
    my $first_col_zero = ($B->col(0) == zero_vector<Scalar>($B->rows));
    if ($make_affine eq undef && $first_col_zero) {
	$make_affine = 1;
    }
    if ($make_affine eq 1) {
	$proj = ones_vector<Scalar>() | $proj;
    }
    my $ret;
    if ($v->isa("Polytope") && $first_col_zero && $make_affine != 0) {
	$ret = new Polytope<Scalar>(POINTS=>$proj);
    } else {
	$ret = new Cone<Scalar>(INPUT_RAYS=>remove_zero_rows($proj));
    }
    if (defined (my $labels = $v->lookup("RAY_LABELS|INPUT_RAY_LABELS"))) {
	$ret->INPUT_RAY_LABELS = $labels;
    }
    return $ret;
}

# @category Producing a vector configuration
# Project a vector configuration V along the subspace with the given basis B.
# The result is still expressed in the original ambient basis.
# If V is a PointConfiguration and the first column of B is zero, the result is a PointConfiguration, else a VectorConfiguration.
# @tparam Scalar coordinate type
# @param VectorConfiguration V
# @param Matrix B a matrix whose rows contain the basis of the space to be projected out
# @return VectorConfiguration
user_function project_out<Scalar>(VectorConfiguration<Scalar> Matrix<Scalar>) { 
    my ($v,$B) = @_;
    my $O = new Matrix<Scalar>($B);
    orthogonalize_subspace($O);
    my $V = new Matrix<Scalar>($v->VECTORS);
    project_to_orthogonal_complement($V, $O);
    my $first_col_zero = ($B->col(0) == zero_vector<Scalar>($B->rows));
    my $ret;
    if ($v->isa("PointConfiguration") && $first_col_zero) {
	$ret = new PointConfiguration<Scalar>(POINTS=>$V);
    } else {
	$ret = new VectorConfiguration<Scalar>(VECTORS=>$V);
    }
    if (defined (my $labels = $v->lookup("LABELS"))) {
	$ret->LABELS = $labels;
    }
    return $ret;
}

# @category Producing a vector configuration
# Project a Cone C along the subspace with the given basis B
# The result is still expressed in the original ambient basis.
# If V is a Polytope and the first column of B is zero, the result is a Polytope, else a Cone.
# @tparam Scalar coordinate type
# @param Cone C
# @param Matrix B a matrix whose rows contain the basis of the space to be projected out
# @return Cone
user_function project_out<Scalar>(Cone<Scalar> Matrix<Scalar>) { 
    my ($v,$B) = @_;
    my $O = new Matrix<Scalar>($B);
    orthogonalize_subspace($O);
    my $V = new Matrix<Scalar>($v->give("RAYS|INPUT_RAYS"));
    project_to_orthogonal_complement($V, $O);
    my $first_col_zero = ($B->col(0) == zero_vector<Scalar>($B->rows));
    my $ret;
    if ($v->isa("Polytope") && $first_col_zero) {
	$ret = new Polytope<Scalar>(POINTS=>$V);
    } else {
	$ret = new Cone<Scalar>(INPUT_RAYS=>$V);
    }
    if (defined (my $labels = $v->lookup("RAY_LABELS|INPUT_RAY_LABELS"))) {
	$ret->INPUT_RAY_LABELS = $labels;
    }
    return $ret;
}
			  

# @category Producing a vector configuration
# Construct a new vector configuration that projects to a given array of vector configurations
# If the n vector configurations are d_1, d_2, ..., d_n-dimensional and all have m vectors,
# the resulting vector configuration is (d_1+...+d_n)-dimensional, has m vectors, and
# the projection to the i-th d_i coordinates gives the i-th input vector configuration.
# @tparam Scalar coordinate type
# @param Array<VectorConfiguration> P_Array
# @return VectorConfiguration
user_function projection_preimage<Scalar>(VectorConfiguration<Scalar> +) {
    my $a = new Array<VectorConfiguration<Scalar>>(@_);
    projection_preimage_impl($a);
}

# @category Producing a vector configuration
# Construct the free sum of two vector configurations.
# 
# @param VectorConfiguration P1
# @param VectorConfiguration 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 VectorConfiguration
user_function free_sum<Scalar>(VectorConfiguration<Scalar> VectorConfiguration<Scalar>; { force_centered=>1, no_coordinates=>0 }) {
   my ($P1, $P2, $options) = @_;
   if (!$P1->isa("PointConfiguration") && $P2->isa("PointConfiguration") ||
       !$P2->isa("PointConfiguration") && $P1->isa("PointConfiguration")) {
       die "free_sum: cannot mix point and vector configurations";
   }
   my $first_coord = ($P1->isa("PointConfiguration") ? 1 : 0);
   free_sum_impl($P1, $P2, "VECTOR", "AFFINE_HULL", $first_coord, $options);
}

# Local Variables:
# mode: perl
# c-basic-offset:4
# End: