• Main Page
  • Related Pages
  • Modules
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

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<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     //int row=cfs[j].row;
00115     double elm=cfs[j].val;
00116     //cout<<row<<"   "<<col<<"  "<<elm<<endl;
00117     if (col>=0) {
00118       newRow.insert(col,elm);
00119     } else if (col==-1) {
00120       rhs = elm;
00121     }
00122   }
00123   // NB no "assembly" of added row yet. Will be done....
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         // Nothing to do
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 {  // use pre-attached solver.
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   // Generate coefficient matrix and right hand side
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   // Treat right hand side as n'th column
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   // Row bounds
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             // Nothing to do
00337             break;                
00338       }
00339     }
00340   }
00341 
00342   // Generate objective function coefficients
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   // Column bounds
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 }

Generated on Thu Jul 29 2010 20:00:14 for FLOPC++ by  doxygen 1.7.1