#!/usr/local/bin/perl

use application "polytope";
use Benchmark qw(:all);

use IO::Uncompress::UnXz qw( $UnXzError );
use strict ;
use warnings ;

# Call with 
# polymake --script massive_cache.pl
#
# This script will walk over all files named
#   "3d3.result.*.xz" (e.g. "3d3.result.000.xz")
# and produce a file 
#   "3d3.result.*.xz.gkset.set"
# containing the massive gkz vectors of the triangulations of the input. The
# output file is in polymake format and can be loaded in polymake. It contains
# a Set<Vector<Rational>>.
#
# Furthermore this script will produce a cache of gkz vectors for simplices and
# store it for further runs of this script.

my $s = simplex(3,3);

# Cache of massive gkz vectors for simplices.
# If there is a file with a stored cache, load it.
my $massive_cache;
if(! -e "saved_massive_cache.map"){
   $massive_cache = new Map<Set<Int>, Vector<Rational>>();
} else {
   $massive_cache = load_data("saved_massive_cache.map");
}

my $pc = new PointConfiguration(POINTS=>$s->LATTICE_POINTS);
my $mult = build_mult($pc);
my $num_triangs = 0;

foreach my $file (glob("3d3.result.*.xz")) {
   print "Entering ", $file, "\n";
   my $t0= Benchmark->new;
   get_gkset($file);
   my $t1=Benchmark->new;
   my $td1=timediff($t1,$t0);
   print "$file: Elapsed time: ", timestr($td1), "\n";
}

sub get_gkset{
   my($file) = @_;

   my ($name) = $file =~ m/([^\/]*)$/;
   if( -e "$name.gkset.set"){
      return;
   }
   my $line ;
   my $gkset = new Set<Vector<Rational>>();
   my $gz = new IO::Uncompress::UnXz $file or die "Cannot open $file: $UnXzError\n" ;
   my $i = 0;
   while(($line = $gz->getline())){
      my $facets;
      if($line =~ m/^\[\[.*\]\]$/){
         $facets = new Array<Set<Int>>(eval $line);
      } else {
         my($triang) = $line =~ m/(\{\{.*\}\})/;
         $triang =~ s/\{/\[/g;
         $triang =~ s/\}/\]/g;
      $facets = new Array<Set<Int>>(eval $triang);
      }
      my $mgkz = eta($pc, $mult, $facets);
      $gkset += $mgkz;
      if($num_triangs % 500000 == 0){
         print $num_triangs," ",$gkset->size(),"\n";
         $| = 1;
      }
      $i++;
      $num_triangs++;
   };
   print "Writing result for file $file containig $i triangulations \n";
   save_data($gkset, "$name.gkset.set");
   $gz->close() ;
}

save_data($massive_cache, "saved_massive_cache.map");

# The multiplicity of a massive face of a triangulation only depends on the
# minimal face of the simplex that it is contained in. This method walks over
# the faces of the simplex and stores the multiplicity for any face, i.e. it
# stores the number of paths from a given face to the top in the Hasse diagram.
sub build_mult {
   my($pc) = @_;
   my $mult = new Map<Int, Map<Set<Int>, Int>>;
   my $d = $pc->DIM;
   my $facets = $pc->CONVEX_HULL->POINTS_IN_FACETS;
   my $all = sequence(0, $pc->POINTS->rows());
   $mult->{$d}->{$all} = 1;
   for(my $i=$d-1; $i>=0; $i--){
      my @oldFaces = keys %{$mult->{$i+1}};
      foreach my $facet (@$facets){
         foreach my $oldFace (@oldFaces){
            my $intersection = $oldFace*$facet;
            if($intersection != $oldFace){
               $mult->{$i}->{$intersection} += $mult->{$i+1}->{$oldFace};
            }
         }
      }
   }
   return $mult;
}

# Build massive gkz vector for simplex.
# Returns a matrix whose rows are the dimension dependent massive gkz vectors
# of the given simplex.
sub eta_simplex {
   my($simplex, $pc, $mult) = @_;
   my $d = $pc->DIM;
   my $facets = $pc->CONVEX_HULL->POINTS_IN_FACETS;
   my $result = new Matrix<Rational>($d+1, $pc->POINTS->rows());
   my $vol = (new Polytope(POINTS=>$pc->POINTS->minor($simplex, All)))->LATTICE_VOLUME;
   foreach my $vertex (@$simplex) {
      $result->elem($d, $vertex) = $vol;
   }
   my @massive = ($simplex);
   for(my $i=$d-1; $i>=0; $i--){
      my $vols = new Map<Set<Int>, Rational>;
      my $containers = new Map<Set<Int>, Set<Int>>;
      my @newFaces = map {
         my $facet = $_;
         map($_*$facet, @massive);
      } @$facets;
      @newFaces = grep {
         my $face = $_;
         my @a = grep($face * $_ == $face, keys %{$mult->{$i}});
         if(@a > 0){
            $containers->{$face} = $a[0];
            $vols->{$face} = (new Polytope(POINTS=>$pc->POINTS->minor($face, All)))->LATTICE_VOLUME;
         }
         @a > 0;
      } @newFaces;
      @newFaces = grep {
         my $dim = (new Polytope(POINTS=>$pc->POINTS->minor($_, All)))->DIM;
         $dim == $i;
      } @newFaces;
      foreach my $face (@newFaces) {
         my $m = new Rational($mult->{$i}->{$containers->{$face}});
         foreach my $vertex (@$face) {
            $result->elem($i, $vertex) += (new Rational($vols->{$face}))/$m;
         }
      }
      @massive = @newFaces;
   }
   return $result;
}

# Ask for massive gkz vector of simplex. If it is not in the cache, compute it
# and store it in the cache.
sub eta_simplex_cached {
   my($pc, $mult,  $simplex) = @_;
   if(! defined $massive_cache->{$simplex}){
      my $A = eta_simplex($simplex, $pc, $mult);
      my $result = new Vector<Rational>($pc->POINTS->rows());
      my $d = $pc->DIM;
      for(my $i=0; $i<=$d; $i++){
         $result += $i%2 == 1 ? $A->row($i) : -$A->row($i);
      }
      $massive_cache->{$simplex} = $result;
   }
   return $massive_cache->{$simplex}
}

# Produce a massive gkz vector of a triangulation using the summation formula.
sub eta {
   my($pc, $mult, $triangulation) = @_;
   my $result = new Vector<Rational>($pc->POINTS->rows());
   foreach my $simplex (@$triangulation) {
      $result += eta_simplex_cached($pc, $mult, $simplex);
   }
   return $result;
}