GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/nl/cad/cdcac_utils.cpp Lines: 133 166 80.1 %
Date: 2021-09-07 Branches: 192 456 42.1 %

Line Exec Source
1
/******************************************************************************
2
 * Top contributors (to current version):
3
 *   Gereon Kremer
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
 * Implements utilities for cdcac.
14
 */
15
16
#include "theory/arith/nl/cad/cdcac_utils.h"
17
18
#ifdef CVC5_POLY_IMP
19
20
#include "theory/arith/nl/cad/projections.h"
21
22
namespace cvc5 {
23
namespace theory {
24
namespace arith {
25
namespace nl {
26
namespace cad {
27
28
using namespace poly;
29
30
9
bool operator==(const CACInterval& lhs, const CACInterval& rhs)
31
{
32
9
  return lhs.d_interval == rhs.d_interval;
33
}
34
bool operator<(const CACInterval& lhs, const CACInterval& rhs)
35
{
36
  return lhs.d_interval < rhs.d_interval;
37
}
38
39
/**
40
 * Induces an ordering on poly intervals that is suitable for redundancy
41
 * removal as implemented in clean_intervals.
42
 */
43
5942
inline bool compareForCleanup(const Interval& lhs, const Interval& rhs)
44
{
45
5942
  const lp_value_t* ll = &(lhs.get_internal()->a);
46
  const lp_value_t* lu =
47
5942
      lhs.get_internal()->is_point ? ll : &(lhs.get_internal()->b);
48
5942
  const lp_value_t* rl = &(rhs.get_internal()->a);
49
  const lp_value_t* ru =
50
5942
      rhs.get_internal()->is_point ? rl : &(rhs.get_internal()->b);
51
52
5942
  int lc = lp_value_cmp(ll, rl);
53
  // Lower bound is smaller
54
5942
  if (lc < 0) return true;
55
  // Lower bound is larger
56
2970
  if (lc > 0) return false;
57
  // Lower bound type is smaller
58
775
  if (!lhs.get_internal()->a_open && rhs.get_internal()->a_open) return true;
59
  // Lower bound type is larger
60
771
  if (lhs.get_internal()->a_open && !rhs.get_internal()->a_open) return false;
61
62
  // Attention: Here it differs from the regular interval ordering!
63
759
  int uc = lp_value_cmp(lu, ru);
64
  // Upper bound is smaller
65
759
  if (uc < 0) return false;
66
  // Upper bound is larger
67
549
  if (uc > 0) return true;
68
  // Upper bound type is smaller
69
405
  if (lhs.get_internal()->b_open && !rhs.get_internal()->b_open) return false;
70
  // Upper bound type is larger
71
405
  if (!lhs.get_internal()->b_open && rhs.get_internal()->b_open) return true;
72
73
  // Identical
74
405
  return false;
75
}
76
77
1481
bool intervalCovers(const Interval& lhs, const Interval& rhs)
78
{
79
1481
  const lp_value_t* ll = &(lhs.get_internal()->a);
80
  const lp_value_t* lu =
81
1481
      lhs.get_internal()->is_point ? ll : &(lhs.get_internal()->b);
82
1481
  const lp_value_t* rl = &(rhs.get_internal()->a);
83
  const lp_value_t* ru =
84
1481
      rhs.get_internal()->is_point ? rl : &(rhs.get_internal()->b);
85
86
1481
  int lc = lp_value_cmp(ll, rl);
87
1481
  int uc = lp_value_cmp(lu, ru);
88
89
  // Lower bound is smaller and upper bound is larger
90
1481
  if (lc < 0 && uc > 0) return true;
91
  // Lower bound is larger or upper bound is smaller
92
748
  if (lc > 0 || uc < 0) return false;
93
94
  // Now both bounds are identical.
95
463
  Assert(lc <= 0 && uc >= 0);
96
463
  Assert(lc == 0 || uc == 0);
97
98
  // Lower bound is the same and the bound type is stricter
99
463
  if (lc == 0 && lhs.get_internal()->a_open && !rhs.get_internal()->a_open)
100
    return false;
101
  // Upper bound is the same and the bound type is stricter
102
463
  if (uc == 0 && lhs.get_internal()->b_open && !rhs.get_internal()->b_open)
103
20
    return false;
104
105
  // Both bounds are weaker
106
443
  return true;
107
}
108
109
281
bool intervalConnect(const Interval& lhs, const Interval& rhs)
110
{
111
281
  Assert(lhs < rhs) << "Can only check for a connection if lhs < rhs.";
112
281
  const lp_value_t* lu = lhs.get_internal()->is_point
113
556
                             ? &(lhs.get_internal()->a)
114
837
                             : &(lhs.get_internal()->b);
115
281
  const lp_value_t* rl = &(rhs.get_internal()->a);
116
281
  int c = lp_value_cmp(lu, rl);
117
281
  if (c < 0)
118
  {
119
60
    Trace("libpoly::interval_connect")
120
30
        << lhs << " and " << rhs << " are separated." << std::endl;
121
30
    return false;
122
  }
123
251
  if (c > 0)
124
  {
125
262
    Trace("libpoly::interval_connect")
126
131
        << lhs << " and " << rhs << " overlap." << std::endl;
127
131
    return true;
128
  }
129
120
  Assert(c == 0);
130
354
  if (lhs.get_internal()->is_point || rhs.get_internal()->is_point
131
218
      || !lhs.get_internal()->b_open || !rhs.get_internal()->a_open)
132
  {
133
76
    Trace("libpoly::interval_connect")
134
38
        << lhs << " and " << rhs
135
38
        << " touch and the intermediate point is covered." << std::endl;
136
38
    return true;
137
  }
138
82
  return false;
139
}
140
141
321
void cleanIntervals(std::vector<CACInterval>& intervals)
142
{
143
  // Simplifies removal of redundancies later on.
144
321
  if (intervals.size() < 2) return;
145
146
  // Sort intervals.
147
264
  std::sort(intervals.begin(),
148
            intervals.end(),
149
5942
            [](const CACInterval& lhs, const CACInterval& rhs) {
150
5942
              return compareForCleanup(lhs.d_interval, rhs.d_interval);
151
5942
            });
152
153
  // Remove intervals that are covered by others.
154
  // Implementation roughly follows
155
  // https://en.cppreference.com/w/cpp/algorithm/remove Find first interval that
156
  // covers the next one.
157
264
  std::size_t first = 0;
158
371
  for (std::size_t n = intervals.size(); first < n - 1; ++first)
159
  {
160
310
    if (intervalCovers(intervals[first].d_interval,
161
310
                       intervals[first + 1].d_interval))
162
    {
163
203
      break;
164
    }
165
  }
166
  // If such an interval exists, remove accordingly.
167
264
  if (first < intervals.size() - 1)
168
  {
169
1374
    for (std::size_t i = first + 2, n = intervals.size(); i < n; ++i)
170
    {
171
1171
      if (!intervalCovers(intervals[first].d_interval, intervals[i].d_interval))
172
      {
173
        // Interval is not covered. Move it to the front and bump front.
174
198
        ++first;
175
198
        intervals[first] = std::move(intervals[i]);
176
      }
177
      // Else: Interval is covered as well.
178
    }
179
    // Erase trailing values
180
1176
    while (intervals.size() > first + 1)
181
    {
182
1176
      intervals.pop_back();
183
    }
184
  }
185
}
186
187
124
std::vector<Node> collectConstraints(const std::vector<CACInterval>& intervals)
188
{
189
124
  std::vector<Node> res;
190
372
  for (const auto& i : intervals)
191
  {
192
248
    res.insert(res.end(), i.d_origins.begin(), i.d_origins.end());
193
  }
194
124
  std::sort(res.begin(), res.end());
195
124
  auto it = std::unique(res.begin(), res.end());
196
124
  res.erase(it, res.end());
197
124
  return res;
198
}
199
200
321
bool sampleOutside(const std::vector<CACInterval>& infeasible, Value& sample)
201
{
202
321
  if (infeasible.empty())
203
  {
204
    // No infeasible region, just take anything: zero
205
20
    sample = poly::Integer();
206
20
    return true;
207
  }
208
301
  if (!is_minus_infinity(get_lower(infeasible.front().d_interval)))
209
  {
210
    // First does not cover -oo, just take sufficiently low value
211
60
    Trace("cdcac") << "Sample before " << infeasible.front().d_interval
212
30
                   << std::endl;
213
30
    const auto* i = infeasible.front().d_interval.get_internal();
214
60
    sample = value_between(
215
90
        Value::minus_infty().get_internal(), true, &i->a, !i->a_open);
216
30
    return true;
217
  }
218
440
  for (std::size_t i = 0, n = infeasible.size(); i < n - 1; ++i)
219
  {
220
    // Search for two subsequent intervals that do not connect
221
281
    if (!intervalConnect(infeasible[i].d_interval,
222
281
                         infeasible[i + 1].d_interval))
223
    {
224
      // Two intervals do not connect, take something from the gap
225
112
      const auto* l = infeasible[i].d_interval.get_internal();
226
112
      const auto* r = infeasible[i + 1].d_interval.get_internal();
227
228
224
      Trace("cdcac") << "Sample between " << infeasible[i].d_interval << " and "
229
112
                     << infeasible[i + 1].d_interval << std::endl;
230
231
112
      if (l->is_point)
232
      {
233
        sample = value_between(&l->a, true, &r->a, !r->a_open);
234
      }
235
      else
236
      {
237
112
        sample = value_between(&l->b, !l->b_open, &r->a, !r->a_open);
238
      }
239
112
      return true;
240
    }
241
    else
242
    {
243
338
      Trace("cdcac") << infeasible[i].d_interval << " and "
244
169
                     << infeasible[i + 1].d_interval << " connect" << std::endl;
245
    }
246
  }
247
159
  if (!is_plus_infinity(get_upper(infeasible.back().d_interval)))
248
  {
249
    // Last does not cover oo, just take something sufficiently large
250
70
    Trace("cdcac") << "Sample above " << infeasible.back().d_interval
251
35
                   << std::endl;
252
35
    const auto* i = infeasible.back().d_interval.get_internal();
253
35
    if (i->is_point)
254
    {
255
10
      sample =
256
20
          value_between(&i->a, true, Value::plus_infty().get_internal(), true);
257
    }
258
    else
259
    {
260
25
      sample = value_between(
261
50
          &i->b, !i->b_open, Value::plus_infty().get_internal(), true);
262
    }
263
35
    return true;
264
  }
265
124
  return false;
266
}
267
268
namespace {
269
/**
270
 * Replace a polynomial at polys[id] with the given pair of polynomials.
271
 * Also update d_mainPolys in the given interval accordingly.
272
 * If one of the factors (from the pair) is from a lower level (has a different
273
 * main variable), push this factor to the d_downPolys.
274
 * The first factor needs to be a proper polynomial (!is_constant(subst.first)),
275
 * but the second factor may be anything.
276
 */
277
void replace_polynomial(PolyVector& polys,
278
                        std::size_t id,
279
                        std::pair<poly::Polynomial, poly::Polynomial> subst,
280
                        CACInterval& interval)
281
{
282
  Assert(polys[id] == subst.first * subst.second);
283
  Assert(!is_constant(subst.first));
284
  // Whether polys[id] has already been replaced
285
  bool replaced = false;
286
  poly::Variable var = main_variable(polys[id]);
287
  // The position within interval.d_mainPolys to be replaced.
288
  auto mit = std::find(
289
      interval.d_mainPolys.begin(), interval.d_mainPolys.end(), polys[id]);
290
  if (main_variable(subst.first) == var)
291
  {
292
    // Replace in polys[id] and *mit
293
    polys[id] = subst.first;
294
    if (mit != interval.d_mainPolys.end())
295
    {
296
      *mit = subst.first;
297
    }
298
    replaced = true;
299
  }
300
  else
301
  {
302
    // Push to d_downPolys
303
    interval.d_downPolys.add(subst.first);
304
  }
305
  // Skip constant poly
306
  if (!is_constant(subst.second))
307
  {
308
    if (main_variable(subst.second) == var)
309
    {
310
      if (replaced)
311
      {
312
        // Append to polys and d_mainPolys
313
        polys.add(subst.second);
314
        interval.d_mainPolys.add(subst.second);
315
      }
316
      else
317
      {
318
        // Replace in polys[id] and *mit
319
        polys[id] = subst.second;
320
        if (mit != interval.d_mainPolys.end())
321
        {
322
          *mit = subst.second;
323
        }
324
        replaced = true;
325
      }
326
    }
327
    else
328
    {
329
      // Push to d_downPolys
330
      interval.d_downPolys.add(subst.second);
331
    }
332
  }
333
  Assert(replaced)
334
      << "At least one of the factors should have the main variable";
335
}
336
}  // namespace
337
338
82
void makeFinestSquareFreeBasis(CACInterval& lhs, CACInterval& rhs)
339
{
340
82
  auto& l = lhs.d_upperPolys;
341
82
  auto& r = rhs.d_lowerPolys;
342
82
  if (l.empty()) return;
343
164
  for (std::size_t i = 0, ln = l.size(); i < ln; ++i)
344
  {
345
164
    for (std::size_t j = 0, rn = r.size(); j < rn; ++j)
346
    {
347
      // All main variables should be the same
348
82
      Assert(main_variable(l[i]) == main_variable(r[j]));
349
82
      if (l[i] == r[j]) continue;
350
148
      Polynomial g = gcd(l[i], r[j]);
351
74
      if (!is_constant(g))
352
      {
353
        auto newl = div(l[i], g);
354
        auto newr = div(r[j], g);
355
        replace_polynomial(l, i, std::make_pair(g, newl), lhs);
356
        replace_polynomial(r, j, std::make_pair(g, newr), rhs);
357
      }
358
    }
359
  }
360
82
  l.reduce();
361
82
  r.reduce();
362
82
  lhs.d_mainPolys.reduce();
363
82
  rhs.d_mainPolys.reduce();
364
82
  lhs.d_downPolys.reduce();
365
82
  rhs.d_downPolys.reduce();
366
}
367
368
}  // namespace cad
369
}  // namespace nl
370
}  // namespace arith
371
}  // namespace theory
372
29502
}  // namespace cvc5
373
374
#endif