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);
}
}