00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "BonMilpRounding.hpp"
00011 #include "CoinHelperFunctions.hpp"
00012 #include "CbcModel.hpp"
00013 #include "BonSubMipSolver.hpp"
00014
00015 #include "CoinTime.hpp"
00016
00017 #include <fstream>
00018
00019 #include <iomanip>
00020
00021 #include "CoinHelperFunctions.hpp"
00022 #include "OsiClpSolverInterface.hpp"
00023
00024
00025
00026 using namespace std;
00027
00028 namespace Bonmin
00029 {
00030 MilpRounding::MilpRounding(BonminSetup * setup)
00031 :
00032 CbcHeuristic(),
00033 setup_(setup),
00034 howOften_(20),
00035 mip_(NULL)
00036 {
00037 Initialize(setup);
00038 }
00039
00040 void
00041 MilpRounding::Initialize(BonminSetup * b){
00042 delete mip_;
00043 mip_ = new SubMipSolver (*b, b->prefix());
00044
00045 }
00046
00047 MilpRounding::MilpRounding(const MilpRounding ©)
00048 :
00049 CbcHeuristic(copy),
00050 setup_(copy.setup_),
00051 howOften_(copy.howOften_),
00052 mip_(new SubMipSolver(*copy.mip_))
00053 {
00054 }
00055
00056 MilpRounding &
00057 MilpRounding::operator=(const MilpRounding & rhs)
00058 {
00059 if(this != &rhs) {
00060 CbcHeuristic::operator=(rhs);
00061 setup_ = rhs.setup_;
00062 howOften_ = rhs.howOften_;
00063 delete mip_;
00064 if(rhs.mip_)
00065 mip_ = new SubMipSolver(*rhs.mip_);
00066 }
00067 return *this;
00068 }
00069
00070 MilpRounding::~MilpRounding(){
00071 delete mip_;
00072 }
00073
00074 struct MatComp{
00075 const int * iRow;
00076 const int * jCol;
00078 bool operator()(int i,int j){
00079 return (jCol[i] < jCol[j]) || (jCol[i] == jCol[j] && iRow[i] < iRow[j]);
00080 }
00081 };
00082
00083
00084 int
00085 MilpRounding::solution(double &solutionValue, double *betterSolution)
00086 {
00087 printf("Entering heuristic\n");
00088 if(model_->getCurrentPassNumber() > 1) return 0;
00089 if ((model_->getNodeCount()%howOften_)!=0||model_->getCurrentPassNumber()>1)
00090 return 0;
00091
00092 int returnCode = 0;
00093
00094 OsiTMINLPInterface * nlp = NULL;
00095 if(setup_->getAlgorithm() == B_BB)
00096 nlp = dynamic_cast<OsiTMINLPInterface *>(model_->solver()->clone());
00097 else
00098 nlp = dynamic_cast<OsiTMINLPInterface *>(setup_->nonlinearSolver()->clone());
00099
00100 TMINLP2TNLP* minlp = nlp->problem();
00101
00102
00103 double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
00104
00105
00106 int n;
00107 int m;
00108 int nnz_jac_g;
00109 int nnz_h_lag;
00110 Ipopt::TNLP::IndexStyleEnum index_style;
00111 minlp->get_nlp_info(n, m, nnz_jac_g,
00112 nnz_h_lag, index_style);
00113
00114 const Bonmin::TMINLP::VariableType* variableType = minlp->var_types();
00115 const double* x_sol = minlp->x_sol();
00116 const double* g_l = minlp->g_l();
00117 const double* g_u = minlp->g_u();
00118
00119 const double * colsol = model_->solver()->getColSolution();
00120
00121
00122
00123 TMINLP* tminlp = nlp->model();
00124 vector<Ipopt::TNLP::LinearityType> c_lin(m);
00125 tminlp->get_constraints_linearity(m, c_lin());
00126 vector<int> c_idx(m);
00127 int n_lin = 0;
00128 for (int i=0;i<m;i++) {
00129 if (c_lin[i]==Ipopt::TNLP::LINEAR)
00130 c_idx[i] = n_lin++;
00131 else
00132 c_idx[i] = -1;
00133 }
00134
00135
00136
00137 vector<int> indexRow(nnz_jac_g);
00138 vector<int> indexCol(nnz_jac_g);
00139 minlp->eval_jac_g(n, x_sol, false,
00140 m, nnz_jac_g,
00141 indexRow(), indexCol(), 0);
00142
00143
00144 vector<double> jac_g(nnz_jac_g);
00145 minlp->eval_jac_g(n, x_sol, false,
00146 m, nnz_jac_g,
00147 NULL, NULL, jac_g());
00148
00149
00150 vector<int> sortedIndex(nnz_jac_g);
00151 CoinIotaN(sortedIndex(), nnz_jac_g, 0);
00152 MatComp c;
00153 c.iRow = indexRow();
00154 c.jCol = indexCol();
00155 std::sort(sortedIndex.begin(), sortedIndex.end(), c);
00156
00157 vector<int> row (nnz_jac_g);
00158 vector<double> value (nnz_jac_g);
00159 vector<int> columnStart(n,0);
00160 vector<int> columnLength(n,0);
00161 int indexCorrection = (index_style == Ipopt::TNLP::C_STYLE) ? 0 : 1;
00162 int iniCol = -1;
00163 int nnz = 0;
00164 for(int i=0; i<nnz_jac_g; i++) {
00165 int thisIndexCol = indexCol[sortedIndex[i]]-indexCorrection;
00166 int thisIndexRow = c_idx[indexRow[sortedIndex[i]] - indexCorrection];
00167 if(thisIndexCol != iniCol) {
00168 iniCol = thisIndexCol;
00169 columnStart[thisIndexCol] = nnz;
00170 columnLength[thisIndexCol] = 0;
00171 }
00172 if(thisIndexRow == -1) continue;
00173 columnLength[thisIndexCol]++;
00174 row[nnz] = thisIndexRow;
00175 value[nnz] = jac_g[i];
00176 nnz++;
00177 }
00178
00179
00180 vector<double> newRowLower(n_lin);
00181 vector<double> newRowUpper(n_lin);
00182 for(int i = 0 ; i < m ; i++){
00183 if(c_idx[i] == -1) continue;
00184 newRowLower[c_idx[i]] = g_l[i];
00185 newRowUpper[c_idx[i]] = g_u[i];
00186 }
00187
00188
00189 vector<double> newSolution(n);
00190 std::copy(x_sol, x_sol + n, newSolution.begin());
00191
00192
00193 CoinPackedMatrix matrix(true,n_lin,n, nnz, value(), row(), columnStart(), columnLength());
00194
00195
00196
00197
00198 double beta = 1;
00199 vector<double> objective(n);
00200 vector<int> idxIntegers;
00201 idxIntegers.reserve(n);
00202 for(int i = 0 ; i < n ; i++){
00203 if(variableType[i] != Bonmin::TMINLP::CONTINUOUS){
00204 idxIntegers.push_back(i);
00205 objective[i] = beta*(1 - 2*colsol[i]);
00206 }
00207 }
00208
00209 #if 0
00210
00211 const double * duals = nlp->getRowPrice() + 2 *n;
00212 vector<double> grad(n, 0);
00213 vector<int> indices(n, 0);
00214 tminlp->eval_grad_f(n, x_sol, false, grad());
00215 for(int i = 0 ; i < m ; i++){
00216 if(c_lin[i] == Ipopt::TNLP::LINEAR) continue;
00217 int nnz;
00218 tminlp->eval_grad_gi(n, x_sol, false, i, nnz, indices(), NULL);
00219 tminlp->eval_grad_gi(n, x_sol, false, i, nnz, NULL, grad());
00220 for(int k = 0 ; k < nnz ; k++){
00221 objective[indices[k]] += alpha *duals[i] * grad[k];
00222 }
00223 }
00224 for(int i = 0 ; i < n ; i++){
00225 if(variableType[i] != Bonmin::TMINLP::CONTINUOUS)
00226 objective[i] += alpha * grad[i];
00227
00228 else objective[i] = 0;
00229 }
00230 std::copy(objective.begin(), objective.end(), std::ostream_iterator<double>(std::cout, " "));
00231 std::cout<<std::endl;
00232 #endif
00233
00234
00235 OsiSolverInterface *si = mip_->solver();
00236 assert(si != NULL);
00237 CoinMessageHandler * handler = model_->messageHandler()->clone();
00238 si->passInMessageHandler(handler);
00239 si->messageHandler()->setLogLevel(1);
00240
00241 si->loadProblem(matrix, model_->solver()->getColLower(), model_->solver()->getColUpper(), objective(),
00242 newRowLower(), newRowUpper());
00243 si->setInteger(idxIntegers(), static_cast<int>(idxIntegers.size()));
00244 si->applyCuts(noGoods);
00245 printf("Done creating mip, start solving\n");
00246
00247 bool hasFractionnal = true;
00248 while(hasFractionnal){
00249 mip_->optimize(DBL_MAX, 0, 60);
00250 printf("Done solving\n");
00251 hasFractionnal = false;
00252 #if 0
00253 bool feasible = false;
00254 if(mip_->getLastSolution()) {
00255 const double* solution = mip_->getLastSolution();
00256 std::copy(solution, solution + n, newSolution.begin());
00257 feasible = true;
00258
00259 }
00260
00261 if(feasible) {
00262
00263
00264 CoinPackedVector v;
00265 double lb = 1;
00266 for (int iColumn=0;iColumn<n;iColumn++) {
00267 if (variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
00268 double value=newSolution[iColumn];
00269 if (fabs(floor(value+0.5)-value)>integerTolerance) {
00270 #ifdef DEBUG_BON_HEURISTIC_DIVE_MIP
00271 cout<<"It should be infeasible because: "<<endl;
00272 cout<<"variable "<<iColumn<<" is not integer"<<endl;
00273 #endif
00274 feasible = false;
00275 break;
00276 }
00277 else {
00278 value=floor(newSolution[iColumn]+0.5);
00279 if(value > 0.5){
00280 v.insert(iColumn, -1);
00281 lb -= value;
00282 }
00283 minlp->SetVariableUpperBound(iColumn, value);
00284 minlp->SetVariableLowerBound(iColumn, value);
00285 }
00286 }
00287 }
00288 }
00289 #endif
00290 }
00291 bool feasible = false;
00292 if(mip_->getLastSolution()) {
00293 const double* solution = mip_->getLastSolution();
00294 std::copy(solution, solution + n, newSolution.begin());
00295 feasible = true;
00296
00297 delete handler;
00298 }
00299
00300
00301 if(feasible) {
00302
00303
00304 printf("We found a solution");
00305 CoinPackedVector v;
00306 double lb = 1;
00307 for (int iColumn=0;iColumn<n;iColumn++) {
00308 if (variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
00309 double value=newSolution[iColumn];
00310 if (fabs(floor(value+0.5)-value)>integerTolerance) {
00311 #ifdef DEBUG_BON_HEURISTIC_DIVE_MIP
00312 cout<<"It should be infeasible because: "<<endl;
00313 cout<<"variable "<<iColumn<<" is not integer"<<endl;
00314 #endif
00315 feasible = false;
00316 break;
00317 }
00318 else {
00319 value=floor(newSolution[iColumn]+0.5);
00320 if(value > 0.5){
00321 v.insert(iColumn, -1);
00322 lb -= value;
00323 }
00324 minlp->SetVariableUpperBound(iColumn, value);
00325 minlp->SetVariableLowerBound(iColumn, value);
00326 }
00327 }
00328 }
00329 OsiRowCut c;
00330 c.setRow(v);
00331 c.setLb(lb);
00332 c.setUb(DBL_MAX);
00333 noGoods.insert(c);
00334 if(feasible) {
00335 nlp->initialSolve();
00336 if(minlp->optimization_status() != Ipopt::SUCCESS) {
00337 printf("Infeasible for NLP");
00338 feasible = false;
00339 }
00340 std::copy(x_sol,x_sol+n, newSolution.begin());
00341 }
00342 }
00343 else {printf("No solution found\n");}
00344 if(feasible) {
00345 double newSolutionValue;
00346 minlp->eval_f(n, newSolution(), true, newSolutionValue);
00347 if(newSolutionValue < solutionValue) {
00348 std::copy(newSolution.begin(), newSolution.end(), betterSolution);
00349 solutionValue = newSolutionValue;
00350 returnCode = 1;
00351 }
00352 }
00353
00354 delete nlp;
00355
00356 #ifdef DEBUG_BON_HEURISTIC_DIVE_MIP
00357 std::cout<<"DiveMIP returnCode = "<<returnCode<<std::endl;
00358 #endif
00359
00360 return returnCode;
00361 }
00362 void
00363 MilpRounding::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions){
00364 roptions->SetRegisteringCategory("Undocumented Heuristics", RegisteredOptions::UndocumentedCategory);
00365 roptions->AddStringOption2(
00366 "MILP_rounding_heuristic",
00367 "if yes runs the heuristic",
00368 "no",
00369 "no", "don't run it",
00370 "yes", "runs the heuristic",
00371 "");
00372 }
00373 }