GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/linear_equality.cpp Lines: 368 804 45.8 %
Date: 2021-09-16 Branches: 621 3298 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
9942
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
9942
  d_trackCallback(this)
66
9942
{}
67
68
9942
LinearEqualityModule::Statistics::Statistics()
69
    : d_statPivots(
70
19884
        smtStatisticsRegistry().registerInt("theory::arith::pivots")),
71
      d_statUpdates(
72
19884
          smtStatisticsRegistry().registerInt("theory::arith::updates")),
73
      d_pivotTime(
74
19884
          smtStatisticsRegistry().registerTimer("theory::arith::pivotTime")),
75
      d_adjTime(
76
19884
          smtStatisticsRegistry().registerTimer("theory::arith::adjTime")),
77
9942
      d_weakeningAttempts(smtStatisticsRegistry().registerInt(
78
19884
          "theory::arith::weakening::attempts")),
79
9942
      d_weakeningSuccesses(smtStatisticsRegistry().registerInt(
80
19884
          "theory::arith::weakening::success")),
81
9942
      d_weakenings(smtStatisticsRegistry().registerInt(
82
19884
          "theory::arith::weakening::total")),
83
9942
      d_weakenTime(smtStatisticsRegistry().registerTimer(
84
19884
          "theory::arith::weakening::time")),
85
      d_forceTime(
86
89478
          smtStatisticsRegistry().registerTimer("theory::arith::forcing::time"))
87
{
88
9942
}
89
90
4882794
void LinearEqualityModule::includeBoundUpdate(ArithVar v, const BoundsInfo& prev){
91
4882794
  Assert(!d_areTracking);
92
93
4882794
  BoundsInfo curr = d_variables.boundsInfo(v);
94
95
4882794
  Assert(prev != curr);
96
4882794
  Tableau::ColIterator basicIter = d_tableau.colIterator(v);
97
52483242
  for(; !basicIter.atEnd(); ++basicIter){
98
23800224
    const Tableau::Entry& entry = *basicIter;
99
23800224
    Assert(entry.getColVar() == v);
100
23800224
    int a_ijSgn = entry.getCoefficient().sgn();
101
102
23800224
    RowIndex ridx = entry.getRowIndex();
103
23800224
    BoundsInfo& counts = d_btracking.get(ridx);
104
23800224
    Debug("includeBoundUpdate") << d_tableau.rowIndexToBasic(ridx) << " " << counts << " to " ;
105
23800224
    counts.addInChange(a_ijSgn, prev, curr);
106
23800224
    Debug("includeBoundUpdate") << counts << " " << a_ijSgn << std::endl;
107
  }
108
4882794
}
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
166034
void LinearEqualityModule::updateUntracked(ArithVar x_i, const DeltaRational& v){
181
166034
  Assert(!d_tableau.isBasic(x_i));
182
166034
  Assert(!d_areTracking);
183
166034
  const DeltaRational& assignment_x_i = d_variables.getAssignment(x_i);
184
166034
  ++(d_statistics.d_statUpdates);
185
186
187
332068
  Debug("arith") <<"update " << x_i << ": "
188
166034
                 << assignment_x_i << "|-> " << v << endl;
189
332068
  DeltaRational diff = v - assignment_x_i;
190
191
166034
  Tableau::ColIterator colIter = d_tableau.colIterator(x_i);
192
3856154
  for(; !colIter.atEnd(); ++colIter){
193
1845060
    const Tableau::Entry& entry = *colIter;
194
1845060
    Assert(entry.getColVar() == x_i);
195
196
1845060
    ArithVar x_j = d_tableau.rowIndexToBasic(entry.getRowIndex());
197
1845060
    const Rational& a_ji = entry.getCoefficient();
198
199
1845060
    const DeltaRational& assignment = d_variables.getAssignment(x_j);
200
3690120
    DeltaRational  nAssignment = assignment+(diff * a_ji);
201
1845060
    d_variables.setAssignment(x_j, nAssignment);
202
203
1845060
    d_basicVariableUpdates(x_j);
204
  }
205
206
166034
  d_variables.setAssignment(x_i, v);
207
208
166034
  if(Debug.isOn("paranoid:check_tableau")){  debugCheckTableau(); }
209
166034
}
210
211
211187
void LinearEqualityModule::updateTracked(ArithVar x_i, const DeltaRational& v){
212
422374
  TimerStat::CodeTimer codeTimer(d_statistics.d_adjTime);
213
214
211187
  Assert(!d_tableau.isBasic(x_i));
215
211187
  Assert(d_areTracking);
216
217
211187
  ++(d_statistics.d_statUpdates);
218
219
422374
  DeltaRational diff =  v - d_variables.getAssignment(x_i);
220
422374
  Debug("arith") <<"update " << x_i << ": "
221
211187
                 << d_variables.getAssignment(x_i) << "|-> " << v << endl;
222
223
224
211187
  BoundCounts before = d_variables.atBoundCounts(x_i);
225
211187
  d_variables.setAssignment(x_i, v);
226
211187
  BoundCounts after = d_variables.atBoundCounts(x_i);
227
228
211187
  bool anyChange = before != after;
229
230
211187
  Tableau::ColIterator colIter = d_tableau.colIterator(x_i);
231
8945703
  for(; !colIter.atEnd(); ++colIter){
232
4367258
    const Tableau::Entry& entry = *colIter;
233
4367258
    Assert(entry.getColVar() == x_i);
234
235
4367258
    RowIndex ridx = entry.getRowIndex();
236
4367258
    ArithVar x_j = d_tableau.rowIndexToBasic(ridx);
237
4367258
    const Rational& a_ji = entry.getCoefficient();
238
239
4367258
    const DeltaRational& assignment = d_variables.getAssignment(x_j);
240
8734516
    DeltaRational  nAssignment = assignment+(diff * a_ji);
241
4367258
    Debug("update") << x_j << " " << a_ji << assignment << " -> " << nAssignment << endl;
242
4367258
    BoundCounts xjBefore = d_variables.atBoundCounts(x_j);
243
4367258
    d_variables.setAssignment(x_j, nAssignment);
244
4367258
    BoundCounts xjAfter = d_variables.atBoundCounts(x_j);
245
246
4367258
    Assert(rowIndexIsTracked(ridx));
247
4367258
    BoundsInfo& next_bc_k = d_btracking.get(ridx);
248
4367258
    if(anyChange){
249
3201510
      next_bc_k.addInAtBoundChange(a_ji.sgn(), before, after);
250
    }
251
4367258
    if(xjBefore != xjAfter){
252
480651
      next_bc_k.addInAtBoundChange(-1, xjBefore, xjAfter);
253
    }
254
255
4367258
    d_basicVariableUpdates(x_j);
256
  }
257
258
211187
  if(Debug.isOn("paranoid:check_tableau")){  debugCheckTableau(); }
259
211187
}
260
261
211187
void LinearEqualityModule::pivotAndUpdate(ArithVar x_i, ArithVar x_j, const DeltaRational& x_i_value){
262
211187
  Assert(x_i != x_j);
263
264
422374
  TimerStat::CodeTimer codeTimer(d_statistics.d_pivotTime);
265
266
  static int instance = 0;
267
268
211187
  if(Debug.isOn("arith::tracking::pre")){
269
    ++instance;
270
    Debug("arith::tracking")  << "pre update #" << instance << endl;
271
    debugCheckTracking();
272
  }
273
274
211187
  if(Debug.isOn("arith::simplex:row")){ debugPivot(x_i, x_j); }
275
276
211187
  RowIndex ridx = d_tableau.basicToRowIndex(x_i);
277
211187
  const Tableau::Entry& entry_ij =  d_tableau.findEntry(ridx, x_j);
278
211187
  Assert(!entry_ij.blank());
279
280
211187
  const Rational& a_ij = entry_ij.getCoefficient();
281
211187
  const DeltaRational& betaX_i = d_variables.getAssignment(x_i);
282
422374
  DeltaRational theta = (x_i_value - betaX_i)/a_ij;
283
422374
  DeltaRational x_j_value = d_variables.getAssignment(x_j) + theta;
284
285
211187
  updateTracked(x_j, x_j_value);
286
287
211187
  if(Debug.isOn("arith::tracking::mid")){
288
    Debug("arith::tracking")  << "postupdate prepivot #" << instance << endl;
289
    debugCheckTracking();
290
  }
291
292
  // Pivots
293
211187
  ++(d_statistics.d_statPivots);
294
295
211187
  d_tableau.pivot(x_i, x_j, d_trackCallback);
296
297
211187
  if(Debug.isOn("arith::tracking::post")){
298
    Debug("arith::tracking")  << "postpivot #" << instance << endl;
299
    debugCheckTracking();
300
  }
301
302
211187
  d_basicVariableUpdates(x_j);
303
304
211187
  if(Debug.isOn("matrix")){
305
    d_tableau.printMatrix();
306
  }
307
211187
}
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
420555
DeltaRational LinearEqualityModule::computeRowBound(RowIndex ridx, bool rowUb, ArithVar skip) const {
427
420555
  DeltaRational sum(0,0);
428
2931735
  for(Tableau::RowIterator i = d_tableau.ridRowIterator(ridx); !i.atEnd(); ++i){
429
2511180
    const Tableau::Entry& entry = (*i);
430
2511180
    ArithVar v = entry.getColVar();
431
2511180
    if(v == skip){ continue; }
432
433
2184107
    const Rational& coeff =  entry.getCoefficient();
434
2184107
    bool vUb = (rowUb == (coeff.sgn() > 0));
435
436
4368214
    const DeltaRational& bound = vUb ?
437
1095902
      d_variables.getUpperBound(v):
438
3272312
      d_variables.getLowerBound(v);
439
440
4368214
    DeltaRational diff = bound * coeff;
441
2184107
    sum = sum + diff;
442
  }
443
420555
  return sum;
444
}
445
446
/**
447
 * Computes the value of a basic variable using the current assignment.
448
 */
449
189972
DeltaRational LinearEqualityModule::computeRowValue(ArithVar x, bool useSafe) const{
450
189972
  Assert(d_tableau.isBasic(x));
451
189972
  DeltaRational sum(0);
452
453
1123948
  for(Tableau::RowIterator i = d_tableau.basicRowIterator(x); !i.atEnd(); ++i){
454
933976
    const Tableau::Entry& entry = (*i);
455
933976
    ArithVar nonbasic = entry.getColVar();
456
933976
    if(nonbasic == x) continue;
457
744004
    const Rational& coeff = entry.getCoefficient();
458
459
744004
    const DeltaRational& assignment = d_variables.getAssignment(nonbasic, useSafe);
460
744004
    sum = sum + (assignment * coeff);
461
  }
462
189972
  return sum;
463
}
464
465
4577852
const Tableau::Entry* LinearEqualityModule::rowLacksBound(RowIndex ridx, bool rowUb, ArithVar skip){
466
4577852
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
467
41499210
  for(; !iter.atEnd(); ++iter){
468
23036785
    const Tableau::Entry& entry = *iter;
469
470
23036785
    ArithVar var = entry.getColVar();
471
23036785
    if(var == skip) { continue; }
472
473
23034093
    int sgn = entry.getCoefficient().sgn();
474
23034093
    bool selectUb = (rowUb == (sgn > 0));
475
46068186
    bool hasBound = selectUb ?
476
11679317
      d_variables.hasUpperBound(var):
477
34388869
      d_variables.hasLowerBound(var);
478
23034093
    if(!hasBound){
479
4576106
      return &entry;
480
    }
481
  }
482
1746
  return NULL;
483
}
484
485
232
void LinearEqualityModule::propagateBasicFromRow(ConstraintP c){
486
232
  Assert(c != NullConstraint);
487
232
  Assert(c->isUpperBound() || c->isLowerBound());
488
232
  Assert(!c->assertedToTheTheory());
489
232
  Assert(!c->hasProof());
490
491
232
  bool upperBound = c->isUpperBound();
492
232
  ArithVar basic = c->getVariable();
493
232
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
494
495
464
  ConstraintCPVec bounds;
496
232
  RationalVectorP coeffs = ARITH_NULLPROOF(new RationalVector());
497
232
  propagateRow(bounds, ridx, upperBound, c, coeffs);
498
232
  c->impliedByFarkas(bounds, coeffs, false);
499
232
  c->tryToPropagate();
500
501
232
  if(coeffs != RationalVectorPSentinel) { delete coeffs; }
502
232
}
503
504
/* An explanation of the farkas coefficients.
505
 *
506
 * We are proving c using the other variables on the row.
507
 * The proof is in terms of the other constraints and the negation of c, ~c.
508
 *
509
 * A row has the form:
510
 *   sum a_i * x_i  = 0
511
 * or
512
 *   sx + sum r y + sum q z = 0
513
 * where r > 0 and q < 0.
514
 *
515
 * If rowUp, we are proving c
516
 *   g = sum r u_y + sum q l_z
517
 *   and c is entailed by -sx <= g
518
 * If !rowUp, we are proving c
519
 *   g = sum r l_y + sum q u_z
520
 *   and c is entailed by -sx >= g
521
 *
522
 *             | s     | c         | ~c       | u_i     | l_i
523
 *   if  rowUp | s > 0 | x >= -g/s | x < -g/s | a_i > 0 | a_i < 0
524
 *   if  rowUp | s < 0 | x <= -g/s | x > -g/s | a_i > 0 | a_i < 0
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
 *
528
 *
529
 * Thus we treat !rowUp as multiplying the row by -1 and rowUp as 1
530
 * for the entire row.
531
 */
532
149744
void LinearEqualityModule::propagateRow(ConstraintCPVec& into, RowIndex ridx, bool rowUp, ConstraintP c, RationalVectorP farkas){
533
149744
  Assert(!c->assertedToTheTheory());
534
149744
  Assert(c->canBePropagated());
535
149744
  Assert(!c->hasProof());
536
537
149744
  if(farkas != RationalVectorPSentinel){
538
70856
    Assert(farkas->empty());
539
70856
    farkas->push_back(Rational(0));
540
  }
541
542
149744
  ArithVar v = c->getVariable();
543
299488
  Debug("arith::propagateRow") << "LinearEqualityModule::propagateRow("
544
149744
                                   << ridx << ", " << rowUp << ", " << v << ") start" << endl;
545
546
149744
  const Rational& multiple = rowUp ? d_one : d_negOne;
547
548
149744
  Debug("arith::propagateRow") << "multiple: " << multiple << endl;
549
550
149744
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
551
2228336
  for(; !iter.atEnd(); ++iter){
552
1039296
    const Tableau::Entry& entry = *iter;
553
1039296
    ArithVar nonbasic = entry.getColVar();
554
1039296
    const Rational& a_ij = entry.getCoefficient();
555
1039296
    int sgn = a_ij.sgn();
556
1039296
    Assert(sgn != 0);
557
1039296
    bool selectUb = rowUp ? (sgn > 0) : (sgn < 0);
558
559
1039296
    Assert(nonbasic != v || (rowUp && a_ij.sgn() > 0 && c->isLowerBound())
560
           || (rowUp && a_ij.sgn() < 0 && c->isUpperBound())
561
           || (!rowUp && a_ij.sgn() > 0 && c->isUpperBound())
562
           || (!rowUp && a_ij.sgn() < 0 && c->isLowerBound()));
563
564
1039296
    if(Debug.isOn("arith::propagateRow")){
565
      if(nonbasic == v){
566
        Debug("arith::propagateRow") << "(target) "
567
                                     << rowUp << " "
568
                                     << a_ij.sgn() << " "
569
                                     << c->isLowerBound() << " "
570
                                     << c->isUpperBound() << endl;
571
572
        Debug("arith::propagateRow") << "(target) ";
573
      }
574
      Debug("arith::propagateRow") << "propagateRow " << a_ij << " * " << nonbasic ;
575
    }
576
577
1039296
    if(nonbasic == v){
578
149744
      if(farkas != RationalVectorPSentinel){
579
70856
        Assert(farkas->front().isZero());
580
141712
        Rational multAij = multiple * a_ij;
581
70856
        Debug("arith::propagateRow") << "(" << multAij << ") ";
582
70856
        farkas->front() = multAij;
583
      }
584
585
149744
      Debug("arith::propagateRow") << c << endl;
586
    }else{
587
588
      ConstraintCP bound = selectUb
589
405537
        ? d_variables.getUpperBoundConstraint(nonbasic)
590
1295089
        : d_variables.getLowerBoundConstraint(nonbasic);
591
592
889552
      if(farkas != RationalVectorPSentinel){
593
811570
        Rational multAij = multiple * a_ij;
594
405785
        Debug("arith::propagateRow") << "(" << multAij << ") ";
595
405785
        farkas->push_back(multAij);
596
      }
597
889552
      Assert(bound != NullConstraint);
598
889552
      Debug("arith::propagateRow") << bound << endl;
599
889552
      into.push_back(bound);
600
    }
601
  }
602
299488
  Debug("arith::propagateRow") << "LinearEqualityModule::propagateRow("
603
149744
                                   << ridx << ", " << rowUp << ", " << v << ") done" << endl;
604
605
149744
}
606
607
973182
ConstraintP LinearEqualityModule::weakestExplanation(bool aboveUpper, DeltaRational& surplus, ArithVar v, const Rational& coeff, bool& anyWeakening, ArithVar basic) const {
608
609
973182
  int sgn = coeff.sgn();
610
973182
  bool ub = aboveUpper?(sgn < 0) : (sgn > 0);
611
612
1946364
  ConstraintP c = ub ?
613
292007
    d_variables.getUpperBoundConstraint(v) :
614
1654357
    d_variables.getLowerBoundConstraint(v);
615
616
  bool weakened;
617
1091894
  do{
618
1091894
    const DeltaRational& bound = c->getValue();
619
620
1091894
    weakened = false;
621
622
2183788
    ConstraintP weaker = ub?
623
369846
      c->getStrictlyWeakerUpperBound(true, true):
624
1813942
      c->getStrictlyWeakerLowerBound(true, true);
625
626
1091894
    if(weaker != NullConstraint){
627
268035
      const DeltaRational& weakerBound = weaker->getValue();
628
629
536070
      DeltaRational diff = aboveUpper ? bound - weakerBound : weakerBound - bound;
630
      //if var == basic,
631
      //  if aboveUpper, weakerBound > bound, multiply by -1
632
      //  if !aboveUpper, weakerBound < bound, multiply by -1
633
268035
      diff = diff * coeff;
634
268035
      if(surplus > diff){
635
118712
        ++d_statistics.d_weakenings;
636
118712
        weakened = true;
637
118712
        anyWeakening = true;
638
118712
        surplus = surplus - diff;
639
640
118712
        Debug("arith::weak") << "found:" << endl;
641
118712
        if(v == basic){
642
12391
          Debug("arith::weak") << "  basic: ";
643
        }
644
237424
        Debug("arith::weak") << "  " << surplus << " "<< diff  << endl
645
118712
                             << "  " << bound << c << endl
646
118712
                             << "  " << weakerBound << weaker << endl;
647
648
118712
        Assert(diff.sgn() > 0);
649
118712
        c = weaker;
650
      }
651
    }
652
  }while(weakened);
653
654
973182
  return c;
655
}
656
657
/* An explanation of the farkas coefficients.
658
 *
659
 * We are proving a conflict on the basic variable x_b.
660
 * If aboveUpper, then the conflict is with the constraint c : x_b <= u_b.
661
 * If !aboveUpper, then the conflict is with the constraint c : x_b >= l_b.
662
 *
663
 * A row has the form:
664
 *   -x_b sum a_i * x_i  = 0
665
 * or
666
 *   -x_b + sum r y + sum q z = 0,
667
 *    x_b = sum r y + sum q z
668
 * where r > 0 and q < 0.
669
 *
670
 *
671
 * If !aboveUp, we are proving ~c: x_b < l_b
672
 *   g = sum r u_y + sum q l_z
673
 *   x_b <= g < l_b
674
 *   and ~c is entailed by x_b <= g
675
 *
676
 * If aboveUp, we are proving ~c : x_b > u_b
677
 *   g = sum r l_y + sum q u_z
678
 *   x_b >= g > u_b
679
 *   and ~c is entailed by x_b >= g
680
 *
681
 *
682
 *               | s     | c         | ~c       | u_i     | l_i
683
 *   if !aboveUp | s > 0 | x >= -g/s | x < -g/s | a_i > 0 | a_i < 0
684
 *   if !aboveUp | s < 0 | x <= -g/s | x > -g/s | a_i > 0 | a_i < 0
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
 *
688
 * Thus we treat aboveUp as multiplying the row by -1 and !aboveUp as 1
689
 * for the entire row.
690
 */
691
71173
ConstraintCP LinearEqualityModule::minimallyWeakConflict(bool aboveUpper, ArithVar basicVar, FarkasConflictBuilder& fcs) const {
692
71173
  Assert(!fcs.underConstruction());
693
142346
  TimerStat::CodeTimer codeTimer(d_statistics.d_weakenTime);
694
695
142346
  Debug("arith::weak") << "LinearEqualityModule::minimallyWeakConflict("
696
71173
                       << aboveUpper <<", "<< basicVar << ", ...) start" << endl;
697
698
71173
  const Rational& adjustSgn = aboveUpper ? d_negOne : d_one;
699
71173
  const DeltaRational& assignment = d_variables.getAssignment(basicVar);
700
142346
  DeltaRational surplus;
701
71173
  if(aboveUpper){
702
31495
    Assert(d_variables.hasUpperBound(basicVar));
703
31495
    Assert(assignment > d_variables.getUpperBound(basicVar));
704
31495
    surplus = assignment - d_variables.getUpperBound(basicVar);
705
  }else{
706
39678
    Assert(d_variables.hasLowerBound(basicVar));
707
39678
    Assert(assignment < d_variables.getLowerBound(basicVar));
708
39678
    surplus = d_variables.getLowerBound(basicVar) - assignment;
709
  }
710
711
71173
  bool anyWeakenings = false;
712
1044355
  for(Tableau::RowIterator i = d_tableau.basicRowIterator(basicVar); !i.atEnd(); ++i){
713
973182
    const Tableau::Entry& entry = *i;
714
973182
    ArithVar v = entry.getColVar();
715
973182
    const Rational& coeff = entry.getCoefficient();
716
973182
    bool weakening = false;
717
973182
    ConstraintP c = weakestExplanation(aboveUpper, surplus, v, coeff, weakening, basicVar);
718
1946364
    Debug("arith::weak") << "weak : " << weakening << " "
719
1946364
                         << c->assertedToTheTheory() << " "
720
973182
                         << d_variables.getAssignment(v) << " "
721
973182
                         << c << endl;
722
973182
    anyWeakenings = anyWeakenings || weakening;
723
724
973182
    fcs.addConstraint(c, coeff, adjustSgn);
725
973182
    if(basicVar == v){
726
71173
      Assert(!c->negationHasProof());
727
71173
      fcs.makeLastConsequent();
728
    }
729
  }
730
71173
  Assert(fcs.consequentIsSet());
731
732
71173
  ConstraintCP conflicted = fcs.commitConflict();
733
734
71173
  ++d_statistics.d_weakeningAttempts;
735
71173
  if(anyWeakenings){
736
53559
    ++d_statistics.d_weakeningSuccesses;
737
  }
738
142346
  Debug("arith::weak") << "LinearEqualityModule::minimallyWeakConflict("
739
71173
                       << aboveUpper <<", "<< basicVar << ", ...) done" << endl;
740
142346
  return conflicted;
741
}
742
743
131314
ArithVar LinearEqualityModule::minVarOrder(ArithVar x, ArithVar y) const {
744
131314
  Assert(x != ARITHVAR_SENTINEL);
745
131314
  Assert(y != ARITHVAR_SENTINEL);
746
131314
  if(x <= y){
747
79759
    return x;
748
  } else {
749
51555
    return y;
750
  }
751
}
752
753
818292
ArithVar LinearEqualityModule::minColLength(ArithVar x, ArithVar y) const {
754
818292
  Assert(x != ARITHVAR_SENTINEL);
755
818292
  Assert(y != ARITHVAR_SENTINEL);
756
818292
  Assert(!d_tableau.isBasic(x));
757
818292
  Assert(!d_tableau.isBasic(y));
758
818292
  uint32_t xLen = d_tableau.getColLength(x);
759
818292
  uint32_t yLen = d_tableau.getColLength(y);
760
818292
  if( xLen > yLen){
761
118299
     return y;
762
699993
  } else if( xLen== yLen ){
763
130797
    return minVarOrder(x,y);
764
  }else{
765
569196
    return x;
766
  }
767
}
768
769
ArithVar LinearEqualityModule::minRowLength(ArithVar x, ArithVar y) const {
770
  Assert(x != ARITHVAR_SENTINEL);
771
  Assert(y != ARITHVAR_SENTINEL);
772
  Assert(d_tableau.isBasic(x));
773
  Assert(d_tableau.isBasic(y));
774
  uint32_t xLen = d_tableau.basicRowLength(x);
775
  uint32_t yLen = d_tableau.basicRowLength(y);
776
  if( xLen > yLen){
777
     return y;
778
  } else if( xLen== yLen ){
779
    return minVarOrder(x,y);
780
  }else{
781
    return x;
782
  }
783
}
784
785
861762
ArithVar LinearEqualityModule::minBoundAndColLength(ArithVar x, ArithVar y) const{
786
861762
  Assert(x != ARITHVAR_SENTINEL);
787
861762
  Assert(y != ARITHVAR_SENTINEL);
788
861762
  Assert(!d_tableau.isBasic(x));
789
861762
  Assert(!d_tableau.isBasic(y));
790
861762
  if(d_variables.hasEitherBound(x) && !d_variables.hasEitherBound(y)){
791
15789
    return y;
792
845973
  }else if(!d_variables.hasEitherBound(x) && d_variables.hasEitherBound(y)){
793
27681
    return x;
794
  }else {
795
818292
    return minColLength(x, y);
796
  }
797
}
798
799
template <bool above>
800
211187
ArithVar LinearEqualityModule::selectSlack(ArithVar x_i, VarPreferenceFunction pref) const{
801
211187
  ArithVar slack = ARITHVAR_SENTINEL;
802
803
5174701
  for(Tableau::RowIterator iter = d_tableau.basicRowIterator(x_i); !iter.atEnd();  ++iter){
804
4963514
    const Tableau::Entry& entry = *iter;
805
4963514
    ArithVar nonbasic = entry.getColVar();
806
4963514
    if(nonbasic == x_i) continue;
807
808
4752327
    const Rational& a_ij = entry.getCoefficient();
809
4752327
    int sgn = a_ij.sgn();
810
4752327
    if(isAcceptableSlack<above>(sgn, nonbasic)){
811
      //If one of the above conditions is met, we have found an acceptable
812
      //nonbasic variable to pivot x_i with.  We can now choose which one we
813
      //prefer the most.
814
1073466
      slack = (slack == ARITHVAR_SENTINEL) ? nonbasic : (this->*pref)(slack, nonbasic);
815
    }
816
  }
817
818
211187
  return slack;
819
}
820
821
const Tableau::Entry* LinearEqualityModule::selectSlackEntry(ArithVar x_i, bool above) const{
822
  for(Tableau::RowIterator iter = d_tableau.basicRowIterator(x_i); !iter.atEnd();  ++iter){
823
    const Tableau::Entry& entry = *iter;
824
    ArithVar nonbasic = entry.getColVar();
825
    if(nonbasic == x_i) continue;
826
827
    const Rational& a_ij = entry.getCoefficient();
828
    int sgn = a_ij.sgn();
829
    if(above && isAcceptableSlack<true>(sgn, nonbasic)){
830
      //If one of the above conditions is met, we have found an acceptable
831
      //nonbasic variable to pivot x_i with.  We can now choose which one we
832
      //prefer the most.
833
      return &entry;
834
    }else if(!above && isAcceptableSlack<false>(sgn, nonbasic)){
835
      return &entry;
836
    }
837
  }
838
839
  return NULL;
840
}
841
842
1837891
void LinearEqualityModule::startTrackingBoundCounts(){
843
1837891
  Assert(!d_areTracking);
844
1837891
  d_areTracking = true;
845
1837891
  if(Debug.isOn("arith::tracking")){
846
    debugCheckTracking();
847
  }
848
1837891
  Assert(d_areTracking);
849
1837891
}
850
851
1837891
void LinearEqualityModule::stopTrackingBoundCounts(){
852
1837891
  Assert(d_areTracking);
853
1837891
  d_areTracking = false;
854
1837891
  if(Debug.isOn("arith::tracking")){
855
    debugCheckTracking();
856
  }
857
1837891
  Assert(!d_areTracking);
858
1837891
}
859
860
861
94986
void LinearEqualityModule::trackRowIndex(RowIndex ridx){
862
94986
  Assert(!rowIndexIsTracked(ridx));
863
94986
  BoundsInfo bi = computeRowBoundInfo(ridx, true);
864
94986
  d_btracking.set(ridx, bi);
865
94986
}
866
867
94986
BoundsInfo LinearEqualityModule::computeRowBoundInfo(RowIndex ridx, bool inQueue) const{
868
94986
  BoundsInfo bi;
869
870
94986
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
871
1028962
  for(; !iter.atEnd();  ++iter){
872
466988
    const Tableau::Entry& entry = *iter;
873
466988
    ArithVar v = entry.getColVar();
874
466988
    const Rational& a_ij = entry.getCoefficient();
875
466988
    bi += (d_variables.selectBoundsInfo(v, inQueue)).multiplyBySgn(a_ij.sgn());
876
  }
877
94986
  return bi;
878
}
879
880
BoundCounts LinearEqualityModule::debugBasicAtBoundCount(ArithVar x_i) const {
881
  return d_btracking[d_tableau.basicToRowIndex(x_i)].atBounds();
882
}
883
884
/**
885
 * If the pivot described in u were performed,
886
 * then the row would qualify as being either at the minimum/maximum
887
 * to the non-basics being at their bounds.
888
 * The minimum/maximum is determined by the direction the non-basic is changing.
889
 */
890
bool LinearEqualityModule::basicsAtBounds(const UpdateInfo& u) const {
891
  Assert(u.describesPivot());
892
893
  ArithVar nonbasic = u.nonbasic();
894
  ArithVar basic = u.leaving();
895
  Assert(basicIsTracked(basic));
896
  int coeffSgn = u.getCoefficient().sgn();
897
  int nbdir = u.nonbasicDirection();
898
899
  ConstraintP c = u.limiting();
900
  int toUB = (c->getType() == UpperBound ||
901
              c->getType() == Equality) ? 1 : 0;
902
  int toLB = (c->getType() == LowerBound ||
903
              c->getType() == Equality) ? 1 : 0;
904
905
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
906
907
  BoundCounts bcs = d_btracking[ridx].atBounds();
908
  // x = c*n + \sum d*m
909
  // 0 = -x + c*n + \sum d*m
910
  // n = 1/c * x + -1/c * (\sum d*m)
911
  BoundCounts nonb = bcs - d_variables.atBoundCounts(nonbasic).multiplyBySgn(coeffSgn);
912
  nonb.addInChange(-1, d_variables.atBoundCounts(basic), BoundCounts(toLB, toUB));
913
  nonb = nonb.multiplyBySgn(-coeffSgn);
914
915
  uint32_t length = d_tableau.basicRowLength(basic);
916
  Debug("basicsAtBounds")
917
    << "bcs " << bcs
918
    << "nonb " << nonb
919
    << "length " << length << endl;
920
  // nonb has nb excluded.
921
  if(nbdir < 0){
922
    return nonb.lowerBoundCount() + 1 == length;
923
  }else{
924
    Assert(nbdir > 0);
925
    return nonb.upperBoundCount() + 1 == length;
926
  }
927
}
928
929
457822
bool LinearEqualityModule::nonbasicsAtLowerBounds(ArithVar basic) const {
930
457822
  Assert(basicIsTracked(basic));
931
457822
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
932
933
457822
  BoundCounts bcs = d_btracking[ridx].atBounds();
934
457822
  uint32_t length = d_tableau.basicRowLength(basic);
935
936
  // return true if excluding the basic is every element is at its "lowerbound"
937
  // The psuedo code is:
938
  //   bcs -= basic.count(basic, basic's sgn)
939
  //   return bcs.lowerBoundCount() + 1 == length
940
  // As basic's sign is always -1, we can pull out the pieces of the count:
941
  //   bcs.lowerBoundCount() - basic.atUpperBoundInd() + 1 == length
942
  // basic.atUpperBoundInd() is either 0 or 1
943
457822
  uint32_t lbc = bcs.lowerBoundCount();
944
1073119
  return  (lbc == length) ||
945
1073119
    (lbc + 1 == length && d_variables.cmpAssignmentUpperBound(basic) != 0);
946
}
947
948
509628
bool LinearEqualityModule::nonbasicsAtUpperBounds(ArithVar basic) const {
949
509628
  Assert(basicIsTracked(basic));
950
509628
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
951
509628
  BoundCounts bcs = d_btracking[ridx].atBounds();
952
509628
  uint32_t length = d_tableau.basicRowLength(basic);
953
509628
  uint32_t ubc = bcs.upperBoundCount();
954
  // See the comment for nonbasicsAtLowerBounds()
955
956
1217646
  return (ubc == length) ||
957
1217646
    (ubc + 1 == length && d_variables.cmpAssignmentLowerBound(basic) != 0);
958
}
959
960
211187
void LinearEqualityModule::trackingMultiplyRow(RowIndex ridx, int sgn) {
961
211187
  Assert(rowIndexIsTracked(ridx));
962
211187
  Assert(sgn != 0);
963
211187
  if(sgn < 0){
964
118410
    BoundsInfo& bi = d_btracking.get(ridx);
965
118410
    bi = bi.multiplyBySgn(sgn);
966
  }
967
211187
}
968
969
56849783
void LinearEqualityModule::trackingCoefficientChange(RowIndex ridx, ArithVar nb, int oldSgn, int currSgn){
970
56849783
  Assert(oldSgn != currSgn);
971
56849783
  BoundsInfo nb_inf = d_variables.boundsInfo(nb);
972
973
56849783
  Assert(rowIndexIsTracked(ridx));
974
975
56849783
  BoundsInfo& row_bi = d_btracking.get(ridx);
976
56849783
  row_bi.addInSgn(nb_inf, oldSgn, currSgn);
977
56849783
}
978
979
ArithVar LinearEqualityModule::minBy(const ArithVarVec& vec, VarPreferenceFunction pf) const{
980
  if(vec.empty()) {
981
    return ARITHVAR_SENTINEL;
982
  }else {
983
    ArithVar sel = vec.front();
984
    ArithVarVec::const_iterator i = vec.begin() + 1;
985
    ArithVarVec::const_iterator i_end = vec.end();
986
    for(; i != i_end; ++i){
987
      sel = (this->*pf)(sel, *i);
988
    }
989
    return sel;
990
  }
991
}
992
993
bool LinearEqualityModule::accumulateBorder(const Tableau::Entry& entry, bool ub){
994
  ArithVar currBasic = d_tableau.rowIndexToBasic(entry.getRowIndex());
995
996
  Assert(basicIsTracked(currBasic));
997
998
  ConstraintP bound = ub ?
999
    d_variables.getUpperBoundConstraint(currBasic):
1000
    d_variables.getLowerBoundConstraint(currBasic);
1001
1002
  if(bound == NullConstraint){ return false; }
1003
  Assert(bound != NullConstraint);
1004
1005
  const Rational& coeff = entry.getCoefficient();
1006
1007
  const DeltaRational& assignment = d_variables.getAssignment(currBasic);
1008
  DeltaRational toBound = bound->getValue() - assignment;
1009
  DeltaRational nbDiff = toBound/coeff;
1010
1011
  // if ub
1012
  // if toUB >= 0
1013
  // then ub >= currBasic
1014
  //   if sgn > 0,
1015
  //   then diff >= 0, so nb must increase for G
1016
  //   else diff <= 0, so nb must decrease for G
1017
  // else ub < currBasic
1018
  //   if sgn > 0,
1019
  //   then diff < 0, so nb must decrease for G
1020
  //   else diff > 0, so nb must increase for G
1021
1022
  int diffSgn = nbDiff.sgn();
1023
1024
  if(diffSgn != 0 && willBeInConflictAfterPivot(entry, nbDiff, ub)){
1025
    return true;
1026
  }else{
1027
    bool areFixing = ub ? (toBound.sgn() < 0 ) : (toBound.sgn() > 0);
1028
    Border border(bound, nbDiff, areFixing, &entry, ub);
1029
    bool increasing =
1030
      (diffSgn > 0) ||
1031
      (diffSgn == 0 && ((coeff.sgn() > 0) == ub));
1032
1033
    // assume diffSgn == 0
1034
    // if coeff > 0,
1035
    //   if ub, inc
1036
    //   else, dec
1037
    // else coeff < 0
1038
    //   if ub, dec
1039
    //   else, inc
1040
1041
    if(increasing){
1042
      Debug("handleBorders") << "push back increasing " << border << endl;
1043
      d_increasing.push_back(border);
1044
    }else{
1045
      Debug("handleBorders") << "push back decreasing " << border << endl;
1046
      d_decreasing.push_back(border);
1047
    }
1048
    return false;
1049
  }
1050
}
1051
1052
bool LinearEqualityModule::willBeInConflictAfterPivot(const Tableau::Entry& entry, const DeltaRational& nbDiff, bool bToUB) const{
1053
  int nbSgn = nbDiff.sgn();
1054
  Assert(nbSgn != 0);
1055
1056
  if(nbSgn > 0){
1057
    if (!d_upperBoundDifference || nbDiff <= *d_upperBoundDifference)
1058
    {
1059
      return false;
1060
    }
1061
  }else{
1062
    if (!d_lowerBoundDifference || nbDiff >= *d_lowerBoundDifference)
1063
    {
1064
      return false;
1065
    }
1066
  }
1067
1068
  // Assume past this point, nb will be in error if this pivot is done
1069
  ArithVar nb = entry.getColVar();
1070
  RowIndex ridx = entry.getRowIndex();
1071
  ArithVar basic = d_tableau.rowIndexToBasic(ridx);
1072
  Assert(rowIndexIsTracked(ridx));
1073
  int coeffSgn = entry.getCoefficient().sgn();
1074
1075
1076
  // if bToUB, then basic is going to be set to its upperbound
1077
  // if not bToUB, then basic is going to be set to its lowerbound
1078
1079
  // Different steps of solving for this:
1080
  // 1) y = a * x + \sum b * z
1081
  // 2) -a * x = -y + \sum b * z
1082
  // 3) x = (-1/a) * ( -y + \sum b * z)
1083
1084
  BoundCounts bc = d_btracking[ridx].atBounds();
1085
1086
  // 1) y = a * x + \sum b * z
1087
  // Get bc(\sum b * z)
1088
  BoundCounts sumOnly = bc - d_variables.atBoundCounts(nb).multiplyBySgn(coeffSgn);
1089
1090
  // y's bounds in the proposed model
1091
  int yWillBeAtUb = (bToUB || d_variables.boundsAreEqual(basic)) ? 1 : 0;
1092
  int yWillBeAtLb = (!bToUB || d_variables.boundsAreEqual(basic)) ? 1 : 0;
1093
  BoundCounts ysBounds(yWillBeAtLb, yWillBeAtUb);
1094
1095
  // 2) -a * x = -y + \sum b * z
1096
  // Get bc(-y + \sum b * z)
1097
  sumOnly.addInChange(-1, d_variables.atBoundCounts(basic), ysBounds);
1098
1099
  // 3) x = (-1/a) * ( -y + \sum b * z)
1100
  // Get bc((-1/a) * ( -y + \sum b * z))
1101
  BoundCounts xsBoundsAfterPivot = sumOnly.multiplyBySgn(-coeffSgn);
1102
1103
  uint32_t length = d_tableau.basicRowLength(basic);
1104
  if(nbSgn > 0){
1105
    // Only check for the upper bound being violated
1106
    return xsBoundsAfterPivot.lowerBoundCount() + 1 == length;
1107
  }else{
1108
    // Only check for the lower bound being violated
1109
    return xsBoundsAfterPivot.upperBoundCount() + 1 == length;
1110
  }
1111
}
1112
1113
UpdateInfo LinearEqualityModule::mkConflictUpdate(const Tableau::Entry& entry, bool ub) const{
1114
  ArithVar currBasic = d_tableau.rowIndexToBasic(entry.getRowIndex());
1115
  ArithVar nb = entry.getColVar();
1116
1117
  ConstraintP bound = ub ?
1118
    d_variables.getUpperBoundConstraint(currBasic):
1119
    d_variables.getLowerBoundConstraint(currBasic);
1120
1121
1122
  const Rational& coeff = entry.getCoefficient();
1123
  const DeltaRational& assignment = d_variables.getAssignment(currBasic);
1124
  DeltaRational toBound = bound->getValue() - assignment;
1125
  DeltaRational nbDiff = toBound/coeff;
1126
1127
  return UpdateInfo::conflict(nb, nbDiff, coeff, bound);
1128
}
1129
1130
UpdateInfo LinearEqualityModule::speculativeUpdate(ArithVar nb, const Rational& focusCoeff, UpdatePreferenceFunction pref){
1131
  Assert(d_increasing.empty());
1132
  Assert(d_decreasing.empty());
1133
  Assert(!d_lowerBoundDifference);
1134
  Assert(!d_upperBoundDifference);
1135
1136
  int focusCoeffSgn = focusCoeff.sgn();
1137
1138
  static int instance = 0;
1139
  ++instance;
1140
  Debug("speculativeUpdate") << "speculativeUpdate " << instance << endl;
1141
  Debug("speculativeUpdate") << "nb " << nb << endl;
1142
  Debug("speculativeUpdate") << "focusCoeff " << focusCoeff << endl;
1143
1144
  if(d_variables.hasUpperBound(nb)){
1145
    ConstraintP ub = d_variables.getUpperBoundConstraint(nb);
1146
    d_upperBoundDifference = ub->getValue() - d_variables.getAssignment(nb);
1147
    Border border(ub, *d_upperBoundDifference, false, NULL, true);
1148
    Debug("handleBorders") << "push back increasing " << border << endl;
1149
    d_increasing.push_back(border);
1150
  }
1151
  if(d_variables.hasLowerBound(nb)){
1152
    ConstraintP lb = d_variables.getLowerBoundConstraint(nb);
1153
    d_lowerBoundDifference = lb->getValue() - d_variables.getAssignment(nb);
1154
    Border border(lb, *d_lowerBoundDifference, false, NULL, false);
1155
    Debug("handleBorders") << "push back decreasing " << border << endl;
1156
    d_decreasing.push_back(border);
1157
  }
1158
1159
  Tableau::ColIterator colIter = d_tableau.colIterator(nb);
1160
  for(; !colIter.atEnd(); ++colIter){
1161
    const Tableau::Entry& entry = *colIter;
1162
    Assert(entry.getColVar() == nb);
1163
1164
    if(accumulateBorder(entry, true)){
1165
      clearSpeculative();
1166
      return mkConflictUpdate(entry, true);
1167
    }
1168
    if(accumulateBorder(entry, false)){
1169
      clearSpeculative();
1170
      return mkConflictUpdate(entry, false);
1171
    }
1172
  }
1173
1174
  UpdateInfo selected;
1175
  BorderHeap& withSgn = focusCoeffSgn > 0 ? d_increasing : d_decreasing;
1176
  BorderHeap& againstSgn = focusCoeffSgn > 0 ? d_decreasing : d_increasing;
1177
1178
  handleBorders(selected, nb, focusCoeff, withSgn, 0, pref);
1179
  int m = 1 - selected.errorsChangeSafe(0);
1180
  handleBorders(selected, nb, focusCoeff, againstSgn, m, pref);
1181
1182
  clearSpeculative();
1183
  return selected;
1184
}
1185
1186
void LinearEqualityModule::clearSpeculative(){
1187
  // clear everything away
1188
  d_increasing.clear();
1189
  d_decreasing.clear();
1190
  d_lowerBoundDifference.reset();
1191
  d_upperBoundDifference.reset();
1192
}
1193
1194
void LinearEqualityModule::handleBorders(UpdateInfo& selected, ArithVar nb, const Rational& focusCoeff, BorderHeap& heap, int minimumFixes, UpdatePreferenceFunction pref){
1195
  Assert(minimumFixes >= 0);
1196
1197
  // The values popped off of the heap
1198
  // should be popped with the values closest to 0
1199
  // being first and larger in absolute value last
1200
1201
1202
  int fixesRemaining = heap.possibleFixes();
1203
1204
  Debug("handleBorders")
1205
    << "handleBorders "
1206
    << "nb " << nb
1207
    << "fc " << focusCoeff
1208
    << "h.e " << heap.empty()
1209
    << "h.dir " << heap.direction()
1210
    << "h.rem " << fixesRemaining
1211
    << "h.0s " << heap.numZeroes()
1212
    << "min " << minimumFixes
1213
    << endl;
1214
1215
  if(heap.empty()){
1216
    // if the heap is empty, return
1217
    return;
1218
  }
1219
1220
  bool zeroesWillDominate = fixesRemaining - heap.numZeroes() < minimumFixes;
1221
1222
  // can the number of fixes ever exceed the minimum?
1223
  // no more than the number of possible fixes can be fixed in total
1224
  // nothing can be fixed before the zeroes are taken care of
1225
  if(minimumFixes > 0 && zeroesWillDominate){
1226
    return;
1227
  }
1228
1229
1230
  int negErrorChange = 0;
1231
  int nbDir = heap.direction();
1232
1233
  // points at the beginning of the heap
1234
  if(zeroesWillDominate){
1235
    heap.dropNonZeroes();
1236
  }
1237
  heap.make_heap();
1238
1239
1240
  // pretend like the previous block had a value of zero.
1241
  // The block that actually has a value of 0 must handle this.
1242
  const DeltaRational zero(0);
1243
  const DeltaRational* prevBlockValue = &zero;
1244
1245
  /** The coefficient changes as the value crosses border. */
1246
  Rational effectiveCoefficient = focusCoeff;
1247
1248
  /* Keeps track of the change to the value of the focus function.*/
1249
  DeltaRational totalFocusChange(0);
1250
1251
  const int focusCoeffSgn = focusCoeff.sgn();
1252
1253
  while(heap.more()  &&
1254
        (fixesRemaining + negErrorChange > minimumFixes ||
1255
         (fixesRemaining + negErrorChange == minimumFixes &&
1256
          effectiveCoefficient.sgn() == focusCoeffSgn))){
1257
    // There are more elements &&
1258
    // we can either fix at least 1 more variable in the error function
1259
    // or we can improve the error function
1260
1261
1262
    int brokenInBlock = 0;
1263
    BorderVec::const_iterator endBlock = heap.end();
1264
1265
    pop_block(heap, brokenInBlock, fixesRemaining, negErrorChange);
1266
1267
    // if endVec == beginVec, block starts there
1268
    // other wise, block starts at endVec
1269
    BorderVec::const_iterator startBlock
1270
      = heap.more() ? heap.end() : heap.begin();
1271
1272
    const DeltaRational& blockValue = (*startBlock).d_diff;
1273
1274
    // if decreasing
1275
    // blockValue < prevBlockValue
1276
    // diff.sgn() = -1
1277
    DeltaRational diff = blockValue - (*prevBlockValue);
1278
    DeltaRational blockChangeToFocus =  diff * effectiveCoefficient;
1279
    totalFocusChange += blockChangeToFocus;
1280
1281
    Debug("handleBorders")
1282
      << "blockValue " << (blockValue)
1283
      << "diff " << diff
1284
      << "blockChangeToFocus " << totalFocusChange
1285
      << "blockChangeToFocus " << totalFocusChange
1286
      << "negErrorChange " << negErrorChange
1287
      << "brokenInBlock " << brokenInBlock
1288
      << "fixesRemaining " << fixesRemaining
1289
      << endl;
1290
1291
    int currFocusChangeSgn = totalFocusChange.sgn();
1292
    for(BorderVec::const_iterator i = startBlock; i != endBlock; ++i){
1293
      const Border& b = *i;
1294
1295
      Debug("handleBorders") << b << endl;
1296
1297
      bool makesImprovement = negErrorChange > 0 ||
1298
        (negErrorChange == 0  && currFocusChangeSgn > 0);
1299
1300
      if(!makesImprovement){
1301
        if(b.ownBorder() || minimumFixes > 0){
1302
          continue;
1303
        }
1304
      }
1305
1306
      UpdateInfo proposal(nb, nbDir);
1307
      if(b.ownBorder()){
1308
        proposal.witnessedUpdate(b.d_diff, b.d_bound, -negErrorChange, currFocusChangeSgn);
1309
      }else{
1310
        proposal.update(b.d_diff, b.getCoefficient(), b.d_bound, -negErrorChange, currFocusChangeSgn);
1311
      }
1312
1313
      if(selected.unbounded() || (this->*pref)(selected, proposal)){
1314
        selected = proposal;
1315
      }
1316
    }
1317
1318
    effectiveCoefficient += updateCoefficient(startBlock, endBlock);
1319
    prevBlockValue = &blockValue;
1320
    negErrorChange -= brokenInBlock;
1321
  }
1322
}
1323
1324
Rational LinearEqualityModule::updateCoefficient(BorderVec::const_iterator startBlock, BorderVec::const_iterator endBlock){
1325
  //update coefficient
1326
  Rational changeToCoefficient(0);
1327
  for(BorderVec::const_iterator i = startBlock; i != endBlock; ++i){
1328
    const Border& curr = *i;
1329
    if(curr.ownBorder()){// breaking its own bound
1330
      if(curr.d_upperbound){
1331
        changeToCoefficient -= 1;
1332
      }else{
1333
        changeToCoefficient += 1;
1334
      }
1335
    }else{
1336
      const Rational& coeff = curr.d_entry->getCoefficient();
1337
      if(curr.d_areFixing){
1338
        if(curr.d_upperbound){// fixing an upper bound
1339
          changeToCoefficient += coeff;
1340
        }else{// fixing a lower bound
1341
          changeToCoefficient -= coeff;
1342
        }
1343
      }else{
1344
        if(curr.d_upperbound){// breaking an upper bound
1345
          changeToCoefficient -= coeff;
1346
        }else{
1347
          // breaking a lower bound
1348
          changeToCoefficient += coeff;
1349
        }
1350
      }
1351
    }
1352
  }
1353
  return changeToCoefficient;
1354
}
1355
1356
void LinearEqualityModule::pop_block(BorderHeap& heap, int& brokenInBlock, int& fixesRemaining, int& negErrorChange){
1357
  Assert(heap.more());
1358
1359
  if(heap.top().d_areFixing){
1360
    fixesRemaining--;
1361
    negErrorChange++;
1362
  }else{
1363
    brokenInBlock++;
1364
  }
1365
  heap.pop_heap();
1366
  const DeltaRational& blockValue = (*heap.end()).d_diff;
1367
1368
  while(heap.more()){
1369
    const Border& top = heap.top();
1370
    if(blockValue == top.d_diff){
1371
      // belongs to the block
1372
      if(top.d_areFixing){
1373
        fixesRemaining--;
1374
        negErrorChange++;
1375
      }else{
1376
        brokenInBlock++;
1377
      }
1378
      heap.pop_heap();
1379
    }else{
1380
      // does not belong to the block
1381
      Assert((heap.direction() > 0) ? (blockValue < top.d_diff)
1382
                                    : (blockValue > top.d_diff));
1383
      break;
1384
    }
1385
  }
1386
}
1387
1388
void LinearEqualityModule::substitutePlusTimesConstant(ArithVar to, ArithVar from, const Rational& mult){
1389
  d_tableau.substitutePlusTimesConstant(to, from, mult, d_trackCallback);
1390
}
1391
void LinearEqualityModule::directlyAddToCoefficient(ArithVar row, ArithVar col, const Rational& mult){
1392
  d_tableau.directlyAddToCoefficient(row, col, mult, d_trackCallback);
1393
}
1394
1395
}  // namespace arith
1396
}  // namespace theory
1397
29577
}  // namespace cvc5