This feed contains pages with tag "optimization".
In an effort to promote interaction between the Computational Geometry and Optimization research communities, we are editing a special issue of Computational Geometry: Theory and Applications (CGTA) on Geometric Aspects of Optimization
In the last 25 years, computational geometry has become very good at dealing with enormous numbers of input objects, but still faces many challenges in high dimensional settings. Techniques from optimization, although typically slower for many objects in small dimensions, are usually more effective in the face of many variables/dimensions.
On the other hand, the geometric and combinatorial approach that is central to computational geometry has led to important insights in optimization, particularly in algorithms for linear programming.
Original papers with a geometric and algorithmic flavour are invited on topics including, but not limited to:
- Interactions between graph theory, geometry and optimization.
- Geometry and combinatorics of linear and semidefinite programming
- Convex Polyhedra, Convex Bodies, Hyperplane Arrangements
- Geometric Optimization Problems
- Oriented Matroids
- Data depth, and related problems.
All submissions will be refereed following the usual high standards of the journal. Papers should be submitted using the CGTA electronic submission system
Under "Issue", please select "GeomOpt10" and, under "Suggesting Editor", please select "SpecialIssueEditors GeomOpt10"
To ensure consideration for the special issue, please submit by July 1, 2010
Every effort will be made to ensure the timeliness of the special issue; in addition papers in the special issue will appear on-line as soon as accepted.
If you have any questions, feel free to contact one of the guest editors:
- David Bremner bremner@unb.ca
- Antoine Deza deza@mcmaster.ca
- Hiroshi Imai imai@is.s.u-tokyo.ac.jp
- Sonoko Moriyama moriso@is.s.u-tokyo.ac.jp
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);
}
}