304 lines · cpp
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 7using namespace mlir;8using namespace presburger;9using namespace mlir::presburger::detail;10 11/// The following are 3 randomly generated vectors with 412/// entries each and define a cone's H-representation13/// using these numbers. We check that the dual contains14/// the same numbers.15/// We do the same in the reverse case.16TEST(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 cone40/// with n inequalities in n variables and check for41/// the determinant being equal to the index.42TEST(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 generated52// (s.t. the cones are unimodular) and the generating functions53// are computed. We check that the results contain the correct54// matrices.55TEST(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-zero90// dot product with all non-null vectors.91TEST(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 the108// coefficients are computed by hand.109// Although the function allows the coefficients of the numerator110// to be arbitrary quasipolynomials, we stick to constants for simplicity,111// as the relevant arithmetic operations on quasipolynomials112// are tested separately.113TEST(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 130TEST(BarvinokTest, computeNumTermsCone) {131 // The following test is taken from132 // Verdoolaege, Sven, et al. "Counting integer points in parametric133 // 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 only156 // either ⌊p/2⌋^2, ⌊p/2⌋, or constants.157 // The total coefficient of ⌊p/2⌋^2 is the sum of coefficients of all158 // terms with 2 affine functions, and159 // the coefficient of total ⌊p/2⌋ is the sum of coefficients of all160 // 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 cuboid175 // 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 either200 // 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 is209 // 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 be232 // (⌊M⌋ + 1)(⌊N⌋ + 1)(⌊P⌋ + 1) =233 // ⌊M⌋⌊N⌋⌊P⌋ + ⌊M⌋⌊N⌋ + ⌊N⌋⌊P⌋ + ⌊M⌋⌊P⌋ + ⌊M⌋ + ⌊N⌋ + ⌊P⌋ + 1.234 for (auto &i : count)235 for (unsigned j = 0; j < 2; j++)236 for (unsigned k = 0; k < 2; k++)237 EXPECT_EQ(i[j][k], 1);238}239 240/// We define some simple polyhedra with unimodular tangent cones and verify241/// that the returned generating functions correspond to those calculated by242/// hand.243TEST(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 ≥ 0280 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