GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/nl/transcendental/sine_solver.cpp Lines: 222 224 99.1 %
Date: 2021-03-22 Branches: 547 1146 47.7 %

Line Exec Source
1
/*********************                                                        */
2
/*! \file sine_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/sine_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
namespace {
37
38
/**
39
 * Ensure a is in the main phase:
40
 *   -pi <= a <= pi
41
 */
42
284
inline Node mkValidPhase(TNode a, TNode pi)
43
{
44
  return mkBounded(
45
568
      NodeManager::currentNM()->mkNode(Kind::MULT, mkRationalNode(-1), pi),
46
      a,
47
852
      pi);
48
}
49
}  // namespace
50
51
4389
SineSolver::SineSolver(TranscendentalState* tstate) : d_data(tstate) {}
52
53
4386
SineSolver::~SineSolver() {}
54
55
142
void SineSolver::doPhaseShift(TNode a, TNode new_a, TNode y)
56
{
57
142
  NodeManager* nm = NodeManager::currentNM();
58
142
  Assert(a.getKind() == Kind::SINE);
59
142
  Trace("nl-ext-tf") << "Basis sine : " << new_a << " for " << a << std::endl;
60
142
  Assert(!d_data->d_pi.isNull());
61
284
  Node shift = nm->mkSkolem("s", nm->integerType(), "number of shifts");
62
  // TODO (cvc4-projects #47) : do not introduce shift here, instead needs model-based
63
  // refinement for constant shifts (cvc4-projects #1284)
64
  Node lem = nm->mkNode(
65
      Kind::AND,
66
284
      mkValidPhase(y, d_data->d_pi),
67
852
      nm->mkNode(Kind::ITE,
68
284
                 mkValidPhase(a[0], d_data->d_pi),
69
284
                 a[0].eqNode(y),
70
568
                 a[0].eqNode(nm->mkNode(Kind::PLUS,
71
                                        y,
72
568
                                        nm->mkNode(Kind::MULT,
73
284
                                                   nm->mkConst(Rational(2)),
74
                                                   shift,
75
142
                                                   d_data->d_pi)))),
76
852
      new_a.eqNode(a));
77
142
  CDProof* proof = nullptr;
78
142
  if (d_data->isProofEnabled())
79
  {
80
24
    proof = d_data->getProof();
81
24
    proof->addStep(lem, PfRule::ARITH_TRANS_SINE_SHIFT, {}, {a[0], y, shift});
82
  }
83
  // note we must do preprocess on this lemma
84
284
  Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : purify : " << lem
85
142
                        << std::endl;
86
142
  d_data->d_im.addPendingLemma(lem, InferenceId::ARITH_NL_T_PURIFY_ARG, proof);
87
142
}
88
89
2576
void SineSolver::checkInitialRefine()
90
{
91
2576
  NodeManager* nm = NodeManager::currentNM();
92
3452
  for (std::pair<const Kind, std::vector<Node> >& tfl : d_data->d_funcMap)
93
  {
94
876
    if (tfl.first != Kind::SINE)
95
    {
96
519
      continue;
97
    }
98
714
    Trace("nl-ext") << "Get initial (sine) refinement lemmas for "
99
357
                       "transcendental functions..."
100
357
                    << std::endl;
101
3241
    for (const Node& t : tfl.second)
102
    {
103
      // initial refinements
104
2884
      if (d_tf_initial_refine.find(t) == d_tf_initial_refine.end())
105
      {
106
250
        d_tf_initial_refine[t] = true;
107
        Node symn = nm->mkNode(Kind::SINE,
108
500
                               nm->mkNode(Kind::MULT, d_data->d_neg_one, t[0]));
109
250
        symn = Rewriter::rewrite(symn);
110
        // Can assume it is its own master since phase is split over 0,
111
        // hence  -pi <= t[0] <= pi implies -pi <= -t[0] <= pi.
112
250
        d_data->d_trMaster[symn] = symn;
113
250
        d_data->d_trSlaves[symn].insert(symn);
114
250
        Assert(d_data->d_trSlaves.find(t) != d_data->d_trSlaves.end());
115
116
        {
117
          // sine bounds: -1 <= sin(t) <= 1
118
          Node lem = nm->mkNode(Kind::AND,
119
500
                                nm->mkNode(Kind::LEQ, t, d_data->d_one),
120
1000
                                nm->mkNode(Kind::GEQ, t, d_data->d_neg_one));
121
250
          CDProof* proof = nullptr;
122
250
          if (d_data->isProofEnabled())
123
          {
124
40
            proof = d_data->getProof();
125
40
            proof->addStep(lem, PfRule::ARITH_TRANS_SINE_BOUNDS, {}, {t[0]});
126
          }
127
250
          d_data->d_im.addPendingLemma(
128
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
129
        }
130
        {
131
          // sine symmetry: sin(t) - sin(-t) = 0
132
500
          Node lem = nm->mkNode(Kind::PLUS, t, symn).eqNode(d_data->d_zero);
133
250
          CDProof* proof = nullptr;
134
250
          if (d_data->isProofEnabled())
135
          {
136
40
            proof = d_data->getProof();
137
40
            proof->addStep(lem, PfRule::ARITH_TRANS_SINE_SYMMETRY, {}, {t[0]});
138
          }
139
250
          d_data->d_im.addPendingLemma(
140
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
141
        }
142
        {
143
          // sine zero tangent:
144
          //   t > 0  =>  sin(t) < t
145
          //   t < 0  =>  sin(t) > t
146
          Node lem =
147
              nm->mkNode(Kind::AND,
148
1000
                         nm->mkNode(Kind::IMPLIES,
149
500
                                    nm->mkNode(Kind::GT, t[0], d_data->d_zero),
150
500
                                    nm->mkNode(Kind::LT, t, t[0])),
151
1000
                         nm->mkNode(Kind::IMPLIES,
152
500
                                    nm->mkNode(Kind::LT, t[0], d_data->d_zero),
153
1500
                                    nm->mkNode(Kind::GT, t, t[0])));
154
250
          CDProof* proof = nullptr;
155
250
          if (d_data->isProofEnabled())
156
          {
157
40
            proof = d_data->getProof();
158
80
            proof->addStep(
159
40
                lem, PfRule::ARITH_TRANS_SINE_TANGENT_ZERO, {}, {t[0]});
160
          }
161
250
          d_data->d_im.addPendingLemma(
162
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
163
        }
164
        {
165
          // sine pi tangent:
166
          //   t > -pi  =>  sin(t) > -pi-t
167
          //   t <  pi  =>  sin(t) <  pi-t
168
          Node lem = nm->mkNode(
169
              Kind::AND,
170
1000
              nm->mkNode(
171
                  Kind::IMPLIES,
172
500
                  nm->mkNode(Kind::GT, t[0], d_data->d_pi_neg),
173
500
                  nm->mkNode(Kind::GT,
174
                             t,
175
500
                             nm->mkNode(Kind::MINUS, d_data->d_pi_neg, t[0]))),
176
1000
              nm->mkNode(
177
                  Kind::IMPLIES,
178
500
                  nm->mkNode(Kind::LT, t[0], d_data->d_pi),
179
500
                  nm->mkNode(Kind::LT,
180
                             t,
181
1500
                             nm->mkNode(Kind::MINUS, d_data->d_pi, t[0]))));
182
250
          CDProof* proof = nullptr;
183
250
          if (d_data->isProofEnabled())
184
          {
185
40
            proof = d_data->getProof();
186
80
            proof->addStep(
187
40
                lem, PfRule::ARITH_TRANS_SINE_TANGENT_PI, {}, {t[0]});
188
          }
189
250
          d_data->d_im.addPendingLemma(
190
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
191
        }
192
        {
193
          Node lem =
194
              nm->mkNode(Kind::AND,
195
                         // sign
196
1000
                         nm->mkNode(Kind::EQUAL,
197
500
                                    nm->mkNode(Kind::LT, t[0], d_data->d_zero),
198
500
                                    nm->mkNode(Kind::LT, t, d_data->d_zero)),
199
                         // zero val
200
1000
                         nm->mkNode(Kind::EQUAL,
201
500
                                    nm->mkNode(Kind::GT, t[0], d_data->d_zero),
202
1500
                                    nm->mkNode(Kind::GT, t, d_data->d_zero)));
203
250
          d_data->d_im.addPendingLemma(lem,
204
                                       InferenceId::ARITH_NL_T_INIT_REFINE);
205
        }
206
      }
207
    }
208
  }
209
2576
}
210
211
1574
void SineSolver::checkMonotonic()
212
{
213
214
1574
  auto it = d_data->d_funcMap.find(Kind::SINE);
215
1574
  if (it == d_data->d_funcMap.end())
216
  {
217
1387
    Trace("nl-ext-exp") << "No sine terms" << std::endl;
218
2774
    return;
219
  }
220
374
  Trace("nl-ext")
221
187
      << "Get monotonicity lemmas for (sine) transcendental functions..."
222
187
      << std::endl;
223
224
  // sort arguments of all transcendentals
225
374
  std::vector<Node> tf_args;
226
374
  std::map<Node, Node> tf_arg_to_term;
227
228
2039
  for (const Node& tf : it->second)
229
  {
230
3704
    Node a = tf[0];
231
3704
    Node mvaa = d_data->d_model.computeAbstractModelValue(a);
232
1852
    if (mvaa.isConst())
233
    {
234
1852
      Trace("nl-ext-tf-mono-debug") << "...tf term : " << a << std::endl;
235
1852
      tf_args.push_back(a);
236
1852
      tf_arg_to_term[a] = tf;
237
    }
238
  }
239
240
187
  if (tf_args.empty())
241
  {
242
    return;
243
  }
244
245
187
  sortByNlModel(
246
187
      tf_args.begin(), tf_args.end(), &d_data->d_model, true, false, true);
247
248
187
  std::vector<Node> mpoints = {d_data->d_pi,
249
187
                               d_data->d_pi_2,
250
187
                               d_data->d_zero,
251
187
                               d_data->d_pi_neg_2,
252
935
                               d_data->d_pi_neg};
253
374
  std::vector<Node> mpoints_vals;
254
255
  // get model values for points
256
1122
  for (const auto& point : mpoints)
257
  {
258
935
    mpoints_vals.emplace_back(d_data->d_model.computeAbstractModelValue(point));
259
935
    Assert(mpoints_vals.back().isConst());
260
  }
261
262
187
  unsigned mdir_index = 0;
263
187
  int monotonic_dir = -1;
264
374
  Node mono_bounds[2];
265
374
  Node targ, targval, t, tval;
266
2039
  for (const auto& sarg : tf_args)
267
  {
268
3704
    Node sargval = d_data->d_model.computeAbstractModelValue(sarg);
269
1852
    Assert(sargval.isConst());
270
3704
    Node s = tf_arg_to_term[sarg];
271
3704
    Node sval = d_data->d_model.computeAbstractModelValue(s);
272
1852
    Assert(sval.isConst());
273
274
    // increment to the proper monotonicity region
275
1852
    bool increment = true;
276
6858
    while (increment && mdir_index < mpoints.size())
277
    {
278
2503
      increment = false;
279
5006
      Node pval = mpoints_vals[mdir_index];
280
2503
      Assert(pval.isConst());
281
2503
      if (sargval.getConst<Rational>() < pval.getConst<Rational>())
282
      {
283
651
        increment = true;
284
1302
        Trace("nl-ext-tf-mono")
285
651
            << "...increment at " << sarg << " since model value is less than "
286
651
            << mpoints[mdir_index] << std::endl;
287
      }
288
2503
      if (increment)
289
      {
290
651
        tval = Node::null();
291
651
        mono_bounds[1] = mpoints[mdir_index];
292
651
        mdir_index++;
293
651
        monotonic_dir = regionToMonotonicityDir(mdir_index);
294
651
        if (mdir_index < mpoints.size())
295
        {
296
651
          mono_bounds[0] = mpoints[mdir_index];
297
        }
298
        else
299
        {
300
          mono_bounds[0] = Node::null();
301
        }
302
      }
303
    }
304
    // store the concavity region
305
1852
    d_data->d_tf_region[s] = mdir_index;
306
3704
    Trace("nl-ext-concavity")
307
1852
        << "Transcendental function " << s << " is in region #" << mdir_index;
308
1852
    Trace("nl-ext-concavity") << ", arg model value = " << sargval << std::endl;
309
310
1852
    if (!tval.isNull())
311
    {
312
1332
      NodeManager* nm = NodeManager::currentNM();
313
2664
      Node mono_lem;
314
1332
      if (monotonic_dir == 1
315
1332
          && sval.getConst<Rational>() > tval.getConst<Rational>())
316
      {
317
216
        mono_lem = nm->mkNode(Kind::IMPLIES,
318
144
                              nm->mkNode(Kind::GEQ, targ, sarg),
319
144
                              nm->mkNode(Kind::GEQ, t, s));
320
      }
321
1260
      else if (monotonic_dir == -1
322
1260
               && sval.getConst<Rational>() < tval.getConst<Rational>())
323
      {
324
96
        mono_lem = nm->mkNode(Kind::IMPLIES,
325
64
                              nm->mkNode(Kind::LEQ, targ, sarg),
326
64
                              nm->mkNode(Kind::LEQ, t, s));
327
      }
328
1332
      if (!mono_lem.isNull())
329
      {
330
104
        if (!mono_bounds[0].isNull())
331
        {
332
104
          Assert(!mono_bounds[1].isNull());
333
312
          mono_lem = nm->mkNode(
334
              Kind::IMPLIES,
335
416
              nm->mkNode(Kind::AND,
336
208
                         mkBounded(mono_bounds[0], targ, mono_bounds[1]),
337
208
                         mkBounded(mono_bounds[0], sarg, mono_bounds[1])),
338
              mono_lem);
339
        }
340
208
        Trace("nl-ext-tf-mono")
341
104
            << "Monotonicity lemma : " << mono_lem << std::endl;
342
343
104
        d_data->d_im.addPendingLemma(mono_lem,
344
                                     InferenceId::ARITH_NL_T_MONOTONICITY);
345
      }
346
    }
347
    // store the previous values
348
1852
    targ = sarg;
349
1852
    targval = sargval;
350
1852
    t = s;
351
1852
    tval = sval;
352
  }
353
}
354
355
32
void SineSolver::doTangentLemma(
356
    TNode e, TNode c, TNode poly_approx, int region, std::uint64_t d)
357
{
358
32
  NodeManager* nm = NodeManager::currentNM();
359
360
  // compute tangent plane
361
  // Figure 3: T( x )
362
  // We use zero slope tangent planes, since the concavity of the Taylor
363
  // approximation cannot be easily established.
364
32
  Convexity convexity = regionToConvexity(region);
365
32
  int mdir = regionToMonotonicityDir(region);
366
32
  bool usec = (mdir == 1) == (convexity == Convexity::CONCAVE);
367
  Node lem = nm->mkNode(
368
      Kind::IMPLIES,
369
128
      nm->mkNode(
370
          Kind::AND,
371
64
          nm->mkNode(
372
64
              Kind::GEQ, e[0], usec ? regionToLowerBound(region) : Node(c)),
373
64
          nm->mkNode(
374
64
              Kind::LEQ, e[0], usec ? Node(c) : regionToUpperBound(region))),
375
64
      nm->mkNode(convexity == Convexity::CONVEX ? Kind::GEQ : Kind::LEQ,
376
                 e,
377
128
                 poly_approx));
378
379
64
  Trace("nl-ext-sine") << "*** Tangent plane lemma (pre-rewrite): " << lem
380
32
                       << std::endl;
381
32
  lem = Rewriter::rewrite(lem);
382
32
  Trace("nl-ext-sine") << "*** Tangent plane lemma : " << lem << std::endl;
383
32
  Assert(d_data->d_model.computeAbstractModelValue(lem) == d_data->d_false);
384
  // Figure 3 : line 9
385
32
  CDProof* proof = nullptr;
386
32
  if (d_data->isProofEnabled())
387
  {
388
8
    proof = d_data->getProof();
389
8
    if (convexity == Convexity::CONVEX)
390
    {
391
4
      if (usec)
392
      {
393
6
        proof->addStep(lem,
394
                       PfRule::ARITH_TRANS_SINE_APPROX_BELOW_NEG,
395
                       {},
396
                       {nm->mkConst<Rational>(2 * d),
397
                        e[0],
398
                        c,
399
                        regionToLowerBound(region),
400
5
                        c});
401
      }
402
      else
403
      {
404
18
        proof->addStep(lem,
405
                       PfRule::ARITH_TRANS_SINE_APPROX_BELOW_NEG,
406
                       {},
407
                       {nm->mkConst<Rational>(2 * d),
408
                        e[0],
409
                        c,
410
                        c,
411
15
                        regionToUpperBound(region)});
412
      }
413
    }
414
    else
415
    {
416
4
      if (usec)
417
      {
418
18
        proof->addStep(lem,
419
                       PfRule::ARITH_TRANS_SINE_APPROX_ABOVE_POS,
420
                       {},
421
                       {nm->mkConst<Rational>(2 * d),
422
                        e[0],
423
                        c,
424
                        regionToLowerBound(region),
425
15
                        c});
426
      }
427
      else
428
      {
429
6
        proof->addStep(lem,
430
                       PfRule::ARITH_TRANS_SINE_APPROX_ABOVE_POS,
431
                       {},
432
                       {nm->mkConst<Rational>(2 * d),
433
                        e[0],
434
                        c,
435
                        c,
436
5
                        regionToUpperBound(region)});
437
      }
438
    }
439
  }
440
32
  d_data->d_im.addPendingLemma(
441
      lem, InferenceId::ARITH_NL_T_TANGENT, proof, true);
442
32
}
443
444
72
void SineSolver::doSecantLemmas(TNode e,
445
                                TNode poly_approx,
446
                                TNode c,
447
                                TNode poly_approx_c,
448
                                unsigned d,
449
                                unsigned actual_d,
450
                                int region)
451
{
452
72
  d_data->doSecantLemmas(getSecantBounds(e, c, d, region),
453
                         poly_approx,
454
                         c,
455
                         poly_approx_c,
456
                         e,
457
                         regionToConvexity(region),
458
                         d,
459
                         actual_d);
460
72
}
461
462
72
std::pair<Node, Node> SineSolver::getSecantBounds(TNode e,
463
                                                  TNode c,
464
                                                  unsigned d,
465
                                                  int region)
466
{
467
72
  std::pair<Node, Node> bounds = d_data->getClosestSecantPoints(e, c, d);
468
469
  // Check if we already have neighboring secant points
470
72
  if (bounds.first.isNull())
471
  {
472
    // lower boundary point for this concavity region
473
72
    bounds.first = regionToLowerBound(region);
474
  }
475
72
  if (bounds.second.isNull())
476
  {
477
    // upper boundary point for this concavity region
478
72
    bounds.second = regionToUpperBound(region);
479
  }
480
72
  return bounds;
481
}
482
483
}  // namespace transcendental
484
}  // namespace nl
485
}  // namespace arith
486
}  // namespace theory
487
26676
}  // namespace CVC4