Commit 11c1a153 authored by Pavel Smirnov's avatar Pavel Smirnov

added multithreading

parent 9652a5f7
......@@ -3,7 +3,6 @@ MINPSC_PATH=tools/minpsc
all:
make build
make test
rebuild:
rm -rf $(GRID_BUILDER_PATH)/CMakeCache.txt
......@@ -16,7 +15,7 @@ rebuild:
make build
build: grid_builder minpsc
build: minpsc
grid_builder:
cmake -B'$(GRID_BUILDER_PATH)' -H'$(GRID_BUILDER_PATH)'
......
......@@ -4,9 +4,6 @@
#include <unordered_map>
#include <omp.h>
#include <ilcplex/ilocplex.h>
using namespace std;
static wval CalcSolutionUBound(const Graph& g) {
......@@ -17,7 +14,7 @@ static wval CalcSolutionUBound(const Graph& g) {
return solution_ubound;
}
static Graph BuildReducedGraph(const Graph& g) {
static void BuildReducedGraph(Graph& g) {
wval solution_ubound = CalcSolutionUBound(g);
vector<wval> lbounds = CalcLB(g);
wval solution_lbound = 0;
......@@ -37,19 +34,18 @@ static Graph BuildReducedGraph(const Graph& g) {
cout << "(" << edges_in_graph << "/" << n*(n-1)/2 << " edges left) " << endl;
for (size_t v = 0; v < n; ++v)
rg.SortEdges(v);
return rg;
g = rg;
}
static pair<Graph, Graph> BuildLBGraph(const Graph& rg) {
size_t n = rg.Size();
Graph ga(n);
Graph gainv(n);
for (size_t v = 0; v < n; ++v) {
const Edge & e = rg.GetEdges(v)[0];
ga.AddOrientedEdge(v, e.to, wval(0));
gainv.AddOrientedEdge(e.to, v, wval(0));
static pair<Graph, Graph> BuildLBGraph(const Graph& reduced_graph) {
Graph graph_a(reduced_graph.Size());
Graph graph_a_inv(reduced_graph.Size());
for (size_t v = 0; v < reduced_graph.Size(); ++v) {
const Edge& e = reduced_graph.GetEdges(v)[0];
graph_a.AddOrientedEdge(v, e.to, wval(0));
graph_a_inv.AddOrientedEdge(e.to, v, wval(0));
}
return {ga, gainv};
return {graph_a, graph_a_inv};
}
static unordered_map<size_t, IloNumVar> AddObjectiveAndVariablesY(
......@@ -246,40 +242,29 @@ static void AddAllOptimizations(
AddOptimization25(env, model, rg, yMapping);
}
wval CplexEx1Connect(const Graph& g, size_t s, unordered_map<string, size_t>& info) {
if (g.Size() <= 1) {
return 0;
}
Graph rg = BuildReducedGraph(g);
size_t edgesLeft = 0;
for (size_t v = 0; v < rg.Size(); ++v) {
edgesLeft += rg.GetEdges(v).size();
}
edgesLeft /= 2;
info["EdgesLeft"] = edgesLeft;
pair<Graph, Graph> lb_graphs = BuildLBGraph(rg); // for opt 23, 24
const auto& ga = lb_graphs.first;
const auto& gainv = lb_graphs.second;
void CplexEx1Solver::Transform(bool optimize) {
BuildReducedGraph(_graph);
IloEnv env;
IloModel model(env);
auto lb_graphs = BuildLBGraph(_graph); // for opt 23, 24
const auto& graph_a = lb_graphs.first;
const auto& graph_a_inv = lb_graphs.second;
auto yMapping = AddObjectiveAndVariablesY(env, model, rg);
AddAllOptimizations(env, model, rg, ga, gainv, yMapping);
auto y_mapping = AddObjectiveAndVariablesY(_env, _model, _graph);
AddAllOptimizations(_env, _model, _graph, graph_a, graph_a_inv, y_mapping);
auto xMapping = AddVariablesX(env, rg);
AddEx1Restrictions(env, model, rg, yMapping, xMapping, s);
auto x_mapping = AddVariablesX(_env, _graph);
AddEx1Restrictions(_env, _model, _graph, y_mapping, x_mapping, _start_vertex);
}
IloCplex cplex(env);
cplex.setOut(env.getNullStream());
wval CplexEx1Solver::Solve() {
IloCplex cplex(_env);
cplex.setOut(_env.getNullStream());
cplex.extract(model);
cplex.extract(_model);
cplex.solve();
wval result = cplex.getObjValue();
env.end();
_env.end();
return result;
}
......@@ -357,48 +342,37 @@ void AddNewEx2Restriction(
}
}
wval CplexEx2Connect(const Graph& g, unordered_map<string, size_t>& info) {
if (g.Size() <= 1) {
return 0;
}
Graph rg = BuildReducedGraph(g);
size_t edgesLeft = 0;
for (size_t v = 0; v < rg.Size(); ++v) {
edgesLeft += rg.GetEdges(v).size();
}
edgesLeft /= 2;
info["EdgesLeft"] = edgesLeft;
pair<Graph, Graph> lb_graphs = BuildLBGraph(rg); // for opt 23, 24
const auto& ga = lb_graphs.first;
const auto& gainv = lb_graphs.second;
void CplexEx2Solver::Transform(bool optimize) {
BuildReducedGraph(_graph);
IloEnv env;
IloModel model(env);
auto lb_graphs = BuildLBGraph(_graph); // for opt 23, 24
const auto& graph_a = lb_graphs.first;
const auto& graph_a_inv = lb_graphs.second;
auto yMapping = AddObjectiveAndVariablesY(env, model, rg);
AddAllOptimizations(env, model, rg, ga, gainv, yMapping);
_y_mapping = AddObjectiveAndVariablesY(_env, _model, _graph);
AddAllOptimizations(_env, _model, _graph, graph_a, graph_a_inv, _y_mapping);
auto zMapping = AddVariablesZ(env, rg);
AddEx2BaseRestrictions(env, model, rg, yMapping, zMapping);
_z_mapping = AddVariablesZ(_env, _graph);
AddEx2BaseRestrictions(_env, _model, _graph, _y_mapping, _z_mapping);
}
IloCplex cplex(env);
cplex.setOut(env.getNullStream());
wval CplexEx2Solver::Solve() {
IloCplex cplex(_env);
cplex.setOut(_env.getNullStream());
for (size_t iterations = 0;; ++iterations) {
cplex.extract(model);
cplex.extract(_model);
cplex.solve();
vector<vector<size_t>> components = FindIsolatedComponents(rg, cplex, yMapping);
vector<vector<size_t>> components = FindIsolatedComponents(_graph, cplex, _y_mapping);
if (components.size() == 1) {
cout << iterations << " restrictions were added" << endl;
break;
}
AddNewEx2Restriction(env, model, rg, components, zMapping);
AddNewEx2Restriction(_env, _model, _graph, components, _z_mapping);
}
wval result = cplex.getObjValue();
env.end();
_env.end();
return result;
}
#ifndef CPLEX_CONNECT_H
#define CPLEX_CONNECT_H
#include "solver.h"
#include "graph.h"
#include "numeric.h"
#include <string>
#include <unordered_map>
#define IL_STD
#include <ilcplex/ilocplex.h>
void ToCplexEx1(const Graph& g, const std::string output_file, size_t s, bool optimize);
wval CplexEx1Connect(const Graph& g, size_t s, std::unordered_map<std::string, size_t>& info);
class CplexSolver : public Solver {
public:
CplexSolver()
: _env()
, _model(_env)
{}
protected:
IloEnv _env;
IloModel _model;
};
class CplexEx1Solver : public CplexSolver {
public:
CplexEx1Solver(size_t start_vertex)
: _start_vertex(start_vertex)
{}
void Transform(bool optimize);
wval Solve();
private:
size_t _start_vertex;
};
class CplexEx2Solver : public CplexSolver {
public:
CplexEx2Solver() {}
wval CplexEx2Connect(const Graph& g, std::unordered_map<std::string, size_t>& info);
void Transform(bool optimize);
wval Solve();
private:
std::unordered_map<size_t, IloNumVar> _y_mapping;
std::unordered_map<size_t, IloNumVar> _z_mapping;
};
#endif
......@@ -4,6 +4,7 @@
#include "bounds.h"
#include <iostream>
#include <thread>
using namespace std;
......@@ -23,6 +24,43 @@ struct DpState {
: color_set(color_set), v(v), q_ind(q_ind) {}
};
struct MemoryForDp {
wval* Dval;
size_t* suffD;
size_t* binsearch_precalc;
};
MemoryForDp AllocateMemoryForDp(size_t max_C, size_t max_pg_size) {
size_t base1 = max_pg_size;
size_t base2 = base1 * max_pg_size;
size_t base3 = base2 * (1LL << max_C);
MemoryForDp memory;
memory.Dval = (wval*)malloc(sizeof(wval) * base3);
if (!memory.Dval) {
throw(std::bad_alloc());
}
memory.suffD = (size_t*)malloc(sizeof(size_t) * base3);
if (!memory.suffD) {
throw(std::bad_alloc());
}
memory.binsearch_precalc = (size_t*)malloc(sizeof(size_t) * max_pg_size * max_pg_size);
if (!memory.binsearch_precalc) {
throw(std::bad_alloc());
}
return memory;
}
void DeallocateMemoryForDp(MemoryForDp& memory) {
free(memory.Dval);
free(memory.suffD);
free(memory.binsearch_precalc);
}
bool CalcDp(
const Graph& pg,
const wval* lbounds,
......@@ -32,7 +70,9 @@ bool CalcDp(
const vector<ColorSet> colset_by_comp,
wval solution_ubound,
bool has_isolated_vertices,
wval& answer
wval& answer,
MemoryForDp& memory,
size_t threads_count
) {
const wval inf(solution_ubound + 1);
......@@ -41,18 +81,13 @@ bool CalcDp(
size_t base3 = base2 * (1LL << C);
vector<bool> Dok(base3);
wval* Dval = (wval*)malloc(sizeof(wval) * base3);
if (!Dval) throw(std::bad_alloc());
wval* Dval = memory.Dval;
size_t* suffD = memory.suffD;
size_t* binsearch_precalc = memory.binsearch_precalc;
DpState * bestdp_for_colset = (DpState*)malloc(sizeof(DpState) * (1LL << C));
DpState* bestdp_for_colset = (DpState*)malloc(sizeof(DpState) * (1LL << C));
if (!bestdp_for_colset) throw(std::bad_alloc());
size_t * suffD = (size_t*)malloc(sizeof(size_t) * base3);
if (!suffD) throw(std::bad_alloc());
size_t* binsearch_precalc = (size_t*)malloc(sizeof(size_t) * pg.Size() * pg.Size());
if (!binsearch_precalc) throw(std::bad_alloc());
for (size_t v = 0; v < pg.Size(); ++v) {
const vector<Edge> & Nv = pg.GetEdges(v);
for (size_t u_ind = 0; u_ind < Nv.size(); ++u_ind) {
......@@ -84,7 +119,15 @@ bool CalcDp(
ColorSet y = color_set + x;
color_set = (((color_set & ~y) / x) >> 1) | y;
}
for (size_t color_set_id = 0; color_set_id < color_sets.size(); ++color_set_id) {
vector<std::thread> threads;
for (size_t thread_id = 0; thread_id < threads_count; ++thread_id) {
// begin threading
threads.emplace_back([&](size_t initial_color_set_id) {
for (
size_t color_set_id = initial_color_set_id;
color_set_id < color_sets.size();
color_set_id += threads_count
) {
size_t color_set = color_sets[color_set_id];
bestdp_for_colset[color_set] = DpState(color_set, 0, 0);
......@@ -196,6 +239,12 @@ bool CalcDp(
}
}
}
}, thread_id);
// end threading
}
for (auto& thread : threads) {
thread.join();
}
}
DpState& best_dps = bestdp_for_colset[(ColorSet(1) << C) - 1];
......@@ -208,9 +257,6 @@ bool CalcDp(
bool good = Dok[base2 * best_dps.color_set + base1 * best_dps.v + best_dps.q_ind];
free(bestdp_for_colset);
free(suffD);
free(Dval);
free(binsearch_precalc);
return good;
}
......@@ -224,14 +270,28 @@ bool RunDp(
const vector<ColorSet> & colset_by_comp,
wval solution_ubound,
bool has_isolated_vertices,
wval& solution
wval& solution,
MemoryForDp& memory,
size_t threads_count
) {
if (C - 1 >= MAX_SET_SIZE)
throw AlgorithmException(string() + "Cannot create color set with " + std::to_string(C - 1) + "indice");
// calc DP
return CalcDp(pg, &lbounds[0], &col[0], C, comp_by_col, colset_by_comp, solution_ubound, has_isolated_vertices, solution);
return CalcDp(
pg,
&lbounds[0],
&col[0],
C,
comp_by_col,
colset_by_comp,
solution_ubound,
has_isolated_vertices,
solution,
memory,
threads_count
);
}
bool NextSubset(vector<size_t> & subset, size_t set_size) {
......@@ -260,137 +320,128 @@ size_t CalcAlgoRunTimes(const vector<size_t> & ci, double eps) {
}
wval DpConnect(
const Graph& g,
const vector<wval>& precalced_lbounds,
double eps,
default_random_engine& engine,
bool preprocess,
unordered_map<string, size_t>& info
) {
vector<wval> lbounds = precalced_lbounds;
size_t n = g.Size();
if (n <= 1)
return 0;
void DpSolver::Transform(bool optimize) {
_optimize = optimize;
_lowerbounds = CalcLB(_graph);
// find obligatory subgraph
Graph og = GetObligatoryGraph(g, lbounds);
Graph obligatory_graph = GetObligatoryGraph(_graph, _lowerbounds);
// find all components
vector<vector<size_t>> comps = og.FindComponents();
vector<size_t> comp_nums(n);
for (size_t comp_num = 0; comp_num < comps.size(); ++comp_num) {
for (size_t v : comps[comp_num])
_comps = obligatory_graph.FindComponents();
vector<size_t> comp_nums(_graph.Size());
for (size_t comp_num = 0; comp_num < _comps.size(); ++comp_num) {
for (size_t v : _comps[comp_num])
comp_nums[v] = comp_num;
}
// calc solution cost lower bound
wval solution_lbound = 0;
for (wval val : lbounds)
for (wval val : _lowerbounds)
solution_lbound += val;
vector<wval> ubounds;
wval solution_ubound;
if (preprocess) {
ubounds = CalcSolutionUBoundPrima(g);
solution_ubound = 0;
if (optimize) {
ubounds = CalcSolutionUBoundPrima(_graph);
_solution_upperbound = 0;
for (wval val : ubounds)
solution_ubound += val;
_solution_upperbound += val;
}
// build padded graph and sort edges
Graph pg(g.Size());
for (size_t v = 0; v < n; ++v) {
for (const Edge & edge : g.GetEdges(v)) {
Graph padded_graph(_graph.Size());
for (size_t v = 0; v < _graph.Size(); ++v) {
for (const Edge & edge : _graph.GetEdges(v)) {
size_t u = edge.to;
if (v < u && comp_nums[v] != comp_nums[u]) {
if (!preprocess ||
GreaterEq(solution_ubound, solution_lbound + (edge.weight - lbounds[v]) + (edge.weight - lbounds[u]))
if (!optimize ||
GreaterEq(
_solution_upperbound,
solution_lbound +
(edge.weight - _lowerbounds[v]) +
(edge.weight - _lowerbounds[u])
)
) {
pg.AddEdge(v, u, edge.weight);
padded_graph.AddEdge(v, u, edge.weight);
}
}
}
pg.AddEdge(v, v, wval(0));
padded_graph.AddEdge(v, v, wval(0));
}
// mapping
wval lost_answer = 0;
if (preprocess) {
vector<size_t> mapping(n, n);
size_t n_mapped = 0;
for (size_t v = 0; v < n; ++v) {
if (pg.Deg(v) > 1) {
mapping[v] = n_mapped++;
}
}
info["VerticesLeft"] = n_mapped;
Graph pg_mapped(n_mapped);
size_t edgesLeft = 0;
for (size_t v = 0; v < n; ++v) {
if (mapping[v] != n) {
for (const Edge& edge : pg.GetEdges(v)) {
_lost_answer = 0;
if (optimize) {
vector<size_t> mapping(padded_graph.Size(), padded_graph.Size());
size_t mapped_graph_size = 0;
for (size_t v = 0; v < padded_graph.Size(); ++v) {
if (padded_graph.Deg(v) > 1) {
mapping[v] = mapped_graph_size++;
}
}
Graph padded_graph_mapped(mapped_graph_size);
for (size_t v = 0; v < padded_graph.Size(); ++v) {
if (mapping[v] != padded_graph.Size()) {
for (const Edge& edge : padded_graph.GetEdges(v)) {
size_t u = edge.to;
if (mapping[v] < mapping[u]) {
pg_mapped.AddEdge(mapping[v], mapping[u], edge.weight);
++edgesLeft;
padded_graph_mapped.AddEdge(mapping[v], mapping[u], edge.weight);
}
}
pg_mapped.AddEdge(mapping[v], mapping[v], wval(0));
padded_graph_mapped.AddEdge(mapping[v], mapping[v], wval(0));
}
}
info["EdgesLeft"] = edgesLeft;
vector<size_t> comp_nums_mapped(n_mapped);
vector<wval> lbounds_mapped(n_mapped);
for (size_t v = 0; v < n; ++v) {
if (mapping[v] != n) {
comp_nums_mapped[mapping[v]] = comp_nums[v];
lbounds_mapped[mapping[v]] = lbounds[v];
vector<wval> lowerbounds_mapped(padded_graph_mapped.Size());
for (size_t v = 0; v < padded_graph.Size(); ++v) {
if (mapping[v] != padded_graph.Size()) {
lowerbounds_mapped[mapping[v]] = _lowerbounds[v];
} else {
lost_answer += lbounds[v];
_lost_answer += _lowerbounds[v];
}
}
vector<vector<size_t>> comps_mapped(comps.size());
vector<vector<size_t>> comps_mapped(_comps.size());
for (size_t comp_num = 0; comp_num < comps_mapped.size(); ++comp_num) {
for (size_t v : comps[comp_num]) {
if (mapping[v] != n) {
for (size_t v : _comps[comp_num]) {
if (mapping[v] != padded_graph.Size()) {
comps_mapped[comp_num].push_back(mapping[v]);
}
}
}
pg = pg_mapped;
n = n_mapped;
swap(comps, comps_mapped);
swap(lbounds, lbounds_mapped);
swap(comp_nums, comp_nums_mapped);
swap(_comps, comps_mapped);
swap(_lowerbounds, lowerbounds_mapped);
_graph = padded_graph_mapped;
} else {
_graph = padded_graph;
}
for (size_t v = 0; v < n; ++v)
pg.SortEdges(v);
for (size_t v = 0; v < _graph.Size(); ++v)
_graph.SortEdges(v);
if (!preprocess) {
ubounds = CalcSolutionUBound(pg, lbounds);
solution_ubound = 0;
if (!optimize) {
ubounds = CalcSolutionUBound(_graph, _lowerbounds);
_solution_upperbound = 0;
for (wval val : ubounds)
solution_ubound += val;
_solution_upperbound += val;
}
}
const wval inf(solution_ubound + 1);
wval DpSolver::Solve() {
const wval inf(_solution_upperbound + 1);
// run dp for random colorings
size_t c = comps.size();
vector<size_t> subset(c); // this stands for (c1, c1 + c2, ... , c1 + c2 + ... + cc), which is a c-size subset of (1, ... , 2 * (c - 1))
size_t c = _comps.size();
vector<size_t> subset(c); // this stands for (c1, c1 + c2, ... , c1 + c2 + ... + cc),
// which is a c-size subset of (1, ... , 2 * (c - 1))
for (size_t i = 0; i < c; ++i)
subset[i] = i + 1;
bool found_solution = false;
wval best_solution;
best_solution = solution_ubound + 1;
best_solution = _solution_upperbound + 1;
MemoryForDp memory = AllocateMemoryForDp(2 * c - 2, _graph.Size());
do {
// restore ci from q
......@@ -398,22 +449,23 @@ wval DpConnect(
bool sizes_correct = true;
for (size_t i = 0; i < c; ++i) {
ci[i] = subset[i] - ((i > 0) ? subset[i - 1] : 0);
sizes_correct = sizes_correct && (ci[i] <= comps[i].size());
sizes_correct = sizes_correct && (ci[i] <= _comps[i].size());
}
if (!sizes_correct) {
continue;
}
size_t times = CalcAlgoRunTimes(ci, eps);
if (times >= INF)
size_t times = CalcAlgoRunTimes(ci, _eps);
if (times >= INF) {
throw new AlgorithmException("DP should be run too many times for this eps.");
}
// repeat dp on random colorings
while (times--) {
// color vertices
vector<ColorSet> colset_by_comp(c);
vector<size_t> comp_by_col(subset.back());
vector<size_t> col(n);
vector<size_t> col(_graph.Size());
for (size_t i = 0; i < c; ++i) {
size_t l = (i > 0) ? subset[i - 1] : 0;
size_t r = subset[i] - 1;
......@@ -423,19 +475,34 @@ wval DpConnect(
for (size_t col_id = l; col_id <= r; ++col_id)
comp_by_col[col_id] = i;
const vector<size_t> & comp = comps[i];
const vector<size_t> & comp = _comps[i];
for (size_t v : comp)
col[v] = distr(engine);
col[v] = distr(_random_engine);
}
// run DP in coloring and relax the answer
wval curr_solution;
if (RunDp(pg, lbounds, col, subset.back(), comp_by_col, colset_by_comp, solution_ubound, preprocess, curr_solution)) {
if (RunDp(
_graph,
_lowerbounds,
col,
subset.back(),
comp_by_col,
colset_by_comp,
_solution_upperbound,
_optimize,
curr_solution,
memory,
_threads_num
)
) {
RelaxMin(best_solution, curr_solution);
}
}
} while (NextSubset(subset, max(c, (c - 1) * 2)));
best_solution += lost_answer;
DeallocateMemoryForDp(memory);
best_solution += _lost_answer;
return best_solution;
}
#ifndef DP_CONNECT_H
#define DP_CONNECT_H
#include "solver.h"
#include "numeric.h"
#include "graph.h"
......@@ -8,13 +9,30 @@
#include <random>
#include <unordered_map>
wval DpConnect(
const Graph& g,
const std::vector<wval>& lbounds,
double eps,
std::default_random_engine& engine,
bool preprocess,
std::unordered_map<std::string, size_t>& info
);
class DpSolver : public Solver {
public:
DpSolver(double eps, size_t seed, size_t threads_num)
: _eps(eps)
, _random_engine(seed)
, _threads_num(threads_num)
, _transformed(false)
, _optimize(false)
{}
void Transform(bool optimize);
wval Solve();
private:
double _eps;
std::default_random_engine _random_engine;
size_t _threads_num;
bool _optimize;
bool _transformed;
std::vector<wval> _lowerbounds;
std::vector<std::vector<size_t>> _comps;
wval _solution_upperbound;
wval _lost_answer;
};
#endif
......@@ -18,6 +18,8 @@ struct Edge {
class Graph {
public:
Graph() : Graph(0) {}
Graph(size_t size);