mlir/unittests/Analysis/Presburger/BarvinokTest.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
   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.
  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.
  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.
  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.
  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.
 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
 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.
 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 }