This feed contains pages with tag "integer_program".

Recently I suggested to some students that they could use the Gnu Linear Programming Toolkit from C++. Shortly afterwards I thought I had better verify that I had not just sent people on a hopeless mission. To test things out, I decided to try using GLPK as part of an ongoing project with Lars Schewe

The basic idea of this example is to use glpk to solve an integer program with row generation.

The main hurdle (assuming you want to actually write object oriented c++) is how to make the glpk callback work in an object oriented way. Luckily glpk provides a pointer "info" that can be passed to the solver, and which is passed back to the callback routine. This can be used to keep track of what object is involved.

Here is the class header

#ifndef GLPSOL_HH #define GLPSOL_HH #include "LP.hh" #include "Vektor.hh" #include "glpk.h" #include "combinat.hh" namespace mpc { class GLPSol : public LP { private: glp_iocp parm; static Vektor<double> get_primal_sol(glp_prob *prob); static void callback(glp_tree *tree, void *info); static int output_handler(void *info, const char *s); protected: glp_prob *root; public: GLPSol(int columns); ~GLPSol() {}; virtual void rowgen(const Vektor<double> &candidate) {}; bool solve(); bool add(const LinearConstraint &cnst); }; } #endif

The class `LP`

is just an abstract base class (like an interface for
java-heads) defining the `add`

method. The method `rowgen`

is virtual
because it is intended to be overridden by a subclass if row
generation is actually required. By default it does nothing.

Notice that the callback method here is static; that means it is
essentially a `C`

function with a funny name. This will be the
function that glpk calls when it wants from help.

#include <assert.h> #include "GLPSol.hh" #include "debug.hh" namespace mpc{ GLPSol::GLPSol(int columns) { // redirect logging to my handler glp_term_hook(output_handler,NULL); // make an LP problem root=glp_create_prob(); glp_add_cols(root,columns); // all of my variables are binary, my objective function is always the same // your milage may vary for (int j=1; j<=columns; j++){ glp_set_obj_coef(root,j,1.0); glp_set_col_kind(root,j,GLP_BV); } glp_init_iocp(&parm); // here is the interesting bit; we pass the address of the current object // into glpk along with the callback function parm.cb_func=GLPSol::callback; parm.cb_info=this; } int GLPSol::output_handler(void *info, const char *s){ DEBUG(1) << s; return 1; } Vektor<double> GLPSol::get_primal_sol(glp_prob *prob){ Vektor<double> sol; assert(prob); for (int i=1; i<=glp_get_num_cols(prob); i++){ sol[i]=glp_get_col_prim(prob,i); } return sol; } // the callback function just figures out what object called glpk and forwards // the call. I happen to decode the solution into a more convenient form, but // you can do what you like void GLPSol::callback(glp_tree *tree, void *info){ GLPSol *obj=(GLPSol *)info; assert(obj); switch(glp_ios_reason(tree)){ case GLP_IROWGEN: obj->rowgen(get_primal_sol(glp_ios_get_prob(tree))); break; default: break; } } bool GLPSol::solve(void) { int ret=glp_simplex(root,NULL); if (ret==0) ret=glp_intopt(root,&parm); if (ret==0) return (glp_mip_status(root)==GLP_OPT); else return false; } bool GLPSol::add(const LinearConstraint&cnst){ int next_row=glp_add_rows(root,1); // for mysterious reasons, glpk wants to index from 1 int indices[cnst.size()+1]; double coeff[cnst.size()+1]; DEBUG(3) << "adding " << cnst << std::endl; int j=1; for (LinearConstraint::const_iterator p=cnst.begin(); p!=cnst.end(); p++){ indices[j]=p->first; coeff[j]=(double)p->second; j++; } int gtype=0; switch(cnst.type()){ case LIN_LEQ: gtype=GLP_UP; break; case LIN_GEQ: gtype=GLP_LO; break; default: gtype=GLP_FX; } glp_set_row_bnds(root,next_row,gtype, (double)cnst.rhs(),(double)cnst.rhs()); glp_set_mat_row(root, next_row, cnst.size(), indices, coeff); return true; } }

All this is a big waste of effort unless we actually do some row generation. I'm not especially proud of the crude rounding I do here, but it shows how to do it, and it does, eventually solve problems.

#include "OMGLPSol.hh" #include "DualGraph.hh" #include "CutIterator.hh" #include "IntSet.hh" namespace mpc{ void OMGLPSol::rowgen(const Vektor<double>&candidate){ if (diameter<=0){ DEBUG(1) << "no path constraints to generate" << std::endl; return; } DEBUG(3) << "Generating paths for " << candidate << std::endl; // this looks like a crude hack, which it is, but motivated by the // following: the boundary complex is determined only by the signs // of the bases, which we here represent as 0 for - and 1 for + Chirotope chi(*this); for (Vektor<double>::const_iterator p=candidate.begin(); p!=candidate.end(); p++){ if (p->second > 0.5) { chi[p->first]=SIGN_POS; } else { chi[p->first]=SIGN_NEG; } } BoundaryComplex bc(chi); DEBUG(3) << chi; DualGraph dg(bc); CutIterator pathins(*this,candidate); int paths_found= dg.all_paths(pathins, IntSet::lex_set(elements(),rank()-1,source_facet), IntSet::lex_set(elements(),rank()-1,sink_facet), diameter-1); DEBUG(1) << "row generation found " << paths_found << " realized paths\n"; DEBUG(1) << "effective cuts: " << pathins.effective() << std::endl; } void OMGLPSol::get_solution(Chirotope &chi) { int nv=glp_get_num_cols(root); for(int i=1;i<=nv;++i) { int val=glp_mip_col_val(root,i); chi[i]=(val==0 ? SIGN_NEG : SIGN_POS); } } }

So ignore the problem specific way I generate constraints, the key
remaining piece of code is `CutIterator`

which filters the generated
constraints to make sure they actually cut off the candidate
solution. This is crucial, because row generation must not add
constraints in the case that it cannot improve the solution, because
glpk assumes that if the user is generating cuts, the solver doesn't
have to.

#ifndef PATH_CONSTRAINT_ITERATOR_HH #define PATH_CONSTRAINT_ITERATOR_HH #include "PathConstraint.hh" #include "CNF.hh" namespace mpc { class CutIterator : public std::iterator<std::output_iterator_tag, void, void, void, void>{ private: LP& _list; Vektor<double> _sol; std::size_t _pcount; std::size_t _ccount; public: CutIterator (LP& list, const Vektor<double>& sol) : _list(list),_sol(sol), _pcount(0), _ccount(0) {} CutIterator& operator=(const Path& p) { PathConstraint pc(p); _ccount+=pc.appendTo(_list,&_sol); _pcount++; if (_pcount %10000==0) { DEBUG(1) << _pcount << " paths generated" << std::endl; } return *this; } CutIterator& operator*() {return *this;} CutIterator& operator++() {return *this;} CutIterator& operator++(int) {return *this;} int effective() { return _ccount; }; }; } #endif

Oh heck, another level of detail; the actual filtering actually
happens in the `appendTo`

method the PathConstraint class. This is
just computing the dot product of two vectors. I would leave it as an
exercise to the readier, but remember some fuzz is neccesary to to
these kinds of comparisons with floating point numbers. Eventually,
the decision is made by the following `feasible`

method of the
`LinearConstraint`

class.

bool feasible(const Vektor<double> & x){ double sum=0; for (const_iterator p=begin();p!=end(); p++){ sum+= p->second*x.at(p->first); } switch (type()){ case LIN_LEQ: return (sum <= _rhs+epsilon); case LIN_GEQ: return (sum >= _rhs-epsilon); default: return (sum <= _rhs+epsilon) && (sum >= _rhs-epsilon); } }

In CS3997 this year, I have decided to have "debates" instead of
presentations. This means that I need to make sure that each student
has has exactly one topic, **and** every topic is assigned to an even
number of students.

Another constraint is that I asked the students to list in order their top 3 preferences. I wanted to maximize (within reason) "student happiness", so I give weight 4 for their first choice, 2 for the second, and 1 for their third.

Finally, the students are numbered `1...26`

in the order they sent me
their preferences. I decided to enforce "first-come first-serve" in
the objective function, so the happiness of student 1 has more weight
than student 26. How much more is a bit of a subjective choice.

If you don't want to look at the solution yet, the students preferences are available separately. 'happy[i,j]' measures how happy student i is being assigned topic j

Almost the real solution is available. In actuallity, I first solved the problem for the first 18 students (so they didn't have to wait), and use the following

printf { i in students, j in topics : x[i,j]=1 } "s.t. fix_%d_%d: x[%d,%d]=1;\n",i,j,i,j;

to print out some constraints, which I then cut and past into the model, and resolved a week or so later when I had all of the data.