Use side effect utility for non-linear lemmas (#3780)

Fixes #3647.

Previously, when doing incremental linearization for transcedental functions, we would add points to the list of secant points at the time when linearization lemmas were generated. However, our strategy has been updated such that lemmas may be abandoned (say in the case that a higher priority lemma is found). Thus, our list of secant points had spurious points corresponding to lemmas that weren't sent. This led to assertion failures, and likely led to gaps in our linearization, hindering our ability to say sat/unsat.

This PR introduces a "lemma side effect" class to ensure that modifications to the state of the nonlinear solver are in sync with the lemmas we send out.

4 files changed, 205 insertions, 89 deletions

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 30cdf15b7..52f46147c 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -298,6 +298,7 @@ libcvc4_add_sources(
   theory/arith/linear_equality.h
   theory/arith/matrix.cpp
   theory/arith/matrix.h
+  theory/arith/nl_lemma_utils.h
   theory/arith/nl_model.cpp
   theory/arith/nl_model.h
   theory/arith/nonlinear_extension.cpp
diff --git a/src/theory/arith/nl_lemma_utils.h b/src/theory/arith/nl_lemma_utils.h
new file mode 100644
index 000000000..74886a1fb
--- /dev/null
+++ b/src/theory/arith/nl_lemma_utils.h
@@ -0,0 +1,53 @@
+/*********************                                                        */
+/*! \file nl_lemma_utils.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Andrew Reynolds
+ ** This file is part of the CVC4 project.
+ ** Copyright (c) 2009-2019 by the authors listed in the file AUTHORS
+ ** in the top-level source directory) and their institutional affiliations.
+ ** All rights reserved.  See the file COPYING in the top-level source
+ ** directory for licensing information.\endverbatim
+ **
+ ** \brief Utilities for processing lemmas from the non-linear solver
+ **/
+
+#ifndef CVC4__THEORY__ARITH__NL_LEMMA_UTILS_H
+#define CVC4__THEORY__ARITH__NL_LEMMA_UTILS_H
+
+#include <tuple>
+#include <vector>
+#include "expr/node.h"
+
+namespace CVC4 {
+namespace theory {
+namespace arith {
+
+/**
+ * A side effect of adding a lemma in the non-linear solver. This is used
+ * to specify how the state of the non-linear solver should update. This
+ * includes:
+ * - A set of secant points to record (for transcendental secant plane
+ * inferences).
+ */
+struct NlLemmaSideEffect
+{
+  NlLemmaSideEffect() {}
+  ~NlLemmaSideEffect() {}
+  /** secant points to add
+   *
+   * A member (tf, d, c) in this vector indicates that point c should be added
+   * to the list of secant points for an application of a transcendental
+   * function tf for Taylor degree d. This is used for incremental linearization
+   * for underapproximation (resp. overapproximations) of convex (resp.
+   * concave) regions of transcendental functions. For details, see
+   * Cimatti et al., CADE 2017.
+   */
+  std::vector<std::tuple<Node, unsigned, Node> > d_secantPoint;
+};
+
+}  // namespace arith
+}  // namespace theory
+}  // namespace CVC4
+
+#endif /* CVC4__THEORY__ARITH__NL_LEMMA_UTILS_H */
diff --git a/src/theory/arith/nonlinear_extension.cpp b/src/theory/arith/nonlinear_extension.cpp
index acae404ba..bcf7cae14 100644
--- a/src/theory/arith/nonlinear_extension.cpp
+++ b/src/theory/arith/nonlinear_extension.cpp
@@ -540,12 +540,20 @@ Node NonlinearExtension::mkMonomialRemFactor(
 }
 
 void NonlinearExtension::sendLemmas(const std::vector<Node>& out,
-                                    bool preprocess)
+                                    bool preprocess,
+                                    std::map<Node, NlLemmaSideEffect>& lemSE)
 {
+  std::map<Node, NlLemmaSideEffect>::iterator its;
   for (const Node& lem : out)
   {
     Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : " << lem << std::endl;
     d_containing.getOutputChannel().lemma(lem, false, preprocess);
+    // process the side effect
+    its = lemSE.find(lem);
+    if (its != lemSE.end())
+    {
+      processSideEffect(its->second);
+    }
     // add to cache if not preprocess
     if (!preprocess)
     {
@@ -554,6 +562,17 @@ void NonlinearExtension::sendLemmas(const std::vector<Node>& out,
   }
 }
 
+void NonlinearExtension::processSideEffect(const NlLemmaSideEffect& se)
+{
+  for (const std::tuple<Node, unsigned, Node>& sp : se.d_secantPoint)
+  {
+    Node tf = std::get<0>(sp);
+    unsigned d = std::get<1>(sp);
+    Node c = std::get<2>(sp);
+    d_secant_points[tf][d].push_back(c);
+  }
+}
+
 unsigned NonlinearExtension::filterLemma(Node lem, std::vector<Node>& out)
 {
   Trace("nl-ext-lemma-debug")
@@ -890,7 +909,8 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
                                       const std::vector<Node>& xts,
                                       std::vector<Node>& lems,
                                       std::vector<Node>& lemsPp,
-                                      std::vector<Node>& wlems)
+                                      std::vector<Node>& wlems,
+                                      std::map<Node, NlLemmaSideEffect>& lemSE)
 {
   d_ms_vars.clear();
   d_ms_proc.clear();
@@ -1262,7 +1282,7 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions,
   }
   if (options::nlExtTfTangentPlanes())
   {
-    lemmas = checkTranscendentalTangentPlanes();
+    lemmas = checkTranscendentalTangentPlanes(lemSE);
     filterLemmas(lemmas, wlems);
   }
   Trace("nl-ext") << "  ...finished with " << wlems.size() << " waiting lemmas."
@@ -1297,8 +1317,8 @@ void NonlinearExtension::check(Theory::Effort e) {
     // If we computed lemmas during collectModelInfo, send them now.
     if (!d_cmiLemmas.empty() || !d_cmiLemmasPp.empty())
     {
-      sendLemmas(d_cmiLemmas);
-      sendLemmas(d_cmiLemmasPp, true);
+      sendLemmas(d_cmiLemmas, false, d_cmiLemmasSE);
+      sendLemmas(d_cmiLemmasPp, true, d_cmiLemmasSE);
       return;
     }
     // Otherwise, we will answer SAT. The values that we approximated are
@@ -1318,8 +1338,10 @@ void NonlinearExtension::check(Theory::Effort e) {
   }
 }
 
-bool NonlinearExtension::modelBasedRefinement(std::vector<Node>& mlems,
-                                              std::vector<Node>& mlemsPp)
+bool NonlinearExtension::modelBasedRefinement(
+    std::vector<Node>& mlems,
+    std::vector<Node>& mlemsPp,
+    std::map<Node, NlLemmaSideEffect>& lemSE)
 {
   // get the assertions
   std::vector<Node> assertions;
@@ -1405,7 +1427,8 @@ bool NonlinearExtension::modelBasedRefinement(std::vector<Node>& mlems,
     if (!false_asserts.empty() || num_shared_wrong_value > 0)
     {
       complete_status = num_shared_wrong_value > 0 ? -1 : 0;
-      checkLastCall(assertions, false_asserts, xts, mlems, mlemsPp, wlems);
+      checkLastCall(
+          assertions, false_asserts, xts, mlems, mlemsPp, wlems, lemSE);
       if (!mlems.empty() || !mlemsPp.empty())
       {
         return true;
@@ -1522,10 +1545,11 @@ void NonlinearExtension::interceptModel(std::map<Node, Node>& arithModel)
   // run a last call effort check
   d_cmiLemmas.clear();
   d_cmiLemmasPp.clear();
+  d_cmiLemmasSE.clear();
   if (!d_builtModel.get())
   {
     Trace("nl-ext") << "interceptModel: do model-based refinement" << std::endl;
-    modelBasedRefinement(d_cmiLemmas, d_cmiLemmasPp);
+    modelBasedRefinement(d_cmiLemmas, d_cmiLemmasPp, d_cmiLemmasSE);
   }
   if (d_builtModel.get())
   {
@@ -3011,7 +3035,8 @@ std::vector<Node> NonlinearExtension::checkTranscendentalMonotonic() {
   return lemmas;
 }
 
-std::vector<Node> NonlinearExtension::checkTranscendentalTangentPlanes()
+std::vector<Node> NonlinearExtension::checkTranscendentalTangentPlanes(
+    std::map<Node, NlLemmaSideEffect>& lemSE)
 {
   std::vector<Node> lemmas;
   Trace("nl-ext") << "Get tangent plane lemmas for transcendental functions..."
@@ -3051,7 +3076,7 @@ std::vector<Node> NonlinearExtension::checkTranscendentalTangentPlanes()
         {
           Trace("nl-ext-tftp") << "- run at degree " << d << "..." << std::endl;
           unsigned prev = lemmas.size();
-          if (!checkTfTangentPlanesFun(tf, d, lemmas))
+          if (!checkTfTangentPlanesFun(tf, d, lemmas, lemSE))
           {
             Trace("nl-ext-tftp")
                 << "...fail, #lemmas = " << (lemmas.size() - prev) << std::endl;
@@ -3069,9 +3094,11 @@ std::vector<Node> NonlinearExtension::checkTranscendentalTangentPlanes()
   return lemmas;
 }
 
-bool NonlinearExtension::checkTfTangentPlanesFun(Node tf,
-                                                 unsigned d,
-                                                 std::vector<Node>& lemmas)
+bool NonlinearExtension::checkTfTangentPlanesFun(
+    Node tf,
+    unsigned d,
+    std::vector<Node>& lemmas,
+    std::map<Node, NlLemmaSideEffect>& lemSE)
 {
   Assert(d_model.isRefineableTfFun(tf));
 
@@ -3257,23 +3284,24 @@ bool NonlinearExtension::checkTfTangentPlanesFun(Node tf,
     Assert(std::find(
                d_secant_points[tf][d].begin(), d_secant_points[tf][d].end(), c)
            == d_secant_points[tf][d].end());
-    // insert into the vector
-    d_secant_points[tf][d].push_back(c);
+    // 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[tf][d];
+    spoints.push_back(c);
+
     // sort
     SortNlModel smv;
     smv.d_nlm = &d_model;
     smv.d_isConcrete = true;
-    std::sort(
-        d_secant_points[tf][d].begin(), d_secant_points[tf][d].end(), smv);
+    std::sort(spoints.begin(), spoints.end(), smv);
     // get the resulting index of c
     unsigned index =
-        std::find(
-            d_secant_points[tf][d].begin(), d_secant_points[tf][d].end(), c)
-        - d_secant_points[tf][d].begin();
+        std::find(spoints.begin(), spoints.end(), c) - spoints.begin();
     // bounds are the next closest upper/lower bound values
     if (index > 0)
     {
-      bounds[0] = d_secant_points[tf][d][index - 1];
+      bounds[0] = spoints[index - 1];
     }
     else
     {
@@ -3289,9 +3317,9 @@ bool NonlinearExtension::checkTfTangentPlanesFun(Node tf,
         bounds[0] = Rewriter::rewrite(nm->mkNode(MINUS, c, d_one));
       }
     }
-    if (index < d_secant_points[tf][d].size() - 1)
+    if (index < spoints.size() - 1)
     {
-      bounds[1] = d_secant_points[tf][d][index + 1];
+      bounds[1] = spoints[index + 1];
     }
     else
     {
@@ -3310,6 +3338,8 @@ bool NonlinearExtension::checkTfTangentPlanesFun(Node tf,
     Trace("nl-ext-tftp-debug2") << "...secant bounds are : " << bounds[0]
                                 << " ... " << bounds[1] << std::endl;
 
+    // the secant plane may be conjunction of 1-2 guarded inequalities
+    std::vector<Node> lemmaConj;
     for (unsigned s = 0; s < 2; s++)
     {
       // compute secant plane
@@ -3370,11 +3400,18 @@ bool NonlinearExtension::checkTfTangentPlanesFun(Node tf,
         lem = Rewriter::rewrite(lem);
         Trace("nl-ext-tftp-lemma") << "*** Secant plane lemma : " << lem
                                    << std::endl;
-        // Figure 3 : line 22
-        lemmas.push_back(lem);
+        lemmaConj.push_back(lem);
         Assert(d_model.computeAbstractModelValue(lem) == d_false);
       }
     }
+    // Figure 3 : line 22
+    Assert(!lemmaConj.empty());
+    Node lem =
+        lemmaConj.size() == 1 ? lemmaConj[0] : nm->mkNode(AND, lemmaConj);
+    lemmas.push_back(lem);
+    // The side effect says that if lem is added, then we should add the
+    // secant point c for (tf,d).
+    lemSE[lem].d_secantPoint.push_back(std::make_tuple(tf, d, c));
   }
   return false;
 }
diff --git a/src/theory/arith/nonlinear_extension.h b/src/theory/arith/nonlinear_extension.h
index 338ae6611..64a601cc5 100644
--- a/src/theory/arith/nonlinear_extension.h
+++ b/src/theory/arith/nonlinear_extension.h
@@ -34,6 +34,7 @@
 #include "context/context.h"
 #include "expr/kind.h"
 #include "expr/node.h"
+#include "theory/arith/nl_lemma_utils.h"
 #include "theory/arith/nl_model.h"
 #include "theory/arith/theory_arith.h"
 #include "theory/uf/equality_engine.h"
@@ -178,9 +179,13 @@ class NonlinearExtension {
    * Otherwise, it returns false. In the latter case, the model object d_model
    * may have information regarding how to construct a model, in the case that
    * we determined the problem is satisfiable.
+   *
+   * The argument lemSE is the "side effect" of the lemmas in mlems and mlemsPp
+   * (for details, see checkLastCall).
    */
   bool modelBasedRefinement(std::vector<Node>& mlems,
-                            std::vector<Node>& mlemsPp);
+                            std::vector<Node>& mlemsPp,
+                            std::map<Node, NlLemmaSideEffect>& lemSE);
   /** returns true if the multiset containing the
    * factors of monomial a is a subset of the multiset
    * containing the factors of monomial b.
@@ -230,13 +235,19 @@ class NonlinearExtension {
    * output channel as a last resort. In other words, only if we are not
    * able to establish SAT via a call to checkModel(...) should wlems be
    * considered. This set typically contains tangent plane lemmas.
+   *
+   * The argument lemSE is the "side effect" of the lemmas from the previous
+   * three calls. If a lemma is mapping to a side effect, it should be
+   * processed via a call to processSideEffect(...) immediately after the
+   * lemma is sent (if it is indeed sent on this call to check).
    */
   int checkLastCall(const std::vector<Node>& assertions,
                     const std::vector<Node>& false_asserts,
                     const std::vector<Node>& xts,
                     std::vector<Node>& lems,
                     std::vector<Node>& lemsPp,
-                    std::vector<Node>& wlems);
+                    std::vector<Node>& wlems,
+                    std::map<Node, NlLemmaSideEffect>& lemSE);
   //---------------------------------------term utilities
   static bool isArithKind(Kind k);
   static Node mkLit(Node a, Node b, int status, bool isAbsolute = false);
@@ -395,7 +406,11 @@ class NonlinearExtension {
   /**
    * Send lemmas in out on the output channel of theory of arithmetic.
    */
-  void sendLemmas(const std::vector<Node>& out, bool preprocess = false);
+  void sendLemmas(const std::vector<Node>& out,
+                  bool preprocess,
+                  std::map<Node, NlLemmaSideEffect>& lemSE);
+  /** Process side effect se */
+  void processSideEffect(const NlLemmaSideEffect& se);
 
   // Returns the NodeMultiset for an existing monomial.
   const NodeMultiset& getMonomialExponentMap(Node monomial) const;
@@ -480,6 +495,8 @@ class NonlinearExtension {
    */
   std::vector<Node> d_cmiLemmas;
   std::vector<Node> d_cmiLemmasPp;
+  /** the side effects of the above lemmas */
+  std::map<Node, NlLemmaSideEffect> d_cmiLemmasSE;
   /**
    * The approximations computed during collectModelInfo. For details, see
    * NlModel::getModelValueRepair.
@@ -872,64 +889,68 @@ class NonlinearExtension {
   std::vector<Node> checkTranscendentalMonotonic();
 
   /** check transcendental tangent planes
-  *
-  * Returns a set of valid theory lemmas, based on
-  * computing an "incremental linearization" of
-  * transcendental functions based on the model values
-  * of transcendental functions and their arguments.
-  * It is based on Figure 3 of "Satisfiability
-  * Modulo Transcendental Functions via Incremental
-  * Linearization" by Cimatti et al., CADE 2017.
-  * This schema is not terminating in general.
-  * It is not enabled by default, and can
-  * be enabled by --nl-ext-tf-tplanes.
-  *
-  * Example:
-  *
-  * Assume we have a term sin(y) where M( y ) = 1 where M is the current model.
-  * Note that:
-  *   sin(1) ~= .841471
-  *
-  * The Taylor series and remainder of sin(y) of degree 7 is
-  *   P_{7,sin(0)}( x ) = x + (-1/6)*x^3 + (1/20)*x^5
-  *   R_{7,sin(0),b}( x ) = (-1/5040)*x^7
-  *
-  * This gives us lower and upper bounds :
-  *   P_u( x ) = P_{7,sin(0)}( x ) + R_{7,sin(0),b}( x )
-  *     ...where note P_u( 1 ) = 4243/5040 ~= .841865
-  *   P_l( x ) = P_{7,sin(0)}( x ) - R_{7,sin(0),b}( x )
-  *     ...where note P_l( 1 ) = 4241/5040 ~= .841468
-  *
-  * Assume that M( sin(y) ) > P_u( 1 ).
-  * Since the concavity of sine in the region 0 < x < PI/2 is -1,
-  * we add a tangent plane refinement.
-  * The tangent plane at the point 1 in P_u is
-  * given by the formula:
-  *   T( x ) = P_u( 1 ) + ((d/dx)(P_u(x)))( 1 )*( x - 1 )
-  * We add the lemma:
-  *   ( 0 < y < PI/2 ) => sin( y ) <= T( y )
-  * which is:
-  *   ( 0 < y < PI/2 ) => sin( y ) <= (391/720)*(y - 2737/1506)
-  *
-  * Assume that M( sin(y) ) < P_u( 1 ).
-  * Since the concavity of sine in the region 0 < x < PI/2 is -1,
-  * we add a secant plane refinement for some constants ( l, u )
-  * such that 0 <= l < M( y ) < u <= PI/2. Assume we choose
-  * l = 0 and u = M( PI/2 ) = 150517/47912.
-  * The secant planes at point 1 for P_l
-  * are given by the formulas:
-  *   S_l( x ) = (x-l)*(P_l( l )-P_l(c))/(l-1) + P_l( l )
-  *   S_u( x ) = (x-u)*(P_l( u )-P_l(c))/(u-1) + P_l( u )
-  * We add the lemmas:
-  *   ( 0 < y < 1 ) => sin( y ) >= S_l( y )
-  *   ( 1 < y < PI/2 ) => sin( y ) >= S_u( y )
-  * which are:
-  *   ( 0 < y < 1 ) => (sin y) >= 4251/5040*y
-  *   ( 1 < y < PI/2 ) => (sin y) >= c1*(y+c2)
-  *     where c1, c2 are rationals (for brevity, omitted here)
-  *     such that c1 ~= .277 and c2 ~= 2.032.
-  */
-  std::vector<Node> checkTranscendentalTangentPlanes();
+   *
+   * Returns a set of valid theory lemmas, based on
+   * computing an "incremental linearization" of
+   * transcendental functions based on the model values
+   * of transcendental functions and their arguments.
+   * It is based on Figure 3 of "Satisfiability
+   * Modulo Transcendental Functions via Incremental
+   * Linearization" by Cimatti et al., CADE 2017.
+   * This schema is not terminating in general.
+   * It is not enabled by default, and can
+   * be enabled by --nl-ext-tf-tplanes.
+   *
+   * Example:
+   *
+   * Assume we have a term sin(y) where M( y ) = 1 where M is the current model.
+   * Note that:
+   *   sin(1) ~= .841471
+   *
+   * The Taylor series and remainder of sin(y) of degree 7 is
+   *   P_{7,sin(0)}( x ) = x + (-1/6)*x^3 + (1/20)*x^5
+   *   R_{7,sin(0),b}( x ) = (-1/5040)*x^7
+   *
+   * This gives us lower and upper bounds :
+   *   P_u( x ) = P_{7,sin(0)}( x ) + R_{7,sin(0),b}( x )
+   *     ...where note P_u( 1 ) = 4243/5040 ~= .841865
+   *   P_l( x ) = P_{7,sin(0)}( x ) - R_{7,sin(0),b}( x )
+   *     ...where note P_l( 1 ) = 4241/5040 ~= .841468
+   *
+   * Assume that M( sin(y) ) > P_u( 1 ).
+   * Since the concavity of sine in the region 0 < x < PI/2 is -1,
+   * we add a tangent plane refinement.
+   * The tangent plane at the point 1 in P_u is
+   * given by the formula:
+   *   T( x ) = P_u( 1 ) + ((d/dx)(P_u(x)))( 1 )*( x - 1 )
+   * We add the lemma:
+   *   ( 0 < y < PI/2 ) => sin( y ) <= T( y )
+   * which is:
+   *   ( 0 < y < PI/2 ) => sin( y ) <= (391/720)*(y - 2737/1506)
+   *
+   * Assume that M( sin(y) ) < P_u( 1 ).
+   * Since the concavity of sine in the region 0 < x < PI/2 is -1,
+   * we add a secant plane refinement for some constants ( l, u )
+   * such that 0 <= l < M( y ) < u <= PI/2. Assume we choose
+   * l = 0 and u = M( PI/2 ) = 150517/47912.
+   * The secant planes at point 1 for P_l
+   * are given by the formulas:
+   *   S_l( x ) = (x-l)*(P_l( l )-P_l(c))/(l-1) + P_l( l )
+   *   S_u( x ) = (x-u)*(P_l( u )-P_l(c))/(u-1) + P_l( u )
+   * We add the lemmas:
+   *   ( 0 < y < 1 ) => sin( y ) >= S_l( y )
+   *   ( 1 < y < PI/2 ) => sin( y ) >= S_u( y )
+   * which are:
+   *   ( 0 < y < 1 ) => (sin y) >= 4251/5040*y
+   *   ( 1 < y < PI/2 ) => (sin y) >= c1*(y+c2)
+   *     where c1, c2 are rationals (for brevity, omitted here)
+   *     such that c1 ~= .277 and c2 ~= 2.032.
+   *
+   * The argument lemSE is the "side effect" of the lemmas in the return
+   * value of this function (for details, see checkLastCall).
+   */
+  std::vector<Node> checkTranscendentalTangentPlanes(
+      std::map<Node, NlLemmaSideEffect>& lemSE);
   /** check transcendental function refinement for tf
    *
    * This method is called by the above method for each refineable
@@ -938,9 +959,13 @@ class NonlinearExtension {
    *
    * This runs Figure 3 of Cimatti et al., CADE 2017 for transcendental
    * function application tf for Taylor degree d. It may add a secant or
-   * tangent plane lemma to lems.
+   * tangent plane lemma to lems and its side effect (if one exists)
+   * to lemSE.
    */
-  bool checkTfTangentPlanesFun(Node tf, unsigned d, std::vector<Node>& lems);
+  bool checkTfTangentPlanesFun(Node tf,
+                               unsigned d,
+                               std::vector<Node>& lems,
+                               std::map<Node, NlLemmaSideEffect>& lemSE);
   //-------------------------------------------- end lemma schemas
 }; /* class NonlinearExtension */