GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/linear_equality.cpp Lines: 376 814 46.2 %
Date: 2021-03-22 Branches: 642 3352 19.2 %

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