[Mlir-commits] [mlir] [MLIR][Presburger] Implement vertex enumeration and chamber decomposition for polytope generating function computation. (PR #78987)
Arjun P
llvmlistbot at llvm.org
Thu Jan 25 15:49:17 PST 2024
================
@@ -147,6 +149,314 @@ 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::solveParametricEquations(FracMatrix equations) {
+ // equations is a d x (d + p + 1) matrix.
+ // Each row represents an equation.
+ unsigned d = equations.getNumRows();
+ unsigned numCols = equations.getNumColumns();
+
+ // If the determinant is zero, there is no unique solution.
+ // Thus we return null.
+ if (FracMatrix(equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/d - 1,
+ /*fromColumn=*/0,
+ /*toColumn=*/d - 1))
+ .determinant() == 0)
+ return std::nullopt;
+
+ // Perform row operations to make each column all zeros except for the
+ // diagonal element, which is made to be one.
+ for (unsigned i = 0; i < d; ++i) {
+ // First ensure that the diagonal element is nonzero, by swapping
+ // it with a row that is non-zero at column i.
+ if (equations(i, i) != 0)
+ continue;
+ for (unsigned j = i + 1; j < d; ++j) {
+ if (equations(j, i) == 0)
+ continue;
+ equations.swapRows(j, i);
+ break;
+ }
+
+ Fraction diagElement = equations(i, i);
+
+ // Apply row operations to make all elements except the diagonal to zero.
+ for (unsigned j = 0; j < d; ++j) {
+ if (i == j)
+ continue;
+ if (equations(j, i) == 0)
+ continue;
+ // Apply row operations to make element (j, i) zero by subtracting the
+ // ith row, appropriately scaled.
+ Fraction currentElement = equations(j, i);
+ equations.addToRow(/*sourceRow=*/i, /*targetRow=*/j,
+ /*scale=*/-currentElement / diagElement);
+ }
+ }
+
+ // Rescale diagonal elements to 1.
+ for (unsigned i = 0; i < d; ++i)
+ equations.scaleRow(i, 1 / equations(i, i));
+
+ // Now we have reduced the equations to the form
+ // x_i + b_1' m_1 + ... + b_p' m_p + c' = 0
+ // i.e. each variable appears exactly once in the system, and has coefficient
+ // one.
+ // Thus we have
+ // x_i = - b_1' m_1 - ... - b_p' m_p - c
+ // and so we return the negation of the last p + 1 columns of the matrix.
+ // We copy these columns and return them.
+ ParamPoint vertex =
+ equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/d - 1,
+ /*fromColumn=*/d, /*toColumn=*/numCols - 1);
+ vertex.negateMatrix();
+ return vertex;
+}
+
+/// This is an implementation of the Clauss-Loechner algorithm for chamber
+/// decomposition.
+/// We maintain a list of pairwise disjoint chambers and their generating
+/// functions. We iterate over the list of regions, each time adding the
+/// generating function to the chambers where the corresponding vertex is
+/// active and separating the chambers where it is not.
+///
+/// Given the region each vertex is active in, for each subset of vertices,
+/// the region that precisely this subset is in, is the intersection of the
+/// regions that these are active in, intersected with the complements of the
+/// remaining regions. The generating function for this region is the sum of
+/// generating functions for the vertices' tangent cones.
+std::vector<std::pair<PresburgerSet, GeneratingFunction>>
+mlir::presburger::detail::computeChamberDecomposition(
+ unsigned numSymbols, ArrayRef<std::pair<PresburgerSet, GeneratingFunction>>
+ regionsAndGeneratingFunctions) {
+ assert(regionsAndGeneratingFunctions.size() > 0 &&
+ "there must be at least one chamber!");
+ // We maintain a list of regions and their associated generating function
+ // initialized with the universe and the empty generating function.
+ std::vector<std::pair<PresburgerSet, GeneratingFunction>> chambers = {
+ {PresburgerSet::getUniverse(
+ PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0)),
+ GeneratingFunction(numSymbols, {}, {}, {})}};
+
+ // We iterate over the region list.
+ // For each activity region R_j (corresponding to a vertex v_j, whose
+ // generating function is gf_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 the generating function is
+ // gf_i + gf_j.
+ // 2. the difference R_i - R_j, where v_j is inactive and the generating
+ // function is gf_i.
+
+ // At each step, we define a new chamber list after considering vertex v_j,
+ // and its generating function, replacing and appending chambers as
+ // discussed above.
+ // The loop has the invariant that the union over all the chambers gives the
+ // universe at every step (modulo lower-dimensional spaces).
+ for (const auto &[region, generatingFunction] :
+ regionsAndGeneratingFunctions) {
+ std::vector<std::pair<PresburgerSet, GeneratingFunction>> newChambers;
+
+ for (auto [currentRegion, currentGeneratingFunction] : chambers) {
+ // First, we check if the intersection of R_j and R_i.
+ PresburgerSet intersection = currentRegion.intersect(region);
+
+ // If the intersection is not full-dimensional, we do not modify
+ // the chamber list.
+ if (!intersection.isFullDim()) {
+ newChambers.emplace_back(currentRegion, currentGeneratingFunction);
+ continue;
+ }
+
+ // If it is, we add the intersection as a chamber where this vertex is
+ // active.
+ newChambers.emplace_back(intersection,
+ currentGeneratingFunction + generatingFunction);
+
+ // We also add the difference if it is full-dimensional.
+ PresburgerSet difference = currentRegion.subtract(region);
+ if (difference.isFullDim())
+ newChambers.emplace_back(difference, currentGeneratingFunction);
+ }
+ chambers = std::move(newChambers);
+ }
+
+ return chambers;
+}
+
+/// 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 satisfiability.
+/// 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<PresburgerSet, GeneratingFunction>>
+mlir::presburger::detail::computePolytopeGeneratingFunction(PolyhedronH poly) {
+ unsigned numVars = poly.getNumRangeVars();
+ unsigned numSymbols = poly.getNumSymbolVars();
+ unsigned numIneqs = poly.getNumInequalities();
+
+ // These vectors store lists of
+ // subsets of inequalities,
+ // the vertices corresponding to them,
+ // the active regions of the vertices, and the generating functions of
+ // the tangent cones of the vertices, in order.
+ // std::vector<IntMatrix> subsets;
+ std::vector<ParamPoint> vertices;
+ std::vector<std::pair<PresburgerSet, GeneratingFunction>>
+ regionsAndGeneratingFunctions;
+
+ FracMatrix a2(numIneqs - numVars, numVars);
+ FracMatrix b2c2(numIneqs - numVars, numSymbols + 1);
+
+ // We iterate over all subsets of inequalities with cardinality numVars,
+ // using bitsets of numIneqs bits to enumerate.
+ // For a given set of numIneqs bits, we consider a subset which contains
+ // the i'th inequality if the i'th bit in the bitset is set.
+ // The largest possible bitset that corresponds to such a subset can be
+ // written as numVar 1's followed by (numIneqs - numVars) 0's.
+ // We start with this and count downwards (in binary), considering only
+ // those numbers whose binary representation has numVars 1's and the
+ // rest 0's.
+ unsigned upperBound = ((1ul << numVars) - 1ul) << (numIneqs - numVars);
+ for (std::bitset<16> indicator(upperBound);
+ indicator.to_ulong() <= upperBound;
+ indicator = std::bitset<16>(indicator.to_ulong() - 1)) {
+
+ if (indicator.count() != numVars)
+ continue;
+
+ // Collect the inequalities corresponding to the bits which are set
+ // and the remaining ones.
+ auto [subset, remainder] = poly.getInequalities().splitByBitset(indicator);
+ // 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.
+ a2 = FracMatrix(
+ remainder.getSubMatrix(0, numIneqs - numVars - 1, 0, numVars - 1));
+ b2c2 = FracMatrix(remainder.getSubMatrix(0, numIneqs - numVars - 1, numVars,
+ numVars + numSymbols));
+
+ // Find the vertex, if any, corresponding to the current subset of
+ // inequalities.
+ std::optional<ParamPoint> vertex =
+ solveParametricEquations(FracMatrix(subset)); // d x (p+1)
+
+ if (!vertex)
+ 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].
+ FracMatrix activeRegion(numIneqs - numVars, numSymbols + 1);
+ for (unsigned i = 0; i < numIneqs - numVars; i++) {
+ activeRegion.setRow(i, vertex->preMultiplyWithRow(a2.getRow(i)));
+ activeRegion.addToRow(i, b2c2.getRow(i), 1);
+ }
+
+ // If any row of activeRegion is all-zero, that means that the
+ // corresponding inequality in `remaining` is *also* satisfied by the
+ // vertex for any values of the parameters. Thus we remove it from
+ // activeRegion and add that inequality to `subset`.
+ // Note that if the row is not all-zero, it is still possible for the
+ // inequality to be satisfied by the corresponding vertex *in some subset of
+ // the parameter space*. However, since this subset is defined by an
+ // equality, it is not full-dimensional and we can therefore ignore it.
+ for (int i = activeRegion.getNumRows() - 1; i >= 0; --i) {
+ if (llvm::any_of(activeRegion.getRow(i),
+ [&](Fraction f) { return f != 0; }))
+ continue;
+ subset.appendExtraRow(remainder.getRow(i));
+ activeRegion.removeRow(i);
+ }
+
+ // We convert the representation of the active region to an integers-only
+ // form so as to store it as a PresburgerSet.
+ IntMatrix activeRegionNorm = activeRegion.normalizeRows();
+ IntegerPolyhedron activeRegionRel = IntegerPolyhedron(
+ PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0));
+ for (unsigned i = 0, e = activeRegion.getNumRows(); i < e; ++i)
+ activeRegionRel.addInequality(activeRegionNorm.getRow(i));
+
+ // Now, we compute the generating function at this vertex.
+ // We collect the inequalities corresponding to each vertex to compute
+ // the tangent cone at that vertex.
+ SmallVector<MPInt> ineq(numVars + 1);
+
+ // We only need the coefficients of the variables (NOT the parameters)
+ // as the generating function only depends on these.
+ // We translate the cones to be pointed at the origin by making the
+ // constant terms zero.
+ ConeH tangentCone = defineHRep(numVars);
+ for (unsigned j = 0, e = subset.getNumRows(); j < e; ++j) {
+ for (unsigned k = 0; k < numVars; ++k)
+ ineq[k] = subset(j, k);
----------------
Superty wrote:
by this point I don't remember what the variable name `subset` refers to, please use a name like activeIneqs, supportingIneqs, tightIneqs, definingIneqs, or whatever
https://github.com/llvm/llvm-project/pull/78987
More information about the Mlir-commits
mailing list