#  Copyright (c) 1997-2006
#  Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
#  http://www.math.tu-berlin.de/polymake,  mailto:polymake@math.tu-berlin.de
#
#  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.
#-----------------------------------------------------------------------------
#  $Project: polymake $$Id: porta.rules 7044 2006-02-21 19:34:28Z gawrilow $


# path to the xporta executable
custom $xporta;

CONFIGURE {
   $xporta=find_via_path("xporta")
   or die "PORTA package not found\n";
}

object RationalPolytope;

# Category: Convex hull computation
# Run the @see external.porta program, implementing the Fourier-Motzkin elimination method.
# The essential drawback of this tool is that it employs a limited-precision arithmetic,
# and therefore can fail on numerically difficult problems.
label porta

sub hextodec {
   my @dec=(0);
   foreach my $h (split(//,shift)) {
      $h="1".$h  if $h=~tr/abcdef/012345/;
      for (my $i=$#dec; $i>=0; $i--) {
	 my $d=$dec[$i]*16+$h;  $dec[$i]=$d%10;  $h=$d/10;
      }
      unshift (@dec,$h%10), $h/=10  if $h;
      unshift (@dec,$h)             if $h;
   }
   join "",@dec
}

sub conv_porta_number {
   my $n=shift;
   $n=~s/\(hex\)([\da-f]+)/&hextodec($1)/eg;
   $n
}

porta.convex_hull.primal: FACETS, AFFINE_HULL : AMBIENT_DIM, DIM, VERTICES | POINTS
WEIGHT 4.20
   my (@conv, @cone);
   foreach (convert_to_rational($this->VERTICES | POINTS)) {
      my ($first, $rest)=split /\s+/,$_,2;
      if ($first eq "0") {
	 push @cone, $rest;
      } else {
	 push @conv, $rest;
      }
   }

   my $tempname=new Poly::Tempfile;
   open P, ">$tempname.poi"
   or die "can't create temporary file $tempname.poi: $!";
   print P "DIM=", $this->AMBIENT_DIM, "\n";
   print P "CONV_SECTION\n", @conv;
   print P "CONE_SECTION\n", @cone  if @cone;
   print P "END\n";
   close P;

   my $command="$xporta -l -T $tempname.poi";
   $command.=" >/dev/null 2>&1" if !$Switches::d;
   system($command) and die "porta exited with error code\n";
   open P, "$tempname.poi.ieq"
   or die "can't read porta results from $tempname.poi.ieq: $!";
   while (<P>) {
      last if $_ eq "INEQUALITIES_SECTION\n";
   }

   my (@AFFINE_HULL, @FACETS);
   while (<P>) {
      next unless /\)(.*)([<=])=/;
      $_=$1;
      my $to= $2 eq "=" ? \@AFFINE_HULL : \@FACETS;
      my $right=$';
      $right=~s/\s//g;
      $right=conv_porta_number($right);
      my $n=1;
      my $x="";
      while (/([+-])(.*?)x(\d+)/) {
	 $_=$';
	 $x.=" 0", $n++  while $n<$3;
	 $x.= $1 eq "-" ? " " : " -";
	 my $coef=$2;
	 $coef=~s/\s//g;
	 $x.= $coef ne "" ? conv_porta_number($coef) : "1";
	 $n++;
      }
      $x.=" 0"x($this->AMBIENT_DIM-$n+1);
      push @$to, "$right$x\n";
   }
   close P;

   my $r=[ ];	# a scratch section
   client("dimension", $this, "-rays", $r);
   if ($r->[0] == $this->DIM) {
      unshift @FACETS, "1".(" 0" x $this->AMBIENT_DIM)."\n";
   }
   $this->FACETS=\@FACETS;
   $this->AFFINE_HULL=\@AFFINE_HULL;
   unlink "porta.log";


porta.convex_hull.dual: VERTICES : AMBIENT_DIM, FACETS | INEQUALITIES, VALID_POINT
WEIGHT 4.10
   my $tempname=new Poly::Tempfile;
   open P, ">$tempname.ieq"
   or die "can't create temporary file $tempname.ieq: $!";
   my $valid_point=$this->VALID_POINT;
   $valid_point =~ s/^(\d+) //;
   print P "DIM=", $this->AMBIENT_DIM,
           "\nVALID\n$valid_point\nINEQUALITIES_SECTION\n";

   foreach my $line (convert_to_rational($this->FACETS | INEQUALITIES)) {
      my ($right, @x)=map /^\d/ ? "+$_" : $_,  split(/\s+/,$line);
      $right=~tr/+-/-+/;
      my $n=0;
      my $left=join("",map(($n++, $_ eq "+0") ? "" : $_."x".$n, @x));
      if ($left ne "") {
	 print P "$left>=$right\n";
      }
      elsif (substr($right,0,1) eq "+") {
	 die "polytope is empty or definition inconsistent";
      }
   }

   if (defined (my $AFFINE_HULL=$this->lookup("AFFINE_HULL | EQUATIONS"))) {
      foreach my $line (convert_to_rational($AFFINE_HULL)) {
	 my ($right, @x)=map /^\d/ ? "+$_" : $_,  split(/\s+/,$line);
	 $right=~tr/+-/-+/;
	 my $n=0;
	 my $left=join("",map(($n++, $_ eq "+0") ? "" : $_."x".$n, @x));
	 if ($left ne "") {
	    print P "$left=$right\n";
	 }
	 elsif ($right ne "-0") {
	    die "polytope is empty or definition inconsistent";
	 }
      }
   }
   print P "END\n";
   close P;

   my $command="$xporta -l -T $tempname.ieq";
   $command.=" >/dev/null 2>&1" if !$Switches::d;
   system($command) and die "xporta exited with error code\n";
   open P, "$tempname.ieq.poi"
   or die "can't read porta results from $tempname.ieq.poi: $!";

   my $first=-1;
   my @vertices;
   while (<P>) {
      if (/CON([VE])_SECTION/) {
	 $first= $1 eq "V" ? 1 : 0;
	 next;
      }
      next if $first<0;
      if (/\)\s*/) {
	 $_=$';
	 s"\s*/\s*"/"g;
         push @vertices, "$first $_";
      }
   }
   close P;

   $this->VERTICES=[ map { join(" ", map { conv_porta_number($_) } split) } @vertices ];
   unlink "porta.log";


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