1#  Copyright (c) 1997-2021
2#  Ewgenij Gawrilow, Michael Joswig, and the polymake team
3#  Technische Universität Berlin, Germany
4#  https://polymake.org
5#
6#  This program is free software; you can redistribute it and/or modify it
7#  under the terms of the GNU General Public License as published by the
8#  Free Software Foundation; either version 2, or (at your option) any
9#  later version: http://www.gnu.org/licenses/gpl.txt.
10#
11#  This program is distributed in the hope that it will be useful,
12#  but WITHOUT ANY WARRANTY; without even the implied warranty of
13#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14#  GNU General Public License for more details.
15#-------------------------------------------------------------------------------
16
17CREDIT latte
18  LattE (Lattice point Enumeration) is a computer software dedicated to the
19  problems of counting lattice points and integration inside convex polytopes.
20  Copyright by Matthias Koeppe, Jesus A. De Loera and others.
21  http://www.math.ucdavis.edu/~latte/
22
23# path to the LattE binaries
24custom $latte_count;
25
26CONFIGURE {
27   find_program($latte_count, "count", {
28                   prompt => "the program `count' from the LattE package",
29                   check => sub { `$_[0] 2>&1` !~ /LattE/ && "this is not LattE's count" }
30                }) or return;
31   `$latte_count 2>&1` =~ /LattE/ or die <<".";
32'$latte_count' is not the count binary from LattE.
33.
34}
35
36# additional parameters for calling LattE's "count". See the manual of LattE for details, recommended options are
37# "--irrational-all-primal": use irrational decomposition instead of polarization
38# "--maxdet=<n>": stop triangulating if determinant of cone < n (needs LiDIA)
39# "--exponential": use exponential substitution for evaluating the generating function
40custom $latte_count_param="";
41
42# additional parameters for calling latte's "count --ehrhart-polynomial"
43# see $latte_count_param for details
44custom $latte_ehrhart_param="";
45
46# @category Lattice points in polytopes
47# Use [[wiki:external_software#latte|LattE]] for counting and detecting lattice points inside convex polytopes.
48label latte
49
50sub write_latte_ineq {
51   my ($this, $tempname, $int)=@_;
52   # daten: ineq / eq
53   my $ineq = $this->give("FACETS | INEQUALITIES");
54   my $n_ineq = $ineq->rows;
55   my $fac = dense(eliminate_denominators_in_rows($ineq));
56   my $n_lines = $n_ineq;
57   my $ah;
58   if (defined (my $AH = $this->lookup("AFFINE_HULL | EQUATIONS"))) {
59      $n_lines = $n_lines + $AH->rows;
60      $ah = dense(eliminate_denominators_in_rows($AH));
61   }
62
63   open(my $P, ">$tempname")
64    or die "can't create temporary file $tempname: $!";
65   #header
66   print $P $n_lines, " ", $this->CONE_AMBIENT_DIM, "\n";
67
68   # facet lines
69   if ($int) {
70      $fac *= $this->CONE_DIM;
71      for (my $i=0; $i<$fac->rows; $i++) {
72        $fac->elem($i,0) -= 1;
73      }
74   }
75   print $P $fac;
76
77   # affine hull lines
78   if ($n_lines > $n_ineq) {
79      print $P $ah;
80      print $P "linearity ", $n_lines-$n_ineq, " ";
81      print $P ++$n_ineq," " while $n_ineq < $n_lines;
82      print $P "\n";
83   }
84   close $P;
85}
86
87sub write_latte_vertices {
88   my ($this, $tempname)=@_;
89   my $vert = dense($this->VERTICES);
90
91   open(my $P, ">$tempname")
92    or die "can't create temporary file $tempname: $!";
93
94   print $P $vert->rows, " ", ($this->CONE_AMBIENT_DIM), "\n";
95   print $P $vert;
96
97   close $P;
98}
99
100sub parse_latte_points {
101   my $P = shift;
102   local $_;
103   while (<$P>) {
104       if (my ($n)=/number of lattice points\D+(\d+)/) {
105           close $P;
106           return $n;
107       } elsif(/Empty polytope/ or /The polytope is empty!/) {
108           close $P;
109           return 0;
110       }
111   }
112   close $P;
113   die "can't parse output from LattE's 'count'\n";
114}
115
116object Polytope<Rational> {
117
118    rule latte.integer_points: N_LATTICE_POINTS : CONE_AMBIENT_DIM, CONE_DIM, FACETS | INEQUALITIES {
119        my $temp_dir=new Tempdir;
120	my $input_file="input.ine";
121        my $interior = 0;
122        write_latte_ineq($this, "$temp_dir/$input_file", $interior);
123
124        # latte 1.2 does not do a full redundancy check by default, so force it when using inequalities
125        # furthermore it might produce wrong output for 0-dim polytopes when given the far face inequality
126        my $check = "--redundancy-check=full-cddlib";
127        if (defined($this->lookup("FACETS")) && $this->CONE_DIM > 1) {
128            $check = "--redundancy-check=none";
129        }
130        my $tdir = new TempChangeDir($temp_dir);
131        if ($Verbose::external) {
132            dbg_print( "running latte's count: cd $temp_dir; $latte_count $check $latte_count_param $input_file" );
133        }
134        open my $P, "$latte_count $check $latte_count_param $input_file 2>&1 |"
135            or die "couldn't run LattE's 'count': $!";
136        $this->N_LATTICE_POINTS = parse_latte_points($P);
137    }
138    precondition : BOUNDED;
139    precondition : FEASIBLE;
140    weight 5.3;
141
142    rule latte.integer_points: N_LATTICE_POINTS : CONE_AMBIENT_DIM , VERTICES {
143	my $temp_dir=new Tempdir;
144        my $input_file="input.poi";
145        write_latte_vertices($this, "$temp_dir/$input_file");
146
147        my $tdir = new TempChangeDir($temp_dir);
148        if ($Verbose::external) {
149            dbg_print( "running latte's count: cd $temp_dir; $latte_count $latte_count_param vrep $input_file" );
150        }
151        open my $P, "$latte_count $latte_count_param vrep $input_file 2>&1 |"
152            or die "couldn't run LattE's 'count': $!";
153        $this->N_LATTICE_POINTS = parse_latte_points($P);
154
155        # FULL_DIM precodition because latte 1.2 fails for some non full-dim polytopes
156    }
157    precondition : BOUNDED;
158    precondition : FEASIBLE;
159    precondition : FULL_DIM;
160    weight 5.3;
161
162
163    # Read [[FACETS]], multiply with POLYTOPE_DIM+1 (or CONE_DIM, respectively!) and substract 1 from the first column.
164    # Call LattE's count on the polytope having this matrix as facet matrix.
165    rule latte.integer_points: N_INTERIOR_LATTICE_POINTS : CONE_AMBIENT_DIM, CONE_DIM, FACETS | INEQUALITIES {
166        my $temp_dir=new Tempdir;
167        my $input_file="input.ine";
168        my $interior = 1;
169        write_latte_ineq($this, "$temp_dir/$input_file", $interior);
170
171        # latte 1.2 does not do a full redundancy check by default, so force it when using inequalities
172        # furthermore it might produce wrong output for 0-dim polytopes when given the far face inequality
173        my $check = "--redundancy-check=full-cddlib";
174        if (defined($this->lookup("FACETS")) && $this->lookup("CONE_DIM") > 1) {
175            $check = "--redundancy-check=none";
176        }
177        my $tdir = new TempChangeDir($temp_dir);
178        if ($Verbose::external) {
179            dbg_print( "running latte's count: cd $temp_dir; $latte_count $check $latte_count_param $input_file" );
180        }
181        open my $P, "$latte_count $check $latte_count_param $input_file 2>&1 |"
182            or die "couldn't run LattE's 'count': $!";
183        $this->N_INTERIOR_LATTICE_POINTS = parse_latte_points($P);
184    }
185    weight 5.3;
186    precondition : BOUNDED;
187    precondition : FEASIBLE;
188
189}
190
191object_specialization Polytope::Lattice {
192
193    rule EHRHART_POLYNOMIAL : {
194	$this->EHRHART_POLYNOMIAL = new UniPolynomial<Rational>(1);
195    }
196    precondition : CONE_DIM { $this->CONE_DIM == 1 }
197
198    rule latte.ehrhartpoly: LATTICE, EHRHART_POLYNOMIAL : CONE_AMBIENT_DIM, FACETS | INEQUALITIES {
199        my $temp_dir=new Tempdir;
200        my $input_file="input.ine";
201        write_latte_ineq($this, "$temp_dir/$input_file");
202
203        my $tdir = new TempChangeDir($temp_dir);
204        if ($Verbose::external) {
205            dbg_print( "running latte's count: cd $temp_dir; $latte_count $latte_ehrhart_param --ehrhart-polynomial $input_file" );
206        }
207        open my $P, "$latte_count $latte_ehrhart_param --ehrhart-polynomial $input_file 2>&1 |"
208            or die "couldn't run LattE's 'count': $!";
209        local $_;
210        while (<$P>) {
211            if(/only implemented for integral polytopes/) {
212                $this->LATTICE = 0;
213                close $P;
214                return;
215            } elsif (/^Ehrhart polynomial: /) {
216                # parse: Ehrhart polynomial:  + 1 * t^0 + 6 * t^1 + 12 * t^2 + 8 * t^3
217                #					+ 1 * t^0 + 8/3 * t^1 + 2 * t^2 + 4/3 * t^3
218		s/ //g;
219                my @list = m/\+?(-?\d+\/?\d*)\*t\^(\d+)/gc;
220                my $coeffs = new Vector<Rational>($list[-1]+1);
221                while (@list) {
222                    $coeffs->[shift(@list)] = shift(@list);
223                }
224		my $exps = new Array<Int>(sequence(0,$coeffs->dim()));
225                $this->EHRHART_POLYNOMIAL = new UniPolynomial($coeffs,$exps);
226                $this->LATTICE=1;
227                close $P;
228                return;
229            } elsif (/^The number of lattice points is 1\./ || /Total number of lattice points: 1/) {
230                $this->LATTICE=1;
231                $this->EHRHART_POLYNOMIAL = new UniPolynomial<Rational>(1);
232                close $P;
233                return;
234            }
235        }
236        close $P;
237        die "could not parse output from latte";
238    }
239    weight 5.8;
240    precondition override : BOUNDED && FEASIBLE;
241    precondition : CONE_AMBIENT_DIM { $this->CONE_AMBIENT_DIM > 0 && $this->CONE_AMBIENT_DIM <= 100 } # upper limits for matrix dimension:
242    precondition : FACETS | INEQUALITIES { ($this->give("FACETS | INEQUALITIES"))->rows <= 50000 }            # default values stored explicitly for cdd used by latte
243    precondition : CONE_DIM { $this->CONE_DIM > 1 }                                                   # latte does not like polytopes that only have a far facet
244
245}
246
247
248# Local Variables:
249# mode: perl
250# cperl-indent-level:4
251# End:
252