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 #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         //int row=cfs[j].row;
00104         double elm=cfs[j].val;
00105         //cout<<row<<"   "<<col<<"  "<<elm<<endl;
00106         if (col>=0) {
00107             newRow.insert(col,elm);
00108         } else if (col==-1) {
00109             rhs = elm;
00110         }
00111     }
00112     // NB no "assembly" of added row yet. Will be done....
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             // Nothing to do
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 {  // use pre-attached solver.
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     // Generate coefficient matrix and right hand side
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     // Treat right hand side as n'th column
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     // Row bounds
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                 // Nothing to do
00341                 break;          
00342         }
00343     }
00344 
00345     // Generate objective function coefficients
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     // Column bounds
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     // Instead of the 2 lines above we should be able to use
00387     // the line below, but due to a bug in OsiGlpk it does not work
00388     // Solver->loadProblem(n, m, Cst, Rnr, Elm, l, u, c, bl, bu);
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         //cout<<"FlopCpp: Problem is primal infeasible."<<endl;
00463     } else if (Solver->isProvenDualInfeasible() == true) {
00464         return mSolverState=MP_model::DUAL_INFEASIBLE;
00465         //cout<<"FlopCpp: Problem is dual infeasible."<<endl;
00466     } else {
00467         return mSolverState=MP_model::ABANDONED;
00468         //cout<<"FlopCpp: Solution process abandoned."<<endl;
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 }

Generated on Sun Nov 6 03:14:50 2011 for FLOPC++ by  doxygen 1.4.7