FLOPC++
|
00001 // ******************** FlopCpp ********************************************** 00002 // File: MP_model.cpp 00003 // $Id$ 00004 // Author: Tim Helge Hultberg (thh@mat.ua.pt) 00005 // Copyright (C) 2003 Tim Helge Hultberg 00006 // All Rights Reserved. 00007 //**************************************************************************** 00008 00009 #include <iostream> 00010 #include <sstream> 00011 using std::cout; 00012 using std::endl; 00013 #include <algorithm> 00014 00015 #include <CoinPackedMatrix.hpp> 00016 #include <OsiSolverInterface.hpp> 00017 #include "MP_model.hpp" 00018 #include "MP_variable.hpp" 00019 #include "MP_constraint.hpp" 00020 #include <CoinTime.hpp> 00021 00022 using namespace flopc; 00023 00024 MP_model& MP_model::default_model = *new MP_model(0); 00025 MP_model* MP_model::current_model = &MP_model::default_model; 00026 MP_model &MP_model::getDefaultModel() { return default_model;} 00027 MP_model *MP_model::getCurrentModel() { return current_model;} 00028 00029 void NormalMessenger::statistics(int bm, int m, int bn, int n, int nz) { 00030 cout<<"FlopCpp: Number of constraint blocks: " <<bm<<endl; 00031 cout<<"FlopCpp: Number of individual constraints: " <<m<<endl; 00032 cout<<"FlopCpp: Number of variable blocks: " <<bn<<endl; 00033 cout<<"FlopCpp: Number of individual variables: " <<n<<endl; 00034 cout<<"FlopCpp: Number of non-zeroes (including rhs): " <<nz<<endl; 00035 } 00036 00037 void NormalMessenger::generationTime(double t) { 00038 cout<<"FlopCpp: Generation time: "<<t<<endl; 00039 } 00040 00041 void VerboseMessenger::constraintDebug(string name, const vector<MP::Coef>& cfs) { 00042 cout<<"FlopCpp: Constraint "<<name<<endl; 00043 for (unsigned int j=0; j<cfs.size(); j++) { 00044 int col=cfs[j].col; 00045 int row=cfs[j].row; 00046 double elm=cfs[j].val; 00047 int stage=cfs[j].stage; 00048 cout<<row<<" "<<col<<" "<<elm<<" "<<stage<<endl; 00049 } 00050 } 00051 00052 void VerboseMessenger::objectiveDebug(const vector<MP::Coef>& cfs) { 00053 cout<<"Objective "<<endl; 00054 for (unsigned int j=0; j<cfs.size(); j++) { 00055 int col=cfs[j].col; 00056 int row=cfs[j].row; 00057 double elm=cfs[j].val; 00058 cout<<row<<" "<<col<<" "<<elm<<endl; 00059 } 00060 } 00061 00062 MP_model::MP_model(OsiSolverInterface* s, Messenger* m) : 00063 messenger(m), Objective(0), Solver(s), 00064 m(0), n(0), nz(0), bl(0), 00065 mSolverState(((s==0)?(MP_model::DETACHED):(MP_model::SOLVER_ONLY))) { 00066 MP_model::current_model = this; 00067 } 00068 00069 MP_model::~MP_model() { 00070 delete messenger; 00071 if (Solver!=0) { 00072 delete Solver; 00073 } 00074 } 00075 00076 00077 MP_model& MP_model::add(MP_constraint& constraint) { 00078 Constraints.insert(&constraint); 00079 return *this; 00080 } 00081 00082 void MP_model::add(MP_constraint* constraint) { 00083 constraint->M = this; 00084 if (constraint->left.isDefined() && constraint->right.isDefined()) { 00085 constraint->offset = m; 00086 m += constraint->size(); 00087 } 00088 } 00089 00090 double MP_model::getInfinity() const { 00091 if (Solver==0) { 00092 return 9.9e+32; 00093 } else { 00094 return Solver->getInfinity(); 00095 } 00096 } 00097 00098 void MP_model::add(MP_variable* v) { 00099 v->M = this; 00100 v->offset = n; 00101 n += v->size(); 00102 } 00103 00104 void MP_model::addRow(const Constraint& constraint) { 00105 vector<MP::Coef> cfs; 00106 vector<Constant> v; 00107 MP::GenerateFunctor f(0,cfs); 00108 constraint->left->generate(MP_domain::getEmpty(),v,f,1.0); 00109 constraint->right->generate(MP_domain::getEmpty(),v,f,-1.0); 00110 CoinPackedVector newRow; 00111 double rhs = 0.0; 00112 for (unsigned int j=0; j<cfs.size(); j++) { 00113 int col=cfs[j].col; 00114 //int row=cfs[j].row; 00115 double elm=cfs[j].val; 00116 //cout<<row<<" "<<col<<" "<<elm<<endl; 00117 if (col>=0) { 00118 newRow.insert(col,elm); 00119 } else if (col==-1) { 00120 rhs = elm; 00121 } 00122 } 00123 // NB no "assembly" of added row yet. Will be done.... 00124 00125 double local_bl = -rhs; 00126 double local_bu = -rhs; 00127 00128 double inf = Solver->getInfinity(); 00129 switch (constraint->sense) { 00130 case LE: 00131 local_bl = - inf; 00132 break; 00133 case GE: 00134 local_bu = inf; 00135 break; 00136 case EQ: 00137 // Nothing to do 00138 break; 00139 } 00140 00141 Solver->addRow(newRow,local_bl,local_bu); 00142 } 00143 00144 void MP_model::setObjective(const MP_expression& o) { 00145 Objective = o; 00146 } 00147 00148 void MP_model::minimize_max(MP_set &s, const MP_expression &obj) { 00149 MP_variable v; 00150 MP_constraint constraint(s); 00151 add(constraint); 00152 constraint(s) = v() >= obj; 00153 minimize(v()); 00154 } 00155 00156 00157 void MP_model::assemble(vector<MP::Coef>& v, vector<MP::Coef>& av) { 00158 std::sort(v.begin(),v.end(),MP::CoefLess()); 00159 int c,r,s; 00160 double val; 00161 std::vector<MP::Coef>::const_iterator i = v.begin(); 00162 while (i!=v.end()) { 00163 c = i->col; 00164 r = i->row; 00165 val = i->val; 00166 s = i->stage; 00167 i++; 00168 while (i!=v.end() && c==i->col && r==i->row) { 00169 val += i->val; 00170 if (i->stage>s) { 00171 s = i->stage; 00172 } 00173 i++; 00174 } 00175 av.push_back(MP::Coef(c,r,val,s)); 00176 } 00177 } 00178 00179 void MP_model::maximize() { 00180 if (Solver!=0) { 00181 attach(Solver); 00182 solve(MP_model::MAXIMIZE); 00183 } else { 00184 cout<<"no solver specified"<<endl; 00185 } 00186 } 00187 00188 void MP_model::maximize(const MP_expression &obj) { 00189 if (Solver!=0) { 00190 Objective = obj; 00191 attach(Solver); 00192 solve(MP_model::MAXIMIZE); 00193 } else { 00194 cout<<"no solver specified"<<endl; 00195 } 00196 } 00197 00198 void MP_model::minimize() { 00199 if (Solver!=0) { 00200 attach(Solver); 00201 solve(MP_model::MINIMIZE); 00202 } else { 00203 cout<<"no solver specified"<<endl; 00204 } 00205 } 00206 00207 void MP_model::minimize(const MP_expression &obj) { 00208 if (Solver!=0) { 00209 Objective = obj; 00210 attach(Solver); 00211 solve(MP_model::MINIMIZE); 00212 } else { 00213 cout<<"no solver specified"<<endl; 00214 } 00215 } 00216 00217 void MP_model::attach(OsiSolverInterface *_solver) { 00218 if (_solver == 0) { 00219 if (Solver == 0) { 00220 mSolverState = MP_model::DETACHED; 00221 return; 00222 } 00223 } else { // use pre-attached solver. 00224 if(Solver && Solver!=_solver) { 00225 detach(); 00226 } 00227 Solver=_solver; 00228 } 00229 double time = CoinCpuTime(); 00230 m=0; 00231 n=0; 00232 vector<MP::Coef> coefs; 00233 vector<MP::Coef> cfs; 00234 00235 typedef std::set<MP_variable* >::iterator varIt; 00236 typedef std::set<MP_constraint* >::iterator conIt; 00237 00238 Objective->insertVariables(Variables); 00239 for (conIt i=Constraints.begin(); i!=Constraints.end(); i++) { 00240 add(*i); 00241 (*i)->insertVariables(Variables); 00242 } 00243 for (varIt j=Variables.begin(); j!=Variables.end(); j++) { 00244 add(*j); 00245 } 00246 00247 // Generate coefficient matrix and right hand side 00248 for (conIt i=Constraints.begin(); i!=Constraints.end(); i++) { 00249 (*i)->coefficients(cfs); 00250 messenger->constraintDebug((*i)->getName(),cfs); 00251 assemble(cfs,coefs); 00252 cfs.erase(cfs.begin(),cfs.end()); 00253 } 00254 nz = coefs.size(); 00255 00256 messenger->statistics(Constraints.size(),m,Variables.size(),n,nz); 00257 00258 if (nz>0) { 00259 Elm = new double[nz]; 00260 Rnr = new int[nz]; 00261 } 00262 Cst = new int[n+2]; 00263 Clg = new int[n+1]; 00264 if (n>0) { 00265 l = new double[n]; 00266 u = new double[n]; 00267 c = new double[n]; 00268 } 00269 if (m>0) { 00270 bl = new double[m]; 00271 bu = new double[m]; 00272 } 00273 const double inf = Solver->getInfinity(); 00274 00275 for (int j=0; j<n; j++) { 00276 Clg[j] = 0; 00277 } 00278 Clg[n] = 0; 00279 00280 // Treat right hand side as n'th column 00281 for (int j=0; j<=n; j++) { 00282 Clg[j] = 0; 00283 } 00284 for (int i=0; i<nz; i++) { 00285 int col = coefs[i].col; 00286 if (col == -1) { 00287 col = n; 00288 } 00289 Clg[col]++; 00290 } 00291 Cst[0]=0; 00292 for (int j=0; j<=n; j++) { 00293 Cst[j+1]=Cst[j]+Clg[j]; 00294 } 00295 for (int i=0; i<=n; i++) { 00296 Clg[i]=0; 00297 } 00298 for (int i=0; i<nz; i++) { 00299 int col = coefs[i].col; 00300 if (col==-1) { 00301 col = n; 00302 } 00303 int row = coefs[i].row; 00304 double elm = coefs[i].val; 00305 Elm[Cst[col]+Clg[col]] = elm; 00306 Rnr[Cst[col]+Clg[col]] = row; 00307 Clg[col]++; 00308 } 00309 00310 // Row bounds 00311 for (int i=0; i<m; i++) { 00312 bl[i] = 0; 00313 bu[i] = 0; 00314 } 00315 for (int j=Cst[n]; j<Cst[n+1]; j++) { 00316 bl[Rnr[j]] = -Elm[j]; 00317 bu[Rnr[j]] = -Elm[j]; 00318 } 00319 00320 for (conIt i=Constraints.begin(); i!=Constraints.end(); i++) { 00321 if ((*i)->left.isDefined() && (*i)->right.isDefined() ) { 00322 int begin = (*i)->offset; 00323 int end = (*i)->offset+(*i)->size(); 00324 switch ((*i)->sense) { 00325 case LE: 00326 for (int k=begin; k<end; k++) { 00327 bl[k] = - inf; 00328 } 00329 break; 00330 case GE: 00331 for (int k=begin; k<end; k++) { 00332 bu[k] = inf; 00333 } 00334 break; 00335 case EQ: 00336 // Nothing to do 00337 break; 00338 } 00339 } 00340 } 00341 00342 // Generate objective function coefficients 00343 vector<Constant> v; 00344 MP::GenerateFunctor f(0,cfs); 00345 coefs.erase(coefs.begin(),coefs.end()); 00346 Objective->generate(MP_domain::getEmpty(), v, f, 1.0); 00347 00348 messenger->objectiveDebug(cfs); 00349 assemble(cfs,coefs); 00350 00351 for (int j=0; j<n; j++) { 00352 c[j] = 0.0; 00353 } 00354 for (size_t i=0; i<coefs.size(); i++) { 00355 int col = coefs[i].col; 00356 double elm = coefs[i].val; 00357 c[col] = elm; 00358 } 00359 00360 // Column bounds 00361 for (int j=0; j<n; j++) { 00362 l[j] = 0.0; 00363 u[j] = inf; 00364 } 00365 00366 for (varIt i=Variables.begin(); i!=Variables.end(); i++) { 00367 for (int k=0; k<(*i)->size(); k++) { 00368 l[(*i)->offset+k] = (*i)->lowerLimit.v[k]; 00369 u[(*i)->offset+k] = (*i)->upperLimit.v[k]; 00370 } 00371 } 00372 00373 Solver->loadProblem(n, m, Cst, Rnr, Elm, l, u, c, bl, bu); 00374 00375 if (nz>0) { 00376 delete [] Elm; 00377 delete [] Rnr; 00378 } 00379 delete [] Cst; 00380 delete [] Clg; 00381 if (n>0) { 00382 delete [] l; 00383 delete [] u; 00384 delete [] c; 00385 } 00386 if (m>0) { 00387 delete [] bl; 00388 delete [] bu; 00389 } 00390 00391 for (varIt i=Variables.begin(); i!=Variables.end(); i++) { 00392 int begin = (*i)->offset; 00393 int end = (*i)->offset+(*i)->size(); 00394 if ((*i)->type == discrete) { 00395 for (int k=begin; k<end; k++) { 00396 Solver->setInteger(k); 00397 } 00398 } 00399 } 00400 mSolverState = MP_model::ATTACHED; 00401 messenger->generationTime(CoinCpuTime()-time); 00402 } 00403 00404 void MP_model::detach() { 00405 assert(Solver); 00406 mSolverState = MP_model::DETACHED; 00407 delete Solver; 00408 Solver = 0; 00409 } 00410 00411 MP_model::MP_status MP_model::solve(const MP_model::MP_direction &dir) { 00412 assert(Solver); 00413 assert(mSolverState != MP_model::DETACHED && 00414 mSolverState != MP_model::SOLVER_ONLY); 00415 Solver->setObjSense(dir); 00416 bool isMIP = false; 00417 for (varIt i=Variables.begin(); i!=Variables.end(); i++) { 00418 if ((*i)->type == discrete) { 00419 isMIP = true; 00420 break; 00421 } 00422 } 00423 if (isMIP == true) { 00424 try { 00425 Solver->branchAndBound(); 00426 } catch (CoinError e) { 00427 cout<<e.message()<<endl; 00428 cout<<"Solving the LP relaxation instead."<<endl; 00429 try { 00430 Solver->initialSolve(); 00431 } catch (CoinError e) { 00432 cout<<e.message()<<endl; 00433 } 00434 } 00435 } else { 00436 try { 00437 Solver->initialSolve(); 00438 } catch (CoinError e) { 00439 cout<<e.message()<<endl; 00440 } 00441 } 00442 00443 if (Solver->isProvenOptimal() == true) { 00444 cout<<"FlopCpp: Optimal obj. value = "<<Solver->getObjValue()<<endl; 00445 cout<<"FlopCpp: Solver(m, n, nz) = "<<Solver->getNumRows()<<" "<< 00446 Solver->getNumCols()<<" "<<Solver->getNumElements()<<endl; 00447 } else if (Solver->isProvenPrimalInfeasible() == true) { 00448 return mSolverState=MP_model::PRIMAL_INFEASIBLE; 00449 cout<<"FlopCpp: Problem is primal infeasible."<<endl; 00450 } else if (Solver->isProvenDualInfeasible() == true) { 00451 return mSolverState=MP_model::DUAL_INFEASIBLE; 00452 cout<<"FlopCpp: Problem is dual infeasible."<<endl; 00453 } else { 00454 return mSolverState=MP_model::ABANDONED; 00455 cout<<"FlopCpp: Solution process abandoned."<<endl; 00456 } 00457 return mSolverState=MP_model::OPTIMAL; 00458 } 00459 00460 namespace flopc { 00461 std::ostream &operator<<(std::ostream &os, 00462 const MP_model::MP_status &condition) { 00463 switch(condition) { 00464 case MP_model::OPTIMAL: 00465 os<<"OPTIMAL"; 00466 break; 00467 case MP_model::PRIMAL_INFEASIBLE: 00468 os<<"PRIMAL_INFEASIBLE"; 00469 break; 00470 case MP_model::DUAL_INFEASIBLE: 00471 os<<"DUAL_INFEASIBLE"; 00472 break; 00473 case MP_model::ABANDONED: 00474 os<<"ABANDONED"; 00475 break; 00476 default: 00477 os<<"UNKNOWN"; 00478 }; 00479 return os; 00480 } 00481 00482 std::ostream &operator<<(std::ostream &os, 00483 const MP_model::MP_direction &direction) { 00484 switch(direction) { 00485 case MP_model::MINIMIZE: 00486 os<<"MINIMIZE"; 00487 break; 00488 case MP_model::MAXIMIZE: 00489 os<<"MAXIMIZE"; 00490 break; 00491 default: 00492 os<<"UNKNOWN"; 00493 }; 00494 return os; 00495 } 00496 }