ClpPredictorCorrector Class Reference

#include <ClpPredictorCorrector.hpp>

Inheritance diagram for ClpPredictorCorrector:

List of all members.

Public Member Functions
Description of algorithm
int	solve ()
Functions used in algorithm
double	findStepLength (const int phase)
	findStepLength.
double	findDirectionVector (const int phase)
	findDirectionVector.
void	createSolution ()
	createSolution. Creates solution from scratch
double	complementarityGap (int &numberComplementarityPairs, const int phase)
	complementarityGap. Computes gap
void	setupForSolve (const int phase)
	setupForSolve.
bool	checkGoodMove (const bool doCorrector)
bool	checkGoodMove2 (const double move)
	: checks for one step size
int	updateSolution ()
	updateSolution. Updates solution at end of iteration
double	affineProduct ()
	Save info on products of affine deltaT*deltaW and deltaS*deltaZ.

---

Detailed Description

This solves LPs using the predictor-corrector method.

It is rather basic as Interior point is not my speciality

It inherits from ClpInterior. It has no data of its own and is never created - only cast from a ClpInterior object at algorithm time.

Definition at line 24 of file ClpPredictorCorrector.hpp.

---

Member Function Documentation

int ClpPredictorCorrector::solve (   )

Primal Dual Predictor Corrector algorithm Method Big TODO Definition at line 38 of file ClpPredictorCorrector.cpp. References affineProduct(), ClpInterior::checkSolution(), complementarityGap(), createSolution(), ClpInterior::createWorkingData(), ClpModel::dualTolerance(), ClpCholeskyBase::factorize(), findDirectionVector(), findStepLength(), CoinMessageHandler::message(), ClpModel::objectiveValue(), ClpCholeskyBase::order(), ClpModel::primalTolerance(), ClpCholeskyBase::rank(), ClpCholeskyBase::resetRowsDropped(), setupForSolve(), ClpCholeskyBase::status(), ClpMatrixBase::times(), ClpMatrixBase::transposeTimes(), and updateSolution(). { problemStatus_=-1; algorithm_=1; //create all regions if (!createWorkingData()) { problemStatus_=4; return 2; } //initializeFeasible(); - this just set fixed flag smallestInfeasibility_=COIN_DBL_MAX; int i; for (i=0;i<LENGTH_HISTORY;i++) historyInfeasibility_[i]=COIN_DBL_MAX; //bool firstTime=true; //firstFactorization(true); cholesky_->order(this); mu_=1.0e10; //set iterations numberIterations_=-1; //initialize solution here createSolution(); //firstFactorization(false); CoinZeroN(dual_,numberRows_); multiplyAdd(solution_+numberColumns_,numberRows_,-1.0,errorRegion_,0.0); matrix_->times(1.0,solution_,errorRegion_); maximumRHSError_=maximumAbsElement(errorRegion_,numberRows_); maximumBoundInfeasibility_=maximumRHSError_; //double maximumDualError_=COIN_DBL_MAX; //initialize actualDualStep_=0.0; actualPrimalStep_=0.0; gonePrimalFeasible_=false; goneDualFeasible_=false; //bool hadGoodSolution=false; diagonalScaleFactor_=1.0; diagonalNorm_=solutionNorm_; mu_=solutionNorm_; int numberFixed=updateSolution(); int numberFixedTotal=numberFixed; //int numberRows_DroppedBefore=0; //double extra=eExtra; //double maximumPerturbation=COIN_DBL_MAX; //constants for infeas interior point const double beta2 = 0.99995; const double tau = 0.00002; double lastComplementarityGap=COIN_DBL_MAX; int lastGoodIteration=0; double bestObjectiveGap=COIN_DBL_MAX; int saveIteration=-1; bool sloppyOptimal=false; double * savePi=NULL; double * savePrimal=NULL; int numberTotal = numberRows_+numberColumns_; while (problemStatus_<0) { complementarityGap_=complementarityGap(numberComplementarityPairs_,0); handler_->message(CLP_BARRIER_ITERATION,messages_) <<numberIterations_ <<primalObjective_*optimizationDirection_- dblParam_[ClpObjOffset] << dualObjective_*optimizationDirection_- dblParam_[ClpObjOffset] <<complementarityGap_ <<numberFixedTotal <<cholesky_->rank() <<CoinMessageEol; double goodGapChange; if (!sloppyOptimal) { goodGapChange=0.93; } else { goodGapChange=0.7; } double gapO; double lastGood=bestObjectiveGap; if (gonePrimalFeasible_&&goneDualFeasible_) { double largestObjective; if (fabs(primalObjective_)>fabs(dualObjective_)) { largestObjective = fabs(primalObjective_); } else { largestObjective = fabs(dualObjective_); } if (largestObjective<1.0) { largestObjective=1.0; } gapO=fabs(primalObjective_-dualObjective_)/largestObjective; handler_->message(CLP_BARRIER_OBJECTIVE_GAP,messages_) <<gapO <<CoinMessageEol; //start saving best if (gapO<bestObjectiveGap) { saveIteration=numberIterations_; bestObjectiveGap=gapO; if (!savePi) { savePi=new double[numberRows_]; savePrimal = new double [numberTotal]; } CoinMemcpyN(dual_,numberRows_,savePi); CoinMemcpyN(solution_,numberTotal,savePrimal); } else if(gapO>2.0*bestObjectiveGap) { //maybe be more sophisticated e.g. re-initialize having //fixed variables and dropped rows //std::cout <<" gap increasing "<<std::endl; } //std::cout <<"could stop"<<std::endl; //gapO=0.0; if (fabs(primalObjective_-dualObjective_)<dualTolerance()) { gapO=0.0; } } else { gapO=COIN_DBL_MAX; if (saveIteration>=0) { handler_->message(CLP_BARRIER_GONE_INFEASIBLE,messages_) <<CoinMessageEol; } } if (gapO<1.0e-7&&!sloppyOptimal) { sloppyOptimal=true; handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL,messages_) <<numberIterations_<<complementarityGap_ <<CoinMessageEol; } if (complementarityGap_>=1.05*lastComplementarityGap) { handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_) <<complementarityGap_<<"increasing" <<CoinMessageEol; if (saveIteration>=0&&sloppyOptimal) { handler_->message(CLP_BARRIER_EXIT2,messages_) <<saveIteration <<CoinMessageEol; break; } else { //lastComplementarityGap=complementarityGap_; } } else if (complementarityGap_<goodGapChange*lastComplementarityGap) { lastGoodIteration=numberIterations_; lastComplementarityGap=complementarityGap_; } else if (numberIterations_-lastGoodIteration>=5&&complementarityGap_<1.0e-3) { handler_->message(CLP_BARRIER_COMPLEMENTARITY,messages_) <<complementarityGap_<<"not decreasing" <<CoinMessageEol; if (gapO>0.75*lastGood) { break; } } if (numberIterations_>maximumBarrierIterations_) { handler_->message(CLP_BARRIER_STOPPING,messages_) <<CoinMessageEol; break; } if (gapO<targetGap_) { problemStatus_=0; handler_->message(CLP_BARRIER_EXIT,messages_) <<" " <<CoinMessageEol; break;//finished } if (complementarityGap_<1.0e-18) { problemStatus_=0; handler_->message(CLP_BARRIER_EXIT,messages_) <<"- small complementarity gap" <<CoinMessageEol; break;//finished } if (complementarityGap_<1.0e-10&&gapO<1.0e-10) { problemStatus_=0; handler_->message(CLP_BARRIER_EXIT,messages_) <<"- objective gap and complementarity gap both small" <<CoinMessageEol; break;//finished } if (gapO<1.0e-9) { double value=gapO*complementarityGap_; value*=actualPrimalStep_; value*=actualDualStep_; //std::cout<<value<<std::endl; if (value<1.0e-17&&numberIterations_>lastGoodIteration) { problemStatus_=0; handler_->message(CLP_BARRIER_EXIT,messages_) <<"- objective gap and complementarity gap both smallish and small steps" <<CoinMessageEol; break;//finished } } bool useAffine=false; bool goodMove=false; bool doCorrector=true; //bool retry=false; worstDirectionAccuracy_=0.0; while (!goodMove) { goodMove=true; int newDropped=0; //Predictor step //Are we going to use the affine direction? if (!useAffine) { //no - normal //prepare for cholesky. Set up scaled diagonal in weights // ** for efficiency may be better if scale factor known before double norm2=0.0; double maximumValue; getNorms(diagonal_,numberColumns_,maximumValue,norm2); diagonalNorm_ = sqrt(norm2/numberComplementarityPairs_); diagonalScaleFactor_=1.0; double maximumAllowable=eScale; //scale so largest is less than allowable ? could do better double factor=0.5; while (maximumValue>maximumAllowable) { diagonalScaleFactor_*=factor; maximumValue*=factor; } /* endwhile */ if (diagonalScaleFactor_!=1.0) { handler_->message(CLP_BARRIER_SCALING,messages_) <<"diagonal"<<diagonalScaleFactor_ <<CoinMessageEol; diagonalNorm_*=diagonalScaleFactor_; } multiplyAdd(NULL,numberTotal,0.0,diagonal_, diagonalScaleFactor_); int * rowsDroppedThisTime = new int [numberRows_]; newDropped=cholesky_->factorize(diagonal_,rowsDroppedThisTime); if (newDropped) { //int newDropped2=cholesky_->factorize(diagonal_,rowsDroppedThisTime); //assert(!newDropped2); if (newDropped<0) { //replace dropped newDropped=-newDropped; //off 1 to allow for reset all newDropped--; //set all bits false cholesky_->resetRowsDropped(); } } delete [] rowsDroppedThisTime; if (cholesky_->status()) { std::cout << "bad cholesky?" <<std::endl; abort(); } } //set up for affine direction setupForSolve(0); double directionAccuracy=findDirectionVector(0); if (directionAccuracy>worstDirectionAccuracy_) { worstDirectionAccuracy_=directionAccuracy; } int phase=0; // predictor, corrector , primal dual // 0 - normal // 1 - second time around i.e. no need for first part // 2 - affine step only // 9 - to exit from while because of error (to go round again) // (9 also used to signal end of while (but then goodMove is true)) int recoveryMode=0; if (!goodMove) { recoveryMode=9; } if (goodMove&&useAffine) { recoveryMode=2; phase=0; } while (recoveryMode!=9) { goodMove=true; if (!recoveryMode) { findStepLength(phase); int nextNumber; //number of complementarity pairs double nextGap=complementarityGap(nextNumber,1); //if (complementarityGap_>1.0e-1&&worstDirectionAccuracy_<1.0e-5) { //wasif (complementarityGap_>1.0e-2*numberComplementarityPairs_) { if (complementarityGap_>1.0e-4*numberComplementarityPairs_) { //std::cout <<"predicted duality gap "<<nextGap<<std::endl; //double part1=nextGap/numberComplementarityPairs_; double part1=nextGap/numberComplementarityPairs_; double part2=nextGap/complementarityGap_; mu_=part1*part2*part2; } else { double phi; if (numberComplementarityPairs_<=500) { phi=pow(numberComplementarityPairs_,2); } else { phi=pow(numberComplementarityPairs_,1.5); if (phi<500.0*500.0) { phi=500.0*500.0; } } mu_=complementarityGap_/phi; //? could use gapO? //if small then should be stopping //if (mu_<1.0e-4/(numberComplementarityPairs_*qqqq)) { //mu_=1.0e-4/(numberComplementarityPairs_*qqqq); //? better to skip corrector? //} } //save information double product=affineProduct(); //can we do corrector step? double xx= complementarityGap_*(beta2-tau) +product; //std::cout<<"gap part "<< //complementarityGap_*(beta2-tau)<<" after adding product = "<<xx; if (xx>0.0) { double saveMu = mu_; double mu2=numberComplementarityPairs_; mu2=xx/mu2; if (mu2>mu_) { //std::cout<<" could increase to "<<mu2<<std::endl; //was mu2=mu2*0.25; mu2=mu2*0.99; if (mu2<mu_) { mu_=mu2; //std::cout<<" changing mu to "<<mu_<<std::endl; } else { //std::cout<<std::endl; } } else { //std::cout<<" should decrease to "<<mu2<<std::endl; mu_=0.5*mu2; std::cout<<" changing mu to "<<mu_<<std::endl; } handler_->message(CLP_BARRIER_MU,messages_) <<saveMu<<mu_ <<CoinMessageEol; } else { //std::cout<<" bad by any standards"<<std::endl; } if (complementarityGap_*(beta2-tau)+product-mu_*numberComplementarityPairs_<0.0) { doCorrector=false; double floatNumber = numberComplementarityPairs_; mu_=complementarityGap_/(floatNumber*floatNumber); //? if small we should be stopping //if (mu_<1.0e-4/(numberComplementarityPairs_*numberComplementarityPairs_)) { //mu_=1.0e-4/(totalVariables*totalVariables); //} handler_->message(CLP_BARRIER_INFO,messages_) <<"no corrector step" <<CoinMessageEol; phase=2; } else { phase=1; } } //set up for next step setupForSolve(phase); double directionAccuracy2=findDirectionVector(phase); if (directionAccuracy2>worstDirectionAccuracy_) { worstDirectionAccuracy_=directionAccuracy2; } double testValue=1.0e2*directionAccuracy; if (1.0e2*projectionTolerance_>testValue) { testValue=1.0e2*projectionTolerance_; } if (primalTolerance()>testValue) { testValue=primalTolerance(); } if (maximumRHSError_>testValue) { testValue=maximumRHSError_; } if (!recoveryMode) { if (directionAccuracy2>testValue&&numberIterations_>=-77) { goodMove=false; useAffine=true;//if bad accuracy doCorrector=false; recoveryMode=9; } } if (goodMove) { findStepLength(phase); if (numberIterations_>=-77) { goodMove=checkGoodMove(doCorrector); } else { goodMove=true; } if (!goodMove) { if (doCorrector) { doCorrector=false; double floatNumber = numberComplementarityPairs_; mu_=complementarityGap_/(floatNumber*floatNumber); handler_->message(CLP_BARRIER_INFO,messages_) <<" no corrector step - original move would be bad" <<CoinMessageEol; phase=2; recoveryMode=1; } else { // if any killed then do zero step and hope for best abort(); } } } //force leave if (goodMove) { recoveryMode=9; } } /* endwhile */ } /* endwhile */ numberFixed=updateSolution(); numberFixedTotal+=numberFixed; } /* endwhile */ if (savePi) { //std::cout<<"Restoring from iteration "<<saveIteration<<std::endl; CoinMemcpyN(savePi,numberRows_,dual_); CoinMemcpyN(savePrimal,numberTotal,solution_); delete [] savePi; delete [] savePrimal; } //recompute slacks // Split out solution CoinZeroN(rowActivity_,numberRows_); CoinMemcpyN(solution_,numberColumns_,columnActivity_); matrix_->times(1.0,columnActivity_,rowActivity_); //unscale objective multiplyAdd(NULL,numberTotal,0.0,cost_,scaleFactor_); multiplyAdd(NULL,numberRows_,0,dual_,scaleFactor_); CoinMemcpyN(cost_,numberColumns_,reducedCost_); matrix_->transposeTimes(-1.0,dual_,reducedCost_); CoinMemcpyN(reducedCost_,numberColumns_,dj_); checkSolution(); handler_->message(CLP_BARRIER_END,messages_) <<sumPrimalInfeasibilities_ <<sumDualInfeasibilities_ <<complementarityGap_ <<objectiveValue() <<CoinMessageEol; //std::cout<<"Absolute primal infeasibility at end "<<sumPrimalInfeasibilities_<<std::endl; //std::cout<<"Absolute dual infeasibility at end "<<sumDualInfeasibilities_<<std::endl; //std::cout<<"Absolute complementarity at end "<<complementarityGap_<<std::endl; //std::cout<<"Primal objective "<<objectiveValue()<<std::endl; //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl; //delete all temporary regions deleteWorkingData(); return problemStatus_; }

---

The documentation for this class was generated from the following files:

• ClpPredictorCorrector.hpp
• ClpPredictorCorrector.cpp

---