GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/linear_equality.cpp Lines: 368 804 45.8 %
Date: 2021-09-30 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
6278
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
6278
  d_trackCallback(this)
66
6278
{}
67
68
6278
LinearEqualityModule::Statistics::Statistics()
69
    : d_statPivots(
70
12556
        smtStatisticsRegistry().registerInt("theory::arith::pivots")),
71
      d_statUpdates(
72
12556
          smtStatisticsRegistry().registerInt("theory::arith::updates")),
73
      d_pivotTime(
74
12556
          smtStatisticsRegistry().registerTimer("theory::arith::pivotTime")),
75
      d_adjTime(
76
12556
          smtStatisticsRegistry().registerTimer("theory::arith::adjTime")),
77
6278
      d_weakeningAttempts(smtStatisticsRegistry().registerInt(
78
12556
          "theory::arith::weakening::attempts")),
79
6278
      d_weakeningSuccesses(smtStatisticsRegistry().registerInt(
80
12556
          "theory::arith::weakening::success")),
81
6278
      d_weakenings(smtStatisticsRegistry().registerInt(
82
12556
          "theory::arith::weakening::total")),
83
6278
      d_weakenTime(smtStatisticsRegistry().registerTimer(
84
12556
          "theory::arith::weakening::time")),
85
      d_forceTime(
86
56502
          smtStatisticsRegistry().registerTimer("theory::arith::forcing::time"))
87
{
88
6278
}
89
90
3004524
void LinearEqualityModule::includeBoundUpdate(ArithVar v, const BoundsInfo& prev){
91
3004524
  Assert(!d_areTracking);
92
93
3004524
  BoundsInfo curr = d_variables.boundsInfo(v);
94
95
3004524
  Assert(prev != curr);
96
3004524
  Tableau::ColIterator basicIter = d_tableau.colIterator(v);
97
29672462
  for(; !basicIter.atEnd(); ++basicIter){
98
13333969
    const Tableau::Entry& entry = *basicIter;
99
13333969
    Assert(entry.getColVar() == v);
100
13333969
    int a_ijSgn = entry.getCoefficient().sgn();
101
102
13333969
    RowIndex ridx = entry.getRowIndex();
103
13333969
    BoundsInfo& counts = d_btracking.get(ridx);
104
13333969
    Debug("includeBoundUpdate") << d_tableau.rowIndexToBasic(ridx) << " " << counts << " to " ;
105
13333969
    counts.addInChange(a_ijSgn, prev, curr);
106
13333969
    Debug("includeBoundUpdate") << counts << " " << a_ijSgn << std::endl;
107
  }
108
3004524
}
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
84111
void LinearEqualityModule::updateUntracked(ArithVar x_i, const DeltaRational& v){
181
84111
  Assert(!d_tableau.isBasic(x_i));
182
84111
  Assert(!d_areTracking);
183
84111
  const DeltaRational& assignment_x_i = d_variables.getAssignment(x_i);
184
84111
  ++(d_statistics.d_statUpdates);
185
186
187
168222
  Debug("arith") <<"update " << x_i << ": "
188
84111
                 << assignment_x_i << "|-> " << v << endl;
189
168222
  DeltaRational diff = v - assignment_x_i;
190
191
84111
  Tableau::ColIterator colIter = d_tableau.colIterator(x_i);
192
1865033
  for(; !colIter.atEnd(); ++colIter){
193
890461
    const Tableau::Entry& entry = *colIter;
194
890461
    Assert(entry.getColVar() == x_i);
195
196
890461
    ArithVar x_j = d_tableau.rowIndexToBasic(entry.getRowIndex());
197
890461
    const Rational& a_ji = entry.getCoefficient();
198
199
890461
    const DeltaRational& assignment = d_variables.getAssignment(x_j);
200
1780922
    DeltaRational  nAssignment = assignment+(diff * a_ji);
201
890461
    d_variables.setAssignment(x_j, nAssignment);
202
203
890461
    d_basicVariableUpdates(x_j);
204
  }
205
206
84111
  d_variables.setAssignment(x_i, v);
207
208
84111
  if(Debug.isOn("paranoid:check_tableau")){  debugCheckTableau(); }
209
84111
}
210
211
136052
void LinearEqualityModule::updateTracked(ArithVar x_i, const DeltaRational& v){
212
272104
  TimerStat::CodeTimer codeTimer(d_statistics.d_adjTime);
213
214
136052
  Assert(!d_tableau.isBasic(x_i));
215
136052
  Assert(d_areTracking);
216
217
136052
  ++(d_statistics.d_statUpdates);
218
219
272104
  DeltaRational diff =  v - d_variables.getAssignment(x_i);
220
272104
  Debug("arith") <<"update " << x_i << ": "
221
136052
                 << d_variables.getAssignment(x_i) << "|-> " << v << endl;
222
223
224
136052
  BoundCounts before = d_variables.atBoundCounts(x_i);
225
136052
  d_variables.setAssignment(x_i, v);
226
136052
  BoundCounts after = d_variables.atBoundCounts(x_i);
227
228
136052
  bool anyChange = before != after;
229
230
136052
  Tableau::ColIterator colIter = d_tableau.colIterator(x_i);
231
4962292
  for(; !colIter.atEnd(); ++colIter){
232
2413120
    const Tableau::Entry& entry = *colIter;
233
2413120
    Assert(entry.getColVar() == x_i);
234
235
2413120
    RowIndex ridx = entry.getRowIndex();
236
2413120
    ArithVar x_j = d_tableau.rowIndexToBasic(ridx);
237
2413120
    const Rational& a_ji = entry.getCoefficient();
238
239
2413120
    const DeltaRational& assignment = d_variables.getAssignment(x_j);
240
4826240
    DeltaRational  nAssignment = assignment+(diff * a_ji);
241
2413120
    Debug("update") << x_j << " " << a_ji << assignment << " -> " << nAssignment << endl;
242
2413120
    BoundCounts xjBefore = d_variables.atBoundCounts(x_j);
243
2413120
    d_variables.setAssignment(x_j, nAssignment);
244
2413120
    BoundCounts xjAfter = d_variables.atBoundCounts(x_j);
245
246
2413120
    Assert(rowIndexIsTracked(ridx));
247
2413120
    BoundsInfo& next_bc_k = d_btracking.get(ridx);
248
2413120
    if(anyChange){
249
1876513
      next_bc_k.addInAtBoundChange(a_ji.sgn(), before, after);
250
    }
251
2413120
    if(xjBefore != xjAfter){
252
279728
      next_bc_k.addInAtBoundChange(-1, xjBefore, xjAfter);
253
    }
254
255
2413120
    d_basicVariableUpdates(x_j);
256
  }
257
258
136052
  if(Debug.isOn("paranoid:check_tableau")){  debugCheckTableau(); }
259
136052
}
260
261
136052
void LinearEqualityModule::pivotAndUpdate(ArithVar x_i, ArithVar x_j, const DeltaRational& x_i_value){
262
136052
  Assert(x_i != x_j);
263
264
272104
  TimerStat::CodeTimer codeTimer(d_statistics.d_pivotTime);
265
266
  static int instance = 0;
267
268
136052
  if(Debug.isOn("arith::tracking::pre")){
269
    ++instance;
270
    Debug("arith::tracking")  << "pre update #" << instance << endl;
271
    debugCheckTracking();
272
  }
273
274
136052
  if(Debug.isOn("arith::simplex:row")){ debugPivot(x_i, x_j); }
275
276
136052
  RowIndex ridx = d_tableau.basicToRowIndex(x_i);
277
136052
  const Tableau::Entry& entry_ij =  d_tableau.findEntry(ridx, x_j);
278
136052
  Assert(!entry_ij.blank());
279
280
136052
  const Rational& a_ij = entry_ij.getCoefficient();
281
136052
  const DeltaRational& betaX_i = d_variables.getAssignment(x_i);
282
272104
  DeltaRational theta = (x_i_value - betaX_i)/a_ij;
283
272104
  DeltaRational x_j_value = d_variables.getAssignment(x_j) + theta;
284
285
136052
  updateTracked(x_j, x_j_value);
286
287
136052
  if(Debug.isOn("arith::tracking::mid")){
288
    Debug("arith::tracking")  << "postupdate prepivot #" << instance << endl;
289
    debugCheckTracking();
290
  }
291
292
  // Pivots
293
136052
  ++(d_statistics.d_statPivots);
294
295
136052
  d_tableau.pivot(x_i, x_j, d_trackCallback);
296
297
136052
  if(Debug.isOn("arith::tracking::post")){
298
    Debug("arith::tracking")  << "postpivot #" << instance << endl;
299
    debugCheckTracking();
300
  }
301
302
136052
  d_basicVariableUpdates(x_j);
303
304
136052
  if(Debug.isOn("matrix")){
305
    d_tableau.printMatrix();
306
  }
307
136052
}
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
234088
DeltaRational LinearEqualityModule::computeRowBound(RowIndex ridx, bool rowUb, ArithVar skip) const {
427
234088
  DeltaRational sum(0,0);
428
1661398
  for(Tableau::RowIterator i = d_tableau.ridRowIterator(ridx); !i.atEnd(); ++i){
429
1427310
    const Tableau::Entry& entry = (*i);
430
1427310
    ArithVar v = entry.getColVar();
431
1427310
    if(v == skip){ continue; }
432
433
1242963
    const Rational& coeff =  entry.getCoefficient();
434
1242963
    bool vUb = (rowUb == (coeff.sgn() > 0));
435
436
2485926
    const DeltaRational& bound = vUb ?
437
629083
      d_variables.getUpperBound(v):
438
1856843
      d_variables.getLowerBound(v);
439
440
2485926
    DeltaRational diff = bound * coeff;
441
1242963
    sum = sum + diff;
442
  }
443
234088
  return sum;
444
}
445
446
/**
447
 * Computes the value of a basic variable using the current assignment.
448
 */
449
103664
DeltaRational LinearEqualityModule::computeRowValue(ArithVar x, bool useSafe) const{
450
103664
  Assert(d_tableau.isBasic(x));
451
103664
  DeltaRational sum(0);
452
453
642700
  for(Tableau::RowIterator i = d_tableau.basicRowIterator(x); !i.atEnd(); ++i){
454
539036
    const Tableau::Entry& entry = (*i);
455
539036
    ArithVar nonbasic = entry.getColVar();
456
539036
    if(nonbasic == x) continue;
457
435372
    const Rational& coeff = entry.getCoefficient();
458
459
435372
    const DeltaRational& assignment = d_variables.getAssignment(nonbasic, useSafe);
460
435372
    sum = sum + (assignment * coeff);
461
  }
462
103664
  return sum;
463
}
464
465
2932015
const Tableau::Entry* LinearEqualityModule::rowLacksBound(RowIndex ridx, bool rowUb, ArithVar skip){
466
2932015
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
467
28238437
  for(; !iter.atEnd(); ++iter){
468
15583480
    const Tableau::Entry& entry = *iter;
469
470
15583480
    ArithVar var = entry.getColVar();
471
15583480
    if(var == skip) { continue; }
472
473
15580788
    int sgn = entry.getCoefficient().sgn();
474
15580788
    bool selectUb = (rowUb == (sgn > 0));
475
31161576
    bool hasBound = selectUb ?
476
7802152
      d_variables.hasUpperBound(var):
477
23359424
      d_variables.hasLowerBound(var);
478
15580788
    if(!hasBound){
479
2930269
      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
83772
void LinearEqualityModule::propagateRow(ConstraintCPVec& into, RowIndex ridx, bool rowUp, ConstraintP c, RationalVectorP farkas){
535
83772
  Assert(!c->assertedToTheTheory());
536
83772
  Assert(c->canBePropagated());
537
83772
  Assert(!c->hasProof());
538
539
83772
  if(farkas != RationalVectorPSentinel){
540
1907
    Assert(farkas->empty());
541
1907
    farkas->push_back(Rational(0));
542
  }
543
544
83772
  ArithVar v = c->getVariable();
545
167544
  Debug("arith::propagateRow") << "LinearEqualityModule::propagateRow("
546
83772
                                   << ridx << ", " << rowUp << ", " << v << ") start" << endl;
547
548
83772
  const Rational& multiple = rowUp ? d_one : d_negOne;
549
550
83772
  Debug("arith::propagateRow") << "multiple: " << multiple << endl;
551
552
83772
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
553
1305942
  for(; !iter.atEnd(); ++iter){
554
611085
    const Tableau::Entry& entry = *iter;
555
611085
    ArithVar nonbasic = entry.getColVar();
556
611085
    const Rational& a_ij = entry.getCoefficient();
557
611085
    int sgn = a_ij.sgn();
558
611085
    Assert(sgn != 0);
559
611085
    bool selectUb = rowUp ? (sgn > 0) : (sgn < 0);
560
561
611085
    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
611085
    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
611085
    if(nonbasic == v){
580
83772
      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
83772
      Debug("arith::propagateRow") << c << endl;
588
    }else{
589
590
      ConstraintCP bound = selectUb
591
224268
        ? d_variables.getUpperBoundConstraint(nonbasic)
592
751581
        : d_variables.getLowerBoundConstraint(nonbasic);
593
594
527313
      if(farkas != RationalVectorPSentinel){
595
20696
        Rational multAij = multiple * a_ij;
596
10348
        Debug("arith::propagateRow") << "(" << multAij << ") ";
597
10348
        farkas->push_back(multAij);
598
      }
599
527313
      Assert(bound != NullConstraint);
600
527313
      Debug("arith::propagateRow") << bound << endl;
601
527313
      into.push_back(bound);
602
    }
603
  }
604
167544
  Debug("arith::propagateRow") << "LinearEqualityModule::propagateRow("
605
83772
                                   << ridx << ", " << rowUp << ", " << v << ") done" << endl;
606
607
83772
}
608
609
765896
ConstraintP LinearEqualityModule::weakestExplanation(bool aboveUpper, DeltaRational& surplus, ArithVar v, const Rational& coeff, bool& anyWeakening, ArithVar basic) const {
610
611
765896
  int sgn = coeff.sgn();
612
765896
  bool ub = aboveUpper?(sgn < 0) : (sgn > 0);
613
614
1531792
  ConstraintP c = ub ?
615
187837
    d_variables.getUpperBoundConstraint(v) :
616
1343955
    d_variables.getLowerBoundConstraint(v);
617
618
  bool weakened;
619
835575
  do{
620
835575
    const DeltaRational& bound = c->getValue();
621
622
835575
    weakened = false;
623
624
1671150
    ConstraintP weaker = ub?
625
235848
      c->getStrictlyWeakerUpperBound(true, true):
626
1435302
      c->getStrictlyWeakerLowerBound(true, true);
627
628
835575
    if(weaker != NullConstraint){
629
175013
      const DeltaRational& weakerBound = weaker->getValue();
630
631
350026
      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
175013
      diff = diff * coeff;
636
175013
      if(surplus > diff){
637
69679
        ++d_statistics.d_weakenings;
638
69679
        weakened = true;
639
69679
        anyWeakening = true;
640
69679
        surplus = surplus - diff;
641
642
69679
        Debug("arith::weak") << "found:" << endl;
643
69679
        if(v == basic){
644
7198
          Debug("arith::weak") << "  basic: ";
645
        }
646
139358
        Debug("arith::weak") << "  " << surplus << " "<< diff  << endl
647
69679
                             << "  " << bound << c << endl
648
69679
                             << "  " << weakerBound << weaker << endl;
649
650
69679
        Assert(diff.sgn() > 0);
651
69679
        c = weaker;
652
      }
653
    }
654
  }while(weakened);
655
656
765896
  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
43018
ConstraintCP LinearEqualityModule::minimallyWeakConflict(bool aboveUpper, ArithVar basicVar, FarkasConflictBuilder& fcs) const {
694
43018
  Assert(!fcs.underConstruction());
695
86036
  TimerStat::CodeTimer codeTimer(d_statistics.d_weakenTime);
696
697
86036
  Debug("arith::weak") << "LinearEqualityModule::minimallyWeakConflict("
698
43018
                       << aboveUpper <<", "<< basicVar << ", ...) start" << endl;
699
700
43018
  const Rational& adjustSgn = aboveUpper ? d_negOne : d_one;
701
43018
  const DeltaRational& assignment = d_variables.getAssignment(basicVar);
702
86036
  DeltaRational surplus;
703
43018
  if(aboveUpper){
704
16988
    Assert(d_variables.hasUpperBound(basicVar));
705
16988
    Assert(assignment > d_variables.getUpperBound(basicVar));
706
16988
    surplus = assignment - d_variables.getUpperBound(basicVar);
707
  }else{
708
26030
    Assert(d_variables.hasLowerBound(basicVar));
709
26030
    Assert(assignment < d_variables.getLowerBound(basicVar));
710
26030
    surplus = d_variables.getLowerBound(basicVar) - assignment;
711
  }
712
713
43018
  bool anyWeakenings = false;
714
808914
  for(Tableau::RowIterator i = d_tableau.basicRowIterator(basicVar); !i.atEnd(); ++i){
715
765896
    const Tableau::Entry& entry = *i;
716
765896
    ArithVar v = entry.getColVar();
717
765896
    const Rational& coeff = entry.getCoefficient();
718
765896
    bool weakening = false;
719
765896
    ConstraintP c = weakestExplanation(aboveUpper, surplus, v, coeff, weakening, basicVar);
720
1531792
    Debug("arith::weak") << "weak : " << weakening << " "
721
1531792
                         << c->assertedToTheTheory() << " "
722
765896
                         << d_variables.getAssignment(v) << " "
723
765896
                         << c << endl;
724
765896
    anyWeakenings = anyWeakenings || weakening;
725
726
765896
    fcs.addConstraint(c, coeff, adjustSgn);
727
765896
    if(basicVar == v){
728
43018
      Assert(!c->negationHasProof());
729
43018
      fcs.makeLastConsequent();
730
    }
731
  }
732
43018
  Assert(fcs.consequentIsSet());
733
734
43018
  ConstraintCP conflicted = fcs.commitConflict();
735
736
43018
  ++d_statistics.d_weakeningAttempts;
737
43018
  if(anyWeakenings){
738
33405
    ++d_statistics.d_weakeningSuccesses;
739
  }
740
86036
  Debug("arith::weak") << "LinearEqualityModule::minimallyWeakConflict("
741
43018
                       << aboveUpper <<", "<< basicVar << ", ...) done" << endl;
742
86036
  return conflicted;
743
}
744
745
122986
ArithVar LinearEqualityModule::minVarOrder(ArithVar x, ArithVar y) const {
746
122986
  Assert(x != ARITHVAR_SENTINEL);
747
122986
  Assert(y != ARITHVAR_SENTINEL);
748
122986
  if(x <= y){
749
76095
    return x;
750
  } else {
751
46891
    return y;
752
  }
753
}
754
755
747309
ArithVar LinearEqualityModule::minColLength(ArithVar x, ArithVar y) const {
756
747309
  Assert(x != ARITHVAR_SENTINEL);
757
747309
  Assert(y != ARITHVAR_SENTINEL);
758
747309
  Assert(!d_tableau.isBasic(x));
759
747309
  Assert(!d_tableau.isBasic(y));
760
747309
  uint32_t xLen = d_tableau.getColLength(x);
761
747309
  uint32_t yLen = d_tableau.getColLength(y);
762
747309
  if( xLen > yLen){
763
92089
     return y;
764
655220
  } else if( xLen== yLen ){
765
122478
    return minVarOrder(x,y);
766
  }else{
767
532742
    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
772861
ArithVar LinearEqualityModule::minBoundAndColLength(ArithVar x, ArithVar y) const{
788
772861
  Assert(x != ARITHVAR_SENTINEL);
789
772861
  Assert(y != ARITHVAR_SENTINEL);
790
772861
  Assert(!d_tableau.isBasic(x));
791
772861
  Assert(!d_tableau.isBasic(y));
792
772861
  if(d_variables.hasEitherBound(x) && !d_variables.hasEitherBound(y)){
793
8284
    return y;
794
764577
  }else if(!d_variables.hasEitherBound(x) && d_variables.hasEitherBound(y)){
795
17268
    return x;
796
  }else {
797
747309
    return minColLength(x, y);
798
  }
799
}
800
801
template <bool above>
802
136052
ArithVar LinearEqualityModule::selectSlack(ArithVar x_i, VarPreferenceFunction pref) const{
803
136052
  ArithVar slack = ARITHVAR_SENTINEL;
804
805
4585933
  for(Tableau::RowIterator iter = d_tableau.basicRowIterator(x_i); !iter.atEnd();  ++iter){
806
4449881
    const Tableau::Entry& entry = *iter;
807
4449881
    ArithVar nonbasic = entry.getColVar();
808
4449881
    if(nonbasic == x_i) continue;
809
810
4313829
    const Rational& a_ij = entry.getCoefficient();
811
4313829
    int sgn = a_ij.sgn();
812
4313829
    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
909421
      slack = (slack == ARITHVAR_SENTINEL) ? nonbasic : (this->*pref)(slack, nonbasic);
817
    }
818
  }
819
820
136052
  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
1202999
void LinearEqualityModule::startTrackingBoundCounts(){
845
1202999
  Assert(!d_areTracking);
846
1202999
  d_areTracking = true;
847
1202999
  if(Debug.isOn("arith::tracking")){
848
    debugCheckTracking();
849
  }
850
1202999
  Assert(d_areTracking);
851
1202999
}
852
853
1202999
void LinearEqualityModule::stopTrackingBoundCounts(){
854
1202999
  Assert(d_areTracking);
855
1202999
  d_areTracking = false;
856
1202999
  if(Debug.isOn("arith::tracking")){
857
    debugCheckTracking();
858
  }
859
1202999
  Assert(!d_areTracking);
860
1202999
}
861
862
863
51832
void LinearEqualityModule::trackRowIndex(RowIndex ridx){
864
51832
  Assert(!rowIndexIsTracked(ridx));
865
51832
  BoundsInfo bi = computeRowBoundInfo(ridx, true);
866
51832
  d_btracking.set(ridx, bi);
867
51832
}
868
869
51832
BoundsInfo LinearEqualityModule::computeRowBoundInfo(RowIndex ridx, bool inQueue) const{
870
51832
  BoundsInfo bi;
871
872
51832
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
873
590868
  for(; !iter.atEnd();  ++iter){
874
269518
    const Tableau::Entry& entry = *iter;
875
269518
    ArithVar v = entry.getColVar();
876
269518
    const Rational& a_ij = entry.getCoefficient();
877
269518
    bi += (d_variables.selectBoundsInfo(v, inQueue)).multiplyBySgn(a_ij.sgn());
878
  }
879
51832
  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
257702
bool LinearEqualityModule::nonbasicsAtLowerBounds(ArithVar basic) const {
932
257702
  Assert(basicIsTracked(basic));
933
257702
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
934
935
257702
  BoundCounts bcs = d_btracking[ridx].atBounds();
936
257702
  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
257702
  uint32_t lbc = bcs.lowerBoundCount();
946
600344
  return  (lbc == length) ||
947
600344
    (lbc + 1 == length && d_variables.cmpAssignmentUpperBound(basic) != 0);
948
}
949
950
318585
bool LinearEqualityModule::nonbasicsAtUpperBounds(ArithVar basic) const {
951
318585
  Assert(basicIsTracked(basic));
952
318585
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
953
318585
  BoundCounts bcs = d_btracking[ridx].atBounds();
954
318585
  uint32_t length = d_tableau.basicRowLength(basic);
955
318585
  uint32_t ubc = bcs.upperBoundCount();
956
  // See the comment for nonbasicsAtLowerBounds()
957
958
767320
  return (ubc == length) ||
959
767320
    (ubc + 1 == length && d_variables.cmpAssignmentLowerBound(basic) != 0);
960
}
961
962
136052
void LinearEqualityModule::trackingMultiplyRow(RowIndex ridx, int sgn) {
963
136052
  Assert(rowIndexIsTracked(ridx));
964
136052
  Assert(sgn != 0);
965
136052
  if(sgn < 0){
966
77560
    BoundsInfo& bi = d_btracking.get(ridx);
967
77560
    bi = bi.multiplyBySgn(sgn);
968
  }
969
136052
}
970
971
44467491
void LinearEqualityModule::trackingCoefficientChange(RowIndex ridx, ArithVar nb, int oldSgn, int currSgn){
972
44467491
  Assert(oldSgn != currSgn);
973
44467491
  BoundsInfo nb_inf = d_variables.boundsInfo(nb);
974
975
44467491
  Assert(rowIndexIsTracked(ridx));
976
977
44467491
  BoundsInfo& row_bi = d_btracking.get(ridx);
978
44467491
  row_bi.addInSgn(nb_inf, oldSgn, currSgn);
979
44467491
}
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
22755
}  // namespace cvc5