GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/linear_equality.cpp Lines: 368 806 45.7 %
Date: 2021-08-06 Branches: 624 3316 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
9853
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
9853
  d_trackCallback(this)
66
9853
{}
67
68
9853
LinearEqualityModule::Statistics::Statistics()
69
    : d_statPivots(
70
19706
        smtStatisticsRegistry().registerInt("theory::arith::pivots")),
71
      d_statUpdates(
72
19706
          smtStatisticsRegistry().registerInt("theory::arith::updates")),
73
      d_pivotTime(
74
19706
          smtStatisticsRegistry().registerTimer("theory::arith::pivotTime")),
75
      d_adjTime(
76
19706
          smtStatisticsRegistry().registerTimer("theory::arith::adjTime")),
77
9853
      d_weakeningAttempts(smtStatisticsRegistry().registerInt(
78
19706
          "theory::arith::weakening::attempts")),
79
9853
      d_weakeningSuccesses(smtStatisticsRegistry().registerInt(
80
19706
          "theory::arith::weakening::success")),
81
9853
      d_weakenings(smtStatisticsRegistry().registerInt(
82
19706
          "theory::arith::weakening::total")),
83
9853
      d_weakenTime(smtStatisticsRegistry().registerTimer(
84
19706
          "theory::arith::weakening::time")),
85
      d_forceTime(
86
88677
          smtStatisticsRegistry().registerTimer("theory::arith::forcing::time"))
87
{
88
9853
}
89
90
4618423
void LinearEqualityModule::includeBoundUpdate(ArithVar v, const BoundsInfo& prev){
91
4618423
  Assert(!d_areTracking);
92
93
4618423
  BoundsInfo curr = d_variables.boundsInfo(v);
94
95
4618423
  Assert(prev != curr);
96
4618423
  Tableau::ColIterator basicIter = d_tableau.colIterator(v);
97
49600369
  for(; !basicIter.atEnd(); ++basicIter){
98
22490973
    const Tableau::Entry& entry = *basicIter;
99
22490973
    Assert(entry.getColVar() == v);
100
22490973
    int a_ijSgn = entry.getCoefficient().sgn();
101
102
22490973
    RowIndex ridx = entry.getRowIndex();
103
22490973
    BoundsInfo& counts = d_btracking.get(ridx);
104
22490973
    Debug("includeBoundUpdate") << d_tableau.rowIndexToBasic(ridx) << " " << counts << " to " ;
105
22490973
    counts.addInChange(a_ijSgn, prev, curr);
106
22490973
    Debug("includeBoundUpdate") << counts << " " << a_ijSgn << std::endl;
107
  }
108
4618423
}
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
166583
void LinearEqualityModule::updateUntracked(ArithVar x_i, const DeltaRational& v){
181
166583
  Assert(!d_tableau.isBasic(x_i));
182
166583
  Assert(!d_areTracking);
183
166583
  const DeltaRational& assignment_x_i = d_variables.getAssignment(x_i);
184
166583
  ++(d_statistics.d_statUpdates);
185
186
187
333166
  Debug("arith") <<"update " << x_i << ": "
188
166583
                 << assignment_x_i << "|-> " << v << endl;
189
333166
  DeltaRational diff = v - assignment_x_i;
190
191
166583
  Tableau::ColIterator colIter = d_tableau.colIterator(x_i);
192
3667843
  for(; !colIter.atEnd(); ++colIter){
193
1750630
    const Tableau::Entry& entry = *colIter;
194
1750630
    Assert(entry.getColVar() == x_i);
195
196
1750630
    ArithVar x_j = d_tableau.rowIndexToBasic(entry.getRowIndex());
197
1750630
    const Rational& a_ji = entry.getCoefficient();
198
199
1750630
    const DeltaRational& assignment = d_variables.getAssignment(x_j);
200
3501260
    DeltaRational  nAssignment = assignment+(diff * a_ji);
201
1750630
    d_variables.setAssignment(x_j, nAssignment);
202
203
1750630
    d_basicVariableUpdates(x_j);
204
  }
205
206
166583
  d_variables.setAssignment(x_i, v);
207
208
166583
  if(Debug.isOn("paranoid:check_tableau")){  debugCheckTableau(); }
209
166583
}
210
211
198030
void LinearEqualityModule::updateTracked(ArithVar x_i, const DeltaRational& v){
212
396060
  TimerStat::CodeTimer codeTimer(d_statistics.d_adjTime);
213
214
198030
  Assert(!d_tableau.isBasic(x_i));
215
198030
  Assert(d_areTracking);
216
217
198030
  ++(d_statistics.d_statUpdates);
218
219
396060
  DeltaRational diff =  v - d_variables.getAssignment(x_i);
220
396060
  Debug("arith") <<"update " << x_i << ": "
221
198030
                 << d_variables.getAssignment(x_i) << "|-> " << v << endl;
222
223
224
198030
  BoundCounts before = d_variables.atBoundCounts(x_i);
225
198030
  d_variables.setAssignment(x_i, v);
226
198030
  BoundCounts after = d_variables.atBoundCounts(x_i);
227
228
198030
  bool anyChange = before != after;
229
230
198030
  Tableau::ColIterator colIter = d_tableau.colIterator(x_i);
231
7615620
  for(; !colIter.atEnd(); ++colIter){
232
3708795
    const Tableau::Entry& entry = *colIter;
233
3708795
    Assert(entry.getColVar() == x_i);
234
235
3708795
    RowIndex ridx = entry.getRowIndex();
236
3708795
    ArithVar x_j = d_tableau.rowIndexToBasic(ridx);
237
3708795
    const Rational& a_ji = entry.getCoefficient();
238
239
3708795
    const DeltaRational& assignment = d_variables.getAssignment(x_j);
240
7417590
    DeltaRational  nAssignment = assignment+(diff * a_ji);
241
3708795
    Debug("update") << x_j << " " << a_ji << assignment << " -> " << nAssignment << endl;
242
3708795
    BoundCounts xjBefore = d_variables.atBoundCounts(x_j);
243
3708795
    d_variables.setAssignment(x_j, nAssignment);
244
3708795
    BoundCounts xjAfter = d_variables.atBoundCounts(x_j);
245
246
3708795
    Assert(rowIndexIsTracked(ridx));
247
3708795
    BoundsInfo& next_bc_k = d_btracking.get(ridx);
248
3708795
    if(anyChange){
249
2682916
      next_bc_k.addInAtBoundChange(a_ji.sgn(), before, after);
250
    }
251
3708795
    if(xjBefore != xjAfter){
252
449291
      next_bc_k.addInAtBoundChange(-1, xjBefore, xjAfter);
253
    }
254
255
3708795
    d_basicVariableUpdates(x_j);
256
  }
257
258
198030
  if(Debug.isOn("paranoid:check_tableau")){  debugCheckTableau(); }
259
198030
}
260
261
198030
void LinearEqualityModule::pivotAndUpdate(ArithVar x_i, ArithVar x_j, const DeltaRational& x_i_value){
262
198030
  Assert(x_i != x_j);
263
264
396060
  TimerStat::CodeTimer codeTimer(d_statistics.d_pivotTime);
265
266
  static int instance = 0;
267
268
198030
  if(Debug.isOn("arith::tracking::pre")){
269
    ++instance;
270
    Debug("arith::tracking")  << "pre update #" << instance << endl;
271
    debugCheckTracking();
272
  }
273
274
198030
  if(Debug.isOn("arith::simplex:row")){ debugPivot(x_i, x_j); }
275
276
198030
  RowIndex ridx = d_tableau.basicToRowIndex(x_i);
277
198030
  const Tableau::Entry& entry_ij =  d_tableau.findEntry(ridx, x_j);
278
198030
  Assert(!entry_ij.blank());
279
280
198030
  const Rational& a_ij = entry_ij.getCoefficient();
281
198030
  const DeltaRational& betaX_i = d_variables.getAssignment(x_i);
282
396060
  DeltaRational theta = (x_i_value - betaX_i)/a_ij;
283
396060
  DeltaRational x_j_value = d_variables.getAssignment(x_j) + theta;
284
285
198030
  updateTracked(x_j, x_j_value);
286
287
198030
  if(Debug.isOn("arith::tracking::mid")){
288
    Debug("arith::tracking")  << "postupdate prepivot #" << instance << endl;
289
    debugCheckTracking();
290
  }
291
292
  // Pivots
293
198030
  ++(d_statistics.d_statPivots);
294
295
198030
  d_tableau.pivot(x_i, x_j, d_trackCallback);
296
297
198030
  if(Debug.isOn("arith::tracking::post")){
298
    Debug("arith::tracking")  << "postpivot #" << instance << endl;
299
    debugCheckTracking();
300
  }
301
302
198030
  d_basicVariableUpdates(x_j);
303
304
198030
  if(Debug.isOn("matrix")){
305
    d_tableau.printMatrix();
306
  }
307
198030
}
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
447445
DeltaRational LinearEqualityModule::computeRowBound(RowIndex ridx, bool rowUb, ArithVar skip) const {
427
447445
  DeltaRational sum(0,0);
428
2814910
  for(Tableau::RowIterator i = d_tableau.ridRowIterator(ridx); !i.atEnd(); ++i){
429
2367465
    const Tableau::Entry& entry = (*i);
430
2367465
    ArithVar v = entry.getColVar();
431
2367465
    if(v == skip){ continue; }
432
433
2063101
    const Rational& coeff =  entry.getCoefficient();
434
2063101
    bool vUb = (rowUb == (coeff.sgn() > 0));
435
436
4126202
    const DeltaRational& bound = vUb ?
437
1046830
      d_variables.getUpperBound(v):
438
3079372
      d_variables.getLowerBound(v);
439
440
4126202
    DeltaRational diff = bound * coeff;
441
2063101
    sum = sum + diff;
442
  }
443
447445
  return sum;
444
}
445
446
/**
447
 * Computes the value of a basic variable using the current assignment.
448
 */
449
184996
DeltaRational LinearEqualityModule::computeRowValue(ArithVar x, bool useSafe) const{
450
184996
  Assert(d_tableau.isBasic(x));
451
184996
  DeltaRational sum(0);
452
453
1056906
  for(Tableau::RowIterator i = d_tableau.basicRowIterator(x); !i.atEnd(); ++i){
454
871910
    const Tableau::Entry& entry = (*i);
455
871910
    ArithVar nonbasic = entry.getColVar();
456
871910
    if(nonbasic == x) continue;
457
686914
    const Rational& coeff = entry.getCoefficient();
458
459
686914
    const DeltaRational& assignment = d_variables.getAssignment(nonbasic, useSafe);
460
686914
    sum = sum + (assignment * coeff);
461
  }
462
184996
  return sum;
463
}
464
465
4106206
const Tableau::Entry* LinearEqualityModule::rowLacksBound(RowIndex ridx, bool rowUb, ArithVar skip){
466
4106206
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
467
36394782
  for(; !iter.atEnd(); ++iter){
468
20248748
    const Tableau::Entry& entry = *iter;
469
470
20248748
    ArithVar var = entry.getColVar();
471
20248748
    if(var == skip) { continue; }
472
473
20246056
    int sgn = entry.getCoefficient().sgn();
474
20246056
    bool selectUb = (rowUb == (sgn > 0));
475
40492112
    bool hasBound = selectUb ?
476
10184812
      d_variables.hasUpperBound(var):
477
30307300
      d_variables.hasLowerBound(var);
478
20246056
    if(!hasBound){
479
4104460
      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
130422
void LinearEqualityModule::propagateRow(ConstraintCPVec& into, RowIndex ridx, bool rowUp, ConstraintP c, RationalVectorP farkas){
533
130422
  Assert(!c->assertedToTheTheory());
534
130422
  Assert(c->canBePropagated());
535
130422
  Assert(!c->hasProof());
536
537
130422
  if(farkas != RationalVectorPSentinel){
538
71653
    Assert(farkas->empty());
539
71653
    farkas->push_back(Rational(0));
540
  }
541
542
130422
  ArithVar v = c->getVariable();
543
260844
  Debug("arith::propagateRow") << "LinearEqualityModule::propagateRow("
544
130422
                                   << ridx << ", " << rowUp << ", " << v << ") start" << endl;
545
546
130422
  const Rational& multiple = rowUp ? d_one : d_negOne;
547
548
130422
  Debug("arith::propagateRow") << "multiple: " << multiple << endl;
549
550
130422
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
551
1999194
  for(; !iter.atEnd(); ++iter){
552
934386
    const Tableau::Entry& entry = *iter;
553
934386
    ArithVar nonbasic = entry.getColVar();
554
934386
    const Rational& a_ij = entry.getCoefficient();
555
934386
    int sgn = a_ij.sgn();
556
934386
    Assert(sgn != 0);
557
934386
    bool selectUb = rowUp ? (sgn > 0) : (sgn < 0);
558
559
934386
    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
934386
    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
934386
    if(nonbasic == v){
578
130422
      if(farkas != RationalVectorPSentinel){
579
71653
        Assert(farkas->front().isZero());
580
143306
        Rational multAij = multiple * a_ij;
581
71653
        Debug("arith::propagateRow") << "(" << multAij << ") ";
582
71653
        farkas->front() = multAij;
583
      }
584
585
130422
      Debug("arith::propagateRow") << c << endl;
586
    }else{
587
588
      ConstraintCP bound = selectUb
589
361210
        ? d_variables.getUpperBoundConstraint(nonbasic)
590
1165174
        : d_variables.getLowerBoundConstraint(nonbasic);
591
592
803964
      if(farkas != RationalVectorPSentinel){
593
828210
        Rational multAij = multiple * a_ij;
594
414105
        Debug("arith::propagateRow") << "(" << multAij << ") ";
595
414105
        farkas->push_back(multAij);
596
      }
597
803964
      Assert(bound != NullConstraint);
598
803964
      Debug("arith::propagateRow") << bound << endl;
599
803964
      into.push_back(bound);
600
    }
601
  }
602
260844
  Debug("arith::propagateRow") << "LinearEqualityModule::propagateRow("
603
130422
                                   << ridx << ", " << rowUp << ", " << v << ") done" << endl;
604
605
130422
}
606
607
883495
ConstraintP LinearEqualityModule::weakestExplanation(bool aboveUpper, DeltaRational& surplus, ArithVar v, const Rational& coeff, bool& anyWeakening, ArithVar basic) const {
608
609
883495
  int sgn = coeff.sgn();
610
883495
  bool ub = aboveUpper?(sgn < 0) : (sgn > 0);
611
612
1766990
  ConstraintP c = ub ?
613
250970
    d_variables.getUpperBoundConstraint(v) :
614
1516020
    d_variables.getLowerBoundConstraint(v);
615
616
  bool weakened;
617
970426
  do{
618
970426
    const DeltaRational& bound = c->getValue();
619
620
970426
    weakened = false;
621
622
1940852
    ConstraintP weaker = ub?
623
306533
      c->getStrictlyWeakerUpperBound(true, true):
624
1634319
      c->getStrictlyWeakerLowerBound(true, true);
625
626
970426
    if(weaker != NullConstraint){
627
216879
      const DeltaRational& weakerBound = weaker->getValue();
628
629
433758
      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
216879
      diff = diff * coeff;
634
216879
      if(surplus > diff){
635
86931
        ++d_statistics.d_weakenings;
636
86931
        weakened = true;
637
86931
        anyWeakening = true;
638
86931
        surplus = surplus - diff;
639
640
86931
        Debug("arith::weak") << "found:" << endl;
641
86931
        if(v == basic){
642
8549
          Debug("arith::weak") << "  basic: ";
643
        }
644
173862
        Debug("arith::weak") << "  " << surplus << " "<< diff  << endl
645
86931
                             << "  " << bound << c << endl
646
86931
                             << "  " << weakerBound << weaker << endl;
647
648
86931
        Assert(diff.sgn() > 0);
649
86931
        c = weaker;
650
      }
651
    }
652
  }while(weakened);
653
654
883495
  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
61696
ConstraintCP LinearEqualityModule::minimallyWeakConflict(bool aboveUpper, ArithVar basicVar, FarkasConflictBuilder& fcs) const {
692
61696
  Assert(!fcs.underConstruction());
693
123392
  TimerStat::CodeTimer codeTimer(d_statistics.d_weakenTime);
694
695
123392
  Debug("arith::weak") << "LinearEqualityModule::minimallyWeakConflict("
696
61696
                       << aboveUpper <<", "<< basicVar << ", ...) start" << endl;
697
698
61696
  const Rational& adjustSgn = aboveUpper ? d_negOne : d_one;
699
61696
  const DeltaRational& assignment = d_variables.getAssignment(basicVar);
700
123392
  DeltaRational surplus;
701
61696
  if(aboveUpper){
702
27889
    Assert(d_variables.hasUpperBound(basicVar));
703
27889
    Assert(assignment > d_variables.getUpperBound(basicVar));
704
27889
    surplus = assignment - d_variables.getUpperBound(basicVar);
705
  }else{
706
33807
    Assert(d_variables.hasLowerBound(basicVar));
707
33807
    Assert(assignment < d_variables.getLowerBound(basicVar));
708
33807
    surplus = d_variables.getLowerBound(basicVar) - assignment;
709
  }
710
711
61696
  bool anyWeakenings = false;
712
945191
  for(Tableau::RowIterator i = d_tableau.basicRowIterator(basicVar); !i.atEnd(); ++i){
713
883495
    const Tableau::Entry& entry = *i;
714
883495
    ArithVar v = entry.getColVar();
715
883495
    const Rational& coeff = entry.getCoefficient();
716
883495
    bool weakening = false;
717
883495
    ConstraintP c = weakestExplanation(aboveUpper, surplus, v, coeff, weakening, basicVar);
718
1766990
    Debug("arith::weak") << "weak : " << weakening << " "
719
1766990
                         << c->assertedToTheTheory() << " "
720
883495
                         << d_variables.getAssignment(v) << " "
721
883495
                         << c << endl;
722
883495
    anyWeakenings = anyWeakenings || weakening;
723
724
883495
    fcs.addConstraint(c, coeff, adjustSgn);
725
883495
    if(basicVar == v){
726
61696
      Assert(!c->negationHasProof());
727
61696
      fcs.makeLastConsequent();
728
    }
729
  }
730
61696
  Assert(fcs.consequentIsSet());
731
732
61696
  ConstraintCP conflicted = fcs.commitConflict();
733
734
61696
  ++d_statistics.d_weakeningAttempts;
735
61696
  if(anyWeakenings){
736
44556
    ++d_statistics.d_weakeningSuccesses;
737
  }
738
123392
  Debug("arith::weak") << "LinearEqualityModule::minimallyWeakConflict("
739
61696
                       << aboveUpper <<", "<< basicVar << ", ...) done" << endl;
740
123392
  return conflicted;
741
}
742
743
130405
ArithVar LinearEqualityModule::minVarOrder(ArithVar x, ArithVar y) const {
744
130405
  Assert(x != ARITHVAR_SENTINEL);
745
130405
  Assert(y != ARITHVAR_SENTINEL);
746
130405
  if(x <= y){
747
79539
    return x;
748
  } else {
749
50866
    return y;
750
  }
751
}
752
753
808219
ArithVar LinearEqualityModule::minColLength(ArithVar x, ArithVar y) const {
754
808219
  Assert(x != ARITHVAR_SENTINEL);
755
808219
  Assert(y != ARITHVAR_SENTINEL);
756
808219
  Assert(!d_tableau.isBasic(x));
757
808219
  Assert(!d_tableau.isBasic(y));
758
808219
  uint32_t xLen = d_tableau.getColLength(x);
759
808219
  uint32_t yLen = d_tableau.getColLength(y);
760
808219
  if( xLen > yLen){
761
114370
     return y;
762
693849
  } else if( xLen== yLen ){
763
129892
    return minVarOrder(x,y);
764
  }else{
765
563957
    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
848667
ArithVar LinearEqualityModule::minBoundAndColLength(ArithVar x, ArithVar y) const{
786
848667
  Assert(x != ARITHVAR_SENTINEL);
787
848667
  Assert(y != ARITHVAR_SENTINEL);
788
848667
  Assert(!d_tableau.isBasic(x));
789
848667
  Assert(!d_tableau.isBasic(y));
790
848667
  if(d_variables.hasEitherBound(x) && !d_variables.hasEitherBound(y)){
791
14428
    return y;
792
834239
  }else if(!d_variables.hasEitherBound(x) && d_variables.hasEitherBound(y)){
793
26020
    return x;
794
  }else {
795
808219
    return minColLength(x, y);
796
  }
797
}
798
799
template <bool above>
800
198030
ArithVar LinearEqualityModule::selectSlack(ArithVar x_i, VarPreferenceFunction pref) const{
801
198030
  ArithVar slack = ARITHVAR_SENTINEL;
802
803
5060812
  for(Tableau::RowIterator iter = d_tableau.basicRowIterator(x_i); !iter.atEnd();  ++iter){
804
4862782
    const Tableau::Entry& entry = *iter;
805
4862782
    ArithVar nonbasic = entry.getColVar();
806
4862782
    if(nonbasic == x_i) continue;
807
808
4664752
    const Rational& a_ij = entry.getCoefficient();
809
4664752
    int sgn = a_ij.sgn();
810
4664752
    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
1047210
      slack = (slack == ARITHVAR_SENTINEL) ? nonbasic : (this->*pref)(slack, nonbasic);
815
    }
816
  }
817
818
198030
  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
1635058
void LinearEqualityModule::startTrackingBoundCounts(){
843
1635058
  Assert(!d_areTracking);
844
1635058
  d_areTracking = true;
845
1635058
  if(Debug.isOn("arith::tracking")){
846
    debugCheckTracking();
847
  }
848
1635058
  Assert(d_areTracking);
849
1635058
}
850
851
1635058
void LinearEqualityModule::stopTrackingBoundCounts(){
852
1635058
  Assert(d_areTracking);
853
1635058
  d_areTracking = false;
854
1635058
  if(Debug.isOn("arith::tracking")){
855
    debugCheckTracking();
856
  }
857
1635058
  Assert(!d_areTracking);
858
1635058
}
859
860
861
92498
void LinearEqualityModule::trackRowIndex(RowIndex ridx){
862
92498
  Assert(!rowIndexIsTracked(ridx));
863
92498
  BoundsInfo bi = computeRowBoundInfo(ridx, true);
864
92498
  d_btracking.set(ridx, bi);
865
92498
}
866
867
92498
BoundsInfo LinearEqualityModule::computeRowBoundInfo(RowIndex ridx, bool inQueue) const{
868
92498
  BoundsInfo bi;
869
870
92498
  Tableau::RowIterator iter = d_tableau.ridRowIterator(ridx);
871
964408
  for(; !iter.atEnd();  ++iter){
872
435955
    const Tableau::Entry& entry = *iter;
873
435955
    ArithVar v = entry.getColVar();
874
435955
    const Rational& a_ij = entry.getCoefficient();
875
435955
    bi += (d_variables.selectBoundsInfo(v, inQueue)).multiplyBySgn(a_ij.sgn());
876
  }
877
92498
  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
411734
bool LinearEqualityModule::nonbasicsAtLowerBounds(ArithVar basic) const {
930
411734
  Assert(basicIsTracked(basic));
931
411734
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
932
933
411734
  BoundCounts bcs = d_btracking[ridx].atBounds();
934
411734
  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
411734
  uint32_t lbc = bcs.lowerBoundCount();
944
962913
  return  (lbc == length) ||
945
962913
    (lbc + 1 == length && d_variables.cmpAssignmentUpperBound(basic) != 0);
946
}
947
948
454261
bool LinearEqualityModule::nonbasicsAtUpperBounds(ArithVar basic) const {
949
454261
  Assert(basicIsTracked(basic));
950
454261
  RowIndex ridx = d_tableau.basicToRowIndex(basic);
951
454261
  BoundCounts bcs = d_btracking[ridx].atBounds();
952
454261
  uint32_t length = d_tableau.basicRowLength(basic);
953
454261
  uint32_t ubc = bcs.upperBoundCount();
954
  // See the comment for nonbasicsAtLowerBounds()
955
956
1077557
  return (ubc == length) ||
957
1077557
    (ubc + 1 == length && d_variables.cmpAssignmentLowerBound(basic) != 0);
958
}
959
960
198030
void LinearEqualityModule::trackingMultiplyRow(RowIndex ridx, int sgn) {
961
198030
  Assert(rowIndexIsTracked(ridx));
962
198030
  Assert(sgn != 0);
963
198030
  if(sgn < 0){
964
110979
    BoundsInfo& bi = d_btracking.get(ridx);
965
110979
    bi = bi.multiplyBySgn(sgn);
966
  }
967
198030
}
968
969
52865117
void LinearEqualityModule::trackingCoefficientChange(RowIndex ridx, ArithVar nb, int oldSgn, int currSgn){
970
52865117
  Assert(oldSgn != currSgn);
971
52865117
  BoundsInfo nb_inf = d_variables.boundsInfo(nb);
972
973
52865117
  Assert(rowIndexIsTracked(ridx));
974
975
52865117
  BoundsInfo& row_bi = d_btracking.get(ridx);
976
52865117
  row_bi.addInSgn(nb_inf, oldSgn, currSgn);
977
52865117
}
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.nothing()
1058
        || nbDiff <= d_upperBoundDifference.value())
1059
    {
1060
      return false;
1061
    }
1062
  }else{
1063
    if (d_lowerBoundDifference.nothing()
1064
        || nbDiff >= d_lowerBoundDifference.value())
1065
    {
1066
      return false;
1067
    }
1068
  }
1069
1070
  // Assume past this point, nb will be in error if this pivot is done
1071
  ArithVar nb = entry.getColVar();
1072
  RowIndex ridx = entry.getRowIndex();
1073
  ArithVar basic = d_tableau.rowIndexToBasic(ridx);
1074
  Assert(rowIndexIsTracked(ridx));
1075
  int coeffSgn = entry.getCoefficient().sgn();
1076
1077
1078
  // if bToUB, then basic is going to be set to its upperbound
1079
  // if not bToUB, then basic is going to be set to its lowerbound
1080
1081
  // Different steps of solving for this:
1082
  // 1) y = a * x + \sum b * z
1083
  // 2) -a * x = -y + \sum b * z
1084
  // 3) x = (-1/a) * ( -y + \sum b * z)
1085
1086
  BoundCounts bc = d_btracking[ridx].atBounds();
1087
1088
  // 1) y = a * x + \sum b * z
1089
  // Get bc(\sum b * z)
1090
  BoundCounts sumOnly = bc - d_variables.atBoundCounts(nb).multiplyBySgn(coeffSgn);
1091
1092
  // y's bounds in the proposed model
1093
  int yWillBeAtUb = (bToUB || d_variables.boundsAreEqual(basic)) ? 1 : 0;
1094
  int yWillBeAtLb = (!bToUB || d_variables.boundsAreEqual(basic)) ? 1 : 0;
1095
  BoundCounts ysBounds(yWillBeAtLb, yWillBeAtUb);
1096
1097
  // 2) -a * x = -y + \sum b * z
1098
  // Get bc(-y + \sum b * z)
1099
  sumOnly.addInChange(-1, d_variables.atBoundCounts(basic), ysBounds);
1100
1101
  // 3) x = (-1/a) * ( -y + \sum b * z)
1102
  // Get bc((-1/a) * ( -y + \sum b * z))
1103
  BoundCounts xsBoundsAfterPivot = sumOnly.multiplyBySgn(-coeffSgn);
1104
1105
  uint32_t length = d_tableau.basicRowLength(basic);
1106
  if(nbSgn > 0){
1107
    // Only check for the upper bound being violated
1108
    return xsBoundsAfterPivot.lowerBoundCount() + 1 == length;
1109
  }else{
1110
    // Only check for the lower bound being violated
1111
    return xsBoundsAfterPivot.upperBoundCount() + 1 == length;
1112
  }
1113
}
1114
1115
UpdateInfo LinearEqualityModule::mkConflictUpdate(const Tableau::Entry& entry, bool ub) const{
1116
  ArithVar currBasic = d_tableau.rowIndexToBasic(entry.getRowIndex());
1117
  ArithVar nb = entry.getColVar();
1118
1119
  ConstraintP bound = ub ?
1120
    d_variables.getUpperBoundConstraint(currBasic):
1121
    d_variables.getLowerBoundConstraint(currBasic);
1122
1123
1124
  const Rational& coeff = entry.getCoefficient();
1125
  const DeltaRational& assignment = d_variables.getAssignment(currBasic);
1126
  DeltaRational toBound = bound->getValue() - assignment;
1127
  DeltaRational nbDiff = toBound/coeff;
1128
1129
  return UpdateInfo::conflict(nb, nbDiff, coeff, bound);
1130
}
1131
1132
UpdateInfo LinearEqualityModule::speculativeUpdate(ArithVar nb, const Rational& focusCoeff, UpdatePreferenceFunction pref){
1133
  Assert(d_increasing.empty());
1134
  Assert(d_decreasing.empty());
1135
  Assert(d_lowerBoundDifference.nothing());
1136
  Assert(d_upperBoundDifference.nothing());
1137
1138
  int focusCoeffSgn = focusCoeff.sgn();
1139
1140
  static int instance = 0;
1141
  ++instance;
1142
  Debug("speculativeUpdate") << "speculativeUpdate " << instance << endl;
1143
  Debug("speculativeUpdate") << "nb " << nb << endl;
1144
  Debug("speculativeUpdate") << "focusCoeff " << focusCoeff << endl;
1145
1146
  if(d_variables.hasUpperBound(nb)){
1147
    ConstraintP ub = d_variables.getUpperBoundConstraint(nb);
1148
    d_upperBoundDifference = ub->getValue() - d_variables.getAssignment(nb);
1149
    Border border(ub, d_upperBoundDifference.value(), false, NULL, true);
1150
    Debug("handleBorders") << "push back increasing " << border << endl;
1151
    d_increasing.push_back(border);
1152
  }
1153
  if(d_variables.hasLowerBound(nb)){
1154
    ConstraintP lb = d_variables.getLowerBoundConstraint(nb);
1155
    d_lowerBoundDifference = lb->getValue() - d_variables.getAssignment(nb);
1156
    Border border(lb, d_lowerBoundDifference.value(), false, NULL, false);
1157
    Debug("handleBorders") << "push back decreasing " << border << endl;
1158
    d_decreasing.push_back(border);
1159
  }
1160
1161
  Tableau::ColIterator colIter = d_tableau.colIterator(nb);
1162
  for(; !colIter.atEnd(); ++colIter){
1163
    const Tableau::Entry& entry = *colIter;
1164
    Assert(entry.getColVar() == nb);
1165
1166
    if(accumulateBorder(entry, true)){
1167
      clearSpeculative();
1168
      return mkConflictUpdate(entry, true);
1169
    }
1170
    if(accumulateBorder(entry, false)){
1171
      clearSpeculative();
1172
      return mkConflictUpdate(entry, false);
1173
    }
1174
  }
1175
1176
  UpdateInfo selected;
1177
  BorderHeap& withSgn = focusCoeffSgn > 0 ? d_increasing : d_decreasing;
1178
  BorderHeap& againstSgn = focusCoeffSgn > 0 ? d_decreasing : d_increasing;
1179
1180
  handleBorders(selected, nb, focusCoeff, withSgn, 0, pref);
1181
  int m = 1 - selected.errorsChangeSafe(0);
1182
  handleBorders(selected, nb, focusCoeff, againstSgn, m, pref);
1183
1184
  clearSpeculative();
1185
  return selected;
1186
}
1187
1188
void LinearEqualityModule::clearSpeculative(){
1189
  // clear everything away
1190
  d_increasing.clear();
1191
  d_decreasing.clear();
1192
  d_lowerBoundDifference.clear();
1193
  d_upperBoundDifference.clear();
1194
}
1195
1196
void LinearEqualityModule::handleBorders(UpdateInfo& selected, ArithVar nb, const Rational& focusCoeff, BorderHeap& heap, int minimumFixes, UpdatePreferenceFunction pref){
1197
  Assert(minimumFixes >= 0);
1198
1199
  // The values popped off of the heap
1200
  // should be popped with the values closest to 0
1201
  // being first and larger in absolute value last
1202
1203
1204
  int fixesRemaining = heap.possibleFixes();
1205
1206
  Debug("handleBorders")
1207
    << "handleBorders "
1208
    << "nb " << nb
1209
    << "fc " << focusCoeff
1210
    << "h.e " << heap.empty()
1211
    << "h.dir " << heap.direction()
1212
    << "h.rem " << fixesRemaining
1213
    << "h.0s " << heap.numZeroes()
1214
    << "min " << minimumFixes
1215
    << endl;
1216
1217
  if(heap.empty()){
1218
    // if the heap is empty, return
1219
    return;
1220
  }
1221
1222
  bool zeroesWillDominate = fixesRemaining - heap.numZeroes() < minimumFixes;
1223
1224
  // can the number of fixes ever exceed the minimum?
1225
  // no more than the number of possible fixes can be fixed in total
1226
  // nothing can be fixed before the zeroes are taken care of
1227
  if(minimumFixes > 0 && zeroesWillDominate){
1228
    return;
1229
  }
1230
1231
1232
  int negErrorChange = 0;
1233
  int nbDir = heap.direction();
1234
1235
  // points at the beginning of the heap
1236
  if(zeroesWillDominate){
1237
    heap.dropNonZeroes();
1238
  }
1239
  heap.make_heap();
1240
1241
1242
  // pretend like the previous block had a value of zero.
1243
  // The block that actually has a value of 0 must handle this.
1244
  const DeltaRational zero(0);
1245
  const DeltaRational* prevBlockValue = &zero;
1246
1247
  /** The coefficient changes as the value crosses border. */
1248
  Rational effectiveCoefficient = focusCoeff;
1249
1250
  /* Keeps track of the change to the value of the focus function.*/
1251
  DeltaRational totalFocusChange(0);
1252
1253
  const int focusCoeffSgn = focusCoeff.sgn();
1254
1255
  while(heap.more()  &&
1256
        (fixesRemaining + negErrorChange > minimumFixes ||
1257
         (fixesRemaining + negErrorChange == minimumFixes &&
1258
          effectiveCoefficient.sgn() == focusCoeffSgn))){
1259
    // There are more elements &&
1260
    // we can either fix at least 1 more variable in the error function
1261
    // or we can improve the error function
1262
1263
1264
    int brokenInBlock = 0;
1265
    BorderVec::const_iterator endBlock = heap.end();
1266
1267
    pop_block(heap, brokenInBlock, fixesRemaining, negErrorChange);
1268
1269
    // if endVec == beginVec, block starts there
1270
    // other wise, block starts at endVec
1271
    BorderVec::const_iterator startBlock
1272
      = heap.more() ? heap.end() : heap.begin();
1273
1274
    const DeltaRational& blockValue = (*startBlock).d_diff;
1275
1276
    // if decreasing
1277
    // blockValue < prevBlockValue
1278
    // diff.sgn() = -1
1279
    DeltaRational diff = blockValue - (*prevBlockValue);
1280
    DeltaRational blockChangeToFocus =  diff * effectiveCoefficient;
1281
    totalFocusChange += blockChangeToFocus;
1282
1283
    Debug("handleBorders")
1284
      << "blockValue " << (blockValue)
1285
      << "diff " << diff
1286
      << "blockChangeToFocus " << totalFocusChange
1287
      << "blockChangeToFocus " << totalFocusChange
1288
      << "negErrorChange " << negErrorChange
1289
      << "brokenInBlock " << brokenInBlock
1290
      << "fixesRemaining " << fixesRemaining
1291
      << endl;
1292
1293
    int currFocusChangeSgn = totalFocusChange.sgn();
1294
    for(BorderVec::const_iterator i = startBlock; i != endBlock; ++i){
1295
      const Border& b = *i;
1296
1297
      Debug("handleBorders") << b << endl;
1298
1299
      bool makesImprovement = negErrorChange > 0 ||
1300
        (negErrorChange == 0  && currFocusChangeSgn > 0);
1301
1302
      if(!makesImprovement){
1303
        if(b.ownBorder() || minimumFixes > 0){
1304
          continue;
1305
        }
1306
      }
1307
1308
      UpdateInfo proposal(nb, nbDir);
1309
      if(b.ownBorder()){
1310
        proposal.witnessedUpdate(b.d_diff, b.d_bound, -negErrorChange, currFocusChangeSgn);
1311
      }else{
1312
        proposal.update(b.d_diff, b.getCoefficient(), b.d_bound, -negErrorChange, currFocusChangeSgn);
1313
      }
1314
1315
      if(selected.unbounded() || (this->*pref)(selected, proposal)){
1316
        selected = proposal;
1317
      }
1318
    }
1319
1320
    effectiveCoefficient += updateCoefficient(startBlock, endBlock);
1321
    prevBlockValue = &blockValue;
1322
    negErrorChange -= brokenInBlock;
1323
  }
1324
}
1325
1326
Rational LinearEqualityModule::updateCoefficient(BorderVec::const_iterator startBlock, BorderVec::const_iterator endBlock){
1327
  //update coefficient
1328
  Rational changeToCoefficient(0);
1329
  for(BorderVec::const_iterator i = startBlock; i != endBlock; ++i){
1330
    const Border& curr = *i;
1331
    if(curr.ownBorder()){// breaking its own bound
1332
      if(curr.d_upperbound){
1333
        changeToCoefficient -= 1;
1334
      }else{
1335
        changeToCoefficient += 1;
1336
      }
1337
    }else{
1338
      const Rational& coeff = curr.d_entry->getCoefficient();
1339
      if(curr.d_areFixing){
1340
        if(curr.d_upperbound){// fixing an upper bound
1341
          changeToCoefficient += coeff;
1342
        }else{// fixing a lower bound
1343
          changeToCoefficient -= coeff;
1344
        }
1345
      }else{
1346
        if(curr.d_upperbound){// breaking an upper bound
1347
          changeToCoefficient -= coeff;
1348
        }else{
1349
          // breaking a lower bound
1350
          changeToCoefficient += coeff;
1351
        }
1352
      }
1353
    }
1354
  }
1355
  return changeToCoefficient;
1356
}
1357
1358
void LinearEqualityModule::pop_block(BorderHeap& heap, int& brokenInBlock, int& fixesRemaining, int& negErrorChange){
1359
  Assert(heap.more());
1360
1361
  if(heap.top().d_areFixing){
1362
    fixesRemaining--;
1363
    negErrorChange++;
1364
  }else{
1365
    brokenInBlock++;
1366
  }
1367
  heap.pop_heap();
1368
  const DeltaRational& blockValue = (*heap.end()).d_diff;
1369
1370
  while(heap.more()){
1371
    const Border& top = heap.top();
1372
    if(blockValue == top.d_diff){
1373
      // belongs to the block
1374
      if(top.d_areFixing){
1375
        fixesRemaining--;
1376
        negErrorChange++;
1377
      }else{
1378
        brokenInBlock++;
1379
      }
1380
      heap.pop_heap();
1381
    }else{
1382
      // does not belong to the block
1383
      Assert((heap.direction() > 0) ? (blockValue < top.d_diff)
1384
                                    : (blockValue > top.d_diff));
1385
      break;
1386
    }
1387
  }
1388
}
1389
1390
void LinearEqualityModule::substitutePlusTimesConstant(ArithVar to, ArithVar from, const Rational& mult){
1391
  d_tableau.substitutePlusTimesConstant(to, from, mult, d_trackCallback);
1392
}
1393
void LinearEqualityModule::directlyAddToCoefficient(ArithVar row, ArithVar col, const Rational& mult){
1394
  d_tableau.directlyAddToCoefficient(row, col, mult, d_trackCallback);
1395
}
1396
1397
}  // namespace arith
1398
}  // namespace theory
1399
29322
}  // namespace cvc5