GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/nl/transcendental/transcendental_solver.cpp Lines: 178 207 86.0 %
Date: 2021-03-22 Branches: 347 829 41.9 %

Line Exec Source
1
/*********************                                                        */
2
/*! \file transcendental_solver.cpp
3
 ** \verbatim
4
 ** Top contributors (to current version):
5
 **   Andrew Reynolds, Gereon Kremer, 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/transcendental_solver.h"
16
17
#include <cmath>
18
#include <set>
19
20
#include "expr/node_algorithm.h"
21
#include "expr/node_builder.h"
22
#include "options/arith_options.h"
23
#include "theory/arith/arith_msum.h"
24
#include "theory/arith/arith_utilities.h"
25
#include "theory/arith/inference_manager.h"
26
#include "theory/arith/nl/nl_model.h"
27
#include "theory/arith/nl/transcendental/taylor_generator.h"
28
#include "theory/rewriter.h"
29
30
using namespace CVC4::kind;
31
32
namespace CVC4 {
33
namespace theory {
34
namespace arith {
35
namespace nl {
36
namespace transcendental {
37
38
4389
TranscendentalSolver::TranscendentalSolver(InferenceManager& im,
39
                                           NlModel& m,
40
                                           ProofNodeManager* pnm,
41
4389
                                           context::UserContext* c)
42
4389
    : d_tstate(im, m, pnm, c), d_expSlv(&d_tstate), d_sineSlv(&d_tstate)
43
{
44
8778
  d_taylor_degree = options::nlExtTfTaylorDegree();
45
4389
}
46
47
4386
TranscendentalSolver::~TranscendentalSolver() {}
48
49
2742
void TranscendentalSolver::initLastCall(const std::vector<Node>& xts)
50
{
51
5380
  std::vector<Node> needsMaster;
52
2742
  d_tstate.init(xts, needsMaster);
53
54
2742
  if (d_tstate.d_im.hasUsed()) {
55
104
    return;
56
  }
57
58
2638
  NodeManager* nm = NodeManager::currentNM();
59
2780
  for (const Node& a : needsMaster)
60
  {
61
    // should not have processed this already
62
142
    Assert(d_tstate.d_trMaster.find(a) == d_tstate.d_trMaster.end());
63
142
    Kind k = a.getKind();
64
142
    Assert(k == Kind::SINE || k == Kind::EXPONENTIAL);
65
    Node y =
66
284
        nm->mkSkolem("y", nm->realType(), "phase shifted trigonometric arg");
67
284
    Node new_a = nm->mkNode(k, y);
68
142
    d_tstate.d_trSlaves[new_a].insert(new_a);
69
142
    d_tstate.d_trSlaves[new_a].insert(a);
70
142
    d_tstate.d_trMaster[a] = new_a;
71
142
    d_tstate.d_trMaster[new_a] = new_a;
72
142
    switch (k)
73
    {
74
142
      case Kind::SINE: d_sineSlv.doPhaseShift(a, new_a, y); break;
75
      case Kind::EXPONENTIAL: d_expSlv.doPurification(a, new_a, y); break;
76
      default: AlwaysAssert(false) << "Unexpected Kind " << k;
77
    }
78
  }
79
}
80
81
222
bool TranscendentalSolver::preprocessAssertionsCheckModel(
82
    std::vector<Node>& assertions)
83
{
84
444
  std::vector<Node> pvars;
85
444
  std::vector<Node> psubs;
86
610
  for (const std::pair<const Node, Node>& tb : d_tstate.d_trMaster)
87
  {
88
388
    pvars.push_back(tb.first);
89
388
    psubs.push_back(tb.second);
90
  }
91
92
  // initialize representation of assertions
93
444
  std::vector<Node> passertions;
94
3152
  for (const Node& a : assertions)
95
96
  {
97
5860
    Node pa = a;
98
2930
    if (!pvars.empty())
99
    {
100
1593
      pa = arithSubstitute(pa, pvars, psubs);
101
1593
      pa = Rewriter::rewrite(pa);
102
    }
103
2930
    if (!pa.isConst() || !pa.getConst<bool>())
104
    {
105
2890
      Trace("nl-ext-cm-assert") << "- assert : " << pa << std::endl;
106
2890
      passertions.push_back(pa);
107
    }
108
  }
109
  // get model bounds for all transcendental functions
110
444
  Trace("nl-ext-cm") << "----- Get bounds for transcendental functions..."
111
222
                     << std::endl;
112
424
  for (std::pair<const Kind, std::vector<Node> >& tfs : d_tstate.d_funcMap)
113
  {
114
532
    for (const Node& tf : tfs.second)
115
    {
116
330
      Trace("nl-ext-cm") << "- Term: " << tf << std::endl;
117
330
      bool success = true;
118
      // tf is Figure 3 : tf( x )
119
660
      std::pair<Node, Node> bounds;
120
330
      if (tfs.first == Kind::PI)
121
      {
122
59
        bounds = {d_tstate.d_pi_bound[0], d_tstate.d_pi_bound[1]};
123
      }
124
      else
125
      {
126
542
        bounds = d_tstate.d_taylor.getTfModelBounds(
127
271
            tf, d_taylor_degree, d_tstate.d_model);
128
271
        if (bounds.first != bounds.second)
129
        {
130
271
          d_tstate.d_model.setUsedApproximate();
131
        }
132
      }
133
330
      if (!bounds.first.isNull() && !bounds.second.isNull())
134
      {
135
        // for each function in the congruence classe
136
660
        for (const Node& ctf : d_tstate.d_funcCongClass[tf])
137
        {
138
          // each term in congruence classes should be master terms
139
330
          Assert(d_tstate.d_trSlaves.find(ctf) != d_tstate.d_trSlaves.end());
140
          // we set the bounds for each slave of tf
141
718
          for (const Node& stf : d_tstate.d_trSlaves[ctf])
142
          {
143
776
            Trace("nl-ext-cm")
144
388
                << "...bound for " << stf << " : [" << bounds.first << ", "
145
388
                << bounds.second << "]" << std::endl;
146
388
            success = d_tstate.d_model.addCheckModelBound(
147
                stf, bounds.first, bounds.second);
148
          }
149
        }
150
      }
151
      else
152
      {
153
        Trace("nl-ext-cm") << "...no bound for " << tf << std::endl;
154
      }
155
330
      if (!success)
156
      {
157
        // a bound was conflicting
158
        Trace("nl-ext-cm") << "...failed to set bound for " << tf << std::endl;
159
        Trace("nl-ext-cm") << "-----" << std::endl;
160
        return false;
161
      }
162
    }
163
  }
164
  // replace the assertions
165
222
  assertions = std::move(passertions);
166
222
  return true;
167
}
168
169
33
void TranscendentalSolver::incrementTaylorDegree() { d_taylor_degree++; }
170
255
unsigned TranscendentalSolver::getTaylorDegree() const
171
{
172
255
  return d_taylor_degree;
173
}
174
175
void TranscendentalSolver::processSideEffect(const NlLemma& se)
176
{
177
  for (const std::tuple<Node, unsigned, Node>& sp : se.d_secantPoint)
178
  {
179
    Node tf = std::get<0>(sp);
180
    unsigned d = std::get<1>(sp);
181
    Node c = std::get<2>(sp);
182
    d_tstate.d_secant_points[tf][d].push_back(c);
183
  }
184
}
185
186
2576
void TranscendentalSolver::checkTranscendentalInitialRefine()
187
{
188
2576
  d_expSlv.checkInitialRefine();
189
2576
  d_sineSlv.checkInitialRefine();
190
2576
}
191
192
1574
void TranscendentalSolver::checkTranscendentalMonotonic()
193
{
194
1574
  d_expSlv.checkMonotonic();
195
1574
  d_sineSlv.checkMonotonic();
196
1574
}
197
198
222
void TranscendentalSolver::checkTranscendentalTangentPlanes()
199
{
200
222
  if (Trace.isOn("nl-ext"))
201
  {
202
    if (!d_tstate.d_funcMap.empty())
203
    {
204
      Trace("nl-ext")
205
          << "Get tangent plane lemmas for transcendental functions..."
206
          << std::endl;
207
    }
208
  }
209
  // this implements Figure 3 of "Satisfiaility Modulo Transcendental Functions
210
  // via Incremental Linearization" by Cimatti et al
211
202
  for (const std::pair<const Kind, std::vector<Node> >& tfs :
212
222
       d_tstate.d_funcMap)
213
  {
214
202
    Kind k = tfs.first;
215
202
    if (k == PI)
216
    {
217
      // We do not use Taylor approximation for PI currently.
218
      // This is because the convergence is extremely slow, and hence an
219
      // initial approximation is superior.
220
59
      continue;
221
    }
222
223
    // we substitute into the Taylor sum P_{n,f(0)}( x )
224
225
414
    for (const Node& tf : tfs.second)
226
    {
227
      // tf is Figure 3 : tf( x )
228
271
      Trace("nl-ext-tftp") << "Compute tangent planes " << tf << std::endl;
229
      // go until max degree is reached, or we don't meet bound criteria
230
1199
      for (unsigned d = 1; d <= d_taylor_degree; d++)
231
      {
232
1131
        Trace("nl-ext-tftp") << "- run at degree " << d << "..." << std::endl;
233
        unsigned prev =
234
1131
            d_tstate.d_im.numPendingLemmas() + d_tstate.d_im.numWaitingLemmas();
235
1131
        if (checkTfTangentPlanesFun(tf, d))
236
        {
237
406
          Trace("nl-ext-tftp") << "...fail, #lemmas = "
238
203
                               << (d_tstate.d_im.numPendingLemmas()
239
406
                                   + d_tstate.d_im.numWaitingLemmas() - prev)
240
203
                               << std::endl;
241
203
          break;
242
        }
243
        else
244
        {
245
928
          Trace("nl-ext-tftp") << "...success" << std::endl;
246
        }
247
      }
248
    }
249
  }
250
222
}
251
252
1131
bool TranscendentalSolver::checkTfTangentPlanesFun(Node tf, unsigned d)
253
{
254
1131
  NodeManager* nm = NodeManager::currentNM();
255
1131
  Kind k = tf.getKind();
256
  // this should only be run on master applications
257
1131
  Assert(d_tstate.d_trSlaves.find(tf) != d_tstate.d_trSlaves.end());
258
259
  // Figure 3 : c
260
2262
  Node c = d_tstate.d_model.computeAbstractModelValue(tf[0]);
261
1131
  int csign = c.getConst<Rational>().sgn();
262
1131
  if (csign == 0)
263
  {
264
    // no secant/tangent plane is necessary
265
    return true;
266
  }
267
1131
  Assert(csign == 1 || csign == -1);
268
269
  // Figure 3: P_l, P_u
270
  // mapped to for signs of c
271
2262
  std::map<int, Node> poly_approx_bounds[2];
272
2262
  TaylorGenerator::ApproximationBounds pbounds;
273
  std::uint64_t actual_d =
274
1131
      d_tstate.d_taylor.getPolynomialApproximationBoundForArg(k, c, d, pbounds);
275
1131
  poly_approx_bounds[0][1] = pbounds.d_lower;
276
1131
  poly_approx_bounds[0][-1] = pbounds.d_lower;
277
1131
  poly_approx_bounds[1][1] = pbounds.d_upperPos;
278
1131
  poly_approx_bounds[1][-1] = pbounds.d_upperNeg;
279
280
  // Figure 3 : v
281
2262
  Node v = d_tstate.d_model.computeAbstractModelValue(tf);
282
283
  // check value of tf
284
2262
  Trace("nl-ext-tftp-debug") << "Process tangent plane refinement for " << tf
285
1131
                             << ", degree " << d << "..." << std::endl;
286
1131
  Trace("nl-ext-tftp-debug") << "  value in model : " << v << std::endl;
287
1131
  Trace("nl-ext-tftp-debug") << "  arg value in model : " << c << std::endl;
288
289
  // compute the concavity
290
1131
  int region = -1;
291
  std::unordered_map<Node, int, NodeHashFunction>::iterator itr =
292
1131
      d_tstate.d_tf_region.find(tf);
293
1131
  if (itr != d_tstate.d_tf_region.end())
294
  {
295
1131
    region = itr->second;
296
1131
    Trace("nl-ext-tftp-debug") << "  region is : " << region << std::endl;
297
  }
298
  // Figure 3 : conc
299
1131
  int concavity = regionToConcavity(k, itr->second);
300
1131
  Trace("nl-ext-tftp-debug") << "  concavity is : " << concavity << std::endl;
301
1131
  if (concavity == 0)
302
  {
303
    // no secant/tangent plane is necessary
304
    return true;
305
  }
306
307
  // Figure 3: P
308
2262
  Node poly_approx;
309
310
  // compute whether this is a tangent refinement or a secant refinement
311
1131
  bool is_tangent = false;
312
1131
  bool is_secant = false;
313
  std::pair<Node, Node> mvb =
314
2262
      d_tstate.d_taylor.getTfModelBounds(tf, d, d_tstate.d_model);
315
  // this is the approximated value of tf(c), which is a value such that:
316
  //    M_A(tf(c)) >= poly_appox_c >= tf(c) or
317
  //    M_A(tf(c)) <= poly_appox_c <= tf(c)
318
  // In other words, it is a better approximation of the true value of tf(c)
319
  // in the case that we add a refinement lemma. We use this value in the
320
  // refinement schemas below.
321
2262
  Node poly_approx_c;
322
3082
  for (unsigned r = 0; r < 2; r++)
323
  {
324
4105
    Node pab = poly_approx_bounds[r][csign];
325
4105
    Node v_pab = r == 0 ? mvb.first : mvb.second;
326
2154
    if (!v_pab.isNull())
327
    {
328
4308
      Trace("nl-trans") << "...model value of " << pab << " is " << v_pab
329
2154
                        << std::endl;
330
331
2154
      Assert(v_pab.isConst());
332
4105
      Node comp = nm->mkNode(r == 0 ? LT : GT, v, v_pab);
333
2154
      Trace("nl-trans") << "...compare : " << comp << std::endl;
334
4105
      Node compr = Rewriter::rewrite(comp);
335
2154
      Trace("nl-trans") << "...got : " << compr << std::endl;
336
2154
      if (compr == d_tstate.d_true)
337
      {
338
203
        poly_approx_c = Rewriter::rewrite(v_pab);
339
        // beyond the bounds
340
203
        if (r == 0)
341
        {
342
108
          poly_approx = poly_approx_bounds[r][csign];
343
108
          is_tangent = concavity == 1;
344
108
          is_secant = concavity == -1;
345
        }
346
        else
347
        {
348
95
          poly_approx = poly_approx_bounds[r][csign];
349
95
          is_tangent = concavity == -1;
350
95
          is_secant = concavity == 1;
351
        }
352
203
        if (Trace.isOn("nl-ext-tftp"))
353
        {
354
          Trace("nl-ext-tftp") << "*** Outside boundary point (";
355
          Trace("nl-ext-tftp") << (r == 0 ? "low" : "high") << ") ";
356
          printRationalApprox("nl-ext-tftp", v_pab);
357
          Trace("nl-ext-tftp") << ", will refine..." << std::endl;
358
          Trace("nl-ext-tftp")
359
              << "    poly_approx = " << poly_approx << std::endl;
360
          Trace("nl-ext-tftp")
361
              << "    is_tangent = " << is_tangent << std::endl;
362
          Trace("nl-ext-tftp") << "    is_secant = " << is_secant << std::endl;
363
        }
364
203
        break;
365
      }
366
      else
367
      {
368
3902
        Trace("nl-ext-tftp")
369
1951
            << "  ...within " << (r == 0 ? "low" : "high") << " bound : ";
370
1951
        printRationalApprox("nl-ext-tftp", v_pab);
371
1951
        Trace("nl-ext-tftp") << std::endl;
372
      }
373
    }
374
  }
375
376
  // Figure 3: P( c )
377
1131
  if (is_tangent || is_secant)
378
  {
379
406
    Trace("nl-trans") << "...poly approximation at c is " << poly_approx_c
380
406
                      << std::endl;
381
  }
382
  else
383
  {
384
    // we may want to continue getting better bounds
385
928
    return false;
386
  }
387
388
203
  if (is_tangent)
389
  {
390
88
    if (k == Kind::EXPONENTIAL)
391
    {
392
56
      d_expSlv.doTangentLemma(tf, c, poly_approx_c, d);
393
    }
394
32
    else if (k == Kind::SINE)
395
    {
396
32
      d_sineSlv.doTangentLemma(tf, c, poly_approx_c, region, d);
397
    }
398
  }
399
115
  else if (is_secant)
400
  {
401
115
    if (k == EXPONENTIAL)
402
    {
403
43
      d_expSlv.doSecantLemmas(tf, poly_approx, c, poly_approx_c, d, actual_d);
404
    }
405
72
    else if (k == Kind::SINE)
406
    {
407
72
      d_sineSlv.doSecantLemmas(
408
          tf, poly_approx, c, poly_approx_c, d, actual_d, region);
409
    }
410
  }
411
203
  return true;
412
}
413
414
1131
int TranscendentalSolver::regionToConcavity(Kind k, int region)
415
{
416
1131
  if (k == EXPONENTIAL)
417
  {
418
749
    if (region == 1)
419
    {
420
749
      return 1;
421
    }
422
  }
423
382
  else if (k == SINE)
424
  {
425
382
    if (region == 1 || region == 2)
426
    {
427
191
      return -1;
428
    }
429
191
    else if (region == 3 || region == 4)
430
    {
431
191
      return 1;
432
    }
433
  }
434
  return 0;
435
}
436
437
}  // namespace transcendental
438
}  // namespace nl
439
}  // namespace arith
440
}  // namespace theory
441
31065
}  // namespace CVC4