Solvers and API for specialized problems, such as solving special linear programming problems with min-cost flow algorithms. It also wraps solvers like semidefinite programming solver Csdp and convex optimization solver Gurobi and lpsolve.
The analysis and background of using primal min-cost flow to solve linear programming problem can be found in the detailed description of class limbo::solvers::MinCostFlow.
Be aware that only Capacity Scaling algorithm supports real value costs, while all other algorithms only support integer costs. All capacity and flow need to be integers.
#ifndef LIMBO_ALGORITHMS_COLORING_SDPCOLORINGCSDP
#define LIMBO_ALGORITHMS_COLORING_SDPCOLORINGCSDP
{
namespace algorithms
{
namespace coloring
{
template <typename GraphType>
class SDPColoringCsdp : public Coloring<GraphType>
{
public:
typedef Coloring<GraphType> base_type;
using typename base_type::graph_type;
using typename base_type::graph_vertex_type;
using typename base_type::graph_edge_type;
using typename base_type::vertex_iterator_type;
using typename base_type::edge_iterator_type;
using typename base_type::edge_weight_type;
typedef typename base_type::EdgeHashType edge_hash_type;
struct FMGainCalcType
{
edge_weight_type
operator()(int32_t v, int8_t origp, int8_t newp, std::vector<int8_t>
const& vPartition)
const
{
typedef typename boost::graph_traits<graph_type>::out_edge_iterator out_edge_iterator_type;
out_edge_iterator_type ei, eie;
edge_weight_type gain = (origp >= 0)? 0 : boost::num_edges(graph)*boost::num_vertices(graph);
for (boost::tie(ei, eie) = boost::out_edges(v, graph); ei != eie; ++ei)
{
graph_vertex_type t = boost::target(*ei, graph);
int8_t pt = vPartition[t];
#ifdef DEBUG_SDPCOLORING
assert((int32_t)t != v);
#endif
if (pt < 0) continue;
edge_weight_type w = boost::get(boost::edge_weight, graph, *ei);
gain += (pt == newp)? -w : (pt == origp)? w : 0;
}
return gain;
}
};
void write_sdp_sol(std::string
const& filename,
struct blockmatrix
const& X)
const;
protected:
void set_sparseblock_entry(
struct sparseblock& block, int32_t entryid, int32_t i, int32_t j,
double value)
const;
void coloring_algos(graph_type
const& g, std::vector<int8_t>& vColor)
const;
virtual void coloring_by_FM(graph_type
const& mg, std::vector<int8_t>& vColor)
const;
};
template <typename GraphType>
: base_type(g)
{
m_rounding_lb = -0.4;
m_rounding_ub = 0.9;
}
template <typename GraphType>
double SDPColoringCsdp<GraphType>::coloring()
{
assert_msg(!this->has_precolored(),
"SDP coloring does not support precolored layout yet");
struct blockmatrix C;
double *b;
struct constraintmatrix *constraints;
struct blockmatrix X,Z;
double *y;
double pobj,dobj;
vertex_iterator_type vi, vie;
edge_iterator_type ei, eie;
uint32_t num_vertices = boost::num_vertices(this->m_graph);
uint32_t num_edges = boost::num_edges(this->m_graph);
uint32_t num_conflict_edges = 0;
uint32_t num_stitch_edges = 0;
for (boost::tie(ei, eie) = boost::edges(this->m_graph); ei != eie; ++ei)
{
if (this->edge_weight(*ei) >= 0)
num_conflict_edges += 1;
else
num_stitch_edges += 1;
}
assert_msg(num_edges > 0 && num_conflict_edges > 0,
"no edges or conflict edges found, no need to solve SDP");
uint32_t num_variables = num_vertices+num_conflict_edges;
uint32_t num_constraints = num_conflict_edges+num_vertices;
C.nblocks = 2;
C.blocks = (struct blockrec *)malloc((C.nblocks+1)*sizeof(struct blockrec));
assert_msg(C.blocks,
"Couldn't allocate storage for C");
construct_objectve_blockrec(C, 1, num_vertices, MATRIX);
for (boost::tie(ei, eie) = boost::edges(this->m_graph); ei != eie; ++ei)
{
graph_vertex_type s = boost::source(*ei, this->m_graph);
graph_vertex_type t = boost::target(*ei, this->m_graph);
edge_weight_type alpha = (this->edge_weight(*ei) >= 0)? -1 : this->stitch_weight();
s += 1; t += 1;
int32_t idx1 = ijtok(s,t,C.blocks[1].blocksize);
int32_t idx2 = ijtok(t,s,C.blocks[1].blocksize);
C.blocks[1].data.mat[idx1] = alpha;
C.blocks[1].data.mat[idx2] = alpha;
}
construct_objectve_blockrec(C, 2, num_conflict_edges, DIAG);
#ifdef DEBUG_SDPCOLORING
print_blockrec("C.blocks[1].data.mat", C.blocks[1]);
print_blockrec("C.blocks[2].data.vec", C.blocks[2]);
#endif
b = (double *)malloc((num_constraints+1)*sizeof(double));
assert_msg(b,
"Failed to allocate storage for right hand side of constraints b");
double beta = -2.0/(this->color_num()-1.0);
for (uint32_t i = 1, ie = num_constraints+1; i != ie; ++i)
{
if (i <= num_conflict_edges)
b[i] = beta;
else
b[i] = 1;
}
constraints=(struct constraintmatrix *)malloc((num_constraints+1)*sizeof(struct constraintmatrix));
assert_msg(constraints,
"Failed to allocate storage for constraints");
uint32_t cnt = 1;
for (boost::tie(ei, eie) = boost::edges(this->m_graph); ei != eie; ++ei)
{
if (this->edge_weight(*ei) >= 0)
{
graph_vertex_type s = boost::source(*ei, this->m_graph);
graph_vertex_type t = boost::target(*ei, this->m_graph);
if (s > t)
std::swap(s, t);
s += 1; t += 1;
struct constraintmatrix& constr = constraints[cnt];
constr.blocks = NULL;
for (uint32_t i = 2; i != 0; --i)
{
struct sparseblock* blockptr;
if (i == 1)
{
blockptr = construct_constraint_sparseblock(i, num_vertices, cnt, 1);
set_sparseblock_entry(*blockptr, 1, s, t, 1);
}
else
{
blockptr = construct_constraint_sparseblock(i, num_conflict_edges, cnt, 1);
set_sparseblock_entry(*blockptr, 1, cnt, cnt, -1);
}
blockptr->next = constr.blocks;
constr.blocks = blockptr;
}
++cnt;
}
}
for (boost::tie(vi, vie) = boost::vertices(this->m_graph); vi != vie; ++vi)
{
graph_vertex_type v = *vi;
v += 1;
struct constraintmatrix& constr = constraints[cnt];
constr.blocks = NULL;
struct sparseblock* blockptr = construct_constraint_sparseblock(1, num_vertices, cnt, 1);
set_sparseblock_entry(*blockptr, 1, v, v, 1);
blockptr->next = constr.blocks;
constr.blocks = blockptr;
++cnt;
}
#ifdef DEBUG_SDPCOLORING
write_prob((char*)"problem.sdpa", num_variables, num_constraints, C, b, constraints);
#endif
initsoln(num_variables, num_constraints, C, b, constraints, &X, &y, &Z);
struct paramstruc params;
int printlevel;
initparams(¶ms, &printlevel);
printlevel = 0;
int ret = limbo::solvers::easy_sdp_ext<int>(num_variables, num_constraints, C, b, constraints, 0.0, &X, &y, &Z, &pobj, &dobj, params, printlevel);
round_sol(X);
free_prob(num_variables, num_constraints, C, b, constraints, X, y, Z);
#ifdef DEBUG_LPCOLORING
#endif
return this->calc_cost(this->m_vColor);
}
template <typename GraphType>
void SDPColoringCsdp<GraphType>::construct_objectve_blockrec(blockmatrix& C, int32_t blocknum, int32_t blocksize, blockcat blockcategory) const
{
struct blockrec& cblock = C.blocks[blocknum];
cblock.blocksize = blocksize;
cblock.blockcategory = blockcategory;
if (blockcategory == MATRIX)
{
cblock.data.mat = (double *)malloc(blocksize*blocksize*sizeof(double));
assert_msg(cblock.data.mat,
"Couldn't allocate storage for cblock.data.mat");
std::fill(cblock.data.mat, cblock.data.mat+blocksize*blocksize, 0);
}
else if (blockcategory == DIAG)
{
cblock.data.vec = (double *)malloc((blocksize+1)*sizeof(double));
assert_msg(cblock.data.vec,
"Couldn't allocate storage for cblock.data.vec");
std::fill(cblock.data.vec, cblock.data.vec+blocksize+1, 0);
}
}
template <typename GraphType>
struct sparseblock* SDPColoringCsdp<GraphType>::construct_constraint_sparseblock(int32_t blocknum, int32_t blocksize, int32_t constraintnum, int32_t entrynum) const
{
struct sparseblock* blockptr = (struct sparseblock *)malloc(sizeof(struct sparseblock));
assert_msg(blockptr,
"Allocation of constraint block failed for blockptr");
blockptr->blocknum = blocknum;
blockptr->blocksize = blocksize;
blockptr->constraintnum = constraintnum;
blockptr->next = NULL;
blockptr->nextbyblock = NULL;
blockptr->entries = (double *) malloc((entrynum+1)*sizeof(double));
assert_msg(blockptr->entries,
"Allocation of constraint block failed for blockptr->entries");
blockptr->iindices = (int *) malloc((entrynum+1)*sizeof(int));
assert_msg(blockptr->iindices,
"Allocation of constraint block failed for blockptr->iindices");
blockptr->jindices = (int *) malloc((entrynum+1)*sizeof(int));
assert_msg(blockptr->jindices,
"Allocation of constraint block failed for blockptr->jindices");
blockptr->numentries = entrynum;
return blockptr;
}
template <typename GraphType>
void SDPColoringCsdp<GraphType>::set_sparseblock_entry(struct sparseblock& block, int32_t entryid, int32_t i, int32_t j, double value) const
{
block.iindices[entryid] = i;
block.jindices[entryid] = j;
block.entries[entryid] = value;
}
template <typename GraphType>
void SDPColoringCsdp<GraphType>::round_sol(struct blockmatrix const& X)
{
#ifdef DEBUG_SDPCOLORING
write_sdp_sol("problem.sol", X);
#endif
std::vector<graph_vertex_type> vParent (boost::num_vertices(this->m_graph));
std::vector<uint32_t> vRank (vParent.size());
disjoint_set_type::SubsetHelper<graph_vertex_type, uint32_t> gp (vParent, vRank);
struct blockrec const& block = X.blocks[1];
assert_msg(block.blockcategory == MATRIX,
"mismatch of block category");
for (int32_t i = 1; i <= block.blocksize; ++i)
for (int32_t j = i+1; j <= block.blocksize; ++j)
{
double ent = block.data.mat[ijtok(i, j, block.blocksize)];
if (ent > m_rounding_ub)
disjoint_set_type::union_set(gp, i-1, j-1);
}
uint32_t mg_count = 0;
for (uint32_t i = 0, ie = vParent.size(); i != ie; ++i)
if (vParent[i] == i)
vG2MG[i] = mg_count++;
for (uint32_t i = 0, ie = vParent.size(); i != ie; ++i)
if (vParent[i] != i)
vG2MG[i] = vG2MG.at(disjoint_set_type::find_set(gp, i));
#ifdef DEBUG_SDPCOLORING
assert(mg_count == disjoint_set_type::count_sets(gp));
#endif
graph_type mg (mg_count);
edge_iterator_type ei, eie;
for (boost::tie(ei, eie) = boost::edges(this->m_graph); ei != eie; ++ei)
{
graph_edge_type const& e = *ei;
graph_vertex_type s = boost::source(e, this->m_graph);
graph_vertex_type t = boost::target(e, this->m_graph);
graph_vertex_type ms = vG2MG.at(s);
graph_vertex_type mt = vG2MG.at(t);
std::pair<graph_edge_type, bool> me = boost::edge(ms, mt, mg);
edge_weight_type w = (this->edge_weight(e) >= 0)? 1 : -this->stitch_weight();
if (me.second)
w += boost::get(boost::edge_weight, mg, me.first);
else
me = boost::add_edge(ms, mt, mg);
boost::put(boost::edge_weight, mg, me.first, w);
#ifdef DEBUG_SDPCOLORING
assert(boost::get(boost::edge_weight, mg, me.first) != 0);
#endif
}
#ifdef DEBUG_SDPCOLORING
this->check_edge_weight(mg, this->stitch_weight()/10, boost::num_edges(this->m_graph));
#endif
std::vector<int8_t> vMColor (mg_count, -1);
coloring_merged_graph(mg, vMColor);
vertex_iterator_type vi, vie;
for (boost::tie(vi, vie) = boost::vertices(this->m_graph); vi != vie; ++vi)
{
graph_vertex_type v = *vi;
this->m_vColor[v] = vMColor.at(vG2MG.at(v));
#ifdef DEBUG_SDPCOLORING
assert(this->m_vColor[v] >= 0 && this->m_vColor[v] < this->color_num());
#endif
}
}
template <typename GraphType>
void SDPColoringCsdp<GraphType>::coloring_merged_graph(graph_type const& mg, std::vector<int8_t>& vMColor) const
{
uint32_t num_vertices = boost::num_vertices(mg);
if (num_vertices <= max_backtrack_num_vertices || num_vertices == boost::num_vertices(this->m_graph))
coloring_algos(mg, vMColor);
else
{
typedef GraphSimplification<graph_type> graph_simplification_type;
graph_simplification_type gs (mg, this->color_num());
gs.simplify(graph_simplification_type::HIDE_SMALL_DEGREE);
std::vector<std::vector<int8_t> > mSubColor (gs.num_component());
std::vector<std::vector<graph_vertex_type> > mSimpl2Orig (gs.num_component());
for (uint32_t sub_comp_id = 0; sub_comp_id < gs.num_component(); ++sub_comp_id)
{
graph_type sg;
std::vector<int8_t>& vSubColor = mSubColor[sub_comp_id];
std::vector<graph_vertex_type>& vSimpl2Orig = mSimpl2Orig[sub_comp_id];
gs.simplified_graph_component(sub_comp_id, sg, vSimpl2Orig);
vSubColor.assign(boost::num_vertices(sg), -1);
#ifdef DEBUG_SDPCOLORING
this->
write_graph(
"initial_merged_graph", sg, vSubColor);
#endif
coloring_algos(sg, vSubColor);
#ifdef DEBUG_SDPCOLORING
this->
write_graph(
"final_merged_graph", sg, vSubColor);
#endif
}
gs.recover(vMColor, mSubColor, mSimpl2Orig);
gs.recover_hide_small_degree(vMColor);
}
}
template <typename GraphType>
void SDPColoringCsdp<GraphType>::coloring_algos(graph_type const& g, std::vector<int8_t>& vColor) const
{
if (boost::num_vertices(g) <= max_backtrack_num_vertices)
coloring_by_backtrack(g, vColor);
else
coloring_by_FM(g, vColor);
}
template <typename GraphType>
void SDPColoringCsdp<GraphType>::coloring_by_backtrack(SDPColoringCsdp<GraphType>::graph_type const& mg, std::vector<int8_t>& vColor) const
{
BacktrackColoring<graph_type> bc (mg);
bc.stitch_weight(1);
bc.color_num(this->color_num());
bc();
for (uint32_t i = 0, ie = vColor.size(); i != ie; ++i)
vColor[i] = bc.color(i);
}
template <typename GraphType>
void SDPColoringCsdp<GraphType>::coloring_by_FM(SDPColoringCsdp<GraphType>::graph_type const& mg, std::vector<int8_t>& vColor) const
{
fmp.set_partitions(vColor.begin(), vColor.end());
fmp();
for (uint32_t i = 0, ie = vColor.size(); i != ie; ++i)
vColor[i] = fmp.partition(i);
}
template <typename GraphType>
void SDPColoringCsdp<GraphType>::write_sdp_sol(std::string const& filename, struct blockmatrix const& X) const
{
uint32_t matrix_size = 0;
for (int32_t blk = 1; blk <= X.nblocks; ++blk)
matrix_size += X.blocks[blk].blocksize;
std::vector<std::vector<double> > mSol (matrix_size, std::vector<double>(matrix_size, 0.0));
int32_t index_offset = 0;
for (int32_t blk = 1; blk <= X.nblocks; ++blk)
{
switch (X.blocks[blk].blockcategory)
{
case DIAG:
for (int32_t i = 1; i <= X.blocks[blk].blocksize; ++i)
{
double ent = X.blocks[blk].data.vec[i];
if (ent != 0.0)
mSol[index_offset+i-1][index_offset+i-1] = ent;
};
break;
case MATRIX:
for (int32_t i = 1; i <= X.blocks[blk].blocksize; ++i)
for (int32_t j = i; j <= X.blocks[blk].blocksize; ++j)
{
double ent = X.blocks[blk].data.mat[ijtok(i, j, X.blocks[blk].blocksize)];
if (ent != 0.0)
mSol[index_offset+i-1][index_offset+j-1] = ent;
};
break;
case PACKEDMATRIX:
}
index_offset += X.blocks[blk].blocksize;
}
std::ofstream out (filename.c_str());
assert_msg(out.good(),
"failed to open file " << filename <<
" for write");
for (std::vector<std::vector<double> >::const_iterator it1 = mSol.begin(), it1e = mSol.end(); it1 != it1e; ++it1)
{
const char* prefix = "";
for (std::vector<double>::const_iterator it2 = it1->begin(), it2e = it1->end(); it2 != it2e; ++it2)
{
out << prefix << *it2;
prefix = ",";
}
out << std::endl;
}
out.close();
}
template <typename GraphType>
void SDPColoringCsdp<GraphType>::print_blockrec(const char* label, blockrec const& block) const
{
printf("%s", label);
if (block.blockcategory == MATRIX)
{
printf("[M]: ");
for (int32_t i = 0, ie = block.blocksize*block.blocksize; i != ie; ++i)
printf("%g ", block.data.mat[i]);
}
else if (block.blockcategory == DIAG)
{
printf("[V]: ");
for (int32_t i = 0; i != block.blocksize; ++i)
printf("%g ", block.data.vec[i+1]);
}
printf("\n");
}
}
}
}
#endif
#include <iostream>
#include <fstream>
#include <string>
#include <boost/graph/graphviz.hpp>
#include <boost/graph/graph_utility.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/undirected_graph.hpp>
#include <boost/graph/erdos_renyi_generator.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/graph/random.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <boost/version.hpp>
#if BOOST_VERSION <= 14601
#include <boost/graph/detail/is_same.hpp>
#else
#include <boost/type_traits/is_same.hpp>
#endif
using std::cout;
using std::endl;
using std::ifstream;
using std::ofstream;
using std::string;
using std::pair;
typedef adjacency_list<vecS, vecS, undirectedS,
property<vertex_index_t, std::size_t, property<vertex_color_t, int> >,
property<edge_index_t, std::size_t, property<edge_weight_t, double> >,
property<graph_name_t, string> > graph_type;
typedef subgraph<graph_type> subgraph_type;
typedef property<vertex_index_t, std::size_t> VertexId;
typedef property<edge_index_t, std::size_t> EdgeID;
typedef graph_traits<graph_type>::vertex_descriptor vertex_descriptor;
typedef graph_traits<graph_type>::edge_descriptor edge_descriptor;
typedef property_map<graph_type, edge_weight_t>::type edge_weight_map_type;
typedef property_map<graph_type, vertex_color_t>::type vertex_color_map_type;
{
graph_type g;
vertex_descriptor a = boost::add_vertex(g);
vertex_descriptor b = boost::add_vertex(g);
vertex_descriptor c = boost::add_vertex(g);
vertex_descriptor d = boost::add_vertex(g);
vertex_descriptor e = boost::add_vertex(g);
boost::add_edge(a, b, g);
boost::add_edge(a, c, g);
boost::add_edge(a, d, g);
boost::add_edge(a, e, g);
boost::add_edge(b, c, g);
boost::add_edge(b, e, g);
boost::add_edge(c, d, g);
boost::add_edge(d, e, g);
BOOST_AUTO(edge_weight_map, get(edge_weight, g));
graph_traits<graph_type>::edge_iterator eit, eit_end;
for (tie(eit, eit_end) = edges(g); eit != eit_end; ++eit)
{
if (source(*eit, g) != a || target(*eit, g) != d)
edge_weight_map[*eit] = 1;
else
edge_weight_map[*eit] = -1;
}
lc.stitch_weight(0.1);
return lc();
}
{
mt19937 gen;
graph_type g;
int N = 40;
std::vector<vertex_descriptor> vertex_set;
std::vector< std::pair<vertex_descriptor, vertex_descriptor> > edge_set;
generate_random_graph(g, N, N * 2, gen,
std::back_inserter(vertex_set),
std::back_inserter(edge_set));
BOOST_AUTO(edge_weight_map, get(edge_weight, g));
unsigned int i = 0;
graph_traits<graph_type>::edge_iterator eit, eit_end;
for (tie(eit, eit_end) = edges(g); eit != eit_end; ++eit, ++i)
{
#if 1
if (i%10 == 0)
edge_weight_map[*eit] = -1;
else
#endif
edge_weight_map[*eit] = 1;
}
lc.stitch_weight(0.1);
return lc();
}
{
ifstream in (filename.c_str());
typedef adjacency_list<vecS, vecS, undirectedS,
property<vertex_index_t, std::size_t, property<vertex_color_t, int, property<vertex_name_t, std::size_t> > >,
property<edge_index_t, std::size_t, property<edge_weight_t, int> >,
property<graph_name_t, string> > tmp_graph_type;
tmp_graph_type tmpg;
dynamic_properties tmpdp;
tmpdp.property("node_id", get(vertex_name, tmpg));
tmpdp.property("label", get(vertex_name, tmpg));
tmpdp.property("weight", get(edge_weight, tmpg));
tmpdp.property("label", get(edge_weight, tmpg));
assert(read_graphviz(in, tmpg, tmpdp, "node_id"));
graph_type g (num_vertices(tmpg));
graph_traits<tmp_graph_type>::vertex_iterator vit, vit_end;
for (tie(vit, vit_end) = vertices(tmpg); vit != vit_end; ++vit)
{
size_t name = get(vertex_name, tmpg, *vit);
int color = get(vertex_color, tmpg, *vit);
put(vertex_color, g, (vertex_descriptor)name, color);
}
graph_traits<tmp_graph_type>::edge_iterator eit, eit_end;
for (tie(eit, eit_end) = edges(tmpg); eit != eit_end; ++eit)
{
size_t s_name = get(vertex_name, tmpg, source(*eit, tmpg));
size_t t_name = get(vertex_name, tmpg, target(*eit, tmpg));
pair<edge_descriptor, bool> pe = add_edge(s_name, t_name, g);
assert(pe.second);
int weight = get(edge_weight, g, *eit);
put(edge_weight, g, pe.first, weight);
}
#ifdef DEBUG_SDPCOLORING
dynamic_properties dp;
dp.property("id", get(vertex_index, g));
dp.property("node_id", get(vertex_index, g));
dp.property("label", get(vertex_index, g));
dp.property("weight", get(edge_weight, g));
dp.property("label", get(edge_weight, g));
ofstream out ("graph_init.gv");
write_graphviz_dp(out, g, dp, string("id"));
out.close();
system("dot -Tpdf graph_init.gv -o graph_init.pdf");
#endif
lc.stitch_weight(0.1);
double cost = lc();
in.close();
return cost;
}
int main(
int argc,
char** argv)
{
double cost;
if (argc < 2)
{
}
cout << "cost = " << cost << endl;
return 0;
}
Compiling and running commands (assuming LIMBO_DIR, and GUROBI_HOME are well defined).
Compiling and running commands (assuming LIMBO_DIR, and LPSOLVE_DIR are well defined).