GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/linear_equality.cpp Lines: 368 804 45.8 %
Date: 2021-09-29 Branches: 620 3296 18.8 %

Line Exec Source
1
/******************************************************************************
2
 * Top contributors (to current version):
3
 *   Tim King, Mathias Preiner, Morgan Deters
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
 * This implements the LinearEqualityModule.
14
 */
15
#include "theory/arith/linear_equality.h"
16
17
#include "base/output.h"
18
#include "smt/smt_statistics_registry.h"
19
#include "theory/arith/constraint.h"
20
21
22
using namespace std;
23
24
namespace cvc5 {
25
namespace theory {
26
namespace arith {
27
28
/* Explicitly instatiate these functions. */
29
30
template ArithVar LinearEqualityModule::selectSlack<true>(ArithVar x_i, VarPreferenceFunction pf) const;
31
template ArithVar LinearEqualityModule::selectSlack<false>(ArithVar x_i, VarPreferenceFunction pf) const;
32
33
template bool LinearEqualityModule::preferWitness<true>(const UpdateInfo& a, const UpdateInfo& b) const;
34
template bool LinearEqualityModule::preferWitness<false>(const UpdateInfo& a, const UpdateInfo& b) const;
35
36
37
void Border::output(std::ostream& out) const{
38
  out << "{Border"
39
      << ", " << d_bound->getVariable()
40
      << ", " << d_bound->getValue()
41
      << ", " << d_diff
42
      << ", " << d_areFixing
43
      << ", " << d_upperbound;
44
  if(ownBorder()){
45
    out << ", ownBorder";
46
  }else{
47
    out << ", " << d_entry->getCoefficient();
48
  }
49
  out << ", " << d_bound
50
      << "}";
51
}
52
53
6271
LinearEqualityModule::LinearEqualityModule(ArithVariables& vars, Tableau& t, BoundInfoMap& boundsTracking, BasicVarModelUpdateCallBack f):
54
  d_variables(vars),
55
  d_tableau(t),
56
  d_basicVariableUpdates(f),
57
  d_increasing(1),
58
  d_decreasing(-1),
59
  d_upperBoundDifference(),
60
  d_lowerBoundDifference(),
61
  d_one(1),
62
  d_negOne(-1),
63
  d_btracking(boundsTracking),
64
  d_areTracking(false),
65
6271
  d_trackCallback(this)
66
6271
{}
67
68
6271
LinearEqualityModule::Statistics::Statistics()
69
    : d_statPivots(
70
12542
        smtStatisticsRegistry().registerInt("theory::arith::pivots")),
71
      d_statUpdates(
72
12542
          smtStatisticsRegistry().registerInt("theory::arith::updates")),
73
      d_pivotTime(
74
12542
          smtStatisticsRegistry().registerTimer("theory::arith::pivotTime")),
75
      d_adjTime(
76
12542
          smtStatisticsRegistry().registerTimer("theory::arith::adjTime")),
77
6271
      d_weakeningAttempts(smtStatisticsRegistry().registerInt(
78
12542
          "theory::arith::weakening::attempts")),
79
6271
      d_weakeningSuccesses(smtStatisticsRegistry().registerInt(
80
12542
          "theory::arith::weakening::success")),
81
6271
      d_weakenings(smtStatisticsRegistry().registerInt(
82
12542
          "theory::arith::weakening::total")),
83
6271
      d_weakenTime(smtStatisticsRegistry().registerTimer(
84
12542
          "theory::arith::weakening::time")),
85
      d_forceTime(
86
56439
          smtStatisticsRegistry().registerTimer("theory::arith::forcing::time"))
87
{
88
6271
}
89
90
2803904
void LinearEqualityModule::includeBoundUpdate(ArithVar v, const BoundsInfo& prev){
91
2803904
  Assert(!d_areTracking);
92
93
2803904
  BoundsInfo curr = d_variables.boundsInfo(v);
94
95
2803904
  Assert(prev != curr);
96
2803904
  Tableau::ColIterator basicIter = d_tableau.colIterator(v);
97
27004658
  for(; !basicIter.atEnd(); ++basicIter){
98
12100377
    const Tableau::Entry& entry = *basicIter;
99
12100377
    Assert(entry.getColVar() == v);
100
12100377
    int a_ijSgn = entry.getCoefficient().sgn();
101
102
12100377
    RowIndex ridx = entry.getRowIndex();
103
12100377
    BoundsInfo& counts = d_btracking.get(ridx);
104
12100377
    Debug("includeBoundUpdate") << d_tableau.rowIndexToBasic(ridx) << " " << counts << " to " ;
105
12100377
    counts.addInChange(a_ijSgn, prev, curr);
106
12100377
    Debug("includeBoundUpdate") << counts << " " << a_ijSgn << std::endl;
107
  }
108
2803904
}
109
110
void LinearEqualityModule::updateMany(const DenseMap<DeltaRational>& many){
111
  for(DenseMap<DeltaRational>::const_iterator i = many.begin(), i_end = many.end(); i != i_end; ++i){
112
    ArithVar nb = *i;
113
    if(!d_tableau.isBasic(nb)){
114
      Assert(!d_tableau.isBasic(nb));
115
      const DeltaRational& newValue = many[nb];
116
      if(newValue != d_variables.getAssignment(nb)){
117
        Trace("arith::updateMany")
118
          << "updateMany:" << nb << " "
119
          << d_variables.getAssignment(nb) << " to "<< newValue << endl;
120
        update(nb, newValue);
121
      }
122
    }
123
  }
124
}
125
126
127
128
129
void LinearEqualityModule::applySolution(const DenseSet& newBasis, const DenseMap<DeltaRational>& newValues){
130
  forceNewBasis(newBasis);
131
  updateMany(newValues);
132
}
133
134
void LinearEqualityModule::forceNewBasis(const DenseSet& newBasis){
135
  TimerStat::CodeTimer codeTimer(d_statistics.d_forceTime);
136
  cout << "force begin" << endl;
137
  DenseSet needsToBeAdded;
138
  for(DenseSet::const_iterator i = newBasis.begin(), i_end = newBasis.end(); i != i_end; ++i){
139
    ArithVar b = *i;
140
    if(!d_tableau.isBasic(b)){
141
      needsToBeAdded.add(b);
142
    }
143
  }
144
145
  while(!needsToBeAdded.empty()){
146
    ArithVar toRemove = ARITHVAR_SENTINEL;
147
    ArithVar toAdd = ARITHVAR_SENTINEL;
148
    DenseSet::const_iterator i = needsToBeAdded.begin(), i_end = needsToBeAdded.end();
149
    for(; toAdd == ARITHVAR_SENTINEL && i != i_end; ++i){
150
      ArithVar v = *i;
151
152
      Tableau::ColIterator colIter = d_tableau.colIterator(v);
153
      for(; !colIter.atEnd(); ++colIter){
154
        const Tableau::Entry& entry = *colIter;
155
        Assert(entry.getColVar() == v);
156
        ArithVar b = d_tableau.rowIndexToBasic(entry.getRowIndex());
157
        if(!newBasis.isMember(b)){
158
          toAdd = v;
159
          if(toRemove == ARITHVAR_SENTINEL ||
160
             d_tableau.basicRowLength(toRemove) > d_tableau.basicRowLength(b)){
161
            toRemove = b;
162
          }
163
        }
164
      }
165
    }
166
    Assert(toRemove != ARITHVAR_SENTINEL);
167
    Assert(toAdd != ARITHVAR_SENTINEL);
168
169
    Trace("arith::forceNewBasis") << toRemove << " " << toAdd << endl;
170
    CVC5Message() << toRemove << " " << toAdd << endl;
171
    d_tableau.pivot(toRemove, toAdd, d_trackCallback);
172
    d_basicVariableUpdates(toAdd);
173
174
    Trace("arith::forceNewBasis") << needsToBeAdded.size() << "to go" << endl;
175
    CVC5Message() << needsToBeAdded.size() << "to go" << endl;
176
    needsToBeAdded.remove(toAdd);
177
  }
178
}
179
180
83541
void LinearEqualityModule::updateUntracked(ArithVar x_i, const DeltaRational& v){
181
83541
  Assert(!d_tableau.isBasic(x_i));
182
83541
  Assert(!d_areTracking);
183
83541
  const DeltaRational& assignment_x_i = d_variables.getAssignment(x_i);
184
83541
  ++(d_statistics.d_statUpdates);
185
186
187
167082
  Debug("arith") <<"update " << x_i << ": "
188
83541
                 << assignment_x_i << "|-> " << v << endl;
189
167082
  DeltaRational diff = v - assignment_x_i;
190
191
83541
  Tableau::ColIterator colIter = d_tableau.colIterator(x_i);
192
1836081
  for(; !colIter.atEnd(); ++colIter){
193
876270
    const Tableau::Entry& entry = *colIter;
194
876270
    Assert(entry.getColVar() == x_i);
195
196
876270
    ArithVar x_j = d_tableau.rowIndexToBasic(entry.getRowIndex());
197
876270
    const Rational& a_ji = entry.getCoefficient();
198
199
876270
    const DeltaRational& assignment = d_variables.getAssignment(x_j);
200
1752540
    DeltaRational  nAssignment = assignment+(diff * a_ji);
201
876270
    d_variables.setAssignment(x_j, nAssignment);
202
203
876270
    d_basicVariableUpdates(x_j);
204
  }
205
206
83541
  d_variables.setAssignment(x_i, v);
207
208
83541
  if(Debug.isOn("paranoid:check_tableau")){  debugCheckTableau(); }
209
83541
}
210
211
134892
void LinearEqualityModule::updateTracked(ArithVar x_i, const DeltaRational& v){
212
269784
  TimerStat::CodeTimer codeTimer(d_statistics.d_adjTime);
213
214
134892
  Assert(!d_tableau.isBasic(x_i));
215
134892
  Assert(d_areTracking);
216
217
134892
  ++(d_statistics.d_statUpdates);
218
219
269784
  DeltaRational diff =  v - d_variables.getAssignment(x_i);
220
269784
  Debug("arith") <<"update " << x_i << ": "
221
134892
                 << d_variables.getAssignment(x_i) << "|-> " << v << endl;
222
223
224
134892
  BoundCounts before = d_variables.atBoundCounts(x_i);
225
134892
  d_variables.setAssignment(x_i, v);
226
134892
  BoundCounts after = d_variables.atBoundCounts(x_i);
227
228
134892
  bool anyChange = before != after;
229
230
134892
  Tableau::ColIterator colIter = d_tableau.colIterator(x_i);
231
4923734
  for(; !colIter.atEnd(); ++colIter){
232
2394421
    const Tableau::Entry& entry = *colIter;
233
2394421
    Assert(entry.getColVar() == x_i);
234
235
2394421
    RowIndex ridx = entry.getRowIndex();
236
2394421
    ArithVar x_j = d_tableau.rowIndexToBasic(ridx);
237
2394421
    const Rational& a_ji = entry.getCoefficient();
238
239
2394421
    const DeltaRational& assignment = d_variables.getAssignment(x_j);
240
4788842
    DeltaRational  nAssignment = assignment+(diff * a_ji);
241
2394421
    Debug("update") << x_j << " " << a_ji << assignment << " -> " << nAssignment << endl;
242
2394421
    BoundCounts xjBefore = d_variables.atBoundCounts(x_j);
243
2394421
    d_variables.setAssignment(x_j, nAssignment);
244
2394421
    BoundCounts xjAfter = d_variables.atBoundCounts(x_j);
245
246
2394421
    Assert(rowIndexIsTracked(ridx));
247
2394421
    BoundsInfo& next_bc_k = d_btracking.get(ridx);
248
2394421
    if(anyChange){
249
1864541
      next_bc_k.addInAtBoundChange(a_ji.sgn(), before, after);
250
    }
251
2394421
    if(xjBefore != xjAfter){
252
277507
      next_bc_k.addInAtBoundChange(-1, xjBefore, xjAfter);
253
    }
254
255
2394421
    d_basicVariableUpdates(x_j);
256
  }
257
258
134892
  if(Debug.isOn("paranoid:check_tableau")){  debugCheckTableau(); }
259
134892
}
260
261
134892
void LinearEqualityModule::pivotAndUpdate(ArithVar x_i, ArithVar x_j, const DeltaRational& x_i_value){
262
134892
  Assert(x_i != x_j);
263
264
269784
  TimerStat::CodeTimer codeTimer(d_statistics.d_pivotTime);
265
266
  static int instance = 0;
267
268
134892
  if(Debug.isOn("arith::tracking::pre")){
269
    ++instance;
270
    Debug("arith::tracking")  << "pre update #" << instance << endl;
271
    debugCheckTracking();
272
  }
273
274
134892
  if(Debug.isOn("arith::simplex:row")){ debugPivot(x_i, x_j); }
275
276
134892
  RowIndex ridx = d_tableau.basicToRowIndex(x_i);
277
134892
  const Tableau::Entry& entry_ij =  d_tableau.findEntry(ridx, x_j);
278
134892
  Assert(!entry_ij.blank());
279
280
134892
  const Rational& a_ij = entry_ij.getCoefficient();
281
134892
  const DeltaRational& betaX_i = d_variables.getAssignment(x_i);
282
269784
  DeltaRational theta = (x_i_value - betaX_i)/a_ij;
283
269784
  DeltaRational x_j_value = d_variables.getAssignment(x_j) + theta;
284
285
134892
  updateTracked(x_j, x_j_value);
286
287
134892
  if(Debug.isOn("arith::tracking::mid")){
288
    Debug("arith::tracking")  << "postupdate prepivot #" << instance << endl;
289
    debugCheckTracking();
290
  }
291
292
  // Pivots
293
134892
  ++(d_statistics.d_statPivots);
294
295
134892
  d_tableau.pivot(x_i, x_j, d_trackCallback);
296
297
134892
  if(Debug.isOn("arith::tracking::post")){
298
    Debug("arith::tracking")  << "postpivot #" << instance << endl;
299
    debugCheckTracking();
300
  }
301
302
134892
  d_basicVariableUpdates(x_j);
303
304
134892
  if(Debug.isOn("matrix")){
305
    d_tableau.printMatrix();
306
  }
307
134892
}
308
309
uint32_t LinearEqualityModule::updateProduct(const UpdateInfo& inf) const {
310
  uint32_t colLen = d_tableau.getColLength(inf.nonbasic());
311
  if(inf.describesPivot()){
312
    Assert(inf.leaving() != inf.nonbasic());
313
    return colLen + d_tableau.basicRowLength(inf.leaving());
314
  }else{
315
    return colLen;
316
  }
317
}
318
319
void LinearEqualityModule::debugCheckTracking(){
320
  Tableau::BasicIterator basicIter = d_tableau.beginBasic(),
321
    endIter = d_tableau.endBasic();
322
  for(; basicIter != endIter; ++basicIter){
323
    ArithVar basic = *basicIter;
324
    Debug("arith::tracking") << "arith::tracking row basic: " << basic << endl;
325
326
    for(Tableau::RowIterator iter = d_tableau.basicRowIterator(basic); !iter.atEnd() && Debug.isOn("arith::tracking"); ++iter){
327
      const Tableau::Entry& entry = *iter;
328
329
      ArithVar var = entry.getColVar();
330
      const Rational& coeff = entry.getCoefficient();
331
      DeltaRational beta = d_variables.getAssignment(var);
332
      Debug("arith::tracking") << var << " " << d_variables.boundsInfo(var)
333
                               << " " << beta << coeff;
334
      if(d_variables.hasLowerBound(var)){
335
        Debug("arith::tracking") << "(lb " << d_variables.getLowerBound(var) << ")";
336
      }
337
      if(d_variables.hasUpperBound(var)){
338
        Debug("arith::tracking") << "(up " << d_variables.getUpperBound(var) << ")";
339
      }
340
      Debug("arith::tracking") << endl;
341
    }
342
    Debug("arith::tracking") << "end row"<< endl;
343
344
    if(basicIsTracked(basic)){
345
      RowIndex ridx = d_tableau.basicToRowIndex(basic);
346
      BoundsInfo computed = computeRowBoundInfo(ridx, false);
347
      Debug("arith::tracking")
348
        << "computed " << computed
349
        << " tracking " << d_btracking[ridx] << endl;
350
      Assert(computed == d_btracking[ridx]);
351
    }
352
  }
353
}
354
355
void LinearEqualityModule::debugPivot(ArithVar x_i, ArithVar x_j){
356
  Debug("arith::pivot") << "debugPivot("<< x_i  <<"|->"<< x_j << ")" << endl;
357
358
  for(Tableau::RowIterator iter = d_tableau.basicRowIterator(x_i); !iter.atEnd(); ++iter){
359
    const Tableau::Entry& entry = *iter;
360
361
    ArithVar var = entry.getColVar();
362
    const Rational& coeff = entry.getCoefficient();
363
    DeltaRational beta = d_variables.getAssignment(var);
364
    Debug("arith::pivot") << var << beta << coeff;
365
    if(d_variables.hasLowerBound(var)){
366
      Debug("arith::pivot") << "(lb " << d_variables.getLowerBound(var) << ")";
367
    }
368
    if(d_variables.hasUpperBound(var)){
369
      Debug("arith::pivot") << "(up " << d_variables.getUpperBound(var) << ")";
370
    }
371
    Debug("arith::pivot") << endl;
372
  }
373
  Debug("arith::pivot") << "end row"<< endl;
374
}
375
376
/**
377
 * This check is quite expensive.
378
 * It should be wrapped in a Debug.isOn() guard.
379
 *   if(Debug.isOn("paranoid:check_tableau")){
380
 *      checkTableau();
381
 *   }
382
 */
383
void LinearEqualityModule::debugCheckTableau(){
384
  Tableau::BasicIterator basicIter = d_tableau.beginBasic(),
385
    endIter = d_tableau.endBasic();
386
  for(; basicIter != endIter; ++basicIter){
387
    ArithVar basic = *basicIter;
388
    DeltaRational sum;
389
    Debug("paranoid:check_tableau") << "starting row" << basic << endl;
390
    Tableau::RowIterator nonbasicIter = d_tableau.basicRowIterator(basic);
391
    for(; !nonbasicIter.atEnd(); ++nonbasicIter){
392
      const Tableau::Entry& entry = *nonbasicIter;
393
      ArithVar nonbasic = entry.getColVar();
394
      if(basic == nonbasic) continue;
395
396
      const Rational& coeff = entry.getCoefficient();
397
      DeltaRational beta = d_variables.getAssignment(nonbasic);
398
      Debug("paranoid:check_tableau") << nonbasic << beta << coeff<<endl;
399
      sum = sum + (beta*coeff);
400
    }
401
    DeltaRational shouldBe = d_variables.getAssignment(basic);
402
    Debug("paranoid:check_tableau") << "ending row" << sum
403
                                    << "," << shouldBe << endl;
404
405
    Assert(sum == shouldBe);
406
  }
407
}
408
bool LinearEqualityModule::debugEntireLinEqIsConsistent(const string& s){
409
  bool result = true;
410
  for(ArithVar var = 0, end = d_tableau.getNumColumns(); var != end; ++var){
411
    //  for(VarIter i = d_variables.begin(), end = d_variables.end(); i != end; ++i){
412
    //ArithVar var = d_arithvarNodeMap.asArithVar(*i);
413
    if(!d_variables.assignmentIsConsistent(var)){
414
      d_variables.printModel(var);
415
      Warning() << s << ":" << "Assignment is not consistent for " << var ;
416
      if(d_tableau.isBasic(var)){
417
        Warning() << " (basic)";
418
      }
419
      Warning() << endl;
420
      result = false;
421
    }
422
  }
423
  return result;
424
}
425
426
230122
DeltaRational LinearEqualityModule::computeRowBound(RowIndex ridx, bool rowUb, ArithVar skip) const {
427
230122
  DeltaRational sum(0,0);
428
1610695
  for(Tableau::RowIterator i = d_tableau.ridRowIterator(ridx); !i.atEnd(); ++i){
429
1380573
    const Tableau::Entry& entry = (*i);
430
1380573
    ArithVar v = entry.getColVar();
431
1380573
    if(v == skip){ continue; }
432
433
1199802
    const Rational& coeff =  entry.getCoefficient();
434
1199802
    bool vUb = (rowUb == (coeff.sgn() > 0));
435
436
2399604
    const DeltaRational& bound = vUb ?
437
605792
      d_variables.getUpperBound(v):
438
1793812
      d_variables.getLowerBound(v);
439
440
2399604
    DeltaRational diff = bound * coeff;
441
1199802
    sum = sum + diff;
442
  }
443
230122
  return sum;
444
}
445
446
/**
447
 * Computes the value of a basic variable using the current assignment.
448
 */
449
101820
DeltaRational LinearEqualityModule::computeRowValue(ArithVar x, bool useSafe) const{
450
101820
  Assert(d_tableau.isBasic(x));
451
101820
  DeltaRational sum(0);
452
453
621356
  for(Tableau::RowIterator i = d_tableau.basicRowIterator(x); !i.atEnd(); ++i){
454
519536
    const Tableau::Entry& entry = (*i);
455
519536
    ArithVar nonbasic = entry.getColVar();
456
519536
    if(nonbasic == x) continue;
457
417716
    const Rational& coeff = entry.getCoefficient();
458
459
417716
    const DeltaRational& assignment = d_variables.getAssignment(nonbasic, useSafe);
460
417716
    sum = sum + (assignment * coeff);
461
  }
462
101820
  return sum;
463
}
464
465
2609583
const Tableau::Entry* LinearEqualityModule::rowLacksBound(RowIndex ridx, bool rowUb, ArithVar skip){
466
2609583
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
467
23150635
  for(; !iter.atEnd(); ++iter){
468
12878363
    const Tableau::Entry& entry = *iter;
469
470
12878363
    ArithVar var = entry.getColVar();
471
12878363
    if(var == skip) { continue; }
472
473
12875671
    int sgn = entry.getCoefficient().sgn();
474
12875671
    bool selectUb = (rowUb == (sgn > 0));
475
25751342
    bool hasBound = selectUb ?
476
6430784
      d_variables.hasUpperBound(var):
477
19320558
      d_variables.hasLowerBound(var);
478
12875671
    if(!hasBound){
479
2607837
      return &entry;
480
    }
481
  }
482
1746
  return NULL;
483
}
484
485
232
void LinearEqualityModule::propagateBasicFromRow(ConstraintP c,
486
                                                 bool produceProofs)
487
{
488
232
  Assert(c != NullConstraint);
489
232
  Assert(c->isUpperBound() || c->isLowerBound());
490
232
  Assert(!c->assertedToTheTheory());
491
232
  Assert(!c->hasProof());
492
493
232
  bool upperBound = c->isUpperBound();
494
232
  ArithVar basic = c->getVariable();
495
232
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
496
497
464
  ConstraintCPVec bounds;
498
232
  RationalVectorP coeffs = produceProofs ? new RationalVector() : nullptr;
499
232
  propagateRow(bounds, ridx, upperBound, c, coeffs);
500
232
  c->impliedByFarkas(bounds, coeffs, false);
501
232
  c->tryToPropagate();
502
503
232
  if(coeffs != RationalVectorPSentinel) { delete coeffs; }
504
232
}
505
506
/* An explanation of the farkas coefficients.
507
 *
508
 * We are proving c using the other variables on the row.
509
 * The proof is in terms of the other constraints and the negation of c, ~c.
510
 *
511
 * A row has the form:
512
 *   sum a_i * x_i  = 0
513
 * or
514
 *   sx + sum r y + sum q z = 0
515
 * where r > 0 and q < 0.
516
 *
517
 * If rowUp, we are proving c
518
 *   g = sum r u_y + sum q l_z
519
 *   and c is entailed by -sx <= g
520
 * If !rowUp, we are proving c
521
 *   g = sum r l_y + sum q u_z
522
 *   and c is entailed by -sx >= g
523
 *
524
 *             | s     | c         | ~c       | u_i     | l_i
525
 *   if  rowUp | s > 0 | x >= -g/s | x < -g/s | a_i > 0 | a_i < 0
526
 *   if  rowUp | s < 0 | x <= -g/s | x > -g/s | a_i > 0 | a_i < 0
527
 *   if !rowUp | s > 0 | x <= -g/s | x > -g/s | a_i < 0 | a_i > 0
528
 *   if !rowUp | s < 0 | x >= -g/s | x < -g/s | a_i < 0 | a_i > 0
529
 *
530
 *
531
 * Thus we treat !rowUp as multiplying the row by -1 and rowUp as 1
532
 * for the entire row.
533
 */
534
80793
void LinearEqualityModule::propagateRow(ConstraintCPVec& into, RowIndex ridx, bool rowUp, ConstraintP c, RationalVectorP farkas){
535
80793
  Assert(!c->assertedToTheTheory());
536
80793
  Assert(c->canBePropagated());
537
80793
  Assert(!c->hasProof());
538
539
80793
  if(farkas != RationalVectorPSentinel){
540
1907
    Assert(farkas->empty());
541
1907
    farkas->push_back(Rational(0));
542
  }
543
544
80793
  ArithVar v = c->getVariable();
545
161586
  Debug("arith::propagateRow") << "LinearEqualityModule::propagateRow("
546
80793
                                   << ridx << ", " << rowUp << ", " << v << ") start" << endl;
547
548
80793
  const Rational& multiple = rowUp ? d_one : d_negOne;
549
550
80793
  Debug("arith::propagateRow") << "multiple: " << multiple << endl;
551
552
80793
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
553
1230585
  for(; !iter.atEnd(); ++iter){
554
574896
    const Tableau::Entry& entry = *iter;
555
574896
    ArithVar nonbasic = entry.getColVar();
556
574896
    const Rational& a_ij = entry.getCoefficient();
557
574896
    int sgn = a_ij.sgn();
558
574896
    Assert(sgn != 0);
559
574896
    bool selectUb = rowUp ? (sgn > 0) : (sgn < 0);
560
561
574896
    Assert(nonbasic != v || (rowUp && a_ij.sgn() > 0 && c->isLowerBound())
562
           || (rowUp && a_ij.sgn() < 0 && c->isUpperBound())
563
           || (!rowUp && a_ij.sgn() > 0 && c->isUpperBound())
564
           || (!rowUp && a_ij.sgn() < 0 && c->isLowerBound()));
565
566
574896
    if(Debug.isOn("arith::propagateRow")){
567
      if(nonbasic == v){
568
        Debug("arith::propagateRow") << "(target) "
569
                                     << rowUp << " "
570
                                     << a_ij.sgn() << " "
571
                                     << c->isLowerBound() << " "
572
                                     << c->isUpperBound() << endl;
573
574
        Debug("arith::propagateRow") << "(target) ";
575
      }
576
      Debug("arith::propagateRow") << "propagateRow " << a_ij << " * " << nonbasic ;
577
    }
578
579
574896
    if(nonbasic == v){
580
80793
      if(farkas != RationalVectorPSentinel){
581
1907
        Assert(farkas->front().isZero());
582
3814
        Rational multAij = multiple * a_ij;
583
1907
        Debug("arith::propagateRow") << "(" << multAij << ") ";
584
1907
        farkas->front() = multAij;
585
      }
586
587
80793
      Debug("arith::propagateRow") << c << endl;
588
    }else{
589
590
      ConstraintCP bound = selectUb
591
206134
        ? d_variables.getUpperBoundConstraint(nonbasic)
592
700237
        : d_variables.getLowerBoundConstraint(nonbasic);
593
594
494103
      if(farkas != RationalVectorPSentinel){
595
20696
        Rational multAij = multiple * a_ij;
596
10348
        Debug("arith::propagateRow") << "(" << multAij << ") ";
597
10348
        farkas->push_back(multAij);
598
      }
599
494103
      Assert(bound != NullConstraint);
600
494103
      Debug("arith::propagateRow") << bound << endl;
601
494103
      into.push_back(bound);
602
    }
603
  }
604
161586
  Debug("arith::propagateRow") << "LinearEqualityModule::propagateRow("
605
80793
                                   << ridx << ", " << rowUp << ", " << v << ") done" << endl;
606
607
80793
}
608
609
750994
ConstraintP LinearEqualityModule::weakestExplanation(bool aboveUpper, DeltaRational& surplus, ArithVar v, const Rational& coeff, bool& anyWeakening, ArithVar basic) const {
610
611
750994
  int sgn = coeff.sgn();
612
750994
  bool ub = aboveUpper?(sgn < 0) : (sgn > 0);
613
614
1501988
  ConstraintP c = ub ?
615
179851
    d_variables.getUpperBoundConstraint(v) :
616
1322137
    d_variables.getLowerBoundConstraint(v);
617
618
  bool weakened;
619
817843
  do{
620
817843
    const DeltaRational& bound = c->getValue();
621
622
817843
    weakened = false;
623
624
1635686
    ConstraintP weaker = ub?
625
225357
      c->getStrictlyWeakerUpperBound(true, true):
626
1410329
      c->getStrictlyWeakerLowerBound(true, true);
627
628
817843
    if(weaker != NullConstraint){
629
170309
      const DeltaRational& weakerBound = weaker->getValue();
630
631
340618
      DeltaRational diff = aboveUpper ? bound - weakerBound : weakerBound - bound;
632
      //if var == basic,
633
      //  if aboveUpper, weakerBound > bound, multiply by -1
634
      //  if !aboveUpper, weakerBound < bound, multiply by -1
635
170309
      diff = diff * coeff;
636
170309
      if(surplus > diff){
637
66849
        ++d_statistics.d_weakenings;
638
66849
        weakened = true;
639
66849
        anyWeakening = true;
640
66849
        surplus = surplus - diff;
641
642
66849
        Debug("arith::weak") << "found:" << endl;
643
66849
        if(v == basic){
644
7069
          Debug("arith::weak") << "  basic: ";
645
        }
646
133698
        Debug("arith::weak") << "  " << surplus << " "<< diff  << endl
647
66849
                             << "  " << bound << c << endl
648
66849
                             << "  " << weakerBound << weaker << endl;
649
650
66849
        Assert(diff.sgn() > 0);
651
66849
        c = weaker;
652
      }
653
    }
654
  }while(weakened);
655
656
750994
  return c;
657
}
658
659
/* An explanation of the farkas coefficients.
660
 *
661
 * We are proving a conflict on the basic variable x_b.
662
 * If aboveUpper, then the conflict is with the constraint c : x_b <= u_b.
663
 * If !aboveUpper, then the conflict is with the constraint c : x_b >= l_b.
664
 *
665
 * A row has the form:
666
 *   -x_b sum a_i * x_i  = 0
667
 * or
668
 *   -x_b + sum r y + sum q z = 0,
669
 *    x_b = sum r y + sum q z
670
 * where r > 0 and q < 0.
671
 *
672
 *
673
 * If !aboveUp, we are proving ~c: x_b < l_b
674
 *   g = sum r u_y + sum q l_z
675
 *   x_b <= g < l_b
676
 *   and ~c is entailed by x_b <= g
677
 *
678
 * If aboveUp, we are proving ~c : x_b > u_b
679
 *   g = sum r l_y + sum q u_z
680
 *   x_b >= g > u_b
681
 *   and ~c is entailed by x_b >= g
682
 *
683
 *
684
 *               | s     | c         | ~c       | u_i     | l_i
685
 *   if !aboveUp | s > 0 | x >= -g/s | x < -g/s | a_i > 0 | a_i < 0
686
 *   if !aboveUp | s < 0 | x <= -g/s | x > -g/s | a_i > 0 | a_i < 0
687
 *   if  aboveUp | s > 0 | x <= -g/s | x > -g/s | a_i < 0 | a_i > 0
688
 *   if  aboveUp | s < 0 | x >= -g/s | x < -g/s | a_i < 0 | a_i > 0
689
 *
690
 * Thus we treat aboveUp as multiplying the row by -1 and !aboveUp as 1
691
 * for the entire row.
692
 */
693
41940
ConstraintCP LinearEqualityModule::minimallyWeakConflict(bool aboveUpper, ArithVar basicVar, FarkasConflictBuilder& fcs) const {
694
41940
  Assert(!fcs.underConstruction());
695
83880
  TimerStat::CodeTimer codeTimer(d_statistics.d_weakenTime);
696
697
83880
  Debug("arith::weak") << "LinearEqualityModule::minimallyWeakConflict("
698
41940
                       << aboveUpper <<", "<< basicVar << ", ...) start" << endl;
699
700
41940
  const Rational& adjustSgn = aboveUpper ? d_negOne : d_one;
701
41940
  const DeltaRational& assignment = d_variables.getAssignment(basicVar);
702
83880
  DeltaRational surplus;
703
41940
  if(aboveUpper){
704
16508
    Assert(d_variables.hasUpperBound(basicVar));
705
16508
    Assert(assignment > d_variables.getUpperBound(basicVar));
706
16508
    surplus = assignment - d_variables.getUpperBound(basicVar);
707
  }else{
708
25432
    Assert(d_variables.hasLowerBound(basicVar));
709
25432
    Assert(assignment < d_variables.getLowerBound(basicVar));
710
25432
    surplus = d_variables.getLowerBound(basicVar) - assignment;
711
  }
712
713
41940
  bool anyWeakenings = false;
714
792934
  for(Tableau::RowIterator i = d_tableau.basicRowIterator(basicVar); !i.atEnd(); ++i){
715
750994
    const Tableau::Entry& entry = *i;
716
750994
    ArithVar v = entry.getColVar();
717
750994
    const Rational& coeff = entry.getCoefficient();
718
750994
    bool weakening = false;
719
750994
    ConstraintP c = weakestExplanation(aboveUpper, surplus, v, coeff, weakening, basicVar);
720
1501988
    Debug("arith::weak") << "weak : " << weakening << " "
721
1501988
                         << c->assertedToTheTheory() << " "
722
750994
                         << d_variables.getAssignment(v) << " "
723
750994
                         << c << endl;
724
750994
    anyWeakenings = anyWeakenings || weakening;
725
726
750994
    fcs.addConstraint(c, coeff, adjustSgn);
727
750994
    if(basicVar == v){
728
41940
      Assert(!c->negationHasProof());
729
41940
      fcs.makeLastConsequent();
730
    }
731
  }
732
41940
  Assert(fcs.consequentIsSet());
733
734
41940
  ConstraintCP conflicted = fcs.commitConflict();
735
736
41940
  ++d_statistics.d_weakeningAttempts;
737
41940
  if(anyWeakenings){
738
32424
    ++d_statistics.d_weakeningSuccesses;
739
  }
740
83880
  Debug("arith::weak") << "LinearEqualityModule::minimallyWeakConflict("
741
41940
                       << aboveUpper <<", "<< basicVar << ", ...) done" << endl;
742
83880
  return conflicted;
743
}
744
745
122822
ArithVar LinearEqualityModule::minVarOrder(ArithVar x, ArithVar y) const {
746
122822
  Assert(x != ARITHVAR_SENTINEL);
747
122822
  Assert(y != ARITHVAR_SENTINEL);
748
122822
  if(x <= y){
749
76085
    return x;
750
  } else {
751
46737
    return y;
752
  }
753
}
754
755
746597
ArithVar LinearEqualityModule::minColLength(ArithVar x, ArithVar y) const {
756
746597
  Assert(x != ARITHVAR_SENTINEL);
757
746597
  Assert(y != ARITHVAR_SENTINEL);
758
746597
  Assert(!d_tableau.isBasic(x));
759
746597
  Assert(!d_tableau.isBasic(y));
760
746597
  uint32_t xLen = d_tableau.getColLength(x);
761
746597
  uint32_t yLen = d_tableau.getColLength(y);
762
746597
  if( xLen > yLen){
763
91844
     return y;
764
654753
  } else if( xLen== yLen ){
765
122314
    return minVarOrder(x,y);
766
  }else{
767
532439
    return x;
768
  }
769
}
770
771
ArithVar LinearEqualityModule::minRowLength(ArithVar x, ArithVar y) const {
772
  Assert(x != ARITHVAR_SENTINEL);
773
  Assert(y != ARITHVAR_SENTINEL);
774
  Assert(d_tableau.isBasic(x));
775
  Assert(d_tableau.isBasic(y));
776
  uint32_t xLen = d_tableau.basicRowLength(x);
777
  uint32_t yLen = d_tableau.basicRowLength(y);
778
  if( xLen > yLen){
779
     return y;
780
  } else if( xLen== yLen ){
781
    return minVarOrder(x,y);
782
  }else{
783
    return x;
784
  }
785
}
786
787
771944
ArithVar LinearEqualityModule::minBoundAndColLength(ArithVar x, ArithVar y) const{
788
771944
  Assert(x != ARITHVAR_SENTINEL);
789
771944
  Assert(y != ARITHVAR_SENTINEL);
790
771944
  Assert(!d_tableau.isBasic(x));
791
771944
  Assert(!d_tableau.isBasic(y));
792
771944
  if(d_variables.hasEitherBound(x) && !d_variables.hasEitherBound(y)){
793
8200
    return y;
794
763744
  }else if(!d_variables.hasEitherBound(x) && d_variables.hasEitherBound(y)){
795
17147
    return x;
796
  }else {
797
746597
    return minColLength(x, y);
798
  }
799
}
800
801
template <bool above>
802
134892
ArithVar LinearEqualityModule::selectSlack(ArithVar x_i, VarPreferenceFunction pref) const{
803
134892
  ArithVar slack = ARITHVAR_SENTINEL;
804
805
4572092
  for(Tableau::RowIterator iter = d_tableau.basicRowIterator(x_i); !iter.atEnd();  ++iter){
806
4437200
    const Tableau::Entry& entry = *iter;
807
4437200
    ArithVar nonbasic = entry.getColVar();
808
4437200
    if(nonbasic == x_i) continue;
809
810
4302308
    const Rational& a_ij = entry.getCoefficient();
811
4302308
    int sgn = a_ij.sgn();
812
4302308
    if(isAcceptableSlack<above>(sgn, nonbasic)){
813
      //If one of the above conditions is met, we have found an acceptable
814
      //nonbasic variable to pivot x_i with.  We can now choose which one we
815
      //prefer the most.
816
907344
      slack = (slack == ARITHVAR_SENTINEL) ? nonbasic : (this->*pref)(slack, nonbasic);
817
    }
818
  }
819
820
134892
  return slack;
821
}
822
823
const Tableau::Entry* LinearEqualityModule::selectSlackEntry(ArithVar x_i, bool above) const{
824
  for(Tableau::RowIterator iter = d_tableau.basicRowIterator(x_i); !iter.atEnd();  ++iter){
825
    const Tableau::Entry& entry = *iter;
826
    ArithVar nonbasic = entry.getColVar();
827
    if(nonbasic == x_i) continue;
828
829
    const Rational& a_ij = entry.getCoefficient();
830
    int sgn = a_ij.sgn();
831
    if(above && isAcceptableSlack<true>(sgn, nonbasic)){
832
      //If one of the above conditions is met, we have found an acceptable
833
      //nonbasic variable to pivot x_i with.  We can now choose which one we
834
      //prefer the most.
835
      return &entry;
836
    }else if(!above && isAcceptableSlack<false>(sgn, nonbasic)){
837
      return &entry;
838
    }
839
  }
840
841
  return NULL;
842
}
843
844
1128716
void LinearEqualityModule::startTrackingBoundCounts(){
845
1128716
  Assert(!d_areTracking);
846
1128716
  d_areTracking = true;
847
1128716
  if(Debug.isOn("arith::tracking")){
848
    debugCheckTracking();
849
  }
850
1128716
  Assert(d_areTracking);
851
1128716
}
852
853
1128716
void LinearEqualityModule::stopTrackingBoundCounts(){
854
1128716
  Assert(d_areTracking);
855
1128716
  d_areTracking = false;
856
1128716
  if(Debug.isOn("arith::tracking")){
857
    debugCheckTracking();
858
  }
859
1128716
  Assert(!d_areTracking);
860
1128716
}
861
862
863
50910
void LinearEqualityModule::trackRowIndex(RowIndex ridx){
864
50910
  Assert(!rowIndexIsTracked(ridx));
865
50910
  BoundsInfo bi = computeRowBoundInfo(ridx, true);
866
50910
  d_btracking.set(ridx, bi);
867
50910
}
868
869
50910
BoundsInfo LinearEqualityModule::computeRowBoundInfo(RowIndex ridx, bool inQueue) const{
870
50910
  BoundsInfo bi;
871
872
50910
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
873
570446
  for(; !iter.atEnd();  ++iter){
874
259768
    const Tableau::Entry& entry = *iter;
875
259768
    ArithVar v = entry.getColVar();
876
259768
    const Rational& a_ij = entry.getCoefficient();
877
259768
    bi += (d_variables.selectBoundsInfo(v, inQueue)).multiplyBySgn(a_ij.sgn());
878
  }
879
50910
  return bi;
880
}
881
882
BoundCounts LinearEqualityModule::debugBasicAtBoundCount(ArithVar x_i) const {
883
  return d_btracking[d_tableau.basicToRowIndex(x_i)].atBounds();
884
}
885
886
/**
887
 * If the pivot described in u were performed,
888
 * then the row would qualify as being either at the minimum/maximum
889
 * to the non-basics being at their bounds.
890
 * The minimum/maximum is determined by the direction the non-basic is changing.
891
 */
892
bool LinearEqualityModule::basicsAtBounds(const UpdateInfo& u) const {
893
  Assert(u.describesPivot());
894
895
  ArithVar nonbasic = u.nonbasic();
896
  ArithVar basic = u.leaving();
897
  Assert(basicIsTracked(basic));
898
  int coeffSgn = u.getCoefficient().sgn();
899
  int nbdir = u.nonbasicDirection();
900
901
  ConstraintP c = u.limiting();
902
  int toUB = (c->getType() == UpperBound ||
903
              c->getType() == Equality) ? 1 : 0;
904
  int toLB = (c->getType() == LowerBound ||
905
              c->getType() == Equality) ? 1 : 0;
906
907
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
908
909
  BoundCounts bcs = d_btracking[ridx].atBounds();
910
  // x = c*n + \sum d*m
911
  // 0 = -x + c*n + \sum d*m
912
  // n = 1/c * x + -1/c * (\sum d*m)
913
  BoundCounts nonb = bcs - d_variables.atBoundCounts(nonbasic).multiplyBySgn(coeffSgn);
914
  nonb.addInChange(-1, d_variables.atBoundCounts(basic), BoundCounts(toLB, toUB));
915
  nonb = nonb.multiplyBySgn(-coeffSgn);
916
917
  uint32_t length = d_tableau.basicRowLength(basic);
918
  Debug("basicsAtBounds")
919
    << "bcs " << bcs
920
    << "nonb " << nonb
921
    << "length " << length << endl;
922
  // nonb has nb excluded.
923
  if(nbdir < 0){
924
    return nonb.lowerBoundCount() + 1 == length;
925
  }else{
926
    Assert(nbdir > 0);
927
    return nonb.upperBoundCount() + 1 == length;
928
  }
929
}
930
931
254123
bool LinearEqualityModule::nonbasicsAtLowerBounds(ArithVar basic) const {
932
254123
  Assert(basicIsTracked(basic));
933
254123
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
934
935
254123
  BoundCounts bcs = d_btracking[ridx].atBounds();
936
254123
  uint32_t length = d_tableau.basicRowLength(basic);
937
938
  // return true if excluding the basic is every element is at its "lowerbound"
939
  // The psuedo code is:
940
  //   bcs -= basic.count(basic, basic's sgn)
941
  //   return bcs.lowerBoundCount() + 1 == length
942
  // As basic's sign is always -1, we can pull out the pieces of the count:
943
  //   bcs.lowerBoundCount() - basic.atUpperBoundInd() + 1 == length
944
  // basic.atUpperBoundInd() is either 0 or 1
945
254123
  uint32_t lbc = bcs.lowerBoundCount();
946
590786
  return  (lbc == length) ||
947
590786
    (lbc + 1 == length && d_variables.cmpAssignmentUpperBound(basic) != 0);
948
}
949
950
314763
bool LinearEqualityModule::nonbasicsAtUpperBounds(ArithVar basic) const {
951
314763
  Assert(basicIsTracked(basic));
952
314763
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
953
314763
  BoundCounts bcs = d_btracking[ridx].atBounds();
954
314763
  uint32_t length = d_tableau.basicRowLength(basic);
955
314763
  uint32_t ubc = bcs.upperBoundCount();
956
  // See the comment for nonbasicsAtLowerBounds()
957
958
756686
  return (ubc == length) ||
959
756686
    (ubc + 1 == length && d_variables.cmpAssignmentLowerBound(basic) != 0);
960
}
961
962
134892
void LinearEqualityModule::trackingMultiplyRow(RowIndex ridx, int sgn) {
963
134892
  Assert(rowIndexIsTracked(ridx));
964
134892
  Assert(sgn != 0);
965
134892
  if(sgn < 0){
966
77082
    BoundsInfo& bi = d_btracking.get(ridx);
967
77082
    bi = bi.multiplyBySgn(sgn);
968
  }
969
134892
}
970
971
44223216
void LinearEqualityModule::trackingCoefficientChange(RowIndex ridx, ArithVar nb, int oldSgn, int currSgn){
972
44223216
  Assert(oldSgn != currSgn);
973
44223216
  BoundsInfo nb_inf = d_variables.boundsInfo(nb);
974
975
44223216
  Assert(rowIndexIsTracked(ridx));
976
977
44223216
  BoundsInfo& row_bi = d_btracking.get(ridx);
978
44223216
  row_bi.addInSgn(nb_inf, oldSgn, currSgn);
979
44223216
}
980
981
ArithVar LinearEqualityModule::minBy(const ArithVarVec& vec, VarPreferenceFunction pf) const{
982
  if(vec.empty()) {
983
    return ARITHVAR_SENTINEL;
984
  }else {
985
    ArithVar sel = vec.front();
986
    ArithVarVec::const_iterator i = vec.begin() + 1;
987
    ArithVarVec::const_iterator i_end = vec.end();
988
    for(; i != i_end; ++i){
989
      sel = (this->*pf)(sel, *i);
990
    }
991
    return sel;
992
  }
993
}
994
995
bool LinearEqualityModule::accumulateBorder(const Tableau::Entry& entry, bool ub){
996
  ArithVar currBasic = d_tableau.rowIndexToBasic(entry.getRowIndex());
997
998
  Assert(basicIsTracked(currBasic));
999
1000
  ConstraintP bound = ub ?
1001
    d_variables.getUpperBoundConstraint(currBasic):
1002
    d_variables.getLowerBoundConstraint(currBasic);
1003
1004
  if(bound == NullConstraint){ return false; }
1005
  Assert(bound != NullConstraint);
1006
1007
  const Rational& coeff = entry.getCoefficient();
1008
1009
  const DeltaRational& assignment = d_variables.getAssignment(currBasic);
1010
  DeltaRational toBound = bound->getValue() - assignment;
1011
  DeltaRational nbDiff = toBound/coeff;
1012
1013
  // if ub
1014
  // if toUB >= 0
1015
  // then ub >= currBasic
1016
  //   if sgn > 0,
1017
  //   then diff >= 0, so nb must increase for G
1018
  //   else diff <= 0, so nb must decrease for G
1019
  // else ub < currBasic
1020
  //   if sgn > 0,
1021
  //   then diff < 0, so nb must decrease for G
1022
  //   else diff > 0, so nb must increase for G
1023
1024
  int diffSgn = nbDiff.sgn();
1025
1026
  if(diffSgn != 0 && willBeInConflictAfterPivot(entry, nbDiff, ub)){
1027
    return true;
1028
  }else{
1029
    bool areFixing = ub ? (toBound.sgn() < 0 ) : (toBound.sgn() > 0);
1030
    Border border(bound, nbDiff, areFixing, &entry, ub);
1031
    bool increasing =
1032
      (diffSgn > 0) ||
1033
      (diffSgn == 0 && ((coeff.sgn() > 0) == ub));
1034
1035
    // assume diffSgn == 0
1036
    // if coeff > 0,
1037
    //   if ub, inc
1038
    //   else, dec
1039
    // else coeff < 0
1040
    //   if ub, dec
1041
    //   else, inc
1042
1043
    if(increasing){
1044
      Debug("handleBorders") << "push back increasing " << border << endl;
1045
      d_increasing.push_back(border);
1046
    }else{
1047
      Debug("handleBorders") << "push back decreasing " << border << endl;
1048
      d_decreasing.push_back(border);
1049
    }
1050
    return false;
1051
  }
1052
}
1053
1054
bool LinearEqualityModule::willBeInConflictAfterPivot(const Tableau::Entry& entry, const DeltaRational& nbDiff, bool bToUB) const{
1055
  int nbSgn = nbDiff.sgn();
1056
  Assert(nbSgn != 0);
1057
1058
  if(nbSgn > 0){
1059
    if (!d_upperBoundDifference || nbDiff <= *d_upperBoundDifference)
1060
    {
1061
      return false;
1062
    }
1063
  }else{
1064
    if (!d_lowerBoundDifference || nbDiff >= *d_lowerBoundDifference)
1065
    {
1066
      return false;
1067
    }
1068
  }
1069
1070
  // Assume past this point, nb will be in error if this pivot is done
1071
  ArithVar nb = entry.getColVar();
1072
  RowIndex ridx = entry.getRowIndex();
1073
  ArithVar basic = d_tableau.rowIndexToBasic(ridx);
1074
  Assert(rowIndexIsTracked(ridx));
1075
  int coeffSgn = entry.getCoefficient().sgn();
1076
1077
1078
  // if bToUB, then basic is going to be set to its upperbound
1079
  // if not bToUB, then basic is going to be set to its lowerbound
1080
1081
  // Different steps of solving for this:
1082
  // 1) y = a * x + \sum b * z
1083
  // 2) -a * x = -y + \sum b * z
1084
  // 3) x = (-1/a) * ( -y + \sum b * z)
1085
1086
  BoundCounts bc = d_btracking[ridx].atBounds();
1087
1088
  // 1) y = a * x + \sum b * z
1089
  // Get bc(\sum b * z)
1090
  BoundCounts sumOnly = bc - d_variables.atBoundCounts(nb).multiplyBySgn(coeffSgn);
1091
1092
  // y's bounds in the proposed model
1093
  int yWillBeAtUb = (bToUB || d_variables.boundsAreEqual(basic)) ? 1 : 0;
1094
  int yWillBeAtLb = (!bToUB || d_variables.boundsAreEqual(basic)) ? 1 : 0;
1095
  BoundCounts ysBounds(yWillBeAtLb, yWillBeAtUb);
1096
1097
  // 2) -a * x = -y + \sum b * z
1098
  // Get bc(-y + \sum b * z)
1099
  sumOnly.addInChange(-1, d_variables.atBoundCounts(basic), ysBounds);
1100
1101
  // 3) x = (-1/a) * ( -y + \sum b * z)
1102
  // Get bc((-1/a) * ( -y + \sum b * z))
1103
  BoundCounts xsBoundsAfterPivot = sumOnly.multiplyBySgn(-coeffSgn);
1104
1105
  uint32_t length = d_tableau.basicRowLength(basic);
1106
  if(nbSgn > 0){
1107
    // Only check for the upper bound being violated
1108
    return xsBoundsAfterPivot.lowerBoundCount() + 1 == length;
1109
  }else{
1110
    // Only check for the lower bound being violated
1111
    return xsBoundsAfterPivot.upperBoundCount() + 1 == length;
1112
  }
1113
}
1114
1115
UpdateInfo LinearEqualityModule::mkConflictUpdate(const Tableau::Entry& entry, bool ub) const{
1116
  ArithVar currBasic = d_tableau.rowIndexToBasic(entry.getRowIndex());
1117
  ArithVar nb = entry.getColVar();
1118
1119
  ConstraintP bound = ub ?
1120
    d_variables.getUpperBoundConstraint(currBasic):
1121
    d_variables.getLowerBoundConstraint(currBasic);
1122
1123
1124
  const Rational& coeff = entry.getCoefficient();
1125
  const DeltaRational& assignment = d_variables.getAssignment(currBasic);
1126
  DeltaRational toBound = bound->getValue() - assignment;
1127
  DeltaRational nbDiff = toBound/coeff;
1128
1129
  return UpdateInfo::conflict(nb, nbDiff, coeff, bound);
1130
}
1131
1132
UpdateInfo LinearEqualityModule::speculativeUpdate(ArithVar nb, const Rational& focusCoeff, UpdatePreferenceFunction pref){
1133
  Assert(d_increasing.empty());
1134
  Assert(d_decreasing.empty());
1135
  Assert(!d_lowerBoundDifference);
1136
  Assert(!d_upperBoundDifference);
1137
1138
  int focusCoeffSgn = focusCoeff.sgn();
1139
1140
  static int instance = 0;
1141
  ++instance;
1142
  Debug("speculativeUpdate") << "speculativeUpdate " << instance << endl;
1143
  Debug("speculativeUpdate") << "nb " << nb << endl;
1144
  Debug("speculativeUpdate") << "focusCoeff " << focusCoeff << endl;
1145
1146
  if(d_variables.hasUpperBound(nb)){
1147
    ConstraintP ub = d_variables.getUpperBoundConstraint(nb);
1148
    d_upperBoundDifference = ub->getValue() - d_variables.getAssignment(nb);
1149
    Border border(ub, *d_upperBoundDifference, false, NULL, true);
1150
    Debug("handleBorders") << "push back increasing " << border << endl;
1151
    d_increasing.push_back(border);
1152
  }
1153
  if(d_variables.hasLowerBound(nb)){
1154
    ConstraintP lb = d_variables.getLowerBoundConstraint(nb);
1155
    d_lowerBoundDifference = lb->getValue() - d_variables.getAssignment(nb);
1156
    Border border(lb, *d_lowerBoundDifference, false, NULL, false);
1157
    Debug("handleBorders") << "push back decreasing " << border << endl;
1158
    d_decreasing.push_back(border);
1159
  }
1160
1161
  Tableau::ColIterator colIter = d_tableau.colIterator(nb);
1162
  for(; !colIter.atEnd(); ++colIter){
1163
    const Tableau::Entry& entry = *colIter;
1164
    Assert(entry.getColVar() == nb);
1165
1166
    if(accumulateBorder(entry, true)){
1167
      clearSpeculative();
1168
      return mkConflictUpdate(entry, true);
1169
    }
1170
    if(accumulateBorder(entry, false)){
1171
      clearSpeculative();
1172
      return mkConflictUpdate(entry, false);
1173
    }
1174
  }
1175
1176
  UpdateInfo selected;
1177
  BorderHeap& withSgn = focusCoeffSgn > 0 ? d_increasing : d_decreasing;
1178
  BorderHeap& againstSgn = focusCoeffSgn > 0 ? d_decreasing : d_increasing;
1179
1180
  handleBorders(selected, nb, focusCoeff, withSgn, 0, pref);
1181
  int m = 1 - selected.errorsChangeSafe(0);
1182
  handleBorders(selected, nb, focusCoeff, againstSgn, m, pref);
1183
1184
  clearSpeculative();
1185
  return selected;
1186
}
1187
1188
void LinearEqualityModule::clearSpeculative(){
1189
  // clear everything away
1190
  d_increasing.clear();
1191
  d_decreasing.clear();
1192
  d_lowerBoundDifference.reset();
1193
  d_upperBoundDifference.reset();
1194
}
1195
1196
void LinearEqualityModule::handleBorders(UpdateInfo& selected, ArithVar nb, const Rational& focusCoeff, BorderHeap& heap, int minimumFixes, UpdatePreferenceFunction pref){
1197
  Assert(minimumFixes >= 0);
1198
1199
  // The values popped off of the heap
1200
  // should be popped with the values closest to 0
1201
  // being first and larger in absolute value last
1202
1203
1204
  int fixesRemaining = heap.possibleFixes();
1205
1206
  Debug("handleBorders")
1207
    << "handleBorders "
1208
    << "nb " << nb
1209
    << "fc " << focusCoeff
1210
    << "h.e " << heap.empty()
1211
    << "h.dir " << heap.direction()
1212
    << "h.rem " << fixesRemaining
1213
    << "h.0s " << heap.numZeroes()
1214
    << "min " << minimumFixes
1215
    << endl;
1216
1217
  if(heap.empty()){
1218
    // if the heap is empty, return
1219
    return;
1220
  }
1221
1222
  bool zeroesWillDominate = fixesRemaining - heap.numZeroes() < minimumFixes;
1223
1224
  // can the number of fixes ever exceed the minimum?
1225
  // no more than the number of possible fixes can be fixed in total
1226
  // nothing can be fixed before the zeroes are taken care of
1227
  if(minimumFixes > 0 && zeroesWillDominate){
1228
    return;
1229
  }
1230
1231
1232
  int negErrorChange = 0;
1233
  int nbDir = heap.direction();
1234
1235
  // points at the beginning of the heap
1236
  if(zeroesWillDominate){
1237
    heap.dropNonZeroes();
1238
  }
1239
  heap.make_heap();
1240
1241
1242
  // pretend like the previous block had a value of zero.
1243
  // The block that actually has a value of 0 must handle this.
1244
  const DeltaRational zero(0);
1245
  const DeltaRational* prevBlockValue = &zero;
1246
1247
  /** The coefficient changes as the value crosses border. */
1248
  Rational effectiveCoefficient = focusCoeff;
1249
1250
  /* Keeps track of the change to the value of the focus function.*/
1251
  DeltaRational totalFocusChange(0);
1252
1253
  const int focusCoeffSgn = focusCoeff.sgn();
1254
1255
  while(heap.more()  &&
1256
        (fixesRemaining + negErrorChange > minimumFixes ||
1257
         (fixesRemaining + negErrorChange == minimumFixes &&
1258
          effectiveCoefficient.sgn() == focusCoeffSgn))){
1259
    // There are more elements &&
1260
    // we can either fix at least 1 more variable in the error function
1261
    // or we can improve the error function
1262
1263
1264
    int brokenInBlock = 0;
1265
    BorderVec::const_iterator endBlock = heap.end();
1266
1267
    pop_block(heap, brokenInBlock, fixesRemaining, negErrorChange);
1268
1269
    // if endVec == beginVec, block starts there
1270
    // other wise, block starts at endVec
1271
    BorderVec::const_iterator startBlock
1272
      = heap.more() ? heap.end() : heap.begin();
1273
1274
    const DeltaRational& blockValue = (*startBlock).d_diff;
1275
1276
    // if decreasing
1277
    // blockValue < prevBlockValue
1278
    // diff.sgn() = -1
1279
    DeltaRational diff = blockValue - (*prevBlockValue);
1280
    DeltaRational blockChangeToFocus =  diff * effectiveCoefficient;
1281
    totalFocusChange += blockChangeToFocus;
1282
1283
    Debug("handleBorders")
1284
      << "blockValue " << (blockValue)
1285
      << "diff " << diff
1286
      << "blockChangeToFocus " << totalFocusChange
1287
      << "blockChangeToFocus " << totalFocusChange
1288
      << "negErrorChange " << negErrorChange
1289
      << "brokenInBlock " << brokenInBlock
1290
      << "fixesRemaining " << fixesRemaining
1291
      << endl;
1292
1293
    int currFocusChangeSgn = totalFocusChange.sgn();
1294
    for(BorderVec::const_iterator i = startBlock; i != endBlock; ++i){
1295
      const Border& b = *i;
1296
1297
      Debug("handleBorders") << b << endl;
1298
1299
      bool makesImprovement = negErrorChange > 0 ||
1300
        (negErrorChange == 0  && currFocusChangeSgn > 0);
1301
1302
      if(!makesImprovement){
1303
        if(b.ownBorder() || minimumFixes > 0){
1304
          continue;
1305
        }
1306
      }
1307
1308
      UpdateInfo proposal(nb, nbDir);
1309
      if(b.ownBorder()){
1310
        proposal.witnessedUpdate(b.d_diff, b.d_bound, -negErrorChange, currFocusChangeSgn);
1311
      }else{
1312
        proposal.update(b.d_diff, b.getCoefficient(), b.d_bound, -negErrorChange, currFocusChangeSgn);
1313
      }
1314
1315
      if(selected.unbounded() || (this->*pref)(selected, proposal)){
1316
        selected = proposal;
1317
      }
1318
    }
1319
1320
    effectiveCoefficient += updateCoefficient(startBlock, endBlock);
1321
    prevBlockValue = &blockValue;
1322
    negErrorChange -= brokenInBlock;
1323
  }
1324
}
1325
1326
Rational LinearEqualityModule::updateCoefficient(BorderVec::const_iterator startBlock, BorderVec::const_iterator endBlock){
1327
  //update coefficient
1328
  Rational changeToCoefficient(0);
1329
  for(BorderVec::const_iterator i = startBlock; i != endBlock; ++i){
1330
    const Border& curr = *i;
1331
    if(curr.ownBorder()){// breaking its own bound
1332
      if(curr.d_upperbound){
1333
        changeToCoefficient -= 1;
1334
      }else{
1335
        changeToCoefficient += 1;
1336
      }
1337
    }else{
1338
      const Rational& coeff = curr.d_entry->getCoefficient();
1339
      if(curr.d_areFixing){
1340
        if(curr.d_upperbound){// fixing an upper bound
1341
          changeToCoefficient += coeff;
1342
        }else{// fixing a lower bound
1343
          changeToCoefficient -= coeff;
1344
        }
1345
      }else{
1346
        if(curr.d_upperbound){// breaking an upper bound
1347
          changeToCoefficient -= coeff;
1348
        }else{
1349
          // breaking a lower bound
1350
          changeToCoefficient += coeff;
1351
        }
1352
      }
1353
    }
1354
  }
1355
  return changeToCoefficient;
1356
}
1357
1358
void LinearEqualityModule::pop_block(BorderHeap& heap, int& brokenInBlock, int& fixesRemaining, int& negErrorChange){
1359
  Assert(heap.more());
1360
1361
  if(heap.top().d_areFixing){
1362
    fixesRemaining--;
1363
    negErrorChange++;
1364
  }else{
1365
    brokenInBlock++;
1366
  }
1367
  heap.pop_heap();
1368
  const DeltaRational& blockValue = (*heap.end()).d_diff;
1369
1370
  while(heap.more()){
1371
    const Border& top = heap.top();
1372
    if(blockValue == top.d_diff){
1373
      // belongs to the block
1374
      if(top.d_areFixing){
1375
        fixesRemaining--;
1376
        negErrorChange++;
1377
      }else{
1378
        brokenInBlock++;
1379
      }
1380
      heap.pop_heap();
1381
    }else{
1382
      // does not belong to the block
1383
      Assert((heap.direction() > 0) ? (blockValue < top.d_diff)
1384
                                    : (blockValue > top.d_diff));
1385
      break;
1386
    }
1387
  }
1388
}
1389
1390
void LinearEqualityModule::substitutePlusTimesConstant(ArithVar to, ArithVar from, const Rational& mult){
1391
  d_tableau.substitutePlusTimesConstant(to, from, mult, d_trackCallback);
1392
}
1393
void LinearEqualityModule::directlyAddToCoefficient(ArithVar row, ArithVar col, const Rational& mult){
1394
  d_tableau.directlyAddToCoefficient(row, col, mult, d_trackCallback);
1395
}
1396
1397
}  // namespace arith
1398
}  // namespace theory
1399
22746
}  // namespace cvc5