1 #include "mlir/Analysis/Presburger/Barvinok.h"
2 #include "./Utils.h"
3 #include "Parser.h"
4 #include <gmock/gmock.h>
5 #include <gtest/gtest.h>
6
7 using namespace mlir;
8 using namespace presburger;
9 using namespace mlir::presburger::detail;
10
11 /// The following are 3 randomly generated vectors with 4
12 /// entries each and define a cone's H-representation
13 /// using these numbers. We check that the dual contains
14 /// the same numbers.
15 /// We do the same in the reverse case.
TEST(BarvinokTest,getDual)16 TEST(BarvinokTest, getDual) {
17 ConeH cone1 = defineHRep(4);
18 cone1.addInequality({1, 2, 3, 4, 0});
19 cone1.addInequality({3, 4, 2, 5, 0});
20 cone1.addInequality({6, 2, 6, 1, 0});
21
22 ConeV dual1 = getDual(cone1);
23
24 EXPECT_EQ_INT_MATRIX(
25 dual1, makeIntMatrix(3, 4, {{1, 2, 3, 4}, {3, 4, 2, 5}, {6, 2, 6, 1}}));
26
27 ConeV cone2 = makeIntMatrix(3, 4, {{3, 6, 1, 5}, {3, 1, 7, 2}, {9, 3, 2, 7}});
28
29 ConeH dual2 = getDual(cone2);
30
31 ConeH expected = defineHRep(4);
32 expected.addInequality({3, 6, 1, 5, 0});
33 expected.addInequality({3, 1, 7, 2, 0});
34 expected.addInequality({9, 3, 2, 7, 0});
35
36 EXPECT_TRUE(dual2.isEqual(expected));
37 }
38
39 /// We randomly generate a nxn matrix to use as a cone
40 /// with n inequalities in n variables and check for
41 /// the determinant being equal to the index.
TEST(BarvinokTest,getIndex)42 TEST(BarvinokTest, getIndex) {
43 ConeV cone = makeIntMatrix(3, 3, {{4, 2, 1}, {5, 2, 7}, {4, 1, 6}});
44 EXPECT_EQ(getIndex(cone), cone.determinant());
45
46 cone = makeIntMatrix(
47 4, 4, {{4, 2, 5, 1}, {4, 1, 3, 6}, {8, 2, 5, 6}, {5, 2, 5, 7}});
48 EXPECT_EQ(getIndex(cone), cone.determinant());
49 }
50
51 // The following cones and vertices are randomly generated
52 // (s.t. the cones are unimodular) and the generating functions
53 // are computed. We check that the results contain the correct
54 // matrices.
TEST(BarvinokTest,unimodularConeGeneratingFunction)55 TEST(BarvinokTest, unimodularConeGeneratingFunction) {
56 ConeH cone = defineHRep(2);
57 cone.addInequality({0, -1, 0});
58 cone.addInequality({-1, -2, 0});
59
60 ParamPoint vertex =
61 makeFracMatrix(2, 3, {{2, 2, 0}, {-1, -Fraction(1, 2), 1}});
62
63 GeneratingFunction gf =
64 computeUnimodularConeGeneratingFunction(vertex, 1, cone);
65
66 EXPECT_EQ_REPR_GENERATINGFUNCTION(
67 gf, GeneratingFunction(
68 2, {1},
69 {makeFracMatrix(3, 2, {{-1, 0}, {-Fraction(1, 2), 1}, {1, 2}})},
70 {{{2, -1}, {-1, 0}}}));
71
72 cone = defineHRep(3);
73 cone.addInequality({7, 1, 6, 0});
74 cone.addInequality({9, 1, 7, 0});
75 cone.addInequality({8, -1, 1, 0});
76
77 vertex = makeFracMatrix(3, 2, {{5, 2}, {6, 2}, {7, 1}});
78
79 gf = computeUnimodularConeGeneratingFunction(vertex, 1, cone);
80
81 EXPECT_EQ_REPR_GENERATINGFUNCTION(
82 gf,
83 GeneratingFunction(
84 1, {1}, {makeFracMatrix(2, 3, {{-83, -100, -41}, {-22, -27, -15}})},
85 {{{8, 47, -17}, {-7, -41, 15}, {1, 5, -2}}}));
86 }
87
88 // The following vectors are randomly generated.
89 // We then check that the output of the function has non-zero
90 // dot product with all non-null vectors.
TEST(BarvinokTest,getNonOrthogonalVector)91 TEST(BarvinokTest, getNonOrthogonalVector) {
92 std::vector<Point> vectors = {Point({1, 2, 3, 4}), Point({-1, 0, 1, 1}),
93 Point({2, 7, 0, 0}), Point({0, 0, 0, 0})};
94 Point nonOrth = getNonOrthogonalVector(vectors);
95
96 for (unsigned i = 0; i < 3; i++)
97 EXPECT_NE(dotProduct(nonOrth, vectors[i]), 0);
98
99 vectors = {Point({0, 1, 3}), Point({-2, -1, 1}), Point({6, 3, 0}),
100 Point({0, 0, -3}), Point({5, 0, -1})};
101 nonOrth = getNonOrthogonalVector(vectors);
102
103 for (const Point &vector : vectors)
104 EXPECT_NE(dotProduct(nonOrth, vector), 0);
105 }
106
107 // The following polynomials are randomly generated and the
108 // coefficients are computed by hand.
109 // Although the function allows the coefficients of the numerator
110 // to be arbitrary quasipolynomials, we stick to constants for simplicity,
111 // as the relevant arithmetic operations on quasipolynomials
112 // are tested separately.
TEST(BarvinokTest,getCoefficientInRationalFunction)113 TEST(BarvinokTest, getCoefficientInRationalFunction) {
114 std::vector<QuasiPolynomial> numerator = {
115 QuasiPolynomial(0, 2), QuasiPolynomial(0, 3), QuasiPolynomial(0, 5)};
116 std::vector<Fraction> denominator = {Fraction(1), Fraction(0), Fraction(4),
117 Fraction(3)};
118 QuasiPolynomial coeff =
119 getCoefficientInRationalFunction(1, numerator, denominator);
120 EXPECT_EQ(coeff.getConstantTerm(), 3);
121
122 numerator = {QuasiPolynomial(0, -1), QuasiPolynomial(0, 4),
123 QuasiPolynomial(0, -2), QuasiPolynomial(0, 5),
124 QuasiPolynomial(0, 6)};
125 denominator = {Fraction(8), Fraction(4), Fraction(0), Fraction(-2)};
126 coeff = getCoefficientInRationalFunction(3, numerator, denominator);
127 EXPECT_EQ(coeff.getConstantTerm(), Fraction(55, 64));
128 }
129
TEST(BarvinokTest,computeNumTermsCone)130 TEST(BarvinokTest, computeNumTermsCone) {
131 // The following test is taken from
132 // Verdoolaege, Sven, et al. "Counting integer points in parametric
133 // polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
134 // 37-66.
135 // It represents a right-angled triangle with right angle at the origin,
136 // with height and base lengths (p/2).
137 GeneratingFunction gf(
138 1, {1, 1, 1},
139 {makeFracMatrix(2, 2, {{0, Fraction(1, 2)}, {0, 0}}),
140 makeFracMatrix(2, 2, {{0, Fraction(1, 2)}, {0, 0}}),
141 makeFracMatrix(2, 2, {{0, 0}, {0, 0}})},
142 {{{-1, 1}, {-1, 0}}, {{1, -1}, {0, -1}}, {{1, 0}, {0, 1}}});
143
144 QuasiPolynomial numPoints = computeNumTerms(gf).collectTerms();
145
146 // First, we make sure that all the affine functions are of the form ⌊p/2⌋.
147 for (const std::vector<SmallVector<Fraction>> &term : numPoints.getAffine()) {
148 for (const SmallVector<Fraction> &aff : term) {
149 EXPECT_EQ(aff.size(), 2u);
150 EXPECT_EQ(aff[0], Fraction(1, 2));
151 EXPECT_EQ(aff[1], Fraction(0, 1));
152 }
153 }
154
155 // Now, we can gather the like terms because we know there's only
156 // either ⌊p/2⌋^2, ⌊p/2⌋, or constants.
157 // The total coefficient of ⌊p/2⌋^2 is the sum of coefficients of all
158 // terms with 2 affine functions, and
159 // the coefficient of total ⌊p/2⌋ is the sum of coefficients of all
160 // terms with 1 affine function,
161 Fraction pSquaredCoeff = 0, pCoeff = 0, constantTerm = 0;
162 SmallVector<Fraction> coefficients = numPoints.getCoefficients();
163 for (unsigned i = 0; i < numPoints.getCoefficients().size(); i++)
164 if (numPoints.getAffine()[i].size() == 2)
165 pSquaredCoeff = pSquaredCoeff + coefficients[i];
166 else if (numPoints.getAffine()[i].size() == 1)
167 pCoeff = pCoeff + coefficients[i];
168
169 // We expect the answer to be (1/2)⌊p/2⌋^2 + (3/2)⌊p/2⌋ + 1.
170 EXPECT_EQ(pSquaredCoeff, Fraction(1, 2));
171 EXPECT_EQ(pCoeff, Fraction(3, 2));
172 EXPECT_EQ(numPoints.getConstantTerm(), Fraction(1, 1));
173
174 // The following generating function corresponds to a cuboid
175 // with length M (x-axis), width N (y-axis), and height P (z-axis).
176 // There are eight terms.
177 gf = GeneratingFunction(
178 3, {1, 1, 1, 1, 1, 1, 1, 1},
179 {makeFracMatrix(4, 3, {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}}),
180 makeFracMatrix(4, 3, {{1, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}}),
181 makeFracMatrix(4, 3, {{0, 0, 0}, {0, 1, 0}, {0, 0, 0}, {0, 0, 0}}),
182 makeFracMatrix(4, 3, {{0, 0, 0}, {0, 0, 0}, {0, 0, 1}, {0, 0, 0}}),
183 makeFracMatrix(4, 3, {{1, 0, 0}, {0, 1, 0}, {0, 0, 0}, {0, 0, 0}}),
184 makeFracMatrix(4, 3, {{1, 0, 0}, {0, 0, 0}, {0, 0, 1}, {0, 0, 0}}),
185 makeFracMatrix(4, 3, {{0, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}}),
186 makeFracMatrix(4, 3, {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}})},
187 {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
188 {{-1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
189 {{1, 0, 0}, {0, -1, 0}, {0, 0, 1}},
190 {{1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
191 {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}},
192 {{-1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
193 {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
194 {{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}}});
195
196 numPoints = computeNumTerms(gf);
197 numPoints = numPoints.collectTerms().simplify();
198
199 // First, we make sure all the affine functions are either
200 // M, N, P, or 1.
201 for (const std::vector<SmallVector<Fraction>> &term : numPoints.getAffine()) {
202 for (const SmallVector<Fraction> &aff : term) {
203 // First, ensure that the coefficients are all nonnegative integers.
204 for (const Fraction &c : aff) {
205 EXPECT_TRUE(c >= 0);
206 EXPECT_EQ(c, c.getAsInteger());
207 }
208 // Now, if the coefficients add up to 1, we can be sure the term is
209 // either M, N, P, or 1.
210 EXPECT_EQ(aff[0] + aff[1] + aff[2] + aff[3], 1);
211 }
212 }
213
214 // We store the coefficients of M, N and P in this array.
215 Fraction count[2][2][2];
216 coefficients = numPoints.getCoefficients();
217 for (unsigned i = 0, e = coefficients.size(); i < e; i++) {
218 unsigned mIndex = 0, nIndex = 0, pIndex = 0;
219 for (const SmallVector<Fraction> &aff : numPoints.getAffine()[i]) {
220 if (aff[0] == 1)
221 mIndex = 1;
222 if (aff[1] == 1)
223 nIndex = 1;
224 if (aff[2] == 1)
225 pIndex = 1;
226 EXPECT_EQ(aff[3], 0);
227 }
228 count[mIndex][nIndex][pIndex] += coefficients[i];
229 }
230
231 // We expect the answer to be
232 // (⌊M⌋ + 1)(⌊N⌋ + 1)(⌊P⌋ + 1) =
233 // ⌊M⌋⌊N⌋⌊P⌋ + ⌊M⌋⌊N⌋ + ⌊N⌋⌊P⌋ + ⌊M⌋⌊P⌋ + ⌊M⌋ + ⌊N⌋ + ⌊P⌋ + 1.
234 for (unsigned i = 0; i < 2; i++)
235 for (unsigned j = 0; j < 2; j++)
236 for (unsigned k = 0; k < 2; k++)
237 EXPECT_EQ(count[i][j][k], 1);
238 }
239
240 /// We define some simple polyhedra with unimodular tangent cones and verify
241 /// that the returned generating functions correspond to those calculated by
242 /// hand.
TEST(BarvinokTest,computeNumTermsPolytope)243 TEST(BarvinokTest, computeNumTermsPolytope) {
244 // A cube of side 1.
245 PolyhedronH poly =
246 parseRelationFromSet("(x, y, z) : (x >= 0, y >= 0, z >= 0, -x + 1 >= 0, "
247 "-y + 1 >= 0, -z + 1 >= 0)",
248 0);
249
250 std::vector<std::pair<PresburgerSet, GeneratingFunction>> count =
251 computePolytopeGeneratingFunction(poly);
252 // There is only one chamber, as it is non-parametric.
253 EXPECT_EQ(count.size(), 9u);
254
255 GeneratingFunction gf = count[0].second;
256 EXPECT_EQ_REPR_GENERATINGFUNCTION(
257 gf,
258 GeneratingFunction(
259 0, {1, 1, 1, 1, 1, 1, 1, 1},
260 {makeFracMatrix(1, 3, {{1, 1, 1}}), makeFracMatrix(1, 3, {{0, 1, 1}}),
261 makeFracMatrix(1, 3, {{0, 1, 1}}), makeFracMatrix(1, 3, {{0, 0, 1}}),
262 makeFracMatrix(1, 3, {{0, 1, 1}}), makeFracMatrix(1, 3, {{0, 0, 1}}),
263 makeFracMatrix(1, 3, {{0, 0, 1}}),
264 makeFracMatrix(1, 3, {{0, 0, 0}})},
265 {{{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
266 {{0, 0, 1}, {-1, 0, 0}, {0, -1, 0}},
267 {{0, 1, 0}, {-1, 0, 0}, {0, 0, -1}},
268 {{0, 1, 0}, {0, 0, 1}, {-1, 0, 0}},
269 {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
270 {{1, 0, 0}, {0, 0, 1}, {0, -1, 0}},
271 {{1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
272 {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}}));
273
274 // A right-angled triangle with side p.
275 poly =
276 parseRelationFromSet("(x, y)[N] : (x >= 0, y >= 0, -x - y + N >= 0)", 0);
277
278 count = computePolytopeGeneratingFunction(poly);
279 // There is only one chamber: p ≥ 0
280 EXPECT_EQ(count.size(), 4u);
281
282 gf = count[0].second;
283 EXPECT_EQ_REPR_GENERATINGFUNCTION(
284 gf, GeneratingFunction(
285 1, {1, 1, 1},
286 {makeFracMatrix(2, 2, {{0, 1}, {0, 0}}),
287 makeFracMatrix(2, 2, {{0, 1}, {0, 0}}),
288 makeFracMatrix(2, 2, {{0, 0}, {0, 0}})},
289 {{{-1, 1}, {-1, 0}}, {{1, -1}, {0, -1}}, {{1, 0}, {0, 1}}}));
290
291 // Cartesian product of a cube with side M and a right triangle with side N.
292 poly = parseRelationFromSet(
293 "(x, y, z, w, a)[M, N] : (x >= 0, y >= 0, z >= 0, -x + M >= 0, -y + M >= "
294 "0, -z + M >= 0, w >= 0, a >= 0, -w - a + N >= 0)",
295 0);
296
297 count = computePolytopeGeneratingFunction(poly);
298
299 EXPECT_EQ(count.size(), 25u);
300
301 gf = count[0].second;
302 EXPECT_EQ(gf.getNumerators().size(), 24u);
303 }
304