Refactor transcendental solver (#5539)

This PR does another round of refactoring on the transcendental solver. It cleans up some variable names, introduces an enum type for Convexity and passes both the intended taylor degree and the actual taylor degree (which will be needed for proofs).

9 files changed, 176 insertions, 93 deletions

diff --git a/src/theory/arith/nl/transcendental/exponential_solver.cpp b/src/theory/arith/nl/transcendental/exponential_solver.cpp
index 7d7d0c23c..482e3bc21 100644
--- a/src/theory/arith/nl/transcendental/exponential_solver.cpp
+++ b/src/theory/arith/nl/transcendental/exponential_solver.cpp
@@ -170,31 +170,41 @@ void ExponentialSolver::doTangentLemma(TNode e, TNode c, TNode poly_approx)
   d_data->d_im.addPendingArithLemma(lem, InferenceId::NL_T_TANGENT, nullptr, true);
 }
 
-void ExponentialSolver::doSecantLemmas(
-    TNode e, TNode poly_approx, TNode c, TNode poly_approx_c, unsigned d)
+void ExponentialSolver::doSecantLemmas(TNode e,
+                                       TNode poly_approx,
+                                       TNode center,
+                                       TNode cval,
+                                       unsigned d,
+                                       unsigned actual_d)
 {
-  d_data->doSecantLemmas(
-      getSecantBounds(e, c, d), poly_approx, c, poly_approx_c, e, d, 1);
+  d_data->doSecantLemmas(getSecantBounds(e, center, d),
+                         poly_approx,
+                         center,
+                         cval,
+                         e,
+                         Convexity::CONVEX,
+                         d,
+                         actual_d);
 }
 
 std::pair<Node, Node> ExponentialSolver::getSecantBounds(TNode e,
-                                                         TNode c,
+                                                         TNode center,
                                                          unsigned d)
 {
-  std::pair<Node, Node> bounds = d_data->getClosestSecantPoints(e, c, d);
+  std::pair<Node, Node> bounds = d_data->getClosestSecantPoints(e, center, d);
 
   // Check if we already have neighboring secant points
   if (bounds.first.isNull())
   {
     // pick c-1
     bounds.first = Rewriter::rewrite(
-        NodeManager::currentNM()->mkNode(Kind::MINUS, c, d_data->d_one));
+        NodeManager::currentNM()->mkNode(Kind::MINUS, center, d_data->d_one));
   }
   if (bounds.second.isNull())
   {
     // pick c+1
     bounds.second = Rewriter::rewrite(
-        NodeManager::currentNM()->mkNode(Kind::PLUS, c, d_data->d_one));
+        NodeManager::currentNM()->mkNode(Kind::PLUS, center, d_data->d_one));
   }
   return bounds;
 }
diff --git a/src/theory/arith/nl/transcendental/exponential_solver.h b/src/theory/arith/nl/transcendental/exponential_solver.h
index b20c23e56..b4d08559a 100644
--- a/src/theory/arith/nl/transcendental/exponential_solver.h
+++ b/src/theory/arith/nl/transcendental/exponential_solver.h
@@ -87,12 +87,16 @@ class ExponentialSolver
   void doTangentLemma(TNode e, TNode c, TNode poly_approx);
 
   /** Sent secant lemmas around c for e */
-  void doSecantLemmas(
-      TNode e, TNode poly_approx, TNode c, TNode poly_approx_c, unsigned d);
+  void doSecantLemmas(TNode e,
+                      TNode poly_approx,
+                      TNode center,
+                      TNode cval,
+                      unsigned d,
+                      unsigned actual_d);
 
  private:
   /** Generate bounds for secant lemmas */
-  std::pair<Node, Node> getSecantBounds(TNode e, TNode c, unsigned d);
+  std::pair<Node, Node> getSecantBounds(TNode e, TNode center, unsigned d);
 
   /** Holds common state for transcendental solvers */
   TranscendentalState* d_data;
diff --git a/src/theory/arith/nl/transcendental/sine_solver.cpp b/src/theory/arith/nl/transcendental/sine_solver.cpp
index 31fd47503..cba85b80e 100644
--- a/src/theory/arith/nl/transcendental/sine_solver.cpp
+++ b/src/theory/arith/nl/transcendental/sine_solver.cpp
@@ -266,19 +266,20 @@ void SineSolver::doTangentLemma(TNode e, TNode c, TNode poly_approx, int region)
   // Figure 3: T( x )
   // We use zero slope tangent planes, since the concavity of the Taylor
   // approximation cannot be easily established.
-  int concavity = regionToConcavity(region);
+  Convexity convexity = regionToConvexity(region);
   int mdir = regionToMonotonicityDir(region);
+  bool usec = (mdir == 1) == (convexity == Convexity::CONCAVE);
   Node lem = nm->mkNode(
       Kind::IMPLIES,
       nm->mkNode(
           Kind::AND,
-          nm->mkNode(Kind::GEQ,
-                     e[0],
-                     mdir == concavity ? Node(c) : regionToLowerBound(region)),
-          nm->mkNode(Kind::LEQ,
-                     e[0],
-                     mdir != concavity ? Node(c) : regionToUpperBound(region))),
-      nm->mkNode(concavity == 1 ? Kind::GEQ : Kind::LEQ, e, poly_approx));
+          nm->mkNode(
+              Kind::GEQ, e[0], usec ? Node(c) : regionToLowerBound(region)),
+          nm->mkNode(
+              Kind::LEQ, e[0], usec ? Node(c) : regionToUpperBound(region))),
+      nm->mkNode(convexity == Convexity::CONVEX ? Kind::GEQ : Kind::LEQ,
+                 e,
+                 poly_approx));
 
   Trace("nl-ext-sine") << "*** Tangent plane lemma (pre-rewrite): " << lem
                        << std::endl;
@@ -294,6 +295,7 @@ void SineSolver::doSecantLemmas(TNode e,
                                 TNode c,
                                 TNode poly_approx_c,
                                 unsigned d,
+                                unsigned actual_d,
                                 int region)
 {
   d_data->doSecantLemmas(getSecantBounds(e, c, d, region),
@@ -301,8 +303,9 @@ void SineSolver::doSecantLemmas(TNode e,
                          c,
                          poly_approx_c,
                          e,
+                         regionToConvexity(region),
                          d,
-                         regionToConcavity(region));
+                         actual_d);
 }
 
 std::pair<Node, Node> SineSolver::getSecantBounds(TNode e,
diff --git a/src/theory/arith/nl/transcendental/sine_solver.h b/src/theory/arith/nl/transcendental/sine_solver.h
index 9f9102a53..e41e6bd4f 100644
--- a/src/theory/arith/nl/transcendental/sine_solver.h
+++ b/src/theory/arith/nl/transcendental/sine_solver.h
@@ -93,6 +93,7 @@ class SineSolver
                       TNode c,
                       TNode poly_approx_c,
                       unsigned d,
+                      unsigned actual_d,
                       int region);
 
  private:
@@ -152,15 +153,15 @@ class SineSolver
       default: return 0;
     }
   }
-  int regionToConcavity(int region)
+  Convexity regionToConvexity(int region)
   {
     switch (region)
     {
       case 1:
-      case 2: return -1;
+      case 2: return Convexity::CONCAVE;
       case 3:
-      case 4: return 1;
-      default: return 0;
+      case 4: return Convexity::CONVEX;
+      default: return Convexity::UNKNOWN;
     }
   }
 
diff --git a/src/theory/arith/nl/transcendental/taylor_generator.cpp b/src/theory/arith/nl/transcendental/taylor_generator.cpp
index a373339e9..1b7257f8f 100644
--- a/src/theory/arith/nl/transcendental/taylor_generator.cpp
+++ b/src/theory/arith/nl/transcendental/taylor_generator.cpp
@@ -212,7 +212,7 @@ void TaylorGenerator::getPolynomialApproximationBounds(
   }
 }
 
-void TaylorGenerator::getPolynomialApproximationBoundForArg(
+unsigned TaylorGenerator::getPolynomialApproximationBoundForArg(
     Kind k, Node c, unsigned d, std::vector<Node>& pbounds)
 {
   getPolynomialApproximationBounds(k, d, pbounds);
@@ -252,7 +252,9 @@ void TaylorGenerator::getPolynomialApproximationBoundForArg(
       getPolynomialApproximationBounds(k, ds, pboundss);
       pbounds[2] = pboundss[2];
     }
+    return ds;
   }
+  return d;
 }
 
 std::pair<Node, Node> TaylorGenerator::getTfModelBounds(Node tf, unsigned d, NlModel& model)
diff --git a/src/theory/arith/nl/transcendental/taylor_generator.h b/src/theory/arith/nl/transcendental/taylor_generator.h
index c647f7fd2..2dbfcccde 100644
--- a/src/theory/arith/nl/transcendental/taylor_generator.h
+++ b/src/theory/arith/nl/transcendental/taylor_generator.h
@@ -79,11 +79,13 @@ class TaylorGenerator
    * polynomials may depend on c. In particular, for P_u+[x] for <k>( c ) where
    * c>0, we return the P_u+[x] from the function above for the minimum degree
    * d' >= d such that (1-c^{2*d'+1}/(2*d'+1)!) is positive.
+   * @return the actual degree of the polynomial approximations (which may be
+   * larger than d).
    */
-  void getPolynomialApproximationBoundForArg(Kind k,
-                                             Node c,
-                                             unsigned d,
-                                             std::vector<Node>& pbounds);
+  unsigned getPolynomialApproximationBoundForArg(Kind k,
+                                                 Node c,
+                                                 unsigned d,
+                                                 std::vector<Node>& pbounds);
 
   /** get transcendental function model bounds
    *
diff --git a/src/theory/arith/nl/transcendental/transcendental_solver.cpp b/src/theory/arith/nl/transcendental/transcendental_solver.cpp
index 8eb69e50b..ad3f49576 100644
--- a/src/theory/arith/nl/transcendental/transcendental_solver.cpp
+++ b/src/theory/arith/nl/transcendental/transcendental_solver.cpp
@@ -233,7 +233,8 @@ bool TranscendentalSolver::checkTfTangentPlanesFun(Node tf, unsigned d)
   // mapped to for signs of c
   std::map<int, Node> poly_approx_bounds[2];
   std::vector<Node> pbounds;
-  d_tstate.d_taylor.getPolynomialApproximationBoundForArg(k, c, d, pbounds);
+  unsigned actual_d =
+      d_tstate.d_taylor.getPolynomialApproximationBoundForArg(k, c, d, pbounds);
   poly_approx_bounds[0][1] = pbounds[0];
   poly_approx_bounds[0][-1] = pbounds[1];
   poly_approx_bounds[1][1] = pbounds[2];
@@ -362,11 +363,12 @@ bool TranscendentalSolver::checkTfTangentPlanesFun(Node tf, unsigned d)
   {
     if (k == EXPONENTIAL)
     {
-      d_expSlv.doSecantLemmas(tf, poly_approx, c, poly_approx_c, d);
+      d_expSlv.doSecantLemmas(tf, poly_approx, c, poly_approx_c, d, actual_d);
     }
     else if (k == Kind::SINE)
     {
-      d_sineSlv.doSecantLemmas(tf, poly_approx, c, poly_approx_c, d, region);
+      d_sineSlv.doSecantLemmas(
+          tf, poly_approx, c, poly_approx_c, d, actual_d, region);
     }
   }
   return true;
diff --git a/src/theory/arith/nl/transcendental/transcendental_state.cpp b/src/theory/arith/nl/transcendental/transcendental_state.cpp
index 41ed2c53d..69678c621 100644
--- a/src/theory/arith/nl/transcendental/transcendental_state.cpp
+++ b/src/theory/arith/nl/transcendental/transcendental_state.cpp
@@ -267,26 +267,26 @@ void TranscendentalState::getCurrentPiBounds()
 }
 
 std::pair<Node, Node> TranscendentalState::getClosestSecantPoints(TNode e,
-                                                                  TNode c,
+                                                                  TNode center,
                                                                   unsigned d)
 {
   // bounds are the minimum and maximum previous secant points
   // should not repeat secant points: secant lemmas should suffice to
   // rule out previous assignment
-  Assert(
-      std::find(d_secant_points[e][d].begin(), d_secant_points[e][d].end(), c)
-      == d_secant_points[e][d].end());
+  Assert(std::find(
+             d_secant_points[e][d].begin(), d_secant_points[e][d].end(), center)
+         == d_secant_points[e][d].end());
   // Insert into the (temporary) vector. We do not update this vector
   // until we are sure this secant plane lemma has been processed. We do
   // this by mapping the lemma to a side effect below.
   std::vector<Node> spoints = d_secant_points[e][d];
-  spoints.push_back(c);
+  spoints.push_back(center);
 
   // sort
   sortByNlModel(spoints.begin(), spoints.end(), &d_model);
   // get the resulting index of c
   unsigned index =
-      std::find(spoints.begin(), spoints.end(), c) - spoints.begin();
+      std::find(spoints.begin(), spoints.end(), center) - spoints.begin();
 
   // bounds are the next closest upper/lower bound values
   return {index > 0 ? spoints[index - 1] : Node(),
@@ -294,24 +294,40 @@ std::pair<Node, Node> TranscendentalState::getClosestSecantPoints(TNode e,
 }
 
 Node TranscendentalState::mkSecantPlane(
-    TNode arg, TNode b, TNode c, TNode approx_b, TNode approx_c)
+    TNode arg, TNode lower, TNode upper, TNode lval, TNode uval)
 {
   NodeManager* nm = NodeManager::currentNM();
   // Figure 3: S_l( x ), S_u( x ) for s = 0,1
-  Node rcoeff_n = Rewriter::rewrite(nm->mkNode(Kind::MINUS, b, c));
+  Node rcoeff_n = Rewriter::rewrite(nm->mkNode(Kind::MINUS, lower, upper));
   Assert(rcoeff_n.isConst());
   Rational rcoeff = rcoeff_n.getConst<Rational>();
   Assert(rcoeff.sgn() != 0);
-  return nm->mkNode(Kind::PLUS,
-                    approx_b,
-                    nm->mkNode(Kind::MULT,
-                               nm->mkNode(Kind::MINUS, approx_b, approx_c),
-                               nm->mkConst(rcoeff.inverse()),
-                               nm->mkNode(Kind::MINUS, arg, b)));
+  Node res =
+      nm->mkNode(Kind::PLUS,
+                 lval,
+                 nm->mkNode(Kind::MULT,
+                            nm->mkNode(Kind::DIVISION,
+                                       nm->mkNode(Kind::MINUS, lval, uval),
+                                       nm->mkNode(Kind::MINUS, lower, upper)),
+                            nm->mkNode(Kind::MINUS, arg, lower)));
+  Trace("nl-trans") << "Creating secant plane for transcendental function of "
+                    << arg << std::endl;
+  Trace("nl-trans") << "\tfrom ( " << lower << " ; " << lval << " ) to ( "
+                    << upper << " ; " << uval << " )" << std::endl;
+  Trace("nl-trans") << "\t" << res << std::endl;
+  Trace("nl-trans") << "\trewritten: " << Rewriter::rewrite(res) << std::endl;
+  return res;
 }
 
-NlLemma TranscendentalState::mkSecantLemma(
-    TNode lower, TNode upper, int concavity, TNode tf, TNode splane)
+NlLemma TranscendentalState::mkSecantLemma(TNode lower,
+                                           TNode upper,
+                                           TNode lapprox,
+                                           TNode uapprox,
+                                           int csign,
+                                           Convexity convexity,
+                                           TNode tf,
+                                           TNode splane,
+                                           unsigned actual_d)
 {
   NodeManager* nm = NodeManager::currentNM();
   // With respect to Figure 3, this is slightly different.
@@ -329,60 +345,73 @@ NlLemma TranscendentalState::mkSecantLemma(
   Node antec_n = nm->mkNode(Kind::AND,
                             nm->mkNode(Kind::GEQ, tf[0], lower),
                             nm->mkNode(Kind::LEQ, tf[0], upper));
+  Trace("nl-trans") << "Bound for secant plane: " << lower << " <= " << tf[0]
+                    << " <= " << upper << std::endl;
+  Trace("nl-trans") << "\t" << antec_n << std::endl;
+  // Convex: actual value is below the secant
+  // Concave: actual value is above the secant
   Node lem = nm->mkNode(
       Kind::IMPLIES,
       antec_n,
-      nm->mkNode(concavity == 1 ? Kind::LEQ : Kind::GEQ, tf, splane));
-  Trace("nl-ext-tftp-debug2")
-      << "*** Secant plane lemma (pre-rewrite) : " << lem << std::endl;
+      nm->mkNode(
+          convexity == Convexity::CONVEX ? Kind::LEQ : Kind::GEQ, tf, splane));
+  Trace("nl-trans-lemma") << "*** Secant plane lemma (pre-rewrite) : " << lem
+                          << std::endl;
   lem = Rewriter::rewrite(lem);
-  Trace("nl-ext-tftp-lemma") << "*** Secant plane lemma : " << lem << std::endl;
+  Trace("nl-trans-lemma") << "*** Secant plane lemma : " << lem << std::endl;
   Assert(d_model.computeAbstractModelValue(lem) == d_false);
   return NlLemma(lem, LemmaProperty::NONE, nullptr, InferenceId::NL_T_SECANT);
 }
 
 void TranscendentalState::doSecantLemmas(const std::pair<Node, Node>& bounds,
                                          TNode poly_approx,
-                                         TNode c,
-                                         TNode approx_c,
+                                         TNode center,
+                                         TNode cval,
                                          TNode tf,
+                                         Convexity convexity,
                                          unsigned d,
-                                         int concavity)
+                                         unsigned actual_d)
 {
-  Trace("nl-ext-tftp-debug2") << "...secant bounds are : " << bounds.first
-                              << " ... " << bounds.second << std::endl;
-  // take the model value of l or u (since may contain PI)
-  // Make secant from bounds.first to c
-  Node lval = d_model.computeAbstractModelValue(bounds.first);
-  Trace("nl-ext-tftp-debug2") << "...model value of bound " << bounds.first
-                              << " is " << lval << std::endl;
-  if (lval != c)
+  int csign = center.getConst<Rational>().sgn();
+  Trace("nl-trans") << "...do secant lemma with center " << center << " val "
+                    << cval << " sign " << csign << std::endl;
+  Trace("nl-trans") << "...secant bounds are : " << bounds.first << " ... "
+                    << bounds.second << std::endl;
+  // take the model value of lower (since may contain PI)
+  // Make secant from bounds.first to center
+  Node lower = d_model.computeAbstractModelValue(bounds.first);
+  Trace("nl-trans") << "...model value of bound " << bounds.first << " is "
+                    << lower << std::endl;
+  if (lower != center)
   {
     // Figure 3 : P(l), P(u), for s = 0
-    Node approx_l = Rewriter::rewrite(
-        poly_approx.substitute(d_taylor.getTaylorVariable(), lval));
-    Node splane = mkSecantPlane(tf[0], lval, c, approx_l, approx_c);
-    NlLemma nlem = mkSecantLemma(lval, c, concavity, tf, splane);
+    Node lval = Rewriter::rewrite(
+        poly_approx.substitute(d_taylor.getTaylorVariable(), lower));
+    Node splane = mkSecantPlane(tf[0], lower, center, lval, cval);
+    NlLemma nlem = mkSecantLemma(
+        lower, center, lval, cval, csign, convexity, tf, splane, actual_d);
     // The side effect says that if lem is added, then we should add the
     // secant point c for (tf,d).
-    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, c));
+    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, center));
     d_im.addPendingArithLemma(nlem, true);
   }
 
-  // Make secant from c to bounds.second
-  Node uval = d_model.computeAbstractModelValue(bounds.second);
-  Trace("nl-ext-tftp-debug2") << "...model value of bound " << bounds.second
-                              << " is " << uval << std::endl;
-  if (c != uval)
+  // take the model value of upper (since may contain PI)
+  // Make secant from center to bounds.second
+  Node upper = d_model.computeAbstractModelValue(bounds.second);
+  Trace("nl-trans") << "...model value of bound " << bounds.second << " is "
+                    << upper << std::endl;
+  if (center != upper)
   {
     // Figure 3 : P(l), P(u), for s = 1
-    Node approx_u = Rewriter::rewrite(
-        poly_approx.substitute(d_taylor.getTaylorVariable(), uval));
-    Node splane = mkSecantPlane(tf[0], c, uval, approx_c, approx_u);
-    NlLemma nlem = mkSecantLemma(c, uval, concavity, tf, splane);
+    Node uval = Rewriter::rewrite(
+        poly_approx.substitute(d_taylor.getTaylorVariable(), upper));
+    Node splane = mkSecantPlane(tf[0], center, upper, cval, uval);
+    NlLemma nlem = mkSecantLemma(
+        center, upper, cval, uval, csign, convexity, tf, splane, actual_d);
     // The side effect says that if lem is added, then we should add the
     // secant point c for (tf,d).
-    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, c));
+    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, center));
     d_im.addPendingArithLemma(nlem, true);
   }
 }
diff --git a/src/theory/arith/nl/transcendental/transcendental_state.h b/src/theory/arith/nl/transcendental/transcendental_state.h
index 0a4e46eca..e1510c3b3 100644
--- a/src/theory/arith/nl/transcendental/transcendental_state.h
+++ b/src/theory/arith/nl/transcendental/transcendental_state.h
@@ -27,6 +27,24 @@ namespace nl {
 namespace transcendental {
 
 /**
+ * This enum indicates whether some function is convex, concave or unknown at
+ * some point.
+ */
+enum class Convexity
+{
+  CONVEX,
+  CONCAVE,
+  UNKNOWN
+};
+inline std::ostream& operator<<(std::ostream& os, Convexity c) {
+  switch (c) {
+    case Convexity::CONVEX: return os << "CONVEX";
+    case Convexity::CONCAVE: return os << "CONCAVE";
+    default: return os << "UNKNOWN";
+  }
+}
+
+/**
  * Holds common state and utilities for transcendental solvers.
  *
  * This includes common lookups and caches as well as generic utilities for
@@ -60,20 +78,24 @@ struct TranscendentalState
    * Get the two closest secant points from the once stored already.
    * "closest" is determined according to the current model.
    * @param e The transcendental term (like (exp t))
-   * @param c The point currently under consideration (probably the model of t)
+   * @param center The point currently under consideration (probably the model
+   * of t)
    * @param d The taylor degree.
    */
-  std::pair<Node, Node> getClosestSecantPoints(TNode e, TNode c, unsigned d);
+  std::pair<Node, Node> getClosestSecantPoints(TNode e,
+                                               TNode center,
+                                               unsigned d);
 
   /**
-   * Construct a secant plane between b and c
+   * Construct a secant plane as function in arg between lower and upper
    * @param arg The argument of the transcendental term
-   * @param b Left secant point
-   * @param c Right secant point
-   * @param approx Approximation for b (not yet substituted)
-   * @param approx_c Approximation for c (already substituted)
+   * @param lower Left secant point
+   * @param upper Right secant point
+   * @param lval Evaluation at lower
+   * @param uval Evaluation at upper
    */
-  Node mkSecantPlane(TNode arg, TNode b, TNode c, TNode approx, TNode approx_c);
+  Node mkSecantPlane(
+      TNode arg, TNode lower, TNode upper, TNode lval, TNode uval);
 
   /**
    * Construct a secant lemma between lower and upper for tf.
@@ -83,26 +105,34 @@ struct TranscendentalState
    * @param tf Current transcendental term
    * @param splane Secant plane as computed by mkSecantPlane()
    */
-  NlLemma mkSecantLemma(
-      TNode lower, TNode upper, int concavity, TNode tf, TNode splane);
+  NlLemma mkSecantLemma(TNode lower,
+                        TNode upper,
+                        TNode lapprox,
+                        TNode uapprox,
+                        int csign,
+                        Convexity convexity,
+                        TNode tf,
+                        TNode splane,
+                        unsigned actual_d);
 
   /**
    * Construct and send secant lemmas (if appropriate)
    * @param bounds Secant bounds
    * @param poly_approx Polynomial approximation
-   * @param c Current point
-   * @param approx_c Approximation for c
+   * @param center Current point
+   * @param cval Evaluation at c
    * @param tf Current transcendental term
    * @param d Current taylor degree
    * @param concavity Concavity in region of c
    */
   void doSecantLemmas(const std::pair<Node, Node>& bounds,
                       TNode poly_approx,
-                      TNode c,
-                      TNode approx_c,
+                      TNode center,
+                      TNode cval,
                       TNode tf,
+                      Convexity convexity,
                       unsigned d,
-                      int concavity);
+                      unsigned actual_d);
 
   Node d_true;
   Node d_false;