[Mlir-commits] [mlir] [MLIR][Presburger][WIP] Implement vertex enumeration and chamber decomposition for polytope generating function computation. (PR #78987)
Kunwar Grover
llvmlistbot at llvm.org
Mon Jan 22 08:01:38 PST 2024
================
@@ -147,6 +148,319 @@ 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 to enumerate.
+ // The largest possible bitset that corresponds to such a subset can be
+ // written as numVar 1's followed by (numIneqs - numVars) 0's.
+ unsigned upperBound = ((1ul << numVars) - 1ul) << (numIneqs - numVars);
+ for (std::bitset<16> indicator(((1ul << numVars) - 1ul)
+ << (numIneqs - numVars));
+ 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.
+ 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, e = vertices.size(); j < e; j++) {
+ newChambers.clear();
+
+ PresburgerRelation r_j = activeRegions[j];
----------------
Groverkss wrote:
I don't think these variable names fit mlir style.
https://github.com/llvm/llvm-project/pull/78987
More information about the Mlir-commits
mailing list