[MLIR][Presburger] Add applyDomain/Range to IntegerRelation

[Groverkss]

Summary

This patch adds support for applying a relation on domain/range of a relation.

[arjunp]

I think we should add a compose member function, it being a natural opreation on relations. You can move applyRange's documentation to compose and just mention in applyRange that it's the same as compose, but provided for uniformity with applyDomain.

mlir/include/mlir/Analysis/Presburger/IntegerRelation.h

525–527

It might be simpler to just say that this returns R2.inverse().compose(R1), and then perhaps provide a concrete example.

532–534

You can just write that it returns R1.compose(R2) and mention that this is provided for uniformity with applyDomain.

[Groverkss]

• Rebase
• Address Arjun's comments

[Groverkss]

mlir/include/mlir/Analysis/Presburger/IntegerRelation.h

525–527

I feel the R2.invese().compose() makes it harder to understand than providing a proper mathematical description.

[arjunp]

I think you can remove the backslashes on exists and in, unless doxygen can somehow render them, it seems unnecessary.

mlir/include/mlir/Analysis/Presburger/IntegerRelation.h

525–527

R2.inverse().compose() is in fact a proper mathematical description. For me personally, a concrete example was more useful to understand what the point of this is, and then for a formal specification the inverse compose thing will show what the general version of that example is.

[Groverkss]

Address Arjun's comments

[arjunp]

mlir/include/mlir/Analysis/Presburger/IntegerRelation.h

524

What do you mean by R1(R2) here? If R3 is A -> C then R3(x) is R2(R1(x)). It seems the standard notation here is R1;R2 (https://en.wikipedia.org/wiki/Composition_of_relations).

526

527

528

533

534

read the original as (i->j), k at first

536

I would prefer this example, it has less "noise" and also shows a little bit happens when the other variables depend on the projected variable

R1 : i -> j : (0 <= i < 2, j = i)
R2: i -> k : (k = i floordiv 2)

R3: k -> j : (0 <= k < 1, 2k <=  j <= 2k + 1)

R1 = {(0, 0), (1, 1)}. R2 maps both 0 and 1 to 0. So R3 = {(0, 0), (0, 1)}

541–548

Not sure if an example is needed for composition but you can do e.g. R2 : j -> k : k = floordiv 2 and R1 from above

543

548

mlir/lib/Analysis/Presburger/IntegerRelation.cpp

2153–2161

I think the current implementation should be thought of as converting R1 to A -> (B x C), converting R2 to B x C, and intersecting

[Groverkss]

mlir/include/mlir/Analysis/Presburger/IntegerRelation.h

524

I meant it like we represent function composition f(g(x)). But this also works for me.

mlir/lib/Analysis/Presburger/IntegerRelation.cpp

2153–2161

good catch!

[Groverkss]

Address Arjun's comments

[arjunp]

mlir/include/mlir/Analysis/Presburger/IntegerRelation.h

524

Right, I mean that we are doing R2(R1(x)) then, not R1(R2).

I guess you can also add that it R3)(x) = R2(R1(x)), might be more understandable

537

539

[Groverkss]

Address Arjun's comments and rebase

mlir/include/mlir/Analysis/Presburger/IntegerRelation.h

public:
   /// to perform much better in the average case. If V is the number of vertices
   /// in the polytope and C is the number of constraints, the algorithm takes
   /// O(VC) time.
   void removeRedundantConstraints();
 
   void removeDuplicateDivs();
 
   /// Converts identifiers of kind srcKind in the range [idStart, idLimit) to
-  /// variables of kind dstKind and placed after all the other variables of kind
-  /// dstKind. The internal ordering among the moved variables is preserved.
+  /// variables of kind dstKind. If `pos` is given, the variables are placed at
+  /// position `pos` of dstKind, otherwise they are placed after all the other
+  /// variables of kind dstKind. The internal ordering among the moved variables
+  /// is preserved.
   void convertIdKind(IdKind srcKind, unsigned idStart, unsigned idLimit,
-                     IdKind dstKind);
+                     IdKind dstKind, unsigned pos);
+  void convertIdKind(IdKind srcKind, unsigned idStart, unsigned idLimit,
+                     IdKind dstKind) {
+    convertIdKind(srcKind, idStart, idLimit, dstKind, getNumIdKind(dstKind));
+  }
   void convertToLocal(IdKind kind, unsigned idStart, unsigned idLimit) {
     convertIdKind(kind, idStart, idLimit, IdKind::Local);
   }
 
   /// Adds additional local ids to the sets such that they both have the union
   /// of the local ids in each set, without changing the set of points that
   /// lie in `this` and `other`.
   ///
public:
   void intersectDomain(const IntegerPolyhedron &poly);
 
   /// Intersect the given `poly` with the range.
   void intersectRange(const IntegerPolyhedron &poly);
 
   /// Get inverse of this relation.
   void inverse();
 
+  /// Given a relation `rel`, apply the relation to the domain of this relation.
+  ///
+  /// Formally, let the relation `this` be R1: A -> B, and `rel` be
+  /// R2: A -> C, returns a relation R3: C -> B, such that, a point (c, b) \in
+  /// R3, iff \exists a, (a, b) \in R1, (a, c) \in R2.
+  void applyDomain(const IntegerRelation &rel);
+
+  /// Given a relation `rel`, apply the relation to the range of this relation.
+  ///
+  /// Formally, let the relation `this` be R1: A -> B, and `rel` be
+  /// R2: B -> C, returns a relation R3: A -> C, such that, a point (a, c) \in
+  /// R3, iff \exists b, (a, b) \in R1, (b, c) \in R2.
+  void applyRange(const IntegerRelation &rel);
+
   void print(raw_ostream &os) const;
   void dump() const;
 
 protected:
   /// Checks all rows of equality/inequality constraints for trivial
   /// contradictions (for example: 1 == 0, 0 >= 1), which may have surfaced
   /// after elimination. Returns true if an invalid constraint is found;
   /// false otherwise.
   bool hasInvalidConstraint() const;
 
   /// Returns the constant lower bound bound if isLower is true, and the upper
   /// bound if isLower is false.
   template <bool isLower>
   Optional<int64_t> computeConstantLowerOrUpperBound(unsigned pos);
 
   /// Eliminates a single identifier at `position` from equality and inequality
   /// constraints. Returns `success` if the identifier was eliminated, and
   /// `failure` otherwise.
   inline LogicalResult gaussianEliminateId(unsigned position) {
     return success(gaussianEliminateIds(position, position + 1) == 1);

mlir/lib/Analysis/Presburger/IntegerRelation.cpp

while (true) {
 
     // Remove the ith equality and the found local variable.
     removeId(j);
     removeEquality(i);
   }
 }
 
 void IntegerRelation::convertIdKind(IdKind srcKind, unsigned idStart,
-                                    unsigned idLimit, IdKind dstKind) {
+                                    unsigned idLimit, IdKind dstKind,
+                                    unsigned pos) {
   assert(idLimit <= getNumIdKind(srcKind) && "Invalid id range");
 
   if (idStart >= idLimit)
     return;
 
   // Append new local variables corresponding to the dimensions to be converted.
-  unsigned newIdsBegin = getIdKindEnd(dstKind);
   unsigned convertCount = idLimit - idStart;
-  appendId(dstKind, convertCount);
+  unsigned newIdsBegin = insertId(dstKind, pos, convertCount);
 
   // Swap the new local variables with dimensions.
   //
   // Essentially, this moves the information corresponding to the specified ids
   // of kind `srcKind` to the `convertCount` newly created ids of kind
   // `dstKind`. In particular, this moves the columns in the constraint
   // matrices, and zeros out the initially occupied columns (because the newly
   // created ids we're swapping with were zero-initialized).

 }
 
 void IntegerRelation::inverse() {
   unsigned numRangeIds = getNumIdKind(IdKind::Range);
   convertIdKind(IdKind::Domain, 0, getIdKindEnd(IdKind::Domain), IdKind::Range);
   convertIdKind(IdKind::Range, 0, numRangeIds, IdKind::Domain);
 }
 
+void IntegerRelation::applyRange(const IntegerRelation &rel) {
+  assert(getRangeSet().getSpace().isCompatible(rel.getDomainSet().getSpace()) &&
+         "Range of `this` should be compatible with Domain of `rel`");
+
+  IntegerRelation copyRel = rel;
+
+  // Let relation `this` be R1: A -> B, and `rel` be R2: B -> C.
+  // We convert R1 to A -> (B -> C) then intersect the range of R1 with R2.
+  // After this, we get R1: A -> C, by projecting out B.
+  unsigned numBIds = getNumRangeIds();
+
+  // Convert R1 from A -> B to A -> (B -> C).
+  appendId(IdKind::Range, copyRel.getNumRangeIds());
+
+  // Convert R2 to a set for intersection.
+  // TODO: Using nested spaces here would help, since we could directly
+  // intersect the range with another relation.
+  copyRel.convertIdKind(IdKind::Domain, 0, numBIds, IdKind::Range, 0);
+
+  // Intersect R2 to range of R1.
+  intersectRange(IntegerPolyhedron(copyRel));
+
+  // Project out B in R1.
+  convertIdKind(IdKind::Range, 0, numBIds, IdKind::Local);
+}
+
+void IntegerRelation::applyDomain(const IntegerRelation &rel) {
+  inverse();
+  applyRange(rel);
+  inverse();
+}
+
 void IntegerRelation::printSpace(raw_ostream &os) const {
   space.print(os);
   os << getNumConstraints() << " constraints\n";
 }
 
 void IntegerRelation::print(raw_ostream &os) const {
   assert(hasConsistentState());
   printSpace(os);

mlir/unittests/Analysis/Presburger/IntegerRelationTest.cpp

IntegerRelation expectedRel = parseRelationFromSet(
         ">= 0, x + y + z floordiv 7 == 0, y >= 0, M - y - 1 >= 0, y + z == 0)",
         1);
 
     IntegerRelation copyRel = rel;
     copyRel.intersectRange(poly);
     EXPECT_TRUE(copyRel.isEqual(expectedRel));
   }
 }
+
+TEST(IntegerRelationTest, applyDomainAndRange) {
+
+  {
+    IntegerRelation map1 = parseRelationFromSet(
+        "(x, y, a, b)[N] : (a - x - N == 0, b - y + N == 0)", 2);
+    IntegerRelation map2 =
+        parseRelationFromSet("(x, y, a)[N] : (a - x - y == 0)", 2);
+
+    map1.applyRange(map2);
+
+    IntegerRelation map3 =
+        parseRelationFromSet("(x, y, a)[N] : (a - x - y == 0)", 2);
+
+    EXPECT_TRUE(map1.isEqual(map3));
+  }
+
+  {
+    IntegerRelation map1 = parseRelationFromSet(
+        "(x, y, a, b)[N] : (a - x + N == 0, b - y - N == 0)", 2);
+    IntegerRelation map2 =
+        parseRelationFromSet("(x, y, a, b)[N] : (a - N == 0, b - N == 0)", 2);
+
+    IntegerRelation map3 =
+        parseRelationFromSet("(x, y, a, b)[N] : (x - N == 0, y - N == 0)", 2);
+
+    map1.applyDomain(map2);
+
+    EXPECT_TRUE(map1.isEqual(map3));
+  }
+}