GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/nl/transcendental/exponential_solver.cpp Lines: 113 121 93.4 %
Date: 2021-03-22 Branches: 248 542 45.8 %

Line Exec Source
1
/*********************                                                        */
2
/*! \file exponential_solver.cpp
3
 ** \verbatim
4
 ** Top contributors (to current version):
5
 **   Gereon Kremer, Andrew Reynolds, Tim King
6
 ** This file is part of the CVC4 project.
7
 ** Copyright (c) 2009-2021 by the authors listed in the file AUTHORS
8
 ** in the top-level source directory and their institutional affiliations.
9
 ** All rights reserved.  See the file COPYING in the top-level source
10
 ** directory for licensing information.\endverbatim
11
 **
12
 ** \brief Implementation of solver for handling transcendental functions.
13
 **/
14
15
#include "theory/arith/nl/transcendental/exponential_solver.h"
16
17
#include <cmath>
18
#include <set>
19
20
#include "expr/node_algorithm.h"
21
#include "expr/node_builder.h"
22
#include "expr/proof.h"
23
#include "options/arith_options.h"
24
#include "theory/arith/arith_msum.h"
25
#include "theory/arith/arith_utilities.h"
26
#include "theory/arith/inference_manager.h"
27
#include "theory/arith/nl/nl_model.h"
28
#include "theory/arith/nl/transcendental/transcendental_state.h"
29
#include "theory/rewriter.h"
30
31
namespace CVC4 {
32
namespace theory {
33
namespace arith {
34
namespace nl {
35
namespace transcendental {
36
37
4389
ExponentialSolver::ExponentialSolver(TranscendentalState* tstate)
38
4389
    : d_data(tstate)
39
{
40
4389
}
41
42
4386
ExponentialSolver::~ExponentialSolver() {}
43
44
void ExponentialSolver::doPurification(TNode a, TNode new_a, TNode y)
45
{
46
  NodeManager* nm = NodeManager::currentNM();
47
  // do both equalities to ensure that new_a becomes a preregistered term
48
  Node lem = nm->mkNode(Kind::AND, a.eqNode(new_a), a[0].eqNode(y));
49
  // note we must do preprocess on this lemma
50
  Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : purify : " << lem
51
                        << std::endl;
52
  d_data->d_im.addPendingLemma(lem, InferenceId::ARITH_NL_T_PURIFY_ARG);
53
}
54
55
2576
void ExponentialSolver::checkInitialRefine()
56
{
57
2576
  NodeManager* nm = NodeManager::currentNM();
58
3452
  for (std::pair<const Kind, std::vector<Node> >& tfl : d_data->d_funcMap)
59
  {
60
876
    if (tfl.first != Kind::EXPONENTIAL)
61
    {
62
715
      continue;
63
    }
64
322
    Trace("nl-ext")
65
161
        << "Get initial (exp) refinement lemmas for transcendental functions..."
66
161
        << std::endl;
67
161
    Assert(tfl.first == Kind::EXPONENTIAL);
68
516
    for (const Node& t : tfl.second)
69
    {
70
      // initial refinements
71
355
      if (d_tf_initial_refine.find(t) == d_tf_initial_refine.end())
72
      {
73
86
        d_tf_initial_refine[t] = true;
74
        {
75
          // exp is always positive: exp(t) > 0
76
172
          Node lem = nm->mkNode(Kind::GT, t, d_data->d_zero);
77
86
          CDProof* proof = nullptr;
78
86
          if (d_data->isProofEnabled())
79
          {
80
18
            proof = d_data->getProof();
81
18
            proof->addStep(lem, PfRule::ARITH_TRANS_EXP_POSITIVITY, {}, {t[0]});
82
          }
83
86
          d_data->d_im.addPendingLemma(
84
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
85
        }
86
        {
87
          // exp at zero: (t = 0) <=> (exp(t) = 1)
88
          Node lem = nm->mkNode(Kind::EQUAL,
89
172
                                t[0].eqNode(d_data->d_zero),
90
344
                                t.eqNode(d_data->d_one));
91
86
          CDProof* proof = nullptr;
92
86
          if (d_data->isProofEnabled())
93
          {
94
18
            proof = d_data->getProof();
95
18
            proof->addStep(lem, PfRule::ARITH_TRANS_EXP_ZERO, {}, {t[0]});
96
          }
97
86
          d_data->d_im.addPendingLemma(
98
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
99
        }
100
        {
101
          // exp on negative values: (t < 0) <=> (exp(t) < 1)
102
          Node lem = nm->mkNode(Kind::EQUAL,
103
172
                                nm->mkNode(Kind::LT, t[0], d_data->d_zero),
104
344
                                nm->mkNode(Kind::LT, t, d_data->d_one));
105
86
          CDProof* proof = nullptr;
106
86
          if (d_data->isProofEnabled())
107
          {
108
18
            proof = d_data->getProof();
109
18
            proof->addStep(lem, PfRule::ARITH_TRANS_EXP_NEG, {}, {t[0]});
110
          }
111
86
          d_data->d_im.addPendingLemma(
112
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
113
        }
114
        {
115
          // exp on positive values: (t <= 0) or (exp(t) > t+1)
116
          Node lem = nm->mkNode(
117
              Kind::OR,
118
172
              nm->mkNode(Kind::LEQ, t[0], d_data->d_zero),
119
172
              nm->mkNode(
120
516
                  Kind::GT, t, nm->mkNode(Kind::PLUS, t[0], d_data->d_one)));
121
86
          CDProof* proof = nullptr;
122
86
          if (d_data->isProofEnabled())
123
          {
124
18
            proof = d_data->getProof();
125
18
            proof->addStep(lem, PfRule::ARITH_TRANS_EXP_SUPER_LIN, {}, {t[0]});
126
          }
127
86
          d_data->d_im.addPendingLemma(
128
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
129
        }
130
      }
131
    }
132
  }
133
2576
}
134
135
1574
void ExponentialSolver::checkMonotonic()
136
{
137
1574
  auto it = d_data->d_funcMap.find(Kind::EXPONENTIAL);
138
1574
  if (it == d_data->d_funcMap.end())
139
  {
140
1482
    Trace("nl-ext-exp") << "No exponential terms" << std::endl;
141
2964
    return;
142
  }
143
144
184
  Trace("nl-ext")
145
92
      << "Get monotonicity lemmas for (exp) transcendental functions..."
146
92
      << std::endl;
147
  // sort arguments of all transcendentals
148
184
  std::vector<Node> tf_args;
149
184
  std::map<Node, Node> tf_arg_to_term;
150
151
269
  for (const Node& tf : it->second)
152
  {
153
354
    Node a = tf[0];
154
354
    Node mvaa = d_data->d_model.computeAbstractModelValue(a);
155
177
    if (mvaa.isConst())
156
    {
157
177
      Trace("nl-ext-exp") << "...tf term : " << a << std::endl;
158
177
      tf_args.push_back(a);
159
177
      tf_arg_to_term[a] = tf;
160
    }
161
  }
162
163
92
  if (tf_args.empty())
164
  {
165
    return;
166
  }
167
168
92
  sortByNlModel(
169
92
      tf_args.begin(), tf_args.end(), &d_data->d_model, true, false, true);
170
171
184
  Node targ, targval, t, tval;
172
269
  for (const auto& sarg : tf_args)
173
  {
174
354
    Node sargval = d_data->d_model.computeAbstractModelValue(sarg);
175
177
    Assert(sargval.isConst());
176
354
    Node s = tf_arg_to_term[sarg];
177
354
    Node sval = d_data->d_model.computeAbstractModelValue(s);
178
177
    Assert(sval.isConst());
179
180
    // store the concavity region
181
177
    d_data->d_tf_region[s] = 1;
182
177
    Trace("nl-ext-concavity") << ", arg model value = " << sargval << std::endl;
183
184
177
    if (!tval.isNull() && sval.getConst<Rational>() > tval.getConst<Rational>())
185
    {
186
13
      NodeManager* nm = NodeManager::currentNM();
187
      Node mono_lem = nm->mkNode(Kind::IMPLIES,
188
26
                                 nm->mkNode(Kind::GEQ, targ, sarg),
189
52
                                 nm->mkNode(Kind::GEQ, t, s));
190
13
      Trace("nl-ext-exp") << "Monotonicity lemma : " << mono_lem << std::endl;
191
192
13
      d_data->d_im.addPendingLemma(mono_lem,
193
                                   InferenceId::ARITH_NL_T_MONOTONICITY);
194
    }
195
    // store the previous values
196
177
    targ = sarg;
197
177
    targval = sargval;
198
177
    t = s;
199
177
    tval = sval;
200
  }
201
}
202
203
56
void ExponentialSolver::doTangentLemma(TNode e,
204
                                       TNode c,
205
                                       TNode poly_approx,
206
                                       std::uint64_t d)
207
{
208
56
  NodeManager* nm = NodeManager::currentNM();
209
  // compute tangent plane
210
  // Figure 3: T( x )
211
  // We use zero slope tangent planes, since the concavity of the Taylor
212
  // approximation cannot be easily established.
213
  // Tangent plane is valid in the interval [c,u).
214
  Node lem = nm->mkNode(Kind::IMPLIES,
215
112
                        nm->mkNode(Kind::GEQ, e[0], c),
216
224
                        nm->mkNode(Kind::GEQ, e, poly_approx));
217
112
  Trace("nl-ext-exp") << "*** Tangent plane lemma (pre-rewrite): " << lem
218
56
                      << std::endl;
219
56
  lem = Rewriter::rewrite(lem);
220
56
  Trace("nl-ext-exp") << "*** Tangent plane lemma : " << lem << std::endl;
221
56
  Assert(d_data->d_model.computeAbstractModelValue(lem) == d_data->d_false);
222
  // Figure 3 : line 9
223
56
  CDProof* proof = nullptr;
224
56
  if (d_data->isProofEnabled())
225
  {
226
5
    proof = d_data->getProof();
227
15
    proof->addStep(lem,
228
                   PfRule::ARITH_TRANS_EXP_APPROX_BELOW,
229
                   {},
230
10
                   {nm->mkConst<Rational>(d), e[0]});
231
  }
232
56
  d_data->d_im.addPendingLemma(
233
      lem, InferenceId::ARITH_NL_T_TANGENT, proof, true);
234
56
}
235
236
43
void ExponentialSolver::doSecantLemmas(TNode e,
237
                                       TNode poly_approx,
238
                                       TNode center,
239
                                       TNode cval,
240
                                       unsigned d,
241
                                       unsigned actual_d)
242
{
243
43
  d_data->doSecantLemmas(getSecantBounds(e, center, d),
244
                         poly_approx,
245
                         center,
246
                         cval,
247
                         e,
248
                         Convexity::CONVEX,
249
                         d,
250
                         actual_d);
251
43
}
252
253
43
std::pair<Node, Node> ExponentialSolver::getSecantBounds(TNode e,
254
                                                         TNode center,
255
                                                         unsigned d)
256
{
257
43
  std::pair<Node, Node> bounds = d_data->getClosestSecantPoints(e, center, d);
258
259
  // Check if we already have neighboring secant points
260
43
  if (bounds.first.isNull())
261
  {
262
    // pick c-1
263
86
    bounds.first = Rewriter::rewrite(
264
129
        NodeManager::currentNM()->mkNode(Kind::MINUS, center, d_data->d_one));
265
  }
266
43
  if (bounds.second.isNull())
267
  {
268
    // pick c+1
269
86
    bounds.second = Rewriter::rewrite(
270
129
        NodeManager::currentNM()->mkNode(Kind::PLUS, center, d_data->d_one));
271
  }
272
43
  return bounds;
273
}
274
275
}  // namespace transcendental
276
}  // namespace nl
277
}  // namespace arith
278
}  // namespace theory
279
26676
}  // namespace CVC4