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

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 {
   eval {
      if($mptopcom && ! -e "$topcom/points2chiro"){
         $topcom = $mptopcom;
      }
   };
   $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, $ignoresym) = @_;
    my $sym_string = "[]";
    if ($p->lookup("GROUP") && !$ignoresym) {
        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 topcom_chirotope($$) {
    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";
    }
    my $error = new TempTextFile;
    my $result = `$cmd < $input 2>$error`;
    if($?){
        my $errormsg = `cat $error`;
        die "Could not run topcom: $errormsg\n";
    }

    $result =~ s/^.*\n\Z//m;   # TOPCOM appends the symmetry generators as text, which we forget here.
    return $result;
}


sub call_topcom_co_or_circuits {
    my ($p, $executable, $ignoresym) = @_;
    my $input = new TempTextFile;
    print $input $p->CHIROTOPE , "\n" , topcom_sym_string($p, $ignoresym) , "\n";
    my $cmd = $executable;
    unless ($DebugLevel) {
       $cmd .= " 2>/dev/null";
    }
    open my $output, "$cmd <$input |"
      or die "can't start $executable: $!\n";
    local $_;
    $_ = <$output>;
    $_ = <$output> if $_ =~ /^\d+,\d+/; # skip dimension
    $_ = <$output> if $_ =~ /^\{/; # opening brace
    my @circuits;
    do {
        if ($_ =~ /\{/) {  # not the closing bracket
            chomp;
            my $pair = new Pair<Set<Int>,Set<Int>>;
            my @e = split /\},\{/, $_;
            $_ = $e[0];
            s/C.*:=? ?//; 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;
        }
    } while (<$output>);
    new Set<Pair<Set,Set>>(\@circuits);
}

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

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

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_on_chiro {
   my($p, $cmd, $opt) = @_;
   my $input = new TempTextFile;
   print $input $p->CHIROTOPE, "\n", topcom_sym_string($p), "\n";
   return compute_triangs_with_topcom($cmd, $input, $opt);
}


sub call_topcom_on_points {
   my($p, $cmd, $opt) = @_;
   my $input = new TempTextFile;
   print $input topcom_input_format($p), "\n";
   return compute_triangs_with_topcom($cmd, $input, $opt);
}



sub call_topcom_chiro2alltriangs($) {
   my $input = new TempTextFile;
   print $input $_[0], "\n";
   return compute_triangs_with_topcom("chiro2alltriangs", $input, "");
}

sub compute_triangs_with_topcom {
   my($cmd, $input, $opt) = @_;
   my $executable = "$topcom/$cmd";
   $executable .= " $opt";
   my $errorex = $executable;
   unless ($DebugLevel) {
      $executable .= " 2>/dev/null";
   }
   open my $output, "$executable < $input |"
     or die "can't start $errorex: $!\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 regular triangulations of a point configuration.
#
# @param PointConfiguration pc or Polytope p input point configuration or polytope
# @return Array<Set<Set<Int>>>
user_function topcom_regular_triangulations($) {
    my $self = shift;
    return new Array<Set<Set<Int>>>(call_topcom_on_points($self, "points2triangs", "--regular"));
}

# @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<Set<Set<Int>>>
user_function topcom_regular_and_connected_triangulations($) {
    my $self = shift;
    return new Array<Set<Set<Int>>>(call_topcom_on_chiro($self, "chiro2triangs", ""));
}

# @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<Set<Set<Int>>>
user_function topcom_fine_and_connected_triangulations($){
    my $self = shift;
    return new Array<Set<Set<Int>>>(call_topcom_on_chiro($self, "chiro2finetriangs", ""));
}

# @category Triangulations, subdivisions and volume
# Computes all fine triangulations (sometimes called “full”) of a chirotope.
# 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<Set<Set<Int>>>
user_function topcom_fine_triangulations($){
    my $self = shift;
    return new Array<Set<Set<Int>>>(call_topcom_on_chiro($self, "chiro2allfinetriangs", ""));
}

# @category Triangulations, subdivisions and volume
# Computes all fine and regular triangulations of a point configuration.
# @param PointConfiguration pc or Polytope p input point configuration or polytope
# @return Array<Set<Set<Int>>>
user_function topcom_fine_and_regular_triangulations($) {
    my $self = shift;
    return new Array<Set<Set<Int>>>(call_topcom_on_points($self, "points2finetriangs", "--regular"));
}

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

# @category Triangulations, subdivisions and volume
# Computes the point configuration of GKZ vectors of a point configuration
# via its chirotope using topcom or mptopcom.
#
# @param PointConfiguration pc input point configuration
# @return PointConfiguration
#
# @example The following [[PointConfiguration]] is called the "mother of all
# examples (moae)". It has two non-regular triangulations, which can be seen
# when comparing the number of points of the output configuration with the
# number of vertices of the convex hull of the output configuration.
# > $moae = new PointConfiguration(POINTS=>[[1,4,0,0],[1,0,4,0],[1,0,0,4],[1,2,1,1],[1,1,2,1],[1,1,1,2]]);
# > $moae = project_full($moae);
# > $SC_moae = secondary_configuration($moae);
# > print $SC_moae -> N_POINTS;
# | 18
# > print $SC_moae -> CONVEX_HULL -> N_VERTICES;
# | 16
user_function secondary_configuration<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 point configuration of GKZ vectors of a point configuration
# via its chirotope using topcom or mptopcom.
#
# @param Polytope pc input polytope
# @return PointConfiguration
#
# @example The square only has two triangulations using its vertices.
# > $square = cube(2,1,0);
# > $SC_square = secondary_configuration($square);
# > print $SC_square -> POINTS;
# | 1 1 2 2 1
# | 1 2 1 1 2
user_function secondary_configuration<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 GKZ secondary polytope of a point configuration via its
# using topcom or mptopcom.
#
# @param PointConfiguration pc input point configuration
# @return Polytope
#
# @example The following [[PointConfiguration]] is called the "mother of all
# examples (moae)". It has two non-regular triangulations, which can be seen
# when comparing the number of points in the secondary configuration with the
# number of vertices of the secondary polytope.
# > $moae = new PointConfiguration(POINTS => [[1,4,0,0],[1,0,4,0],[1,0,0,4],[1,2,1,1],[1,1,2,1],[1,1,1,2]]);
# > $moae = project_full($moae);
# > $SC_moae = secondary_configuration($moae);
# > $SP_moae = secondary_polytope($moae);
# > print $SC_moae -> N_POINTS;
# | 18
# > print $SP_moae -> N_VERTICES;
# | 16
user_function secondary_polytope<Scalar>(PointConfiguration<Scalar>) {
	my $self=shift;
    my $v = new Matrix<Scalar>( map { gkz_vector<Scalar>($self->POINTS, $_) } @{topcom_regular_triangulations($self)} );
    my $vh = ones_vector<Scalar>() | $v;
    return new Polytope<Scalar>(VERTICES => $vh, LINEALITY_SPACE => []);
}

# @category Triangulations, subdivisions and volume
# Computes the GKZ secondary polytope of a point configuration via its
# using topcom or mptopcom.
#
# @param Polytope pc input polytope
# @return Polytope
#
# @example The square only has two triangulations using its vertices.
# > $square = cube(2,1,0);
# > $SP_square = secondary_polytope($square);
# > print $SP_square -> VERTICES;
# | 1 1 2 2 1
# | 1 2 1 1 2
user_function secondary_polytope<Scalar>(Polytope<Scalar>) {
    my $self=shift;
    my $v = new Matrix<Scalar>( map { gkz_vector<Scalar>($self->VERTICES, $_) } @{topcom_regular_triangulations($self)} );
    my $vh = ones_vector<Scalar>() | $v;
    return new Polytope<Scalar>(VERTICES => $vh, LINEALITY_SPACE => []);
}

# @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<Rational> {

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

}


object Polytope {

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

rule topcom.cocircuits : COCIRCUITS : CHIROTOPE {
   $this->COCIRCUITS=call_topcom_cocircuits($this, true);
}
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<Rational> {

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

}

object VectorConfiguration {

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

rule topcom.cocircuits : COCIRCUITS : CHIROTOPE {
   $this->COCIRCUITS=call_topcom_cocircuits($this, true);
}
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: