#  Copyright (c) 1997-2020
#  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.
#-------------------------------------------------------------------------------

CREDIT topcom
  TOPCOM is a package for computing Triangulations Of Point Configurations and Oriented Matroids.
  Copyright by Jörg Rambau.
  http://www.rambau.wm.uni-bayreuth.de/TOPCOM/

# path to the programs from the TOPCOM package
custom $topcom;

CONFIGURE {
   if($mptopcom){
      $topcom = $mptopcom."/points2chiro";
   } else {
      $topcom =~ s{(?<=\S)$}{/points2chiro};
   }
   my $path=find_program($topcom, "points2chiro", { prompt => "the `points2chiro' program from the TOPCOM package" }) or return;
   ($topcom) = $path =~ $directory_re;
}

# @category Triangulation and volume
# Use the [[wiki:external_software#TOPCOM]] package for computing polytope triangulations.

label topcom

#Converts an Array<Array<Int>> of symmetries into a text that TOPCOM understands
sub symmetry_to_text {
   return "[[".join("],[", map { join(",", @{$_}) } @{ $_[0]})."]]";
}

# extracts the symmetry generators of a polymake object into a format that TOPCOM understands
sub topcom_sym_string {
    my $p = shift;
    my $sym_string = "[]";
    if ($p->lookup("GROUP")) {
        if (defined( my $a = $p->isa("PointConfiguration") ? $p->lookup("GROUP.POINTS_ACTION") : $p->lookup("GROUP.RAYS_ACTION") )) {
            $sym_string = symmetry_to_text($a->GENERATORS);
        }
    }
    return $sym_string;
}

# @category Triangulations, subdivisions and volume
# This converts a polytope, cone or point configuration into a format that topcom understands
# @param Cone P (or PointConfiguration)
# @return String
# @example To convert a 2-cube without symmetries into topcom format, type
# > print topcom_input_format(cube(2));
# | [[1,-1,-1],[1,1,-1],[1,-1,1],[1,1,1]]
# | []
# If you also want the symmetry group, you can type
# > print topcom_input_format(cube(2,group=>1));
# | [[1,-1,-1],[1,1,-1],[1,-1,1],[1,1,1]]
# | [[1,0,3,2],[0,2,1,3]]
user_function topcom_input_format($) {
    my $p = shift;
    my $V = $p->isa("PointConfiguration") ? $p->POINTS : $p->RAYS;
    my $topcom_string = "[" . join(",", map { "[".join(",", @{$_})."]" } @{$V}) . "]";
    return $topcom_string . "\n" . topcom_sym_string($p) . "\n";
}    

# @category Triangulations, subdivisions and volume
# This computes the chirotope of a point or vector configuration.
# @param Matrix V The points or vectors, given as rows.
# @option Array<Array<Int>> symmetry A list of generators of a symmetry group of the points, given as permutations on the indices. If specified, the chirotope is only computed up to symmetry.
# @return String
sub chirotope_impl($$) {
    my ($p, $V) = @_;
    my $input = new TempTextFile;
    print $input "[", join(",", map { "[".join(",", @{$_})."]" } @{$V}), "]" . "\n" . topcom_sym_string($p) . "\n";
    my $cmd = "$topcom/points2chiro";
    unless (defined($p->lookup("GROUP"))) {
       $cmd .= " --nosymmetries";
    }
    unless ($DebugLevel) {
       $cmd .= " 2>/dev/null";
    }
    my $result = `$cmd <$input`;
    $result =~ s/^.*\n\Z//m;   # TOPCOM appends the symmetry generators as text, which we forget here.
    return $result;
}

# @category Triangulations, subdivisions and volume
# Compute the chirotope of a polytope using TOPCOM.
# @param Polytope P
# @return String
user_function chirotope<Scalar> [!is_qe(Scalar)] (Polytope<Scalar>) {
    my $p = shift;
    return chirotope_impl($p, $p->VERTICES);
}
                                  
# @category Triangulations, subdivisions and volume
# Compute the chirotope of a point or vector configuration using TOPCOM.
# @param VectorConfiguration P
# @return String
user_function chirotope<Scalar> [!is_qe(Scalar)] (VectorConfiguration<Scalar>) {
    my $p = shift;
    return chirotope_impl($p, $p->VECTORS);
}

# @category Triangulations, subdivisions and volume
# Compute the chirotope of a polytope using polymake's native implementation.
# @param Polytope P
# @return String
user_function chirotope (Polytope<QuadraticExtension>) {
    my $p = shift;
    return chirotope_impl_native($p->VERTICES);
}

# @category Triangulations, subdivisions and volume
# Compute the chirotope of a point or vector configuration using polymake's native implementation.
# @param VectorConfiguration P
# @return String
user_function chirotope (VectorConfiguration<QuadraticExtension>) {
    my $p = shift;
    return chirotope_impl_native($p->VECTORS);
}


sub call_topcom_co_or_circuits {
    my ($p, $executable) = @_;
    my $input = new TempTextFile;
    print $input $p->CHIROTOPE , "\n" , topcom_sym_string($p) , "\n";
    my $cmd = $executable;
    unless ($DebugLevel) {
       $cmd .= " 2>/dev/null";
    }
    open my $output, "$cmd <$input |"
      or die "can't start $executable: $!\n";
    local $_;
    <$output>; # dimensions
    <$output>; # opening bracket
    my @circuits;
    while (<$output>) {
        if ($_ =~ /\{/) {  # not the closing bracket
            chomp;
            my $pair = new Pair<Set<Int>,Set<Int>>;
            my @e = split /\},\{/, $_;
            $_ = $e[0];
            s/\{//g; s/\[//g; s/\]//g; s/\}//g;
            my @x=( split /,/, $_ );
            $pair->first=\@x;

            $_ = $e[1];
            s/\{//g; s/\[//g; s/\]//g; s/\}//g;
            @x = split /,/, $_;
            $pair->second=\@x;

            push @circuits, $pair;
        }
    }
    new Set<Pair<Set,Set>>(\@circuits);
}

sub call_topcom_circuits($)
{
    call_topcom_co_or_circuits($_[0], "$topcom/chiro2circuits");
}

sub call_topcom_cocircuits($)
{
   call_topcom_co_or_circuits($_[0], "$topcom/chiro2cocircuits");
}

sub call_topcom_chiro2placingtriang($) {
   my $input = new TempTextFile;
   print $input $_[0], "\n";
   my $cmd = "$topcom/chiro2placingtriang";
   unless ($DebugLevel) {
      $cmd .= " 2>/dev/null";
   }
   local $_ = `$cmd <$input`;
   chomp;
   s/^\{//; s/\}$//; s/\},\{/}\n{/g; tr/,/ /;
   [ split /\n/, $_ ];
}

sub call_topcom_chiro2triangs {
    my ($p) = @_;
    my $input = new TempTextFile;
    print $input $p->CHIROTOPE, "\n", topcom_sym_string($p), "\n";
    my $cmd = "$topcom/chiro2triangs";
    unless ($DebugLevel) {
       $cmd .= " 2>/dev/null";
    }
    open my $output, "$cmd <$input |"
      or die "can't run $topcom/chiro2triangs: $!\n";
    local $_;
    my @triangs=();
    while (<$output>) {
        s/T.*://; s/\];//; s/^\{//; s/\}$//; s/\},\{/}\n{/g; tr/,/ /;
        my @e = split /\n/, $_;
        push @triangs, \@e;
    }
    \@triangs;
}


sub call_topcom_chiro2finetriangs($) {
    my ($p) = @_;
    my $input = new TempTextFile;
    print $input $p->CHIROTOPE, "\n", topcom_sym_string($p), "\n";
    my $cmd = "$topcom/chiro2finetriangs";
    unless ($DebugLevel) {
       $cmd .= " 2>/dev/null";
    }
    open my $output, "$cmd <$input |"
      or die "can't start $topcom/chiro2finetriangs: $!\n";
    local $_;
    my @triangs;
    while (<$output>) {
       s/T.*://; s/\];//; s/^\{//; s/\}$//; s/\},\{/}\n{/g; tr/,/ /;
       my @e = split /\n/, $_;
       push @triangs, \@e;
    }
    \@triangs;
}

sub call_topcom_chiro2alltriangs($) {
   my $input = new TempTextFile;
   print $input $_[0], "\n";
   my $cmd = "$topcom/chiro2alltriangs";
   unless ($DebugLevel) {
      $cmd .= " 2>/dev/null";
   }
   open my $output, "$cmd <$input |"
     or die "can't start $topcom/chiro2alltriangs: $!\n";
   local $_;
   my @triangs;
   while (<$output>) {
      s/T.*://; s/\];//; s/^\{//; s/\}$//; s/\},\{/}\n{/g; tr/,/ /;
      my @e = split /\n/, $_;
      push @triangs, \@e;
   }
   \@triangs;
}

# @category Triangulations, subdivisions and volume
# Computes all triangulations of a point configuration that are connected by bistellar
# flips to the regular triangulations. The triangulations are computed via the chirotope.
# If the input point configuration or polytope has a symmetry group, only triangulations up to symmetry will be computed.
# @param PointConfiguration pc or Polytope p input point configuration or polytope
# @return Array<Array<Set<Int>>>
user_function topcom_regular_and_connected_triangulations($) {
    my $self = shift;
    return new Array<Array<Set<Int>>>(call_topcom_chiro2triangs($self));
}

# @category Triangulations, subdivisions and volume
# Computes all fine triangulations of a point configuration that are connected by bistellar flips
# to a fine seed triangulation. The triangulations are computed via the chirotope.
# If the input point configuration or polytope has a symmetry group, only fine triangulations up to symmetry will be computed.
# @param PointConfiguration pc or Polytope p input point configuration or polytope
# @return Array<Array<Set<Int>>>
user_function topcom_fine_and_connected_triangulations($){
    my $self = shift;
    return new Array<Array<Set<Int>>>(call_topcom_chiro2finetriangs($self));
}

# @category Triangulations, subdivisions and volume
# Computes all triangulations of a point configuration via its chirotope.
# @param PointConfiguration pc input point configuration
# @return Array<Array<Set<Int>>>
user_function topcom_all_triangulations {
   my $self=shift;
   return new Array<Array<Set<Int>>>(call_topcom_chiro2alltriangs($self->CHIROTOPE));
}

# @category Triangulations, subdivisions and volume
# Computes the GKZ secondary configuration of a point configuration via its chirotope.
# @param PointConfiguration pc input point configuration
# @return PointConfiguration
user_function secondary_polytope<Scalar>(PointConfiguration<Scalar>) {
    my $self=shift;
    my $v = new Matrix<Scalar>( map { gkz_vector<Scalar>($self->POINTS, $_) } @{topcom_all_triangulations($self)} );
    my $vh = ones_vector<Scalar>() | $v;
    return new PointConfiguration<Scalar>(POINTS=>$vh);
}

# @category Triangulations, subdivisions and volume
# Computes the GKZ secondary configuration of a point configuration via its chirotope.
# @param Polytope pc input polytope
# @return PointConfiguration
user_function secondary_polytope<Scalar>(Polytope<Scalar>) {
    my $self=shift;
    my $v = new Matrix<Scalar>( map { gkz_vector<Scalar>($self->VERTICES, $_) } @{topcom_all_triangulations($self)} );
    my $vh = ones_vector<Scalar>() | $v;
    return new PointConfiguration<Scalar>(POINTS=>$vh);
}

# @category Triangulations, subdivisions and volume
# Computes the fiber polytope of a projection of point configurations P->Q via the GKZ secondary configuration.
# @param PointConfiguration pc (or Polytope) source point configuration or polytope
# @param PointConfiguration pc target point configuration
# @return PointConfiguration
user_function fiber_polytope<Scalar>($ PointConfiguration<Scalar>) {
    my ($P, $Q) = @_;
    my $V = $P->isa("PointConfiguration") ? $P->POINTS : $P->VERTICES;
    if ($V->rows() != $Q->N_POINTS) {
        croak("fiber_polytope: The original and target configurations must have the same number of points");
    }
    my $W = new Matrix<Scalar>( map { gkz_vector<Scalar>($Q->POINTS, $_) * $V } @{topcom_all_triangulations($Q)} );
    my $Wh = ones_vector<Scalar>() | $W;
    return new PointConfiguration<Scalar>(POINTS=>$Wh);
}

# @category Triangulations, subdivisions and volume
# Computes the fiber polytope of a projection of point configurations P->Q via the GKZ secondary configuration.
# @param PointConfiguration pc (or Polytope) source point configuration or polytope
# @param Polytope pc target polytope
# @return PointConfiguration
user_function fiber_polytope<Scalar>($ Polytope<Scalar>) {
    my ($P, $Q) = @_;
    my $V = $P->isa("PointConfiguration") ? $P->POINTS : $P->VERTICES;
    if ($V->rows() != $Q->N_VERTICES) {
        croak("fiber_polytope: The original and target configurations must have the same number of points");
    }
    my $W = new Matrix<Scalar>( map { gkz_vector<Scalar>($Q->VERTICES, $_) * $V } @{topcom_all_triangulations($Q)} );
    my $Wh = ones_vector<Scalar>() | $W;
    return new PointConfiguration<Scalar>(POINTS=>$Wh);
}
                                                       
# @category Triangulations, subdivisions and volume
# Computes the fiber polytope of a projection of point configurations P -pi-> Q via the GKZ secondary configuration.
# @param PointConfiguration P (or Polytope) source point configuration or polytope
# @param Matrix pi the projection acting on P
# @return PointConfiguration
user_function fiber_polytope<Scalar>($ Matrix<Scalar>) {
    my ($P, $pi) = @_;
    my $V = $P->isa("PointConfiguration") ? $P->POINTS : $P->VERTICES;
    my $Q = new PointConfiguration<Scalar>(POINTS=>$V * $pi);
    return fiber_polytope($P, $Q);
}

# @category Triangulations, subdivisions and volume
# returns all sets of points that form a
# circuit with the given list of points
# @param Polytope or PointConfiguration P
# @param Set<Int> S subset of point indices
# @return Set<Set<Int>> A list of point sets that intersect positively the set S
user_function positive_circuits {
   my $self=shift;
   my $set=new Set<Int>(num_sorted_uniq(sort {$b <=> $a} @_));
   my $pos_circuits = new Set<Set<Int>>;
   foreach (@{$self->CIRCUITS}) {
     if ( $_->first == $set ) {
	$pos_circuits += $_->second;
     } else {
       if ( $_->second == $set ) {
        $pos_circuits += $_->first;
       }
     }
   }
   return $pos_circuits;
}


object Polytope {

rule CHIROTOPE : VERTICES {
   $this->CHIROTOPE=chirotope($this);
}
weight 6.10;

rule topcom.circuits : CIRCUITS : CHIROTOPE {
   $this->CIRCUITS=call_topcom_circuits($this);
}
weight 6.10;
precondition : FULL_DIM;

rule topcom.cocircuits : COCIRCUITS : CHIROTOPE {
   $this->COCIRCUITS=call_topcom_cocircuits($this);
}
weight 6.10;
precondition : FULL_DIM;

rule topcom.triangulation.poly: TRIANGULATION(new).FACETS : CHIROTOPE {
    $this->TRIANGULATION->FACETS = call_topcom_chiro2placingtriang($this->CHIROTOPE);
}

} # end object Polytope


object VectorConfiguration {

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

rule topcom.circuits : CIRCUITS : CHIROTOPE {
   $this->CIRCUITS=call_topcom_circuits($this);
}
weight 6.10;
precondition : FULL_DIM;

rule topcom.cocircuits : COCIRCUITS : CHIROTOPE {
   $this->COCIRCUITS=call_topcom_cocircuits($this);
}
weight 6.10;
precondition : FULL_DIM;

} # end object VectorConfiguration


object PointConfiguration {

rule topcom.triangulation.pc : TRIANGULATION(new).FACETS : CHIROTOPE {
    $this->TRIANGULATION->FACETS = call_topcom_chiro2placingtriang($this->CHIROTOPE);
}
precondition : FULL_DIM;

} # end object PointConfiguration



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