FLOPC++
MP_model.cpp
Go to the documentation of this file.
1 // ******************** FlopCpp **********************************************
2 // File: MP_model.cpp
3 // $Id$
4 // Author: Tim Helge Hultberg (thh@mat.ua.pt)
5 // Copyright (C) 2003 Tim Helge Hultberg
6 // All Rights Reserved.
7 //****************************************************************************
8 
9 #include <iostream>
10 #include <sstream>
11 using std::cout;
12 using std::endl;
13 #include <algorithm>
14 
15 #include <CoinPackedMatrix.hpp>
16 #include <OsiSolverInterface.hpp>
17 #include "MP_model.hpp"
18 #include "MP_variable.hpp"
19 #include "MP_constraint.hpp"
20 #include <CoinTime.hpp>
21 
22 using namespace flopc;
23 
28 
29 void NormalMessenger::statistics(int bm, int m, int bn, int n, int nz) {
30  cout<<"FlopCpp: Number of constraint blocks: " <<bm<<endl;
31  cout<<"FlopCpp: Number of individual constraints: " <<m<<endl;
32  cout<<"FlopCpp: Number of variable blocks: " <<bn<<endl;
33  cout<<"FlopCpp: Number of individual variables: " <<n<<endl;
34  cout<<"FlopCpp: Number of non-zeroes (including rhs): " <<nz<<endl;
35 }
36 
38  cout<<"FlopCpp: Generation time: "<<t<<endl;
39 }
40 
41 void VerboseMessenger::constraintDebug(string name, const vector<MP::Coef>& cfs) {
42  cout<<"FlopCpp: Constraint "<<name<<endl;
43  for (unsigned int j=0; j<cfs.size(); j++) {
44  int col=cfs[j].col;
45  int row=cfs[j].row;
46  double elm=cfs[j].val;
47  int stage=cfs[j].stage;
48  cout<<row<<" "<<col<<" "<<elm<<" "<<stage<<endl;
49  }
50 }
51 
52 void VerboseMessenger::objectiveDebug(const vector<MP::Coef>& cfs) {
53  cout<<"Objective "<<endl;
54  for (unsigned int j=0; j<cfs.size(); j++) {
55  int col=cfs[j].col;
56  int row=cfs[j].row;
57  double elm=cfs[j].val;
58  cout<<row<<" "<<col<<" "<<elm<<endl;
59  }
60 }
61 
62 MP_model::MP_model(OsiSolverInterface* s, Messenger* m) :
63  messenger(m), Objective(0), Solver(s),
64  m(0), n(0), nz(0), bl(0),
66  MP_model::current_model = this;
67 }
68 
70  delete messenger;
71  if (Solver!=0) {
72  delete Solver;
73  }
74 }
75 
76 
78  Constraints.insert(&constraint);
79  return *this;
80 }
81 
82 void MP_model::add(MP_constraint* constraint) {
83  constraint->M = this;
84  if (constraint->left.isDefined() && constraint->right.isDefined()) {
85  constraint->offset = m;
86  m += constraint->size();
87  }
88 }
89 
90 double MP_model::getInfinity() const {
91  if (Solver==0) {
92  return 9.9e+32;
93  } else {
94  return Solver->getInfinity();
95  }
96 }
97 
99  v->M = this;
100  v->offset = n;
101  n += v->size();
102 }
103 
104 void MP_model::addRow(const Constraint& constraint) {
105  vector<MP::Coef> cfs;
106  vector<Constant> v;
107  MP::GenerateFunctor f(0,cfs);
108  constraint->left->generate(MP_domain::getEmpty(),v,f,1.0);
109  constraint->right->generate(MP_domain::getEmpty(),v,f,-1.0);
110  CoinPackedVector newRow;
111  double rhs = 0.0;
112  for (unsigned int j=0; j<cfs.size(); j++) {
113  int col=cfs[j].col;
114  //int row=cfs[j].row;
115  double elm=cfs[j].val;
116  //cout<<row<<" "<<col<<" "<<elm<<endl;
117  if (col>=0) {
118  newRow.insert(col,elm);
119  } else if (col==-1) {
120  rhs = elm;
121  }
122  }
123  // NB no "assembly" of added row yet. Will be done....
124 
125  double local_bl = -rhs;
126  double local_bu = -rhs;
127 
128  double inf = Solver->getInfinity();
129  switch (constraint->sense) {
130  case LE:
131  local_bl = - inf;
132  break;
133  case GE:
134  local_bu = inf;
135  break;
136  case EQ:
137  // Nothing to do
138  break;
139  }
140 
141  Solver->addRow(newRow,local_bl,local_bu);
142 }
143 
145  Objective = o;
146 }
147 
149  MP_variable v;
150  MP_constraint constraint(s);
151  add(constraint);
152  constraint(s) = v() >= obj;
153  minimize(v());
154 }
155 
156 
157 void MP_model::assemble(vector<MP::Coef>& v, vector<MP::Coef>& av) {
158  std::sort(v.begin(),v.end(),MP::CoefLess());
159  int c,r,s;
160  double val;
161  std::vector<MP::Coef>::const_iterator i = v.begin();
162  while (i!=v.end()) {
163  c = i->col;
164  r = i->row;
165  val = i->val;
166  s = i->stage;
167  i++;
168  while (i!=v.end() && c==i->col && r==i->row) {
169  val += i->val;
170  if (i->stage>s) {
171  s = i->stage;
172  }
173  i++;
174  }
175  av.push_back(MP::Coef(c,r,val,s));
176  }
177 }
178 
180  if (Solver!=0) {
181  attach(Solver);
183  } else {
184  cout<<"no solver specified"<<endl;
185  }
186 }
187 
189  if (Solver!=0) {
190  Objective = obj;
191  attach(Solver);
193  } else {
194  cout<<"no solver specified"<<endl;
195  }
196 }
197 
199  if (Solver!=0) {
200  attach(Solver);
202  } else {
203  cout<<"no solver specified"<<endl;
204  }
205 }
206 
208  if (Solver!=0) {
209  Objective = obj;
210  attach(Solver);
212  } else {
213  cout<<"no solver specified"<<endl;
214  }
215 }
216 
217 void MP_model::attach(OsiSolverInterface *_solver) {
218  if (_solver == 0) {
219  if (Solver == 0) {
221  return;
222  }
223  } else { // use pre-attached solver.
224  if(Solver && Solver!=_solver) {
225  detach();
226  }
227  Solver=_solver;
228  }
229  double time = CoinCpuTime();
230  m=0;
231  n=0;
232  vector<MP::Coef> coefs;
233  vector<MP::Coef> cfs;
234 
235  typedef std::set<MP_variable* >::iterator varIt;
236  typedef std::set<MP_constraint* >::iterator conIt;
237 
238  Objective->insertVariables(Variables);
239  for (conIt i=Constraints.begin(); i!=Constraints.end(); i++) {
240  add(*i);
241  (*i)->insertVariables(Variables);
242  }
243  for (varIt j=Variables.begin(); j!=Variables.end(); j++) {
244  add(*j);
245  }
246 
247  // Generate coefficient matrix and right hand side
248  for (conIt i=Constraints.begin(); i!=Constraints.end(); i++) {
249  (*i)->coefficients(cfs);
250  messenger->constraintDebug((*i)->getName(),cfs);
251  assemble(cfs,coefs);
252  cfs.erase(cfs.begin(),cfs.end());
253  }
254  nz = coefs.size();
255 
256  messenger->statistics(Constraints.size(),m,Variables.size(),n,nz);
257 
258  if (nz>0) {
259  Elm = new double[nz];
260  Rnr = new int[nz];
261  }
262  Cst = new int[n+2];
263  Clg = new int[n+1];
264  if (n>0) {
265  l = new double[n];
266  u = new double[n];
267  c = new double[n];
268  }
269  if (m>0) {
270  bl = new double[m];
271  bu = new double[m];
272  }
273  const double inf = Solver->getInfinity();
274 
275  for (int j=0; j<n; j++) {
276  Clg[j] = 0;
277  }
278  Clg[n] = 0;
279 
280  // Treat right hand side as n'th column
281  for (int j=0; j<=n; j++) {
282  Clg[j] = 0;
283  }
284  for (int i=0; i<nz; i++) {
285  int col = coefs[i].col;
286  if (col == -1) {
287  col = n;
288  }
289  Clg[col]++;
290  }
291  Cst[0]=0;
292  for (int j=0; j<=n; j++) {
293  Cst[j+1]=Cst[j]+Clg[j];
294  }
295  for (int i=0; i<=n; i++) {
296  Clg[i]=0;
297  }
298  for (int i=0; i<nz; i++) {
299  int col = coefs[i].col;
300  if (col==-1) {
301  col = n;
302  }
303  int row = coefs[i].row;
304  double elm = coefs[i].val;
305  Elm[Cst[col]+Clg[col]] = elm;
306  Rnr[Cst[col]+Clg[col]] = row;
307  Clg[col]++;
308  }
309 
310  // Row bounds
311  for (int i=0; i<m; i++) {
312  bl[i] = 0;
313  bu[i] = 0;
314  }
315  for (int j=Cst[n]; j<Cst[n+1]; j++) {
316  bl[Rnr[j]] = -Elm[j];
317  bu[Rnr[j]] = -Elm[j];
318  }
319 
320  for (conIt i=Constraints.begin(); i!=Constraints.end(); i++) {
321  if ((*i)->left.isDefined() && (*i)->right.isDefined() ) {
322  int begin = (*i)->offset;
323  int end = (*i)->offset+(*i)->size();
324  switch ((*i)->sense) {
325  case LE:
326  for (int k=begin; k<end; k++) {
327  bl[k] = - inf;
328  }
329  break;
330  case GE:
331  for (int k=begin; k<end; k++) {
332  bu[k] = inf;
333  }
334  break;
335  case EQ:
336  // Nothing to do
337  break;
338  }
339  }
340  }
341 
342  // Generate objective function coefficients
343  vector<Constant> v;
344  MP::GenerateFunctor f(0,cfs);
345  coefs.erase(coefs.begin(),coefs.end());
346  Objective->generate(MP_domain::getEmpty(), v, f, 1.0);
347 
349  assemble(cfs,coefs);
350 
351  for (int j=0; j<n; j++) {
352  c[j] = 0.0;
353  }
354  for (size_t i=0; i<coefs.size(); i++) {
355  int col = coefs[i].col;
356  double elm = coefs[i].val;
357  c[col] = elm;
358  }
359 
360  // Column bounds
361  for (int j=0; j<n; j++) {
362  l[j] = 0.0;
363  u[j] = inf;
364  }
365 
366  for (varIt i=Variables.begin(); i!=Variables.end(); i++) {
367  for (int k=0; k<(*i)->size(); k++) {
368  l[(*i)->offset+k] = (*i)->lowerLimit.v[k];
369  u[(*i)->offset+k] = (*i)->upperLimit.v[k];
370  }
371  }
372 
373  Solver->loadProblem(n, m, Cst, Rnr, Elm, l, u, c, bl, bu);
374 
375  if (nz>0) {
376  delete [] Elm;
377  delete [] Rnr;
378  }
379  delete [] Cst;
380  delete [] Clg;
381  if (n>0) {
382  delete [] l;
383  delete [] u;
384  delete [] c;
385  }
386  if (m>0) {
387  delete [] bl;
388  delete [] bu;
389  }
390 
391  for (varIt i=Variables.begin(); i!=Variables.end(); i++) {
392  int begin = (*i)->offset;
393  int end = (*i)->offset+(*i)->size();
394  if ((*i)->type == discrete) {
395  for (int k=begin; k<end; k++) {
396  Solver->setInteger(k);
397  }
398  }
399  }
401  messenger->generationTime(CoinCpuTime()-time);
402 }
403 
405  assert(Solver);
407  delete Solver;
408  Solver = 0;
409 }
410 
412  assert(Solver);
413  assert(mSolverState != MP_model::DETACHED &&
415  Solver->setObjSense(dir);
416  bool isMIP = false;
417  for (varIt i=Variables.begin(); i!=Variables.end(); i++) {
418  if ((*i)->type == discrete) {
419  isMIP = true;
420  break;
421  }
422  }
423  if (isMIP == true) {
424  try {
425  Solver->branchAndBound();
426  } catch (CoinError e) {
427  cout<<e.message()<<endl;
428  cout<<"Solving the LP relaxation instead."<<endl;
429  try {
430  Solver->initialSolve();
431  } catch (CoinError e) {
432  cout<<e.message()<<endl;
433  }
434  }
435  } else {
436  try {
437  Solver->initialSolve();
438  } catch (CoinError e) {
439  cout<<e.message()<<endl;
440  }
441  }
442 
443  if (Solver->isProvenOptimal() == true) {
444  cout<<"FlopCpp: Optimal obj. value = "<<Solver->getObjValue()<<endl;
445  cout<<"FlopCpp: Solver(m, n, nz) = "<<Solver->getNumRows()<<" "<<
446  Solver->getNumCols()<<" "<<Solver->getNumElements()<<endl;
447  } else if (Solver->isProvenPrimalInfeasible() == true) {
449  cout<<"FlopCpp: Problem is primal infeasible."<<endl;
450  } else if (Solver->isProvenDualInfeasible() == true) {
452  cout<<"FlopCpp: Problem is dual infeasible."<<endl;
453  } else {
455  cout<<"FlopCpp: Solution process abandoned."<<endl;
456  }
458 }
459 
460 namespace flopc {
461  std::ostream &operator<<(std::ostream &os,
462  const MP_model::MP_status &condition) {
463  switch(condition) {
464  case MP_model::OPTIMAL:
465  os<<"OPTIMAL";
466  break;
468  os<<"PRIMAL_INFEASIBLE";
469  break;
471  os<<"DUAL_INFEASIBLE";
472  break;
473  case MP_model::ABANDONED:
474  os<<"ABANDONED";
475  break;
476  default:
477  os<<"UNKNOWN";
478  };
479  return os;
480  }
481 
482  std::ostream &operator<<(std::ostream &os,
483  const MP_model::MP_direction &direction) {
484  switch(direction) {
485  case MP_model::MINIMIZE:
486  os<<"MINIMIZE";
487  break;
488  case MP_model::MAXIMIZE:
489  os<<"MAXIMIZE";
490  break;
491  default:
492  os<<"UNKNOWN";
493  };
494  return os;
495  }
496 }
friend class CoefLess
void attach(OsiSolverInterface *solver=0)
attaches the symantic representation of a model and data to a particular OsiSolverInterface ...
Definition: MP_model.cpp:217
static void assemble(vector< MP::Coef > &v, vector< MP::Coef > &av)
Definition: MP_model.cpp:157
MP_model & add(MP_constraint &constraint)
Adds a constrataint block to the model.
Definition: MP_model.cpp:77
virtual void constraintDebug(string name, const vector< MP::Coef > &cfs)
Definition: MP_model.hpp:41
virtual void constraintDebug(string name, const vector< MP::Coef > &cfs)
Definition: MP_model.cpp:41
Symbolic representation of a linear expression.This is one of the main public interface classes...
Messenger * messenger
Definition: MP_model.hpp:241
MP_model(OsiSolverInterface *s, Messenger *m=new NormalMessenger)
Constructs an MP_model from an OsiSolverInterface *.
Definition: MP_model.cpp:62
std::ostream & operator<<(std::ostream &os, const MP_model::MP_status &condition)
allows print of result from call to solve();
Definition: MP_model.cpp:461
double getInfinity() const
Definition: MP_model.cpp:90
static MP_model * getCurrentModel()
Definition: MP_model.cpp:27
std::set< MP_constraint *>::iterator conIt
Definition: MP_model.hpp:235
std::set< MP_variable *>::iterator varIt
Definition: MP_model.hpp:234
Inteface for hooking up to internal flopc++ message handling.In more advanced use of FlopC++...
Definition: MP_model.hpp:36
virtual void statistics(int bm, int m, int bn, int n, int nz)
Definition: MP_model.cpp:29
virtual void generationTime(double t)
Definition: MP_model.hpp:44
void setObjective(const MP_expression &o)
sets the "current objective" to the parameter o
Definition: MP_model.cpp:144
MP_expression Objective
Definition: MP_model.hpp:246
static MP_model * current_model
Definition: MP_model.hpp:237
double * Elm
Definition: MP_model.hpp:259
void addRow(const Constraint &constraint)
Adds a constraint to the MP_model.
Definition: MP_model.cpp:104
MP_status mSolverState
Definition: MP_model.hpp:265
A solver is attached, but not yet solved.
Definition: MP_model.hpp:112
MP_status
Reflects the state of the solution from solve()
Definition: MP_model.hpp:98
This is the anchor point for all constructs in a FlopC++ model.The constructors take an OsiSolverInte...
Definition: MP_model.hpp:90
virtual void statistics(int bm, int m, int bn, int n, int nz)
Definition: MP_model.hpp:43
All flopc++ code is contained within the flopc namespace.
Definition: flopc.cpp:11
OsiSolverInterface * Solver
Definition: MP_model.hpp:251
if solve is called and solver finds model primal infeasible.
Definition: MP_model.hpp:102
MP_model::MP_status solve(const MP_model::MP_direction &dir)
Definition: MP_model.cpp:411
set< MP_variable * > Variables
Definition: MP_model.hpp:248
set< MP_constraint * > Constraints
Definition: MP_model.hpp:247
static const MP_domain & getEmpty()
returns a reference to the "empty" set.
Definition: MP_domain.cpp:48
void minimize_max(MP_set &d, const MP_expression &obj)
Definition: MP_model.cpp:148
Symantic representation of a variable.This is one of the main public interface classes. It should be directly declared by clients of the FlopC++. The parametersof construction are MP_set s which specify the indexes over which the variable is defined.
Definition: MP_variable.hpp:35
if solve is called and solver finds the model dual infeasible.
Definition: MP_model.hpp:104
Representation of a set for indexing into some other construct.This is one of the main public interfa...
Definition: MP_set.hpp:79
virtual void objectiveDebug(const vector< MP::Coef > &cfs)
Definition: MP_model.hpp:42
static MP_model & default_model
Definition: MP_model.hpp:236
static MP_model & getDefaultModel()
Definition: MP_model.cpp:26
No solver is attached.
Definition: MP_model.hpp:114
int size() const
if the solve method is called and the optimal solution found.
Definition: MP_model.hpp:100
virtual void objectiveDebug(const vector< MP::Coef > &cfs)
Definition: MP_model.cpp:52
MP_direction
used when calling the solve() method.
Definition: MP_model.hpp:94
void detach()
detaches an OsiSolverInterface object from the model. In essence, this will clean up any intermediate...
Definition: MP_model.cpp:404
virtual void generationTime(double t)
Definition: MP_model.cpp:37
Semantic representation of a linear constraint.This is one of the main public interface classes...