GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/linear_equality.cpp Lines: 368 791 46.5 %
Date: 2021-11-07 Branches: 620 3246 19.1 %

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