[Mlir-commits] [mlir] [MLIR][Presburger][WIP] Implement vertex enumeration and chamber decomposition for polytope generating function computation. (PR #78987)

llvmlistbot at llvm.org llvmlistbot at llvm.org
Mon Jan 22 06:57:34 PST 2024


llvmbot wrote:


<!--LLVM PR SUMMARY COMMENT-->

@llvm/pr-subscribers-mlir-presburger

Author: None (Abhinav271828)

<details>
<summary>Changes</summary>

We implement a function to compute the generating function corresponding to a full-dimensional parametric polytope whose tangent cones are all unimodular.

---
Full diff: https://github.com/llvm/llvm-project/pull/78987.diff


5 Files Affected:

- (modified) mlir/include/mlir/Analysis/Presburger/Barvinok.h (+12) 
- (modified) mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h (+1-1) 
- (modified) mlir/include/mlir/Analysis/Presburger/IntegerRelation.h (+7) 
- (modified) mlir/lib/Analysis/Presburger/Barvinok.cpp (+312) 
- (modified) mlir/lib/Analysis/Presburger/IntegerRelation.cpp (+42) 


``````````diff
diff --git a/mlir/include/mlir/Analysis/Presburger/Barvinok.h b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
index b70ec33b693235f..a5bc49d79408c6c 100644
--- a/mlir/include/mlir/Analysis/Presburger/Barvinok.h
+++ b/mlir/include/mlir/Analysis/Presburger/Barvinok.h
@@ -27,6 +27,7 @@
 #include "mlir/Analysis/Presburger/GeneratingFunction.h"
 #include "mlir/Analysis/Presburger/IntegerRelation.h"
 #include "mlir/Analysis/Presburger/Matrix.h"
+#include "mlir/Analysis/Presburger/PresburgerRelation.h"
 #include "mlir/Analysis/Presburger/QuasiPolynomial.h"
 #include <optional>
 
@@ -84,6 +85,17 @@ ConeH getDual(ConeV cone);
 GeneratingFunction unimodularConeGeneratingFunction(ParamPoint vertex, int sign,
                                                     ConeH cone);
 
+/// Find the solution of a set of equations that express affine constraints
+/// between a set of variables and a set of parameters. The solution expresses
+/// each variable as an affine function of the parameters.
+/// If there is no solution, return std::null.
+std::optional<ParamPoint> findVertex(IntMatrix equations);
+
+/// Compute the generating function corresponding to a polytope.
+/// All tangent cones of the polytope must be unimodular.
+std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
+polytopeGeneratingFunction(PolyhedronH poly);
+
 /// Find a vector that is not orthogonal to any of the given vectors,
 /// i.e., has nonzero dot product with those of the given vectors
 /// that are not null.
diff --git a/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h b/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
index c38eab6efd0fc16..35f15319b4568eb 100644
--- a/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
+++ b/mlir/include/mlir/Analysis/Presburger/GeneratingFunction.h
@@ -72,7 +72,7 @@ class GeneratingFunction {
     return denominators;
   }
 
-  GeneratingFunction operator+(GeneratingFunction &gf) const {
+  GeneratingFunction operator+(GeneratingFunction gf) const {
     assert(numParam == gf.getNumParams() &&
            "two generating functions with different numbers of parameters "
            "cannot be added!");
diff --git a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
index c476a022a482725..141dc1f8fa97ed2 100644
--- a/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
+++ b/mlir/include/mlir/Analysis/Presburger/IntegerRelation.h
@@ -711,6 +711,13 @@ class IntegerRelation {
   /// return `this \ set`.
   PresburgerRelation subtract(const PresburgerRelation &set) const;
 
+  // Remove equalities which have only zero coefficients.
+  void removeTrivialEqualities();
+
+  // Verify whether the relation is full-dimensional, i.e.,
+  // has the same number of dimensions as the number of variables.
+  bool isFullDim();
+
   void print(raw_ostream &os) const;
   void dump() const;
 
diff --git a/mlir/lib/Analysis/Presburger/Barvinok.cpp b/mlir/lib/Analysis/Presburger/Barvinok.cpp
index e0fd0dd8caa4d31..3afc9ad8638ccd2 100644
--- a/mlir/lib/Analysis/Presburger/Barvinok.cpp
+++ b/mlir/lib/Analysis/Presburger/Barvinok.cpp
@@ -10,6 +10,7 @@
 #include "mlir/Analysis/Presburger/Utils.h"
 #include "llvm/ADT/Sequence.h"
 #include <algorithm>
+#include <bitset>
 
 using namespace mlir;
 using namespace presburger;
@@ -147,6 +148,317 @@ GeneratingFunction mlir::presburger::detail::unimodularConeGeneratingFunction(
                             std::vector({denominator}));
 }
 
+/// We use Gaussian elimination to find the solution to a set of d equations
+/// of the form
+/// a_1 x_1 + ... + a_d x_d + b_1 m_1 + ... + b_p m_p + c = 0
+/// where x_i are variables,
+/// m_i are parameters and
+/// a_i, b_i, c are rational coefficients.
+/// The solution expresses each x_i as an affine function of the m_i, and is
+/// therefore represented as a matrix of size d x (p+1).
+/// If there is no solution, we return null.
+std::optional<ParamPoint>
+mlir::presburger::detail::findVertex(IntMatrix equations) {
+  // equations is a d x (d + p + 1) matrix.
+  // Each row represents an equation.
+
+  unsigned numEqs = equations.getNumRows();
+  unsigned numCols = equations.getNumColumns();
+
+  // First, we check that the system has a solution, and return
+  // null if not.
+  IntMatrix coeffs(numEqs, numEqs);
+  for (unsigned i = 0; i < numEqs; i++)
+    for (unsigned j = 0; j < numEqs; j++)
+      coeffs(i, j) = equations(i, j);
+
+  if (coeffs.determinant() == 0)
+    return std::nullopt;
+
+  // We work with rational numbers.
+  FracMatrix equationsF(equations);
+
+  for (unsigned i = 0; i < numEqs; ++i) {
+    // First ensure that the diagonal element is nonzero, by swapping
+    // it with a nonzero row.
+    if (equationsF(i, i) == 0) {
+      for (unsigned j = i + 1; j < numEqs; ++j) {
+        if (equationsF(j, i) != 0) {
+          equationsF.swapRows(j, i);
+          break;
+        }
+      }
+    }
+
+    Fraction b = equationsF(i, i);
+
+    // Set all elements except the diagonal to zero.
+    for (unsigned j = 0; j < numEqs; ++j) {
+      if (equationsF(j, i) == 0 || j == i)
+        continue;
+      // Set element (j, i) to zero
+      // by subtracting the ith row,
+      // appropriately scaled.
+      Fraction a = equationsF(j, i);
+      equationsF.addToRow(j, equationsF.getRow(i), -a / b);
+    }
+  }
+
+  // Rescale diagonal elements to 1.
+  for (unsigned i = 0; i < numEqs; ++i) {
+    Fraction a = equationsF(i, i);
+    for (unsigned j = 0; j < numCols; ++j)
+      equationsF(i, j) = equationsF(i, j) / a;
+  }
+
+  // We copy the last p+1 columns of the matrix as the values of x_i.
+  // We shift the parameter terms to the RHS, and so flip their sign.
+  ParamPoint vertex(numEqs, numCols - numEqs);
+  for (unsigned i = 0; i < numEqs; ++i)
+    for (unsigned j = 0; j < numCols - numEqs; ++j)
+      vertex(i, j) = -equationsF(i, numEqs + j);
+
+  return vertex;
+}
+
+/// For a polytope expressed as a set of inequalities, compute the generating
+/// function corresponding to the number of lattice points present. This
+/// algorithm has three main steps:
+/// 1. Enumerate the vertices, by iterating over subsets of inequalities and
+///    checking for solubility.
+/// 2. For each vertex, identify the tangent cone and compute the generating
+///    function corresponding to it. The sum of these GFs is the GF of the
+///    polytope.
+/// 3. [Clauss-Loechner decomposition] Identify the regions in parameter space
+///    (chambers) where each vertex is active, and accordingly compute the
+///    GF of the polytope in each chamber.
+///
+/// Verdoolaege, Sven, et al. "Counting integer points in parametric
+/// polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
+/// 37-66.
+std::vector<std::pair<PresburgerRelation, GeneratingFunction>>
+mlir::presburger::detail::polytopeGeneratingFunction(PolyhedronH poly) {
+  unsigned numVars = poly.getNumRangeVars();
+  unsigned numParams = poly.getNumSymbolVars();
+  unsigned numIneqs = poly.getNumInequalities();
+
+  // The generating function of the polytope is computed as a set of generating
+  // functions, each one associated with a region in parameter space (chamber).
+  std::vector<std::pair<PresburgerRelation, GeneratingFunction>> gf({});
+
+  // The active region will be defined as activeRegionCoeffs @ p +
+  // activeRegionConstant ≥ 0. The active region is a polyhedron in parameter
+  // space.
+  FracMatrix activeRegion(numIneqs - numVars, numParams + 1);
+
+  // These vectors store lists of
+  // subsets of inequalities,
+  // the vertices corresponding to them, and
+  // the active regions of the vertices, in order.
+  std::vector<IntMatrix> subsets;
+  std::vector<ParamPoint> vertices;
+  std::vector<PresburgerRelation> activeRegions;
+
+  FracMatrix a2(numIneqs - numVars, numVars);
+  FracMatrix b2c2(numIneqs - numVars, numParams + 1);
+
+  // We iterate over all subsets of inequalities with cardinality numVars,
+  // using bitsets up to 2^numIneqs to enumerate.
+  for (std::bitset<16> indicator(((1ul << numVars) - 1ul)
+                                 << (numIneqs - numVars));
+       indicator.to_ulong() <=
+       ((1ul << numVars) - 1ul)
+           << (numIneqs - numVars); // d 1's followed by n-numVars 0's
+       indicator = std::bitset<16>(indicator.to_ulong() - 1)) {
+
+    if (indicator.count() != numVars)
+      continue;
+
+    // Collect the inequalities corresponding to the bits which are set.
+    IntMatrix subset(numVars, numVars + numParams + 1);
+    unsigned j1 = 0, j2 = 0;
+    for (unsigned i = 0; i < numIneqs; i++)
+      if (indicator.test(i))
+        subset.setRow(j1++, poly.getInequality(i));
+
+      else {
+        // All other inequalities are stored in a2 and b2c2.
+        // These are column-wise splits of the inequalities;
+        // a2 stores the coefficients of the variables, and
+        // b2c2 stores the coefficients of the parameters and the constant term.
+        for (unsigned k = 0; k < numVars; k++)
+          a2(j2, k) = poly.atIneq(i, k);
+        for (unsigned k = numVars; k < numVars + numParams + 1; k++)
+          b2c2(j2, k - numVars) = poly.atIneq(i, k);
+        j2++;
+      }
+
+    // Find the vertex, if any, corresponding to the current subset of
+    // inequalities.
+    std::optional<ParamPoint> vertex = findVertex(subset); // d x (p+1)
+
+    if (vertex == std::nullopt)
+      continue;
+    // If this subset corresponds to a vertex, store it.
+    vertices.push_back(*vertex);
+    subsets.push_back(subset);
+
+    // Let the current vertex be [X | y], where
+    // X represents the coefficients of the parameters and
+    // y represents the constant term.
+
+    // The region (in parameter space) where this vertex is active is given
+    // by substituting the vertex into the *remaining* inequalities of the
+    // polytope (those which were not collected into `subset`), i.e.,
+    // [A2 | B2 | c2].
+    // Thus, the coefficients of the parameters after substitution become
+    // (A2 • X + B2)
+    // and the constant terms become
+    // (A2 • y + c2).
+    // The region is therefore given by
+    // (A2 • X + B2) p + (A2 • y + c2) ≥ 0
+    // This is equivalent to A2 • [X | y] + [B2 | c2]
+    // Thus we premultiply [X | y] with each row of A2
+    // and add each row of [B2 | c2].
+    for (unsigned i = 0; i < numIneqs - numVars; i++) {
+      activeRegion.setRow(i, (*vertex).preMultiplyWithRow(a2.getRow(i)));
+      activeRegion.addToRow(i, b2c2.getRow(i), 1);
+    }
+
+    // We convert the representation of the active region to an integers-only
+    // form so as to store it as an PresburgerRelation.
+    // We do this by taking the LCM of the denominators of all the coefficients
+    // and multiplying by it throughout.
+    IntMatrix activeRegionNorm = IntMatrix(numIneqs - numVars, numParams + 1);
+    IntegerRelation activeRegionRel =
+        IntegerRelation(PresburgerSpace::getRelationSpace(0, numParams, 0, 0));
+    MPInt lcmDenoms = MPInt(1);
+    for (unsigned i = 0; i < numIneqs - numVars; i++) {
+      for (unsigned j = 0; j < numParams + 1; j++)
+        lcmDenoms = lcm(lcmDenoms, activeRegion(i, j).den);
+      for (unsigned j = 0; j < numParams + 1; j++)
+        activeRegionNorm(i, j) =
+            (activeRegion(i, j) * lcmDenoms).getAsInteger();
+
+      activeRegionRel.addInequality(activeRegionNorm.getRow(i));
+    }
+
+    activeRegions.push_back(PresburgerRelation(activeRegionRel));
+  }
+
+  // Now, we use Clauss-Loechner decomposition to identify regions in parameter
+  // space where each vertex is active. These regions (chambers) have the
+  // property that no two of them have a full-dimensional intersection, i.e.,
+  // they may share "faces" or "edges", but their intersection can only have
+  // up to numVars-1 dimensions.
+
+  // We maintain a list of regions and their associated vertex sets,
+  // initialized with the first vertex and its corresponding activity region.
+  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> chambers = {
+      std::make_pair(activeRegions[0], std::vector({0u}))};
+  // Note that instead of storing lists of actual vertices, we store lists
+  // of indices. Thus the set {2, 3, 4} represents the vertex set
+  // {vertices[2], vertices[3], vertices[4]}.
+
+  std::vector<std::pair<PresburgerRelation, std::vector<unsigned>>> newChambers;
+
+  // We iterate over the vertex set.
+  // For each vertex v_j and its activity region R_j,
+  // we examine all the current chambers R_i.
+  // If R_j has a full-dimensional intersection with an existing chamber R_i,
+  // then that chamber is replaced by two new ones:
+  // 1. the intersection R_i \cap R_j, where v_j is active;
+  // 2. the difference R_i - R_j, where v_j is inactive.
+  // Once we have examined all R_i, we add a final chamber
+  // R_j - (union of all existing chambers),
+  // in which only v_j is active.
+  for (unsigned j = 1; j < vertices.size(); j++) {
+    newChambers.clear();
+
+    PresburgerRelation r_j = activeRegions[j];
+    ParamPoint v_j = vertices[j];
+
+    for (unsigned i = 0; i < chambers.size(); i++) {
+      auto [r_i, v_i] = chambers[i];
+
+      // First, we check if the intersection of R_j and R_i.
+      // It is a disjoint union of convex regions in the parameter space,
+      // and so we know that it is full-dimensional if any of the disjuncts
+      // is full-dimensional.
+      PresburgerRelation intersection = r_i.intersect(r_j);
+      bool isFullDim = false;
+      for (auto disjunct : intersection.getAllDisjuncts())
+        if (disjunct.isFullDim()) {
+          isFullDim = true;
+          break;
+        }
+      isFullDim = (numParams == 0) || isFullDim;
+
+      // If the intersection is not full-dimensional, we do not modify
+      // the chamber list.
+      if (!isFullDim)
+        newChambers.push_back(chambers[i]);
+      else {
+        // If it is, we add the intersection and the difference as new chambers.
+        PresburgerRelation subtraction = r_i.subtract(r_j);
+        newChambers.push_back(std::make_pair(subtraction, v_i));
+
+        v_i.push_back(j);
+        newChambers.push_back(std::make_pair(intersection, v_i));
+      }
+    }
+
+    // Finally we compute the chamber where only v_j is active by subtracting
+    // all existing chambers from R_j.
+    for (auto chamber : newChambers)
+      r_j = r_j.subtract(chamber.first);
+    newChambers.push_back(std::make_pair(r_j, std::vector({j})));
+
+    // We filter `chambers` to remove empty regions.
+    chambers.clear();
+    for (auto chamber : newChambers) {
+      bool empty = true;
+      for (auto disjunct : chamber.first.getAllDisjuncts())
+        if (!disjunct.isEmpty()) {
+          empty = false;
+          break;
+        }
+      if (!empty)
+        chambers.push_back(chamber);
+    }
+  }
+
+  // Now, we compute the generating function. For each chamber, we iterate over
+  // the vertices active in it, and compute the generating function for each
+  // of them. The sum of these generating functions is the GF corresponding to
+  // the entire polytope.
+  SmallVector<MPInt> ineq(numVars + 1);
+  for (auto chamber : chambers) {
+    GeneratingFunction chamberGf(numParams, {}, {}, {});
+    for (unsigned i : chamber.second) {
+      // We collect the inequalities corresponding to each vertex.
+      // We only need the coefficients of the variables (NOT the parameters)
+      // as the generating function only depends on these.
+      ConeH tgtCone = defineHRep(numVars);
+      for (unsigned j = 0; j < numVars; j++) {
+        for (unsigned k = 0; k < numVars; k++)
+          ineq[k] = subsets[i](j, k);
+        ineq[numVars] = subsets[i](j, numVars + numParams);
+        tgtCone.addInequality(ineq);
+      }
+      // We assume that the tangent cone is unimodular.
+      SmallVector<std::pair<int, ConeH>, 4> unimodCones = {
+          std::make_pair(1, tgtCone)};
+      for (auto signedCone : unimodCones)
+        chamberGf =
+            chamberGf + unimodularConeGeneratingFunction(
+                            vertices[i], signedCone.first, signedCone.second);
+    }
+    gf.push_back(std::make_pair(chamber.first, chamberGf));
+  }
+  return gf;
+}
+
 /// We use an iterative procedure to find a vector not orthogonal
 /// to a given set, ignoring the null vectors.
 /// Let the inputs be {x_1, ..., x_k}, all vectors of length n.
diff --git a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
index 7d2a63d17676f57..b05075b2892ace0 100644
--- a/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
+++ b/mlir/lib/Analysis/Presburger/IntegerRelation.cpp
@@ -2498,6 +2498,48 @@ void IntegerRelation::printSpace(raw_ostream &os) const {
   os << getNumConstraints() << " constraints\n";
 }
 
+void IntegerRelation::removeTrivialEqualities() {
+  bool flag;
+  for (unsigned i = 0; i < getNumEqualities(); i++) {
+    flag = true;
+    for (unsigned j = 0; j < getNumVars() + 1; j++)
+      if (atEq(i, j) != 0)
+        flag = false;
+    if (flag)
+      removeEquality(i);
+  }
+}
+
+bool IntegerRelation::isFullDim() {
+  removeTrivialEqualities();
+
+  // If there is a non-trivial equality, the space cannot be full-dimensional.
+  if (getNumEqualities() > 0)
+    return false;
+
+  // If along the direction of any of the inequalities, the upper and lower
+  // optima are the same, then the region is not full-dimensional.
+  Simplex simplex(*this);
+  for (unsigned i = 0; i < getNumInequalities(); i++) {
+    auto ineq = inequalities.getRow(i);
+    auto upOpt = simplex.computeOptimum(Simplex::Direction::Up, ineq);
+    auto downOpt = simplex.computeOptimum(Simplex::Direction::Down, ineq);
+
+    if (upOpt.getKind() == OptimumKind::Unbounded ||
+        downOpt.getKind() == OptimumKind::Unbounded)
+      continue;
+
+    // Check if the upper and lower optima are equal.
+    if (upOpt.getKind() == OptimumKind::Bounded &&
+        downOpt.getKind() == OptimumKind::Bounded && *upOpt == *downOpt)
+      return false;
+  }
+  // If none of the inequalities were such that the upper and lower optima
+  // along their direction were equal, then we conclude that the region is full
+  // dimensional.
+  return true;
+}
+
 void IntegerRelation::print(raw_ostream &os) const {
   assert(hasConsistentState());
   printSpace(os);

``````````

</details>


https://github.com/llvm/llvm-project/pull/78987


More information about the Mlir-commits mailing list