Commit 770aa6c4 authored by Gauthier Quesnel's avatar Gauthier Quesnel
Browse files

initializer: reuse same algorithms for both solver/optimizer

parent 38bb194f
Pipeline #5758 passed with stage
in 3 minutes and 12 seconds
......@@ -255,6 +255,122 @@ is_signbit_change(Float lhs, Float rhs) noexcept
return std::signbit(lhs) != std::signbit(rhs);
}
template<typename Cost, typename Mode>
void
init_with_bastert(bit_array& x,
const Cost& c,
const int variables,
const int value_if_0) noexcept
{
for (int i = 0; i != variables; ++i)
if (init_x<Mode>(c[i], value_if_0))
x.set(i);
else
x.unset(i);
}
inline void
init_with_random(bit_array& x,
random_engine& rng,
const int variables,
const double init_ramdom) noexcept
{
std::bernoulli_distribution dist(init_ramdom);
for (int i = 0; i != variables; ++i)
if (dist(rng))
x.set(i);
else
x.unset(i);
}
template<typename Cost, typename Mode>
void
init_with_pre_solve(bit_array& x_pessimistic,
bit_array& x_optimistic,
random_engine& rng,
const Cost& c,
const std::vector<merged_constraint>& constraints) noexcept
{
int max_length = 0;
for (const auto& cst : constraints)
if (length(cst.elements) > max_length)
max_length = length(cst.elements);
struct reduced_cost
{
float value;
int factor;
int id;
};
std::vector<reduced_cost> R(max_length);
for (const auto& cst : constraints) {
R.resize(cst.elements.size());
const int r_size = length(cst.elements);
for (int i = 0; i != r_size; ++i) {
R[i].value = static_cast<float>(c[cst.elements[i].variable_index]);
R[i].factor = cst.elements[i].factor;
R[i].id = cst.elements[i].variable_index;
}
std::shuffle(std::begin(R), std::end(R), rng);
std::sort(
std::begin(R), std::end(R), [](const auto& lhs, const auto& rhs) {
return is_better_solution<Mode>(lhs.value, rhs.value);
});
if (!x_pessimistic.empty()) {
int sum = 0;
int best = -2;
for (int i = -1; i < r_size; ++i) {
if (cst.min <= sum && sum <= cst.max) {
best = i;
break;
}
if (i + 1 < r_size)
sum += R[i + 1].factor;
}
int i = 0;
for (; i <= best; ++i)
x_pessimistic.set(R[i].id);
for (; i != r_size; ++i) {
x_pessimistic.unset(R[i].id);
}
}
if (!x_optimistic.empty()) {
int sum = 0;
int best = -2;
for (int i = -1; i < r_size; ++i) {
if (cst.min <= sum && sum <= cst.max)
best = i;
if (best != -2 && i + 1 < r_size &&
stop_iterating<Mode>(R[i + 1].value))
break;
if (i + 1 < r_size)
sum += R[i + 1].factor;
}
int i = 0;
for (; i <= best; ++i)
x_optimistic.set(R[i].id);
for (; i != r_size; ++i)
x_optimistic.unset(R[i].id);
}
}
}
/// The bkmin and bkmax constraint bounds are not equal and can be assigned to
/// -infinity or +infinity. We have to scan the r vector and search a value j
/// such as b(0, k) <= Sum A(k, R[j]) < b(1, k).
......
......@@ -36,6 +36,7 @@
#include <utility>
#include <vector>
#include "itm-common.hpp"
#include "private.hpp"
#include "result.hpp"
#include "sparse-matrix.hpp"
......@@ -44,122 +45,6 @@
namespace baryonyx {
namespace itm {
template<typename Cost, typename Mode>
void
init_with_bastert(bit_array& x,
const Cost& c,
const int variables,
const int value_if_0) noexcept
{
for (int i = 0; i != variables; ++i)
if (init_x<Mode>(c[i], value_if_0))
x.set(i);
else
x.unset(i);
}
inline void
init_with_random(bit_array& x,
random_engine& rng,
const int variables,
const double init_ramdom) noexcept
{
std::bernoulli_distribution dist(init_ramdom);
for (int i = 0; i != variables; ++i)
if (dist(rng))
x.set(i);
else
x.unset(i);
}
template<typename Cost, typename Mode>
void
init_with_pre_solve(bit_array& x_pessimistic,
bit_array& x_optimistic,
random_engine& rng,
const Cost& c,
const std::vector<merged_constraint>& constraints) noexcept
{
int max_length = 0;
for (const auto& cst : constraints)
if (length(cst.elements) > max_length)
max_length = length(cst.elements);
struct reduced_cost
{
float value;
int factor;
int id;
};
std::vector<reduced_cost> R(max_length);
for (const auto& cst : constraints) {
R.resize(cst.elements.size());
const int r_size = length(cst.elements);
for (int i = 0; i != r_size; ++i) {
R[i].value = static_cast<float>(c[cst.elements[i].variable_index]);
R[i].factor = cst.elements[i].factor;
R[i].id = cst.elements[i].variable_index;
}
std::shuffle(std::begin(R), std::end(R), rng);
std::sort(
std::begin(R), std::end(R), [](const auto& lhs, const auto& rhs) {
return is_better_solution<Mode>(lhs.value, rhs.value);
});
{
int sum = 0;
int best = -2;
for (int i = -1; i < r_size; ++i) {
if (cst.min <= sum && sum <= cst.max) {
best = i;
break;
}
if (i + 1 < r_size)
sum += R[i + 1].factor;
}
int i = 0;
for (; i <= best; ++i)
x_pessimistic.set(R[i].id);
for (; i != r_size; ++i) {
x_pessimistic.unset(R[i].id);
}
}
{
int sum = 0;
int best = -2;
for (int i = -1; i < r_size; ++i) {
if (cst.min <= sum && sum <= cst.max)
best = i;
if (best != -2 && i + 1 < r_size &&
stop_iterating<Mode>(R[i + 1].value))
break;
if (i + 1 < r_size)
sum += R[i + 1].factor;
}
int i = 0;
for (; i <= best; ++i)
x_optimistic.set(R[i].id);
for (; i != r_size; ++i)
x_optimistic.unset(R[i].id);
}
}
}
template<typename Cost, typename Mode>
struct storage
{
......@@ -616,12 +501,12 @@ struct optimize_functor
m_rng, length(constraints), variables, norm_costs, constraints);
compute_order compute(p.order, variables);
compute.init(slv, x);
bool is_a_solution = false;
while (!stop_task.load()) {
auto kappa = static_cast<Float>(best_recorder.reinit(
m_thread_id, is_a_solution, p.kappa_min, p.kappa_max, x));
compute.init(slv, x);
auto best_remaining = INT_MAX;
is_a_solution = false;
......
......@@ -31,6 +31,7 @@
#include <utility>
#include <vector>
#include "itm-common.hpp"
#include "observer.hpp"
#include "private.hpp"
#include "sparse-matrix.hpp"
......@@ -105,12 +106,35 @@ struct solver_functor
Solver slv(
m_rng, length(constraints), variables, norm_costs, constraints);
// // solver_initializer<Solver, Float, Mode> initializer(
// // slv, p.init_policy, p.init_policy_random, p.init_random);
// initializer.init(slv, x);
compute_order compute(p.order, variables);
compute.init(slv, x);
{
std::bernoulli_distribution choose_mutation(
m_ctx->parameters.init_random);
bit_array empty_x;
switch (p.init_policy) {
case solver_parameters::init_policy_type::pessimistic_solve:
init_with_pre_solve<Cost, Mode>(
x, empty_x, m_rng, original_costs, constraints);
break;
case solver_parameters::init_policy_type::optimistic_solve:
init_with_pre_solve<Cost, Mode>(
empty_x, x, m_rng, original_costs, constraints);
break;
case solver_parameters::init_policy_type::bastert:
case solver_parameters::init_policy_type::cycle:
case solver_parameters::init_policy_type::crossover_cycle:
init_with_bastert<Cost, Mode>(x, original_costs, variables, 0);
break;
}
for (int i = 0, e = x.size(); i != e; ++i)
if (choose_mutation(m_rng))
x.invert(i);
}
bool start_push = false;
auto kappa = static_cast<Float>(p.kappa_min);
......@@ -120,6 +144,8 @@ struct solver_functor
m_begin = std::chrono::steady_clock::now();
m_end = std::chrono::steady_clock::now();
compute.init(slv, x);
for (int i = 0; i != p.limit; ++i) {
auto remaining = compute.run(slv, x, m_rng, kappa, delta, theta);
obs.make_observation(slv.ap, slv.P.get(), slv.pi.get());
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment