qbpp_easy_solver.hpp

Go to the documentation of this file.

 
#ifndef QBPP_EASY_SOLVER_HPP
#define QBPP_EASY_SOLVER_HPP
 
#include <tbb/blocked_range.h>
#include <tbb/parallel_for.h>
 
#include <boost/circular_buffer.hpp>
#include <limits>
#include <random>
#include <set>
#include <thread>
#include <vector>
 
#include "qbpp.hpp"
#include "qbpp_misc.hpp"
 
namespace qbpp {
 
namespace easy_solver {
 
 
class RandomSet;
class MinSet;
class RandomMinSet;
class Tabu;
class SolDelta;
class TabuSolDelta;
class PosMinSolDelta;
class EasySolver;
 
 
class RandomSet {
  std::vector<vindex_t> position_;
 
  std::vector<vindex_t> variables_;
 
 public:
  explicit RandomSet(vindex_t size) : position_(size, vindex_limit) {
    variables_.reserve(size);
  }
 
  vindex_t var_count() const {
    return static_cast<vindex_t>(variables_.size());
  }
 
  void insert(vindex_t index) {
    if (position_[index] != vindex_limit) {
      throw std::runtime_error(THROW_MESSAGE("RandomSet: Insert variable (",
                                             index, ") already in set"));
    }
    variables_.push_back(index);
    position_[index] = static_cast<vindex_t>(variables_.size() - 1);
  }
 
  void swap(vindex_t pos1, vindex_t pos2) {
    if (pos1 == pos2) return;
    std::swap(variables_[pos1], variables_[pos2]);
    position_[variables_[pos1]] = pos1;
    position_[variables_[pos2]] = pos2;
  }
 
  void erase(vindex_t index) {
    if (position_[index] == vindex_limit) {
      throw std::runtime_error(
          THROW_MESSAGE("RandomSet: Erase variable (", index, ") not in set"));
    }
    swap(position_[index], var_count() - 1);
    variables_.pop_back();
    position_[index] = vindex_limit;
  }
 
  bool has(vindex_t index) const { return position_[index] != vindex_limit; }
 
  vindex_t select_at_random() const {
    return variables_[qbpp::misc::RandomGenerator::gen(var_count())];
  }
 
  void print(const std::string &prefix = "") const {
    std::cout << prefix << " Size = " << var_count() << " : ";
    for (vindex_t i = 0; i < var_count(); ++i) {
      std::cout << " " << variables_[i];
    }
    std::cout << std::endl;
  }
};
 
class MinSet {
  using ValIndexMap = std::set<std::pair<energy_t, vindex_t>>;
 
  ValIndexMap set_;
 
 public:
  MinSet() = default;
 
  vindex_t var_count() const { return static_cast<vindex_t>(set_.size()); }
 
  void insert(vindex_t index, energy_t delta) {
    auto tuple_val = std::make_pair(delta, index);
    auto result = set_.insert(tuple_val);
    if (!result.second) {
      throw std::runtime_error(THROW_MESSAGE("MinSet: Insert variable (", index,
                                             ") already in set"));
    }
  }
 
  void erase(vindex_t index, energy_t delta) {
    auto pair_val = std::make_pair(delta, index);
    auto result = set_.erase(pair_val);
    if (result == 0) {
      throw std::runtime_error(
          THROW_MESSAGE("MinSet: Erase variable (", index, ") not in set"));
    }
  }
 
  vindex_t get_first() const { return (*set_.begin()).second; }
 
  void print(const std::string &prefix) const {
    std::cout << prefix;
    for (auto &[val, i] : set_) {
      std::cout << "(" << i << "," << val << ")";
    }
    std::cout << std::endl;
  }
};
 
class RandomMinSet {
  MinSet min_set_;
  RandomSet random_set_;
 
 public:
  RandomMinSet(vindex_t size) : random_set_(size) {};
 
  vindex_t var_count() const {
    return static_cast<vindex_t>(min_set_.var_count());
  }
 
  void insert(vindex_t index, energy_t delta) {
    min_set_.insert(index, delta);
    random_set_.insert(index);
  }
 
  void erase(vindex_t index, energy_t delta) {
    min_set_.erase(index, delta);
    random_set_.erase(index);
  }
 
  bool has(vindex_t index) const { return random_set_.has(index); }
 
  vindex_t select_at_random() const { return random_set_.select_at_random(); }
 
  vindex_t get_first() const { return min_set_.get_first(); }
 
  void print(const std::string &prefix) const {
    min_set_.print(prefix + " MIN:");
    random_set_.print(prefix + "RANDOM:");
  }
};
 
class Tabu {
  impl::BitVector variables_;
 
  boost::circular_buffer<vindex_t> tabu_list;
 
 public:
  vindex_t tabu_size() const { return static_cast<vindex_t>(tabu_list.size()); }
 
  Tabu(vindex_t size, vindex_t tabu_size)
      : variables_(size), tabu_list(tabu_size) {};
 
  bool has(vindex_t index) const { return variables_.get(index); }
 
  void insert(vindex_t index) {
    if (has(index)) {
      auto it = std::find(tabu_list.begin(), tabu_list.end(), index);
      tabu_list.erase(it);
    } else if (tabu_list.full()) {
      variables_.set(tabu_list.front(), false);
      tabu_list.pop_front();
    }
    tabu_list.push_back(index);
    variables_.set(index, true);
  }
 
  vindex_t non_tabu_random() {
    vindex_t index;
    do {
      index = qbpp::misc::RandomGenerator::gen(variables_.size());
    } while (has(index));
    return index;
  }
 
  void print() {
    std::cout << "Tabu list:";
    for (const auto &i : tabu_list) {
      std::cout << " " << i;
    }
    std::cout << std::endl;
  }
};
 
class SolDelta : public Sol {
 protected:
  const EasySolver &easy_solver_;
 
  std::vector<energy_t> delta_;
 
  RandomMinSet neg_set_;
 
  std::shared_ptr<qbpp::SolHolder> sol_holder_ptr_;
 
 public:
  explicit SolDelta(const EasySolver &easy_solver);
 
  virtual ~SolDelta() = default;
 
  energy_t get_delta(vindex_t i) const { return delta_[i]; }
 
  virtual void flip(vindex_t index);
 
  std::optional<double> set_if_better(const std::string &solver_name);
 
  vindex_t neg_count() const { return neg_set_.var_count(); }
 
  vindex_t neg_random() { return neg_set_.select_at_random(); }
 
  vindex_t neg_min() const { return neg_set_.get_first(); }
 
  virtual void before_delta_updated(
      [[maybe_unused]] vindex_t flip_index,
      [[maybe_unused]] vindex_t update_delta_index) {}
 
  virtual void after_delta_updated(
      [[maybe_unused]] vindex_t flip_index,
      [[maybe_unused]] vindex_t update_delta_index) {}
 
  void greedy();
 
  void random_flip(size_t iteration);
 
  void move_to(qbpp::Sol destination);
};
 
class TabuSolDelta : public SolDelta {
 protected:
  const vindex_t tabu_size_;
 
  Tabu tabu_;
 
 public:
  TabuSolDelta(const EasySolver &easy_solver, vindex_t tabu_size)
      : SolDelta(easy_solver),
        tabu_size_(tabu_size),
        tabu_(var_count(), tabu_size) {}
 
  virtual ~TabuSolDelta() = default;
 
  void flip(vindex_t index) override {
    SolDelta::flip(index);
    tabu_.insert(index);
  }
 
  bool tabu_has(vindex_t index) const { return tabu_.has(index); }
 
  vindex_t non_tabu_random() { return tabu_.non_tabu_random(); }
};
 
class PosMinSolDelta : public TabuSolDelta {
  MinSet pos_set_;
 
 public:
  explicit PosMinSolDelta(const EasySolver &easy_solver);
 
  virtual ~PosMinSolDelta() = default;
 
  vindex_t pos_count() const { return pos_set_.var_count(); }
 
  vindex_t pos_min() const { return pos_set_.get_first(); }
 
  void before_delta_updated(vindex_t , vindex_t k) override {
    if (delta_[k] > 0) {
      pos_set_.erase(k, delta_[k]);
    }
  }
 
  void after_delta_updated(vindex_t , vindex_t k) override {
    if (delta_[k] > 0) {
      pos_set_.insert(k, delta_[k]);
    }
  }
 
  void search(size_t iteration);
 
  void print() const {
    for (vindex_t i = 0; i < var_count(); ++i) {
      std::cout << "(" << i << "," << delta_[i] << ")";
    }
    std::cout << std::endl;
    neg_set_.print("NEG:");
    pos_set_.print("POS:");
  }
};
 
class EasySolver {
 protected:
  const qbpp::QuadModel quad_model_;
 
  std::optional<uint32_t> time_limit_;
 
  size_t thread_count_{std::thread::hardware_concurrency() < 4
                           ? 4
                           : std::thread::hardware_concurrency()};
 
  std::optional<energy_t> target_energy_;
 
  const std::shared_ptr<qbpp::SolHolder> sol_holder_ptr_;
 
  mutable std::mutex callback_mutex_;
 
  bool enable_default_callback_ = false;
 
  void single_search(size_t thread_id);
 
  double start_time_;
 
  mutable std::optional<energy_t> prev_energy_ = std::nullopt;
 
 public:
 
  explicit EasySolver(const qbpp::QuadModel &quad_model,
                      std::shared_ptr<qbpp::SolHolder> sol_holder_ptr = nullptr)
      : quad_model_(quad_model),
        sol_holder_ptr_(sol_holder_ptr
                            ? sol_holder_ptr
                            : std::make_shared<qbpp::SolHolder>(quad_model_)) {}
 
 
  virtual ~EasySolver() = default;
 
  const qbpp::QuadModel &get_quad_model() const { return quad_model_; }
 
  const std::optional<uint32_t> get_time_limit() const { return time_limit_; }
 
  const std::optional<energy_t> get_target_energy() const {
    return target_energy_;
  }
 
  const std::shared_ptr<qbpp::SolHolder> get_sol_holder_ptr() const {
    return sol_holder_ptr_;
  }
 
  vindex_t var_count() const { return quad_model_.var_count(); }
 
  void set_time_limit(uint32_t limit) { time_limit_ = limit; }
 
  void set_target_energy(energy_t energy) { target_energy_ = energy; }
 
  void set_thread_count(unsigned int count) { thread_count_ = count; }
 
  size_t get_thread_count() { return thread_count_; }
 
  qbpp::Sol get_sol() const { return sol_holder_ptr_->get_sol(); }
 
  double get_tts() const { return sol_holder_ptr_->get_tts(); }
 
  void enable_default_callback() { enable_default_callback_ = true; }
 
  virtual void callback(const qbpp::Sol &sol, double tts) const {
    std::lock_guard<std::mutex> lock(callback_mutex_);
    if (enable_default_callback_) {
      if (!prev_energy_.has_value() ||
          sol.get_energy() < prev_energy_.value()) {
        std::cout << "TTS = " << std::fixed << std::setprecision(3)
                  << std::setfill('0') << tts
                  << "s Energy = " << sol.get_energy() << std::endl;
        prev_energy_ = sol.get_energy();
      }
    }
  }
 
  qbpp::Sol search();
 
};
 
 
inline SolDelta::SolDelta(const EasySolver &easy_solver)
    : Sol(easy_solver.get_quad_model()),
      easy_solver_(easy_solver),
      delta_(var_count()),
      neg_set_(var_count()),
      sol_holder_ptr_(easy_solver.get_sol_holder_ptr()) {
  for (vindex_t i = 0; i < var_count(); ++i) {
    delta_[i] = quad_model_.get_linear(i);
  }
}
 
inline void SolDelta::flip(vindex_t index) {
  if (index >= var_count()) {
    throw std::out_of_range(
        THROW_MESSAGE("Flip index (", index, ") is out of range."));
  }
 
  for (vindex_t j = 0; j < quad_model_.get_degree(index); ++j) {
    auto [k, coeff] = quad_model_.get_quadratic(index, j);
 
    if (delta_[k] < 0) {
      neg_set_.erase(k, delta_[k]);
    }
    before_delta_updated(index, k);
 
    delta_[k] +=
        static_cast<energy_t>((2 * get(index) - 1) * (2 * get(k) - 1) * coeff);
 
    if (delta_[k] < 0) {
      neg_set_.insert(k, delta_[k]);
    }
    after_delta_updated(index, k);
  }
 
  if (delta_[index] < 0) {
    neg_set_.erase(index, delta_[index]);
  }
  before_delta_updated(index, index);
 
  delta_[index] = -delta_[index];
 
  if (delta_[index] < 0) {
    neg_set_.insert(index, delta_[index]);
  }
 
  after_delta_updated(index, index);
 
  if (energy_.has_value()) {
    *energy_ -= delta_[index];
  } else {
    energy_ = comp_energy();
  }
 
  bit_vector_.flip(index);
}
 
inline void SolDelta::greedy() {
  while (neg_count() > 0) {
    vindex_t min_index = neg_min();
    flip(min_index);
  }
  set_if_better("Easy");
}
 
inline void SolDelta::random_flip(size_t iteration) {
  for (vindex_t i = 0; i < iteration; ++i) {
    vindex_t flip_index = qbpp::misc::RandomGenerator::gen(this->var_count());
    flip(flip_index);
    set_if_better("Easy");
  }
}
 
inline void SolDelta::move_to(qbpp::Sol destination) {
  MinSet to_be_flipped;
  for (vindex_t i = 0; i < var_count(); ++i) {
    if (get(i) != destination.get(i)) {
      to_be_flipped.insert(i, delta_[i]);
    }
  }
  while (to_be_flipped.var_count() > 0) {
    vindex_t min_index = to_be_flipped.get_first();
    for (vindex_t j = 0; j < quad_model_.get_degree(min_index); ++j) {
      auto pair = quad_model_.get_quadratic(min_index, j);
      if (get(pair.first) != destination.get(pair.first)) {
        to_be_flipped.erase(pair.first, delta_[pair.first]);
      }
    }
    to_be_flipped.erase(min_index, delta_[min_index]);
    flip(min_index);
    for (vindex_t j = 0; j < quad_model_.get_degree(min_index); ++j) {
      auto pair = quad_model_.get_quadratic(min_index, j);
      if (get(pair.first) != destination.get(pair.first))
        to_be_flipped.insert(pair.first, delta_[pair.first]);
    }
    set_if_better("Easy");
  }
}
 
inline std::optional<double> SolDelta::set_if_better(
    const std::string &solver_name) {
  std::optional<double> tts =
      sol_holder_ptr_->set_if_better(*this, solver_name);
  if (tts.has_value()) {
    easy_solver_.callback(*this, tts.value());
  }
  return tts;
}
 
 
inline PosMinSolDelta::PosMinSolDelta(const EasySolver &easy_solver)
    : TabuSolDelta(easy_solver,
                   std::min(easy_solver.get_quad_model().var_count() / 2,
                            static_cast<vindex_t>(5))) {
  for (vindex_t i = 0; i < var_count(); ++i) {
    if (delta_[i] < 0) {
      neg_set_.insert(i, delta_[i]);
    } else if (delta_[i] > 0) {
      pos_set_.insert(i, delta_[i]);
    }
  }
}
 
inline void PosMinSolDelta::search(size_t iteration) {
  greedy();
  for (vindex_t i = 0; i < iteration; ++i) {
    std::optional<vindex_t> flip_index = std::nullopt;
    vindex_t min_index = neg_min();
    if (neg_count() >= 1 &&
        get_energy() + get_delta(min_index) < sol_holder_ptr_->get_energy()) {
      flip_index = min_index;
    } else {
      for (vindex_t j = 0; j < tabu_size_ * 2; ++j) {
        vindex_t candidate = vindex_limit;
        if (pos_count() == 0 || neg_count() == 0) {
          candidate = qbpp::misc::RandomGenerator::gen(var_count());
        } else if (qbpp::misc::RandomGenerator::gen(neg_count() + 1) == 0) {
          candidate = pos_min();
        } else {
          candidate = neg_random();
        }
        if (candidate == vindex_limit) {
          throw std::runtime_error(THROW_MESSAGE("Unexpected error."));
        }
        if (!tabu_.has(candidate)) {
          flip_index = candidate;
          break;
        }
      }
      if (!flip_index.has_value()) {
        flip_index = tabu_.non_tabu_random();
      }
    }
    flip(flip_index.value());
    tabu_.insert(flip_index.value());
    set_if_better("Easy");
    if (easy_solver_.get_target_energy().has_value() &&
        sol_holder_ptr_->get_energy() <= easy_solver_.get_target_energy())
      break;
  }
  greedy();
}
 
inline void EasySolver::single_search(size_t thread_id) {
  const size_t min_flip_count = 100;
  PosMinSolDelta pos_min_sol_delta(*this);
  size_t max_flip_count = (var_count() / 2);
  if (thread_count_ > 1) {
    double k = std::pow(static_cast<double>(max_flip_count) / min_flip_count,
                        1.0 / static_cast<double>(thread_count_ - 1));
    max_flip_count =
        static_cast<size_t>(min_flip_count * std::pow(k, thread_id));
    if (max_flip_count < min_flip_count) max_flip_count = min_flip_count;
  }
  energy_t prev_energy = sol_holder_ptr_->get_energy();
  size_t random_flip_count = 1;
  while (!time_limit_.has_value() ||
         qbpp::get_time() - start_time_ < time_limit_.value()) {
    pos_min_sol_delta.move_to(sol_holder_ptr_->get_sol());
    pos_min_sol_delta.random_flip(random_flip_count);
    pos_min_sol_delta.search(max_flip_count);
    if (target_energy_.has_value() &&
        sol_holder_ptr_->get_energy() <= target_energy_.value())
      break;
    if (prev_energy == sol_holder_ptr_->get_energy()) {
      random_flip_count = (random_flip_count + 1) % (max_flip_count / 2);
      if (random_flip_count == 0) random_flip_count = 1;
    } else {
      random_flip_count = 1;
    }
    prev_energy = sol_holder_ptr_->get_energy();
  }
}
 
inline qbpp::Sol EasySolver::search() {
  if (quad_model_.term_count() == 0) return sol_holder_ptr_->get_sol();
  start_time_ = qbpp::get_time();
  tbb::parallel_for(size_t(0), static_cast<size_t>(thread_count_),
                    [this](size_t thread_id) { single_search(thread_id); });
  return sol_holder_ptr_->get_sol();
}
 
}  
}  
 
#endif