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