GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/nl/cad/cdcac_utils.cpp Lines: 131 166 78.9 %
Date: 2021-05-22 Branches: 187 458 40.8 %

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
4194
inline bool compareForCleanup(const Interval& lhs, const Interval& rhs)
44
{
45
4194
  const lp_value_t* ll = &(lhs.get_internal()->a);
46
  const lp_value_t* lu =
47
4194
      lhs.get_internal()->is_point ? ll : &(lhs.get_internal()->b);
48
4194
  const lp_value_t* rl = &(rhs.get_internal()->a);
49
  const lp_value_t* ru =
50
4194
      rhs.get_internal()->is_point ? rl : &(rhs.get_internal()->b);
51
52
4194
  int lc = lp_value_cmp(ll, rl);
53
  // Lower bound is smaller
54
4194
  if (lc < 0) return true;
55
  // Lower bound is larger
56
1974
  if (lc > 0) return false;
57
  // Lower bound type is smaller
58
535
  if (!lhs.get_internal()->a_open && rhs.get_internal()->a_open) return true;
59
  // Lower bound type is larger
60
535
  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
525
  int uc = lp_value_cmp(lu, ru);
64
  // Upper bound is smaller
65
525
  if (uc < 0) return false;
66
  // Upper bound is larger
67
411
  if (uc > 0) return true;
68
  // Upper bound type is smaller
69
297
  if (lhs.get_internal()->b_open && !rhs.get_internal()->b_open) return false;
70
  // Upper bound type is larger
71
297
  if (!lhs.get_internal()->b_open && rhs.get_internal()->b_open) return true;
72
73
  // Identical
74
297
  return false;
75
}
76
77
970
bool intervalCovers(const Interval& lhs, const Interval& rhs)
78
{
79
970
  const lp_value_t* ll = &(lhs.get_internal()->a);
80
  const lp_value_t* lu =
81
970
      lhs.get_internal()->is_point ? ll : &(lhs.get_internal()->b);
82
970
  const lp_value_t* rl = &(rhs.get_internal()->a);
83
  const lp_value_t* ru =
84
970
      rhs.get_internal()->is_point ? rl : &(rhs.get_internal()->b);
85
86
970
  int lc = lp_value_cmp(ll, rl);
87
970
  int uc = lp_value_cmp(lu, ru);
88
89
  // Lower bound is smaller and upper bound is larger
90
970
  if (lc < 0 && uc > 0) return true;
91
  // Lower bound is larger or upper bound is smaller
92
431
  if (lc > 0 || uc < 0) return false;
93
94
  // Now both bounds are identical.
95
271
  Assert(lc <= 0 && uc >= 0);
96
271
  Assert(lc == 0 || uc == 0);
97
98
  // Lower bound is the same and the bound type is stricter
99
271
  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
271
  if (uc == 0 && lhs.get_internal()->b_open && !rhs.get_internal()->b_open)
103
10
    return false;
104
105
  // Both bounds are weaker
106
261
  return true;
107
}
108
109
148
bool intervalConnect(const Interval& lhs, const Interval& rhs)
110
{
111
148
  Assert(lhs < rhs) << "Can only check for a connection if lhs < rhs.";
112
148
  const lp_value_t* lu = lhs.get_internal()->is_point
113
294
                             ? &(lhs.get_internal()->a)
114
442
                             : &(lhs.get_internal()->b);
115
148
  const lp_value_t* rl = &(rhs.get_internal()->a);
116
148
  int c = lp_value_cmp(lu, rl);
117
148
  if (c < 0)
118
  {
119
40
    Trace("libpoly::interval_connect")
120
20
        << lhs << " and " << rhs << " are separated." << std::endl;
121
20
    return false;
122
  }
123
128
  if (c > 0)
124
  {
125
110
    Trace("libpoly::interval_connect")
126
55
        << lhs << " and " << rhs << " overlap." << std::endl;
127
55
    return true;
128
  }
129
73
  Assert(c == 0);
130
217
  if (lhs.get_internal()->is_point || rhs.get_internal()->is_point
131
138
      || !lhs.get_internal()->b_open || !rhs.get_internal()->a_open)
132
  {
133
48
    Trace("libpoly::interval_connect")
134
24
        << lhs << " and " << rhs
135
24
        << " touch and the intermediate point is covered." << std::endl;
136
24
    return true;
137
  }
138
49
  return false;
139
}
140
141
135
void cleanIntervals(std::vector<CACInterval>& intervals)
142
{
143
  // Simplifies removal of redundancies later on.
144
135
  if (intervals.size() < 2) return;
145
146
  // Sort intervals.
147
125
  std::sort(intervals.begin(),
148
            intervals.end(),
149
4194
            [](const CACInterval& lhs, const CACInterval& rhs) {
150
4194
              return compareForCleanup(lhs.d_interval, rhs.d_interval);
151
4194
            });
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
125
  std::size_t first = 0;
158
169
  for (std::size_t n = intervals.size(); first < n - 1; ++first)
159
  {
160
147
    if (intervalCovers(intervals[first].d_interval,
161
147
                       intervals[first + 1].d_interval))
162
    {
163
103
      break;
164
    }
165
  }
166
  // If such an interval exists, remove accordingly.
167
125
  if (first < intervals.size() - 1)
168
  {
169
926
    for (std::size_t i = first + 2, n = intervals.size(); i < n; ++i)
170
    {
171
823
      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
126
        ++first;
175
126
        intervals[first] = std::move(intervals[i]);
176
      }
177
      // Else: Interval is covered as well.
178
    }
179
    // Erase trailing values
180
800
    while (intervals.size() > first + 1)
181
    {
182
800
      intervals.pop_back();
183
    }
184
  }
185
}
186
187
42
std::vector<Node> collectConstraints(const std::vector<CACInterval>& intervals)
188
{
189
42
  std::vector<Node> res;
190
134
  for (const auto& i : intervals)
191
  {
192
92
    res.insert(res.end(), i.d_origins.begin(), i.d_origins.end());
193
  }
194
42
  std::sort(res.begin(), res.end());
195
42
  auto it = std::unique(res.begin(), res.end());
196
42
  res.erase(it, res.end());
197
42
  return res;
198
}
199
200
135
bool sampleOutside(const std::vector<CACInterval>& infeasible, Value& sample)
201
{
202
135
  if (infeasible.empty())
203
  {
204
    // No infeasible region, just take anything: zero
205
    sample = poly::Integer();
206
    return true;
207
  }
208
135
  if (!is_minus_infinity(get_lower(infeasible.front().d_interval)))
209
  {
210
    // First does not cover -oo, just take sufficiently low value
211
12
    Trace("cdcac") << "Sample before " << infeasible.front().d_interval
212
6
                   << std::endl;
213
6
    const auto* i = infeasible.front().d_interval.get_internal();
214
12
    sample = value_between(
215
18
        Value::minus_infty().get_internal(), true, &i->a, !i->a_open);
216
6
    return true;
217
  }
218
208
  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
148
    if (!intervalConnect(infeasible[i].d_interval,
222
148
                         infeasible[i + 1].d_interval))
223
    {
224
      // Two intervals do not connect, take something from the gap
225
69
      const auto* l = infeasible[i].d_interval.get_internal();
226
69
      const auto* r = infeasible[i + 1].d_interval.get_internal();
227
228
138
      Trace("cdcac") << "Sample between " << infeasible[i].d_interval << " and "
229
69
                     << infeasible[i + 1].d_interval << std::endl;
230
231
69
      if (l->is_point)
232
      {
233
        sample = value_between(&l->a, true, &r->a, !r->a_open);
234
      }
235
      else
236
      {
237
69
        sample = value_between(&l->b, !l->b_open, &r->a, !r->a_open);
238
      }
239
69
      return true;
240
    }
241
    else
242
    {
243
158
      Trace("cdcac") << infeasible[i].d_interval << " and "
244
79
                     << infeasible[i + 1].d_interval << " connect" << std::endl;
245
    }
246
  }
247
60
  if (!is_plus_infinity(get_upper(infeasible.back().d_interval)))
248
  {
249
    // Last does not cover oo, just take something sufficiently large
250
36
    Trace("cdcac") << "Sample above " << infeasible.back().d_interval
251
18
                   << std::endl;
252
18
    const auto* i = infeasible.back().d_interval.get_internal();
253
18
    if (i->is_point)
254
    {
255
4
      sample =
256
8
          value_between(&i->a, true, Value::plus_infty().get_internal(), true);
257
    }
258
    else
259
    {
260
14
      sample = value_between(
261
28
          &i->b, !i->b_open, Value::plus_infty().get_internal(), true);
262
    }
263
18
    return true;
264
  }
265
42
  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
28
void makeFinestSquareFreeBasis(CACInterval& lhs, CACInterval& rhs)
339
{
340
28
  auto& l = lhs.d_upperPolys;
341
28
  auto& r = rhs.d_lowerPolys;
342
28
  if (l.empty()) return;
343
56
  for (std::size_t i = 0, ln = l.size(); i < ln; ++i)
344
  {
345
56
    for (std::size_t j = 0, rn = r.size(); j < rn; ++j)
346
    {
347
      // All main variables should be the same
348
28
      Assert(main_variable(l[i]) == main_variable(r[j]));
349
28
      if (l[i] == r[j]) continue;
350
48
      Polynomial g = gcd(l[i], r[j]);
351
24
      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
28
  l.reduce();
361
28
  r.reduce();
362
28
  lhs.d_mainPolys.reduce();
363
28
  rhs.d_mainPolys.reduce();
364
28
  lhs.d_downPolys.reduce();
365
28
  rhs.d_downPolys.reduce();
366
}
367
368
}  // namespace cad
369
}  // namespace nl
370
}  // namespace arith
371
}  // namespace theory
372
28191
}  // namespace cvc5
373
374
#endif