Changeset 92


Ignore:
Timestamp:
Jan 13, 2003 10:05:41 AM (17 years ago)
Author:
forrest
Message:

Free variables done better (I hope)

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpPrimalColumnSteepest.cpp

    r89 r92  
    235235        case ClpSimplex::isFree:
    236236        case ClpSimplex::superBasic:
    237           if (fabs(value)>tolerance) {
    238             // we are going to bias towards free
    239             value *= 10.0;
     237          if (fabs(value)>1.0e2*tolerance) {
     238            // we are going to bias towards free (but only if reasonable)
     239            value *= 1000.0;
    240240            // store square in list
    241241            if (infeas[iSequence+addSequence])
     
    295295      case ClpSimplex::isFree:
    296296      case ClpSimplex::superBasic:
    297         if (fabs(value)>tolerance) {
    298           // we are going to bias towards free
    299           value *= 10.0;
     297        if (fabs(value)>1.0e2*tolerance) {
     298          // we are going to bias towards free (but only if reasonable)
     299          value *= 1000.0;
    300300          // store square in list
    301301          if (infeas[sequenceOut])
     
    360360      case ClpSimplex::isFree:
    361361        if (fabs(value)>tolerance) {
    362           // we are going to bias towards free
    363           value *= 10.0;
     362          // we are going to bias towards free (but only if reasonable)
     363          value *= 1000.0;
    364364          // store square in list
    365365          if (infeas[iSequence+addSequence])
     
    705705      case ClpSimplex::isFree:
    706706      case ClpSimplex::superBasic:
    707         if (fabs(value)>tolerance) {
     707        if (fabs(value)>1.0e2*tolerance) {
     708          // we are going to bias towards free (but only if reasonable)
    708709          // store square in list
    709           infeasible_->quickAdd(iSequence,value*value);
     710          infeasible_->quickAdd(iSequence,1.0e6*value*value);
    710711        }
    711712        break;
  • trunk/ClpSimplex.cpp

    r85 r92  
    100100  numberFake_(0),
    101101  progressFlag_(0),
     102  firstFree_(-1),
    102103  sumOfRelaxedDualInfeasibilities_(0.0),
    103104  sumOfRelaxedPrimalInfeasibilities_(0.0)
     
    484485    solveType -= 10;
    485486  }
     487#ifdef CLP_DEBUG
     488  if (solveType==1) {
     489    int numberFreeIn=0,numberFreeOut=0;
     490    double biggestDj=0.0;
     491    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     492      switch(getColumnStatus(iColumn)) {
     493       
     494      case basic:
     495        if (columnLower_[iColumn]<-largeValue_
     496            &&columnUpper_[iColumn]>largeValue_)
     497          numberFreeIn++;
     498        break;
     499      default:
     500        if (columnLower_[iColumn]<-largeValue_
     501            &&columnUpper_[iColumn]>largeValue_) {
     502          numberFreeOut++;
     503          biggestDj = max(fabs(dj_[iColumn]),biggestDj);
     504        }
     505        break;
     506      }
     507    }
     508    if (numberFreeIn+numberFreeOut)
     509      printf("%d in basis, %d out - largest dj %g\n",
     510             numberFreeIn,numberFreeOut,biggestDj);
     511  }
     512#endif
    486513  if (!solveType) {
    487514    if (!valuesPass) {
     
    566593        }
    567594        totalSlacks=numberBasic;
     595
    568596        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    569597          switch(getColumnStatus(iColumn)) {
     
    641669        }
    642670      } else {
    643         //#define TESTFREE
    644671        // all slack basis
    645672        int numberBasic=0;
     
    663690            rowActivityWork_[iRow]=0.0;
    664691          }
    665 #ifdef TESTFREE
    666692          setRowStatus(iRow,basic);
    667693          numberBasic++;
    668 #endif
    669694        }
    670695        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     
    681706            }
    682707          } else {
    683 #ifndef TESTFREE
    684             numberBasic++;
    685             setColumnStatus(iColumn,basic);
    686 #else
    687708            setColumnStatus(iColumn,isFree);
    688 #endif
    689709            columnActivityWork_[iColumn]=0.0;
    690710          }
     
    11331153  numberFake_(0),
    11341154  progressFlag_(0),
     1155  firstFree_(-1),
    11351156  sumOfRelaxedDualInfeasibilities_(0.0),
    11361157  sumOfRelaxedPrimalInfeasibilities_(0.0)
     
    12211242  numberFake_(0),
    12221243  progressFlag_(0),
     1244  firstFree_(-1),
    12231245  sumOfRelaxedDualInfeasibilities_(0.0),
    12241246  sumOfRelaxedPrimalInfeasibilities_(0.0)
     
    13491371  numberFake_ = rhs.numberFake_;
    13501372  progressFlag_ = rhs.progressFlag_;
     1373  firstFree_ = rhs.firstFree_;
    13511374  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    13521375  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
     
    14271450
    14281451  objectiveValue_ = 0.0;
     1452  firstFree_ = -1;
    14291453  // now look at primal solution
    14301454  columnPrimalInfeasibility_=0.0;
     
    14451469
    14461470  for (iRow=0;iRow<numberRows_;iRow++) {
     1471    //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
    14471472    double infeasibility=0.0;
    14481473    objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];
     
    14721497  solution = columnActivityWork_;
    14731498  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1499    //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
    14741500    double infeasibility=0.0;
    14751501    objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];
     
    15021528{
    15031529
    1504   double * solution;
    15051530  int iRow,iColumn;
    15061531  sumDualInfeasibilities_=0.0;
     
    15111536  rowDualInfeasibility_=0.0;
    15121537  rowDualSequence_=-1;
     1538  firstFree_ = -1;
    15131539  primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_,
    15141540                                   columnPrimalInfeasibility_);
    15151541  remainingDualInfeasibility_=0.0;
    1516   solution = rowActivityWork_;
    15171542  double relaxedTolerance=dualTolerance_;
    15181543  // we can't really trust infeasibilities if there is dual error
     
    15221547  sumOfRelaxedDualInfeasibilities_ = 0.0;
    15231548
     1549  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1550    if (getColumnStatus(iColumn) != basic) {
     1551      // not basic
     1552      double value = reducedCostWork_[iColumn];
     1553      double distanceUp = columnUpperWork_[iColumn]-
     1554        columnActivityWork_[iColumn];
     1555      double distanceDown = columnActivityWork_[iColumn] -
     1556        columnLowerWork_[iColumn];
     1557      if (distanceUp>primalTolerance_) {
     1558        // Check if "free"
     1559        if (firstFree_<0&&distanceDown>primalTolerance_) {
     1560          if (algorithm_>0) {
     1561            // primal
     1562            firstFree_ = iColumn;
     1563          } else if (fabs(value)>1.0e2*relaxedTolerance) {
     1564            // dual with dj
     1565            firstFree_ = iColumn;
     1566          }
     1567        }
     1568        // should not be negative
     1569        if (value<0.0) {
     1570          value = - value;
     1571          if (value>columnDualInfeasibility_) {
     1572            columnDualInfeasibility_=value;
     1573            columnDualSequence_=iColumn;
     1574          }
     1575          if (value>dualTolerance_) {
     1576            if (getColumnStatus(iColumn) != isFree) {
     1577              numberDualInfeasibilitiesWithoutFree_ ++;
     1578              sumDualInfeasibilities_ += value-dualTolerance_;
     1579              if (value>relaxedTolerance)
     1580                sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
     1581              numberDualInfeasibilities_ ++;
     1582            } else {
     1583              // free so relax a lot
     1584              value *= 0.01;
     1585              if (value>dualTolerance_) {
     1586                sumDualInfeasibilities_ += value-dualTolerance_;
     1587                if (value>relaxedTolerance)
     1588                  sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
     1589                numberDualInfeasibilities_ ++;
     1590              }
     1591            }
     1592            // maybe we can make feasible by increasing tolerance
     1593            if (distanceUp<largeValue_) {
     1594              if (distanceUp>primalToleranceToGetOptimal_)
     1595                primalToleranceToGetOptimal_=distanceUp;
     1596            } else {
     1597              //gap too big for any tolerance
     1598              remainingDualInfeasibility_=
     1599                max(remainingDualInfeasibility_,value);
     1600            }
     1601          }
     1602        }
     1603      }
     1604      if (distanceDown>primalTolerance_) {
     1605        // should not be positive
     1606        if (value>0.0) {
     1607          if (value>columnDualInfeasibility_) {
     1608            columnDualInfeasibility_=value;
     1609            columnDualSequence_=iColumn;
     1610          }
     1611          if (value>dualTolerance_) {
     1612            sumDualInfeasibilities_ += value-dualTolerance_;
     1613            if (value>relaxedTolerance)
     1614              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
     1615            numberDualInfeasibilities_ ++;
     1616            if (getColumnStatus(iColumn) != isFree)
     1617              numberDualInfeasibilitiesWithoutFree_ ++;
     1618            // maybe we can make feasible by increasing tolerance
     1619            if (distanceDown<largeValue_&&
     1620                distanceDown>primalToleranceToGetOptimal_)
     1621              primalToleranceToGetOptimal_=distanceDown;
     1622          }
     1623        }
     1624      }
     1625    }
     1626  }
    15241627  for (iRow=0;iRow<numberRows_;iRow++) {
    15251628    if (getRowStatus(iRow) != basic) {
    15261629      // not basic
    15271630      double value = rowReducedCost_[iRow];
    1528       double distance;
    1529       distance = rowUpperWork_[iRow]-solution[iRow];
    1530       if (distance>primalTolerance_) {
     1631      double distanceUp = rowUpperWork_[iRow]-rowActivityWork_[iRow];
     1632      double distanceDown = rowActivityWork_[iRow] -rowLowerWork_[iRow];
     1633      if (distanceUp>primalTolerance_) {
     1634        // Check if "free"
     1635        if (firstFree_<0&&distanceDown>primalTolerance_) {
     1636          if (algorithm_>0) {
     1637            // primal
     1638            firstFree_ = iRow+numberColumns_;
     1639          } else if (fabs(value)>1.0e2*relaxedTolerance) {
     1640            // dual with dj
     1641            firstFree_ = iRow+numberColumns_;
     1642          }
     1643        }
    15311644        // should not be negative
    15321645        if (value<0.0) {
     
    15441657              numberDualInfeasibilitiesWithoutFree_ ++;
    15451658            // maybe we can make feasible by increasing tolerance
    1546             if (distance<largeValue_) {
    1547               if (distance>primalToleranceToGetOptimal_)
    1548                 primalToleranceToGetOptimal_=distance;
     1659            if (distanceUp<largeValue_) {
     1660              if (distanceUp>primalToleranceToGetOptimal_)
     1661                primalToleranceToGetOptimal_=distanceUp;
    15491662            } else {
    15501663              //gap too big for any tolerance
     
    15551668        }
    15561669      }
    1557       distance = solution[iRow] -rowLowerWork_[iRow];
    1558       if (distance>primalTolerance_) {
     1670      if (distanceDown>primalTolerance_) {
    15591671        // should not be positive
    15601672        if (value>0.0) {
     
    15711683              numberDualInfeasibilitiesWithoutFree_ ++;
    15721684            // maybe we can make feasible by increasing tolerance
    1573             if (distance<largeValue_&&
    1574                 distance>primalToleranceToGetOptimal_)
    1575               primalToleranceToGetOptimal_=distance;
    1576           }
    1577         }
    1578       }
    1579     }
    1580   }
    1581   solution = columnActivityWork_;
    1582   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1583     if (getColumnStatus(iColumn) != basic) {
    1584       // not basic
    1585       double value = reducedCostWork_[iColumn];
    1586       double distance;
    1587       distance = columnUpperWork_[iColumn]-solution[iColumn];
    1588       if (distance>primalTolerance_) {
    1589         // should not be negative
    1590         if (value<0.0) {
    1591           value = - value;
    1592           if (value>columnDualInfeasibility_) {
    1593             columnDualInfeasibility_=value;
    1594             columnDualSequence_=iColumn;
    1595           }
    1596           if (value>dualTolerance_) {
    1597             sumDualInfeasibilities_ += value-dualTolerance_;
    1598             if (value>relaxedTolerance)
    1599               sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    1600             numberDualInfeasibilities_ ++;
    1601             if (getColumnStatus(iColumn) != isFree)
    1602               numberDualInfeasibilitiesWithoutFree_ ++;
    1603             // maybe we can make feasible by increasing tolerance
    1604             if (distance<largeValue_) {
    1605               if (distance>primalToleranceToGetOptimal_)
    1606                 primalToleranceToGetOptimal_=distance;
    1607             } else {
    1608               //gap too big for any tolerance
    1609               remainingDualInfeasibility_=
    1610                 max(remainingDualInfeasibility_,value);
    1611             }
    1612           }
    1613         }
    1614       }
    1615       distance = solution[iColumn] -columnLowerWork_[iColumn];
    1616       if (distance>primalTolerance_) {
    1617         // should not be positive
    1618         if (value>0.0) {
    1619           if (value>columnDualInfeasibility_) {
    1620             columnDualInfeasibility_=value;
    1621             columnDualSequence_=iColumn;
    1622           }
    1623           if (value>dualTolerance_) {
    1624             sumDualInfeasibilities_ += value-dualTolerance_;
    1625             if (value>relaxedTolerance)
    1626               sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    1627             numberDualInfeasibilities_ ++;
    1628             if (getColumnStatus(iColumn) != isFree)
    1629               numberDualInfeasibilitiesWithoutFree_ ++;
    1630             // maybe we can make feasible by increasing tolerance
    1631             if (distance<largeValue_&&
    1632                 distance>primalToleranceToGetOptimal_)
    1633               primalToleranceToGetOptimal_=distance;
     1685            if (distanceDown<largeValue_&&
     1686                distanceDown>primalToleranceToGetOptimal_)
     1687              primalToleranceToGetOptimal_=distanceDown;
    16341688          }
    16351689        }
  • trunk/ClpSimplexDual.cpp

    r90 r92  
    234234    factorization_->slackValue(-1.0);
    235235    factorization_->zeroTolerance(1.0e-13);
    236    
     236    // might want to do first time - factorization_->pivotTolerance(0.99);
     237
     238    // Do Initial factorization
    237239    int factorizationStatus = internalFactorize(0);
     240
    238241    if (factorizationStatus<0)
    239242      return 1; // some error
     
    248251      perturb();
    249252   
    250     double objectiveChange;
    251253    // for dual we will change bounds using dualBound_
    252254    // for this we need clean basis so it is after factorize
    253255    gutsOfSolution(rowActivityWork_,columnActivityWork_);
    254256   
     257    double objectiveChange;
    255258    numberFake_ =0; // Number of variables at fake bounds
    256259    changeBounds(true,NULL,objectiveChange);
     
    325328  perturbation_ = savePerturbation;
    326329  dualBound_ = saveDualBound_;
     330  factorization_->relaxAccuracyCheck(1.0);
    327331  return problemStatus_;
    328332}
     
    341345  int debugIteration=-1;
    342346#endif
     347  // if can't trust much and long way from optimal then relax
     348  if (largestPrimalError_>10.0)
     349    factorization_->relaxAccuracyCheck(min(1.0e2,largestPrimalError_/10.0));
     350  else
     351    factorization_->relaxAccuracyCheck(1.0);
    343352  // status stays at -1 while iterating, >=0 finished, -2 to invert
    344353  // status -3 to go to top without an invert
     
    452461          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
    453462#endif
     463        double checkValue=1.0e-7;
     464        // if can't trust much and long way from optimal then relax
     465        if (largestPrimalError_>10.0)
     466          checkValue = min(1.0e-4,1.0e-8*largestPrimalError_);
    454467        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
    455             fabs(btranAlpha-alpha_)>1.0e-7*(1.0+fabs(alpha_))) {
     468            fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
    456469          handler_->message(CLP_DUAL_CHECK,messages_)
    457470            <<btranAlpha
    458471            <<alpha_
    459472            <<CoinMessageEol;
    460           dualRowPivot_->unrollWeights();
    461473          if (factorization_->pivots()) {
     474            dualRowPivot_->unrollWeights();
    462475            problemStatus_=-2; // factorize now
    463476            rowArray_[0]->clear();
     
    470483            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
    471484                fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
     485              dualRowPivot_->unrollWeights();
    472486              // need to reject something
    473487              char x = isColumn(sequenceOut_) ? 'C' :'R';
     
    516530          assert(dualOut_<=oldOut);
    517531#endif
    518           if(dualOut_<0.0&&factorization_->pivots()) {
     532          if(dualOut_<0.0&&factorization_->pivots()&&
     533             getStatus(sequenceIn_)!=isFree) {
    519534            // going backwards - factorize
    520535            dualRowPivot_->unrollWeights();
     
    543558          }
    544559        }
    545         assert(fabs(dualOut_)<1.0e30);
     560        assert(fabs(dualOut_)<1.0e50);
    546561        // if stable replace in basis
    547562        int updateStatus = factorization_->replaceColumn(rowArray_[2],
     
    962977          upper_[iSequence]=newUpperValue;
    963978          if (newLowerValue > lowerValue) {
    964             if (newUpperValue < upperValue)
     979            if (newUpperValue < upperValue) {
    965980              setFakeBound(iSequence,ClpSimplexDual::bothFake);
    966             else
     981            } else {
    967982              setFakeBound(iSequence,ClpSimplexDual::lowerFake);
     983            }
    968984          } else {
    969985            if (newUpperValue < upperValue)
     
    10221038            setFakeBound(iSequence,ClpSimplexDual::bothFake);
    10231039          }
     1040        } else {
     1041          // set non basic free variables to fake bounds
     1042          // I don't think we should ever get here
     1043          assert(!("should not be here"));
     1044          lower_[iSequence]=-0.5*dualBound_;
     1045          upper_[iSequence]= 0.5*dualBound_;
     1046          setFakeBound(iSequence,ClpSimplexDual::bothFake);
     1047          setStatus(iSequence,atUpperBound);
     1048          solution_[iSequence]=0.5*dualBound_;
    10241049        }
    10251050      }
    10261051    }
     1052
    10271053    return 1;
    10281054  }
     
    11691195            keep = 2;
    11701196        } else {
    1171           if (alpha>=acceptablePivot)
     1197          if (fabs(alpha)>10.0*acceptablePivot)
    11721198            keep = 2;
    1173           else if (-alpha>=acceptablePivot)
    1174             keep = 2;
     1199
     1200
    11751201        }
    11761202        break;
     
    12131239  marker[0][0] = numberRemaining;
    12141240
    1215   if (!numberRemaining)
     1241  if (!numberRemaining&&sequenceIn_<0)
    12161242    return; // Looks infeasible
    12171243
     
    15251551      oldValue = dj_[sequenceIn_];
    15261552      theta_ = oldValue/alpha;
    1527 #if 0
    1528       if (numberIterations_>2000)
    1529         exit(99);
    1530       if (numberIterations_>2000-20)
    1531         handler_->setLogLevel(63);
    1532       if (numberIterations_>2000-20)
    1533         printf("theta %g oldValue %g tol %g %g\n",theta_,oldValue,dualTolerance_,
    1534                newTolerance);
    1535 #endif
    15361553      if (theta_<MINIMUMTHETA) {
    15371554        // can't pivot to zero
     
    15601577            alpha=workColumn[iSequence];
    15611578          double value = dj_[iSequence]-theta_*alpha;
    1562 #if 0
    1563           if (numberIterations_>2000-20)
    1564             printf("%d alpha %g value %g\n",iSequence,alpha,value);
    1565 #endif
    15661579           
    15671580          // can't be free here
     
    15801593              dj_[iSequence] += modification;
    15811594              cost_[iSequence] +=  modification;
    1582 #if 0
    1583               if (numberIterations_>2000-20)
    1584                 printf("%d acost %g mod %g\n",iSequence,cost_[iSequence],
    1585                        modification);
    1586 #endif
    15871595#endif
    15881596            }
     
    15991607              dj_[iSequence] += modification;
    16001608              cost_[iSequence] +=  modification;
    1601 #if 0
    1602               if (numberIterations_>2000-20)
    1603                 printf("%d bcost %g mod %g\n",iSequence,cost_[iSequence],
    1604                        modification);
    1605 #endif
    16061609#endif
    16071610            }
     
    18481851  handler_->printing(numberDualInfeasibilitiesWithoutFree_
    18491852                     <numberDualInfeasibilities_)
    1850                        <<numberDualInfeasibilities_-
    1851     numberDualInfeasibilitiesWithoutFree_;
     1853                       <<numberDualInfeasibilitiesWithoutFree_;
    18521854  handler_->message()<<CoinMessageEol;
    18531855
     1856  // dual bound coming in
     1857  //double saveDualBound = dualBound_;
    18541858  while (problemStatus_<=-3) {
    18551859    bool cleanDuals=situationChanged;
  • trunk/ClpSimplexPrimal.cpp

    r66 r92  
    718718  handler_->printing(numberDualInfeasibilitiesWithoutFree_
    719719                     <numberDualInfeasibilities_)
    720                        <<numberDualInfeasibilities_-
    721     numberDualInfeasibilitiesWithoutFree_;
     720                       <<numberDualInfeasibilitiesWithoutFree_;
    722721  handler_->message()<<CoinMessageEol;
    723722  assert (primalFeasible());
  • trunk/Test/ClpMain.cpp

    r86 r92  
    915915    ClpDualRowSteepest steep;
    916916    models[0].setDualRowPivotAlgorithm(steep);
    917     models[0].setPrimalTolerance(1.0e-8);
     917    models[0].setPrimalTolerance(1.0e-7);
    918918    ClpPrimalColumnSteepest steepP;
    919919    models[0].setPrimalColumnPivotAlgorithm(steepP);
  • trunk/include/ClpSimplex.hpp

    r80 r92  
    726726  /// Progress flag - at present 0 bit says artificials out
    727727  int progressFlag_;
     728  /// First free/super-basic variable (-1 if none)
     729  int firstFree_;
    728730  /// Sum of Dual infeasibilities using tolerance based on error in duals
    729731  double sumOfRelaxedDualInfeasibilities_;
Note: See TracChangeset for help on using the changeset viewer.