GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/nl/transcendental/transcendental_state.cpp Lines: 181 193 93.8 %
Date: 2021-03-22 Branches: 493 1064 46.3 %

Line Exec Source
1
/*********************                                                        */
2
/*! \file transcendental_state.cpp
3
 ** \verbatim
4
 ** Top contributors (to current version):
5
 **   Gereon Kremer, Andrew Reynolds
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/transcendental_state.h"
16
17
#include "expr/proof.h"
18
#include "theory/arith/arith_utilities.h"
19
#include "theory/arith/inference_manager.h"
20
#include "theory/arith/nl/nl_model.h"
21
#include "theory/arith/nl/transcendental/taylor_generator.h"
22
#include "theory/rewriter.h"
23
24
namespace CVC4 {
25
namespace theory {
26
namespace arith {
27
namespace nl {
28
namespace transcendental {
29
30
4389
TranscendentalState::TranscendentalState(InferenceManager& im,
31
                                         NlModel& model,
32
                                         ProofNodeManager* pnm,
33
4389
                                         context::UserContext* c)
34
4389
    : d_im(im), d_model(model), d_pnm(pnm), d_ctx(c)
35
{
36
4389
  d_true = NodeManager::currentNM()->mkConst(true);
37
4389
  d_false = NodeManager::currentNM()->mkConst(false);
38
4389
  d_zero = NodeManager::currentNM()->mkConst(Rational(0));
39
4389
  d_one = NodeManager::currentNM()->mkConst(Rational(1));
40
4389
  d_neg_one = NodeManager::currentNM()->mkConst(Rational(-1));
41
4389
  if (d_pnm != nullptr)
42
  {
43
401
    d_proof.reset(new CDProofSet<CDProof>(d_pnm, d_ctx, "nl-trans"));
44
401
    d_proofChecker.reset(new TranscendentalProofRuleChecker());
45
401
    d_proofChecker->registerTo(pnm->getChecker());
46
  }
47
4389
}
48
49
2204
bool TranscendentalState::isProofEnabled() const
50
{
51
2204
  return d_proof.get() != nullptr;
52
}
53
54
327
CDProof* TranscendentalState::getProof()
55
{
56
327
  Assert(isProofEnabled());
57
327
  return d_proof->allocateProof(d_ctx);
58
}
59
60
2742
void TranscendentalState::init(const std::vector<Node>& xts,
61
                               std::vector<Node>& needsMaster)
62
{
63
2742
  d_funcCongClass.clear();
64
2742
  d_funcMap.clear();
65
2742
  d_tf_region.clear();
66
67
2742
  bool needPi = false;
68
  // for computing congruence
69
5484
  std::map<Kind, ArgTrie> argTrie;
70
23167
  for (std::size_t i = 0, xsize = xts.size(); i < xsize; ++i)
71
  {
72
    // Ignore if it is not a transcendental
73
20425
    if (!isTranscendentalKind(xts[i].getKind()))
74
    {
75
14201
      continue;
76
    }
77
12448
    Node a = xts[i];
78
6224
    Kind ak = a.getKind();
79
6224
    bool consider = true;
80
    // if we've already computed master for a
81
6224
    if (d_trMaster.find(a) != d_trMaster.end())
82
    {
83
      // a master has at least one slave
84
5753
      consider = (d_trSlaves.find(a) != d_trSlaves.end());
85
    }
86
    else
87
    {
88
471
      if (ak == Kind::SINE)
89
      {
90
        // always not a master
91
310
        consider = false;
92
      }
93
      else
94
      {
95
251
        for (const Node& ac : a)
96
        {
97
90
          if (isTranscendentalKind(ac.getKind()))
98
          {
99
            consider = false;
100
            break;
101
          }
102
        }
103
      }
104
471
      if (!consider)
105
      {
106
        // wait to assign a master below
107
310
        needsMaster.push_back(a);
108
      }
109
      else
110
      {
111
161
        d_trMaster[a] = a;
112
161
        d_trSlaves[a].insert(a);
113
      }
114
    }
115
6224
    if (ak == Kind::EXPONENTIAL || ak == Kind::SINE)
116
    {
117
5769
      needPi = needPi || (ak == Kind::SINE);
118
      // if we didn't indicate that it should be purified above
119
11538
      if (consider)
120
      {
121
3765
        ensureCongruence(a, argTrie);
122
      }
123
    }
124
455
    else if (ak == Kind::PI)
125
    {
126
455
      Assert(consider);
127
455
      needPi = true;
128
455
      d_funcMap[ak].push_back(a);
129
455
      d_funcCongClass[a].push_back(a);
130
    }
131
  }
132
  // initialize pi if necessary
133
2742
  if (needPi && d_pi.isNull())
134
  {
135
73
    mkPi();
136
73
    getCurrentPiBounds();
137
  }
138
139
2742
  if (Trace.isOn("nl-ext-mv"))
140
  {
141
    Trace("nl-ext-mv") << "Arguments of trancendental functions : "
142
                       << std::endl;
143
    for (std::pair<const Kind, std::vector<Node> >& tfl : d_funcMap)
144
    {
145
      Kind k = tfl.first;
146
      if (k == Kind::SINE || k == Kind::EXPONENTIAL)
147
      {
148
        for (const Node& tf : tfl.second)
149
        {
150
          Node v = tf[0];
151
          d_model.computeConcreteModelValue(v);
152
          d_model.computeAbstractModelValue(v);
153
          d_model.printModelValue("nl-ext-mv", v);
154
        }
155
      }
156
    }
157
  }
158
2742
}
159
160
3765
void TranscendentalState::ensureCongruence(TNode a,
161
                                           std::map<Kind, ArgTrie>& argTrie)
162
{
163
3765
  NodeManager* nm = NodeManager::currentNM();
164
7530
  std::vector<Node> repList;
165
7530
  for (const Node& ac : a)
166
  {
167
7530
    Node r = d_model.computeConcreteModelValue(ac);
168
3765
    repList.push_back(r);
169
  }
170
7530
  Node aa = argTrie[a.getKind()].add(a, repList);
171
3765
  if (aa != a)
172
  {
173
    // apply congruence to pairs of terms that are disequal and congruent
174
326
    Assert(aa.getNumChildren() == a.getNumChildren());
175
652
    Node mvaa = d_model.computeAbstractModelValue(a);
176
652
    Node mvaaa = d_model.computeAbstractModelValue(aa);
177
326
    if (mvaa != mvaaa)
178
    {
179
112
      std::vector<Node> exp;
180
112
      for (unsigned j = 0, size = a.getNumChildren(); j < size; j++)
181
      {
182
56
        exp.push_back(a[j].eqNode(aa[j]));
183
      }
184
112
      Node expn = exp.size() == 1 ? exp[0] : nm->mkNode(Kind::AND, exp);
185
112
      Node cong_lemma = expn.impNode(a.eqNode(aa));
186
56
      d_im.addPendingLemma(cong_lemma, InferenceId::ARITH_NL_CONGRUENCE);
187
    }
188
  }
189
  else
190
  {
191
    // new representative of congruence class
192
3439
    d_funcMap[a.getKind()].push_back(a);
193
  }
194
  // add to congruence class
195
3765
  d_funcCongClass[aa].push_back(a);
196
3765
}
197
198
73
void TranscendentalState::mkPi()
199
{
200
73
  NodeManager* nm = NodeManager::currentNM();
201
73
  if (d_pi.isNull())
202
  {
203
73
    d_pi = nm->mkNullaryOperator(nm->realType(), Kind::PI);
204
146
    d_pi_2 = Rewriter::rewrite(
205
219
        nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(1) / Rational(2))));
206
146
    d_pi_neg_2 = Rewriter::rewrite(
207
219
        nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(-1) / Rational(2))));
208
146
    d_pi_neg = Rewriter::rewrite(
209
219
        nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(-1))));
210
    // initialize bounds
211
73
    d_pi_bound[0] = nm->mkConst(Rational(103993) / Rational(33102));
212
73
    d_pi_bound[1] = nm->mkConst(Rational(104348) / Rational(33215));
213
  }
214
73
}
215
216
73
void TranscendentalState::getCurrentPiBounds()
217
{
218
73
  NodeManager* nm = NodeManager::currentNM();
219
  Node pi_lem = nm->mkNode(Kind::AND,
220
146
                           nm->mkNode(Kind::GEQ, d_pi, d_pi_bound[0]),
221
292
                           nm->mkNode(Kind::LEQ, d_pi, d_pi_bound[1]));
222
73
  CDProof* proof = nullptr;
223
73
  if (isProofEnabled())
224
  {
225
16
    proof = getProof();
226
48
    proof->addStep(
227
32
        pi_lem, PfRule::ARITH_TRANS_PI, {}, {d_pi_bound[0], d_pi_bound[1]});
228
  }
229
73
  d_im.addPendingLemma(pi_lem, InferenceId::ARITH_NL_T_PI_BOUND, proof);
230
73
}
231
232
115
std::pair<Node, Node> TranscendentalState::getClosestSecantPoints(TNode e,
233
                                                                  TNode center,
234
                                                                  unsigned d)
235
{
236
  // bounds are the minimum and maximum previous secant points
237
  // should not repeat secant points: secant lemmas should suffice to
238
  // rule out previous assignment
239
115
  Assert(std::find(
240
             d_secant_points[e][d].begin(), d_secant_points[e][d].end(), center)
241
         == d_secant_points[e][d].end());
242
  // Insert into the (temporary) vector. We do not update this vector
243
  // until we are sure this secant plane lemma has been processed. We do
244
  // this by mapping the lemma to a side effect below.
245
230
  std::vector<Node> spoints = d_secant_points[e][d];
246
115
  spoints.push_back(center);
247
248
  // sort
249
115
  sortByNlModel(spoints.begin(), spoints.end(), &d_model);
250
  // get the resulting index of c
251
  unsigned index =
252
115
      std::find(spoints.begin(), spoints.end(), center) - spoints.begin();
253
254
  // bounds are the next closest upper/lower bound values
255
230
  return {index > 0 ? spoints[index - 1] : Node(),
256
230
          index < spoints.size() - 1 ? spoints[index + 1] : Node()};
257
}
258
259
230
Node TranscendentalState::mkSecantPlane(
260
    TNode arg, TNode lower, TNode upper, TNode lval, TNode uval)
261
{
262
230
  NodeManager* nm = NodeManager::currentNM();
263
  // Figure 3: S_l( x ), S_u( x ) for s = 0,1
264
460
  Node rcoeff_n = Rewriter::rewrite(nm->mkNode(Kind::MINUS, lower, upper));
265
230
  Assert(rcoeff_n.isConst());
266
460
  Rational rcoeff = rcoeff_n.getConst<Rational>();
267
230
  Assert(rcoeff.sgn() != 0);
268
  Node res =
269
      nm->mkNode(Kind::PLUS,
270
                 lval,
271
920
                 nm->mkNode(Kind::MULT,
272
920
                            nm->mkNode(Kind::DIVISION,
273
460
                                       nm->mkNode(Kind::MINUS, lval, uval),
274
460
                                       nm->mkNode(Kind::MINUS, lower, upper)),
275
690
                            nm->mkNode(Kind::MINUS, arg, lower)));
276
460
  Trace("nl-trans") << "Creating secant plane for transcendental function of "
277
230
                    << arg << std::endl;
278
460
  Trace("nl-trans") << "\tfrom ( " << lower << " ; " << lval << " ) to ( "
279
230
                    << upper << " ; " << uval << " )" << std::endl;
280
230
  Trace("nl-trans") << "\t" << res << std::endl;
281
230
  Trace("nl-trans") << "\trewritten: " << Rewriter::rewrite(res) << std::endl;
282
460
  return res;
283
}
284
285
230
NlLemma TranscendentalState::mkSecantLemma(TNode lower,
286
                                           TNode upper,
287
                                           TNode lapprox,
288
                                           TNode uapprox,
289
                                           int csign,
290
                                           Convexity convexity,
291
                                           TNode tf,
292
                                           TNode splane,
293
                                           unsigned actual_d)
294
{
295
230
  NodeManager* nm = NodeManager::currentNM();
296
  // With respect to Figure 3, this is slightly different.
297
  // In particular, we chose b to be the model value of bounds[s],
298
  // which is a constant although bounds[s] may not be (e.g. if it
299
  // contains PI).
300
  // To ensure that c...b does not cross an inflection point,
301
  // we guard with the symbolic version of bounds[s].
302
  // This leads to lemmas e.g. of this form:
303
  //   ( c <= x <= PI/2 ) => ( sin(x) < ( P( b ) - P( c ) )*( x -
304
  //   b ) + P( b ) )
305
  // where b = (PI/2)^M, the current value of PI/2 in the model.
306
  // This is sound since we are guarded by the symbolic
307
  // representation of PI/2.
308
  Node antec_n = nm->mkNode(Kind::AND,
309
460
                            nm->mkNode(Kind::GEQ, tf[0], lower),
310
920
                            nm->mkNode(Kind::LEQ, tf[0], upper));
311
460
  Trace("nl-trans") << "Bound for secant plane: " << lower << " <= " << tf[0]
312
230
                    << " <= " << upper << std::endl;
313
230
  Trace("nl-trans") << "\t" << antec_n << std::endl;
314
  // Convex: actual value is below the secant
315
  // Concave: actual value is above the secant
316
  Node lem = nm->mkNode(
317
      Kind::IMPLIES,
318
      antec_n,
319
460
      nm->mkNode(
320
460
          convexity == Convexity::CONVEX ? Kind::LEQ : Kind::GEQ, tf, splane));
321
460
  Trace("nl-trans-lemma") << "*** Secant plane lemma (pre-rewrite) : " << lem
322
230
                          << std::endl;
323
230
  lem = Rewriter::rewrite(lem);
324
230
  Trace("nl-trans-lemma") << "*** Secant plane lemma : " << lem << std::endl;
325
230
  Assert(d_model.computeAbstractModelValue(lem) == d_false);
326
230
  CDProof* proof = nullptr;
327
230
  if (isProofEnabled())
328
  {
329
42
    proof = getProof();
330
42
    if (tf.getKind() == Kind::EXPONENTIAL)
331
    {
332
18
      if (csign == 1)
333
      {
334
80
        proof->addStep(
335
            lem,
336
            PfRule::ARITH_TRANS_EXP_APPROX_ABOVE_POS,
337
            {},
338
64
            {nm->mkConst<Rational>(2 * actual_d), tf[0], lower, upper});
339
      }
340
      else
341
      {
342
10
        proof->addStep(
343
            lem,
344
            PfRule::ARITH_TRANS_EXP_APPROX_ABOVE_NEG,
345
            {},
346
8
            {nm->mkConst<Rational>(2 * actual_d), tf[0], lower, upper});
347
      }
348
    }
349
24
    else if (tf.getKind() == Kind::SINE)
350
    {
351
24
      if (convexity == Convexity::CONCAVE)
352
      {
353
84
        proof->addStep(lem,
354
                       PfRule::ARITH_TRANS_SINE_APPROX_BELOW_POS,
355
                       {},
356
                       {nm->mkConst<Rational>(2 * actual_d),
357
                        tf[0],
358
                        lower,
359
                        upper,
360
                        lapprox,
361
                        uapprox
362
363
72
                       });
364
      }
365
      else
366
      {
367
84
        proof->addStep(lem,
368
                       PfRule::ARITH_TRANS_SINE_APPROX_ABOVE_NEG,
369
                       {},
370
                       {nm->mkConst<Rational>(2 * actual_d),
371
                        tf[0],
372
                        lower,
373
                        upper,
374
                        lapprox,
375
72
                        uapprox});
376
      }
377
    }
378
  }
379
  return NlLemma(
380
460
      InferenceId::ARITH_NL_T_SECANT, lem, LemmaProperty::NONE, proof);
381
}
382
383
115
void TranscendentalState::doSecantLemmas(const std::pair<Node, Node>& bounds,
384
                                         TNode poly_approx,
385
                                         TNode center,
386
                                         TNode cval,
387
                                         TNode tf,
388
                                         Convexity convexity,
389
                                         unsigned d,
390
                                         unsigned actual_d)
391
{
392
115
  int csign = center.getConst<Rational>().sgn();
393
230
  Trace("nl-trans") << "...do secant lemma with center " << center << " val "
394
115
                    << cval << " sign " << csign << std::endl;
395
230
  Trace("nl-trans") << "...secant bounds are : " << bounds.first << " ... "
396
115
                    << bounds.second << std::endl;
397
  // take the model value of lower (since may contain PI)
398
  // Make secant from bounds.first to center
399
230
  Node lower = d_model.computeAbstractModelValue(bounds.first);
400
230
  Trace("nl-trans") << "...model value of bound " << bounds.first << " is "
401
115
                    << lower << std::endl;
402
115
  if (lower != center)
403
  {
404
    // Figure 3 : P(l), P(u), for s = 0
405
    Node lval = Rewriter::rewrite(
406
230
        poly_approx.substitute(d_taylor.getTaylorVariable(), lower));
407
230
    Node splane = mkSecantPlane(tf[0], lower, center, lval, cval);
408
    NlLemma nlem = mkSecantLemma(
409
230
        lower, center, lval, cval, csign, convexity, tf, splane, actual_d);
410
    // The side effect says that if lem is added, then we should add the
411
    // secant point c for (tf,d).
412
115
    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, center));
413
115
    d_im.addPendingLemma(nlem, true);
414
  }
415
416
  // take the model value of upper (since may contain PI)
417
  // Make secant from center to bounds.second
418
230
  Node upper = d_model.computeAbstractModelValue(bounds.second);
419
230
  Trace("nl-trans") << "...model value of bound " << bounds.second << " is "
420
115
                    << upper << std::endl;
421
115
  if (center != upper)
422
  {
423
    // Figure 3 : P(l), P(u), for s = 1
424
    Node uval = Rewriter::rewrite(
425
230
        poly_approx.substitute(d_taylor.getTaylorVariable(), upper));
426
230
    Node splane = mkSecantPlane(tf[0], center, upper, cval, uval);
427
    NlLemma nlem = mkSecantLemma(
428
230
        center, upper, cval, uval, csign, convexity, tf, splane, actual_d);
429
    // The side effect says that if lem is added, then we should add the
430
    // secant point c for (tf,d).
431
115
    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, center));
432
115
    d_im.addPendingLemma(nlem, true);
433
  }
434
115
}
435
436
}  // namespace transcendental
437
}  // namespace nl
438
}  // namespace arith
439
}  // namespace theory
440
26676
}  // namespace CVC4