00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "BonHeuristicFPump.hpp"
00011 #include "CoinHelperFunctions.hpp"
00012 #include "CbcModel.hpp"
00013
00014 #include "OsiAuxInfo.hpp"
00015
00016 #include "CoinTime.hpp"
00017
00018 #include <fstream>
00019
00020 #include <iomanip>
00021
00022 using namespace std;
00023
00024
00025
00026 namespace Bonmin
00027 {
00028 class score_sorter {
00029 public:
00031 score_sorter(const vector<double>& score):
00032 score_(score) {}
00033
00034 bool operator() (const int x, const int y) const {
00035 return score_[x]>score_[y];
00036 }
00037
00038 private:
00039 const vector<double>& score_;
00040 };
00041
00042
00043 HeuristicFPump::HeuristicFPump()
00044 :
00045 CbcHeuristic(),
00046 setup_(NULL),
00047 objective_norm_(1),
00048 enableAdvanced_(false)
00049 {}
00050
00051 HeuristicFPump::HeuristicFPump(BonminSetup * setup)
00052 :
00053 CbcHeuristic(),
00054 setup_(setup),
00055 objective_norm_(1),
00056 enableAdvanced_(false)
00057 {
00058 Initialize(setup->options());
00059 }
00060
00061 HeuristicFPump::HeuristicFPump(const HeuristicFPump ©)
00062 :
00063 CbcHeuristic(copy),
00064 setup_(copy.setup_),
00065 objective_norm_(copy.objective_norm_),
00066 enableAdvanced_(copy.enableAdvanced_)
00067 {}
00068
00069 HeuristicFPump &
00070 HeuristicFPump::operator=(const HeuristicFPump & rhs)
00071 {
00072 if(this != &rhs) {
00073 CbcHeuristic::operator=(rhs);
00074 setup_ = rhs.setup_;
00075 objective_norm_ = rhs.objective_norm_;
00076 enableAdvanced_ = rhs.enableAdvanced_;
00077 }
00078 return *this;
00079 }
00080
00081 int
00082 HeuristicFPump::solution(double &solutionValue, double *betterSolution)
00083 {
00084 if(model_->getNodeCount() || model_->getCurrentPassNumber() > 1) return 0;
00085
00086 bool integerSolutionAlreadyExists = false;
00087 if(model_->getSolutionCount()) {
00088
00089 integerSolutionAlreadyExists = true;
00090 if(!enableAdvanced_)
00091 return 0;
00092 assert(solutionValue < 1.0e50);
00093 }
00094
00095 const int maxNumberIterations = 200;
00096 const double toleranceObjectiveFP = 1.0e-5;
00097
00098 int returnCode = 0;
00099
00100 OsiTMINLPInterface * nlp = NULL;
00101 if(setup_->getAlgorithm() == B_BB)
00102 nlp = dynamic_cast<OsiTMINLPInterface *>(model_->solver()->clone());
00103 else
00104 nlp = dynamic_cast<OsiTMINLPInterface *>(setup_->nonlinearSolver()->clone());
00105
00106 TMINLP2TNLP* minlp = nlp->problem();
00107
00108
00109 double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
00110 double primalTolerance;
00111 #if 0
00112 OsiSolverInterface * solver = model_->solver();
00113 solver->getDblParam(OsiPrimalTolerance,primalTolerance);
00114 #endif
00115 primalTolerance=1.0e-6;
00116
00117 int numberColumns;
00118 int numberRows;
00119 int nnz_jac_g;
00120 int nnz_h_lag;
00121 TNLP::IndexStyleEnum index_style;
00122 minlp->get_nlp_info(numberColumns, numberRows, nnz_jac_g,
00123 nnz_h_lag, index_style);
00124
00125 const Bonmin::TMINLP::VariableType* variableType = minlp->var_types();
00126 const double* x_sol = minlp->x_sol();
00127 const double* x_l = minlp->x_l();
00128 const double* x_u = minlp->x_u();
00129
00130 #ifdef DEBUG_BON_HEURISTIC_FPUMP
00131 const double* g_sol = minlp->g_sol();
00132 const double* g_l = minlp->g_l();
00133 const double* g_u = minlp->g_u();
00134
00135 for(int i=0; i<numberColumns; i++)
00136 cout<<"x_l["<<i<<"]= "<<x_l[i]<<" "
00137 <<"x_sol["<<i<<"]= "<<x_sol[i]<<" "
00138 <<"x_u["<<i<<"]= "<<x_u[i]<<" "
00139 <<"variableType = "<<variableType[i]<<endl;
00140 for(int i=0; i<numberRows; i++)
00141 cout<<"g_l["<<i<<"]= "<<g_l[i]<<" "
00142 <<"g_sol["<<i<<"]= "<<g_sol[i]<<" "
00143 <<"g_u["<<i<<"]= "<<g_u[i]<<endl;
00144
00145 cout<<"obj_value = "<<minlp->obj_value()<<endl;
00146
00147 cout<<"optimization_status = "<<minlp->optimization_status()<<endl;
00148 #endif
00149
00150
00151
00152 if(minlp->optimization_status() != SUCCESS)
00153 return returnCode;
00154
00155
00156 double* newSolution = new double [numberColumns];
00157 memcpy(newSolution,x_sol,numberColumns*sizeof(double));
00158 double* new_g_sol = new double [numberRows];
00159
00160
00161 vector<int> integerColumns;
00162 int numberFractionalVariables = 0;
00163 for (int iColumn=0;iColumn<numberColumns;iColumn++) {
00164 if (variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
00165 integerColumns.push_back(iColumn);
00166 double value=newSolution[iColumn];
00167 if (fabs(floor(value+0.5)-value)>integerTolerance) {
00168 numberFractionalVariables++;
00169 }
00170 }
00171 }
00172 int numberIntegerColumns = (int) integerColumns.size();
00173
00174
00175 const int numberOldSolutionsStored = 4;
00176 double ** oldSolution = new double * [numberOldSolutionsStored];
00177 for (int j=0;j<numberOldSolutionsStored;j++) {
00178 oldSolution[j]= new double[numberIntegerColumns];
00179 for (int i=0;i<numberIntegerColumns;i++)
00180 oldSolution[j][i]=-COIN_DBL_MAX;
00181 }
00182
00183 RoundingFPump roundObj(minlp);
00184
00185 bool stopDueToAlmostZeroObjective = false;
00186 double* x_bar = new double[numberIntegerColumns];
00187 int* indexes_x_bar = new int[numberIntegerColumns];
00188 double* copy_newSolution = new double[numberColumns];
00189 int iteration = 0;
00190 while(numberFractionalVariables) {
00191 iteration++;
00192 if(iteration > maxNumberIterations) {
00193 break;
00194 }
00195 memcpy(copy_newSolution, newSolution, numberColumns*sizeof(double));
00196 roundObj.round(integerTolerance, primalTolerance, copy_newSolution);
00197 bool flip = true;
00198 for(int iIntCol=0; iIntCol<numberIntegerColumns; iIntCol++) {
00199 int iColumn = integerColumns[iIntCol];
00200 double value=copy_newSolution[iColumn];
00201 #if 0
00202 double value=newSolution[iColumn];
00203 if (fabs(floor(value+0.5)-value)>integerTolerance) {
00204 value = floor(value+0.5);
00205
00206 if(value < x_l[iColumn]-primalTolerance)
00207 value++;
00208 else if(value > x_u[iColumn]+primalTolerance)
00209 value--;
00210 }
00211 #endif
00212 x_bar[iIntCol]=value;
00213 indexes_x_bar[iIntCol]=iColumn;
00214 if(flip &&
00215 fabs(x_bar[iIntCol]-oldSolution[0][iIntCol])>integerTolerance)
00216 flip = false;
00217 }
00218
00219 #ifdef DEBUG_BON_HEURISTIC_FPUMP
00220 cout<<"iteration= "<<iteration<<", flip= "<<flip<<endl;
00221 #endif
00222
00223
00224
00225 if(flip) {
00226 vector<int> sortedIntegerColumns(numberIntegerColumns);
00227 vector<double> score(numberIntegerColumns);
00228 for(int iIntCol=0; iIntCol<numberIntegerColumns; iIntCol++) {
00229 int iColumn = integerColumns[iIntCol];
00230 sortedIntegerColumns[iIntCol] = iIntCol;
00231 double value=newSolution[iColumn];
00232 score[iIntCol] = fabs(value-oldSolution[0][iIntCol]);
00233 }
00234 sort(sortedIntegerColumns.begin(),sortedIntegerColumns.end(),
00235 score_sorter(score));
00236
00237 int maxNumberToMove = 1;
00238 int numberMoved = 0;
00239 for(int i=0; i<numberIntegerColumns; i++) {
00240 int iIntCol = sortedIntegerColumns[i];
00241 if(score[iIntCol] > 0.00) {
00242 int iColumn = integerColumns[iIntCol];
00243 double value=newSolution[iColumn];
00244 if(value-oldSolution[0][iIntCol]>0.0)
00245 value = oldSolution[0][iIntCol]+1.0;
00246 else
00247 value = oldSolution[0][iIntCol]-1.0;
00248
00249 if(value < x_l[iColumn]-primalTolerance)
00250 value++;
00251 else if(value > x_u[iColumn]+primalTolerance)
00252 value--;
00253 assert(fabs(floor(value+0.5)-value)<=integerTolerance);
00254 x_bar[iIntCol]=value;
00255 numberMoved++;
00256 } else
00257 break;
00258 if(numberMoved >= maxNumberToMove)
00259 break;
00260 }
00261
00262
00263
00264
00265 bool matched;
00266 for (int k = numberOldSolutionsStored-1; k > 0; k--) {
00267 double * b = oldSolution[k];
00268 matched = true;
00269 for(int iIntCol=0; iIntCol<numberIntegerColumns; iIntCol++) {
00270 if (fabs(x_bar[iIntCol]-b[iIntCol])>integerTolerance) {
00271 matched=false;
00272 break;
00273 }
00274 }
00275 if (matched) break;
00276 }
00277
00278 #ifdef DEBUG_BON_HEURISTIC_FPUMP
00279 cout<<"matched= "<<matched<<endl;
00280 #endif
00281
00282 if (matched) {
00283
00284 for(int iIntCol=0; iIntCol<numberIntegerColumns; iIntCol++) {
00285 int iColumn = integerColumns[iIntCol];
00286 double value=newSolution[iColumn];
00287 double random = max(0.0,CoinDrand48()-0.3);
00288 double difference = fabs(value-oldSolution[0][iIntCol]);
00289 if(difference+random>0.5) {
00290 if(value-oldSolution[0][iIntCol]>0.0)
00291 value = oldSolution[0][iIntCol]+1.0;
00292 else
00293 value = oldSolution[0][iIntCol]-1.0;
00294
00295 if(value < x_l[iColumn]-primalTolerance)
00296 value++;
00297 else if(value > x_u[iColumn]+primalTolerance)
00298 value--;
00299 assert(fabs(floor(value+0.5)-value)<=integerTolerance);
00300 } else {
00301
00302 value = oldSolution[0][iIntCol];
00303 }
00304 x_bar[iIntCol]=value;
00305 }
00306 }
00307 }
00308
00309 for (int j=numberOldSolutionsStored-1;j>0;j--) {
00310 for (int i = 0; i < numberIntegerColumns; i++)
00311 oldSolution[j][i]=oldSolution[j-1][i];
00312 }
00313 for (int j = 0; j < numberIntegerColumns; j++)
00314 oldSolution[0][j] = x_bar[j];
00315
00316
00317
00318 double obj_nlp;
00319 if(integerSolutionAlreadyExists)
00320
00321 obj_nlp = nlp->solveFeasibilityProblem(numberIntegerColumns,
00322 x_bar,indexes_x_bar,
00323 objective_norm_, solutionValue);
00324 else
00325 obj_nlp = nlp->solveFeasibilityProblem(numberIntegerColumns,
00326 x_bar,indexes_x_bar,
00327 1,0,objective_norm_);
00328
00329
00330 #ifdef DEBUG_BON_HEURISTIC_FPUMP
00331 cout<<"obj_nlp= "<<obj_nlp<<endl;
00332 #endif
00333
00334 memcpy(newSolution,x_sol,numberColumns*sizeof(double));
00335
00336 if(obj_nlp < toleranceObjectiveFP) {
00337 stopDueToAlmostZeroObjective = true;
00338 break;
00339 }
00340
00341
00342 numberFractionalVariables = 0;
00343 for(int iIntCol=0; iIntCol<numberIntegerColumns; iIntCol++) {
00344 int iColumn = integerColumns[iIntCol];
00345 double value=newSolution[iColumn];
00346 if (fabs(floor(value+0.5)-value)>integerTolerance)
00347 numberFractionalVariables++;
00348 }
00349
00350 }
00351
00352 for (int j=0;j<numberOldSolutionsStored;j++)
00353 delete [] oldSolution[j];
00354 delete [] oldSolution;
00355 delete [] x_bar;
00356 delete [] indexes_x_bar;
00357
00358
00359
00360 for(int iIntCol=0; iIntCol<numberIntegerColumns; iIntCol++) {
00361 int iColumn = integerColumns[iIntCol];
00362 double value=floor(newSolution[iColumn]+0.5);
00363 minlp->SetVariableUpperBound(iColumn, floor(value));
00364 minlp->SetVariableLowerBound(iColumn, ceil(value));
00365 }
00366 nlp->initialSolve();
00367 bool feasible = true;
00368 if(minlp->optimization_status() != SUCCESS) {
00369 feasible = false;
00370
00371
00372 }
00373 memcpy(newSolution,x_sol,numberColumns*sizeof(double));
00374
00375 if(feasible) {
00376 double newSolutionValue;
00377 minlp->eval_f(numberColumns, newSolution, true, newSolutionValue);
00378 if(newSolutionValue < solutionValue) {
00379 memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
00380 solutionValue = newSolutionValue;
00381 returnCode = 1;
00382 }
00383 }
00384
00385 #ifdef DEBUG_BON_HEURISTIC_FPUMP
00386 cout<<"returnCode= "<<returnCode<<endl;
00387 #endif
00388
00389 #if 0
00390 delete [] indexRow;
00391 delete [] indexCol;
00392 delete [] row;
00393 delete [] columnStart;
00394 delete [] columnLength;
00395 #endif
00396 delete [] newSolution;
00397 delete [] new_g_sol;
00398 delete [] copy_newSolution;
00399
00400 return returnCode;
00401 }
00402
00403 void
00404 HeuristicFPump::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions){
00405 roptions->SetRegisteringCategory("MINLP Heuristics", RegisteredOptions::BonminCategory);
00406 roptions->AddBoundedIntegerOption("feasibility_pump_objective_norm","Norm of feasibility pump objective function",
00407 1, 2, 1,"");
00408 roptions->setOptionExtraInfo("feasibility_pump_objective_norm", 63);
00409 roptions->AddStringOption2("heuristic_feasibility_pump", "whether the heuristic feasibility pump should be used",
00410 "no", "no", "don't use it", "yes", "use it", "");
00411 roptions->setOptionExtraInfo("heuristic_feasibility_pump", 63);
00412
00413 roptions->SetRegisteringCategory("Test Heuristics", RegisteredOptions::UndocumentedCategory);
00414 roptions->AddStringOption2("unstable_fp","use at your own risks",
00415 "no",
00416 "no", "",
00417 "yes", "","");
00418 roptions->setOptionExtraInfo("unstable_fp", 63);
00419 }
00420
00421 void
00422 HeuristicFPump::Initialize(Ipopt::SmartPtr<Bonmin::OptionsList> options){
00423 options->GetIntegerValue("feasibility_pump_objective_norm", objective_norm_, "bonmin.");
00424 options->GetEnumValue("unstable_fp", enableAdvanced_, "bonmin.");
00425 }
00426
00427 RoundingFPump::RoundingFPump(TMINLP2TNLP* minlp)
00428 :
00429 minlp_(minlp)
00430 {
00431 gutsOfConstructor();
00432 }
00433
00434 RoundingFPump::~RoundingFPump()
00435 {
00436 delete [] col_and_jac_g_;
00437 }
00438
00439 void
00440 RoundingFPump::gutsOfConstructor()
00441 {
00442
00443 int nnz_jac_g;
00444 int nnz_h_lag;
00445 TNLP::IndexStyleEnum index_style;
00446 minlp_->get_nlp_info(numberColumns_, numberRows_, nnz_jac_g,
00447 nnz_h_lag, index_style);
00448
00449 const double* x_sol = minlp_->x_sol();
00450
00451
00452
00453
00454 int* indexRow = new int[nnz_jac_g];
00455 int* indexCol = new int[nnz_jac_g];
00456 minlp_->eval_jac_g(numberColumns_, x_sol, false,
00457 numberRows_, nnz_jac_g,
00458 indexRow, indexCol, 0);
00459
00460
00461 double* jac_g = new double [nnz_jac_g];
00462 double* zero_sol = new double [numberColumns_];
00463 minlp_->get_starting_point(numberColumns_, 1, zero_sol, 0, NULL, NULL, numberRows_, 0, NULL);
00464
00465 minlp_->eval_jac_g(numberColumns_, zero_sol, true,
00466 numberRows_, nnz_jac_g,
00467 0, 0, jac_g);
00468
00469 col_and_jac_g_ = new vector<pair<int, int> >[numberRows_];
00470
00471 int indexCorrection = (index_style == TNLP::C_STYLE) ? 0 : 1;
00472 for(int i=0; i<nnz_jac_g; i++) {
00473 int thisIndexRow = indexRow[i]-indexCorrection;
00474 int thisIndexCol = indexCol[i]-indexCorrection;
00475 pair<int, int> value(thisIndexCol, static_cast<int>(jac_g[i]));
00476 col_and_jac_g_[thisIndexRow].push_back(value);
00477 }
00478
00479 delete [] indexRow;
00480 delete [] indexCol;
00481 delete [] jac_g;
00482 delete [] zero_sol;
00483 }
00484
00485 void
00486 RoundingFPump::round(const double integerTolerance,
00487 const double primalTolerance,
00488 double* solution)
00489 {
00490 const Bonmin::TMINLP::VariableType* variableType = minlp_->var_types();
00491 const double* x_l = minlp_->x_l();
00492 const double* x_u = minlp_->x_u();
00493 const double* g_l = minlp_->g_l();
00494 const double* g_u = minlp_->g_u();
00495
00496
00497 for(int iRow=0; iRow<numberRows_; iRow++) {
00498 if(g_l[iRow] == 1.0 && g_u[iRow] == 1.0) {
00499 bool sosConstraint = true;
00500 double weightedSum = 0.0;
00501 int counter = 1;
00502 vector<pair<int, int> > jac_g = col_and_jac_g_[iRow];
00503 for (unsigned int j=0; j<jac_g.size(); j++) {
00504 int iColumn = jac_g[j].first;
00505 if (solution[iColumn]>=1.0-integerTolerance ||
00506 jac_g[j].second != 1.0 ||
00507 variableType[iColumn] == Bonmin::TMINLP::CONTINUOUS) {
00508 sosConstraint = false;
00509 break;
00510 }
00511 else {
00512 weightedSum += counter * solution[iColumn];
00513 counter++;
00514 }
00515 }
00516 #ifdef DEBUG_BON_HEURISTIC_FPUMP
00517 if(sosConstraint) {
00518 cout<<"weightedSum= "<<weightedSum
00519 <<", numberColumns_= "<<numberColumns_<<endl;
00520 }
00521 #endif
00522
00523 if(sosConstraint) {
00524 double fl = floor(weightedSum + 0.5);
00525 int indexColumnSelected = static_cast<int>(fl) - 1;
00526 if(indexColumnSelected < 0){
00527 continue;
00528 }
00529 assert(indexColumnSelected < jac_g.size());
00530 for (unsigned int j=0; j<jac_g.size(); j++) {
00531 int iColumn = jac_g[j].first;
00532 if(j == indexColumnSelected)
00533 solution[iColumn] = 1.0;
00534 else
00535 solution[iColumn] = 0.0;
00536 }
00537 }
00538 }
00539 }
00540
00541 for(int iColumn=0; iColumn<numberColumns_; iColumn++) {
00542 if(variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
00543 double value=solution[iColumn];
00544 if (fabs(floor(value+0.5)-value)>integerTolerance) {
00545 value = floor(value+0.5);
00546
00547 if(value < x_l[iColumn]-primalTolerance)
00548 value++;
00549 else if(value > x_u[iColumn]+primalTolerance)
00550 value--;
00551 solution[iColumn] = value;
00552 }
00553 }
00554 }
00555 }
00556 }