SMT-RAT  24.02
Toolbox for Strategic and Parallel Satisfiability-Modulo-Theories Solving
Tableau.tpp
Go to the documentation of this file.
1 /**
2  * @file Tableau.hpp
3  * @author Florian Corzilius <corzilius@cs.rwth-aachen.de>
4  * @since 2012-04-05
5  * @version 2014-10-01
6  */
7 
8 #pragma once
9 
10 #include "Tableau.h"
11 #include "TableauSettings.h"
12 
13 // #define DEBUG_METHODS_TABLEAU
14 // #define DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
15 // #define LRA_PEDANTIC_CORRECTNESS_CHECKS
16 // #define DEBUG_SOI_SIMPLEX
17 
18 namespace smtrat
19 {
20  namespace lra
21  {
22  template<class Settings, typename T1, typename T2>
23  Tableau<Settings,T1,T2>::Tableau( ModuleInput::iterator _defaultBoundPosition ):
24  mRowsCompressed( true ),
25  mWidth( 0 ),
26  mPivotingSteps( 0 ),
27  mMaxPivotsWithoutBlandsRule( 0 ),
28  mVarIDCounter( 1 ),
29  mDefaultBoundPosition( _defaultBoundPosition ),
30  mUnusedIDs(),
31  mVariableIdAllocator(),
32  mRows(),
33  mColumns(),
34  mNonActiveBasics(),
35  mConflictingRows(),
36  mCurDelta( 0 ),
37  mOriginalVars(),
38  mSlackVars(),
39  mConstraintToBound(),
40  mLearnedLowerBounds(),
41  mLearnedUpperBounds(),
42  mNewLearnedBounds()
43  {
44  mpEntries = new std::vector< TableauEntry<T1,T2> >();
45  mpEntries->push_back( TableauEntry<T1,T2>() );
46  mpTheta = new Value<T1>();
47  }
48 
49  template<class Settings, typename T1, typename T2>
50  Tableau<Settings,T1,T2>::~Tableau()
51  {
52  while( !mConstraintToBound.empty() )
53  {
54  std::vector< const Bound<T1,T2>* >* toDelete = mConstraintToBound.begin()->second;
55  mConstraintToBound.erase( mConstraintToBound.begin() );
56  if( toDelete != NULL ) delete toDelete;
57  }
58  while( !mOriginalVars.empty() )
59  {
60  Variable<T1,T2>* varToDelete = mOriginalVars.begin()->second;
61  mOriginalVars.erase( mOriginalVars.begin() );
62  delete varToDelete;
63  }
64  while( !mSlackVars.empty() )
65  {
66  Variable<T1,T2>* varToDelete = mSlackVars.begin()->second;
67  mSlackVars.erase( mSlackVars.begin() );
68  delete varToDelete;
69  }
70  delete mpEntries;
71  delete mpTheta;
72  }
73 
74  template<class Settings, typename T1, typename T2>
75  EntryID Tableau<Settings,T1,T2>::newTableauEntry( const T2& _content )
76  {
77  if( mUnusedIDs.empty() )
78  {
79  mpEntries->push_back( TableauEntry<T1,T2>( LAST_ENTRY_ID, LAST_ENTRY_ID, LAST_ENTRY_ID, LAST_ENTRY_ID, 0, 0, _content ) );
80  return ( ( mpEntries->size() ) - 1);
81  }
82  else
83  {
84  EntryID id = mUnusedIDs.top();
85  mUnusedIDs.pop();
86  (*mpEntries)[id].rContent() = _content;
87  return id;
88  }
89  }
90 
91  template<class Settings, typename T1, typename T2>
92  void Tableau<Settings,T1,T2>::removeEntry( EntryID _entryID )
93  {
94  TableauEntry<T1,T2>& entry = (*mpEntries)[_entryID];
95  const EntryID& up = entry.vNext( false );
96  const EntryID& down = entry.vNext( true );
97  if( up != LAST_ENTRY_ID )
98  {
99  (*mpEntries)[up].setVNext( true, down );
100  }
101  if( down != LAST_ENTRY_ID )
102  {
103  (*mpEntries)[down].setVNext( false, up );
104  }
105  else
106  {
107  entry.columnVar()->rStartEntry() = up;
108  }
109  const EntryID& left = entry.hNext( true );
110  const EntryID& right = entry.hNext( false );
111  if( left != LAST_ENTRY_ID )
112  {
113  (*mpEntries)[left].setHNext( false, right );
114  }
115  else
116  {
117  entry.rowVar()->rStartEntry() = right;
118  }
119  if( right != LAST_ENTRY_ID )
120  {
121  (*mpEntries)[right].setHNext( true, left );
122  }
123  --(entry.rowVar()->rSize());
124  --(entry.columnVar()->rSize());
125  mUnusedIDs.push( _entryID );
126  }
127 
128  template<class Settings, typename T1, typename T2>
129  void Tableau<Settings,T1,T2>::activateBound( const Bound<T1,T2>* _bound, const FormulaT& _formula )
130  {
131  _bound->pOrigins()->push_back( _formula );
132  const Variable<T1,T2>& var = _bound->variable();
133  if( !var.hasBound() && var.isBasic() && !var.isOriginal() )
134  activateBasicVar( _bound->pVariable() );
135  if( _bound->isUpperBound() )
136  {
137  if( *var.pSupremum() > *_bound )
138  {
139  _bound->pVariable()->setSupremum( _bound );
140  if( !(*var.pInfimum() > _bound->limit() && !_bound->deduced()) && !var.isBasic() && (*var.pSupremum() < var.assignment()) )
141  {
142  updateBasicAssignments( var.position(), Value<T1>( (*var.pSupremum()).limit() - var.assignment() ) );
143  _bound->pVariable()->rAssignment() = (*var.pSupremum()).limit();
144  }
145  }
146  }
147  if( _bound->isLowerBound() )
148  {
149  if( *var.pInfimum() < *_bound )
150  {
151  _bound->pVariable()->setInfimum( _bound );
152  if( !(*var.pSupremum() < _bound->limit() && !_bound->deduced()) && !var.isBasic() && (*var.pInfimum() > var.assignment()) )
153  {
154  updateBasicAssignments( var.position(), Value<T1>( (*var.pInfimum()).limit() - var.assignment() ) );
155  _bound->pVariable()->rAssignment() = (*var.pInfimum()).limit();
156  }
157  }
158  }
159  }
160 
161  template<class Settings, typename T1, typename T2>
162  Variable<T1,T2>* Tableau<Settings,T1,T2>::getVariable( const Poly& _lhs, T1& _factor, T1& _boundValue )
163  {
164  Variable<T1, T2>* result;
165  if( _lhs.nr_terms() == 1 || ( _lhs.nr_terms() == 2 && _lhs.has_constant_term() ) )
166  {
167  // TODO: do not store the expanded polynomial, but use the coefficient and coprimeCoefficients
168  const typename Poly::PolyType& expandedPoly = _lhs.polynomial();
169  auto term = expandedPoly.begin();
170  for( ; term != expandedPoly.end(); ++term )
171  if( !term->is_constant() ) break;
172  carl::Variable var = term->monomial()->begin()->first;
173  _factor = T1( term->coeff() ) * _lhs.coefficient();
174  _boundValue = T1( -_lhs.constant_part() )/_factor;
175  auto basicIter = mOriginalVars.find( var );
176  // constraint not found, add new nonbasic variable
177  if( basicIter == mOriginalVars.end() )
178  {
179  typename Poly::PolyType* varPoly = new typename Poly::PolyType( var );
180  result = newNonbasicVariable( varPoly, var.type() == carl::VariableType::VT_INT );
181  mOriginalVars.insert( std::pair<carl::Variable, Variable<T1, T2>*>( var, result ) );
182  }
183  else
184  {
185  result = basicIter->second;
186  }
187  }
188  else
189  {
190  T1 constantPart( _lhs.constant_part() );
191  bool negative = (_lhs.lterm().coeff() < typename Poly::CoeffType(T1( 0 )));
192  typename Poly::PolyType* linearPart;
193  if( negative )
194  linearPart = new typename Poly::PolyType( -_lhs + (Rational)constantPart );
195  else
196  linearPart = new typename Poly::PolyType( _lhs - (Rational)constantPart );
197  _factor = linearPart->coprime_factor();
198  assert( _factor > 0 );
199  constantPart *= _factor;
200  (*linearPart) *= _factor;
201 // linearPart->makeOrdered();
202  _boundValue = (negative ? constantPart : T1(-constantPart));
203  typename carl::FastPointerMap<typename Poly::PolyType, Variable<T1, T2>*>::iterator slackIter = mSlackVars.find( linearPart );
204  if( slackIter == mSlackVars.end() )
205  {
206  result = newBasicVariable( linearPart, _lhs.integer_valued() );
207  mSlackVars.insert( std::pair<const typename Poly::PolyType*, Variable<T1, T2>*>( linearPart, result ) );
208  }
209  else
210  {
211  delete linearPart;
212  result = slackIter->second;
213  }
214  if( negative )
215  _factor = T1(-_factor);
216  }
217  return result;
218  }
219 
220  template<class Settings, typename T1, typename T2>
221  Variable<T1,T2>* Tableau<Settings,T1,T2>::getObjectiveVariable( const Poly& _lhs )
222  {
223  return newBasicVariable( new typename Poly::PolyType( _lhs ), _lhs.integer_valued() );
224  }
225 
226  template<class Settings, typename T1, typename T2>
227  std::pair<const Bound<T1,T2>*, bool> Tableau<Settings,T1,T2>::newBound( const FormulaT& _constraint )
228  {
229  auto ctbIter = mConstraintToBound.find( _constraint );
230  if( ctbIter != mConstraintToBound.end() )
231  return std::make_pair( *ctbIter->second->begin(), false );
232  assert( _constraint.type() == carl::FormulaType::CONSTRAINT );
233  const ConstraintT& constraint = _constraint.constraint();
234  assert( constraint.is_consistent() == 2 );
235  T1 factor( 0 );
236  T1 boundValue( 0 );
237  Variable<T1, T2>* newVar = getVariable( constraint.lhs(), factor, boundValue );
238  bool negative = (factor < T1(0));
239  std::pair<const Bound<T1,T2>*, bool> result;
240  switch( constraint.relation() )
241  {
242  case carl::Relation::EQ:
243  {
244  // TODO: Take value from an allocator to assure the values are located close to each other in the memory.
245  Value<T1>* value = new Value<T1>( boundValue );
246  result = newVar->addEqualBound( value, mDefaultBoundPosition, _constraint );
247  std::vector< const Bound<T1,T2>* >* boundVector = new std::vector< const Bound<T1,T2>* >();
248  result.first->boundExists();
249  boundVector->push_back( result.first );
250  assert( mConstraintToBound.find( _constraint ) == mConstraintToBound.end() );
251  mConstraintToBound[_constraint] = boundVector;
252  break;
253  }
254  case carl::Relation::LEQ:
255  {
256  Value<T1>* value = new Value<T1>( boundValue );
257  result = negative ? newVar->addLowerBound( value, mDefaultBoundPosition, _constraint ) : newVar->addUpperBound( value, mDefaultBoundPosition, _constraint );
258  std::vector< const Bound<T1,T2>* >* boundVector = new std::vector< const Bound<T1,T2>* >();
259  boundVector->push_back( result.first );
260  assert( mConstraintToBound.find( _constraint ) == mConstraintToBound.end() );
261  mConstraintToBound[_constraint] = boundVector;
262  result.first->boundExists();
263  // create the complement
264  Value<T1>* vc = constraint.integer_valued() ? new Value<T1>( boundValue + (negative ? T1( -1 ) : T1( 1 )) ) : new Value<T1>( boundValue, (negative ? T1( -1 ) : T1( 1 )) );
265  FormulaT complConstr( _constraint.constraint().lhs(), carl::inverse( _constraint.constraint().relation() ) );
266  const Bound<T1,T2>* complement = negative ? newVar->addUpperBound( vc, mDefaultBoundPosition, complConstr ).first : newVar->addLowerBound( vc, mDefaultBoundPosition, complConstr ).first;
267  auto ctbInsertRes = mConstraintToBound.insert( std::make_pair( complConstr, nullptr ) );
268  if( ctbInsertRes.second )
269  {
270  std::vector< const Bound<T1,T2>* >* compBoundVector = new std::vector< const Bound<T1,T2>* >();
271  compBoundVector->push_back( complement );
272  ctbInsertRes.first->second = compBoundVector;
273  result.first->setComplement( complement );
274  complement->setComplement( result.first );
275  }
276  break;
277  }
278  case carl::Relation::GEQ:
279  {
280  Value<T1>* value = new Value<T1>( boundValue );
281  result = negative ? newVar->addUpperBound( value, mDefaultBoundPosition, _constraint ) : newVar->addLowerBound( value, mDefaultBoundPosition, _constraint );
282  std::vector< const Bound<T1,T2>* >* boundVector = new std::vector< const Bound<T1,T2>* >();
283  boundVector->push_back( result.first );
284  assert( mConstraintToBound.find( _constraint ) == mConstraintToBound.end() );
285  mConstraintToBound[_constraint] = boundVector;
286  result.first->boundExists();
287  // create the complement
288  Value<T1>* vc = constraint.integer_valued() ? new Value<T1>( boundValue + (negative ? T1( 1 ) : T1( -1 )) ) : new Value<T1>( boundValue, (negative ? T1( 1 ) : T1( -1 ) ) );
289  FormulaT complConstr( _constraint.constraint().lhs(), carl::inverse( _constraint.constraint().relation() ) );
290  const Bound<T1,T2>* complement = negative ? newVar->addLowerBound( vc, mDefaultBoundPosition, complConstr ).first : newVar->addUpperBound( vc, mDefaultBoundPosition, complConstr ).first;
291  auto ctbInsertRes = mConstraintToBound.insert( std::make_pair( complConstr, nullptr ) );
292  if( ctbInsertRes.second )
293  {
294  std::vector< const Bound<T1,T2>* >* compBoundVector = new std::vector< const Bound<T1,T2>* >();
295  compBoundVector->push_back( complement );
296  ctbInsertRes.first->second = compBoundVector;
297  result.first->setComplement( complement );
298  complement->setComplement( result.first );
299  }
300  break;
301  }
302  case carl::Relation::LESS:
303  {
304  Value<T1>* value = new Value<T1>( boundValue, (negative ? T1( 1 ) : T1( -1 ) ) );
305  result = negative ? newVar->addLowerBound( value, mDefaultBoundPosition, _constraint ) : newVar->addUpperBound( value, mDefaultBoundPosition, _constraint );
306  std::vector< const Bound<T1,T2>* >* boundVector = new std::vector< const Bound<T1,T2>* >();
307  boundVector->push_back( result.first );
308  assert( mConstraintToBound.find( _constraint ) == mConstraintToBound.end() );
309  mConstraintToBound[_constraint] = boundVector;
310  result.first->boundExists();
311  // create the complement
312  Value<T1>* vc = new Value<T1>( boundValue );
313  FormulaT complConstr( _constraint.constraint().lhs(), carl::inverse( _constraint.constraint().relation() ) );
314  const Bound<T1,T2>* complement = negative ? newVar->addUpperBound( vc, mDefaultBoundPosition, complConstr ).first : newVar->addLowerBound( vc, mDefaultBoundPosition, complConstr ).first;
315  auto ctbInsertRes = mConstraintToBound.insert( std::make_pair( complConstr, nullptr ) );
316  if( ctbInsertRes.second )
317  {
318  std::vector< const Bound<T1,T2>* >* compBoundVector = new std::vector< const Bound<T1,T2>* >();
319  compBoundVector->push_back( complement );
320  ctbInsertRes.first->second = compBoundVector;
321  result.first->setComplement( complement );
322  complement->setComplement( result.first );
323  }
324  break;
325  }
326  case carl::Relation::GREATER:
327  {
328  Value<T1>* value = new Value<T1>( boundValue, (negative ? T1( -1 ) : T1( 1 )) );
329  result = negative ? newVar->addUpperBound( value, mDefaultBoundPosition, _constraint ) : newVar->addLowerBound( value, mDefaultBoundPosition, _constraint );
330  std::vector< const Bound<T1,T2>* >* boundVector = new std::vector< const Bound<T1,T2>* >();
331  boundVector->push_back( result.first );
332  assert( mConstraintToBound.find( _constraint ) == mConstraintToBound.end() );
333  mConstraintToBound[_constraint] = boundVector;
334  result.first->boundExists();
335  // create the complement
336  Value<T1>* vc = new Value<T1>( boundValue );
337  FormulaT complConstr( _constraint.constraint().lhs(), carl::inverse( _constraint.constraint().relation() ) );
338  const Bound<T1,T2>* complement = negative ? newVar->addLowerBound( vc, mDefaultBoundPosition, complConstr ).first : newVar->addUpperBound( vc, mDefaultBoundPosition, complConstr ).first;
339  auto ctbInsertRes = mConstraintToBound.insert( std::make_pair( complConstr, nullptr ) );
340  if( ctbInsertRes.second )
341  {
342  std::vector< const Bound<T1,T2>* >* compBoundVector = new std::vector< const Bound<T1,T2>* >();
343  compBoundVector->push_back( complement );
344  ctbInsertRes.first->second = compBoundVector;
345  result.first->setComplement( complement );
346  complement->setComplement( result.first );
347  }
348  break;
349  }
350  case carl::Relation::NEQ:
351  {
352  FormulaT constraintLess = FormulaT( smtrat::ConstraintT( constraint.lhs(), carl::Relation::LESS ) );
353  Value<T1>* valueA = constraint.integer_valued() ? new Value<T1>( boundValue - T1( 1 ) ) : new Value<T1>( boundValue, (negative ? T1( 1 ) : T1( -1 ) ) );
354  auto resultLess = negative ? newVar->addLowerBound( valueA, mDefaultBoundPosition, constraintLess ) : newVar->addUpperBound( valueA, mDefaultBoundPosition, constraintLess );
355  auto ctbInsertRes = mConstraintToBound.insert( std::make_pair( constraintLess, nullptr ) );
356  if( ctbInsertRes.second )
357  {
358  std::vector< const Bound<T1,T2>* >* boundVector = new std::vector< const Bound<T1,T2>* >();
359  boundVector->push_back( resultLess.first );
360  ctbInsertRes.first->second = boundVector;
361  }
362 // if( !constraint.integer_valued() )
363  resultLess.first->setNeqRepresentation( _constraint );
364 
365  std::vector< const Bound<T1,T2>* >* boundVectorB = new std::vector< const Bound<T1,T2>* >();
366  boundVectorB->push_back( resultLess.first );
367 
368  FormulaT constraintLeq = FormulaT( smtrat::ConstraintT( constraint.lhs(), carl::Relation::LEQ ) );
369  Value<T1>* valueB = new Value<T1>( boundValue );
370  auto resultLeq = negative ? newVar->addLowerBound( valueB, mDefaultBoundPosition, constraintLeq ) : newVar->addUpperBound( valueB, mDefaultBoundPosition, constraintLeq );
371  ctbInsertRes = mConstraintToBound.insert( std::make_pair( constraintLeq, nullptr ) );
372  if( ctbInsertRes.second )
373  {
374  std::vector< const Bound<T1,T2>* >* boundVector = new std::vector< const Bound<T1,T2>* >();
375  boundVector->push_back( resultLeq.first );
376  ctbInsertRes.first->second = boundVector;
377  }
378  resultLeq.first->setNeqRepresentation( _constraint );
379 
380  boundVectorB->push_back( resultLeq.first );
381 
382  FormulaT constraintGeq = FormulaT( smtrat::ConstraintT( constraint.lhs(), carl::Relation::GEQ ) );
383  Value<T1>* valueC = new Value<T1>( boundValue );
384  auto resultGeq = negative ? newVar->addUpperBound( valueC, mDefaultBoundPosition, constraintGeq ) : newVar->addLowerBound( valueC, mDefaultBoundPosition, constraintGeq );
385  ctbInsertRes = mConstraintToBound.insert( std::make_pair( constraintGeq, nullptr ) );
386  if( ctbInsertRes.second )
387  {
388  std::vector< const Bound<T1,T2>* >* boundVector = new std::vector< const Bound<T1,T2>* >();
389  boundVector->push_back( resultGeq.first );
390  ctbInsertRes.first->second = boundVector;
391  }
392  resultGeq.first->setNeqRepresentation( _constraint );
393 
394  boundVectorB->push_back( resultGeq.first );
395 
396  FormulaT constraintGreater = FormulaT( smtrat::ConstraintT( constraint.lhs(), carl::Relation::GREATER ) );
397  Value<T1>* valueD = constraint.integer_valued() ? new Value<T1>( boundValue + T1( 1 ) ) : new Value<T1>( boundValue, (negative ? T1( -1 ) : T1( 1 )) );
398  auto resultGreater = negative ? newVar->addUpperBound( valueD, mDefaultBoundPosition, constraintGreater ) : newVar->addLowerBound( valueD, mDefaultBoundPosition, constraintGreater );
399  ctbInsertRes = mConstraintToBound.insert( std::make_pair( constraintGreater, nullptr ) );
400  if( ctbInsertRes.second )
401  {
402  std::vector< const Bound<T1,T2>* >* boundVector = new std::vector< const Bound<T1,T2>* >();
403  boundVector->push_back( resultGreater.first );
404  ctbInsertRes.first->second = boundVector;
405  }
406 // if( !constraint.integer_valued() )
407  resultGreater.first->setNeqRepresentation( _constraint );
408 
409  boundVectorB->push_back( resultGreater.first );
410  assert( mConstraintToBound.find( _constraint ) == mConstraintToBound.end() );
411  mConstraintToBound[_constraint] = boundVectorB;
412  resultLess.first->setComplement( resultGeq.first );
413  resultLeq.first->setComplement( resultGreater.first );
414  resultGeq.first->setComplement( resultLess.first );
415  resultGreater.first->setComplement( resultLeq.first );
416  result = resultGreater; //??
417  break;
418  }
419  }
420  return result;
421  }
422 
423  template<class Settings, typename T1, typename T2>
424  void Tableau<Settings,T1,T2>::removeBound( const FormulaT& _constraint )
425  {
426  auto iter = mConstraintToBound.find( _constraint );
427  if( iter != mConstraintToBound.end() )
428  {
429  std::vector< const Bound<T1,T2>* >* boundVector = iter->second;
430  Variable<T1, T2>* boundVar = boundVector->back()->pVariable();
431  if( _constraint.constraint().relation() == carl::Relation::NEQ )
432  {
433  assert( boundVector->size() == 4 );
434  for( const Bound<T1,T2>* bound : *boundVector )
435  {
436  bound->resetNeqRepresentation();
437  if( bound->markedAsDeleted() )
438  {
439  auto iterB = mConstraintToBound.find( bound->asConstraint() );
440  assert( iterB != mConstraintToBound.end() );
441  std::vector< const Bound<T1,T2>* >* boundVectorB = iter->second;
442  assert( boundVectorB->size() == 1 );
443  const Bound<T1,T2>* boundB = boundVectorB->back();
444  assert( !boundB->isActive() );
445  assert( boundVar == boundB->pVariable() );
446  boundVar->removeBound( boundB );
447  boundVectorB->pop_back();
448  assert( !boundB->isActive() );
449  assert( boundVar->pInfimum() != boundB );
450  assert( boundVar->pSupremum() != boundB );
451  delete boundB;
452  delete boundVectorB;
453  mConstraintToBound.erase( iterB );
454  }
455  }
456  boundVector->clear();
457  delete boundVector;
458  }
459  else
460  {
461  while( boundVector->size() > 1 )
462  {
463  const Bound<T1,T2>* toDel = boundVector->back();
464  boundVar->removeBound( toDel );
465  boundVector->pop_back();
466  delete toDel;
467  }
468  const Bound<T1,T2>* bound = boundVector->back();
469  assert(!bound->isActive());
470  if( !bound->neqRepresentation().is_true() )
471  {
472  bound->markAsDeleted();
473  }
474  else
475  {
476  assert( !bound->isActive() );
477  assert( boundVar->pInfimum() != bound );
478  assert( boundVar->pSupremum() != bound );
479  boundVar->removeBound( bound );
480  boundVector->pop_back();
481  delete bound;
482  delete boundVector;
483  }
484  }
485  mConstraintToBound.erase( iter );
486  if( !boundVar->isOriginal() && boundVar->isBasic() && boundVar->lowerbounds().size() == 1 && boundVar->upperbounds().size() == 1 )
487  {
488  deleteVariable( boundVar );
489  }
490  }
491  }
492 
493  template<class Settings, typename T1, typename T2>
494  void Tableau<Settings,T1,T2>::deleteVariable( Variable<T1, T2>* _variable, bool _optimizationVar )
495  {
496  assert( !_variable->isOriginal() && _variable->isBasic() && _variable->lowerbounds().size() == 1 && _variable->upperbounds().size() == 1 );
497  if( _variable->positionInNonActives() == mNonActiveBasics.end() )
498  {
499  deactivateBasicVar( _variable );
500  compressRows();
501  }
502  assert( !_variable->hasBound() );
503  mNonActiveBasics.erase( _variable->positionInNonActives() );
504  _variable->setPositionInNonActives( mNonActiveBasics.end() );
505  if( !_optimizationVar )
506  mSlackVars.erase( _variable->pExpression() );
507  assert( _variable->isBasic() );
508  mVariableIdAllocator.free( _variable->id() );
509  delete _variable;
510  }
511 
512  template<class Settings, typename T1, typename T2>
513  Variable<T1, T2>* Tableau<Settings,T1,T2>::newNonbasicVariable( const typename Poly::PolyType* _poly, bool _isInteger )
514  {
515  Variable<T1, T2>* var = new Variable<T1, T2>( mWidth++, _poly, mDefaultBoundPosition, _isInteger, mVariableIdAllocator.get() );
516  mColumns.push_back( var );
517  return var;
518  }
519 
520  template<class Settings, typename T1, typename T2>
521  Variable<T1, T2>* Tableau<Settings,T1,T2>::newBasicVariable( const typename Poly::PolyType* _poly, bool _isInteger )
522  {
523  mNonActiveBasics.emplace_front();
524  Variable<T1, T2>* var = new Variable<T1, T2>( mNonActiveBasics.begin(), _poly, mDefaultBoundPosition, _isInteger, mVariableIdAllocator.get() );
525  for( auto term = _poly->begin(); term != _poly->end(); ++term )
526  {
527  assert( !term->is_constant() );
528  assert( carl::is_integer( term->coeff() ) );
529  carl::Variable var = term->monomial()->begin()->first;
530  Variable<T1, T2>* nonBasic;
531  auto nonBasicIter = mOriginalVars.find( var );
532  if( mOriginalVars.end() == nonBasicIter )
533  {
534  typename Poly::PolyType* varPoly = new typename Poly::PolyType( var );
535  nonBasic = newNonbasicVariable( varPoly, var.type() == carl::VariableType::VT_INT );
536  mOriginalVars.insert( std::pair<carl::Variable, Variable<T1, T2>*>( var, nonBasic ) );
537  }
538  else
539  {
540  nonBasic = nonBasicIter->second;
541  }
542  mNonActiveBasics.front().emplace_back( nonBasic, T2( carl::get_num( term->coeff() ) ) );
543  }
544  return var;
545  }
546 
547  template<class Settings, typename T1, typename T2>
548  void Tableau<Settings,T1,T2>::activateBasicVar( Variable<T1, T2>* _var )
549  {
550  if( _var->positionInNonActives() == mNonActiveBasics.end() )
551  return;
552  assert( _var->isBasic() );
553  assert( !_var->isOriginal() );
554  assert( !_var->hasBound() );
555  compressRows();
556  std::map<size_t,T2> coeffs;
557  for( auto lravarCoeffPair = _var->positionInNonActives()->begin(); lravarCoeffPair != _var->positionInNonActives()->end(); ++lravarCoeffPair )
558  {
559  Variable<T1, T2>* lravar = lravarCoeffPair->first;
560  if( lravar->isBasic() )
561  {
562  if( lravar->positionInNonActives() != mNonActiveBasics.end() && !lravar->isOriginal() )
563  {
564  if( Settings::omit_division )
565  {
566  T2 l = carl::lcm( lravarCoeffPair->second, lravar->factor() );
567  assert( l > 0 );
568  if( lravarCoeffPair->second < 0 && lravar->factor() < 0 )
569  l *= T2( -1 );
570  T2 ca = carl::div( l, lravar->factor() );
571  T2 cb = carl::div( l, lravarCoeffPair->second );
572  _var->rFactor() *= cb;
573  for( auto iter = coeffs.begin(); iter != coeffs.end(); ++iter )
574  {
575  iter->second *= cb;
576  }
577  auto iterB = lravarCoeffPair;
578  ++iterB;
579  for( ; iterB != _var->positionInNonActives()->end(); ++iterB )
580  {
581  iterB->second *= cb;
582  }
583  for( auto lravarCoeffPairB = lravar->positionInNonActives()->begin(); lravarCoeffPairB != lravar->positionInNonActives()->end(); ++lravarCoeffPairB )
584  {
585  _var->positionInNonActives()->emplace_back( lravarCoeffPairB->first, ca*lravarCoeffPairB->second );
586  }
587  }
588  else
589  {
590  for( auto lravarCoeffPairB = lravar->positionInNonActives()->begin(); lravarCoeffPairB != lravar->positionInNonActives()->end(); ++lravarCoeffPairB )
591  {
592  _var->positionInNonActives()->emplace_back( lravarCoeffPairB->first, lravarCoeffPair->second*lravarCoeffPairB->second );
593  }
594  }
595  }
596  else
597  {
598  if( Settings::omit_division )
599  {
600  T2 l = carl::lcm( lravarCoeffPair->second, lravar->factor() );
601  assert( l > 0 );
602  if( lravarCoeffPair->second < 0 && lravar->factor() < 0 )
603  l *= T2( -1 );
604  T2 ca = carl::div( l, lravar->factor() );
605  T2 cb = carl::div( l, lravarCoeffPair->second );
606  _var->rFactor() *= cb;
607  for( auto iter = coeffs.begin(); iter != coeffs.end(); ++iter )
608  {
609  iter->second *= cb;
610  }
611  auto iterB = lravarCoeffPair;
612  ++iterB;
613  for( ; iterB != _var->positionInNonActives()->end(); ++iterB )
614  {
615  iterB->second *= cb;
616  }
617  Iterator rowIter = Iterator( lravar->startEntry(), mpEntries );
618  while( true )
619  {
620  coeffs[(*rowIter).columnVar()->position()] += ca*(*rowIter).content();
621  if( rowIter.hEnd( false ) ) break;
622  else rowIter.hMove( false );
623  }
624  }
625  else
626  {
627  Iterator rowIter = Iterator( lravar->startEntry(), mpEntries );
628  while( true )
629  {
630  coeffs[(*rowIter).columnVar()->position()] += lravarCoeffPair->second*(*rowIter).content();
631  if( rowIter.hEnd( false ) ) break;
632  else rowIter.hMove( false );
633  }
634  }
635  }
636  }
637  else
638  {
639  coeffs[lravar->position()] += lravarCoeffPair->second;
640  }
641  }
642  mNonActiveBasics.erase( _var->positionInNonActives() );
643  _var->setPositionInNonActives( mNonActiveBasics.end() );
644 
645  T2 g = carl::abs( _var->factor() );
646  for( auto iter = coeffs.begin(); iter != coeffs.end(); ++iter )
647  {
648  if( g == 1 ) break;
649  if( iter->second != T2( 0 ) )
650  carl::gcd_assign( g, iter->second );
651  }
652  if( !(g == 1) )
653  {
654  assert( g > 0 );
655  for( auto iter = coeffs.begin(); iter != coeffs.end(); ++iter )
656  {
657  if( iter->second != T2( 0 ) )
658  carl::div_assign( iter->second, g );
659  }
660  carl::div_assign( _var->rFactor(), g );
661  }
662 
663  _var->setPosition( mRows.size() );
664  mRows.push_back( _var );
665  _var->rAssignment() = Value<T1>( 0 );
666  EntryID lastInsertedEntry = LAST_ENTRY_ID;
667  _var->rSize() = 0;
668  for( const auto& coeff : coeffs )
669  {
670  if( coeff.second == T2( 0 ) )
671  continue;
672  ++(_var->rSize());
673  EntryID entryID = newTableauEntry( coeff.second );
674  TableauEntry<T1,T2>& entry = (*mpEntries)[entryID];
675  // Fix the position.
676  entry.setColumnVar( mColumns[coeff.first] );
677  entry.setRowVar( _var );
678  EntryID& columnStart = mColumns[coeff.first]->rStartEntry();
679  // Set it as column end.
680  if( columnStart != LAST_ENTRY_ID )
681  {
682  (*mpEntries)[columnStart].setVNext( true, entryID );
683  }
684  entry.setVNext( false, columnStart );
685  columnStart = entryID;
686  ++(mColumns[coeff.first]->rSize());
687  entry.setVNext( true, LAST_ENTRY_ID );
688  // Put it in the row.
689  if( lastInsertedEntry == LAST_ENTRY_ID )
690  {
691  _var->rStartEntry() = entryID;
692  entry.setHNext( true, lastInsertedEntry );
693  }
694  else
695  {
696  Iterator rowIter = Iterator( lastInsertedEntry, mpEntries );
697  (*rowIter).setHNext( false, entryID );
698  entry.setHNext( true, rowIter.entryID() );
699  }
700  // For now, the entry will be the rightmost in this row.
701  entry.setHNext( false, LAST_ENTRY_ID );
702  lastInsertedEntry = entryID;
703  _var->rAssignment() += mColumns[coeff.first]->assignment() * coeff.second;
704  }
705  if( Settings::omit_division )
706  {
707  _var->rAssignment() /= _var->factor();
708  }
709  assert( checkCorrectness() == mRows.size() );
710  assert( _var->positionInNonActives() == mNonActiveBasics.end() );
711  }
712 
713  template<class Settings, typename T1, typename T2>
714  void Tableau<Settings,T1,T2>::deactivateBasicVar( Variable<T1, T2>* _var )
715  {
716  assert( _var->isBasic() );
717  assert( !_var->isOriginal() );
718  if( Settings::pivot_into_local_conflict )
719  {
720  auto crIter = mConflictingRows.begin();
721  for( ; crIter != mConflictingRows.end(); ++crIter )
722  if( (*crIter) == _var ) break;
723  if( crIter != mConflictingRows.end() )
724  {
725  mConflictingRows.erase( crIter );
726  }
727  }
728  mNonActiveBasics.emplace_front();
729  EntryID entryToRemove = _var->startEntry();
730  while( entryToRemove != LAST_ENTRY_ID )
731  {
732  TableauEntry<T1,T2>& entry = (*mpEntries)[entryToRemove];
733  Variable<T1, T2>* nbVar = entry.columnVar();
734  mNonActiveBasics.front().emplace_back( nbVar, entry.content() );
735  EntryID up = entry.vNext( false );
736  EntryID down = entry.vNext( true );
737  if( up != LAST_ENTRY_ID )
738  {
739  (*mpEntries)[up].setVNext( true, down );
740  }
741  if( down != LAST_ENTRY_ID )
742  {
743  (*mpEntries)[down].setVNext( false, up );
744  }
745  else
746  {
747  nbVar->rStartEntry() = up;
748  }
749  --(nbVar->rSize());
750  entry.setRowVar( NULL ); //Not necessary but cleaner.
751  mUnusedIDs.push( entryToRemove );
752  entryToRemove = entry.hNext( false );
753  }
754  mRows[_var->position()] = NULL;
755  _var->rStartEntry() = LAST_ENTRY_ID;
756  _var->rSize() = 0;
757  _var->rPosition() = 0;
758  _var->setPositionInNonActives( mNonActiveBasics.begin() );
759  mRowsCompressed = false;
760  }
761 
762  template<class Settings, typename T1, typename T2>
763  void Tableau<Settings,T1,T2>::storeAssignment()
764  {
765  for( Variable<T1, T2>* basicVar : mRows )
766  {
767  if( basicVar != NULL )
768  basicVar->storeAssignment();
769  }
770  for( Variable<T1, T2>* nonbasicVar : mColumns )
771  nonbasicVar->storeAssignment();
772  }
773 
774  template<class Settings, typename T1, typename T2>
775  void Tableau<Settings,T1,T2>::resetAssignment()
776  {
777  for( Variable<T1, T2>* basicVar : mRows )
778  {
779  if( basicVar != NULL )
780  {
781  basicVar->resetAssignment();
782  }
783  }
784  for( Variable<T1, T2>* nonbasicVar : mColumns )
785  nonbasicVar->resetAssignment();
786  }
787 
788  template<class Settings, typename T1, typename T2>
789  RationalAssignment Tableau<Settings,T1,T2>::getRationalAssignment() const
790  {
791  T1 minDelta = -1;
792  mCurDelta = T1(0);
793  // For all non-basic variables find the minimum of all (c2-c1)/(k1-k2), where ..
794  for( auto originalVar : mColumns )
795  {
796  adaptDelta( *originalVar, false, minDelta );
797  adaptDelta( *originalVar, true, minDelta );
798  }
799  // For all basic variables find the minimum of all (c2-c1)/(k1-k2), where ..
800  for( auto var : mRows )
801  {
802  if( var != NULL )
803  {
804  adaptDelta( *var, false, minDelta );
805  adaptDelta( *var, true, minDelta );
806  }
807  }
808  mCurDelta = minDelta < 0 ? T1(1) : minDelta;
809  RationalAssignment result;
810  // Calculate the rational assignment of all original variables.
811  for( auto var : mColumns )
812  {
813  if( var->isOriginal() )
814  {
815  T1 value = var->assignment().mainPart();
816  value += (var->assignment().deltaPart() * mCurDelta);
817  result.insert( std::pair<const carl::Variable,T1>( var->expression().single_variable(), value ) );
818  }
819  }
820  for( auto var : mRows )
821  {
822  if( var != NULL && var->isOriginal() )
823  {
824  T1 value = var->assignment().mainPart();
825  value += (var->assignment().deltaPart() * mCurDelta);
826  result.insert( std::pair<const carl::Variable,T1>( var->expression().single_variable(), value ) );
827  }
828  }
829  return result;
830  }
831 
832  template<class Settings, typename T1, typename T2>
833  void Tableau<Settings,T1,T2>::adaptDelta( const Variable<T1,T2>& _variable, bool _upperBound, T1& _minDelta ) const
834  {
835  const Value<T1>& assValue = _variable.assignment();
836  const Bound<T1,T2>& bound = _upperBound ? _variable.supremum() : _variable.infimum();
837  if( !bound.isInfinite() )
838  {
839  // .. the bound limit is c1+k1*delta, the variable assignment is c2+k2*delta, then ..
840  if( (_upperBound && bound.limit().mainPart() > assValue.mainPart() && bound.limit().deltaPart() < assValue.deltaPart()) // .. c1>c2 and k1<k2
841  || (!_upperBound && bound.limit().mainPart() < assValue.mainPart() && bound.limit().deltaPart() > assValue.deltaPart()) ) // .. c1<c2 and k1>k2
842  {
843  if( _upperBound )
844  mCurDelta = ( bound.limit().mainPart() - assValue.mainPart() ) / ( assValue.deltaPart() - bound.limit().deltaPart() );
845  else
846  mCurDelta = ( assValue.mainPart() - bound.limit().mainPart() ) / ( bound.limit().deltaPart() - assValue.deltaPart() );
847  if( _minDelta < 0 || mCurDelta < _minDelta )
848  _minDelta = mCurDelta;
849  }
850  }
851  }
852 
853  template<class Settings, typename T1, typename T2>
854  void Tableau<Settings,T1,T2>::compressRows()
855  {
856  if( mRowsCompressed ) return;
857  std::deque<size_t> emptyPositions;
858  size_t curPos = 0;
859  while( curPos < mRows.size() )
860  {
861  if( mRows[curPos] == NULL )
862  {
863  emptyPositions.push_back( curPos );
864  }
865  else if( !emptyPositions.empty() )
866  {
867  size_t emptyPos = emptyPositions.front();
868  emptyPositions.pop_front();
869  mRows[emptyPos] = mRows[curPos];
870  mRows[emptyPos]->rPosition() = emptyPos;
871  emptyPositions.push_back( curPos );
872  }
873  ++curPos;
874  }
875  mRows.resize( mRows.size() - emptyPositions.size() );
876  mRowsCompressed = true;
877  }
878 
879  template<class Settings, typename T1, typename T2>
880  bool Tableau<Settings,T1,T2>::usedBlandsRule(){
881  return !(mPivotingSteps < mMaxPivotsWithoutBlandsRule);
882  }
883 
884  template<class Settings, typename T1, typename T2>
885  std::vector<T2> Tableau<Settings,T1,T2>::getInfeasibilityRow(){
886  std::vector<T1> result;
887 
888  for( Variable<T1, T2>* nonbasicVar : mColumns ){
889  T2 runningSum = 0;
890  Iterator columnIter = Iterator( nonbasicVar->startEntry(), mpEntries ); // startEntry points to last row -> has to iterate upwards
891 
892  if(nonbasicVar->size() == 0){
893  result.push_back(T1 (0.0));
894  continue;
895  }
896  while(true){
897  const Variable<T1, T2>& basicVar = *(*columnIter).rowVar();
898  assert( (*columnIter).rowVar() != NULL );
899  T2 factor = 1;
900  if(Settings::omit_division){
901  factor = basicVar.factor();
902  }
903  // check direction:
904  if(basicVar.supremum() < basicVar.assignment()){
905  runningSum += (*columnIter).rContent()/factor;
906  }
907  if(basicVar.infimum() > basicVar.assignment()){
908  runningSum -= (*columnIter).rContent()/factor;
909  }
910  if( columnIter.vEnd( false ) )
911  break;
912  else
913  columnIter.vMove( false );
914  }
915  result.push_back(runningSum);
916  }
917  return result;
918  }
919 
920  template<class Settings, typename T1, typename T2>
921  void Tableau<Settings,T1,T2>::setInfeasibilityRow(){
922  mInfeasibilityRow = getInfeasibilityRow();
923  }
924 
925  /**
926  * Compute the sum of violations over all variables.
927  */
928  template<class Settings, typename T1, typename T2>
929  Value<T1> Tableau<Settings,T1,T2>::violationSum(){
930  Value<T1> sum (0.0);
931  // sufficient to iterate over basic variables, nonbasic variables fulfill their bounds.
932  for( Variable<T1, T2>* basicVar : mRows ){
933  if(basicVar->supremum() < basicVar->assignment()){ //supremum (s1,s2) , assignment (a1,a2) -> s1 < a1 || ( s1 == a1 && s2 < a2 )
934  auto diff = basicVar->assignment() - basicVar->supremum().limit();
935  sum += diff;
936  }
937  if(basicVar->infimum() > basicVar->assignment()){
938  auto diff = basicVar->infimum().limit() - basicVar->assignment();
939  sum += diff;
940  }
941  }
942 
943  // Repeat for columns -> needed as VioSum is also called from dVioSum, where a nonbasic viarable might be set out of bounds
944  for( Variable<T1, T2>* nonbasicVar : mColumns ){
945  if(nonbasicVar->supremum() < nonbasicVar->assignment()){ //supremum (s1,s2) , assignment (a1,a2) -> s1 < a1 || ( s1 == a1 && s2 < a2 )
946  auto diff = nonbasicVar->assignment() - nonbasicVar->supremum().limit();
947  sum += diff;
948  }
949  if(nonbasicVar->infimum() > nonbasicVar->assignment()){
950  auto diff = nonbasicVar->infimum().limit() - nonbasicVar->assignment();
951  sum += diff;
952  }
953  }
954 
955  return sum;
956  }
957 
958  /**
959  * Naive implementation of computing all differences in sum when updating nVar with update. Todo: more efficient implementation.
960  */
961  template<class Settings, typename T1, typename T2>
962  Value<T1> Tableau<Settings,T1,T2>::dViolationSum(const Variable<T1, T2>* nVar, const Value<T1>& update){
963  // iterate over rows
964  assert(!nVar->isBasic()); // cant update basic variable
965 
966 
967  #ifdef DEBUG_SOI_SIMPLEX
968  SMTRAT_LOG_DEBUG("smtrat","Update: " << nVar->position() << " by " << update);
969  std::cout << "lower bound: " << nVar->infimum() << " assignment: " << nVar->assignment() << " upper bound: " << nVar->supremum() << '\n';
970  #endif
971 
972  #ifdef DEBUG_SOI_SIMPLEX
973  Value<T1> previousAssignment = nVar->assignment();
974  #endif
975  Value<T1> savedSum = violationSum();
976  EntryID column_index = nVar->startEntry();
977  updateNonbasicVariable(column_index, update);
978  Value<T1> updated_dVio = violationSum() - savedSum;
979  updateNonbasicVariable(column_index, update*T1(-1));
980 
981  #ifdef DEBUG_SOI_SIMPLEX
982  assert(previousAssignment == nVar->assignment()); // assert that update was reverted
983  #endif
984 
985  return updated_dVio;
986  }
987 
988  /**
989  * Checks if a conflicting pair exists. In a conflict is found, the cell and false is returned. If not, LAST_ENTRY_ID and true is returned.
990  */
991  template<class Settings, typename T1, typename T2>
992  std::pair<EntryID,bool> Tableau<Settings,T1,T2>::hasConflict(){
993  for( const Variable<T1, T2>* basicVar : mRows ){
994  assert( basicVar != NULL );
995  const Variable<T1,T2>& bVar = *basicVar;
996  Value<T1> thetaB = Value<T1>( 0 );
997  bool upperBoundViolated = false;
998  bool lowerBoundViolated = false;
999  if( bVar.supremum() < bVar.assignment() )
1000  {
1001  thetaB = bVar.supremum().limit() - bVar.assignment();
1002  upperBoundViolated = true;
1003  }
1004  else if( bVar.infimum() > bVar.assignment() )
1005  {
1006  thetaB = bVar.infimum().limit() - bVar.assignment();
1007  lowerBoundViolated = true;
1008  }
1009  if( upperBoundViolated || lowerBoundViolated )
1010  {
1011  EntryID result = isSuitableConflictDetection( bVar, lowerBoundViolated );
1012  if( result == LAST_ENTRY_ID )
1013  {
1014  // Found a conflicting row.
1015  return std::pair<EntryID,bool>( bVar.startEntry(), false );
1016  }
1017  }
1018  }
1019  // no conflict found
1020  return std::pair<EntryID,bool>( LAST_ENTRY_ID, true );
1021  }
1022 
1023  /**
1024  * Function to update a nonbasic variable, _pivotingElement needs to be in the column of the nonbasic variable.
1025  */
1026  template<class Settings, typename T1, typename T2>
1027  void Tableau<Settings,T1,T2>::updateNonbasicVariable(EntryID _pivotingElement){
1028 
1029  Variable<T1, T2>* columnVar = (*mpEntries)[_pivotingElement].columnVar();
1030  assert(!columnVar->isBasic()); // dummy check
1031 
1032  SMTRAT_LOG_DEBUG("smtrat", "Updated from" << columnVar -> rAssignment());
1033  columnVar -> rAssignment() = columnVar -> assignment() + (*mpTheta);
1034  SMTRAT_LOG_DEBUG("smtrat", "To " << columnVar -> rAssignment());
1035 
1036  assert(columnVar-> assignment() >= columnVar -> rAssignment() && columnVar-> infimum() <= columnVar -> rAssignment());
1037 
1038  Iterator rowIter = Iterator( columnVar->startEntry(), mpEntries );
1039  Variable<T1, T2>* rowVar;
1040 
1041  while(true){
1042  rowVar = (*mpEntries)[rowIter.entryID()].rowVar();
1043 
1044  if(Settings::omit_division){
1045  rowVar -> rAssignment() += (*mpTheta) * (*rowIter).rContent()/rowVar->factor(); // update nonBasic variable
1046  }
1047  else{
1048  rowVar -> rAssignment() += (*mpTheta) * (*rowIter).rContent();
1049  }
1050 
1051  if( rowIter.vEnd( false ) )
1052  {
1053  break;
1054  }
1055  else{
1056  rowIter.vMove( false );
1057  }
1058  }
1059 
1060  onlyUpdate = false;
1061  }
1062 
1063  /**
1064  * Function to update a nonbasic variable, _pivotingElement needs to be in the column of the nonbasic variable.
1065  * Used in the dViolationSum function to update die nonbasic variable.
1066  * For this reason, the assertion checks on the variable assignments are omitted (the nonbasic variable might be pivoted into the basic).
1067  */
1068  template<class Settings, typename T1, typename T2>
1069  void Tableau<Settings,T1,T2>::updateNonbasicVariable(EntryID _pivotingElement, Value<T1>update){
1070 
1071  Variable<T1, T2>* columnVar = (*mpEntries)[_pivotingElement].columnVar();
1072  assert(!columnVar->isBasic()); // dummy check
1073 
1074  SMTRAT_LOG_DEBUG("smtrat", "Updated from" << columnVar -> rAssignment());
1075  columnVar -> rAssignment() = columnVar -> assignment() + update;
1076  SMTRAT_LOG_DEBUG("smtrat", "To " << columnVar -> rAssignment());
1077 
1078  Iterator rowIter = Iterator( columnVar->startEntry(), mpEntries );
1079  Variable<T1, T2>* rowVar;
1080 
1081  while(true){
1082  rowVar = (*mpEntries)[rowIter.entryID()].rowVar();
1083 
1084  if(Settings::omit_division){
1085  rowVar -> rAssignment() += update * (*rowIter).rContent()/rowVar->factor(); // update nonBasic variable
1086  }
1087  else{
1088  rowVar -> rAssignment() += update * (*rowIter).rContent();
1089  }
1090 
1091  if( rowIter.vEnd( false ) )
1092  {
1093  break;
1094  }
1095  else{
1096  rowIter.vMove( false );
1097  }
1098  }
1099  }
1100 
1101  template<class Settings, typename T1, typename T2>
1102  void Tableau<Settings,T1,T2>::computeLeavingCandidates(const std::size_t i, std::vector< std::pair< Value<T1>, Variable<T1,T2>* > >& leaving_candidates){
1103  const Variable<T1,T2>& nVar = *mColumns[i];
1104 
1105  // consider possible leaving variables
1106 
1107  // consider all basic vars with coeff != 0
1108  Iterator rowIter = Iterator( nVar.startEntry(), mpEntries );
1109  EntryID row_id;
1110  Value<T1> cell_content;
1111  // first element is head of row
1112  while(true){
1113  row_id = rowIter.entryID();
1114  cell_content = (*mpEntries)[row_id].content();
1115  if(cell_content != T1(0.0)){
1116  Variable<T1,T2>& bVar = *(*mpEntries)[rowIter.entryID()].rowVar();
1117  Value<T1> assignment = bVar.assignment();
1118 
1119  // computed delta values
1120  if(! bVar.supremum().isInfinite()){
1121  T1 content = 1/(*mpEntries)[row_id].content();
1122  if(Settings::omit_division){
1123  content = 1/( (*mpEntries)[row_id].content() / bVar.factor() );
1124  }
1125 
1126  Value<T1> needed_update(bVar.supremum().limit());
1127  needed_update = (needed_update - assignment) * content;
1128 
1129  leaving_candidates.emplace_back(needed_update, &bVar);// lower bound
1130  }
1131 
1132  if(! bVar.infimum().isInfinite()){
1133  T1 content = 1/(*mpEntries)[row_id].content();
1134 
1135  if(Settings::omit_division){
1136  content = 1/( (*mpEntries)[row_id].content() / bVar.factor() );
1137  }
1138 
1139  Value<T1> needed_update(bVar.infimum().limit());
1140  needed_update = (needed_update - assignment) * content;
1141 
1142  leaving_candidates.emplace_back(needed_update, &bVar);// upper bound
1143  }
1144 
1145  }
1146  // moves upwards
1147  if( rowIter.vEnd( false ) )
1148  {
1149  break;
1150  }
1151  else{
1152  rowIter.vMove( false );
1153  }
1154  }
1155 
1156  // nonbasic variable is flexible
1157  // to test whether nonbasic variable should only be updated
1158  // Transfer from paper: variable shall just be updated to its bounds, the coefficient is one (no pivot step)
1159 
1160  // computed delta values
1161  Value<T1> assignment = nVar.assignment();
1162  if(! nVar.supremum().isInfinite()){
1163  Value<T1> needed_update(nVar.supremum().limit());
1164  needed_update = needed_update - assignment; // content of cell is 1 -> only update operation
1165  // one could skip those 0 updates, BUT for the faster computation of dVio they are important: Influence the slope of the error function
1166  leaving_candidates.emplace_back(needed_update, mColumns[i]);// lower bound
1167  }
1168 
1169  if(! nVar.infimum().isInfinite()){
1170  Value<T1> needed_update(nVar.infimum().limit());
1171  needed_update = needed_update - assignment; // content of cell is 1 -> only update operation
1172  leaving_candidates.emplace_back(needed_update, mColumns[i]);// upper bound
1173  }
1174  }
1175 
1176  template<class Settings, typename T1, typename T2>
1177  std::map<Value<T1>, Value<T1>> Tableau<Settings,T1,T2>::compute_dVio(const std::vector< std::pair< Value<T1>, Variable<T1,T2>* > >& candidates, const Variable<T1,T2>& nVar, bool positive){
1178  std::map<Value<T1>, Value<T1>> dVio_map; // add 0 dVIo by default for 0 update
1179  dVio_map.insert(std::pair<Value<T1>, Value<T1>> (Value<T1>(0), Value<T1>(0)));
1180  Value<T1> delta_old = Value<T1>(0.0);
1181  Value<T1> delta_new;
1182  Value<T1> new_dVio;
1183  T1 beta = mInfeasibilityRow[nVar.position()];
1184  std::size_t candidate_index = 0;
1185  if(candidates.empty())
1186  return dVio_map;
1187  // beta_0 has to be adjusted for the delta_0 values
1188  while(candidates[candidate_index].first == Value<T1>(0.0) && candidate_index < candidates.size()){
1189  int state;
1190  Variable<T1,T2>& updated_var = *(candidates[candidate_index].second);
1191  Iterator rowIter = Iterator( nVar.startEntry(), mpEntries );
1192  bool was_set = false;
1193  T1 coeff;
1194  // if updated var is not basic, then the current factor is 1
1195  if(updated_var.isBasic()){
1196  while(true){
1197  Variable<T1,T2>& bVar = *((*rowIter).rowVar());
1198  if(bVar.expression() == updated_var.expression()){
1199  coeff = (*mpEntries)[rowIter.entryID()].content();
1200  if(Settings::omit_division){
1201  coeff /= bVar.factor();
1202  }
1203  was_set = true;
1204  break;
1205  }
1206 
1207  // moves upwards
1208  if( rowIter.vEnd( false ) )
1209  {
1210  break;
1211  }
1212  else{
1213  rowIter.vMove( false );
1214  }
1215  }
1216  assert(was_set);
1217  }
1218  else{
1219  coeff = 1;
1220  }
1221 
1222  if( (!positive && coeff > 0) || (positive && coeff < 0) ){
1223  // decreasing part
1224  if(!updated_var.infimum().isInfinite() && !updated_var.supremum().isInfinite() && updated_var.infimum().limit() == updated_var.supremum().limit()){
1225  state = -1; // if one bound is set to directly the upper or lower bound -> the both bounds concide
1226  candidate_index++;
1227  }
1228  else{
1229  if(!updated_var.supremum().isInfinite() && updated_var.supremum() == updated_var.assignment()){
1230  // enters bounds
1231  state = 0;
1232  }
1233  else if(!updated_var.infimum().isInfinite() && updated_var.infimum() == updated_var.assignment()){
1234  // leaves bounds
1235  state = -1;
1236  }
1237  else{
1238  assert(false);
1239  }
1240  }
1241  }
1242  else{
1243  // increasing part
1244  if(!updated_var.infimum().isInfinite() && !updated_var.supremum().isInfinite() && updated_var.infimum().limit() == updated_var.supremum().limit()){
1245  state = 1;
1246  candidate_index++;
1247  }
1248  else{
1249  if(!updated_var.supremum().isInfinite() && updated_var.supremum() == updated_var.assignment()){
1250  // leaves bounds
1251  state = 1;
1252  }
1253  else if(!updated_var.infimum().isInfinite() && updated_var.infimum() == updated_var.assignment()){
1254  // enters bounds
1255  state = 0;
1256  }
1257  else{
1258  assert(false);
1259  }
1260  }
1261  }
1262  beta += state * coeff;
1263  candidate_index++;
1264  }
1265  while(candidate_index < candidates.size()){
1266  delta_new = candidates[candidate_index].first;
1267  // though the paper assumes strict greater, the equal values are handles in the do-while part
1268  SMTRAT_LOG_DEBUG("smtrat", "beta" << beta << " delta_new " << delta_new << " delta_old " << delta_old);
1269  new_dVio += (delta_new - delta_old) * beta;
1270  T1 summed_slopes;
1271 
1272  do{ // handle the k_i elements (elements with same coefficient)
1273  Variable<T1,T2>& updated_var = *(candidates[candidate_index].second);
1274  T1 coeff;
1275  // if updated var is not basic, then the current factor is 1
1276  if(updated_var.isBasic()){
1277  // find the coefficient between the nonbasic var and the var to be updated
1278  Iterator rowIter = Iterator( nVar.startEntry(), mpEntries );
1279  bool was_set = false;
1280  while(true){
1281  Variable<T1,T2>& bVar = *((*rowIter).rowVar());
1282  if(bVar.expression() == updated_var.expression()){
1283  coeff = (*mpEntries)[rowIter.entryID()].content();
1284  if(Settings::omit_division){
1285  coeff /= bVar.factor();
1286  }
1287  was_set = true;
1288  break;
1289  }
1290 
1291  // moves upwards
1292  if( rowIter.vEnd( false ) )
1293  {
1294  break;
1295  }
1296  else{
1297  rowIter.vMove( false );
1298  }
1299  }
1300  assert(was_set);
1301  }
1302  else{
1303  coeff = 1; // var is non-basic
1304  }
1305  int state;
1306  if(delta_new * coeff < 0){ // next to state the difference between the new direction and the previous direction is denoted for better understanding
1307  // decreasing direction
1308  if(!updated_var.infimum().isInfinite() && !updated_var.supremum().isInfinite() && updated_var.infimum().limit() == updated_var.supremum().limit()){
1309  state = -1; // if one bound is set to directly the upper or lower bound -> the both bounds concide
1310  candidate_index++;
1311  }
1312  if(!updated_var.supremum().isInfinite() && delta_new * coeff + updated_var.assignment() == updated_var.supremum().limit()){
1313  state = -1; // to supremum -> 0 - 1
1314  }
1315  else if(!updated_var.infimum().isInfinite() && delta_new * coeff + updated_var.assignment() == updated_var.infimum().limit()){
1316  state = -1; // to infimum -> -1 - 0
1317  }
1318  else{
1319  assert(false);
1320  }
1321  }
1322  else if(delta_new * coeff > 0){
1323  // increasing direction
1324  if(!updated_var.infimum().isInfinite() && !updated_var.supremum().isInfinite() && updated_var.infimum().limit() == updated_var.supremum().limit()){
1325  state = 1;
1326  candidate_index++;
1327  }
1328  if(!updated_var.supremum().isInfinite() && delta_new * coeff + updated_var.assignment() == updated_var.supremum().limit()){
1329  state = 1; // to supremum -> 1 - 0
1330  }
1331  else if(!updated_var.infimum().isInfinite() && delta_new * coeff + updated_var.assignment() == updated_var.infimum().limit()){
1332  state = 1; // to infimum -> 0 - (-1)
1333  }
1334  else{
1335  assert(false); // not allowed to occur
1336  }
1337  }
1338  else{
1339  assert(false); // not allowed to occur
1340  }
1341  SMTRAT_LOG_DEBUG("smtrat", "added " << state * coeff << " to " << summed_slopes);
1342 
1343  summed_slopes += state * coeff;
1344  SMTRAT_LOG_DEBUG("smtrat", "summed_slopes" << summed_slopes);
1345 
1346  candidate_index++;
1347  if(candidate_index >= candidates.size()){
1348  break;
1349  }
1350 
1351  }
1352  while(candidates[candidate_index].first == delta_new);
1353  delta_old = delta_new;
1354  beta += summed_slopes;
1355  assert(dVio_map.find(delta_new) == dVio_map.end());
1356  dVio_map.insert(std::pair<Value<T1>, Value<T1>> (delta_new, new_dVio));
1357  }
1358  return dVio_map;
1359  }
1360 
1361  template<class Settings, typename T1, typename T2>
1362  std::pair<EntryID,bool> Tableau<Settings,T1,T2>::nextPivotingElementInfeasibilities(){ // entry: position in tableau to be pivoted, bool: true if NO conflict found
1363  if( Settings::use_pivoting_strategy && mPivotingSteps < mMaxPivotsWithoutBlandsRule )
1364  {
1365  #ifdef DEBUG_SOI_SIMPLEX
1366  assert(mInfeasibilityRow == getInfeasibilityRow());
1367  #endif
1368  // leaving rule: minimizes d_Violation
1369  std::vector< std::tuple <std::size_t, Value<T1>, Variable<T1,T2>* > > update_candidates;
1370 
1371  // check conflicts
1372  std::pair<EntryID,bool> conflictPair = hasConflict();
1373  if(!conflictPair.second){
1374  SMTRAT_LOG_DEBUG("smtrat", "conflict found");
1375  // reset values for assertions
1376  mOldVioSum = Value<T1>(-1);
1377  mOld_dVioSum = Value<T1>(1);
1378  // conflict found: return conflict pair
1379  return conflictPair;
1380  }
1381 
1382  SMTRAT_LOG_DEBUG("smtrat", "Vio sum" << violationSum() << ", oldVioSum " << mOldVioSum << ", old_dVio" << mOld_dVioSum);
1383 
1384  #ifdef DEBUG_SOI_SIMPLEX
1385  mOldVioSum = violationSum();
1386  #endif
1387 
1388  #ifdef DEBUG_SOI_SIMPLEX
1389  std::stringstream s;
1390  for(auto e : mInfeasibilityRow){
1391  s << e << " ";
1392  }
1393  SMTRAT_LOG_DEBUG("smtrat", "mInfeasibilityRow row: "<< s.str());
1394  s.clear();
1395  s.str(""); //clear content
1396  #endif
1397 
1398  // for all flexible variables
1399  SMTRAT_LOG_DEBUG("smtrat", "Iteration over "<<mColumns.size()<<" columns");
1400 
1401  bool found_improvement;
1402  for(std::size_t column = 0; column < mColumns.size(); column++){
1403  assert( mColumns[column] != NULL ); // non_basic variable must exist
1404 
1405  const Variable<T1,T2>& nVar = *mColumns[column];
1406  assert(!nVar.isBasic());
1407  // check if nonbasis variable is flexible
1408  bool isFlexible = false;
1409  if(mInfeasibilityRow[column] > 0 && nVar.infimum() < nVar.assignment()){
1410  isFlexible = true;
1411  }
1412  else if (mInfeasibilityRow[column] < 0 && nVar.supremum() > nVar.assignment()){
1413  isFlexible = true;
1414  }
1415  if(! isFlexible){
1416  // nonbasic variable is not flexible
1417  SMTRAT_LOG_DEBUG("smtrat", "column " << nVar.position() << " is not flexible");
1418  continue;
1419  }
1420 
1421  SMTRAT_LOG_DEBUG("smtrat", "column "<< nVar.position() <<" is flexible");
1422 
1423  std::vector< std::pair< Value<T1>, Variable<T1,T2>* > > leaving_candidates;
1424  computeLeavingCandidates(column, leaving_candidates); // function to generate leaving candidates for the given column
1425  std::stable_sort( leaving_candidates.begin(), leaving_candidates.end() );// sorts by default by first element; IMPORTANT to use stable_sort -> neighbouring elements need to remain neighboured
1426  std::vector< std::pair< Value<T1>, Variable<T1,T2>* > > positive_candidates;
1427  std::vector< std::pair< Value<T1>, Variable<T1,T2>* > > negative_candidates;
1428 
1429  for(auto element : leaving_candidates){
1430  if(element.first > Value<T1>(0.0)){
1431  positive_candidates.push_back(element);
1432  }
1433  else if(element.first < Value<T1>(0.0)){
1434  negative_candidates.push_back(element);
1435  }
1436  else{
1437  positive_candidates.push_back(element);
1438  negative_candidates.push_back(element);
1439  }
1440  }
1441  std::reverse(negative_candidates.begin(),negative_candidates.end());
1442 
1443  #ifdef DEBUG_SOI_SIMPLEX
1444  SMTRAT_LOG_DEBUG("smtrat", "Positive Candidates:");
1445  for(auto h : positive_candidates){
1446  SMTRAT_LOG_DEBUG("smtrat", h.first << " " << h.second->expression());
1447  }
1448  SMTRAT_LOG_DEBUG("smtrat", "Negative Candidates:");
1449  for(auto h : negative_candidates){
1450  SMTRAT_LOG_DEBUG("smtrat", h.first << " " << h.second->expression());
1451  }
1452  #endif
1453 
1454  std::map<Value<T1>, Value<T1>> leaving_dVio;
1455  // according to the paper is suffices to consider only one direction
1456  if(mInfeasibilityRow[column] < 0){
1457  leaving_dVio = compute_dVio(positive_candidates, nVar, true);
1458  }
1459  else{
1460  leaving_dVio = compute_dVio(negative_candidates, nVar, false);
1461  }
1462 
1463  #ifdef DEBUG_SOI_SIMPLEX
1464  std::stringstream s;
1465  for(auto e : leaving_candidates){
1466  s << "(" << e.first << "," << e.second->position() << ")";
1467  }
1468  SMTRAT_LOG_DEBUG("smtrat", "Candidates for "<< nVar.position() << "(" << nVar.expression() << ")" << ": "<<s.str());
1469  s.clear();
1470  s.str("");
1471  #endif
1472 
1473  if(leaving_candidates.empty())
1474  continue;
1475 
1476  // select (delta, k) pair to minimize dVio (line 22)
1477  assert(!leaving_candidates.empty());
1478 
1479  std::pair<Value<T1>,Variable<T1,T2>*>* min_pair;
1480  Value<T1> min_val = Value<T1>(std::numeric_limits<double>::max());
1481  SMTRAT_LOG_DEBUG("smtrat", "Size leaving candidates is "<< leaving_candidates.size());
1482  for(std::size_t candidate_index = 0; candidate_index < leaving_candidates.size(); candidate_index++){ // can be optimized to only consider positive/negative values
1483  std::pair< Value<T1>, Variable<T1,T2>* >* candidate = &leaving_candidates[candidate_index]; // rename for easier access
1484  if(leaving_dVio.find(candidate->first) == leaving_dVio.end()){
1485  // doesnt need to be considered
1486  continue;
1487  }
1488  assert(leaving_dVio.find(candidate->first) != leaving_dVio.end());
1489  Value<T1> helper_diff = leaving_dVio[candidate->first];;
1490  SMTRAT_LOG_DEBUG("smtrat", "candidate is" << "(" << candidate->first << "," << candidate->second->position() << ")" << " id " << candidate->second->id());
1491 
1492  // tie-breaking in the size of variable numbers
1493  if( helper_diff < min_val){
1494  min_val = helper_diff;
1495  min_pair = candidate;
1496  SMTRAT_LOG_DEBUG("smtrat","Min_pair set to: "<< "(" << min_pair->first << "," << min_pair->second->position() << ")" );
1497  }
1498  else if(helper_diff == min_val && candidate->first.abs() > min_pair->first.abs() ){
1499  min_val = helper_diff;
1500  min_pair = candidate;
1501  SMTRAT_LOG_DEBUG("smtrat", "Min_pair set to: "<< "(" << min_pair->first << "," << min_pair->second->position() << ")" );
1502  }
1503  else if((helper_diff == min_val && candidate->first.abs() == min_pair->first.abs() && min_pair->second->id() > candidate->second->id() )){
1504  SMTRAT_LOG_DEBUG("smtrat", "ID: " << min_pair->second->id() << " " << candidate->second->id() );
1505  min_val = helper_diff;
1506  min_pair = candidate;
1507  SMTRAT_LOG_DEBUG("smtrat", "Min_pair set to: "<< "(" << min_pair->first << "," << min_pair->second->position() << ")" );
1508  }
1509 
1510  SMTRAT_LOG_DEBUG("smtrat", "dVio updated" << min_val);
1511  }
1512 
1513  SMTRAT_LOG_DEBUG("smtrat", "min dVio " << min_val);
1514  assert(min_pair != NULL);
1515  update_candidates.emplace_back(column, min_pair->first, min_pair->second);
1516  SMTRAT_LOG_DEBUG("smtrat", "Added touple ("<< column <<"," << min_pair->first <<"," << min_pair->second->position() << ") to update_candidates");
1517 
1518  if(min_val < Value<T1>(0.0)){
1519  SMTRAT_LOG_DEBUG("smtrat","Found Element with dVio < 0, continue now");
1520  found_improvement = true;
1521  }
1522 
1523  // Abort computation of candidates after finding one with progress (atleast mFullCandidateSearch times all canidates are evaluated, then the amount increases depending on state of computation)
1524  if(mPivotingSteps > mFullCandidateSearch && (double)column > (double)mPivotingSteps/5.0){
1525  if(found_improvement){
1526  break;
1527  }
1528  }
1529  }
1530 
1531  if(update_candidates.empty()){// represents case that there is no flexible variable
1532  mOldVioSum = Value<T1>(-1);
1533  mOld_dVioSum = Value<T1>(1);
1534  std::vector<Variable<T1, T2>*> errorVars;
1535  for(std::size_t row = 0; row < mRows.size(); row++){
1536  assert( mRows[row] != NULL ); // non_basic variable must exist
1537 
1538  const Variable<T1,T2>& bVar = *mRows[row];
1539  assert(bVar.isBasic());
1540 
1541  if(bVar.infimum() > bVar.assignment() || bVar.supremum() < bVar.assignment()){
1542  errorVars.push_back(mRows[row]);
1543  }
1544  }
1545 
1546  if(errorVars.empty()){
1547  // sat case
1548  SMTRAT_LOG_DEBUG("smtrat", "In SAT case.");
1549  return std::pair<EntryID,bool>( lra::LAST_ENTRY_ID, true );
1550  }
1551  else{
1552  // Multiline conflict detection case (conflict!= 0 and E != 0 )
1553  SMTRAT_LOG_DEBUG("smtrat", "Multiline conflict case, size: " << errorVars.size());
1554  //assert(false);
1555  return std::pair<EntryID,bool>( lra::LAST_ENTRY_ID, false );
1556  }
1557  }
1558 
1559  assert(!update_candidates.empty());
1560 
1561  //entering rule: select minimizing update rule
1562  // line 24 in algorithm
1563  #ifdef DEBUG_SOI_SIMPLEX
1564  SMTRAT_LOG_DEBUG("smtrat", "Size update_candidates: " << update_candidates.size());
1565  for(auto e : update_candidates){
1566  s << "(" << std::get<0>(e) << "," << std::get<1>(e) << "," << std::get<2>(e)->position() << ")";
1567  s << ":" << mInfeasibilityRow[std::get<0>(e)] << ",";
1568  }
1569  SMTRAT_LOG_DEBUG("smtrat", "Update candidates: " <<s.str());
1570  s.clear();
1571  s.str("");
1572  #endif
1573 
1574  std::tuple <std::size_t, Value<T1>, Variable<T1,T2>* >* min_pair;
1575  Value<T1> min_val = Value<T1>(std::numeric_limits<double>::max());
1576  for(std::size_t candidate_index = 0; candidate_index < update_candidates.size(); candidate_index++){
1577  std::tuple<std::size_t, Value<T1>, Variable<T1,T2>*>* candidate = &update_candidates[candidate_index];
1578  // compute sign(dVio)* abs(mInfeasibilityRow[v_j])
1579  const Variable<T1,T2>& nVar = *mColumns[std::get<0>(*candidate)];
1580  Value<T1> vio_sum = dViolationSum(&nVar, std::get<1>(*candidate));
1581  Value<T1> weighted_dVio;
1582  if(vio_sum == Value<T1>(0)){
1583  weighted_dVio = Value<T1>(0);
1584  }
1585  else{
1586  assert(vio_sum != Value<T1>(0));
1587  weighted_dVio = vio_sum < Value<T1>(0) ? Value<T1>(-1) : Value<T1>(1); // non-zero part of sign function
1588  }
1589  weighted_dVio *= abs(mInfeasibilityRow[std::get<0>(*candidate)]);
1590  if(min_val > weighted_dVio) {
1591  min_val = weighted_dVio;
1592  min_pair = candidate;
1593  }
1594  else if(min_val == weighted_dVio && std::get<1>(*candidate) < std::get<1>(*min_pair)){ // tie-break in second argument
1595  min_val = weighted_dVio;
1596  min_pair = candidate;
1597  }
1598  }
1599  assert(min_pair != NULL);
1600  assert(min_val <= Value<T1>(0));
1601 
1602  // set update value
1603  *mpTheta = std::get<1>(*min_pair);
1604 
1605  // min_pair: pivot 1 with 3 (1 is nB and 3 is basic)
1606  // if 3 is nB, must be same as 1
1607  // 1 == 3 -> only update operation (no pivot)
1608  EntryID returnPosition = LAST_ENTRY_ID;
1609 
1610  SMTRAT_LOG_DEBUG("smtrat", "Weighted min dVio" << min_val);
1611  mOld_dVioSum = dViolationSum(mColumns[std::get<0>(*min_pair)], std::get<1>(*min_pair));
1612 
1613  if(!std::get<2>(*min_pair)->isBasic()){
1614  // Only update case
1615  assert(std::get<2>(*min_pair)->position() == std::get<0>(*min_pair));
1616  SMTRAT_LOG_DEBUG("smtrat", "Apply update without pivot on var " << std::get<0>(*min_pair) << " by "<< std::get<1>(*min_pair));
1617 
1618  Iterator columnIter = Iterator( mColumns[std::get<0>(*min_pair)]->startEntry(), mpEntries );
1619  returnPosition = columnIter.entryID();
1620  onlyUpdate = true;
1621  }
1622  else{
1623  SMTRAT_LOG_DEBUG("smtrat", "Pivoting pairs of " << *mpTheta << " found, at candidate ("<< std::get<0>(*min_pair) << "," << std::get<1>(*min_pair)<< "," << std::get<2>(*min_pair)->position() << ")");
1624 
1625  assert(std::get<2>(*min_pair)->isBasic());
1626  // apply update with min_pair
1627  // compute pivot entry (via nonbasic variable)
1628  Iterator columnIter = Iterator( mColumns[std::get<0>(*min_pair)]->startEntry(), mpEntries );
1629  while(true){
1630  const Variable<T1, T2>* basicVar = (*columnIter).rowVar();
1631  SMTRAT_LOG_DEBUG("smtrat", "Basic var" << basicVar->position() << "min pair" << std::get<2>(*min_pair)->position());
1632  if(basicVar == std::get<2>(*min_pair)){
1633  returnPosition = columnIter.entryID();
1634  break;
1635  }
1636  if( columnIter.vEnd( false ) ) {
1637  assert(false); // make sure that pivot element is always found
1638  }
1639  else{
1640  columnIter.vMove( false );
1641  }
1642  }
1643  assert(returnPosition != LAST_ENTRY_ID);
1644  }
1645 
1646  // entry: position in tableau to be pivoted, bool: true if NO conflict found
1647  return std::pair<EntryID,bool>( returnPosition, true );
1648  }
1649  else // Bland's rule
1650  {
1651  const Variable<T1, T2>* bestBasicVar = nullptr;
1652  std::pair<EntryID,bool> bestResult( LAST_ENTRY_ID, true );
1653  for( const Variable<T1, T2>* basicVar : mRows )
1654  {
1655  assert( basicVar != NULL );
1656  const Variable<T1,T2>& bVar = *basicVar;
1657  Value<T1> thetaB = Value<T1>( 0 );
1658  bool upperBoundViolated = false;
1659  bool lowerBoundViolated = false;
1660  if( bVar.supremum() < bVar.assignment() )
1661  {
1662  thetaB = bVar.supremum().limit() - bVar.assignment();
1663  upperBoundViolated = true;
1664  }
1665  else if( bVar.infimum() > bVar.assignment() )
1666  {
1667  thetaB = bVar.infimum().limit() - bVar.assignment();
1668  lowerBoundViolated = true;
1669  }
1670  if( upperBoundViolated || lowerBoundViolated )
1671  {
1672  EntryID result = isSuitable( bVar, lowerBoundViolated );
1673  if( result == LAST_ENTRY_ID )
1674  {
1675  // Found a conflicting row.
1676  return std::pair<EntryID,bool>( bVar.startEntry(), false );
1677  }
1678  else
1679  {
1680  if( bestBasicVar == nullptr || *basicVar < *bestBasicVar )
1681  {
1682  bestBasicVar = basicVar;
1683  // Found a pivoting element
1684  *mpTheta = thetaB;
1685  if( Settings::omit_division )
1686  (*mpTheta) *= bVar.factor();
1687  (*mpTheta) /= (*mpEntries)[result].content();
1688  bestResult = std::pair<EntryID,bool>( result, true );
1689  }
1690  }
1691  }
1692  }
1693  return bestResult;
1694  }
1695  SMTRAT_LOG_FATAL("smtrat", "No return statement called!");
1696  assert(false);
1697  }
1698 
1699  template<class Settings, typename T1, typename T2>
1700  bool Tableau<Settings,T1,T2>::hasMultilineConflict(){
1701  for(std::size_t column = 0; column < mColumns.size(); column++){
1702  assert( mColumns[column] != NULL ); // non_basic variable must exist
1703 
1704  const Variable<T1,T2>& nVar = *mColumns[column];
1705  assert(!nVar.isBasic());
1706  // check if nonbasis variable is flexible
1707  bool isFlexible = false;
1708  if(mInfeasibilityRow[column] > 0 && nVar.infimum() < nVar.assignment()){
1709  isFlexible = true;
1710  }
1711  else if (mInfeasibilityRow[column] < 0 && nVar.supremum() > nVar.assignment()){
1712  isFlexible = true;
1713  }
1714 
1715  if(isFlexible){
1716  return false;
1717  }
1718  }
1719 
1720  // rebuild set of error variables
1721  std::vector<Variable<T1, T2>*> errorVars;
1722  for(std::size_t row = 0; row < mRows.size(); row++){
1723  assert( mRows[row] != NULL ); // non_basic variable must exist
1724 
1725  const Variable<T1,T2>& bVar = *mRows[row];
1726  assert(bVar.isBasic());
1727 
1728  if(bVar.infimum() > bVar.assignment() || bVar.supremum() < bVar.assignment()){
1729  errorVars.push_back(mRows[row]);
1730  }
1731  }
1732 
1733  if(errorVars.empty()){
1734  return false;
1735  }
1736 
1737  return true;
1738  }
1739 
1740  /**
1741  *
1742  * @return In SUM-Simplex a conflict is possibly created with multiple rows.
1743  * This function returns the single conflict of possible and otherwise construct the multiline conflict.
1744  */
1745  template<class Settings, typename T1, typename T2>
1746  std::vector< const Bound<T1, T2>* > Tableau<Settings,T1,T2>::getMultilineConflict(){
1747  SMTRAT_LOG_DEBUG("smtrat", "Searching one line conflict");
1748  std::pair<EntryID,bool> conflictPair = hasConflict();
1749  if(!conflictPair.second){
1750  SMTRAT_LOG_DEBUG("smtrat", "Simple conflict found in multiline conflict search.");
1751  return getConflict(conflictPair.first);
1752  }
1753  SMTRAT_LOG_DEBUG("smtrat", "Constructing multiline conflict.");
1754  std::vector< const Bound<T1, T2>* > conflict;
1755 
1756  std::vector<T2> infeas_row = getInfeasibilityRow();
1757  //add lower or upper bounds depending on infeasibility row
1758  for(std::size_t nonbasic_index = 0; nonbasic_index<infeas_row.size(); nonbasic_index++){
1759  if(infeas_row[nonbasic_index] > 0 && !mColumns[nonbasic_index]->pInfimum()->isInfinite()){
1760  conflict.push_back(mColumns[nonbasic_index]->pInfimum());
1761  }
1762  if(infeas_row[nonbasic_index] < 0 && !mColumns[nonbasic_index]->pSupremum()->isInfinite()){
1763  conflict.push_back(mColumns[nonbasic_index]->pSupremum());
1764  }
1765  }
1766 
1767  // rebuild set of error variables
1768  std::vector<Variable<T1, T2>*> errorVars;
1769  for(std::size_t row = 0; row < mRows.size(); row++){
1770  assert( mRows[row] != NULL ); // non_basic variable must exist
1771 
1772  const Variable<T1,T2>& bVar = *mRows[row];
1773  assert(bVar.isBasic());
1774 
1775  if(bVar.infimum() > bVar.assignment() || bVar.supremum() < bVar.assignment()){
1776  errorVars.push_back(mRows[row]);
1777  }
1778  }
1779 
1780  for(Variable<T1, T2>* v: errorVars){
1781  int direction;
1782  if(v->supremum() < v->assignment()){
1783  direction = 1;
1784  }
1785  else if(v->infimum() > v->assignment()){
1786  direction = -1;
1787  }
1788  else{
1789  direction = 0;
1790  }
1791  assert(direction != 0);
1792 
1793  if( v->infimum() > v->assignment()){
1794  // lower bound is broken
1795  conflict.push_back(v->pInfimum());
1796  }
1797  else{
1798  conflict.push_back(v->pSupremum());
1799  }
1800  }
1801  SMTRAT_LOG_DEBUG("smtrat", "Found multiline-conflict:");
1802  for([[maybe_unused]] auto c: conflict){
1803  SMTRAT_LOG_DEBUG("smtrat", *c->origins().begin());
1804  }
1805  return conflict;
1806  }
1807 
1808  template<class Settings, typename T1, typename T2>
1809  std::pair<EntryID,bool> Tableau<Settings,T1,T2>::nextPivotingElement()
1810  {
1811  // Dynamic strategy for a fixed number of steps
1812  if( Settings::use_pivoting_strategy && mPivotingSteps < mMaxPivotsWithoutBlandsRule )
1813  {
1814  FindPivot:
1815  EntryID bestTableauEntry = LAST_ENTRY_ID;
1816  EntryID beginOfFirstConflictRow = LAST_ENTRY_ID;
1817  Value<T1> bestDiff = Value<T1>( 0 );
1818  Value<T1> bestThetaB = Value<T1>( 0 );
1819  bool initialSearch = mConflictingRows.empty();
1820  std::vector<Variable<T1,T2>*>& rowsToConsider = Settings::pivot_into_local_conflict && initialSearch ? mRows : mConflictingRows;
1821  // TODO: instead of running through all rows, just go through those which got conflicting
1822  typename std::vector<Variable<T1,T2>*>::iterator bestVar = rowsToConsider.end();
1823  for( auto basicVar = rowsToConsider.begin(); basicVar != rowsToConsider.end(); )
1824  {
1825  assert( *basicVar != NULL );
1826  Variable<T1,T2>& bVar = **basicVar;
1827  Value<T1> diff = Value<T1>( 0 );
1828  Value<T1> thetaB = Value<T1>( 0 );
1829  bool upperBoundViolated = false;
1830  bool lowerBoundViolated = false;
1831  if( bVar.supremum() < bVar.assignment() )
1832  {
1833  thetaB = bVar.supremum().limit() - bVar.assignment();
1834  diff = thetaB * T2(-1);
1835  upperBoundViolated = true;
1836  }
1837  else if( bVar.infimum() > bVar.assignment() )
1838  {
1839  thetaB = bVar.infimum().limit() - bVar.assignment();
1840  diff = thetaB;
1841  lowerBoundViolated = true;
1842  }
1843  else
1844  {
1845  if( Settings::pivot_into_local_conflict && !initialSearch )
1846  {
1847  bool resetBestVarToEnd = bestVar == mConflictingRows.end();
1848  basicVar = mConflictingRows.erase( basicVar );
1849  if( resetBestVarToEnd ) bestVar = mConflictingRows.end();
1850  if( mConflictingRows.empty() )
1851  {
1852  goto FindPivot;
1853  }
1854  }
1855  else
1856  {
1857  ++basicVar;
1858  }
1859  continue;
1860  }
1861  if( Settings::use_theta_based_pivot_strategy )
1862  {
1863  if( diff <= bestDiff )
1864  {
1865  ++basicVar;
1866  continue;
1867  }
1868  }
1869  else if( Settings::use_activity_based_pivot_strategy && bestVar != rowsToConsider.end() )
1870  {
1871  if( (*basicVar)->conflictActivity() < (*bestVar)->conflictActivity()
1872  || ((*basicVar)->conflictActivity() == (*bestVar)->conflictActivity() && diff <= bestDiff) )
1873  {
1874  ++basicVar;
1875  continue;
1876  }
1877  }
1878  if( upperBoundViolated || lowerBoundViolated )
1879  {
1880  EntryID result = isSuitable( bVar, lowerBoundViolated );
1881  if( result == LAST_ENTRY_ID )
1882  {
1883  bestTableauEntry = LAST_ENTRY_ID;
1884  // Found a conflicting row.
1885  if( beginOfFirstConflictRow == LAST_ENTRY_ID )
1886  {
1887  beginOfFirstConflictRow = bVar.startEntry();
1888  bestVar = basicVar;
1889  break;
1890  }
1891  }
1892  else
1893  {
1894  if( bestVar == rowsToConsider.end() )
1895  {
1896  bestTableauEntry = result;
1897  bestVar = basicVar;
1898  bestDiff = diff;
1899  bestThetaB = thetaB;
1900  }
1901  else
1902  {
1903  assert( result != LAST_ENTRY_ID );
1904  assert( bestVar != rowsToConsider.end() );
1905  assert( bestTableauEntry != LAST_ENTRY_ID );
1906  if( Settings::prefer_equations )
1907  {
1908  if( !(*bestVar)->involvesEquation() && bVar.involvesEquation() )
1909  {
1910  bestTableauEntry = result;
1911  bestVar = basicVar;
1912  }
1913  else if( (*bestVar)->involvesEquation() || !bVar.involvesEquation() )
1914  {
1915  bestTableauEntry = result;
1916  bestThetaB = thetaB;
1917  bestDiff = diff;
1918  if( Settings::pivot_into_local_conflict && initialSearch && (*bestVar)->involvesEquation() )
1919  mConflictingRows.push_back( *bestVar );
1920  bestVar = basicVar;
1921  }
1922  }
1923  else
1924  {
1925 
1926  bestTableauEntry = result;
1927  bestThetaB = thetaB;
1928  bestDiff = diff;
1929  bestVar = basicVar;
1930  }
1931  }
1932  }
1933  }
1934  ++basicVar;
1935  }
1936  if( bestTableauEntry == LAST_ENTRY_ID && beginOfFirstConflictRow != LAST_ENTRY_ID )
1937  {
1938  // Found a conflict
1939  if( Settings::pivot_into_local_conflict )
1940  mConflictingRows.clear();
1941  return std::pair<EntryID,bool>( beginOfFirstConflictRow, false );
1942  }
1943  else if( bestTableauEntry != LAST_ENTRY_ID )
1944  {
1945  // The best pivoting element found
1946  *mpTheta = bestThetaB;
1947  if( Settings::omit_division )
1948  (*mpTheta) *= (*bestVar)->factor();
1949  (*mpTheta) /= (*mpEntries)[bestTableauEntry].content();
1950  if( Settings::pivot_into_local_conflict && !initialSearch )
1951  {
1952  assert( bestVar != mConflictingRows.end() );
1953  mConflictingRows.erase( bestVar );
1954  }
1955  return std::pair<EntryID,bool>( bestTableauEntry, true );
1956  }
1957  else
1958  {
1959  // Found no pivoting element, that is no variable violates its bounds.
1960  assert( mConflictingRows.empty() );
1961  return std::pair<EntryID,bool>( LAST_ENTRY_ID, true );
1962  }
1963  }
1964  else // Bland's rule
1965  {
1966  const Variable<T1, T2>* bestBasicVar = nullptr;
1967  std::pair<EntryID,bool> bestResult( LAST_ENTRY_ID, true );
1968  for( const Variable<T1, T2>* basicVar : mRows )
1969  {
1970  assert( basicVar != NULL );
1971  const Variable<T1,T2>& bVar = *basicVar;
1972  Value<T1> thetaB = Value<T1>( 0 );
1973  bool upperBoundViolated = false;
1974  bool lowerBoundViolated = false;
1975  if( bVar.supremum() < bVar.assignment() )
1976  {
1977  thetaB = bVar.supremum().limit() - bVar.assignment();
1978  upperBoundViolated = true;
1979  }
1980  else if( bVar.infimum() > bVar.assignment() )
1981  {
1982  thetaB = bVar.infimum().limit() - bVar.assignment();
1983  lowerBoundViolated = true;
1984  }
1985  if( upperBoundViolated || lowerBoundViolated )
1986  {
1987  EntryID result = isSuitable( bVar, lowerBoundViolated );
1988  if( result == LAST_ENTRY_ID )
1989  {
1990  // Found a conflicting row.
1991  return std::pair<EntryID,bool>( bVar.startEntry(), false );
1992  }
1993  else
1994  {
1995  if( bestBasicVar == nullptr || *basicVar < *bestBasicVar )
1996  {
1997  bestBasicVar = basicVar;
1998  // Found a pivoting element
1999  *mpTheta = thetaB;
2000  if( Settings::omit_division )
2001  (*mpTheta) *= bVar.factor();
2002  (*mpTheta) /= (*mpEntries)[result].content();
2003  bestResult = std::pair<EntryID,bool>( result, true );
2004  }
2005  }
2006  }
2007  }
2008  return bestResult;
2009  }
2010  }
2011 
2012  template<class Settings, typename T1, typename T2>
2013  std::pair<EntryID,bool> Tableau<Settings,T1,T2>::optimizeIndependentNonbasics( const Variable<T1, T2>& _objective )
2014  {
2015  EntryID firstColumnToCheck = LAST_ENTRY_ID;
2016  Iterator objectiveIter = Iterator( _objective.startEntry(), mpEntries );
2017  while( true )
2018  {
2019  Variable<T1, T2>& varForMinimizaton = *(*objectiveIter).columnVar();
2020  if( (*mpEntries)[varForMinimizaton.startEntry()].vNext( false ) == LAST_ENTRY_ID ) // Non-basic variable only occurs in objective function
2021  {
2022  bool increaseVar = Settings::omit_division ?
2023  ((*objectiveIter).content() < 0 && _objective.factor() > 0) || ((*objectiveIter).content() > 0 && _objective.factor() < 0) :
2024  (*objectiveIter).content() < 0;
2025  if( increaseVar )
2026  {
2027  if( varForMinimizaton.supremum().isInfinite() )
2028  {
2029  return std::make_pair( LAST_ENTRY_ID, true );
2030  }
2031  else if( varForMinimizaton.supremum() > varForMinimizaton.assignment() )
2032  {
2033  updateBasicAssignments( varForMinimizaton.position(), varForMinimizaton.supremum().limit() - varForMinimizaton.assignment() );
2034  varForMinimizaton.rAssignment() = varForMinimizaton.supremum().limit();
2035  }
2036  }
2037  else
2038  {
2039  if( varForMinimizaton.infimum().isInfinite() )
2040  {
2041  return std::make_pair( LAST_ENTRY_ID, true );
2042  }
2043  else if( varForMinimizaton.infimum() < varForMinimizaton.assignment() )
2044  {
2045  updateBasicAssignments( varForMinimizaton.position(), varForMinimizaton.assignment() - varForMinimizaton.infimum().limit() );
2046  varForMinimizaton.rAssignment() = varForMinimizaton.infimum().limit();
2047  }
2048  }
2049  }
2050  else if( firstColumnToCheck == LAST_ENTRY_ID )
2051  {
2052  firstColumnToCheck = objectiveIter.entryID();
2053  }
2054  if( objectiveIter.hEnd( false ) )
2055  break;
2056  else
2057  objectiveIter.hMove( false );
2058  }
2059  return std::make_pair( firstColumnToCheck, firstColumnToCheck != LAST_ENTRY_ID );
2060  }
2061 
2062  template<class Settings, typename T1, typename T2>
2063  std::pair<EntryID,bool> Tableau<Settings,T1,T2>::nextPivotingElementForOptimizing( const Variable<T1, T2>& _objective )
2064  {
2065  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2066  std::cout << "Find next pivoting element to minimize ";
2067  _objective.print();
2068  std::cout << " in:" << std::endl;
2069  print();
2070  #endif
2071  assert( _objective.isBasic() );
2072  assert( _objective.positionInNonActives() == mNonActiveBasics.end() );
2073  std::pair<EntryID,bool> ret = optimizeIndependentNonbasics( _objective );
2074  if( ret.first == LAST_ENTRY_ID )
2075  return ret;
2076  Iterator objectiveIter = Iterator( ret.first, mpEntries );
2077  assert( _objective.infimum().isInfinite() );
2078  // Default value for infinity.
2079  Value<T1> infinityValue = Value<T1>(T1(-1));
2080  // The best entry to pivot at to achieve the best improvement (bestImprovement) on the objective function.
2081  EntryID bestPivotingEntry = LAST_ENTRY_ID;
2082  Value<T1> bestImprovement = Value<T1>(T1(0));
2083  // Go through all columns of the objective functions row.
2084  while( true )
2085  {
2086  const Variable<T1, T2>& columnVar = *(*objectiveIter).columnVar();
2087  bool increaseColumnVar = Settings::omit_division ?
2088  ((*objectiveIter).content() < 0 && _objective.factor() > 0) || ((*objectiveIter).content() > 0 && _objective.factor() < 0) :
2089  (*objectiveIter).content() < 0;
2090  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2091  std::cout << "Check the non-basic variable "; columnVar.print(); std::cout << std::endl;
2092  std::cout << " " << (increaseColumnVar ? "Increase" : "Decrease") << " non-basic variable's assignment." << std::endl;
2093  #endif
2094  // The margin of the column variable according to its bounds (-1 if it is infinite)
2095  Value<T1> columnVarMargin = increaseColumnVar ?
2096  (columnVar.supremum().isInfinite() ? infinityValue : (columnVar.supremum().limit() - columnVar.assignment())) :
2097  (columnVar.infimum().isInfinite() ? infinityValue : (columnVar.assignment() - columnVar.infimum().limit()));
2098  if( !columnVarMargin.is_zero() )
2099  {
2100  // Calculate the change we minimally need on the column variable in order to improve the objective
2101  // more than with the currently best found pivoting entry.
2102  Value<T1> minNeededColumnVarChange = bestImprovement;
2103  assert( bestImprovement >= T1(0) );
2104  if( Settings::omit_division )
2105  minNeededColumnVarChange *= _objective.factor();
2106  minNeededColumnVarChange /= (*objectiveIter).content();
2107  minNeededColumnVarChange.abs_();
2108  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2109  std::cout << " We need more change on non-basic variable than: " << minNeededColumnVarChange << std::endl;
2110  #endif
2111  // This column variable allows more improvement than we could gain so far.
2112  if( columnVarMargin == infinityValue || columnVarMargin > minNeededColumnVarChange )
2113  {
2114  // Search the value, for which we change this column variable without violating any bound of the row variables in it's column.
2115  Value<T1> minColumnVarChange = infinityValue;
2116  EntryID criticalColumnEntry = LAST_ENTRY_ID;
2117  Iterator columnIter = Iterator( columnVar.startEntry(), mpEntries );
2118  assert( columnIter.vEnd( true ) ); // Is the lowest row.
2119  assert( !columnIter.vEnd( false ) ); // There should be at least one more row containing this column variable.
2120  // Skip the objective function's row.
2121  columnIter.vMove( false );
2122  while( true )
2123  {
2124  Variable<T1, T2>& rowVar = *((*columnIter).rowVar());
2125  assert( rowVar != _objective );
2126  bool entryNegative = Settings::omit_division ?
2127  ((*columnIter).content() < 0 && rowVar.factor() > 0) || ((*columnIter).content() > 0 && rowVar.factor() < 0) :
2128  (*columnIter).content() < 0;
2129  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2130  std::cout << " Check the basic variable "; rowVar.print(); std::cout << std::endl;
2131  std::cout << " " << (increaseColumnVar != entryNegative ? "Increase" : "Decrease") << " basic variable's assignment." << std::endl;
2132  #endif
2133  Value<T1> changeOnColumnVar = (increaseColumnVar == entryNegative) ?
2134  (rowVar.infimum().isInfinite() ? infinityValue : rowVar.assignment() - rowVar.infimum().limit()) :
2135  (rowVar.supremum().isInfinite() ? infinityValue : rowVar.supremum().limit() - rowVar.assignment());
2136  if( changeOnColumnVar == infinityValue )
2137  {
2138  // If this is the first entry to be checked, take it as currently most critical entry (which has actually no constraint).
2139  if( criticalColumnEntry == LAST_ENTRY_ID )
2140  {
2141  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2142  std::cout << " Take basic variable allowing infinite change." << std::endl;
2143  #endif
2144  criticalColumnEntry = columnIter.entryID();
2145  }
2146  }
2147  else if( changeOnColumnVar == T1(0) )
2148  {
2149  // No change allowed making this column variable not suitable at all.
2150  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2151  std::cout << " The non-basic variable has no margin." << std::endl;
2152  #endif
2153  break;
2154  }
2155  else
2156  {
2157  if( Settings::omit_division )
2158  changeOnColumnVar *= rowVar.factor();
2159  changeOnColumnVar /= (*columnIter).content();
2160  changeOnColumnVar.abs_();
2161  if( columnVarMargin != infinityValue && changeOnColumnVar > columnVarMargin )
2162  changeOnColumnVar = columnVarMargin;
2163  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2164  std::cout << " Possible change on non-basic variable: " << changeOnColumnVar << std::endl;
2165  #endif
2166  if( changeOnColumnVar > minNeededColumnVarChange )
2167  {
2168  if( minColumnVarChange == infinityValue || changeOnColumnVar < minColumnVarChange )
2169  {
2170  // Found a row which allows less change on the column variable.
2171  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2172  std::cout << " Take basic variable with stricter margin." << std::endl;
2173  #endif
2174  minColumnVarChange = changeOnColumnVar;
2175  criticalColumnEntry = columnIter.entryID();
2176  }
2177  }
2178  else
2179  {
2180  // This column cannot improve the best improvement found yet .. go to the next column.
2181  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2182  std::cout << " The non-basic variable cannot improve the best found one." << std::endl;
2183  #endif
2184  break;
2185  }
2186  }
2187  if( columnIter.vEnd( false ) )
2188  {
2189  // All rows are inspected.
2190  assert( criticalColumnEntry != LAST_ENTRY_ID );
2191  if( minColumnVarChange == infinityValue )
2192  {
2193  // No row variable constraints the column variable in the desired direction.
2194  if( columnVarMargin == infinityValue )
2195  {
2196  // Infinite change possible, which cannot be improved.
2197  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2198  std::cout << " Take this non-basic variable having no limits." << std::endl;
2199  #endif
2200  return std::make_pair( LAST_ENTRY_ID, true );
2201  }
2202  else
2203  {
2204  // We can use the whole margin of the column variable.
2205  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2206  std::cout << " Use the whole margin of the non-basic variable." << std::endl;
2207  #endif
2208  minColumnVarChange = columnVarMargin;
2209  }
2210  }
2211  assert( minColumnVarChange > minNeededColumnVarChange );
2212  // Calculate improvement on objective.
2213  minColumnVarChange *= (*objectiveIter).content();
2214  if( Settings::omit_division )
2215  minColumnVarChange /= _objective.factor();
2216  minColumnVarChange.abs_();
2217  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2218  std::cout << " Take this non-basic variable improving the objective by " << minColumnVarChange << std::endl;
2219  #endif
2220  // Set it as the new best improvement on objective.
2221  assert( minColumnVarChange > bestImprovement );
2222  bestImprovement = minColumnVarChange;
2223  // Take the corresponding entry as the best for pivoting.
2224  bestPivotingEntry = criticalColumnEntry;
2225  break;
2226  }
2227  else
2228  {
2229  columnIter.vMove( false );
2230  }
2231  }
2232  }
2233  }
2234  if( objectiveIter.hEnd( false ) )
2235  break;
2236  else
2237  {
2238  objectiveIter.hMove( false );
2239  }
2240  }
2241  if( bestPivotingEntry == LAST_ENTRY_ID )
2242  {
2243  // We could not find any suitable column variable for pivoting -> Minimum reached.
2244  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2245  std::cout << "No non-basic variable is suitable to actually gain something." << std::endl;
2246  #endif
2247  return nextZeroPivotingElementForOptimizing( _objective );
2248  }
2249  assert( bestImprovement != infinityValue );
2250  assert( bestImprovement > T1(0) );
2251  // Calculate theta (how much do we change the assignment of the column/non-basic variable when pivoting).
2252  *mpTheta = bestImprovement * T1(-1); // Change on objective (negative as we are minimizing).
2253  const Variable<T1, T2>& bestColumnVar = *((*mpEntries)[bestPivotingEntry].columnVar());
2254  Iterator columnIter = Iterator( bestColumnVar.startEntry(), mpEntries );
2255  assert( *(*columnIter).rowVar() == _objective );
2256  if( Settings::omit_division )
2257  *mpTheta *= _objective.factor();
2258  *mpTheta /= (*columnIter).content();
2259  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2260  std::cout << "Found the non-basic variable " << (*mpEntries)[bestPivotingEntry].columnVar()->expression() << " optimizing the objective by " << bestImprovement;
2261  std::cout << " if we pivot it with " << (*mpEntries)[bestPivotingEntry].rowVar()->expression();
2262  std::cout << " where theta (change on non-basic variable) is " << *mpTheta << "." << std::endl;
2263  #endif
2264  return std::make_pair( bestPivotingEntry, true );
2265  }
2266 
2267  template<class Settings, typename T1, typename T2>
2268  std::pair<EntryID,bool> Tableau<Settings,T1,T2>::nextZeroPivotingElementForOptimizing( const Variable<T1, T2>& _objective ) const
2269  {
2270  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2271  std::cout << "Try to find a pivoting candidate, such that the non-basic (column) variable might improve the objective later." << std::endl;
2272  #endif
2273  EntryID smallestPivotingElement = LAST_ENTRY_ID;
2274  // find the smallest non-basic variable, which allows us to gain optimization
2275  const Variable<T1, T2>* bestNonBasicVar = nullptr;
2276  bool increaseBestNonbasicVar = true;
2277  Iterator objectiveIter = Iterator( _objective.startEntry(), mpEntries );
2278  while( true )
2279  {
2280  const Variable<T1, T2>* varForMinimizaton = (*objectiveIter).columnVar();
2281  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2282  std::cout << "Checking " << std::endl; varForMinimizaton->print(); std::cout << " .. " << std::endl;
2283  #endif
2284  if( (*mpEntries)[varForMinimizaton->startEntry()].vNext( false ) != LAST_ENTRY_ID ) // Non-basic variable only occurs in objective function
2285  {
2286  if( bestNonBasicVar == nullptr || *varForMinimizaton < *bestNonBasicVar )
2287  {
2288  bool increaseVar = Settings::omit_division ?
2289  ((*objectiveIter).content() < 0 && _objective.factor() > 0) || ((*objectiveIter).content() > 0 && _objective.factor() < 0) :
2290  (*objectiveIter).content() < 0;
2291  if( increaseVar )
2292  {
2293  if( varForMinimizaton->supremum().isInfinite() || varForMinimizaton->supremum() > varForMinimizaton->assignment() )
2294  {
2295  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2296  std::cout << "Is smaller!" << std::endl;
2297  #endif
2298  bestNonBasicVar = varForMinimizaton;
2299  increaseBestNonbasicVar = increaseVar;
2300  }
2301  }
2302  else
2303  {
2304  if( varForMinimizaton->infimum().isInfinite() || varForMinimizaton->infimum() < varForMinimizaton->assignment() )
2305  {
2306  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2307  std::cout << "Is smaller!" << std::endl;
2308  #endif
2309  bestNonBasicVar = varForMinimizaton;
2310  increaseBestNonbasicVar = increaseVar;
2311  }
2312  }
2313  }
2314  }
2315  if( objectiveIter.hEnd( false ) )
2316  break;
2317  else
2318  objectiveIter.hMove( false );
2319  }
2320  if( bestNonBasicVar != nullptr )
2321  {
2322  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2323  std::cout << "Smallest non-basic variable which might allow to improve the objective later: ";
2324  bestNonBasicVar->print(); std::cout << std::endl;
2325  std::cout << "Search for smallest basic variable, which is on its bound, to pivot with." << std::endl;
2326  #endif
2327  Iterator columnIter = Iterator( bestNonBasicVar->startEntry(), mpEntries );
2328  const Variable<T1, T2>* bestRowVar = nullptr;
2329  assert( columnIter.vEnd( true ) ); // Is the lowest row.
2330  assert( !columnIter.vEnd( false ) ); // There should be at least one more row containing this column variable.
2331  // Skip the objective function's row.
2332  columnIter.vMove( false );
2333  while( true )
2334  {
2335  const Variable<T1, T2>* rowVar = (*columnIter).rowVar();
2336  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2337  std::cout << "Checking " << std::endl; rowVar->print(); std::cout << " .. " << std::endl;
2338  #endif
2339  if( bestRowVar == nullptr || *rowVar < *bestRowVar )
2340  {
2341  if( ((increaseBestNonbasicVar == entryIsNegative( *columnIter )) ?
2342  (rowVar->infimum() == rowVar->assignment()) :
2343  (rowVar->supremum() == rowVar->assignment())) )
2344  {
2345  #ifdef DEBUG_NEXT_PIVOT_FOR_OPTIMIZATION
2346  std::cout << "Is smaller!" << std::endl;
2347  #endif
2348  smallestPivotingElement = columnIter.entryID();
2349  bestRowVar = rowVar;
2350  }
2351  }
2352  if( columnIter.vEnd( false ) )
2353  break;
2354  else
2355  columnIter.vMove( false );
2356  }
2357 
2358  }
2359  if( smallestPivotingElement != LAST_ENTRY_ID )
2360  *mpTheta = T1(0);
2361  return std::make_pair( smallestPivotingElement, smallestPivotingElement != LAST_ENTRY_ID );
2362  }
2363 
2364  /**
2365  * Checks if the basic variable is pivotable.
2366  * Returns LAST_ENTRY_ID as EntryID if a conflict was found. In case of no conflict, the according to the strategy best pivot nonbasic variable is returned.
2367  */
2368  template<class Settings, typename T1, typename T2>
2369  EntryID Tableau<Settings,T1,T2>::isSuitable( const Variable<T1, T2>& _basicVar, bool _forIncreasingAssignment ) const
2370  {
2371  EntryID bestEntry = LAST_ENTRY_ID;
2372  // Check all entries in the row / nonbasic variables
2373  Iterator rowIter = Iterator( _basicVar.startEntry(), mpEntries );
2374  while( true )
2375  {
2376  const Variable<T1, T2>& nonBasicVar = *(*rowIter).columnVar();
2377  if( (_forIncreasingAssignment || entryIsNegative( *rowIter )) && (!_forIncreasingAssignment || entryIsPositive( *rowIter )) )
2378  {
2379  if( nonBasicVar.supremum() > nonBasicVar.assignment() && betterEntry( rowIter.entryID(), bestEntry ) )
2380  bestEntry = rowIter.entryID(); // Nonbasic variable suitable
2381  }
2382  else
2383  {
2384  if( nonBasicVar.infimum() < nonBasicVar.assignment() && betterEntry( rowIter.entryID(), bestEntry ) )
2385  bestEntry = rowIter.entryID(); // Nonbasic variable suitable
2386  }
2387  if( rowIter.hEnd( false ) )
2388  break;
2389  else
2390  rowIter.hMove( false );
2391  }
2392  return bestEntry;
2393  }
2394 
2395  /**
2396  * Checks if the basic variable is pivotable.
2397  * Returns LAST_ENTRY_ID as EntryID if a conflict was found. In difference to isSuitable returns this function an EntryID as soon as the EntryID bestEntry
2398  * is unequal to LAST_ENTRY_ID. The assertions are still checked.
2399  */
2400  template<class Settings, typename T1, typename T2>
2401  EntryID Tableau<Settings,T1,T2>::isSuitableConflictDetection( const Variable<T1, T2>& _basicVar, bool _forIncreasingAssignment ) const
2402  {
2403  EntryID bestEntry = LAST_ENTRY_ID;
2404  // Check all entries in the row / nonbasic variables
2405  Iterator rowIter = Iterator( _basicVar.startEntry(), mpEntries );
2406  while( true )
2407  {
2408  const Variable<T1, T2>& nonBasicVar = *(*rowIter).columnVar();
2409  if( (_forIncreasingAssignment || entryIsNegative( *rowIter )) && (!_forIncreasingAssignment || entryIsPositive( *rowIter )) )
2410  {
2411  if( nonBasicVar.supremum() > nonBasicVar.assignment() && betterEntry( rowIter.entryID(), bestEntry ) ){
2412  assert( rowIter.entryID() != LAST_ENTRY_ID );
2413  return rowIter.entryID(); // Nonbasic variable suitable
2414  }
2415  }
2416  else
2417  {
2418  if( nonBasicVar.infimum() < nonBasicVar.assignment() && betterEntry( rowIter.entryID(), bestEntry ) ){
2419  assert( rowIter.entryID() != LAST_ENTRY_ID );
2420  return rowIter.entryID(); // Nonbasic variable suitable
2421  }
2422  }
2423  if( rowIter.hEnd( false ) )
2424  break;
2425  else
2426  rowIter.hMove( false );
2427  }
2428  return bestEntry;
2429  }
2430 
2431  template<class Settings, typename T1, typename T2>
2432  bool Tableau<Settings,T1,T2>::betterEntry( EntryID _isBetter, EntryID _than ) const
2433  {
2434  assert( _isBetter != LAST_ENTRY_ID );
2435  if( _than == LAST_ENTRY_ID ) return true;
2436  const Variable<T1,T2>& isBetterNbVar = *((*mpEntries)[_isBetter].columnVar());
2437  const Variable<T1,T2>& thanColumnNbVar = *((*mpEntries)[_than].columnVar());
2438  if( !Settings::use_pivoting_strategy || mPivotingSteps >= mMaxPivotsWithoutBlandsRule )
2439  return isBetterNbVar < thanColumnNbVar;
2440  if( Settings::use_activity_based_pivot_strategy )
2441  {
2442  if( isBetterNbVar.conflictActivity() < thanColumnNbVar.conflictActivity() )
2443  return true;
2444  else if( isBetterNbVar.conflictActivity() == thanColumnNbVar.conflictActivity() )
2445  {
2446  size_t valueA = boundedVariables( isBetterNbVar );
2447  size_t valueB = boundedVariables( thanColumnNbVar, valueA );
2448  if( valueA < valueB ) return true;
2449  else if( valueA == valueB )
2450  {
2451  if( isBetterNbVar.conflictActivity() < thanColumnNbVar.conflictActivity() )
2452  return true;
2453  else if( isBetterNbVar.conflictActivity() == thanColumnNbVar.conflictActivity() && isBetterNbVar.size() < thanColumnNbVar.size() )
2454  return true;
2455  }
2456  }
2457  }
2458  else
2459  {
2460  switch( Settings::nonbasic_var_choice_strategy )
2461  {
2462  case NBCS::LESS_BOUNDED_VARIABLES:
2463  {
2464  size_t valueA = boundedVariables( isBetterNbVar );
2465  size_t valueB = boundedVariables( thanColumnNbVar, valueA );
2466  if( valueA < valueB ) return true;
2467  else if( valueA == valueB )
2468  {
2469  if( isBetterNbVar.size() < thanColumnNbVar.size() ) return true;
2470  }
2471  break;
2472  }
2473  case NBCS::LESS_COLUMN_ENTRIES:
2474  {
2475  if( isBetterNbVar.size() < thanColumnNbVar.size() )
2476  return true;
2477  else if( isBetterNbVar.size() == thanColumnNbVar.size() )
2478  {
2479  size_t valueA = boundedVariables( isBetterNbVar );
2480  size_t valueB = boundedVariables( thanColumnNbVar, valueA );
2481  if( valueA < valueB ) return true;
2482  }
2483  break;
2484  }
2485  default:
2486  assert( false );
2487  }
2488  }
2489  return false;
2490  }
2491 
2492  template<class Settings, typename T1, typename T2>
2493  std::vector< const Bound<T1, T2>* > Tableau<Settings,T1,T2>::getConflict( EntryID _rowEntry ) const
2494  {
2495  assert( _rowEntry != LAST_ENTRY_ID );
2496  const Variable<T1,T2>& basicVar = *((*mpEntries)[_rowEntry].rowVar());
2497  // Upper bound is violated
2498  std::vector< const Bound<T1, T2>* > conflict;
2499  if( basicVar.supremum() < basicVar.assignment() )
2500  {
2501  auto iter = basicVar.upperbounds().rbegin();
2502  while( *iter != basicVar.pSupremum() && iter != basicVar.upperbounds().rend() )
2503  {
2504  if( (*iter)->isActive() && **iter < basicVar.assignment() )
2505  break;
2506  ++iter;
2507  }
2508  conflict.push_back( *iter );
2509  // Check all entries in the row / basic variables
2510  Iterator rowIter = Iterator( basicVar.startEntry(), mpEntries );
2511  while( true )
2512  {
2513  if( (!Settings::omit_division || ((*rowIter).content() < 0 && basicVar.factor() > 0) || ((*rowIter).content() > 0 && basicVar.factor() < 0))
2514  && (Settings::omit_division || (*rowIter).content() < 0) )
2515  {
2516  // In case of omit_division = true: check on same bound type
2517  // if division not omitted: factor set to 1, just check on negative content
2518  assert( !((*rowIter).columnVar()->supremum() > (*rowIter).columnVar()->assignment()) );
2519  conflict.push_back( (*rowIter).columnVar()->pSupremum() );
2520  }
2521  else
2522  {
2523  assert( !((*rowIter).columnVar()->infimum() < (*rowIter).columnVar()->assignment()) );
2524  conflict.push_back( (*rowIter).columnVar()->pInfimum() );
2525  }
2526  if( rowIter.hEnd( false ) )
2527  {
2528  break;
2529  }
2530  else
2531  {
2532  rowIter.hMove( false );
2533  }
2534  }
2535  }
2536  // Lower bound is violated
2537  else
2538  {
2539  assert( basicVar.infimum() > basicVar.assignment() );
2540  auto iter = basicVar.lowerbounds().begin();
2541  while( *iter != basicVar.pInfimum() && iter != basicVar.lowerbounds().end() )
2542  {
2543  if( (*iter)->isActive() && **iter > basicVar.assignment() )
2544  break;
2545  ++iter;
2546  }
2547  conflict.push_back( *iter );
2548  // Check all entries in the row / basic variables
2549  Iterator rowIter = Iterator( basicVar.startEntry(), mpEntries );
2550  while( true )
2551  {
2552  if( (!Settings::omit_division || ((*rowIter).content() > 0 && basicVar.factor() > 0) || ((*rowIter).content() < 0 && basicVar.factor() < 0))
2553  && (Settings::omit_division || (*rowIter).content() > 0) )
2554  {
2555  assert( !((*rowIter).columnVar()->supremum() > (*rowIter).columnVar()->assignment()) );
2556  conflict.push_back( (*rowIter).columnVar()->pSupremum() );
2557  }
2558  else
2559  {
2560  assert( !((*rowIter).columnVar()->infimum() < (*rowIter).columnVar()->assignment()) );
2561  conflict.push_back( (*rowIter).columnVar()->pInfimum() );
2562  }
2563  if( rowIter.hEnd( false ) )
2564  {
2565  break;
2566  }
2567  else
2568  {
2569  rowIter.hMove( false );
2570  }
2571  }
2572  }
2573  return conflict;
2574  }
2575 
2576  template<class Settings, typename T1, typename T2>
2577  std::vector< std::vector< const Bound<T1, T2>* > > Tableau<Settings,T1,T2>::getConflictsFrom( EntryID _rowEntry ) const
2578  {
2579  std::vector< std::vector< const Bound<T1, T2>* > > conflicts;
2580  const Variable<T1,T2>* firstConflictingVar = (*mpEntries)[_rowEntry].rowVar();
2581  bool posOfFirstConflictFound = false;
2582  for( Variable<T1,T2>* rowElement : mRows )
2583  {
2584  if( !posOfFirstConflictFound )
2585  {
2586  if( rowElement == firstConflictingVar )
2587  posOfFirstConflictFound = true;
2588  else
2589  continue;
2590  }
2591  assert( rowElement != NULL );
2592  // Upper bound is violated
2593  const Variable<T1,T2>& basicVar = *rowElement;
2594  if( basicVar.supremum() < basicVar.assignment() )
2595  {
2596  conflicts.emplace_back();
2597  conflicts.back().push_back( basicVar.pSupremum() );
2598  // Check all entries in the row / basic variables
2599  Iterator rowIter = Iterator( basicVar.startEntry(), mpEntries );
2600  while( true )
2601  {
2602  if( (!Settings::omit_division || ( (*rowIter).content() < 0 && basicVar.factor() > 0) || ((*rowIter).content() > 0 && basicVar.factor() < 0))
2603  && (Settings::omit_division || (*rowIter).content() < 0) )
2604  {
2605  if( (*rowIter).columnVar()->supremum() > (*rowIter).columnVar()->assignment() )
2606  {
2607  // Not a conflict.
2608  conflicts.pop_back();
2609  break;
2610  }
2611  else
2612  {
2613  conflicts.back().push_back( (*rowIter).columnVar()->pSupremum() );
2614  }
2615  }
2616  else
2617  {
2618  if( (*rowIter).columnVar()->infimum() < (*rowIter).columnVar()->assignment() )
2619  {
2620  // Not a conflict.
2621  conflicts.pop_back();
2622  break;
2623  }
2624  else
2625  {
2626  conflicts.back().push_back( (*rowIter).columnVar()->pInfimum() );
2627  }
2628  }
2629  if( rowIter.hEnd( false ) )
2630  {
2631  break;
2632  }
2633  else
2634  {
2635  rowIter.hMove( false );
2636  }
2637  }
2638  }
2639  // Lower bound is violated
2640  else if( basicVar.infimum() > basicVar.assignment() )
2641  {
2642  conflicts.emplace_back();
2643  conflicts.back().push_back( basicVar.pInfimum() );
2644  // Check all entries in the row / basic variables
2645  Iterator rowIter = Iterator( basicVar.startEntry(), mpEntries );
2646  while( true )
2647  {
2648  if( (!Settings::omit_division || ((*rowIter).content() > 0 && basicVar.factor() > 0) || ((*rowIter).content() < 0 && basicVar.factor() < 0))
2649  && (Settings::omit_division || (*rowIter).content() > 0) )
2650  {
2651  if( (*rowIter).columnVar()->supremum() > (*rowIter).columnVar()->assignment() )
2652  {
2653  // Not a conflict.
2654  conflicts.pop_back();
2655  break;
2656  }
2657  else
2658  {
2659  conflicts.back().push_back( (*rowIter).columnVar()->pSupremum() );
2660  }
2661  }
2662  else
2663  {
2664  if( (*rowIter).columnVar()->infimum() < (*rowIter).columnVar()->assignment() )
2665  {
2666  // Not a conflict.
2667  conflicts.pop_back();
2668  break;
2669  }
2670  else
2671  {
2672  conflicts.back().push_back( (*rowIter).columnVar()->pInfimum() );
2673  }
2674  }
2675  if( rowIter.hEnd( false ) )
2676  {
2677  break;
2678  }
2679  else
2680  {
2681  rowIter.hMove( false );
2682  }
2683  }
2684  }
2685  }
2686  return conflicts;
2687  }
2688 
2689  template<class Settings, typename T1, typename T2>
2690  void Tableau<Settings,T1,T2>::updateBasicAssignments( size_t _column, const Value<T1>& _change )
2691  {
2692  Variable<T1,T2>& nonbasicVar = *mColumns[_column];
2693  if( nonbasicVar.size() > 0 )
2694  {
2695  Iterator columnIter = Iterator( nonbasicVar.startEntry(), mpEntries );
2696  while( true )
2697  {
2698  Variable<T1, T2>& basic = *((*columnIter).rowVar());
2699  if( Settings::omit_division )
2700  basic.rAssignment() += (_change * (*columnIter).content())/basic.factor();
2701  else
2702  basic.rAssignment() += (_change * (*columnIter).content());
2703  if( columnIter.vEnd( false ) )
2704  {
2705  break;
2706  }
2707  else
2708  {
2709  columnIter.vMove( false );
2710  }
2711  }
2712  }
2713  }
2714 
2715  template<class Settings, typename T1, typename T2>
2716  Variable<T1, T2>* Tableau<Settings,T1,T2>::pivot( EntryID _pivotingElement, bool _optimizing )
2717  {
2718  if(onlyUpdate){
2719  ++mPivotingSteps; // alsoc count here to accredit only update operations as pivoting steps
2720  SMTRAT_LOG_DEBUG("smtrat", "Apply only update on "<<_pivotingElement << " by " << *mpTheta);
2721  updateNonbasicVariable(_pivotingElement);
2722  return (*mpEntries)[_pivotingElement].columnVar(); // Not sure what return value is used for (in current implementation it is never used)
2723  }
2724  SMTRAT_LOG_DEBUG("smtrat", "Pivot and update on "<< _pivotingElement << " by " << *mpTheta);
2725  // Find all columns having "a nonzero entry in the pivoting row"**, update this entry and store it.
2726  // First the column with ** left to the pivoting column until the leftmost column with **.
2727  std::vector<Iterator> pivotingRowLeftSide;
2728  TableauEntry<T1,T2>& pivotEntry = (*mpEntries)[_pivotingElement];
2729  T2& pivotContent = pivotEntry.rContent();
2730  Iterator iterTemp = Iterator( _pivotingElement, mpEntries );
2731  Variable<T1, T2>* rowVar = pivotEntry.rowVar();
2732  Variable<T1, T2>* columnVar = pivotEntry.columnVar();
2733  pivotEntry.setRowVar( columnVar );
2734  Iterator colIter = Iterator( columnVar->startEntry(), mpEntries );
2735  while( true )
2736  {
2737  (*colIter).setColumnVar( rowVar );
2738  if( colIter.vEnd( false ) )
2739  break;
2740  colIter.vMove( false );
2741  }
2742  while( !iterTemp.hEnd( true ) )
2743  {
2744  iterTemp.hMove( true );
2745  (*iterTemp).setRowVar( columnVar );
2746  if( Settings::omit_division )
2747  (*iterTemp).rContent() = -(*iterTemp).content();
2748  else
2749  (*iterTemp).rContent() /= -pivotContent;
2750  pivotingRowLeftSide.push_back( iterTemp );
2751  }
2752  // Then the column with ** right to the pivoting column until the rightmost column with **.
2753  std::vector<Iterator> pivotingRowRightSide = std::vector<Iterator>();
2754  iterTemp = Iterator( _pivotingElement, mpEntries );
2755  while( !iterTemp.hEnd( false ) )
2756  {
2757  iterTemp.hMove( false );
2758  (*iterTemp).setRowVar( columnVar );
2759  if( Settings::omit_division )
2760  (*iterTemp).rContent() = -(*iterTemp).content();
2761  else
2762  (*iterTemp).rContent() /= -pivotContent;
2763  pivotingRowRightSide.push_back( iterTemp );
2764  }
2765  // Swap the variables
2766  mRows[rowVar->position()] = columnVar;
2767  mColumns[columnVar->position()] = rowVar;
2768  // Update the assignments of the pivoting variables
2769  if( Settings::omit_division )
2770  {
2771  rowVar->rAssignment() += ((*mpTheta) * pivotContent) / rowVar->factor();
2772  }
2773  else
2774  rowVar->rAssignment() += (*mpTheta) * pivotContent;
2775  assert( rowVar->supremum() > rowVar->assignment() || rowVar->supremum() == rowVar->assignment() );
2776  assert( rowVar->infimum() < rowVar->assignment() || rowVar->infimum() == rowVar->assignment() );
2777  columnVar->rAssignment() += (*mpTheta);
2778  // Adapt both variables.
2779  Variable<T1, T2>& basicVar = *columnVar;
2780  Variable<T1, T2>& nonbasicVar = *rowVar;
2781  size_t tmpPosition = basicVar.position();
2782  basicVar.rPosition() = nonbasicVar.position();
2783  nonbasicVar.rPosition() = tmpPosition;
2784  size_t tmpSize = basicVar.size();
2785  basicVar.rSize() = nonbasicVar.size();
2786  nonbasicVar.rSize() = tmpSize;
2787  basicVar.setBasic( true );
2788  nonbasicVar.setBasic( false );
2789  EntryID tmpStartEntry = basicVar.startEntry();
2790  basicVar.rStartEntry() = nonbasicVar.startEntry();
2791  nonbasicVar.rStartEntry() = tmpStartEntry;
2792  // Update the content of the pivoting entry
2793  if( Settings::omit_division )
2794  {
2795  basicVar.rFactor() = pivotContent;
2796  pivotContent = nonbasicVar.factor();
2797  nonbasicVar.rFactor() = Rational(1);
2798  }
2799  else
2800  {
2801  pivotContent = carl::div( T2(1), pivotContent );
2802  }
2803  if( !_optimizing && Settings::use_refinement && basicVar.hasBound() )
2804  {
2805  rowRefinement( columnVar ); // Note, we have swapped the variables, so the current basic var is now corresponding to what we have stored in columnVar.
2806  }
2807  // Let (p_r,p_c,p_e) be the pivoting entry, where p_r is the row number, p_c the column number and p_e the content.
2808  // For all rows R having a nonzero entry in the pivoting column:
2809  // For all columns C having a nonzero entry (r_r,r_c,r_e) in the pivoting row:
2810  // Update the entry (t_r,t_c,t_e) of the intersection of R and C to (t_r,t_c,t_e+r_e).
2811  if( pivotEntry.vNext( false ) == LAST_ENTRY_ID )
2812  {
2813  update( true, _pivotingElement, pivotingRowLeftSide, pivotingRowRightSide, _optimizing );
2814  }
2815  else if( pivotEntry.vNext( true ) == LAST_ENTRY_ID )
2816  {
2817  update( false, _pivotingElement, pivotingRowLeftSide, pivotingRowRightSide, _optimizing );
2818  }
2819  else
2820  {
2821  update( true, _pivotingElement, pivotingRowLeftSide, pivotingRowRightSide, _optimizing );
2822  update( false, _pivotingElement, pivotingRowLeftSide, pivotingRowRightSide, _optimizing );
2823  }
2824  ++mPivotingSteps;
2825  if( !basicVar.hasBound() && !basicVar.isOriginal() )
2826  {
2827  deactivateBasicVar( columnVar );
2828  compressRows();
2829  if( basicVar.lowerbounds().size() == 1 && basicVar.upperbounds().size() == 1 )
2830  {
2831  mSlackVars.erase( basicVar.pExpression() );
2832  mNonActiveBasics.erase( basicVar.positionInNonActives() );
2833  basicVar.setPositionInNonActives( mNonActiveBasics.end() );
2834  assert( columnVar->isBasic() );
2835  mVariableIdAllocator.free( columnVar->id() );
2836  delete columnVar;
2837  }
2838  }
2839  else
2840  {
2841  assert( basicVar.supremum() >= basicVar.assignment() || basicVar.infimum() <= basicVar.assignment() );
2842  // assert( nonbasicVar.supremum() == nonbasicVar.assignment() || nonbasicVar.infimum() == nonbasicVar.assignment() );
2843  assert( checkCorrectness() == mRows.size() );
2844  }
2845  return columnVar;
2846  }
2847 
2848  template<class Settings, typename T1, typename T2>
2849  void Tableau<Settings,T1,T2>::update( bool _downwards, EntryID _pivotingElement, std::vector<Iterator>& _pivotingRowLeftSide, std::vector<Iterator>& _pivotingRowRightSide, bool _optimizing )
2850  {
2851  std::vector<Iterator> leftColumnIters = std::vector<Iterator>( _pivotingRowLeftSide );
2852  std::vector<Iterator> rightColumnIters = std::vector<Iterator>( _pivotingRowRightSide );
2853  Iterator pivotingColumnIter = Iterator( _pivotingElement, mpEntries );
2854  const T2& pivotingRowFactor = Settings::omit_division ? (*mpEntries)[_pivotingElement].rowVar()->factor() : T2(1);
2855  while( true )
2856  {
2857  if( !pivotingColumnIter.vEnd( _downwards ) )
2858  {
2859  pivotingColumnIter.vMove( _downwards );
2860  }
2861  else
2862  {
2863  break;
2864  }
2865  // Update the assignment of the basic variable corresponding to this row
2866  Variable<T1,T2>& currBasicVar = *((*pivotingColumnIter).rowVar());
2867  if( Settings::omit_division )
2868  currBasicVar.rAssignment() += ((*mpTheta) * (*pivotingColumnIter).content())/currBasicVar.factor();
2869  else
2870  currBasicVar.rAssignment() += (*mpTheta) * (*pivotingColumnIter).content();
2871  assert( !_optimizing || currBasicVar.infimum() < currBasicVar.assignment() || currBasicVar.infimum() == currBasicVar.assignment() );
2872  assert( !_optimizing || currBasicVar.supremum() > currBasicVar.assignment() || currBasicVar.supremum() == currBasicVar.assignment() );
2873  // Update the row
2874  Iterator currentRowIter = pivotingColumnIter;
2875  T2 ca = T2(1);
2876  T2 g = T2(1);
2877  if( Settings::omit_division )
2878  {
2879  T2 l = carl::lcm( (*pivotingColumnIter).content(), pivotingRowFactor );
2880  assert( l > 0 );
2881  if( (*pivotingColumnIter).content() < 0 && pivotingRowFactor < 0 )
2882  l *= T2( -1 );
2883  ca = carl::div( l, pivotingRowFactor );
2884  T2 cb = carl::div( l, (*pivotingColumnIter).content() );
2885  currBasicVar.rFactor() *= cb;
2886  Iterator rowIter = Iterator( currBasicVar.startEntry(), mpEntries );
2887  while( true )
2888  {
2889  (*rowIter).rContent() *= cb;
2890  if( rowIter.hEnd( false ) ) break;
2891  rowIter.hMove( false );
2892  }
2893  g = carl::abs( currBasicVar.factor() );
2894  }
2895  auto pivotingRowIter = _pivotingRowLeftSide.begin();
2896  for( auto currentColumnIter = leftColumnIters.begin(); currentColumnIter != leftColumnIters.end(); ++currentColumnIter )
2897  {
2898  assert( pivotingRowIter != _pivotingRowLeftSide.end() );
2899  while( (_downwards && !(*currentColumnIter).vEnd( _downwards ) && (**currentColumnIter).rowVar()->position() < (*pivotingColumnIter).rowVar()->position())
2900  || (!_downwards && !(*currentColumnIter).vEnd( _downwards ) && (**currentColumnIter).rowVar()->position() > (*pivotingColumnIter).rowVar()->position()) )
2901  {
2902  (*currentColumnIter).vMove( _downwards );
2903  }
2904  while( !currentRowIter.hEnd( true ) && (*currentRowIter).columnVar()->position() > (**currentColumnIter).columnVar()->position() )
2905  {
2906  currentRowIter.hMove( true );
2907  }
2908  if( Settings::omit_division )
2909  addToEntry( ca * (**pivotingRowIter).content(), currentRowIter, true, *currentColumnIter, _downwards );
2910  else
2911  addToEntry( (*pivotingColumnIter).content() * (**pivotingRowIter).content(), currentRowIter, true, *currentColumnIter, _downwards );
2912  ++pivotingRowIter;
2913  }
2914  currentRowIter = pivotingColumnIter;
2915  pivotingRowIter = _pivotingRowRightSide.begin();
2916  for( auto currentColumnIter = rightColumnIters.begin(); currentColumnIter != rightColumnIters.end(); ++currentColumnIter )
2917  {
2918  assert( pivotingRowIter != _pivotingRowRightSide.end() );
2919  while( (_downwards && !(*currentColumnIter).vEnd( _downwards ) && (**currentColumnIter).rowVar()->position() < (*pivotingColumnIter).rowVar()->position())
2920  || (!_downwards && !(*currentColumnIter).vEnd( _downwards ) && (**currentColumnIter).rowVar()->position() > (*pivotingColumnIter).rowVar()->position()) )
2921  {
2922  (*currentColumnIter).vMove( _downwards );
2923  }
2924  while( !currentRowIter.hEnd( false ) && (*currentRowIter).columnVar()->position() < (**currentColumnIter).columnVar()->position() )
2925  {
2926  currentRowIter.hMove( false );
2927  }
2928  if( Settings::omit_division )
2929  addToEntry( ca * (**pivotingRowIter).content(), currentRowIter, false, *currentColumnIter, _downwards );
2930  else
2931  addToEntry( (*pivotingColumnIter).content() * (**pivotingRowIter).content(), currentRowIter, false, *currentColumnIter, _downwards );
2932  ++pivotingRowIter;
2933  }
2934  if( Settings::omit_division )
2935  {
2936  (*pivotingColumnIter).rContent() = ca * (*mpEntries)[_pivotingElement].content();
2937  Iterator rowIter = Iterator( currBasicVar.startEntry(), mpEntries );
2938  while( !(g == 1) )
2939  {
2940  carl::gcd_assign( g, (*rowIter).content() );
2941  if( rowIter.hEnd( false ) ) break;
2942  rowIter.hMove( false );
2943  }
2944  if( !(g == 1) )
2945  {
2946  assert( g > 0 );
2947  rowIter = Iterator( currBasicVar.startEntry(), mpEntries );
2948  while( true )
2949  {
2950  carl::div_assign( (*rowIter).rContent(), g );
2951  if( rowIter.hEnd( false ) ) break;
2952  else rowIter.hMove( false );
2953  }
2954  carl::div_assign( currBasicVar.rFactor(), g );
2955  }
2956  }
2957  else
2958  (*pivotingColumnIter).rContent() *= (*mpEntries)[_pivotingElement].content();
2959  if( !_optimizing && (currBasicVar.supremum() > currBasicVar.assignment() || currBasicVar.infimum() < currBasicVar.assignment()) )
2960  {
2961  if( Settings::pivot_into_local_conflict )
2962  {
2963  mConflictingRows.push_back( (*pivotingColumnIter).rowVar() );
2964  }
2965  if( Settings::use_refinement )
2966  {
2967  rowRefinement( (*pivotingColumnIter).rowVar() );
2968  }
2969  }
2970  }
2971  }
2972 
2973  template<class Settings, typename T1, typename T2>
2974  void Tableau<Settings,T1,T2>::addToEntry( const T2& _toAdd, Iterator& _horiIter, bool _horiIterLeftFromVertIter, Iterator& _vertIter, bool _vertIterBelowHoriIter )
2975  {
2976  if( _horiIter == _vertIter )
2977  {
2978  // Entry already exists, so update it only and maybe remove it.
2979  T2& currentRowContent = (*_horiIter).rContent();
2980  currentRowContent += _toAdd;
2981  if( currentRowContent == 0 )
2982  {
2983  EntryID toRemove = _horiIter.entryID();
2984  _vertIter.vMove( !_vertIterBelowHoriIter );
2985  _horiIter.hMove( !_horiIterLeftFromVertIter );
2986  removeEntry( toRemove );
2987  }
2988  }
2989  else
2990  {
2991  EntryID entryID = newTableauEntry( _toAdd );
2992  TableauEntry<T1,T2>& entry = (*mpEntries)[entryID];
2993  // Set the position.
2994  Variable<T1,T2>* basicVar = (*mpEntries)[_horiIter.entryID()].rowVar();
2995  Variable<T1,T2>* nonbasicVar = (*mpEntries)[_vertIter.entryID()].columnVar();
2996  entry.setRowVar( basicVar );
2997  entry.setColumnVar( nonbasicVar );
2998  if( (_vertIterBelowHoriIter && (*_vertIter).rowVar()->position() > (*_horiIter).rowVar()->position())
2999  || (!_vertIterBelowHoriIter && (*_vertIter).rowVar()->position() < (*_horiIter).rowVar()->position()) )
3000  {
3001  // Entry vertically between two entries.
3002  EntryID upperEntryID = (*_vertIter).vNext( !_vertIterBelowHoriIter );
3003  if( upperEntryID != LAST_ENTRY_ID )
3004  {
3005  (*mpEntries)[upperEntryID].setVNext( _vertIterBelowHoriIter, entryID );
3006  }
3007  (*_vertIter).setVNext( !_vertIterBelowHoriIter, entryID );
3008  entry.setVNext( !_vertIterBelowHoriIter, upperEntryID );
3009  entry.setVNext( _vertIterBelowHoriIter, _vertIter.entryID() );
3010  }
3011  else
3012  {
3013  // Entry will be the lowest in this column.
3014  (*_vertIter).setVNext( _vertIterBelowHoriIter, entryID );
3015  entry.setVNext( !_vertIterBelowHoriIter, _vertIter.entryID() );
3016  entry.setVNext( _vertIterBelowHoriIter, LAST_ENTRY_ID );
3017  if( _vertIterBelowHoriIter )
3018  nonbasicVar->rStartEntry() = entryID;
3019  }
3020  if( (_horiIterLeftFromVertIter && (*_horiIter).columnVar()->position() < (*_vertIter).columnVar()->position())
3021  || (!_horiIterLeftFromVertIter && (*_horiIter).columnVar()->position() > (*_vertIter).columnVar()->position()) )
3022  {
3023  // Entry horizontally between two entries.
3024  EntryID rightEntryID = (*_horiIter).hNext( !_horiIterLeftFromVertIter );
3025  if( rightEntryID != LAST_ENTRY_ID )
3026  {
3027  (*mpEntries)[rightEntryID].setHNext( _horiIterLeftFromVertIter, entryID );
3028  }
3029  (*_horiIter).setHNext( !_horiIterLeftFromVertIter, entryID );
3030  entry.setHNext( !_horiIterLeftFromVertIter, rightEntryID );
3031  entry.setHNext( _horiIterLeftFromVertIter, _horiIter.entryID() );
3032  }
3033  else
3034  {
3035  // Entry will be the leftmost in this row.
3036  (*_horiIter).setHNext( _horiIterLeftFromVertIter, entryID );
3037  entry.setHNext( !_horiIterLeftFromVertIter, _horiIter.entryID() );
3038  entry.setHNext( _horiIterLeftFromVertIter, LAST_ENTRY_ID );
3039  if( _horiIterLeftFromVertIter )
3040  basicVar->rStartEntry() = entryID;
3041  }
3042  // Set the content of the entry.
3043  ++(basicVar->rSize());
3044  ++(nonbasicVar->rSize());
3045  }
3046  }
3047 
3048  template<class Settings, typename T1, typename T2>
3049  void Tableau<Settings,T1,T2>::rowRefinement( Variable<T1,T2>* _basicVar )
3050  {
3051  // Collect the bounds which form an upper resp. lower refinement.
3052  Variable<T1,T2>& basicVar = *_basicVar;
3053  if( basicVar.size() > 128 ) return;
3054  std::vector<const Bound<T1, T2>*>* uPremise = new std::vector<const Bound<T1, T2>*>();
3055  std::vector<const Bound<T1, T2>*>* lPremise = new std::vector<const Bound<T1, T2>*>();
3056  Iterator rowEntry = Iterator( basicVar.startEntry(), mpEntries );
3057  const T2& rowFactor = Settings::omit_division ? basicVar.factor() : T2(1);
3058  while( true )
3059  {
3060  if( (!Settings::omit_division || ((*rowEntry).content() > 0 && rowFactor > 0) || ((*rowEntry).content() < 0 && rowFactor < 0))
3061  && (Settings::omit_division || (*rowEntry).content() > 0) )
3062  {
3063  if( uPremise != NULL )
3064  {
3065  const Bound<T1, T2>* sup = (*rowEntry).columnVar()->pSupremum();
3066  if( sup->pLimit() != NULL )
3067  {
3068  uPremise->push_back( sup );
3069  }
3070  else
3071  {
3072  delete uPremise;
3073  uPremise = NULL;
3074  if( lPremise == NULL ) return;
3075  }
3076  }
3077  if( lPremise != NULL )
3078  {
3079  const Bound<T1, T2>* inf = (*rowEntry).columnVar()->pInfimum();
3080  if( inf->pLimit() != NULL )
3081  {
3082  lPremise->push_back( inf );
3083  }
3084  else
3085  {
3086  delete lPremise;
3087  lPremise = NULL;
3088  if( uPremise == NULL ) return;
3089  }
3090  }
3091  }
3092  else
3093  {
3094  if( uPremise != NULL )
3095  {
3096  const Bound<T1, T2>* inf = (*rowEntry).columnVar()->pInfimum();
3097  if( inf->pLimit() != NULL )
3098  {
3099  uPremise->push_back( inf );
3100  }
3101  else
3102  {
3103  delete uPremise;
3104  uPremise = NULL;
3105  if( lPremise == NULL ) return;
3106  }
3107  }
3108  if( lPremise != NULL )
3109  {
3110  const Bound<T1, T2>* sup = (*rowEntry).columnVar()->pSupremum();
3111  if( sup->pLimit() != NULL )
3112  {
3113  lPremise->push_back( sup );
3114  }
3115  else
3116  {
3117  delete lPremise;
3118  lPremise = NULL;
3119  if( uPremise == NULL ) return;
3120  }
3121  }
3122  }
3123  if( rowEntry.hEnd( false ) ) break;
3124  else rowEntry.hMove( false );
3125  }
3126  if( uPremise != NULL )
3127  {
3128  // Found an upper refinement.
3129  Value<T1>* newlimit = new Value<T1>();
3130  typename std::vector< const Bound<T1, T2>* >::iterator bound = uPremise->begin();
3131  Iterator rowEntry = Iterator( basicVar.startEntry(), mpEntries );
3132  while( true )
3133  {
3134  *newlimit += (*bound)->limit() * (*rowEntry).content();
3135  ++bound;
3136  if( !rowEntry.hEnd( false ) ) rowEntry.hMove( false );
3137  else break;
3138  }
3139  // Learn that the strongest weaker upper bound should be activated.
3140  const typename Bound<T1, T2>::BoundSet& upperBounds = basicVar.upperbounds();
3141  auto ubound = upperBounds.begin();
3142  while( ubound != upperBounds.end() )
3143  {
3144  if( (!Settings::omit_division || (**ubound > (*newlimit)/rowFactor && (*ubound)->type() != Bound<T1, T2>::EQUAL && !(*ubound)->deduced()))
3145  && (Settings::omit_division || (**ubound > *newlimit && (*ubound)->type() != Bound<T1, T2>::EQUAL && !(*ubound)->deduced())) )
3146  {
3147  break;
3148  }
3149  if( *ubound == basicVar.pSupremum() )
3150  {
3151  delete newlimit;
3152  delete uPremise;
3153  goto CheckLowerPremise;
3154  }
3155  ++ubound;
3156  }
3157  if( ubound != --upperBounds.end() )
3158  {
3159  assert( ((*ubound)->type() != Bound<T1, T2>::EQUAL) );
3160  LearnedBound learnedBound( *newlimit, NULL, ubound, std::move( *uPremise ) );
3161  delete uPremise;
3162  if( Settings::introduce_new_constraint_in_refinement )
3163  {
3164  if( (!Settings::omit_division || newlimit->mainPart() > (*ubound)->limit().mainPart()*rowFactor || (*ubound)->limit().deltaPart() == 0)
3165  && (Settings::omit_division || newlimit->mainPart() > (*ubound)->limit().mainPart() || (*ubound)->limit().deltaPart() == 0) )
3166  {
3167  typename Poly::PolyType lhs = Settings::omit_division ?
3168  (*ubound)->variable().expression()*(Rational)rowFactor - (Rational)newlimit->mainPart() :
3169  (*ubound)->variable().expression() - (Rational)newlimit->mainPart();
3170  carl::Relation rel = newlimit->deltaPart() != 0 ? carl::Relation::LESS : carl::Relation::LEQ;
3171  FormulaT constraint = FormulaT( smtrat::ConstraintT( lhs, rel ) );
3172  learnedBound.newBound = basicVar.addUpperBound( newlimit, mDefaultBoundPosition, constraint, true ).first;
3173  }
3174  else
3175  {
3176  learnedBound.newBound = NULL;
3177  }
3178  }
3179  else
3180  {
3181  delete newlimit;
3182  learnedBound.newBound = NULL;
3183  }
3184  auto iter = mLearnedUpperBounds.find( _basicVar );
3185  if( iter != mLearnedUpperBounds.end() )
3186  {
3187  if( *learnedBound.nextWeakerBound < *iter->second.nextWeakerBound )
3188  {
3189  iter->second.nextWeakerBound = learnedBound.nextWeakerBound;
3190  iter->second.premise = learnedBound.premise;
3191  mNewLearnedBounds.push_back( iter );
3192  }
3193  }
3194  else
3195  {
3196  iter = mLearnedUpperBounds.emplace( _basicVar, std::move(learnedBound) ).first;
3197  mNewLearnedBounds.push_back( iter );
3198  }
3199  }
3200  else
3201  {
3202  delete newlimit;
3203  delete uPremise;
3204  }
3205  }
3206  CheckLowerPremise:
3207  if( lPremise != NULL )
3208  {
3209  // Found an lower refinement.
3210  Value<T1>* newlimit = new Value<T1>();
3211  typename std::vector< const Bound<T1, T2>* >::iterator bound = lPremise->begin();
3212  Iterator rowEntry = Iterator( basicVar.startEntry(), mpEntries );
3213  while( true )
3214  {
3215  *newlimit += (*bound)->limit() * (*rowEntry).content();
3216  ++bound;
3217  if( !rowEntry.hEnd( false ) ) rowEntry.hMove( false );
3218  else break;
3219  }
3220  // Learn that the strongest weaker lower bound should be activated.
3221  const typename Bound<T1, T2>::BoundSet& lowerBounds = basicVar.lowerbounds();
3222  assert( !lowerBounds.empty() );
3223  auto lbound = --lowerBounds.end();
3224  while( true )
3225  {
3226  if( (!Settings::omit_division || (**lbound < (*newlimit)/rowFactor && (*lbound)->type() != Bound<T1, T2>::EQUAL && !(*lbound)->deduced()))
3227  && (Settings::omit_division || (**lbound < *newlimit && (*lbound)->type() != Bound<T1, T2>::EQUAL && !(*lbound)->deduced())) )
3228  {
3229  break;
3230  }
3231  if( *lbound == basicVar.pInfimum() )
3232  {
3233  delete newlimit;
3234  delete lPremise;
3235  return;
3236  }
3237  if( lbound == lowerBounds.begin() )
3238  break;
3239  --lbound;
3240  }
3241 
3242  if( lbound != lowerBounds.begin() )
3243  {
3244  assert( ((*lbound)->type() != Bound<T1, T2>::EQUAL) );
3245  LearnedBound learnedBound( *newlimit, NULL, lbound, std::move( *lPremise ) );
3246  delete lPremise;
3247  if( Settings::introduce_new_constraint_in_refinement )
3248  {
3249  if( (!Settings::omit_division || newlimit->mainPart() > (*lbound)->limit().mainPart()*rowFactor || (*lbound)->limit().deltaPart() == 0)
3250  && (Settings::omit_division || newlimit->mainPart() > (*lbound)->limit().mainPart() || (*lbound)->limit().deltaPart() == 0) )
3251  {
3252  typename Poly::PolyType lhs = Settings::omit_division ?
3253  (*lbound)->variable().expression()*(Rational)rowFactor - (Rational)newlimit->mainPart() :
3254  (*lbound)->variable().expression() - (Rational)newlimit->mainPart();
3255  carl::Relation rel = newlimit->deltaPart() != 0 ? carl::Relation::GREATER : carl::Relation::GEQ;
3256  FormulaT constraint = FormulaT( smtrat::ConstraintT( lhs, rel ) );
3257  learnedBound.newBound = basicVar.addLowerBound( newlimit, mDefaultBoundPosition, constraint, true ).first;
3258  }
3259  else
3260  {
3261  learnedBound.newBound = NULL;
3262  }
3263  }
3264  else
3265  {
3266  delete newlimit;
3267  learnedBound.newBound = NULL;
3268  }
3269  auto iter = mLearnedLowerBounds.find( _basicVar );
3270  if( iter != mLearnedLowerBounds.end() )
3271  {
3272  if( *learnedBound.nextWeakerBound > *iter->second.nextWeakerBound )
3273  {
3274  iter->second.nextWeakerBound = learnedBound.nextWeakerBound;
3275  iter->second.premise = learnedBound.premise;
3276  mNewLearnedBounds.push_back( iter );
3277  }
3278  }
3279  else
3280  {
3281  iter = mLearnedLowerBounds.emplace( _basicVar, std::move(learnedBound) ).first;
3282  mNewLearnedBounds.push_back( iter );
3283  }
3284  }
3285  else
3286  {
3287  delete newlimit;
3288  delete lPremise;
3289  }
3290  }
3291  }
3292 
3293  template<class Settings, typename T1, typename T2>
3294  size_t Tableau<Settings,T1,T2>::unboundedVariables( const Variable<T1,T2>& _var, size_t _stopCriterium ) const
3295  {
3296  if( _var.startEntry() == LAST_ENTRY_ID )
3297  {
3298  return 0;
3299  }
3300  else
3301  {
3302  size_t unboundedVars = 0;
3303  if( _var.isBasic() )
3304  {
3305  Iterator rowEntry = Iterator( _var.startEntry(), mpEntries );
3306  while( true )
3307  {
3308  if( (*rowEntry).columnVar()->infimum().isInfinite() || (*rowEntry).columnVar()->supremum().isInfinite() )
3309  {
3310  ++unboundedVars;
3311  if( _stopCriterium != 0 && unboundedVars >= _stopCriterium )
3312  return _stopCriterium + 1;
3313  }
3314  if( rowEntry.hEnd( false ) )
3315  break;
3316  rowEntry.hMove( false );
3317  }
3318  }
3319  else
3320  {
3321  Iterator columnEntry = Iterator( _var.startEntry(), mpEntries );
3322  while( true )
3323  {
3324  if( (*columnEntry).rowVar()->infimum().isInfinite() || (*columnEntry).rowVar()->supremum().isInfinite() )
3325  {
3326  ++unboundedVars;
3327  if( _stopCriterium != 0 && unboundedVars >= _stopCriterium )
3328  return _stopCriterium + 1;
3329  }
3330  if( columnEntry.vEnd( false ) )
3331  break;
3332  columnEntry.vMove( false );
3333  }
3334  }
3335  return unboundedVars;
3336  }
3337  }
3338 
3339  template<class Settings, typename T1, typename T2>
3340  size_t Tableau<Settings,T1,T2>::boundedVariables( const Variable<T1,T2>& _var, size_t _stopCriterium ) const
3341  {
3342  if( _var.startEntry() == LAST_ENTRY_ID )
3343  {
3344  return 0;
3345  }
3346  else
3347  {
3348  size_t boundedVars = 0;
3349  if( _var.isBasic() )
3350  {
3351  Iterator rowEntry = Iterator( _var.startEntry(), mpEntries );
3352  while( true )
3353  {
3354  if( !(*rowEntry).columnVar()->infimum().isInfinite() || !(*rowEntry).columnVar()->supremum().isInfinite() )
3355  {
3356  ++boundedVars;
3357  if( _stopCriterium != 0 && boundedVars >= _stopCriterium )
3358  return _stopCriterium+1;
3359  }
3360  if( rowEntry.hEnd( false ) )
3361  break;
3362  rowEntry.hMove( false );
3363  }
3364  }
3365  else
3366  {
3367  Iterator columnEntry = Iterator( _var.startEntry(), mpEntries );
3368  while( true )
3369  {
3370  if( !(*columnEntry).rowVar()->infimum().isInfinite() || !(*columnEntry).rowVar()->supremum().isInfinite() )
3371  {
3372  ++boundedVars;
3373  if( _stopCriterium != 0 && boundedVars >= _stopCriterium )
3374  return _stopCriterium+1;
3375  }
3376  if( columnEntry.vEnd( false ) )
3377  break;
3378  columnEntry.vMove( false );
3379  }
3380  }
3381  return boundedVars;
3382  }
3383  }
3384 
3385  template<class Settings, typename T1, typename T2>
3386  size_t Tableau<Settings,T1,T2>::checkCorrectness() const
3387  {
3388  #ifdef LRA_PEDANTIC_CORRECTNESS_CHECKS
3389  size_t rowNumber = 0;
3390  for( ; rowNumber < mRows.size(); ++rowNumber )
3391  {
3392  if( !rowCorrect( rowNumber ) ) return rowNumber;
3393  }
3394  return rowNumber;
3395  #else
3396  return mRows.size();
3397  #endif
3398  }
3399 
3400  template<class Settings, typename T1, typename T2>
3401  bool Tableau<Settings,T1,T2>::rowCorrect( size_t _rowNumber ) const
3402  {
3403  if( mRows[_rowNumber] == NULL ) return false;
3404  if( _rowNumber != mRows[_rowNumber]->position() ) return false;
3405  size_t numOfRowElements = 0;
3406  typename Poly::PolyType sumOfNonbasics;
3407  Iterator rowEntry = Iterator( mRows[_rowNumber]->startEntry(), mpEntries );
3408  while( !rowEntry.hEnd( false ) )
3409  {
3410  sumOfNonbasics += (*((*rowEntry).columnVar()->pExpression())) * typename Poly::PolyType( (*rowEntry).content() );
3411  ++numOfRowElements;
3412  rowEntry.hMove( false );
3413  }
3414  ++numOfRowElements;
3415  if( numOfRowElements != mRows[_rowNumber]->size() )
3416  {
3417  return false;
3418  }
3419  sumOfNonbasics += (*((*rowEntry).columnVar()->pExpression())) * typename Poly::PolyType( (*rowEntry).content() );
3420  if( Settings::omit_division )
3421  sumOfNonbasics += -(*mRows[_rowNumber]->pExpression()) * typename Poly::PolyType( mRows[_rowNumber]->factor() );
3422  else
3423  sumOfNonbasics += -(*mRows[_rowNumber]->pExpression());
3424  if( !carl::is_zero(sumOfNonbasics) )
3425  {
3426  return false;
3427  }
3428  return true;
3429  }
3430 
3431  template<class Settings, typename T1, typename T2>
3432  bool Tableau<Settings,T1,T2>::isConflicting() const
3433  {
3434  for( Variable<T1,T2>* rowVar : mRows )
3435  {
3436  if( rowVar != NULL )
3437  {
3438  if( rowVar->isConflicting() ) return true;
3439  }
3440  }
3441  for( Variable<T1,T2>* columnVar : mColumns )
3442  {
3443  if( columnVar->isConflicting() ) return true;
3444  }
3445  return false;
3446  }
3447 
3448  enum GOMORY_SET
3449  {
3450  J_PLUS,
3451  J_MINUS,
3452  K_PLUS,
3453  K_MINUS
3454  };
3455 
3456  template<class Settings, typename T1, typename T2>
3457  const typename Poly::PolyType* Tableau<Settings,T1,T2>::gomoryCut( const T2& _ass, Variable<T1,T2>* _rowVar )
3458  {
3459  typename Poly::PolyType* sum = new typename Poly::PolyType();
3460  Iterator row_iterator = Iterator( _rowVar->startEntry(), mpEntries );
3461  std::vector<GOMORY_SET> splitting;
3462  // Check, whether the premises for a Gomory Cut are satisfied
3463  while( true )
3464  {
3465  const Variable<T1, T2>& nonBasicVar = *(*row_iterator).columnVar();
3466  if( !nonBasicVar.infimum().isInfinite() && nonBasicVar.infimum() == nonBasicVar.assignment() )
3467  {
3468  if( ((*row_iterator).content() < 0 && _rowVar->factor() > 0) || ((*row_iterator).content() > 0 && _rowVar->factor() < 0) )
3469  {
3470  splitting.push_back( J_MINUS );
3471  }
3472  else
3473  {
3474  splitting.push_back( J_PLUS );
3475  }
3476  }
3477  else if( !nonBasicVar.supremum().isInfinite() && nonBasicVar.supremum() == nonBasicVar.assignment() )
3478  {
3479  if( ((*row_iterator).content() < 0 && _rowVar->factor() > 0) || ((*row_iterator).content() > 0 && _rowVar->factor() < 0) )
3480  {
3481  splitting.push_back( K_MINUS );
3482  }
3483  else
3484  {
3485  splitting.push_back( K_PLUS );
3486  }
3487  }
3488  else
3489  {
3490  return sum;
3491  }
3492  if( row_iterator.hEnd( false ) )
3493  {
3494  break;
3495  }
3496  row_iterator.hMove( false );
3497  }
3498  // A Gomory Cut can be constructed
3499  T2 coeff;
3500  #ifdef LRA_DEBUG_GOMORY_CUT
3501  std::cout << "Create Cut for: " << _rowVar->expression() << std::endl;
3502  std::cout << "_ass = " << _ass << std::endl;
3503  #endif
3504  T2 f_zero = _ass - T2(carl::floor( (Rational)_ass ));
3505  #ifdef LRA_DEBUG_GOMORY_CUT
3506  std::cout << "f_zero = " << f_zero << std::endl;
3507  #endif
3508  // Construction of the Gomory Cut
3509  std::vector<GOMORY_SET>::const_iterator vec_iter = splitting.begin();
3510  row_iterator = Iterator( _rowVar->startEntry(), mpEntries );
3511  while( true )
3512  {
3513  const Variable<T1, T2>& nonBasicVar = *(*row_iterator).columnVar();
3514  if( (*vec_iter) == J_MINUS )
3515  {
3516  #ifdef LRA_DEBUG_GOMORY_CUT
3517  std::cout << "nonBasicVar.expression() = " << nonBasicVar.expression() << std::endl;
3518  std::cout << "(*row_iterator).content() = " << (*row_iterator).content() << std::endl;
3519  std::cout << "_rowVar->factor() = " << _rowVar->factor() << std::endl;
3520  std::cout << "f_zero = " << f_zero << std::endl;
3521  std::cout << "(_rowVar->factor() * f_zero) = " << (_rowVar->factor() * f_zero) << std::endl;
3522  #endif
3523  coeff = -((*row_iterator).content() / (_rowVar->factor() * f_zero));
3524  #ifdef LRA_DEBUG_GOMORY_CUT
3525  std::cout << "A: coeff = " << coeff << std::endl;
3526  #endif
3527  *sum += (Rational)coeff*( nonBasicVar.expression() - (Rational)nonBasicVar.infimum().limit().mainPart() );
3528  #ifdef LRA_DEBUG_GOMORY_CUT
3529  std::cout << "(Rational)coeff*( nonBasicVar.expression() - (Rational)nonBasicVar.infimum().limit().mainPart() ) = " << ((Rational)coeff*( nonBasicVar.expression() - (Rational)nonBasicVar.infimum().limit().mainPart() )) << std::endl;
3530  #endif
3531  }
3532  else if( (*vec_iter) == J_PLUS )
3533  {
3534  #ifdef LRA_DEBUG_GOMORY_CUT
3535  std::cout << "nonBasicVar.expression() = " << nonBasicVar.expression() << std::endl;
3536  std::cout << "(*row_iterator).content() = " << (*row_iterator).content() << std::endl;
3537  std::cout << "_rowVar->factor() = " << _rowVar->factor() << std::endl;
3538  std::cout << "f_zero = " << f_zero << std::endl;
3539  std::cout << "(_rowVar->factor() * ( (Rational)1 - f_zero )) = " << (_rowVar->factor() * ( (Rational)1 - f_zero )) << std::endl;
3540  #endif
3541  coeff = (*row_iterator).content()/(_rowVar->factor() * Rational(Rational(1) - Rational(f_zero)));
3542  #ifdef LRA_DEBUG_GOMORY_CUT
3543  std::cout << "B: coeff = " << coeff << std::endl;
3544  #endif
3545  *sum += (Rational)coeff*( nonBasicVar.expression() - (Rational)nonBasicVar.infimum().limit().mainPart() );
3546  #ifdef LRA_DEBUG_GOMORY_CUT
3547  std::cout << "(Rational)coeff*( nonBasicVar.expression() - (Rational)nonBasicVar.infimum().limit().mainPart() ) = " << ((Rational)coeff*( nonBasicVar.expression() - (Rational)nonBasicVar.infimum().limit().mainPart() )) << std::endl;
3548  #endif
3549  }
3550  else if( (*vec_iter) == K_MINUS )
3551  {
3552  #ifdef LRA_DEBUG_GOMORY_CUT
3553  std::cout << "nonBasicVar.expression() = " << nonBasicVar.expression() << std::endl;
3554  std::cout << "(*row_iterator).content() = " << (*row_iterator).content() << std::endl;
3555  std::cout << "_rowVar->factor() = " << _rowVar->factor() << std::endl;
3556  std::cout << "f_zero = " << f_zero << std::endl;
3557  std::cout << "(_rowVar->factor() * ( (Rational)1 - f_zero )) = " << (_rowVar->factor() * (Rational(1) - Rational(f_zero))) << std::endl;
3558  #endif
3559  coeff = -((*row_iterator).content()/(_rowVar->factor() * (Rational(1) - Rational(f_zero))));
3560  #ifdef LRA_DEBUG_GOMORY_CUT
3561  std::cout << "C: coeff = " << coeff << std::endl;
3562  #endif
3563  *sum += ((Rational)-coeff) * ( nonBasicVar.expression() - (Rational)nonBasicVar.supremum().limit().mainPart() );
3564  #ifdef LRA_DEBUG_GOMORY_CUT
3565  std::cout << "(Rational)coeff * ( (Rational)nonBasicVar.supremum().limit().mainPart() - nonBasicVar.expression() ) = " << ((Rational)coeff * ( (Rational)nonBasicVar.supremum().limit().mainPart() - nonBasicVar.expression() )) << std::endl;
3566  #endif
3567  }
3568  else// if( (*vec_iter) == K_PLUS )
3569  {
3570  #ifdef LRA_DEBUG_GOMORY_CUT
3571  std::cout << "nonBasicVar.expression() = " << nonBasicVar.expression() << std::endl;
3572  std::cout << "(*row_iterator).content() = " << (*row_iterator).content() << std::endl;
3573  std::cout << "_rowVar->factor() = " << _rowVar->factor() << std::endl;
3574  std::cout << "f_zero = " << f_zero << std::endl;
3575  std::cout << "(_rowVar->factor() * f_zero) = " << (_rowVar->factor() * f_zero) << std::endl;
3576  #endif
3577  coeff = (*row_iterator).content()/(_rowVar->factor() * f_zero);
3578  #ifdef LRA_DEBUG_GOMORY_CUT
3579  std::cout << "D: coeff = " << coeff << std::endl;
3580  #endif
3581  *sum += ((Rational)-coeff) * (nonBasicVar.expression() - (Rational)nonBasicVar.supremum().limit().mainPart());
3582  #ifdef LRA_DEBUG_GOMORY_CUT
3583  std::cout << "(Rational)coeff * ( (Rational)nonBasicVar.supremum().limit().mainPart() - nonBasicVar.expression() ) = " << ((Rational)coeff * ( (Rational)nonBasicVar.supremum().limit().mainPart() - nonBasicVar.expression() )) << std::endl;
3584  #endif
3585  }
3586  if( row_iterator.hEnd( false ) )
3587  {
3588  break;
3589  }
3590  row_iterator.hMove( false );
3591  ++vec_iter;
3592  }
3593  *sum -= (Rational)1;
3594  #ifdef LRA_DEBUG_GOMORY_CUT
3595  std::cout << "sum = " << sum << std::endl;
3596  #endif
3597  //const Constraint* gomory_constr = newConstraint( *sum , carl::Relation::GEQ );
3598  //newBound(gomory_constr);
3599  // TODO: check whether there is already a basic variable with this polynomial (psum, cf. LRAModule::initialize(..))
3600  return sum;
3601  }
3602 
3603  template<class Settings, typename T1, typename T2>
3604  size_t Tableau<Settings,T1,T2>::getNumberOfEntries( Variable<T1,T2>* _rowVar )
3605  {
3606  size_t result = 0;
3607  Iterator row_iterator = Iterator( _rowVar->startEntry(), mpEntries );
3608  while( true )
3609  {
3610  ++result;
3611  if( row_iterator.hEnd( false ) )
3612  {
3613  row_iterator.hMove( false );
3614  }
3615  else
3616  {
3617  return result;
3618  }
3619  }
3620  }
3621 
3622  template<class Settings, typename T1, typename T2>
3623  void Tableau<Settings,T1,T2>::collect_premises( const Variable<T1,T2>* _rowVar, FormulasT& premises ) const
3624  {
3625  Iterator row_iterator = Iterator( _rowVar->startEntry(), mpEntries );
3626  while( true )
3627  {
3628  const Variable<T1, T2>& nonBasicVar = *(*row_iterator).columnVar();
3629  if( nonBasicVar.infimum() == nonBasicVar.assignment() )
3630  {
3631  premises.push_back( (*row_iterator).columnVar()->infimum().origins().front() );
3632  }
3633  else if( nonBasicVar.supremum() == nonBasicVar.assignment() )
3634  {
3635  premises.push_back( (*row_iterator).columnVar()->supremum().origins().front() );
3636  }
3637  if( !row_iterator.hEnd( false ) )
3638  {
3639  row_iterator.hMove( false );
3640  }
3641  else
3642  {
3643  return;
3644  }
3645  }
3646  }
3647 
3648  #ifdef DEBUG_METHODS_TABLEAU
3649 
3650  template<class Settings, typename T1, typename T2>
3651  void Tableau<Settings,T1,T2>::printHeap( std::ostream& _out, int _maxEntryLength, const std::string _init ) const
3652  {
3653  for( EntryID pos = 1; pos < mpEntries->size(); ++pos )
3654  {
3655  std::cout << _init;
3656  printEntry( pos, _out, _maxEntryLength );
3657  _out << std::endl;
3658  }
3659  }
3660 
3661  template<class Settings, typename T1, typename T2>
3662  void Tableau<Settings,T1,T2>::printEntry( EntryID _entry, std::ostream& _out, int _maxEntryLength ) const
3663  {
3664  _out << std::setw( 4 ) << _entry << ": ";
3665  std::stringstream out;
3666  if( _entry != LAST_ENTRY_ID )
3667  {
3668  out << (*mpEntries)[_entry].content();
3669  _out << std::setw( _maxEntryLength ) << out.str();
3670  }
3671  else
3672  {
3673  _out << std::setw( _maxEntryLength ) << "NULL";
3674  }
3675  _out << " at (";
3676  _out << std::setw( 4 ) << ((_entry == 0 || (*mpEntries)[_entry].rowVar() == NULL) ? 0 : (*mpEntries)[_entry].rowVar()->position());
3677  _out << ",";
3678  _out << std::setw( 4 ) << (_entry == 0 ? 0 : (*mpEntries)[_entry].columnVar()->position());
3679  _out << ") [up:";
3680  _out << std::setw( 4 ) << (*mpEntries)[_entry].vNext( false );
3681  _out << ", down:";
3682  _out << std::setw( 4 ) << (*mpEntries)[_entry].vNext( true );
3683  _out << ", left:";
3684  _out << std::setw( 4 ) << (*mpEntries)[_entry].hNext( true );
3685  _out << ", right:";
3686  _out << std::setw( 4 ) << (*mpEntries)[_entry].hNext( false );
3687  _out << "]";
3688  if( _entry != 0 && (*mpEntries)[_entry].rowVar() != NULL )
3689  {
3690  _out << " (" << *(*mpEntries)[_entry].rowVar()->pExpression() << ", " << *(*mpEntries)[_entry].columnVar()->pExpression() << ")";
3691  }
3692  }
3693 
3694  template<class Settings, typename T1, typename T2>
3695  void Tableau<Settings,T1,T2>::printVariables( bool _allBounds, std::ostream& _out, const std::string _init ) const
3696  {
3697  _out << _init << "Basic variables:" << std::endl;
3698  for( Variable<T1,T2>* rowVar : mRows )
3699  {
3700  if( rowVar != NULL )
3701  {
3702  _out << _init << " ";
3703  rowVar->print( _out );
3704  _out << "(" << unboundedVariables( *rowVar ) << ")" << std::endl;
3705  if( _allBounds ) rowVar->printAllBounds( _out, _init + " " );
3706  }
3707  }
3708  _out << _init << "Nonbasic variables:" << std::endl;
3709  for( Variable<T1,T2>* columnVar : mColumns )
3710  {
3711  _out << _init << " ";
3712  columnVar->print( _out );
3713  _out << "(" << unboundedVariables( *columnVar ) << ")" << std::endl;
3714  if( _allBounds ) columnVar->printAllBounds( _out, _init + " " );
3715  }
3716  }
3717 
3718  template<class Settings, typename T1, typename T2>
3719  void Tableau<Settings,T1,T2>::printLearnedBounds( const std::string _init, std::ostream& _out ) const
3720  {
3721  for( auto learnedBound = mLearnedLowerBounds.begin(); learnedBound != mLearnedLowerBounds.end(); ++learnedBound )
3722  {
3723  printLearnedBound( *learnedBound->first, learnedBound->second, _init, _out );
3724  }
3725  for( auto learnedBound = mLearnedUpperBounds.begin(); learnedBound != mLearnedUpperBounds.end(); ++learnedBound )
3726  {
3727  printLearnedBound( *learnedBound->first, learnedBound->second, _init, _out );
3728  }
3729  }
3730 
3731  template<class Settings, typename T1, typename T2>
3732  void Tableau<Settings,T1,T2>::printLearnedBound( const Variable<T1,T2>& _var, const LearnedBound& _learnedBound, const std::string _init, std::ostream& _out ) const
3733  {
3734  for( auto premiseBound = _learnedBound.premise->begin(); premiseBound != _learnedBound.premise->end(); ++premiseBound )
3735  {
3736  _out << _init;
3737  _out << *(*premiseBound)->variable().pExpression();
3738  (*premiseBound)->print( true, _out, true );
3739  _out << std::endl;
3740  }
3741  _out << _init << " | " << std::endl;
3742  _out << _init << " V " << std::endl;
3743  _out << _init << *_var.pExpression();
3744  _learnedBound.nextWeakerBound->print( true, _out, true );
3745  _out << std::endl;
3746  _out << std::endl;
3747  if( Settings::introduce_new_constraint_in_refinement )
3748  {
3749  _out << _init << *_var.pExpression();
3750  _learnedBound.newBound->print( true, _out, true );
3751  _out << std::endl << std::endl;
3752  }
3753  }
3754 
3755  template<class Settings, typename T1, typename T2>
3756  void Tableau<Settings,T1,T2>::print( EntryID _pivotingElement, std::ostream& _out, const std::string _init, bool _friendlyNames, bool _withZeroColumns ) const
3757  {
3758  char frameSign = '-';
3759  char separator = '|';
3760  char pivoting_separator = '#';
3761  size_t pivotingRow = 0;
3762  size_t pivotingColumn = 0;
3763  size_t basic_var_assign_width = 1;
3764  size_t basic_var_infimum_width = 1;
3765  size_t basic_var_supremum_width = 1;
3766  size_t basic_var_name_width = 1;
3767  std::vector<size_t> columnWidths;
3768  // Set the widths.
3769  for( Variable<T1,T2>* rowVar : mRows )
3770  {
3771  if( rowVar != NULL )
3772  {
3773  std::stringstream outA;
3774  if( Settings::omit_division )
3775  {
3776  if( !(rowVar->factor() == 1) )
3777  outA << rowVar->factor();
3778  }
3779  outA << rowVar->expression();
3780  size_t rowVarNameSize = outA.str().size();
3781  if( Settings::omit_division )
3782  {
3783  if( !(rowVar->factor() == 1) )
3784  rowVarNameSize += 5;
3785  }
3786  if( rowVarNameSize > basic_var_name_width )
3787  basic_var_name_width = rowVarNameSize;
3788  if( rowVar->assignment().toString().size() > basic_var_assign_width )
3789  basic_var_assign_width = rowVar->assignment().toString().size();
3790  if( rowVar->infimum().toString().size() > basic_var_infimum_width )
3791  basic_var_infimum_width = rowVar->infimum().toString().size();
3792  if( rowVar->supremum().toString().size() > basic_var_supremum_width )
3793  basic_var_supremum_width = rowVar->supremum().toString().size();
3794  }
3795  }
3796  size_t table_width = basic_var_assign_width + basic_var_infimum_width + basic_var_supremum_width + basic_var_name_width + 8;
3797  for( Variable<T1,T2>* columnVar : mColumns )
3798  {
3799  if( columnVar->size() == 0 )
3800  {
3801  columnWidths.push_back( 0 );
3802  continue;
3803  }
3804  std::stringstream ss;
3805  ss << columnVar->expression();
3806  size_t column_width = ss.str().size();
3807  if( columnVar->assignment().toString().size() > column_width )
3808  column_width = columnVar->assignment().toString().size();
3809  if( columnVar->infimum().toString().size()+2 > column_width )
3810  column_width = columnVar->infimum().toString().size()+2;
3811  if( columnVar->supremum().toString().size()+2 > column_width )
3812  column_width = columnVar->supremum().toString().size()+2;
3813  Iterator columnIter = Iterator( columnVar->startEntry(), mpEntries );
3814  while( true )
3815  {
3816  std::stringstream outA;
3817  outA << (*columnIter).content();
3818  if( outA.str().size() > column_width )
3819  column_width = outA.str().size();
3820  if( columnIter.vEnd( false ) )
3821  {
3822  break;
3823  }
3824  else
3825  {
3826  columnIter.vMove( false );
3827  }
3828  }
3829  table_width += column_width + 3;
3830  columnWidths.push_back( column_width );
3831  }
3832  // Find the row and column number of the pivoting element.
3833  if( _pivotingElement != LAST_ENTRY_ID )
3834  {
3835  pivotingRow = (*mpEntries)[_pivotingElement].rowVar()->position();
3836  pivotingColumn = (*mpEntries)[_pivotingElement].columnVar()->position();
3837  }
3838  _out << _init << std::setw( (int) table_width ) << std::setfill( frameSign ) << "" << std::endl;
3839  _out << _init << std::setw( (int) basic_var_name_width ) << std::setfill( ' ' ) << "";
3840  _out << " " << separator;
3841  for( Variable<T1,T2>* columnVar : mColumns )
3842  {
3843  if( columnVar->size() == 0 && !_withZeroColumns )
3844  continue;
3845  _out << " ";
3846  std::stringstream out;
3847  out << columnVar->expression();
3848  _out << std::setw( (int) columnWidths[columnVar->position()] ) << out.str();
3849  if( _pivotingElement != LAST_ENTRY_ID && pivotingColumn == columnVar->position() )
3850  _out << " " << pivoting_separator;
3851  else
3852  _out << " " << separator;
3853  }
3854  _out << std::endl;
3855  _out << _init << std::setw( (int) table_width ) << std::setfill( '-' ) << "" << std::endl;
3856  _out << std::setfill( ' ' );
3857  for( Variable<T1,T2>* rowVar : mRows )
3858  {
3859  size_t columnNumber = 0;
3860  _out << _init;
3861  if( rowVar != NULL )
3862  {
3863  std::stringstream out;
3864  if( Settings::omit_division )
3865  {
3866  if( !(rowVar->factor() == 1) )
3867  out << "(" << rowVar->factor() << ")*(";
3868  }
3869  out << rowVar->expression();
3870  if( Settings::omit_division )
3871  {
3872  if( !(rowVar->factor() == 1) )
3873  out << ")";
3874  }
3875  _out << std::setw( (int) basic_var_name_width ) << out.str();
3876  if( _pivotingElement != LAST_ENTRY_ID && pivotingRow == rowVar->position() )
3877  _out << " " << pivoting_separator;
3878  else
3879  _out << " " << separator;
3880  Iterator rowIter = Iterator( rowVar->startEntry(), mpEntries );
3881  size_t currentColumn = 0;
3882  while( true )
3883  {
3884  for( size_t i = currentColumn; i < (*rowIter).columnVar()->position(); ++i )
3885  {
3886  assert( columnNumber < mColumns.size() );
3887  if( mColumns[columnNumber]->size() != 0 )
3888  {
3889  _out << " ";
3890  _out << std::setw( (int) columnWidths[columnNumber] ) << "0";
3891  if( _pivotingElement != LAST_ENTRY_ID && (pivotingRow == rowVar->position() || pivotingColumn == columnNumber) )
3892  _out << " " << pivoting_separator;
3893  else
3894  _out << " " << separator;
3895  }
3896  ++columnNumber;
3897  }
3898  _out << " ";
3899  std::stringstream out;
3900  out << (*rowIter).content();
3901  _out << std::setw( (int) columnWidths[columnNumber] ) << out.str();
3902  if( _pivotingElement != LAST_ENTRY_ID && (pivotingRow == rowVar->position() || pivotingColumn == columnNumber) )
3903  _out << " " << pivoting_separator;
3904  else
3905  _out << " " << separator;
3906  ++columnNumber;
3907  currentColumn = (*rowIter).columnVar()->position()+1;
3908  if( rowIter.hEnd( false ) )
3909  {
3910  for( size_t i = currentColumn; i < mWidth; ++i )
3911  {
3912  if( mColumns[columnNumber]->size() != 0 )
3913  {
3914  _out << " ";
3915  _out << std::setw( (int) columnWidths[columnNumber] ) << "0";
3916  if( _pivotingElement != LAST_ENTRY_ID && (pivotingRow == rowVar->position() || pivotingColumn == columnNumber) )
3917  _out << " " << pivoting_separator;
3918  else
3919  _out << " " << separator;
3920  }
3921  ++columnNumber;
3922  }
3923  break;
3924  }
3925  rowIter.hMove( false );
3926  }
3927  _out << " ";
3928  _out << std::setw( (int) basic_var_assign_width ) << rowVar->assignment();
3929  _out << " [";
3930  _out << std::setw( (int) basic_var_infimum_width ) << rowVar->infimum();
3931  _out << ", ";
3932  _out << std::setw( (int) basic_var_supremum_width ) << rowVar->supremum();
3933  _out << "]";
3934  _out << std::endl;
3935  }
3936  }
3937  _out << _init << std::setw( (int) table_width ) << std::setfill( frameSign ) << "" << std::endl;
3938  _out << _init << std::setw( (int) basic_var_name_width ) << std::setfill( ' ' ) << "";
3939  _out << " " << separator;
3940  for( Variable<T1,T2>* columnVar : mColumns )
3941  {
3942  if( columnVar->size() == 0 && !_withZeroColumns )
3943  continue;
3944  _out << " ";
3945  _out << std::setw( (int) columnWidths[columnVar->position()] ) << columnVar->assignment().toString();
3946  if( _pivotingElement != LAST_ENTRY_ID && pivotingColumn == columnVar->position() )
3947  _out << " " << pivoting_separator;
3948  else
3949  _out << " " << separator;
3950  }
3951  _out << std::endl;
3952  _out << _init << std::setw( (int) basic_var_name_width ) << std::setfill( ' ' ) << "";
3953  _out << " " << separator;
3954  for( Variable<T1,T2>* columnVar : mColumns )
3955  {
3956  if( columnVar->size() == 0 && !_withZeroColumns )
3957  continue;
3958  _out << " ";
3959  std::stringstream outB;
3960  outB << "[" << columnVar->infimum().toString() << ",";
3961  _out << std::left << std::setw( (int) columnWidths[columnVar->position()] ) << outB.str();
3962  _out << std::right << "";
3963  if( _pivotingElement != LAST_ENTRY_ID && pivotingColumn == columnVar->position() )
3964  _out << " " << pivoting_separator;
3965  else
3966  _out << " " << separator;
3967  }
3968  _out << std::endl;
3969  _out << _init << std::setw( (int) basic_var_name_width ) << std::setfill( ' ' ) << "";
3970  _out << " " << separator;
3971  for( Variable<T1,T2>* columnVar : mColumns )
3972  {
3973  if( columnVar->size() == 0 && !_withZeroColumns )
3974  continue;
3975  _out << " ";
3976  std::stringstream outB;
3977  outB << " " << columnVar->supremum().toString() << "]";
3978  _out << std::setw( (int) columnWidths[columnVar->position()] ) << outB.str();
3979  if( _pivotingElement != LAST_ENTRY_ID && pivotingColumn == columnVar->position() )
3980  _out << " " << pivoting_separator;
3981  else
3982  _out << " " << separator;
3983  }
3984  _out << std::endl;
3985  _out << std::setfill( ' ' );
3986  }
3987  #endif
3988  } // end namspace lra
3989 } // end namspace smtrat