libs/cassowary/ClSimplexSolver.cc

   1 // $Id$
   2 //
   3 // Cassowary Incremental Constraint Solver
   4 // Original Smalltalk Implementation by Alan Borning
   5 // This C++ Implementation by Greg J. Badros, <gjb@cs.washington.edu>
   6 // http://www.cs.washington.edu/homes/gjb
   7 // (C) 1998, 1999 Greg J. Badros and Alan Borning
   8 // See ../LICENSE for legal details regarding this software
   9 //
  10 // ClSimplexSolver.cc
  11
  12 using namespace std;
  13
  14 #include <cassowary/debug.h>
  15 #include <cassowary/ClSimplexSolver.h>
  16 #include <cassowary/ClErrors.h>
  17 #include <cassowary/ClVariable.h>
  18 #include <cassowary/ClPoint.h>
  19 #include <cassowary/ClSlackVariable.h>
  20 #include <cassowary/ClObjectiveVariable.h>
  21 #include <cassowary/ClDummyVariable.h>
  22 #include <cassowary/cl_auto_ptr.h>
  23
  24 #include <algorithm>
  25 #include <float.h>
  26 #include <sstream>
  27 #include <queue>
  28
  29 #ifdef HAVE_CONFIG_H
  30 #include <config.h>
  31 #define CONFIG_H_INCLUDED
  32 #endif
  33
  34 // Need to delete all expressions
  35 // and all slack and dummy variables
  36 // See NewExpression -- all allocation is done in there
  37 ClSimplexSolver::~ClSimplexSolver()
  38 {
  39 #ifdef CL_SOLVER_STATS
  40   cerr << "_slackCounter == " << _slackCounter
  41        << "\n_artificialCounter == " << _artificialCounter
  42        << "\n_dummyCounter == " << _dummyCounter << endl;
  43   cerr << "stayMinusErrorVars " << _stayMinusErrorVars.size() << ", "
  44        << "stayPlusErrorVars " << _stayPlusErrorVars.size() << ", "
  45        << "errorVars " << _errorVars.size() << ", "
  46        << "markerVars " << _markerVars.size() << endl;
  47 #endif
  48   // Cannot print *this here, since local ClVariable-s may have been
  49   // destructed already
  50 }
  51
  52 // Add the constraint cn to the tableau
  53 ClSimplexSolver &
  54 ClSimplexSolver::AddConstraint(ClConstraint *const pcn)
  55 {
  56 #ifdef CL_TRACE
  57   Tracer TRACER(__FUNCTION__);
  58   cerr << "(" << *pcn << ")" << endl;
  59 #endif
  60
  61   if (!pcn->FIsOkayForSimplexSolver()) {
  62     throw ExCLTooDifficultSpecial("SimplexSolver cannot handle this constraint object");
  63   }
  64
  65   if (pcn->IsStrictInequality()) {
  66     // cannot handle strict inequalities
  67     throw ExCLStrictInequalityNotAllowed();
  68   }
  69
  70   if (pcn->ReadOnlyVars().size() > 0) {
  71     // cannot handle read-only vars
  72     throw ExCLReadOnlyNotAllowed();
  73   }
  74
  75   if (pcn->IsEditConstraint())
  76     {
  77     ClEditConstraint *pcnEdit = dynamic_cast<ClEditConstraint *>(pcn);
  78     const ClVariable &v = pcnEdit->variable();
  79     if (!v.IsExternal() ||
  80         (!FIsBasicVar(v) && !ColumnsHasKey(v)))
  81       {
  82       // we could try to make this case work,
  83       // but it'd be unnecessarily inefficient --
  84       // and probably easier for the client application
  85       // to deal with
  86       throw ExCLEditMisuse("(ExCLEditMisuse) Edit constraint on variable not in tableau.");
  87       }
  88     ClEditInfo *pcei = PEditInfoFromClv(v);
  89     if (pcei)
  90       {
  91       // we need to only add a partial _editInfoList entry for this
  92       // edit constraint since the variable is already being edited.
  93       // otherwise a more complete entry is added later in this function
  94       _editInfoList.push_back(new ClEditInfo(v, NULL, clvNil, clvNil, 0));
  95       return *this;
  96       }
  97     }
  98
  99   ClVariable clvEplus, clvEminus;
 100   Number prevEConstant;
 101   ClLinearExpression *pexpr = NewExpression(pcn, /* output to: */
 102                                             clvEplus,clvEminus,
 103                                             prevEConstant);
 104   bool fAddedOkDirectly = false;
 105
 106   try
 107     {
 108     // If possible Add expr directly to the appropriate tableau by
 109     // choosing a subject for expr (a variable to become basic) from
 110     // among the current variables in expr.  If this doesn't work use an
 111     // artificial variable.  After adding expr re-Optimize.
 112     fAddedOkDirectly = TryAddingDirectly(*pexpr);
 113     }
 114   catch (ExCLRequiredFailure &error)
 115     {
 116 #ifdef CL_TRACE
 117     cerr << "could not Add directly -- caught ExCLRequiredFailure error" << endl;
 118 #endif
 119     RemoveConstraintInternal(pcn);
 120     throw;
 121     }
 122
 123   if (!fAddedOkDirectly)
 124     { // could not Add directly
 125     ExCLRequiredFailureWithExplanation e;
 126     if (!AddWithArtificialVariable(*pexpr, e))
 127       {
 128 #ifdef CL_DEBUG_FAILURES
 129       cerr << "Failed solve! Could not Add constraint.\n"
 130            << *this << endl;
 131 #endif
 132       RemoveConstraintInternal(pcn);
 133       if (FIsExplaining())
 134         throw e;
 135       else
 136         throw ExCLRequiredFailure();
 137       }
 138     }
 139
 140   _fNeedsSolving = true;
 141
 142   if (pcn->IsEditConstraint())
 143     {
 144     ClEditConstraint *pcnEdit = dynamic_cast<ClEditConstraint *>(pcn);
 145     ClVariable clv = pcnEdit->variable();
 146     _editInfoList.push_back(new ClEditInfo(clv, pcnEdit, clvEplus, clvEminus,
 147                                        prevEConstant));
 148     }
 149
 150   if (_fAutosolve)
 151     {
 152     Optimize(_objective);
 153     SetExternalVariables();
 154     }
 155
 156   pcn->addedTo(*this);
 157   return *this;
 158 }
 159
 160 // Add weak stays to the x and y parts of each point. These have
 161 // increasing weights so that the solver will try to satisfy the x
 162 // and y stays on the same point, rather than the x stay on one and
 163 // the y stay on another.
 164 ClSimplexSolver &
 165 ClSimplexSolver::AddPointStays(const vector<const ClPoint *> &listOfPoints,
 166                                const ClStrength &strength)
 167 {
 168 #ifdef CL_TRACE
 169   Tracer TRACER(__FUNCTION__);
 170 #endif
 171
 172   vector<const ClPoint *>::const_iterator it = listOfPoints.begin();
 173   double weight = 1.0;
 174   static const double multiplier = 2.0;
 175   for ( ; it != listOfPoints.end(); ++it )
 176     {
 177     AddPointStay((*it)->X(),(*it)->Y(),strength,weight);
 178     weight *= multiplier;
 179     }
 180   return *this;
 181 }
 182
 183 ClSimplexSolver &
 184 ClSimplexSolver::AddPointStay(const ClPoint &clp, const ClStrength &strength, double weight)
 185 {
 186   AddPointStay(clp.X(),clp.Y(),strength,weight);
 187   return *this;
 188 }
 189
 190
 191 ClSimplexSolver &
 192 ClSimplexSolver::RemoveEditVarsTo(unsigned int n)
 193 {
 194   queue<ClVariable> qclv;
 195   ClVarSet sclvStillEditing; // Set of edit variables that we need to *not* remove
 196 #ifdef DEBUG_NESTED_EDITS
 197   cerr << __FUNCTION__ << " " << n << endl;
 198 #endif
 199   unsigned int i = 0;
 200   for ( ClEditInfoList::const_iterator it = _editInfoList.begin();
 201         (it != _editInfoList.end() && _editInfoList.size() != static_cast<unsigned int>(n));
 202         ++it, ++i )
 203     {
 204     const ClEditInfo *pcei = (*it);
 205     assert(pcei);
 206 #ifdef DEBUG_NESTED_EDITS
 207     cerr << __FUNCTION__ << "Checking " << pcei->_clv
 208          << ", index = " << i << endl;
 209 #endif
 210     if (i >= n)
 211       qclv.push(pcei->_clv);
 212     else
 213       sclvStillEditing.insert(pcei->_clv);
 214     }
 215   while (!qclv.empty())
 216     {
 217     ClVariable clv = qclv.front();
 218     // only remove the variable if it's not in the set of variable
 219     // from a previous nested outer edit
 220     // e.g., if I do:
 221     // Edit x,y
 222     // Edit w,h,x,y
 223     // EndEdit
 224     // The end edit needs to only get rid of the edits on w,h
 225     // not the ones on x,y
 226     if (sclvStillEditing.find(clv) == sclvStillEditing.end())
 227       {
 228 #ifdef DEBUG_NESTED_EDITS
 229       cerr << __FUNCTION__ << ": Removing " << clv << endl;
 230 #endif
 231       RemoveEditVar(clv);
 232       }
 233 #ifdef DEBUG_NESTED_EDITS
 234     else
 235       {
 236       cerr << __FUNCTION__ << ": Not removing " << clv << endl;
 237       }
 238 #endif
 239     qclv.pop();
 240     }
 241   while (_editInfoList.size() > n) {
 242     _editInfoList.pop_back();
 243   }
 244
 245   return *this;
 246 }
 247
 248
 249 /* A predicate used for remove_if */
 250 class VarInVarSet : public unary_function<ClVariable,bool> {
 251 public:
 252   VarInVarSet(const ClVarSet &clvset) :
 253       _set(clvset),
 254       _setEnd(clvset.end())
 255     { }
 256
 257   bool operator ()(ClVariable clv) const {
 258     return (_set.find(clv) != _setEnd);
 259   }
 260
 261 private:
 262   const ClVarSet &_set;
 263   const ClVarSet::iterator _setEnd;
 264 };
 265
 266
 267
 268 // Remove the constraint cn from the tableau
 269 // Also remove any error variable associated with cn
 270 ClSimplexSolver &
 271 ClSimplexSolver::RemoveConstraintInternal(const ClConstraint *const pcn)
 272 {
 273 #ifdef CL_TRACE
 274   Tracer TRACER(__FUNCTION__);
 275   cerr << "(" << *pcn << ")" << endl;
 276 #endif
 277
 278   // We are about to remove a constraint.  There may be some stay
 279   // constraints that were unsatisfied previously -- if we just
 280   // removed the constraint these could come into play.  Instead,
 281   // Reset all of the stays so that things should stay where they are
 282   // at the moment.
 283   _fNeedsSolving = true;
 284
 285   ResetStayConstants();
 286
 287   // remove any error variables from the objective function
 288   ClLinearExpression *pzRow = RowExpression(_objective);
 289
 290 #ifdef CL_TRACE
 291   cerr << _errorVars << endl << endl;
 292 #endif
 293
 294   ClConstraintToVarSetMap::iterator
 295     it_eVars = _errorVars.find(pcn);
 296   bool fFoundErrorVar = (it_eVars != _errorVars.end());
 297
 298   if (fFoundErrorVar)
 299     {
 300     ClVarSet &eVars = (*it_eVars).second;
 301     ClVarSet::iterator it = eVars.begin();
 302     for ( ; it != eVars.end(); ++it )
 303       {
 304       const ClLinearExpression *pexpr = RowExpression(*it);
 305       if (pexpr == NULL )
 306         {
 307         pzRow->AddVariable(*it,-pcn->weight() * pcn->strength().symbolicWeight().AsDouble(),
 308                            _objective,*this);
 309         }
 310       else
 311         { // the error variable was in the basis
 312         pzRow->AddExpression(*pexpr,-pcn->weight() * pcn->strength().symbolicWeight().AsDouble(),
 313                              _objective,*this);
 314         }
 315       }
 316     }
 317
 318   ClConstraintToVarMap::iterator
 319     it_marker = _markerVars.find(pcn);
 320   if (it_marker == _markerVars.end())
 321     { // could not find the constraint
 322     throw ExCLConstraintNotFound();
 323     }
 324   // try to make the marker variable basic if it isn't already
 325   const ClVariable marker = (*it_marker).second;
 326   _markerVars.erase(it_marker);
 327   _constraintsMarked.erase(marker);
 328 #ifdef CL_TRACE
 329   cerr << "Looking to remove var " << marker << endl;
 330 #endif
 331   if (!FIsBasicVar(marker))
 332     { // not in the basis, so need to do some work
 333     // first choose which variable to move out of the basis
 334     // only consider restricted basic variables
 335     ClVarSet &col = _columns[marker];
 336     ClVarSet::iterator it_col = col.begin();
 337 #ifdef CL_TRACE
 338     cerr << "Must Pivot -- columns are " << col << endl;
 339 #endif
 340
 341     ClVariable exitVar = clvNil;
 342     bool fExitVarSet = false;
 343     double minRatio = 0.0;
 344     for ( ; it_col != col.end(); ++it_col)
 345       {
 346       const ClVariable &v = *it_col;
 347       if (v.IsRestricted() )
 348         {
 349         const ClLinearExpression *pexpr = RowExpression(v);
 350         assert(pexpr != NULL );
 351         Number coeff = pexpr->CoefficientFor(marker);
 352 #ifdef CL_TRACE
 353         cerr << "Marker " << marker << "'s coefficient in " << *pexpr << " is "
 354              << coeff << endl;
 355 #endif
 356         // only consider negative coefficients
 357         if (coeff < 0.0)
 358           {
 359           Number r = - pexpr->Constant() / coeff;
 360           if (!fExitVarSet || r < minRatio)
 361             {
 362             minRatio = r;
 363             exitVar = v;
 364             fExitVarSet = true;
 365             }
 366           }
 367         }
 368       }
 369     // if we didn't set exitvar above, then either the marker
 370     // variable has a positive coefficient in all equations, or it
 371     // only occurs in equations for unrestricted variables.  If it
 372     // does occur in an equation for a restricted variable, pick the
 373     // equation that gives the smallest ratio.  (The row with the
 374     // marker variable will become infeasible, but all the other rows
 375     // will still be feasible; and we will be dropping the row with
 376     // the marker variable.  In effect we are removing the
 377     // non-negativity restriction on the marker variable.)
 378     if (!fExitVarSet)
 379       {
 380 #ifdef CL_TRACE
 381       cerr << "exitVar did not get set" << endl;
 382 #endif
 383       it_col = col.begin();
 384       for ( ; it_col != col.end(); ++it_col)
 385         {
 386         ClVariable v = *it_col;
 387         if (v.IsRestricted() )
 388           {
 389           const ClLinearExpression *pexpr = RowExpression(v);
 390           assert(pexpr != NULL);
 391           Number coeff = pexpr->CoefficientFor(marker);
 392           Number r = pexpr->Constant() / coeff;
 393           if (!fExitVarSet || r < minRatio)
 394             {
 395             minRatio = r;
 396             exitVar = v;
 397             fExitVarSet = true;
 398             }
 399           }
 400         }
 401       }
 402
 403     if (!fExitVarSet)
 404       { // exitVar is still nil
 405       // If col is empty, then exitVar doesn't occur in any equations,
 406       // so just remove it.  Otherwise pick an exit var from among the
 407       // unrestricted variables whose equation involves the marker var
 408       if (col.size() == 0)
 409         {
 410         RemoveColumn(marker);
 411         }
 412       else
 413         {
 414           // A. Beurive' Tue Sep 14 18:26:05 CEST 1999
 415           // Don't pick the objective, or it will be removed!
 416           it_col = col.begin();
 417           for ( ; it_col != col.end(); ++it_col)
 418             {
 419               ClVariable v = *it_col;
 420               if (v != _objective)
 421                 {
 422                   exitVar = v;
 423                   fExitVarSet = true;
 424                   break;
 425                 }
 426             }
 427           assert(fExitVarSet == true);
 428         }
 429       }
 430
 431     if (fExitVarSet)
 432       {
 433       Pivot(marker,exitVar);
 434       }
 435     }
 436
 437   if (FIsBasicVar(marker))
 438     {
 439     ClLinearExpression *pexpr = RemoveRow(marker);
 440 #ifdef CL_TRACE
 441     cerr << "delete@ " << pexpr << endl;
 442 #endif
 443     delete pexpr;
 444     }
 445
 446   // Delete any error variables.  If cn is an inequality, it also
 447   // contains a slack variable; but we use that as the marker variable
 448   // and so it has been deleted when we removed its row.
 449   if (fFoundErrorVar)
 450     {
 451     ClVarSet &eVars = (*it_eVars).second;
 452     ClVarSet::iterator it = eVars.begin();
 453     for ( ; it != eVars.end(); ++it )
 454       {
 455       ClVariable v = (*it);
 456       if (v != marker)
 457         {
 458         RemoveColumn(v);
 459         }
 460       }
 461     }
 462
 463   if (pcn->isStayConstraint())
 464     {
 465     // iterate over the stay{Plus,Minus}ErrorVars and remove those
 466     // variables v in those vectors that are also in set eVars
 467     if (fFoundErrorVar)
 468       {
 469       ClVarSet &eVars = (*it_eVars).second;
 470       _stayPlusErrorVars
 471         .erase(remove_if(_stayPlusErrorVars.begin(),_stayPlusErrorVars.end(),
 472                          VarInVarSet(eVars)),
 473                _stayPlusErrorVars.end());
 474       _stayMinusErrorVars
 475         .erase(remove_if(_stayMinusErrorVars.begin(),_stayMinusErrorVars.end(),
 476                          VarInVarSet(eVars)),
 477                _stayMinusErrorVars.end());
 478       }
 479     }
 480   else if (pcn->IsEditConstraint())
 481     {
 482     const ClEditConstraint *pcnEdit = dynamic_cast<const ClEditConstraint *>(pcn);
 483     const ClVariable clv = pcnEdit->variable();
 484     ClEditInfo *pcei = PEditInfoFromClv(clv);
 485     assert(pcei);
 486     ClVariable clvEditMinus = pcei->_clvEditMinus;
 487     RemoveColumn(clvEditMinus);  // clvEditPlus is a marker var and gets removed later
 488     delete pcei;
 489     _editInfoList.remove(pcei);
 490     }
 491
 492   if (fFoundErrorVar)
 493     {
 494     // This code is not needed since the variables are deleted
 495     // when they are removed from the row --
 496     // leaving it in results in double deletions
 497     // delete the constraint's error variables
 498     //    ClVarSet &evars_set = (*it_eVars).second;
 499     //    ClVarSet::const_iterator it_set = evars_set.begin();
 500     //    for ( ; it_set != evars_set.end(); ++it_set)
 501     //      {
 502     //      delete *it_set;
 503     //      }
 504     _errorVars.erase((*it_eVars).first);
 505     }
 506
 507   if (_fAutosolve)
 508     {
 509     Optimize(_objective);
 510     SetExternalVariables();
 511     }
 512
 513   return *this;
 514 }
 515
 516
 517 // Re-initialize this solver from the original constraints, thus
 518 // getting rid of any accumulated numerical problems.  (Actually,
 519 // Alan hasn't observed any such problems yet, but here's the method
 520 // anyway.)
 521 void
 522 ClSimplexSolver::Reset()
 523 {
 524 #ifdef CL_TRACE
 525   Tracer TRACER(__FUNCTION__);
 526   cerr << "()" << endl;
 527 #endif
 528   // FIXGJB  -- can postpone writing this for a while
 529   // gotta be careful, though, as it's a likely place for
 530   // a memory leak to sneak in
 531   assert(false);
 532 }
 533
 534
 535 // Re-solve the cuurent collection of constraints, given the new
 536 // values for the edit variables that have already been
 537 // suggested (see SuggestValue() method)
 538 void
 539 ClSimplexSolver::Resolve()
 540 { // CODE DUPLICATED ABOVE
 541 #ifdef CL_TRACE
 542   Tracer TRACER(__FUNCTION__);
 543 #endif
 544   DualOptimize();
 545   SetExternalVariables();
 546   _infeasibleRows.clear();
 547   if (_fResetStayConstantsAutomatically)
 548     ResetStayConstants();
 549 }
 550
 551 ClSimplexSolver &
 552 ClSimplexSolver::SuggestValue(ClVariable v, Number x)
 553 {
 554 #ifdef CL_TRACE
 555   Tracer TRACER(__FUNCTION__);
 556 #endif
 557   ClEditInfo *pcei = PEditInfoFromClv(v);
 558   if (NULL == pcei)
 559     {
 560 #ifndef CL_NO_IO
 561         std::stringstream ss;
 562     ss << "SuggestValue for variable " << v << ", but var is not an edit variable" << ends;
 563     throw ExCLEditMisuse(ss.str().c_str());
 564 #else
 565     throw ExCLEditMisuse(v.Name().c_str());
 566 #endif
 567     }
 568   ClVariable clvEditPlus = pcei->_clvEditPlus;
 569   ClVariable clvEditMinus = pcei->_clvEditMinus;
 570   Number delta = x - pcei->_prevEditConstant;
 571   pcei->_prevEditConstant = x;
 572   DeltaEditConstant(delta,clvEditPlus,clvEditMinus);
 573   return *this;
 574 }
 575
 576 // Re-solve the cuurent collection of constraints, given the new
 577 // values for the edit variables that have already been
 578 // suggested (see SuggestValue() method)
 579 // This is not guaranteed to work if you remove an edit constraint
 580 // from the middle of the edit constraints you added
 581 // (e.g., edit A, edit B, edit C, remove B -> this will fail!)
 582 // DEPRECATED
 583 void
 584 ClSimplexSolver::Resolve(const vector<Number> &newEditConstants)
 585 {
 586   ClEditInfoList::iterator it = _editInfoList.begin();
 587   unsigned int i = 0;
 588   for (; i < newEditConstants.size() && it != _editInfoList.end(); ++it, ++i)
 589     {
 590     ClEditInfo *pcei = (*it);
 591     SuggestValue(pcei->_clv,newEditConstants[i]);
 592     }
 593   Resolve();
 594 }
 595
 596
 597 //// protected
 598
 599 // Add the constraint expr=0 to the inequality tableau using an
 600 // artificial variable.  To do this, create an artificial variable
 601 // av and Add av=expr to the inequality tableau, then make av be 0.
 602 // (Raise an exception if we can't attain av=0 -- and prepare explanation)
 603 bool
 604 ClSimplexSolver::AddWithArtificialVariable(ClLinearExpression &expr,
 605                                            ExCLRequiredFailureWithExplanation &e)
 606 {
 607 #ifdef CL_TRACE
 608   Tracer TRACER(__FUNCTION__);
 609   cerr << "(" << expr << ")" << endl;
 610 #endif
 611
 612   // Allocate the objects on the heap because the objects
 613   // will remain in the tableau if we throw an exception,
 614   // and that will result in the destructor cleaning up
 615   // after us
 616   ClSlackVariable *pav = new ClSlackVariable(++_artificialCounter,"a");
 617   ClObjectiveVariable *paz = new ClObjectiveVariable("az");
 618   ClLinearExpression *pazRow = new ClLinearExpression(expr);
 619   // the artificial objective is av, which we know is equal to expr
 620   // (which contains only parametric variables)
 621
 622 #ifdef CL_FIND_LEAK
 623   cerr << "aC = " << _artificialCounter
 624        << "\nDeletes = " << _cArtificialVarsDeleted << endl;
 625 #endif
 626 #ifdef CL_TRACE
 627   cerr << __FUNCTION__ << " before addRow-s:\n"
 628        << (*this) << endl;
 629 #endif
 630
 631   // the artificial objective is av, which we know is equal to expr
 632   // (which contains only parametric variables)
 633
 634   // objective is treated as a row in the tableau,
 635   // so do the substitution for its value (we are minimizing
 636   // the artificial variable)
 637   // this row will be removed from the tableau after optimizing
 638   addRow(*paz,*pazRow);
 639
 640   // now Add the normal row to the tableau -- when artifical
 641   // variable is minimized to 0 (if possible)
 642   // this row remains in the tableau to maintain the constraint
 643   // we are trying to Add
 644   addRow(*pav,expr);
 645
 646 #ifdef CL_TRACE
 647   cerr << __FUNCTION__ << " after addRow-s:\n"
 648        << (*this) << endl;
 649 #endif
 650
 651   // try to Optimize az to 0
 652   // note we are *not* optimizing the real objective, but optimizing
 653   // the artificial objective to see if the error in the constraint
 654   // we are adding can be set to 0
 655   Optimize(*paz);
 656
 657   // Careful, we want to get the Expression that is in
 658   // the tableau, not the one we initialized it with!
 659   ClLinearExpression *pazTableauRow = RowExpression(*paz);
 660 #ifdef CL_TRACE
 661   cerr << "pazTableauRow->Constant() == " << pazTableauRow->Constant() << endl;
 662 #endif
 663
 664   // Check that we were able to make the objective value 0
 665   // If not, the original constraint was not satisfiable
 666   if (!ClApprox(pazTableauRow->Constant(),0.0))
 667     {
 668     BuildExplanation(e, paz, pazTableauRow);
 669     // remove the artificial objective row that we just
 670     // added temporarily
 671     delete RemoveRow(*paz);
 672     // and delete the artificial objective variable that we also added above
 673     delete paz;
 674     return false;
 675     }
 676
 677   // see if av is a basic variable
 678   const ClLinearExpression *pe = RowExpression(*pav);
 679   if (pe != NULL)
 680     {
 681     // FIXGJB: do we ever even get here?
 682     // Find another variable in this row and Pivot, so that av becomes parametric
 683     // If there isn't another variable in the row then
 684     // the tableau contains the equation av = 0  -- just delete av's row
 685     if (pe->IsConstant())
 686       {
 687       // FIXGJB: do we ever get here?
 688       assert(ClApprox(pe->Constant(),0.0));
 689       delete RemoveRow(*pav);
 690       // remove the temporary objective function
 691       // FIXGJB may need this too: delete RemoveRow(*paz);
 692       delete pav;
 693 #ifdef CL_FIND_LEAK
 694       ++_cArtificialVarsDeleted;
 695 #endif
 696       return true;
 697       }
 698     ClVariable entryVar = pe->AnyPivotableVariable();
 699     if (entryVar.IsNil())
 700       {
 701       BuildExplanation(e, *pav, pe);
 702       return false; /* required failure */
 703       }
 704     Pivot(entryVar, *pav);
 705     }
 706   // now av should be parametric
 707   assert(RowExpression(*pav) == NULL);
 708   RemoveColumn(*pav);
 709   delete pav;
 710 #ifdef CL_FIND_LEAK
 711   ++_cArtificialVarsDeleted;
 712 #endif
 713   // remove the temporary objective function
 714   delete RemoveRow(*paz);
 715   delete paz;
 716   return true;
 717 }
 718
 719
 720 // Using the given equation (av = cle) build an explanation which
 721 // implicates all constraints used to construct the equation. That
 722 // is, everything for which the variables in the equation are markers.
 723 void ClSimplexSolver::BuildExplanation(ExCLRequiredFailureWithExplanation &e,
 724                                        ClVariable av,
 725                                        const ClLinearExpression *pcle)
 726 {
 727   ClVarToConstraintMap::iterator it_cn;
 728   it_cn = _constraintsMarked.find(av);
 729   if (it_cn != _constraintsMarked.end())
 730     {
 731       e.AddConstraint((*it_cn).second);
 732     }
 733
 734   assert(pcle != NULL);
 735
 736   const ClVarToNumberMap & terms = pcle->Terms();
 737   ClVarToNumberMap::const_iterator it_term;
 738   for (it_term = terms.begin(); it_term != terms.end(); it_term++)
 739     {
 740     it_cn = _constraintsMarked.find((*it_term).first);
 741     if (it_cn != _constraintsMarked.end())
 742       {
 743       e.AddConstraint((*it_cn).second);
 744       }
 745     }
 746 }
 747
 748
 749
 750 // We are trying to Add the constraint expr=0 to the appropriate
 751 // tableau.  Try to Add expr directly to the tableaus without
 752 // creating an artificial variable.  Return true if successful and
 753 // false if not.
 754 bool
 755 ClSimplexSolver::TryAddingDirectly(ClLinearExpression &expr)
 756 {
 757 #ifdef CL_TRACE
 758   Tracer TRACER(__FUNCTION__);
 759   cerr << "(" << expr << ")" << endl;
 760 #endif
 761   ClVariable subject = ChooseSubject(expr);
 762   if (subject.get_pclv() == NULL )
 763     {
 764 #ifdef CL_TRACE
 765     cerr << "- returning false" << endl;
 766 #endif
 767     return false;
 768     }
 769   expr.NewSubject(subject);
 770   if (ColumnsHasKey(subject))
 771     {
 772     SubstituteOut(subject,expr);
 773     }
 774   addRow(subject,expr);
 775 #ifdef CL_TRACE
 776   cerr << "- returning true" << endl;
 777 #endif
 778   return true; // successfully added directly
 779 }
 780
 781
 782 // We are trying to Add the constraint expr=0 to the tableaux.  Try
 783 // to choose a subject (a variable to become basic) from among the
 784 // current variables in expr.  If expr contains any unrestricted
 785 // variables, then we must choose an unrestricted variable as the
 786 // subject.  Also, if the subject is new to the solver we won't have
 787 // to do any substitutions, so we prefer new variables to ones that
 788 // are currently noted as parametric.  If expr contains only
 789 // restricted variables, if there is a restricted variable with a
 790 // negative coefficient that is new to the solver we can make that
 791 // the subject.  Otherwise we can't find a subject, so return nil.
 792 // (In this last case we have to Add an artificial variable and use
 793 // that variable as the subject -- this is done outside this method
 794 // though.)
 795 //
 796 // Note: in checking for variables that are new to the solver, we
 797 // ignore whether a variable occurs in the objective function, since
 798 // new slack variables are added to the objective function by
 799 // 'NewExpression:', which is called before this method.
 800 ClVariable
 801 ClSimplexSolver::ChooseSubject(ClLinearExpression &expr)
 802 {
 803 #ifdef CL_TRACE
 804   Tracer TRACER(__FUNCTION__);
 805   cerr << "(" << expr << ")" << endl;
 806 #endif
 807   ClVariable subject(clvNil); // the current best subject, if any
 808
 809   // true iff we have found a subject that is an unrestricted variable
 810   bool foundUnrestricted = false;
 811
 812   // true iff we have found a restricted variable that is new to the
 813   // solver (except for being in the obj. function) and that has a
 814   // negative coefficient
 815   bool foundNewRestricted = false;
 816
 817   const ClVarToNumberMap &terms = expr.Terms();
 818   ClVarToNumberMap::const_iterator it = terms.begin();
 819   for ( ; it != terms.end(); ++it )
 820     {
 821     ClVariable v = (*it).first;
 822     Number c = (*it).second;
 823
 824     if (foundUnrestricted)
 825       {
 826       // We have already found an unrestricted variable.  The only
 827       // time we will want to use v instead of the current choice
 828       // 'subject' is if v is unrestricted and new to the solver and
 829       // 'subject' isn't new.  If this is the case just pick v
 830       // immediately and return.
 831       if (!v.IsRestricted())
 832         {
 833         if (!ColumnsHasKey(v))
 834           return v;
 835         }
 836       }
 837     else
 838       { // we haven't found an restricted variable yet
 839       if (v.IsRestricted())
 840         {
 841         // v is restricted.  If we have already found a suitable
 842         // restricted variable just stick with that.  Otherwise, if v
 843         // is new to the solver and has a negative coefficient pick
 844         // it.  Regarding being new to the solver -- if the variable
 845         // occurs only in the objective function we regard it as being
 846         // new to the solver, since error variables are added to the
 847         // objective function when we make the Expression.  We also
 848         // never pick a dummy variable here.
 849         if (!foundNewRestricted && !v.IsDummy() && c < 0.0)
 850           {
 851           const ClTableauColumnsMap &col = Columns();
 852           ClTableauColumnsMap::const_iterator it_col = col.find(v);
 853           if (it_col == col.end() ||
 854               ( col.size() == 1 && ColumnsHasKey(_objective) ) )
 855             {
 856             subject = v;
 857             foundNewRestricted = true;
 858             }
 859           }
 860         }
 861       else
 862         {
 863         // v is unrestricted.
 864         // If v is also new to the solver just pick it now
 865         subject = v;
 866         foundUnrestricted = true;
 867         }
 868       }
 869     }
 870   if (!subject.IsNil())
 871     return subject;
 872
 873   // subject is nil.
 874   // Make one last check -- if all of the variables in expr are dummy
 875   // variables, then we can pick a dummy variable as the subject
 876   Number coeff = 0;
 877   it = terms.begin();
 878   for ( ; it != terms.end(); ++it )
 879     {
 880     ClVariable v = (*it).first;
 881     Number c = (*it).second;
 882     if (!v.IsDummy())
 883       return clvNil; // nope, no luck
 884     // if v is new to the solver, tentatively make it the subject
 885     if (!ColumnsHasKey(v))
 886       {
 887       subject = v;
 888       coeff = c;
 889       }
 890     }
 891
 892   // If we get this far, all of the variables in the Expression should
 893   // be dummy variables.  If the Constant is nonzero we are trying to
 894   // Add an unsatisfiable required constraint.  (Remember that dummy
 895   // variables must take on a value of 0.)  Otherwise, if the Constant
 896   // is Zero, multiply by -1 if necessary to make the coefficient for
 897   // the subject negative."
 898   if (!ClApprox(expr.Constant(),0.0))
 899     {
 900 #ifdef CL_DEBUG_FAILURES
 901     cerr << "required failure in choose subject:\n"
 902          << *this << endl;
 903 #endif
 904     if (FIsExplaining())
 905       {
 906       ExCLRequiredFailureWithExplanation e;
 907       BuildExplanation(e, clvNil, &expr);
 908       throw e;
 909       }
 910     else
 911       throw ExCLRequiredFailure();
 912     }
 913   if (coeff > 0.0)
 914     {
 915     expr.MultiplyMe(-1);
 916     }
 917   return subject;
 918 }
 919
 920 // Each of the non-required edits will be represented by an equation
 921 // of the form
 922 //    v = c + eplus - eminus
 923 // where v is the variable with the edit, c is the previous edit
 924 // value, and eplus and eminus are slack variables that hold the
 925 // error in satisfying the edit constraint.  We are about to change
 926 // something, and we want to fix the constants in the equations
 927 // representing the edit constraints.  If one of eplus and eminus is
 928 // basic, the other must occur only in the Expression for that basic
 929 // error variable.  (They can't both be basic.)  Fix the Constant in
 930 // this Expression.  Otherwise they are both nonbasic.  Find all of
 931 // the expressions in which they occur, and fix the constants in
 932 // those.  See the UIST paper for details.
 933 // (This comment was for resetEditConstants(), but that is now
 934 // gone since it was part of the screwey vector-based interface
 935 // to resolveing. --02/15/99 gjb)
 936 void
 937 ClSimplexSolver::DeltaEditConstant(Number delta,
 938                                    ClVariable plusErrorVar,
 939                                    ClVariable minusErrorVar)
 940 {
 941 #ifdef CL_TRACE
 942   Tracer TRACER(__FUNCTION__);
 943   cerr << "(" << delta << ", " << plusErrorVar << ", " << minusErrorVar << ")" << endl;
 944 #endif
 945   // first check if the plusErrorVar is basic
 946   ClLinearExpression *pexprPlus = RowExpression(plusErrorVar);
 947   if (pexprPlus != NULL )
 948     {
 949     pexprPlus->IncrementConstant(delta);
 950     // error variables are always restricted
 951     // so the row is infeasible if the Constant is negative
 952     if (pexprPlus->Constant() < 0.0)
 953       {
 954       _infeasibleRows.insert(plusErrorVar);
 955       }
 956     return;
 957     }
 958   // check if minusErrorVar is basic
 959   ClLinearExpression *pexprMinus = RowExpression(minusErrorVar);
 960   if (pexprMinus != NULL)
 961     {
 962     pexprMinus->IncrementConstant(-delta);
 963     if (pexprMinus->Constant() < 0.0)
 964       {
 965       _infeasibleRows.insert(minusErrorVar);
 966       }
 967     return;
 968     }
 969   // Neither is basic.  So they must both be nonbasic, and will both
 970   // occur in exactly the same expressions.  Find all the expressions
 971   // in which they occur by finding the column for the minusErrorVar
 972   // (it doesn't matter whether we look for that one or for
 973   // plusErrorVar).  Fix the constants in these expressions.
 974   ClVarSet &columnVars = _columns[minusErrorVar];
 975   ClVarSet::iterator it = columnVars.begin();
 976   for (; it != columnVars.end(); ++it)
 977     {
 978     ClVariable basicVar = *it;
 979     ClLinearExpression *pexpr = RowExpression(basicVar);
 980     assert(pexpr != NULL );
 981     double c = pexpr->CoefficientFor(minusErrorVar);
 982     pexpr->IncrementConstant(c*delta);
 983     if (basicVar.IsRestricted() && pexpr->Constant() < 0.0)
 984       {
 985       _infeasibleRows.insert(basicVar);
 986       }
 987     }
 988 }
 989
 990 // We have set new values for the constants in the edit constraints.
 991 // Re-Optimize using the dual simplex algorithm.
 992 void
 993 ClSimplexSolver::DualOptimize()
 994 {
 995 #ifdef CL_TRACE
 996   Tracer TRACER(__FUNCTION__);
 997   cerr << "()" << endl;
 998 #endif
 999   const ClLinearExpression *pzRow = RowExpression(_objective);
1000   // need to handle infeasible rows
1001   while (!_infeasibleRows.empty())
1002     {
1003     ClVarSet::iterator it_exitVar = _infeasibleRows.begin();
1004     ClVariable exitVar = *it_exitVar;
1005     _infeasibleRows.erase(it_exitVar);
1006     ClVariable entryVar;
1007     // exitVar might have become basic after some other pivoting
1008     // so allow for the case of its not being there any longer
1009     ClLinearExpression *pexpr = RowExpression(exitVar);
1010     if (pexpr != NULL )
1011       {
1012       // make sure the row is still not feasible
1013       if (pexpr->Constant() < 0.0)
1014         {
1015         double ratio = DBL_MAX;
1016         double r;
1017         ClVarToNumberMap &terms = pexpr->Terms();
1018         ClVarToNumberMap::iterator it = terms.begin();
1019         for ( ; it != terms.end(); ++it )
1020           {
1021           ClVariable v = (*it).first;
1022           Number c = (*it).second;
1023           if (c > 0.0 && v.IsPivotable())
1024             {
1025             Number zc = pzRow->CoefficientFor(v);
1026             r = zc/c; // FIXGJB r:= zc/c or Zero, as ClSymbolicWeight-s
1027             if (r < ratio)
1028               {
1029               entryVar = v;
1030               ratio = r;
1031               }
1032             }
1033           }
1034         if (ratio == DBL_MAX)
1035           {
1036           stringstream ss;
1037           ss << "ratio == nil (DBL_MAX)" << ends;
1038           throw ExCLInternalError(ss.str().c_str());
1039           }
1040         Pivot(entryVar,exitVar);
1041         }
1042       }
1043     }
1044 }
1045
1046 // Make a new linear Expression representing the constraint cn,
1047 // replacing any basic variables with their defining expressions.
1048 // Normalize if necessary so that the Constant is non-negative.  If
1049 // the constraint is non-required give its error variables an
1050 // appropriate weight in the objective function.
1051 ClLinearExpression *
1052 ClSimplexSolver::NewExpression(const ClConstraint *pcn,
1053                                /* output to */
1054                                ClVariable &clvEplus,
1055                                ClVariable &clvEminus,
1056                                Number &prevEConstant)
1057 {
1058 #ifdef CL_TRACE
1059   Tracer TRACER(__FUNCTION__);
1060   cerr << "(" << *pcn << ")" << endl;
1061   cerr << "cn.IsInequality() == " << pcn->IsInequality() << endl;
1062   cerr << "cn.IsRequired() == " << pcn->IsRequired() << endl;
1063 #endif
1064   const ClLinearExpression &cnExpr = pcn->Expression();
1065   cl_auto_ptr<ClLinearExpression> pexpr ( new ClLinearExpression(cnExpr.Constant()) );
1066   cl_auto_ptr<ClSlackVariable> pslackVar;
1067   cl_auto_ptr<ClDummyVariable> pdummyVar;
1068   cl_auto_ptr<ClSlackVariable> peminus(0);
1069   cl_auto_ptr<ClSlackVariable> peplus(0);
1070   const ClVarToNumberMap &cnTerms = cnExpr.Terms();
1071   ClVarToNumberMap::const_iterator it = cnTerms.begin();
1072   for ( ; it != cnTerms.end(); ++it)
1073     {
1074     ClVariable v = (*it).first;
1075     Number c = (*it).second;
1076     const ClLinearExpression *pe = RowExpression(v);
1077     if (pe == NULL)
1078       {
1079       pexpr->AddVariable(v,c);
1080       }
1081     else
1082       {
1083       pexpr->AddExpression(*pe,c);
1084       }
1085     }
1086
1087   // Add slack and error variables as needed
1088   if (pcn->IsInequality())
1089     {
1090     // cn is an inequality, so Add a slack variable.  The original
1091     // constraint is expr>=0, so that the resulting equality is
1092     // expr-slackVar=0.  If cn is also non-required Add a negative
1093     // error variable, giving
1094     //    expr-slackVar = -errorVar, in other words
1095     //    expr-slackVar+errorVar=0.
1096     // Since both of these variables are newly created we can just Add
1097     // them to the Expression (they can't be basic).
1098     ++_slackCounter;
1099     ReinitializeAutoPtr(pslackVar,new ClSlackVariable (_slackCounter, "s"));
1100     pexpr->setVariable(*pslackVar,-1);
1101     // index the constraint under its slack variable and vice-versa
1102     _markerVars[pcn] = pslackVar.get();
1103     _constraintsMarked[pslackVar.get()] = pcn;
1104
1105     if (!pcn->IsRequired())
1106       {
1107       ++_slackCounter;
1108       ReinitializeAutoPtr(peminus,new ClSlackVariable (_slackCounter, "em"));
1109       pexpr->setVariable(*peminus,1.0);
1110       // Add emnius to the objective function with the appropriate weight
1111       ClLinearExpression *pzRow = RowExpression(_objective);
1112       // FIXGJB: pzRow->AddVariable(eminus,pcn->strength().symbolicWeight() * pcn->weight());
1113       ClSymbolicWeight sw = pcn->strength().symbolicWeight().Times(pcn->weight());
1114       pzRow->setVariable(*peminus,sw.AsDouble());
1115       _errorVars[pcn].insert(peminus.get());
1116       NoteAddedVariable(*peminus,_objective);
1117       }
1118     }
1119   else
1120     { // cn is an equality
1121     if (pcn->IsRequired())
1122       {
1123       // Add a dummy variable to the Expression to serve as a marker
1124       // for this constraint.  The dummy variable is never allowed to
1125       // enter the basis when pivoting.
1126       ++_dummyCounter;
1127       ReinitializeAutoPtr(pdummyVar,new ClDummyVariable (_dummyCounter, "d"));
1128       pexpr->setVariable(*pdummyVar,1.0);
1129       _markerVars[pcn] = pdummyVar.get();
1130       _constraintsMarked[pdummyVar.get()] = pcn;
1131 #ifdef CL_TRACE
1132       cerr << "Adding dummyVar == d" << _dummyCounter << endl;
1133 #endif
1134       }
1135     else
1136       {
1137       // cn is a non-required equality.  Add a positive and a negative
1138       // error variable, making the resulting constraint
1139       //       expr = eplus - eminus,
1140       // in other words:  expr-eplus+eminus=0
1141       ++_slackCounter;
1142       ReinitializeAutoPtr(peplus,new ClSlackVariable (_slackCounter, "ep"));
1143       ReinitializeAutoPtr(peminus,new ClSlackVariable (_slackCounter, "em"));
1144
1145       pexpr->setVariable(*peplus,-1.0);
1146       pexpr->setVariable(*peminus,1.0);
1147       // index the constraint under one of the error variables
1148       _markerVars[pcn] = peplus.get();
1149       _constraintsMarked[peplus.get()] = pcn;
1150
1151       ClLinearExpression *pzRow = RowExpression(_objective);
1152       // FIXGJB: pzRow->AddVariable(eplus,pcn->strength().symbolicWeight() * pcn->weight());
1153       ClSymbolicWeight sw = pcn->strength().symbolicWeight().Times(pcn->weight());
1154       double swCoeff = sw.AsDouble();
1155 #ifdef CL_TRACE
1156       if (swCoeff == 0)
1157         {
1158         cerr << "sw == " << sw << endl
1159              << "cn == " << *pcn << endl;
1160         cerr << "adding " << *peplus << " and " << *peminus
1161              << " with swCoeff == " << swCoeff << endl;
1162         }
1163 #endif
1164       pzRow->setVariable(*peplus,swCoeff);
1165       NoteAddedVariable(*peplus,_objective);
1166       // FIXGJB: pzRow->AddVariable(eminus,pcn->strength().symbolicWeight() * pcn->weight());
1167       pzRow->setVariable(*peminus,swCoeff);
1168       NoteAddedVariable(*peminus,_objective);
1169       _errorVars[pcn].insert(peminus.get());
1170       _errorVars[pcn].insert(peplus.get());
1171       if (pcn->isStayConstraint())
1172         {
1173         _stayPlusErrorVars.push_back(peplus.get());
1174         _stayMinusErrorVars.push_back(peminus.get());
1175         }
1176       else if (pcn->IsEditConstraint())
1177         {
1178         clvEplus = peplus.get();
1179         clvEminus = peminus.get();
1180         prevEConstant = cnExpr.Constant();
1181         }
1182       }
1183     }
1184
1185   // the Constant in the Expression should be non-negative.
1186   // If necessary normalize the Expression by multiplying by -1
1187   if (pexpr->Constant() < 0)
1188     {
1189 #ifdef CL_TRACE
1190     cerr << "NewExpression's Constant is " << pexpr->Constant() << ", < 0, so flipping" << endl;
1191 #endif
1192     pexpr->MultiplyMe(-1);
1193     }
1194 #ifdef CL_TRACE
1195   cerr << "- returning " << *pexpr << endl;
1196 #endif
1197   // Terrible Name -- release() does *not* delete the object,
1198   // only makes sure that the destructor won't delete the object
1199   // (it releases the cl_auto_ptr from the responsibility of deleting the object)
1200   pslackVar.release();
1201   pdummyVar.release();
1202   peminus.release();
1203   peplus.release();
1204   return pexpr.release();
1205 }
1206
1207 // Minimize the value of the objective.  (The tableau should already
1208 // be feasible.)
1209 void
1210 ClSimplexSolver::Optimize(ClVariable zVar)
1211 {
1212 #ifdef CL_TRACE
1213   Tracer TRACER(__FUNCTION__);
1214   cerr << "(" << zVar << ")\n"
1215        << *this << endl;
1216 #endif
1217   ClLinearExpression *pzRow = RowExpression(zVar);
1218   assert(pzRow != NULL);
1219   ClVariable entryVar = clvNil;
1220   ClVariable exitVar = clvNil;
1221   while (true)
1222     {
1223     Number objectiveCoeff = 0;
1224     // Find the most negative coefficient in the objective function
1225     // (ignoring the non-pivotable dummy variables).  If all
1226     // coefficients are positive we're done
1227     ClVarToNumberMap &terms = pzRow->Terms();
1228     ClVarToNumberMap::iterator it = terms.begin();
1229     for (; it != terms.end(); ++it)
1230       {
1231       ClVariable v = (*it).first;
1232       Number c = (*it).second;
1233       if (v.IsPivotable() && c < objectiveCoeff)
1234         {
1235         objectiveCoeff = c;
1236         entryVar = v;
1237         // A. Beurive' Tue Jul 13 23:03:05 CEST 1999 Why the most
1238         // negative?  I encountered unending cycles of pivots!
1239         break;
1240         }
1241       }
1242     // if all coefficients were positive (or if the objective
1243     // function has no pivotable variables)
1244     // we are at an optimum
1245     if (objectiveCoeff >= -_epsilon)
1246       return;
1247 #ifdef CL_TRACE
1248     cerr << "entryVar == " << entryVar << ", "
1249          << "objectiveCoeff == " << objectiveCoeff
1250          << endl;
1251 #endif
1252
1253     // choose which variable to move out of the basis
1254     // Only consider pivotable basic variables
1255     // (i.e. restricted, non-dummy variables)
1256     double minRatio = DBL_MAX;
1257     ClVarSet &columnVars = _columns[entryVar];
1258     ClVarSet::iterator it_rowvars = columnVars.begin();
1259     Number r = 0.0;
1260     for (; it_rowvars != columnVars.end(); ++it_rowvars)
1261       {
1262       ClVariable v = *it_rowvars;
1263 #ifdef CL_TRACE
1264       cerr << "Checking " << v << endl;
1265 #endif
1266       if (v.IsPivotable())
1267         {
1268         const ClLinearExpression *pexpr = RowExpression(v);
1269         Number coeff = pexpr->CoefficientFor(entryVar);
1270         // only consider negative coefficients
1271         if (coeff < 0.0)
1272           {
1273           r = - pexpr->Constant() / coeff;
1274           if (r < minRatio)
1275             {
1276 #ifdef CL_TRACE
1277             cerr << "New minRatio == " << r << endl;
1278 #endif
1279             minRatio = r;
1280             exitVar = v;
1281             }
1282           }
1283         }
1284       }
1285     // If minRatio is still nil at this point, it means that the
1286     // objective function is unbounded, i.e. it can become
1287     // arbitrarily negative.  This should never happen in this
1288     // application.
1289     if (minRatio == DBL_MAX)
1290       {
1291       stringstream ss;
1292       ss << "objective function is unbounded!" << ends;
1293       throw ExCLInternalError(ss.str().c_str());
1294       }
1295     Pivot(entryVar, exitVar);
1296 #ifdef CL_TRACE
1297     cerr << "After Optimize:\n"
1298          << *this << endl;
1299 #endif
1300     }
1301 }
1302
1303 // Do a Pivot.  Move entryVar into the basis (i.e. make it a basic variable),
1304 // and move exitVar out of the basis (i.e., make it a parametric variable)
1305 void
1306 ClSimplexSolver::Pivot(ClVariable entryVar, ClVariable exitVar)
1307 {
1308 #ifdef CL_TRACE
1309   Tracer TRACER(__FUNCTION__);
1310   cerr << "(" << entryVar << ", " << exitVar << ")" << endl;
1311 #endif
1312
1313   // the entryVar might be non-pivotable if we're doing a RemoveConstraint --
1314   // otherwise it should be a pivotable variable -- enforced at call sites,
1315   // hopefully
1316
1317   // expr is the Expression for the exit variable (about to leave the basis) --
1318   // so that the old tableau includes the equation:
1319   //   exitVar = expr
1320   ClLinearExpression *pexpr = RemoveRow(exitVar);
1321
1322   // Compute an Expression for the entry variable.  Since expr has
1323   // been deleted from the tableau we can destructively modify it to
1324   // build this Expression.
1325   pexpr->ChangeSubject(exitVar,entryVar);
1326   SubstituteOut(entryVar,*pexpr);
1327
1328   if (entryVar.IsExternal())
1329     {
1330     // entry var is no longer a parametric variable since we're moving
1331     // it into the basis
1332     _externalParametricVars.erase(entryVar);
1333     }
1334   addRow(entryVar,*pexpr);
1335 }
1336
1337
1338
1339 // Each of the non-required stays will be represented by an equation
1340 // of the form
1341 //     v = c + eplus - eminus
1342 // where v is the variable with the stay, c is the previous value of
1343 // v, and eplus and eminus are slack variables that hold the error
1344 // in satisfying the stay constraint.  We are about to change
1345 // something, and we want to fix the constants in the equations
1346 // representing the stays.  If both eplus and eminus are nonbasic
1347 // they have value 0 in the current solution, meaning the previous
1348 // stay was exactly satisfied.  In this case nothing needs to be
1349 // changed.  Otherwise one of them is basic, and the other must
1350 // occur only in the Expression for that basic error variable.
1351 // Reset the Constant in this Expression to 0.
1352 void
1353 ClSimplexSolver::ResetStayConstants()
1354 {
1355 #ifdef CL_TRACE
1356   Tracer TRACER(__FUNCTION__);
1357   cerr << "()" << endl;
1358 #endif
1359   ClVarVector::const_iterator
1360     itStayPlusErrorVars = _stayPlusErrorVars.begin();
1361   ClVarVector::const_iterator
1362     itStayMinusErrorVars = _stayMinusErrorVars.begin();
1363
1364   for ( ; itStayPlusErrorVars != _stayPlusErrorVars.end();
1365         ++itStayPlusErrorVars, ++itStayMinusErrorVars )
1366     {
1367     ClLinearExpression *pexpr = RowExpression(*itStayPlusErrorVars);
1368     if (pexpr == NULL )
1369       {
1370       pexpr = RowExpression(*itStayMinusErrorVars);
1371       }
1372     if (pexpr != NULL)
1373       {
1374       pexpr->Set_constant(0.0);
1375       }
1376     }
1377 }
1378
1379 // Set the external variables known to this solver to their appropriate values.
1380 // Set each external basic variable to its value, and set each
1381 // external parametric variable to 0.  (It isn't clear that we will
1382 // ever have external parametric variables -- every external
1383 // variable should either have a stay on it, or have an equation
1384 // that defines it in terms of other external variables that do have
1385 // stays.  For the moment I'll put this in though.)  Variables that
1386 // are internal to the solver don't actually store values -- their
1387 // values are just implicit in the tableu -- so we don't need to set
1388 // them."
1389 void
1390 ClSimplexSolver::SetExternalVariables()
1391 {
1392 #ifdef CL_TRACE
1393   Tracer TRACER(__FUNCTION__);
1394   cerr << "()\n"
1395        << *this << endl;
1396 #endif
1397
1398   // FIXGJB -- oughta check some invariants here
1399
1400   // Set external parametric variables first
1401   // in case I've screwed up
1402   ClVarSet::iterator itParVars = _externalParametricVars.begin();
1403   for ( ; itParVars != _externalParametricVars.end(); ++itParVars )
1404     {
1405     ClVariable v = *itParVars;
1406 #ifndef NDEBUG
1407     // defensively skip it if it is basic -- ChangeValue is virtual
1408     // so don't want to call it twice;  this should never
1409     // happen
1410     if (FIsBasicVar(v))
1411       {
1412 #ifndef CL_NO_IO
1413       // WARNING
1414       cerr << __FUNCTION__ << "Error: variable " << v
1415            << " in _externalParametricVars is basic" << endl;
1416       cerr << "Row is: " << *RowExpression(v) << endl;
1417 #endif
1418       continue;
1419       }
1420 #endif
1421     ChangeClv(v,0.0);
1422     }
1423
1424   // Only iterate over the rows w/ external variables
1425   ClVarSet::iterator itRowVars = _externalRows.begin();
1426   for ( ; itRowVars != _externalRows.end() ; ++itRowVars )
1427     {
1428     ClVariable v = *itRowVars;
1429     ClLinearExpression *pexpr = RowExpression(v);
1430     ChangeClv(v,pexpr->Constant());
1431     }
1432
1433   _fNeedsSolving = false;
1434   if (_pfnResolveCallback)
1435     _pfnResolveCallback(this);
1436 }
1437
1438 #ifndef CL_NO_IO
1439 ostream &
1440 PrintTo(ostream &xo, const ClVarVector &varlist)
1441 {
1442   ClVarVector::const_iterator it = varlist.begin();
1443   xo << varlist.size() << ":" << "[ ";
1444   if (it != varlist.end())
1445     {
1446     xo << *it;
1447     ++it;
1448     }
1449   for (; it != varlist.end(); ++it)
1450     {
1451     xo << ", " << *it;
1452     }
1453   xo << " ]";
1454   return xo;
1455 }
1456
1457 ostream &operator<<(ostream &xo, const ClVarVector &varlist)
1458 { return PrintTo(xo,varlist); }
1459
1460
1461 ostream &
1462 PrintTo(ostream &xo, const ClConstraintToVarSetMap &mapCnToVarSet)
1463 {
1464   ClConstraintToVarSetMap::const_iterator it = mapCnToVarSet.begin();
1465   for ( ; it != mapCnToVarSet.end(); ++it) {
1466     const ClConstraint *pcn = (*it).first;
1467     const ClVarSet &set = (*it).second;
1468     xo << "CN: " << pcn << *pcn << ":: " << set << endl;
1469   }
1470   return xo;
1471 }
1472
1473 ostream &operator <<(ostream &xo, const ClConstraintToVarSetMap &mapCnToVarSet)
1474 { return PrintTo(xo,mapCnToVarSet); }
1475
1476
1477
1478 ostream &
1479 ClSimplexSolver::PrintOn(ostream &xo) const
1480 {
1481   ClTableau::PrintOn(xo);
1482
1483   xo << "_stayPlusErrorVars: "
1484      << _stayPlusErrorVars << endl;
1485   xo << "_stayMinusErrorVars: "
1486      << _stayMinusErrorVars << endl;
1487   xo << "_editInfoList:\n"
1488      << _editInfoList << endl;
1489   return xo;
1490 }
1491
1492
1493 ostream &
1494 ClSimplexSolver::PrintInternalInfo(ostream &xo) const
1495 {
1496   ClTableau::PrintInternalInfo(xo);
1497   xo << "; edvars: " << _editInfoList.size();
1498   xo << endl;
1499   printExternalVariablesTo(xo);
1500   return xo;
1501 }
1502
1503 ostream &operator<<(ostream &xo, const ClSimplexSolver &clss)
1504 {
1505   return clss.PrintOn(xo);
1506 }
1507
1508 #endif
1509
1510 bool
1511 ClSimplexSolver::FIsConstraintSatisfied(const ClConstraint *const pcn) const
1512 {
1513   ClConstraintToVarMap::const_iterator it_marker = _markerVars.find(pcn);
1514   if (it_marker == _markerVars.end())
1515     { // could not find the constraint
1516     throw ExCLConstraintNotFound();
1517     }
1518
1519 #ifndef CL_NO_IO
1520   bool fCnsays = pcn->FIsSatisfied();
1521 #endif
1522
1523   ClConstraintToVarSetMap::const_iterator it_eVars = _errorVars.find(pcn);
1524
1525   if (it_eVars != _errorVars.end())
1526     {
1527     const ClVarSet &eVars = (*it_eVars).second;
1528     ClVarSet::const_iterator it = eVars.begin();
1529     for ( ; it != eVars.end(); ++it )
1530       {
1531       const ClLinearExpression *pexpr = RowExpression(*it);
1532       if (pexpr != NULL && !ClApprox(pexpr->Constant(),0.0))
1533         {
1534 #ifndef CL_NO_IO
1535         if (fCnsays)
1536           cerr << __FUNCTION__ << ": constraint says satisfiable, but solver does not" << endl;
1537 #endif
1538         return false;
1539         }
1540       }
1541     }
1542
1543 #ifndef CL_NO_IO
1544   if (!fCnsays)
1545     cerr << __FUNCTION__ << ": solver says satisfiable, but constraint does not" << endl;
1546 #endif
1547   return true;
1548 }
1549
1550
1551
1552 #ifndef CL_NO_ID
1553
1554 ostream &PrintTo(ostream &xo, const ClSimplexSolver::ClEditInfoList &listPEditInfo)
1555 {
1556   ClSimplexSolver::ClEditInfoList::const_iterator it = listPEditInfo.begin();
1557   for ( ; it != listPEditInfo.end(); ++it) {
1558     const ClSimplexSolver::ClEditInfo *pcei = (*it);
1559     xo << *pcei << endl;
1560   }
1561   return xo;
1562 }
1563
1564
1565 ostream &operator<<(ostream &xo, const ClSimplexSolver::ClEditInfoList &listPEditInfo)
1566 { return PrintTo(xo,listPEditInfo); }
1567
1568 #endif
1569
1570 // A. Beurive' Tue Jul  6 17:03:32 CEST 1999
1571 void
1572 ClSimplexSolver::ChangeStrengthAndWeight(ClConstraint *pcn, const ClStrength &strength, double weight)
1573 {
1574   ClConstraintToVarSetMap::iterator it_eVars = _errorVars.find(pcn);
1575   // Only for constraints that already have error variables (i.e. non-required constraints)
1576   assert(it_eVars != _errorVars.end());
1577
1578   ClLinearExpression *pzRow = RowExpression(_objective);
1579
1580   Number old_coeff = pcn->weight() * pcn->strength().symbolicWeight().AsDouble();
1581   pcn->setStrength(strength);
1582   pcn->setWeight(weight);
1583   Number new_coeff = pcn->weight() * pcn->strength().symbolicWeight().AsDouble();
1584
1585   if (new_coeff != old_coeff)
1586     {
1587 #ifdef CL_TRACE
1588       cerr << "Changing strength and/or weight for constraint: " << endl << *pcn << endl;
1589       cerr << "Updating objective row from:" << endl << *pzRow << endl;
1590 #endif
1591       ClVarSet &eVars = (*it_eVars).second;
1592       ClVarSet::iterator it = eVars.begin();
1593       for ( ; it != eVars.end(); ++it )
1594         {
1595           const ClLinearExpression *pexpr = RowExpression(*it);
1596           if (pexpr == NULL )
1597             {
1598               pzRow->AddVariable(*it,-old_coeff,_objective,*this);
1599               pzRow->AddVariable(*it,new_coeff,_objective,*this);
1600             }
1601           else
1602             {
1603               pzRow->AddExpression(*pexpr,-old_coeff,_objective,*this);
1604               pzRow->AddExpression(*pexpr,new_coeff,_objective,*this);
1605             }
1606         }
1607 #ifdef CL_TRACE
1608       cerr << "to: " << endl << *pzRow << endl;
1609 #endif
1610
1611       if (_fAutosolve)
1612         {
1613           Optimize(_objective);
1614           SetExternalVariables();
1615         }
1616     }
1617 }
1618
1619 // A. Beurive' Tue Jul  6 17:03:42 CEST 1999
1620 void
1621 ClSimplexSolver::ChangeStrength(ClConstraint *pcn, const ClStrength &strength)
1622 {
1623   ChangeStrengthAndWeight(pcn,strength,pcn->weight());
1624 }
1625
1626 // A. Beurive' Tue Jul  6 17:03:42 CEST 1999
1627 void
1628 ClSimplexSolver::ChangeWeight(ClConstraint *pcn, double weight)
1629 {
1630   ChangeStrengthAndWeight(pcn,pcn->strength(),weight);
1631 }