MP_model.cpp

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

Generated on Fri Aug 26 03:02:58 2011 for FLOPC++ by  doxygen 1.4.7