Limbo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Limbo.Solvers

Table of Contents

Introduction

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.

Examples

Primal Min-Cost Flow Solvers

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.

See documented version: test/solvers/test_MinCostFlow.cpp

#include <iostream>
#include <limbo/solvers/MinCostFlow.h>
void test(std::string const& filename, int alg)
{
typedef limbo::solvers::LinearModel<double, int> model_type;
model_type optModel;
optModel.read(filename);
// print problem
optModel.print(std::cout);
// solve
limbo::solvers::MinCostFlowSolver<double, int>* minCostFlowSolver = NULL;
switch (alg)
{
case 0:
minCostFlowSolver = new limbo::solvers::CostScaling<double, int>();
break;
case 1:
minCostFlowSolver = new limbo::solvers::CapacityScaling<double, int>();
break;
case 2:
minCostFlowSolver = new limbo::solvers::NetworkSimplex<double, int>();
break;
case 3:
default:
minCostFlowSolver = new limbo::solvers::CycleCanceling<double, int>();
break;
}
limbo::solvers::MinCostFlow<double, int> solver (&optModel);
limbo::solvers::SolverProperty status = solver(minCostFlowSolver);
//limbo::solvers::SolverProperty status = solver();
std::cout << "Problem solved " << limbo::solvers::toString(status) << "\n";
delete minCostFlowSolver;
// print solutions
optModel.printSolution(std::cout);
// print problem
optModel.print(std::cout);
// print graph with solution information
solver.printGraph(true);
}
int main(int argc, char** argv)
{
if (argc > 1)
{
int alg = 0;
if (argc > 2)
alg = atoi(argv[2]);
// test file API
test(argv[1], alg);
}
else
std::cout << "at least 1 argument required\n";
return 0;
}

Compiling and running commands (assuming LIMBO_DIR, BOOST_DIR and LEMON_DIR are well defined). Limbo.Parsers.LpParser is required for limbo::solvers::MinCostFlow to read input files in .lp format.

1 g++ -o test_MinCostFlow test_MinCostFlow.cpp -I $LIMBO_DIR/include -I $BOOST_DIR/include -I $LEMON_DIR/include -L $LEMON_DIR/lib -lemon -L $LIMBO_DIR/lib -llpparser
2 # test primal min-cost flow for linear programming problem
3 ./test_MinCostFlow lpmcf/benchmarks/mcf.lp

Output

1 # debug.lgf
2 @nodes
3 label supply name potential
4 0 0 node_constr_0 -4
5 1 0 node_constr_1 -4
6 2 -1 subg_constr_0 -1
7 3 -1 subg_constr_1 0
8 4 1 weight_constr -4
9 5 1 st -8
10 @arcs
11  label name capacity_lower capacity_upper cost flow
12 0 2 0 x0_0 0 1 5 0
13 0 3 1 x0_1 0 1 4 0
14 1 2 2 x1_0 0 1 3 1
15 1 3 3 x1_1 0 1 7 0
16 5 3 4 xd1 0 1 8 1
17 5 2 5 xd0 0 1 8 0
18 4 0 6 x0 0 1 0 0
19 4 1 7 x1 0 1 0 1
1 Minimize
2 5.437 x0_0 + 4.689 x0_1 + 3.223 x1_0 + 7.654 x1_1
3  + 8.7657 xd1 + 8.32079 xd0
4 
5 Subject To
6 node_constr_0: 1 x0_0 + 1 x0_1 + -1 x0 = 0
7 node_constr_1: 1 x1_0 + 1 x1_1 + -1 x1 = 0
8 subg_constr_0: 1 x0_0 + 1 x1_0 + 1 xd0 = 1
9 subg_constr_1: 1 x0_1 + 1 x1_1 + 1 xd1 = 1
10 weight_constr: 1 x0 + 1 x1 = 1
11 Bounds
12 0 <= x0_0 <= 1
13 0 <= x0_1 <= 1
14 0 <= x1_0 <= 1
15 0 <= x1_1 <= 1
16 0 <= xd1 <= 1
17 0 <= xd0 <= 1
18 0 <= x0 <= 1
19 0 <= x1 <= 1
20 Generals
21 End
22 Problem solved OPTIMAL
23 # Objective 11.9887
24 x0_0 0
25 x0_1 0
26 x1_0 1
27 x1_1 0
28 xd1 1
29 xd0 0
30 x0 0
31 x1 1

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.

Dual Min-Cost Flow Solvers

The analysis and background of using dual min-cost flow to solve linear programming problem can be found in the detailed description of class limbo::solvers::DualMinCostFlow and limbo::solvers::lpmcf::LpDualMcf.

See documented version: test/solvers/test_DualMinCostFlow.cpp

#include <iostream>
#include <limbo/solvers/DualMinCostFlow.h>
void test(std::string const& filename, int alg)
{
typedef limbo::solvers::LinearModel<int, int> model_type;
model_type optModel;
optModel.read(filename);
// print problem
optModel.print(std::cout);
// solve
limbo::solvers::MinCostFlowSolver<int, int>* minCostFlowSolver = NULL;
switch (alg)
{
case 0:
minCostFlowSolver = new limbo::solvers::CostScaling<int, int>();
break;
case 1:
minCostFlowSolver = new limbo::solvers::CapacityScaling<int, int>();
break;
case 2:
minCostFlowSolver = new limbo::solvers::NetworkSimplex<int, int>();
break;
case 3:
default:
minCostFlowSolver = new limbo::solvers::CycleCanceling<int, int>();
break;
}
limbo::solvers::DualMinCostFlow<int, int> solver (&optModel);
limbo::solvers::SolverProperty status = solver(minCostFlowSolver);
//limbo::solvers::SolverProperty status = solver();
std::cout << "Problem solved " << limbo::solvers::toString(status) << "\n";
delete minCostFlowSolver;
// print solutions
optModel.printSolution(std::cout);
// print problem
optModel.print(std::cout);
// print graph with solution information
solver.printGraph(true);
}
int main(int argc, char** argv)
{
if (argc > 1)
{
int alg = 0;
if (argc > 2)
alg = atoi(argv[2]);
// test file API
test(argv[1], alg);
}
else
std::cout << "at least 1 argument required\n";
return 0;
}

Compiling and running commands (assuming LIMBO_DIR, BOOST_DIR and LEMON_DIR are well defined). Limbo.Parsers.LpParser is required for limbo::solvers::DualMinCostFlow to read input files in .lp format.

1 g++ -o test_DualMinCostFlow test_DualMinCostFlow.cpp -I $LIMBO_DIR/include -I $BOOST_DIR/include -I $LEMON_DIR/include -L $LEMON_DIR/lib -lemon -L $LIMBO_DIR/lib -llpparser
2 # test dual min-cost flow for linear programming problem
3 ./test_DualMinCostFlow lpmcf/benchmarks/problem.lp

Output

1 # debug.lgf
2 @nodes
3 label supply name potential
4 0 -4 R1 -2
5 1 -5 L1 -5
6 2 -5 R2 -2
7 3 -1 L2 -2
8 4 3 x2 -3
9 5 6 x1 -7
10 6 3 x3 -4
11 7 3 additional -7
12 @arcs
13  label capacity_upper cost flow
14 5 4 0 3 4 1
15 5 1 1 3 2 2
16 5 0 2 3 2 3
17 4 1 3 3 1 0
18 4 0 4 3 1 1
19 4 3 5 3 1 1
20 4 2 6 3 1 2
21 6 3 7 3 2 0
22 6 2 8 3 2 3
23 0 7 9 3 10 0
24 7 1 10 3 2 0
25 7 1 11 3 2 3
26 2 7 12 3 10 0
27 3 7 13 3 10 0
28 4 7 14 3 10 0
29 5 7 15 3 10 0
30 6 7 16 3 10 0
1 Minimize
2 2 R1 + -2 L1 + 1 R2 + -1 L2
3 
4 
5 Subject To
6 C0: 1 x2 + -1 x1 >= 4
7 C1: -1 L1 + 1 x1 >= -2
8 C2: 1 R1 + -1 x1 >= 2
9 C3: -1 L1 + 1 x2 >= -1
10 C4: 1 R1 + -1 x2 >= 1
11 C5: -1 L2 + 1 x2 >= -1
12 C6: 1 R2 + -1 x2 >= 1
13 C7: -1 L2 + 1 x3 >= -2
14 C8: 1 R2 + -1 x3 >= 2
15 Bounds
16 R1 >= -10
17 2 <= L1 <= 2
18 R2 >= -10
19 L2 >= -10
20 x2 >= -10
21 x1 >= -10
22 x3 >= -10
23 Generals
24 R1
25 L1
26 R2
27 L2
28 x2
29 x1
30 x3
31 End
32 
33 Problem solved OPTIMAL
34 # Objective 6
35 R1 5
36 L1 2
37 R2 5
38 L2 5
39 x2 4
40 x1 0
41 x3 3

Csdp Solver

See documented version: limbo/algorithms/coloring/SDPColoringCsdp.h

#ifndef LIMBO_ALGORITHMS_COLORING_SDPCOLORINGCSDP
#define LIMBO_ALGORITHMS_COLORING_SDPCOLORINGCSDP
#include <limbo/string/String.h>
#include <limbo/algorithms/coloring/GraphSimplification.h>
#include <limbo/algorithms/coloring/Coloring.h>
#include <limbo/algorithms/coloring/BacktrackColoring.h>
#include <limbo/algorithms/partition/FMMultiWay.h>
#include <limbo/containers/DisjointSet.h>
// as the original csdp easy_sdp api is not very flexible to printlevel
// I made small modification to support that
#include <limbo/solvers/api/CsdpEasySdpApi.h>
namespace limbo
{
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;
using typename base_type::ColorNumType;
typedef typename base_type::EdgeHashType edge_hash_type;
struct FMGainCalcType
{
typedef edge_weight_type value_type;
graph_type const& graph;
FMGainCalcType(graph_type const& g) : graph(g) {}
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;
// if v is not partitioned, then any partition will have preassigned large gain
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
// skip unpartitioned vertex
if (pt < 0) continue;
edge_weight_type w = boost::get(boost::edge_weight, graph, *ei);
// assume origp != newp, pt >= 0
gain += (pt == newp)? -w : (pt == origp)? w : 0;
//gain += w * ((vPartition[t] == origp && origp >= 0) - (vPartition[t] == newp));
}
return gain;
}
};
SDPColoringCsdp(graph_type const& g);
virtual ~SDPColoringCsdp() {}
void write_sdp_sol(std::string const& filename, struct blockmatrix const& X) const;
void print_blockrec(const char* label, blockrec const& block) const;
protected:
virtual double coloring();
void construct_objectve_blockrec(blockmatrix& C, int32_t blocknum, int32_t blocksize, blockcat blockcategory) const;
struct sparseblock* construct_constraint_sparseblock(int32_t blocknum, int32_t blocksize, int32_t constraintnum, int32_t entrynum) const;
void set_sparseblock_entry(struct sparseblock& block, int32_t entryid, int32_t i, int32_t j, double value) const;
void round_sol(struct blockmatrix const& X);
void coloring_merged_graph(graph_type const& mg, std::vector<int8_t>& vMColor) const;
void coloring_algos(graph_type const& g, std::vector<int8_t>& vColor) const;
virtual void coloring_by_backtrack(graph_type const& mg, std::vector<int8_t>& vColor) const;
virtual void coloring_by_FM(graph_type const& mg, std::vector<int8_t>& vColor) const;
double m_rounding_lb;
double m_rounding_ub;
const static uint32_t max_backtrack_num_vertices = 7;
};
template <typename GraphType>
SDPColoringCsdp<GraphType>::SDPColoringCsdp(SDPColoringCsdp<GraphType>::graph_type const& g)
: base_type(g)
{
m_rounding_lb = -0.4;
m_rounding_ub = 0.9;
}
template <typename GraphType>
double SDPColoringCsdp<GraphType>::coloring()
{
// Since Csdp is written in C, the api here is also in C
// Please refer to the documation of Csdp for different notations
// basically, X is primal variables, C, b, constraints and pobj are all for primal
// y, Z, and dobj are for dual problem
//
// Csdp has very complex storage structure for matrix
// I still do not have a full understanding about the block concept, especially blocks.blocksize
// with some reverse engineering, for the coloring problem here, matrices in C, b, and constraints mainly consists of 2 blocks
// the first block is for vertex variables, and the second block is for slack variables introduced to resolve '>=' operators in the constraints
assert_msg(!this->has_precolored(), "SDP coloring does not support precolored layout yet");
// The problem and solution data.
struct blockmatrix C; // objective matrix
double *b; // right hand side of constraints
struct constraintmatrix *constraints; // constraint matrices
// Storage for the initial and final solutions.
struct blockmatrix X,Z;
double *y;
double pobj,dobj;
// iterators used to traverse through the graph
vertex_iterator_type vi, vie;
edge_iterator_type ei, eie;
// compute total number of vertices and edges
uint32_t num_vertices = boost::num_vertices(this->m_graph);
uint32_t num_edges = boost::num_edges(this->m_graph);
// compute total number of conflict edges and stitch edges
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) // conflict edge
num_conflict_edges += 1;
else // stitch edge
num_stitch_edges += 1;
}
assert_msg(num_edges > 0 && num_conflict_edges > 0, "no edges or conflict edges found, no need to solve SDP");
// compute total number of variables and constraints
uint32_t num_variables = num_vertices+num_conflict_edges;
uint32_t num_constraints = num_conflict_edges+num_vertices;
// setup blockmatrix C
C.nblocks = 2;
C.blocks = (struct blockrec *)malloc((C.nblocks+1)*sizeof(struct blockrec));
assert_msg(C.blocks, "Couldn't allocate storage for C");
// C.blocks[0] is not used according to the example of Csdp
// block 1 for vertex variables
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);
// 1 for conflict edge, -alpha for stitch edge
// add unary negative operator, because Csdp solves maximization problem
// but we are solving minimization problem
edge_weight_type alpha = (this->edge_weight(*ei) >= 0)? -1 : this->stitch_weight();
// variable starts from 1 instead of 0 in Csdp
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;
}
// block 2 for slack variables
// this block is all 0s, so we use diagonal format to represent
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
// setup right hand side of constraints b
// the order is first for conflict edges and then for vertices
// the order matters for constraint matrices
b = (double *)malloc((num_constraints+1)*sizeof(double));
assert_msg(b, "Failed to allocate storage for right hand side of constraints b");
// -1/(k-1) according to Bei Yu's DAC2014 paper
// consider in the constraints, xij+xji >= beta, so beta should be -2/(k-1)
double beta = -2.0/(this->color_num()-1.0); // right hand side of constraints for conflict edges
for (uint32_t i = 1, ie = num_constraints+1; i != ie; ++i)
{
if (i <= num_conflict_edges) // slack for each conflict edge, xij >= -0.5
b[i] = beta;
else // slack for each vertex, xii = 1
b[i] = 1;
}
// setup constraint matrix constraints
// the order should be the same as right hand side b
constraints=(struct constraintmatrix *)malloc((num_constraints+1)*sizeof(struct constraintmatrix));
assert_msg(constraints, "Failed to allocate storage for constraints");
// for conflict edges, xij
uint32_t cnt = 1;
for (boost::tie(ei, eie) = boost::edges(this->m_graph); ei != eie; ++ei)
{
if (this->edge_weight(*ei) >= 0) // conflict edge
{
graph_vertex_type s = boost::source(*ei, this->m_graph);
graph_vertex_type t = boost::target(*ei, this->m_graph);
if (s > t) // due to symmetry, only need to create constraint matrices for upper-matrix
std::swap(s, t);
// variable starts from 1 instead of 0 in Csdp
s += 1; t += 1;
struct constraintmatrix& constr = constraints[cnt];
// Terminate the linked list with a NULL pointer.
constr.blocks = NULL;
// inverse order to initialize blocks, because linked list will reverse the order
// first set block 2 for diagonal values and then block 1 for upper-matrix values
for (uint32_t i = 2; i != 0; --i)
{
struct sparseblock* blockptr;
if (i == 1) // block 1, vertex variables
{
blockptr = construct_constraint_sparseblock(i, num_vertices, cnt, 1);
set_sparseblock_entry(*blockptr, 1, s, t, 1);
}
else // block 2, slack variables
{
blockptr = construct_constraint_sparseblock(i, num_conflict_edges, cnt, 1);
set_sparseblock_entry(*blockptr, 1, cnt, cnt, -1);
}
// insert block to linked list
blockptr->next = constr.blocks;
constr.blocks = blockptr;
}
++cnt;
}
}
// for vertices, xii
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];
// Terminate the linked list with a NULL pointer.
constr.blocks = NULL;
struct sparseblock* blockptr = construct_constraint_sparseblock(1, num_vertices, cnt, 1);
set_sparseblock_entry(*blockptr, 1, v, v, 1);
// insert block to linked list
blockptr->next = constr.blocks;
constr.blocks = blockptr;
++cnt;
}
#ifdef DEBUG_SDPCOLORING
write_prob((char*)"problem.sdpa", num_variables, num_constraints, C, b, constraints);
#endif
// after all configuration ready
// start solving sdp
// set initial solution
initsoln(num_variables, num_constraints, C, b, constraints, &X, &y, &Z);
// use default params
// only set printlevel to zero
struct paramstruc params;
int printlevel;
initparams(&params, &printlevel);
//#ifndef DEBUG_SDPCOLORING
printlevel = 0;
//#endif
// A return code for the call to easy_sdp().
// solve sdp
// objective value is (dobj+pobj)/2
//int ret = easy_sdp(num_variables, num_constraints, C, b, constraints, 0.0, &X, &y, &Z, &pobj, &dobj);
int ret = limbo::solvers::easy_sdp_ext<int>(num_variables, num_constraints, C, b, constraints, 0.0, &X, &y, &Z, &pobj, &dobj, params, printlevel);
assert_msg(ret == 0, "SDP failed");
// round result to get colors
round_sol(X);
// Free storage allocated for the problem and return.
free_prob(num_variables, num_constraints, C, b, constraints, X, y, Z);
#ifdef DEBUG_LPCOLORING
this->write_graph("final_output");
#endif
// return objective value
//return (dobj+pobj)/2;
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");
// initialize to all 0s
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");
// initialize to all 0s
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
// merge vertices with SDP solution with disjoint set
std::vector<graph_vertex_type> vParent (boost::num_vertices(this->m_graph));
std::vector<uint32_t> vRank (vParent.size());
typedef limbo::containers::DisjointSet disjoint_set_type;
disjoint_set_type::SubsetHelper<graph_vertex_type, uint32_t> gp (vParent, vRank);
// check SDP solution in X
// we are only interested in block 1
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) // merge vertices if SDP solution rounded to 1.0
disjoint_set_type::union_set(gp, i-1, j-1); // Csdp array starts from 1 instead of 0
}
// construct merged graph
// for vertices in merged graph
std::vector<graph_vertex_type> vG2MG (vParent.size(), std::numeric_limits<graph_vertex_type>::max()); // mapping from graph to merged graph
uint32_t mg_count = 0; // count number of vertices in merged graph
for (uint32_t i = 0, ie = vParent.size(); i != ie; ++i) // check subset elements and compute mg_count
if (vParent[i] == i)
vG2MG[i] = mg_count++;
for (uint32_t i = 0, ie = vParent.size(); i != ie; ++i) // check other elements
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); // merged graph
// for edges in merged graph
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);
// need to consider if this setting is still reasonable when stitch is on
edge_weight_type w = (this->edge_weight(e) >= 0)? 1 : -this->stitch_weight();
if (me.second) // already exist, update weight
w += boost::get(boost::edge_weight, mg, me.first);
else // not exist, add edge
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->print_edge_weight(this->m_graph);
//this->check_edge_weight(this->m_graph, this->stitch_weight()/10, 4);
//this->print_edge_weight(mg);
this->check_edge_weight(mg, this->stitch_weight()/10, boost::num_edges(this->m_graph));
#endif
// coloring for merged graph
std::vector<int8_t> vMColor (mg_count, -1); // coloring solution for merged graph
coloring_merged_graph(mg, vMColor);
// apply coloring solution from merged graph to graph
// first, map colors to subsets
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 small number of vertices or no vertex merged, no need to simplify graph
if (num_vertices <= max_backtrack_num_vertices || num_vertices == boost::num_vertices(this->m_graph))
coloring_algos(mg, vMColor);
else
{
// simplify merged graph
typedef GraphSimplification<graph_type> graph_simplification_type;
graph_simplification_type gs (mg, this->color_num());
gs.simplify(graph_simplification_type::HIDE_SMALL_DEGREE);
// in order to recover color from articulation points
// we have to record all components and mappings
// but graph is not necessary
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
// solve coloring
coloring_algos(sg, vSubColor);
#ifdef DEBUG_SDPCOLORING
this->write_graph("final_merged_graph", sg, vSubColor);
#endif
}
// recover color assignment according to the simplification level set previously
// HIDE_SMALL_DEGREE needs to be recovered manually for density balancing
gs.recover(vMColor, mSubColor, mSimpl2Orig);
// recover colors for simplified vertices without balanced assignment
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
{
// currently backtrack coloring is used
// TO DO: add faster coloring approach like FM partition based
BacktrackColoring<graph_type> bc (mg);
bc.stitch_weight(1); // already scaled in edge weights
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
{
limbo::algorithms::partition::FMMultiWay<FMGainCalcType> fmp (FMGainCalcType(mg), boost::num_vertices(mg), this->color_num());
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
{
// compute dimensions of matrix
uint32_t matrix_size = 0;
for (int32_t blk = 1; blk <= X.nblocks; ++blk)
matrix_size += X.blocks[blk].blocksize;
// collect data from X and store to mSol
std::vector<std::vector<double> > mSol (matrix_size, std::vector<double>(matrix_size, 0.0));
// as i and j starts from 1, set index_offset to -1
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:
default: assert_msg(0, "Invalid Block Type");
}
index_offset += X.blocks[blk].blocksize;
}
// write to file
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]: "); // show it is a matrix
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]: "); // show it is a vector
for (int32_t i = 0; i != block.blocksize; ++i)
printf("%g ", block.data.vec[i+1]);
}
printf("\n");
}
} // namespace coloring
} // namespace algorithms
} // namespace limbo
#endif

See documented version: test/algorithms/test_SDPColoring.cpp

#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 <limbo/algorithms/coloring/SDPColoringCsdp.h>
#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;
using namespace boost;
// do not use setS, it does not compile for subgraph
// do not use custom property tags, it does not compile for most utilities
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;
double simple_graph()
{
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;
}
//test relaxed LP based coloring
limbo::algorithms::coloring::SDPColoringCsdp<graph_type> lc (g);
lc.stitch_weight(0.1);
// THREE or FOUR
lc.color_num(limbo::algorithms::coloring::SDPColoringCsdp<graph_type>::THREE);
return lc();
}
double random_graph()
{
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) // generate stitch
edge_weight_map[*eit] = -1;
else // generate conflict
#endif
edge_weight_map[*eit] = 1;
}
//test relaxed LP based coloring
limbo::algorithms::coloring::SDPColoringCsdp<graph_type> lc (g);
lc.stitch_weight(0.1);
// THREE or FOUR
lc.color_num(limbo::algorithms::coloring::SDPColoringCsdp<graph_type>::THREE);
return lc();
}
double real_graph(string const& filename)
{
ifstream in (filename.c_str());
// the graphviz reader in boost cannot specify vertex_index_t
// I have to create a temporary graph and then construct the real graph
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"));
// real graph
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
//test relaxed LP based coloring
limbo::algorithms::coloring::SDPColoringCsdp<graph_type> lc (g);
lc.stitch_weight(0.1);
// THREE or FOUR
lc.color_num(limbo::algorithms::coloring::SDPColoringCsdp<graph_type>::THREE);
double cost = lc();
in.close();
return cost;
}
int main(int argc, char** argv)
{
double cost;
if (argc < 2)
{
cost = simple_graph();
//cost = random_graph();
}
else cost = real_graph(argv[1]);
cout << "cost = " << cost << endl;
return 0;
}

Compiling and running commands (assuming LIMBO_DIR, BOOST_DIR and LEMON_DIR are well defined). Limbo.Parsers.LpParser is required for limbo::solvers::lpmcf::LpDualMcf to read input files in .lp format.

1 g++ -o test_lpmcf compare.cpp -I $LIMBO_DIR/include -I $BOOST_DIR/include -I $LEMON_DIR/include -L $LEMON_DIR/lib -lemon -L $LIMBO_DIR/lib -llpparser
2 # test: min-cost flow for network graph
3 ./test_lpmcf benchmarks/graph.lgf

Generalized Linear Programming API

A generalized API to solve problem described by limbo::solvers::LinearModel with Gurobi and lpsolve solvers.

See documented version: test/solvers/test_GurobiApi.cpp and test/solvers/test_LPSolveApi.cpp

For Gurobi

#include <iostream>
#include <limbo/solvers/api/GurobiApi.h>
int main()
{
// ILP model
typedef limbo::solvers::LinearModel<double, long> model_type;
model_type optModel;
// create variables
model_type::variable_type var1 = optModel.addVariable(0, 1, limbo::solvers::CONTINUOUS, "x1");
model_type::variable_type var2 = optModel.addVariable(0, 1, limbo::solvers::CONTINUOUS, "x2");
model_type::variable_type var3 = optModel.addVariable(0, 1, limbo::solvers::CONTINUOUS, "x3");
model_type::variable_type var4 = optModel.addVariable(0, 1, limbo::solvers::CONTINUOUS, "x4");
// create objective
optModel.setObjective(var1+var2+var3+var4);
optModel.setOptimizeType(limbo::solvers::MIN);
// create constraints
optModel.addConstraint(var1 - var2 >= 0.5, "c1");
optModel.addConstraint(var4 - var3 >= 0.1, "c2");
optModel.addConstraint(var2 - var3 >= 0.2, "c3");
// solve by Gurobi
typedef limbo::solvers::GurobiLinearApi<model_type::coefficient_value_type, model_type::variable_value_type> solver_type;
solver_type solver (&optModel);
limbo::solvers::GurobiParameters gurobiParams;
gurobiParams.setNumThreads(1);
gurobiParams.setOutputFlag(1);
limbo::solvers::SolverProperty optStatus = solver(&gurobiParams);
std::cout << "optStatus = " << optStatus << std::endl;
return 0;
}

Compiling and running commands (assuming LIMBO_DIR, and GUROBI_HOME are well defined).

1 g++ -o test_GurobiApi test_GurobiApi.cpp -I $LIMBO_DIR/include -I $GUROBI_HOME/include -L $GUROBI_HOME/lib -lgurobi60
2 # test Gurobi API for linear programming problem
3 ./test_GurobiApi

Output

1 Optimize a model with 3 rows, 4 columns and 6 nonzeros
2 Coefficient statistics:
3  Matrix range [1e+00, 1e+00]
4  Objective range [1e+00, 1e+00]
5  Bounds range [1e+00, 1e+00]
6  RHS range [1e-01, 5e-01]
7 Presolve removed 3 rows and 4 columns
8 Presolve time: 0.00s
9 Presolve: All rows and columns removed
10 Iteration Objective Primal Inf. Dual Inf. Time
11  0 1.0000000e+00 0.000000e+00 0.000000e+00 0s
12 
13 Solved in 0 iterations and 0.00 seconds
14 Optimal objective 1.000000000e+00
15 optStatus = 5

For lpsolve (same model construction as Gurobi)

#include <iostream>
#include <limbo/solvers/api/LPSolveApi.h>
int main()
{
// ILP model
typedef limbo::solvers::LinearModel<double, long> model_type;
model_type optModel;
// create variables
model_type::variable_type var1 = optModel.addVariable(0, 1, limbo::solvers::CONTINUOUS, "x1");
model_type::variable_type var2 = optModel.addVariable(0, 1, limbo::solvers::CONTINUOUS, "x2");
model_type::variable_type var3 = optModel.addVariable(0, 1, limbo::solvers::CONTINUOUS, "x3");
model_type::variable_type var4 = optModel.addVariable(0, 1, limbo::solvers::CONTINUOUS, "x4");
// create objective
optModel.setObjective(var1+var2+var3+var4);
optModel.setOptimizeType(limbo::solvers::MIN);
// create constraints
optModel.addConstraint(var1 - var2 >= 0.5, "c1");
optModel.addConstraint(var4 - var3 >= 0.1, "c2");
optModel.addConstraint(var2 - var3 >= 0.2, "c3");
// solve by LPSolve
typedef limbo::solvers::LPSolveLinearApi<model_type::coefficient_value_type, model_type::variable_value_type> solver_type;
solver_type solver (&optModel);
limbo::solvers::LPSolveParameters lpsolveParams;
lpsolveParams.setVerbose(FULL);
limbo::solvers::SolverProperty optStatus = solver(&lpsolveParams);
std::cout << "optStatus = " << optStatus << std::endl;
return 0;
}

Compiling and running commands (assuming LIMBO_DIR, and LPSOLVE_DIR are well defined).

1 g++ -o test_LPSolveApi test_LPSolveApi.cpp -I $LIMBO_DIR/include -I $LPSOLVE_DIR -L $LPSOLVE_DIR -llpsolve55
2 # test lpsolve API for linear programming problem
3 ./test_LPSolveApi

Output

1 Model name: 'LPSolveLinearApi' - run #1
2 Objective: Minimize(R0)
3 
4 SUBMITTED
5 Model size: 3 constraints, 4 variables, 6 non-zeros.
6 Sets: 0 GUB, 0 SOS.
7 
8 PRESOLVE Elimination loops performed.......... O2:M2:I2
9  3 bounds............................... TIGHTENED.
10  [ +0 < Z < +4 ]
11 
12 REDUCED
13 Model size: 3 constraints, 4 variables, 6 non-zeros.
14 Sets: 0 GUB, 0 SOS.
15 Row-types: 0 LE, 3 GE, 0 EQ.
16 
17 
18 CONSTRAINT CLASSES
19 General REAL 3
20 
21 Using DUAL simplex for phase 1 and PRIMAL simplex for phase 2.
22 The primal and dual simplex pricing strategy set to 'Devex'.
23 
24 Objective value -0.9 at iter 2.
25 Optimal solution with dual simplex at iter 3.
26 
27 Primal objective:
28 
29  Column name Value Objective Min Max
30  --------------------------------------------------------------------------
31  x1 1 0.7 0 1e+30
32  x2 1 0.2 -1 1e+30
33  x3 1 0 -3 1e+30
34  x4 1 0.1 0 1e+30
35 
36 Primal variables:
37 
38  Column name Value Slack Min Max
39  --------------------------------------------------------------------------
40  x1 0.7 0 -1e+30 1e+30
41  x2 0.2 0 -1e+30 1e+30
42  x3 0 4 -0.1 0.3
43  x4 0.1 0 -1e+30 1e+30
44 
45 Dual variables:
46 
47  Row name Value Slack Min Max
48  --------------------------------------------------------------------------
49  c1 1 0.5 -0.2 0.8
50  c2 1 0.1 0 1
51  c3 2 0.2 0 0.5
52 
53 
54 Optimal solution 1 after 3 iter.
55 
56 Relative numeric accuracy ||*|| = 3.70074e-17
57 
58  MEMO: lp_solve version 5.5.2.5 for 64 bit OS, with 64 bit REAL variables.
59  In the total iteration count 3, 0 (0.0%) were bound flips.
60  There were 0 refactorizations, 0 triggered by time and 0 by density.
61  ... on average 3.0 major pivots per refactorization.
62  The largest [LUSOL v2.2.1.0] fact(B) had 4 NZ entries, 1.0x largest basis.
63  The constraint matrix inf-norm is 1, with a dynamic range of 1.
64  Time to load data was 0.000 seconds, presolve used 0.000 seconds,
65  ... 0.001 seconds in simplex solver, in total 0.001 seconds.
66 optStatus = 5

All Examples

Possible dependencies: Boost, Lemon.

References

Generated on Wed Jul 3 2019 09:58:21 for Limbo by   doxygen 1.8.8