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<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
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 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
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 {
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
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
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
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
00337 break;
00338 }
00339 }
00340 }
00341
00342
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
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 }