GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/nl/transcendental/transcendental_solver.cpp Lines: 178 207 86.0 %
Date: 2021-09-07 Branches: 345 821 42.0 %

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