00001
00002
00003
00004
00005
00006
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<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<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 solution(0), messenger(m), Objective(0), Solver(s),
00064 m(0), n(0), nz(0), bl(0),
00065 mSolverState(((s==NULL)?(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& c) {
00078 Constraints.insert(&c);
00079 return *this;
00080 }
00081
00082 void MP_model::add(MP_constraint* c) {
00083 c->M = this;
00084 if (c->left.isDefined() && c->right.isDefined()) {
00085 c->offset = m;
00086 m += c->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& c) {
00105 vector<Coef> cfs;
00106 vector<Constant> v;
00107 ObjectiveGenerateFunctor f(cfs);
00108 c.left->generate(MP_domain::getEmpty(),v,f,1.0);
00109 c.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
00115 double elm=cfs[j].val;
00116
00117 if (col>=0) {
00118 newRow.insert(col,elm);
00119 } else if (col==-1) {
00120 rhs = elm;
00121 }
00122 }
00123
00124
00125 double bl = -rhs;
00126 double bu = -rhs;
00127
00128 double inf = Solver->getInfinity();
00129 switch (c.sense) {
00130 case LE:
00131 bl = - inf;
00132 break;
00133 case GE:
00134 bu = inf;
00135 break;
00136 case EQ:
00137
00138 break;
00139 }
00140
00141 Solver->addRow(newRow,bl,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 c(s);
00151 add(c);
00152 c(s) = v() >= obj;
00153 minimize(v());
00154 }
00155
00156
00157 class flopc::CoefLess {
00158 public:
00159 bool operator() (const Coef& a, const Coef& b) const {
00160 if (a.col < b.col) {
00161 return true;
00162 } else if (a.col == b.col && a.row < b.row) {
00163 return true;
00164 } else {
00165 return false;
00166 }
00167 }
00168 };
00169
00170 void MP_model::assemble(vector<Coef>& v, vector<Coef>& av) {
00171 std::sort(v.begin(),v.end(),CoefLess());
00172 int c,r,s;
00173 double val;
00174 std::vector<Coef>::const_iterator i = v.begin();
00175 while (i!=v.end()) {
00176 c = i->col;
00177 r = i->row;
00178 val = i->val;
00179 s = i->stage;
00180 i++;
00181 while (i!=v.end() && c==i->col && r==i->row) {
00182 val += i->val;
00183 if (i->stage>s) {
00184 s = i->stage;
00185 }
00186 i++;
00187 }
00188 av.push_back(Coef(c,r,val,s));
00189 }
00190 }
00191
00192 void MP_model::maximize() {
00193 if (Solver!=0) {
00194 attach(Solver);
00195 solve(MP_model::MAXIMIZE);
00196 } else {
00197 cout<<"no solver specified"<<endl;
00198 }
00199 }
00200
00201 void MP_model::maximize(const MP_expression &obj) {
00202 if (Solver!=0) {
00203 Objective = obj;
00204 attach(Solver);
00205 solve(MP_model::MAXIMIZE);
00206 } else {
00207 cout<<"no solver specified"<<endl;
00208 }
00209 }
00210
00211 void MP_model::minimize() {
00212 if (Solver!=0) {
00213 attach(Solver);
00214 solve(MP_model::MINIMIZE);
00215 } else {
00216 cout<<"no solver specified"<<endl;
00217 }
00218 }
00219
00220 void MP_model::minimize(const MP_expression &obj) {
00221 if (Solver!=0) {
00222 Objective = obj;
00223 attach(Solver);
00224 solve(MP_model::MINIMIZE);
00225 } else {
00226 cout<<"no solver specified"<<endl;
00227 }
00228 }
00229
00230 void MP_model::attach(OsiSolverInterface *_solver) {
00231 if (_solver==NULL) {
00232 if (Solver==NULL) {
00233 mSolverState = MP_model::DETACHED;
00234 return;
00235 }
00236 } else {
00237 if(Solver && Solver!=_solver) {
00238 detach();
00239 }
00240 Solver=_solver;
00241 }
00242 double time = CoinCpuTime();
00243 m=0;
00244 n=0;
00245 vector<Coef> coefs;
00246 vector<Coef> cfs;
00247
00248 typedef std::set<MP_variable* >::iterator varIt;
00249 typedef std::set<MP_constraint* >::iterator conIt;
00250
00251 Objective->insertVariables(Variables);
00252 for (conIt i=Constraints.begin(); i!=Constraints.end(); i++) {
00253 add(*i);
00254 (*i)->insertVariables(Variables);
00255 }
00256 for (varIt j=Variables.begin(); j!=Variables.end(); j++) {
00257 add(*j);
00258 }
00259
00260
00261 bool doAssemble = true;
00262 if (doAssemble == true) {
00263 GenerateFunctor f(cfs);
00264 for (conIt i=Constraints.begin(); i!=Constraints.end(); i++) {
00265 (*i)->coefficients(f);
00266 messenger->constraintDebug((*i)->getName(),cfs);
00267 assemble(cfs,coefs);
00268 cfs.erase(cfs.begin(),cfs.end());
00269 }
00270 } else {
00271 GenerateFunctor f(coefs);
00272 for (conIt i=Constraints.begin(); i!=Constraints.end(); i++) {
00273 (*i)->coefficients(f);
00274 }
00275 }
00276 nz = coefs.size();
00277
00278 messenger->statistics(Constraints.size(),m,Variables.size(),n,nz);
00279
00280 if (nz>0) {
00281 Elm = new double[nz];
00282 Rnr = new int[nz];
00283 }
00284 Cst = new int[n+2];
00285 Clg = new int[n+1];
00286 if (n>0) {
00287 l = new double[n];
00288 u = new double[n];
00289 }
00290 if (m>0) {
00291 bl = new double[m];
00292 bu = new double[m];
00293 }
00294 const double inf = Solver->getInfinity();
00295
00296 for (int j=0; j<n; j++) {
00297 Clg[j] = 0;
00298 }
00299 Clg[n] = 0;
00300
00301
00302 for (int j=0; j<=n; j++) {
00303 Clg[j] = 0;
00304 }
00305 for (int i=0; i<nz; i++) {
00306 int col = coefs[i].col;
00307 if (col == -1) {
00308 col = n;
00309 }
00310 Clg[col]++;
00311 }
00312 Cst[0]=0;
00313 for (int j=0; j<=n; j++) {
00314 Cst[j+1]=Cst[j]+Clg[j];
00315 }
00316 for (int i=0; i<=n; i++) {
00317 Clg[i]=0;
00318 }
00319 for (int i=0; i<nz; i++) {
00320 int col = coefs[i].col;
00321 if (col==-1) {
00322 col = n;
00323 }
00324 int row = coefs[i].row;
00325 double elm = coefs[i].val;
00326 Elm[Cst[col]+Clg[col]] = elm;
00327 Rnr[Cst[col]+Clg[col]] = row;
00328 Clg[col]++;
00329 }
00330
00331
00332 for (int i=0; i<m; i++) {
00333 bl[i] = 0;
00334 bu[i] = 0;
00335 }
00336 for (int j=Cst[n]; j<Cst[n+1]; j++) {
00337 bl[Rnr[j]] = -Elm[j];
00338 bu[Rnr[j]] = -Elm[j];
00339 }
00340
00341 for (conIt i=Constraints.begin(); i!=Constraints.end(); i++) {
00342 if ((*i)->left.isDefined() && (*i)->right.isDefined() ) {
00343 int begin = (*i)->offset;
00344 int end = (*i)->offset+(*i)->size();
00345 switch ((*i)->sense) {
00346 case LE:
00347 for (int k=begin; k<end; k++) {
00348 bl[k] = - inf;
00349 }
00350 break;
00351 case GE:
00352 for (int k=begin; k<end; k++) {
00353 bu[k] = inf;
00354 }
00355 break;
00356 case EQ:
00357
00358 break;
00359 }
00360 }
00361 }
00362
00363
00364 vector<Constant> v;
00365 if (doAssemble == true) {
00366 ObjectiveGenerateFunctor f(cfs);
00367 coefs.erase(coefs.begin(),coefs.end());
00368 Objective->generate(MP_domain::getEmpty(), v, f, 1.0);
00369
00370 messenger->objectiveDebug(cfs);
00371 assemble(cfs,coefs);
00372 } else {
00373 ObjectiveGenerateFunctor f(coefs);
00374 coefs.erase(coefs.begin(),coefs.end());
00375 Objective->generate(MP_domain::getEmpty(), v, f, 1.0);
00376 }
00377
00378 c = new double[n];
00379 for (int j=0; j<n; j++) {
00380 c[j] = 0.0;
00381 }
00382 for (size_t i=0; i<coefs.size(); i++) {
00383 int col = coefs[i].col;
00384 double elm = coefs[i].val;
00385 c[col] = elm;
00386 }
00387
00388
00389 for (int j=0; j<n; j++) {
00390 l[j] = 0.0;
00391 u[j] = inf;
00392 }
00393
00394 for (varIt i=Variables.begin(); i!=Variables.end(); i++) {
00395 for (int k=0; k<(*i)->size(); k++) {
00396 l[(*i)->offset+k] = (*i)->lowerLimit.v[k];
00397 u[(*i)->offset+k] = (*i)->upperLimit.v[k];
00398 }
00399 }
00400
00401 Solver->loadProblem(n, m, Cst, Rnr, Elm, l, u, c, bl, bu);
00402
00403 if (nz>0) {
00404 delete [] Elm;
00405 delete [] Rnr;
00406 }
00407 delete [] Cst;
00408 delete [] Clg;
00409 if (n>0) {
00410 delete [] l;
00411 delete [] u;
00412 }
00413 if (m>0) {
00414 delete [] bl;
00415 delete [] bu;
00416 }
00417
00418 for (varIt i=Variables.begin(); i!=Variables.end(); i++) {
00419 int begin = (*i)->offset;
00420 int end = (*i)->offset+(*i)->size();
00421 if ((*i)->type == discrete) {
00422 for (int k=begin; k<end; k++) {
00423 Solver->setInteger(k);
00424 }
00425 }
00426 }
00427 mSolverState = MP_model::ATTACHED;
00428 messenger->generationTime(time-CoinCpuTime());
00429 }
00430
00431 void MP_model::detach() {
00432 assert(Solver);
00433 mSolverState=MP_model::DETACHED;
00435 delete Solver;
00436 Solver=NULL;
00437 }
00438
00439 MP_model::MP_status MP_model::solve(const MP_model::MP_direction &dir) {
00440 assert(Solver);
00441 assert(mSolverState != MP_model::DETACHED &&
00442 mSolverState != MP_model::SOLVER_ONLY);
00443 Solver->setObjSense(dir);
00444 bool isMIP = false;
00445 for (varIt i=Variables.begin(); i!=Variables.end(); i++) {
00446 if ((*i)->type == discrete) {
00447 isMIP = true;
00448 break;
00449 }
00450 }
00451 if (isMIP == true) {
00452 try {
00453 Solver->branchAndBound();
00454 } catch (CoinError e) {
00455 cout<<e.message()<<endl;
00456 cout<<"Solving the LP relaxation instead."<<endl;
00457 try {
00458 Solver->initialSolve();
00459 } catch (CoinError e) {
00460 cout<<e.message()<<endl;
00461 }
00462 }
00463 } else {
00464 try {
00465 Solver->initialSolve();
00466 } catch (CoinError e) {
00467 cout<<e.message()<<endl;
00468 }
00469 }
00470
00471 if (Solver->isProvenOptimal() == true) {
00472 cout<<"FlopCpp: Optimal obj. value = "<<Solver->getObjValue()<<endl;
00473 cout<<"FlopCpp: Solver(m, n, nz) = "<<Solver->getNumRows()<<" "<<
00474 Solver->getNumCols()<<" "<<Solver->getNumElements()<<endl;
00475 solution = Solver->getColSolution();
00476 reducedCost = Solver->getReducedCost();
00477 rowPrice = Solver->getRowPrice();
00478 rowActivity = Solver->getRowActivity();
00479 } else if (Solver->isProvenPrimalInfeasible() == true) {
00480 return mSolverState=MP_model::PRIMAL_INFEASIBLE;
00481
00482 } else if (Solver->isProvenDualInfeasible() == true) {
00483 return mSolverState=MP_model::DUAL_INFEASIBLE;
00484
00485 } else {
00486 return mSolverState=MP_model::ABANDONED;
00487
00488 }
00489 return mSolverState=MP_model::OPTIMAL;
00490 }
00491
00492 namespace flopc {
00493 std::ostream &operator<<(std::ostream &os,
00494 const MP_model::MP_status &condition) {
00495 switch(condition) {
00496 case MP_model::OPTIMAL:
00497 os<<"OPTIMAL";
00498 break;
00499 case MP_model::PRIMAL_INFEASIBLE:
00500 os<<"PRIMAL_INFEASIBLE";
00501 break;
00502 case MP_model::DUAL_INFEASIBLE:
00503 os<<"DUAL_INFEASIBLE";
00504 break;
00505 case MP_model::ABANDONED:
00506 os<<"ABANDONED";
00507 break;
00508 default:
00509 os<<"UNKNOWN";
00510 };
00511 return os;
00512 }
00513
00514 std::ostream &operator<<(std::ostream &os,
00515 const MP_model::MP_direction &direction) {
00516 switch(direction) {
00517 case MP_model::MINIMIZE:
00518 os<<"MINIMIZE";
00519 break;
00520 case MP_model::MAXIMIZE:
00521 os<<"MAXIMIZE";
00522 break;
00523 default:
00524 os<<"UNKNOWN";
00525 };
00526 return os;
00527 }
00528 }