feat: use preprocessor to infer simple equality bounds in the solver …

…ensuring correct explanation

dlinear/solver/QsoptexTheorySolver.cpp

@@ -276,7 +276,7 @@ void QsoptexTheorySolver::UpdateExplanation([[maybe_unused]] LiteralSet &explana
     // For each free variable in the literal, add their bounds to the explanation
     for (const auto &col_var : predicate_abstractor_[lit.var].GetFreeVariables()) {
       const int theory_col = var_to_theory_col_.at(col_var.get_id());
-      TheoryBoundsToExplanation(theory_col, true, explanation);
+      TheoryBoundsToExplanation(theory_col, explanation);
     }
   }
 }

dlinear/solver/SoplexTheorySolver.cpp

@@ -274,7 +274,7 @@ void SoplexTheorySolver::UpdateExplanation(LiteralSet &explanation) {
     // Add all the active bounds for the free variables in the row to the explanation
     for (const auto &col_var : predicate_abstractor_[var].GetFreeVariables()) {
       const int theory_col = var_to_theory_col_.at(col_var.get_id());
-      TheoryBoundsToExplanation(theory_col, true, explanation);
+      TheoryBoundsToExplanation(theory_col, explanation);
     }
   }
 }
@@ -289,6 +289,14 @@ void SoplexTheorySolver::EnableSPXVarBound() {
   for (int theory_col = 0; theory_col < static_cast<int>(theory_bounds_.size()); theory_col++) {
     spx_.changeBoundsRational(theory_col, theory_bounds_[theory_col].active_lower_bound().get_mpq_t(),
                               theory_bounds_[theory_col].active_upper_bound().get_mpq_t());
+    if (theory_col == static_cast<int>(theory_bounds_.size() - 1)) continue;
+    DLINEAR_TRACE_FMT("EnableSPXVarBound: {} = [{}, {}]", theory_col_to_var_[theory_col],
+                      theory_bounds_[theory_col].active_lower_bound(), theory_bounds_[theory_col].active_upper_bound());
+  }
+  for (const auto &[var, value] : preprocessor_.env()) {
+    const int theory_col = var_to_theory_col_.at(var.get_id());
+    DLINEAR_TRACE_FMT("EnableSPXVarBound: {} = {} (theory_col: {})", var, value, theory_col);
+    spx_.changeBoundsRational(theory_col, value.get_mpq_t(), value.get_mpq_t());
   }
 }
 

dlinear/solver/TheorySolver.cpp

@@ -138,13 +138,11 @@ void TheorySolver::TheoryBoundsToExplanations(Violation violation, int theory_ro
     explanations.insert(explanation);
   }
 }
-void TheorySolver::TheoryBoundsToExplanation(const int theory_col, const bool active, LiteralSet &explanation) const {
-  if (active) {
-    theory_bounds_.at(theory_col).GetActiveExplanation(theory_row_to_lit_, explanation);
+void TheorySolver::TheoryBoundsToExplanation(const int theory_col, LiteralSet &explanation) {
+  if (!config_.disable_theory_preprocessor() && preprocessor_.env().contains(theory_col_to_var_[theory_col])) {
+    preprocessor_.GetActiveExplanation(theory_col, explanation);
   } else {
-    for (const auto &bound : theory_bounds_[theory_col].bounds()) {
-      explanation.insert(theory_row_to_lit_[bound.idx]);
-    }
+    theory_bounds_.at(theory_col).GetActiveExplanation(theory_row_to_lit_, explanation);
   }
 }
 void TheorySolver::TheoryBoundsToExplanation(const int theory_col, const mpq_class &value,
@@ -178,9 +176,12 @@ void TheorySolver::TheoryBoundsToBoundIdxs(const int theory_col, const mpq_class
   for (auto it = theory_bounds_[theory_col].GetActiveBound(value); it; ++it) bound_idxs.insert(it->idx);
 }
 void TheorySolver::TheoryBoundsToBoundIdxs(int theory_col, const TheorySolver::BoundViolationType type,
-                                           std::set<int> &bound_idxs) const {
+                                           std::set<int> &bound_idxs) {
   if (type == BoundViolationType::NO_BOUND_VIOLATION) return;
   const TheorySolverBoundVector &bound = theory_bounds_[theory_col];
+  if (!config_.disable_theory_preprocessor() && preprocessor_.env().contains(theory_col_to_var_[theory_col])) {
+    return preprocessor_.GetActiveBoundIdxs(theory_col, bound_idxs);
+  }
   return TheoryBoundsToBoundIdxs(
       theory_col,
       type == BoundViolationType::LOWER_BOUND_VIOLATION ? bound.active_lower_bound() : bound.active_upper_bound(),

dlinear/solver/TheorySolver.h

@@ -276,10 +276,9 @@ class TheorySolver {
   /**
    * Gather the bounds of the @p theory_col and produce an explanation for the SAT solver.
    * @param theory_col column of the theory solver the bounds are associated with
-   * @param active whether to only consider the active bounds (true) or include the inactive ones as well (false)
    * @param[out] explanation set of literals that correspond to the conflicting bounds
    */
-  void TheoryBoundsToExplanation(int theory_col, bool active, LiteralSet &explanation) const;
+  void TheoryBoundsToExplanation(int theory_col, LiteralSet &explanation);
   /**
    * Gather the bounds that enforced @p value on @p theory_col and produce an explanation for the SAT solver.
    * @param theory_col theory column the bounds are associated with
@@ -322,7 +321,7 @@ class TheorySolver {
    * @param type type of violation the bound is associated with
    * @param[out] bound_idxs set of indexes of the bounds
    */
-  void TheoryBoundsToBoundIdxs(int theory_col, BoundViolationType type, std::set<int> &bound_idxs) const;
+  void TheoryBoundsToBoundIdxs(int theory_col, BoundViolationType type, std::set<int> &bound_idxs);
   /**
    * Generate a tuple (var, type, value) that represents a bound on the variable.
    *

dlinear/solver/TheorySolverBoundPreprocessor.cpp

@@ -114,6 +114,14 @@ bool TheorySolverBoundPreprocessor::AddConstraint(const int theory_row, const Fo
   return true;
 }
 
+void TheorySolverBoundPreprocessor::GetActiveExplanation(int theory_col, LiteralSet& explanation) {
+  GetExplanation(theory_cols_[theory_col], explanation);
+}
+
+void TheorySolverBoundPreprocessor::GetActiveBoundIdxs(const int theory_col, std::set<int>& bound_idxs) {
+  GetExplanationBoundIdxs(theory_cols_[theory_col], bound_idxs);
+}
+
 TheorySolverBoundPreprocessor::Explanations TheorySolverBoundPreprocessor::Process(
     const std::vector<int>& enabled_theory_rows) {
   Explanations explanations;
@@ -162,6 +170,7 @@ void TheorySolverBoundPreprocessor::SetEnvironmentFromBounds() {
     if (active_bound == nullptr) continue;
     const Variable& var = theory_cols_[theory_col];
     env_[var] = *active_bound;
+    DLINEAR_TRACE_FMT("TheorySolverBoundPreprocessor::SetEnvironmentFromBounds: {} = {}", var, *active_bound);
   }
 }
 
@@ -170,18 +179,29 @@ TheorySolverBoundPreprocessor::PropagateEqBinomialResult TheorySolverBoundPrepro
   DLINEAR_TRACE_FMT("TheorySolverBoundPreprocessor::PropagateEqBinomial({})", theory_row);
 
   const mpq_class& coeff = row_to_eq_binomial_edge_coefficients_.at(theory_row);
-  const auto& [from, to] = ExtractBoundEdge(predicate_abstractor_[theory_rows_[theory_row].var]);
+  const Formula& formula = predicate_abstractor_[theory_rows_[theory_row].var];
+  const auto& [from, to] = ExtractBoundEdge(formula);
+
+  DLINEAR_ASSERT(is_equal_to(formula), "Formula should be an equality relation");
+  DLINEAR_ASSERT(formula.GetFreeVariables().size() == 2, "Formula should have exactly two free variables");
+  DLINEAR_ASSERT(formula.IsFlattened(), "The formula must be flattened");
+
+  const Expression& lhs = get_lhs_expression(formula);
+  const std::map<Expression, mpq_class>& map = get_expr_to_coeff_map_in_addition(lhs);
+  DLINEAR_ASSERT(map.size() == 2, "Expression should have exactly two variables");
+  DLINEAR_ASSERT(get_constant_value(get_rhs_expression(formula)) == 0, "The right-hand side must be 0");
+  const bool from_is_first = get_variable(map.cbegin()->first).equal_to(from);
 
   Environment::const_iterator updater_it;
   Variable to_update_variable;
   mpq_class new_value;
   if ((updater_it = env_.find(from)) != env_.end()) {
     DLINEAR_TRACE_FMT("TheorySolverBoundPreprocessor::PropagateEqBinomial: {} --{} --> {}", to, theory_row, from);
-    new_value = updater_it->second * coeff;
+    new_value = from_is_first ? mpq_class{updater_it->second / coeff} : mpq_class{updater_it->second * coeff};
     to_update_variable = to;
   } else if ((updater_it = env_.find(to)) != env_.end()) {
     DLINEAR_TRACE_FMT("TheorySolverBoundPreprocessor::PropagateEqBinomial: {} --{} --> {}", from, theory_row, to);
-    new_value = updater_it->second / coeff;
+    new_value = from_is_first ? mpq_class{updater_it->second * coeff} : mpq_class{updater_it->second} / coeff;
     to_update_variable = from;
   } else {
     // Neither of the two variables is in the environment. Can't propagate
@@ -332,7 +352,7 @@ bool TheorySolverBoundPreprocessor::ShouldEvaluate(const Formula& formula) const
   DLINEAR_TRACE_FMT("TheorySolverBoundPreprocessor::ShouldEvaluate({})", formula);
   // All free variables must be in the environment
   return std::all_of(formula.GetFreeVariables().begin(), formula.GetFreeVariables().end(),
-                     [this](const Variable& v) { return env_.find(v) != env_.end(); });
+                     [this](const Variable& v) { return env_.contains(v); });
 }
 
 bool TheorySolverBoundPreprocessor::ShouldPropagateEqBinomial(const Literal& lit) const {
@@ -487,6 +507,47 @@ void TheorySolverBoundPreprocessor::GetExplanation(const Variable& var, LiteralS
   }
 }
 
+void TheorySolverBoundPreprocessor::AddPathToExplanationBoundIdxs(const Variable& from, const Variable& to,
+                                                                  const TheorySolverBoundVector& from_bounds,
+                                                                  const TheorySolverBoundVector& to_bounds,
+                                                                  std::set<int>& explanation) {
+  DLINEAR_ASSERT_FMT(to_bounds.GetActiveEqualityBound() != nullptr, "The ending variable {} must have an eq bound", to);
+  graph_.AllPaths(from, to, [&](std::vector<Variable>& path) {
+    // Insert start and end bounds to the explanation
+    if (from_bounds.GetActiveEqualityBound() != nullptr) {
+      for (auto it = from_bounds.GetActiveBound(); it; ++it) explanation.insert(it->idx);
+    }
+    if (to_bounds.GetActiveEqualityBound() != nullptr) {
+      for (auto it = to_bounds.GetActiveBound(); it; ++it) explanation.insert(it->idx);
+    }
+    // Insert all rows from the edges in the path to the explanation
+    for (std::size_t i = 1; i < path.size(); i++) {
+      DLINEAR_ASSERT(graph_.GetEdgeWeight(path[i - 1], path[i]) != nullptr, "Edge must exist");
+      const int theory_row = *graph_.GetEdgeWeight(path[i - 1], path[i]);
+      explanation.insert(theory_row);
+    }
+    return VisitResult::STOP;
+  });
+}
+
+void TheorySolverBoundPreprocessor::GetExplanationBoundIdxs(const dlinear::drake::symbolic::Variable& var,
+                                                            std::set<int>& bound_idxs) {
+  const TheorySolverBoundVector& bounds = theory_bounds_.at(var_to_cols_.at(var.get_id()));
+  // If the variable has its bounds set directly by some literals, simply add them to the explanation
+  if (bounds.GetActiveEqualityBound() != nullptr) {
+    for (auto it = bounds.GetActiveBound(); it; ++it) bound_idxs.insert(it->idx);
+  } else {  // else we need to find the variable responsible for the bound propagation from the bound_graph
+    graph_.BFS(var, [&](const Variable& from, const Variable& to, const Weight&) {
+      if (to.equal_to(from)) return VisitResult::CONTINUE;
+      const TheorySolverBoundVector& to_bounds = theory_bounds_.at(var_to_cols_.at(to.get_id()));
+      if (to_bounds.GetActiveEqualityBound() == nullptr) return VisitResult::CONTINUE;
+      AddPathToExplanationBoundIdxs(var, to, bounds, to_bounds, bound_idxs);
+      return VisitResult::SKIP;
+    });
+  }
+}
+
+#if DEBUGGING_PREPROCESSOR
 std::vector<Variable> TheorySolverBoundPreprocessor::GetExplanationOrigins(const Variable& var) {
   const TheorySolverBoundVector& bounds = theory_bounds_.at(var_to_cols_.at(var.get_id()));
   // If the variable has its bounds set directly by some literals, it is its own origin
@@ -505,14 +566,14 @@ std::vector<Variable> TheorySolverBoundPreprocessor::GetExplanationOrigins(const
   DLINEAR_ASSERT(!origins.empty(), "There must be at least one origin");
   return origins;
 }
+#endif
 
 std::ostream& operator<<(std::ostream& os, const TheorySolverBoundPreprocessor& preprocessor) {
-  return os << "TheorySolverBoundPreprocessor{"
-            << "env_ = " << preprocessor.env() << ", "
+  return os << "TheorySolverBoundPreprocessor{" << "env_ = " << preprocessor.env() << ", "
             << "theory_cols_ = " << preprocessor.theory_cols() << ", "
             << "theory_rows_ = " << preprocessor.theory_rows() << ", "
-            << "theory_bounds_ = " << preprocessor.theory_bounds() << ", "
-            << "graph_ = " << preprocessor.bound_graph() << "}";
+            << "theory_bounds_ = " << preprocessor.theory_bounds() << ", " << "graph_ = " << preprocessor.bound_graph()
+            << "}";
 }
 
 }  // namespace dlinear

dlinear/solver/TheorySolverBoundPreprocessor.h

@@ -56,6 +56,9 @@ class TheorySolverBoundPreprocessor {
   Explanations Process(const std::vector<int>& enabled_theory_rows = {});
   void Process(const std::vector<int>& enabled_theory_rows, Explanations& explanations);
 
+  void GetActiveExplanation(int theory_col, LiteralSet& explanation);
+  void GetActiveBoundIdxs(int theory_col, std::set<int>& bound_idxs);
+
   void Clear();
 
   [[nodiscard]] const Config& config() const { return config_; }
@@ -104,6 +107,9 @@ class TheorySolverBoundPreprocessor {
   void AddPathToExplanation(const Variable& from, const Variable& to, LiteralSet& explanation);
   void AddPathToExplanation(const Variable& from, const Variable& to, const TheorySolverBoundVector& from_bounds,
                             const TheorySolverBoundVector& to_bounds, LiteralSet& explanation);
+  void AddPathToExplanationBoundIdxs(const Variable& from, const Variable& to,
+                                     const TheorySolverBoundVector& from_bounds,
+                                     const TheorySolverBoundVector& to_bounds, std::set<int>& explanation);
 
   std::pair<Variable, Variable> ExtractBoundEdge(const Formula& formula) const;
   /**
@@ -118,7 +124,11 @@ class TheorySolverBoundPreprocessor {
   mpq_class ExtractEqBoundCoefficient(const Formula& formula) const;
 
   void GetExplanation(const Variable& var, LiteralSet& explanation);
+  void GetExplanationBoundIdxs(const Variable& var, std::set<int>& bound_idxs);
+
+#if DEBUGGING_PREPROCESSOR
   std::vector<Variable> GetExplanationOrigins(const Variable& var);
+#endif
 
  private:
   const Config& config_;

test/solver/smt2/bignum_lra1.smt2

@@ -10,6 +10,14 @@
 (declare-fun x4 () Real)
 (declare-fun x5 () Real)
 (declare-fun x6 () Real)
-(assert (and (or (= x1 (/ 1 100000)) (= x1 (/ 1 100002))) (or (= x2 (* (/ 1 120030) x1)) (= x2 (* (/ 1 10003) x1))) (or (= x3 (* (/ 1 100310) x2)) (= x3 (* (/ 1 199900) x2))) (or (= x4 (* (/ 1 400000) x3)) (= x4 (* (/ 1 800000) x3))) (or (= x5 (* (/ (- 1) 40000) x4)) (= x5 (* (/ (- 1) 8000) x4))) (or (= x6 (* (/ (- 1) 3000) x5)) (= x6 (* (/ (- 1) 2000) x5))) (> x6 0)))
+(assert (and 
+    (or (= x1 (/ 1 100000)) (= x1 (/ 1 100002))) 
+    (or (= x2 (* (/ 1 120030) x1)) (= x2 (* (/ 1 10003) x1))) 
+    (or (= x3 (* (/ 1 100310) x2)) (= x3 (* (/ 1 199900) x2))) 
+    (or (= x4 (* (/ 1 400000) x3)) (= x4 (* (/ 1 800000) x3))) 
+    (or (= x5 (* (/ (- 1) 40000) x4)) (= x5 (* (/ (- 1) 8000) x4))) 
+    (or (= x6 (* (/ (- 1) 3000) x5)) (= x6 (* (/ (- 1) 2000) x5))) 
+    (> x6 0)
+))
 (check-sat)
 (exit)

test/solver/smt2/preprocessor.smt2

@@ -0,0 +1,19 @@
+(set-info :smt-lib-version 2.6)
+(set-logic QF_LRA)
+(set-info :category "industrial")
+(set-info :status unsat)
+(declare-fun x1 () Real)
+(declare-fun x2 () Real)
+(declare-fun x3 () Real)
+(declare-fun x4 () Real)
+(declare-fun u () Real)
+(declare-fun g () Real)
+(assert (= x1 1))
+(assert (<= g 0))
+(assert (= x2 (+ 1 (* x1 2))))
+(assert (= x3 (+ 2 (* x2 3))))
+(assert (= x4 (+ 7 (* x3 5))))
+(assert (= u (+ 2 (* x1 5))))
+(assert (< x4 g))
+(check-sat)
+(exit)