adaptive_diagonal_curves.h

Go to the documentation of this file.

/*
 * Copyright 2021 MusicScience37 (Kenta Kabashima)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <iterator>
#include <limits>
#include <memory>
#include <queue>
#include <string_view>
#include <utility>
#include <vector>
 
#include <hash_tables/maps/multi_open_address_map_st.h>
 
#include "num_collect/base/concepts/real_scalar.h"
#include "num_collect/base/exception.h"
#include "num_collect/base/index_type.h"
#include "num_collect/base/precondition.h"
#include "num_collect/logging/iterations/iteration_logger.h"
#include "num_collect/logging/log_tag_view.h"
#include "num_collect/logging/logging_macros.h"
#include "num_collect/opt/concepts/multi_variate_objective_function.h"
#include "num_collect/opt/concepts/objective_function.h"
#include "num_collect/opt/impl/ternary_vector.h"
#include "num_collect/opt/optimizer_base.h"
#include "num_collect/util/assert.h"
#include "num_collect/util/safe_cast.h"
 
namespace num_collect::opt {
 
constexpr auto adaptive_diagonal_curves_tag =
    logging::log_tag_view("num_collect::opt::adaptive_diagonal_curves");
 
namespace impl {
 
template <concepts::objective_function ObjectiveFunction>
class adc_sample_dict;
 
template <concepts::multi_variate_objective_function ObjectiveFunction>
class adc_sample_dict<ObjectiveFunction> {
public:
    using objective_function_type = ObjectiveFunction;
 
    using variable_type = typename objective_function_type::variable_type;
 
    using value_type = typename objective_function_type::value_type;
 
    explicit adc_sample_dict(
        const objective_function_type& obj_fun = objective_function_type())
        : obj_fun_(obj_fun) {
        constexpr std::size_t initial_space = 10000;
        value_dict_.reserve_approx(initial_space);
    }
 
    void change_objective_function(const objective_function_type& obj_fun) {
        obj_fun_ = obj_fun;
    }
 
    void init(const variable_type& lower, const variable_type& upper) {
        NUM_COLLECT_PRECONDITION(lower.size() == upper.size(),
            "Element-wise limits must have the same size.");
        NUM_COLLECT_PRECONDITION((lower.array() < upper.array()).all(),
            "Element-wise limits must satisfy lower < upper for each element.");
 
        lower_ = lower;
        width_ = upper - lower;
        dim_ = lower.size();
        value_dict_.clear();
    }
 
    [[nodiscard]] auto operator()(const ternary_vector& point) -> value_type {
        return value_dict_.get_or_create_with_factory(
            point, [this, &point] { return evaluate_on(point); });
    }
 
    [[nodiscard]] auto dim() const -> index_type { return dim_; }
 
    [[nodiscard]] auto opt_variable() const -> const variable_type& {
        return opt_variable_;
    }
 
    [[nodiscard]] auto opt_point() const -> const ternary_vector& {
        return opt_point_;
    }
 
    [[nodiscard]] auto opt_value() const -> const value_type& {
        return opt_value_;
    }
 
    [[nodiscard]] auto evaluations() const noexcept -> index_type {
        return static_cast<index_type>(value_dict_.size());
    }
 
private:
    [[nodiscard]] auto evaluate_on(const ternary_vector& point) -> value_type {
        NUM_COLLECT_DEBUG_ASSERT(point.dim() == dim_);
        auto var = variable_type(dim_);
        for (index_type i = 0; i < dim_; ++i) {
            var(i) = lower_(i) +
                width_(i) * point.elem_as<typename variable_type::Scalar>(i);
        }
        obj_fun_.evaluate_on(var);
 
        if (value_dict_.empty() || obj_fun_.value() < opt_value_) {
            opt_point_ = point;
            opt_variable_ = var;
            opt_value_ = obj_fun_.value();
        }
 
        return obj_fun_.value();
    }
 
    objective_function_type obj_fun_;
 
    variable_type lower_{};
 
    variable_type width_{};
 
    index_type dim_{0};
 
    hash_tables::maps::multi_open_address_map_st<impl::ternary_vector,
        value_type>
        value_dict_{};
 
    ternary_vector opt_point_{};
 
    variable_type opt_variable_{};
 
    value_type opt_value_{};
};
 
template <base::concepts::real_scalar Value>
class adc_rectangle {
public:
    using value_type = Value;
 
    adc_rectangle(
        const impl::ternary_vector& vertex, const value_type& ave_value)
        : vertex_(vertex), ave_value_(ave_value) {}
 
    [[nodiscard]] auto vertex() const -> const impl::ternary_vector& {
        return vertex_;
    }
 
    [[nodiscard]] auto ave_value() const -> const value_type& {
        return ave_value_;
    }
 
    [[nodiscard]] auto sample_points() const
        -> std::pair<ternary_vector, ternary_vector> {
        return determine_sample_points(vertex_);
    }
 
    [[nodiscard]] auto dist() const -> value_type {
        auto squared_sum = static_cast<value_type>(0);
        for (index_type i = 0; i < vertex_.dim(); ++i) {
            using std::pow;
            squared_sum +=
                pow(static_cast<value_type>(3), -2 * (vertex_.digits(i) - 1));
        }
        using std::sqrt;
        const auto half = static_cast<value_type>(0.5);
        return half * sqrt(squared_sum);
    }
 
    [[nodiscard]] static auto determine_sample_points(
        const ternary_vector& lowest_vertex)
        -> std::pair<ternary_vector, ternary_vector> {
        auto res = std::make_pair(lowest_vertex, lowest_vertex);
        const auto dim = lowest_vertex.dim();
        for (index_type i = 0; i < dim; ++i) {
            const auto digits = lowest_vertex.digits(i);
            NUM_COLLECT_DEBUG_ASSERT(digits > 0);
            std::uint_fast32_t one_count = 0;
            for (index_type j = 0; j < digits; ++j) {
                if (lowest_vertex(i, j) == ternary_vector::digit_type{1}) {
                    ++one_count;
                }
            }
 
            auto last_digit =  // NOLINTNEXTLINE: false positive
                static_cast<std::int_fast32_t>(lowest_vertex(i, digits - 1));
            ++last_digit;
            constexpr std::uint_fast32_t odd_mask = 1;
            if ((one_count & odd_mask) == odd_mask) {
                res.first(i, digits - 1) =
                    static_cast<ternary_vector::digit_type>(last_digit);
            } else {
                res.second(i, digits - 1) =
                    static_cast<ternary_vector::digit_type>(last_digit);
            }
        }
        normalize_point(res.first);
        normalize_point(res.second);
        return res;
    }
 
private:
    static void normalize_point(ternary_vector& point) {
        for (index_type i = 0; i < point.dim(); ++i) {
            for (index_type j = point.digits(i) - 1; j > 0; --j) {
                if (point(i, j) == ternary_vector::digit_type{3}) {
                    point(i, j) = 0;
                    std::int_fast32_t temp =
                        point(i, j - 1);  // NOLINT: false positive
                    ++temp;
                    point(i, j - 1) =
                        static_cast<ternary_vector::digit_type>(temp);
                }
            }
        }
    }
 
    impl::ternary_vector vertex_;
 
    value_type ave_value_;
};
 
template <base::concepts::real_scalar Value>
class adc_group {
public:
    using value_type = Value;
 
    using rectangle_type = adc_rectangle<value_type>;
 
    using rectangle_pointer_type = std::shared_ptr<rectangle_type>;
 
    explicit adc_group(value_type dist) : dist_(dist) {}
 
    void push(rectangle_pointer_type rect) {
        NUM_COLLECT_DEBUG_ASSERT(rect);
        rects_.push(std::move(rect));
    }
 
    [[nodiscard]] auto min_rect() const -> const rectangle_pointer_type& {
        NUM_COLLECT_DEBUG_ASSERT(!rects_.empty());
        return rects_.top();
    }
 
    [[nodiscard]] auto empty() const -> bool { return rects_.empty(); }
 
    [[nodiscard]] auto pop() -> rectangle_pointer_type {
        NUM_COLLECT_DEBUG_ASSERT(!rects_.empty());
        auto rect = rects_.top();
        rects_.pop();
        return rect;
    }
 
    [[nodiscard]] auto dist() const -> const value_type& { return dist_; }
 
private:
    struct greater {
        [[nodiscard]] auto operator()(const rectangle_pointer_type& left,
            const rectangle_pointer_type& right) const -> bool {
            return left->ave_value() > right->ave_value();
        }
    };
 
    using queue_type = std::priority_queue<rectangle_pointer_type,
        std::vector<rectangle_pointer_type>, greater>;
 
    queue_type rects_{};
 
    value_type dist_;
};
 
}  // namespace impl
 
template <concepts::objective_function ObjectiveFunction>
class adaptive_diagonal_curves
    : public optimizer_base<adaptive_diagonal_curves<ObjectiveFunction>> {
public:
    using this_type = adaptive_diagonal_curves<ObjectiveFunction>;
 
    using objective_function_type = ObjectiveFunction;
 
    using variable_type = typename objective_function_type::variable_type;
 
    using value_type = typename objective_function_type::value_type;
 
    enum class state_type : std::uint8_t {
        none,         
        local,        
        local_last,   
        global,       
        global_last,  
    };
 
    [[nodiscard]] static auto state_name(state_type state) -> std::string_view {
        switch (state) {
        case state_type::none:
            return "none";
        case state_type::local:
            return "local";
        case state_type::local_last:
            return "local (last)";
        case state_type::global:
            return "global";
        case state_type::global_last:
            return "global (last)";
        default:
            return "invalid process";
        }
    }
 
    explicit adaptive_diagonal_curves(
        const objective_function_type& obj_fun = objective_function_type())
        : optimizer_base<adaptive_diagonal_curves<ObjectiveFunction>>(
              adaptive_diagonal_curves_tag),
          value_dict_(obj_fun) {}
 
    void change_objective_function(const objective_function_type& obj_fun) {
        value_dict_.change_objective_function(obj_fun);
    }
 
    void init(const variable_type& lower, const variable_type& upper) {
        value_dict_.init(lower, upper);
        groups_.clear();
        iterations_ = 0;
        state_ = state_type::none;
 
        // Initialize groups_, optimal_value_, optimal_group_index_.
        create_first_rectangle();
 
        // prec_optimal_value_, prec_optimal_group_index, and
        // iterations_in_current_phase_ are initialized in switch_state().
    }
 
    void iterate() {
        switch_state();
 
        switch (state_) {
        case state_type::local:
            iterate_locally();
            break;
        case state_type::local_last:
            iterate_locally_last();
            break;
        case state_type::global:
            iterate_globally();
            break;
        case state_type::global_last:
            iterate_globally_last();
            break;
        default:
            NUM_COLLECT_LOG_AND_THROW(algorithm_failure,
                "invalid state (bug in adaptive_diagonal_curve class)");
        }
 
        ++iterations_;
    }
 
    [[nodiscard]] auto is_stop_criteria_satisfied() const -> bool {
        return evaluations() >= max_evaluations_;
    }
 
    void configure_iteration_logger(
        logging::iterations::iteration_logger<this_type>& iteration_logger)
        const {
        iteration_logger.template append<index_type>(
            "Iter.", &this_type::iterations);
        iteration_logger.template append<index_type>(
            "Eval.", &this_type::evaluations);
        iteration_logger.template append<value_type>(
            "Value", &this_type::opt_value);
        constexpr index_type state_width = 15;
        iteration_logger
            .template append<std::string_view>(
                "State", &this_type::last_state_name)
            ->width(state_width);
    }
 
    [[nodiscard]] auto opt_variable() const -> const variable_type& {
        return value_dict_.opt_variable();
    }
 
    [[nodiscard]] auto opt_value() const -> const value_type& {
        return value_dict_.opt_value();
    }
 
    [[nodiscard]] auto iterations() const noexcept -> index_type {
        return iterations_;
    }
 
    [[nodiscard]] auto evaluations() const noexcept -> index_type {
        return value_dict_.evaluations();
    }
 
    [[nodiscard]] auto last_state() const noexcept -> state_type {
        return state_;
    }
 
    [[nodiscard]] auto last_state_name() const noexcept -> std::string_view {
        return state_name(last_state());
    }
 
    auto max_evaluations(index_type value) -> adaptive_diagonal_curves& {
        NUM_COLLECT_PRECONDITION(value > 0, this->logger(),
            "Maximum number of function evaluations must be a positive "
            "integer.");
 
        max_evaluations_ = value;
        return *this;
    }
 
    auto min_rate_imp(value_type value) -> adaptive_diagonal_curves& {
        NUM_COLLECT_PRECONDITION(value > static_cast<value_type>(0),
            this->logger(),
            "Minimum rate of improvement in the function value required for "
            "potentially optimal rectangles must be a positive value.");
        min_rate_imp_ = value;
        return *this;
    }
 
    auto decrease_rate_bound(value_type value) -> adaptive_diagonal_curves& {
        NUM_COLLECT_PRECONDITION(value > static_cast<value_type>(0),
            this->logger(),
            "Rate of function value used to check whether the function value "
            "decreased in the current phase must be a positive value.");
        decrease_rate_bound_ = value;
        return *this;
    }
 
private:
    using dict_type = impl::adc_sample_dict<objective_function_type>;
 
    using group_type = impl::adc_group<value_type>;
 
    using rectangle_type = typename group_type::rectangle_type;
 
    void create_first_rectangle() {
        const index_type dim = value_dict_.dim();
        auto point = impl::ternary_vector{dim};
        for (index_type i = 0; i < dim; ++i) {
            point.push_back(i, 0);
        }
 
        const auto [lower_vertex, upper_vertex] =
            rectangle_type::determine_sample_points(point);
        const auto lower_vertex_value = value_dict_(lower_vertex);
        const auto upper_vertex_value = value_dict_(upper_vertex);
        const auto ave_value = half * (lower_vertex_value + upper_vertex_value);
        const auto rect = std::make_shared<rectangle_type>(point, ave_value);
 
        groups_.emplace_back(rect->dist());
        groups_.front().push(rect);
        using std::min;
        optimal_value_ = min(lower_vertex_value, upper_vertex_value);
        optimal_group_index_ = 0;
    }
 
    void switch_state() {
        using std::abs;
        switch (state_) {
        case state_type::local:
            ++iterations_in_current_phase_;
            if (iterations_in_current_phase_ > value_dict_.dim()) {
                state_ = state_type::local_last;
            }
            return;
        case state_type::local_last:
            switch_state_on_local_last();
            return;
        case state_type::global:
            if (optimal_value_ <= prec_optimal_value_ -
                    decrease_rate_bound_ * abs(prec_optimal_value_)) {
                state_ = state_type::local;
                prec_optimal_value_ = optimal_value_;
                iterations_in_current_phase_ = 1;
                prec_optimal_group_index_ = optimal_group_index_;
                return;
            }
            ++iterations_in_current_phase_;
            if (iterations_in_current_phase_ >
                util::safe_cast<index_type>((static_cast<std::size_t>(2)
                    << util::safe_cast<std::size_t>(value_dict_.dim())))) {
                state_ = state_type::global_last;
                return;
            }
            return;
        case state_type::global_last:
            if (optimal_value_ <= prec_optimal_value_ -
                    decrease_rate_bound_ * abs(prec_optimal_value_)) {
                state_ = state_type::local;
                prec_optimal_value_ = optimal_value_;
                iterations_in_current_phase_ = 1;
                prec_optimal_group_index_ = optimal_group_index_;
                return;
            }
            state_ = state_type::global;
            iterations_in_current_phase_ = 1;
            prec_optimal_group_index_ = optimal_group_index_;
            return;
        default:
            state_ = state_type::local;
            prec_optimal_value_ = optimal_value_;
            iterations_in_current_phase_ = 1;
            prec_optimal_group_index_ = optimal_group_index_;
        }
    }
 
    void switch_state_on_local_last() {
        using std::abs;
        if (optimal_value_ <= prec_optimal_value_ -
                decrease_rate_bound_ * abs(prec_optimal_value_)) {
            state_ = state_type::local;
            prec_optimal_value_ = optimal_value_;
            iterations_in_current_phase_ = 1;
            prec_optimal_group_index_ = optimal_group_index_;
            return;
        }
 
        const bool is_optimal_smallest =
            (optimal_group_index_ == groups_.size() - 1);
        const bool is_all_smallest =
            std::all_of(std::begin(groups_), std::end(groups_) - 1,
                [](const group_type& group) { return group.empty(); });
        if ((!is_optimal_smallest) || is_all_smallest) {
            state_ = state_type::local;
            iterations_in_current_phase_ = 1;
            prec_optimal_group_index_ = optimal_group_index_;
            return;
        }
 
        state_ = state_type::global;
        prec_optimal_value_ = optimal_value_;
        iterations_in_current_phase_ = 1;
        prec_optimal_group_index_ = optimal_group_index_;
    }
 
    void iterate_locally() {
        const std::size_t max_group_index = std::max<std::size_t>(
            std::max<std::size_t>(prec_optimal_group_index_, 1) - 1,
            min_nonempty_group_index());
        divide_nondominated_rectangles(
            min_nonempty_group_index(), max_group_index);
    }
 
    void iterate_locally_last() {
        const std::size_t max_group_index = std::max<std::size_t>(
            prec_optimal_group_index_, min_nonempty_group_index());
        divide_nondominated_rectangles(
            min_nonempty_group_index(), max_group_index);
    }
 
    void iterate_globally() {
        const std::size_t max_group_index =
            (std::max<std::size_t>(
                 prec_optimal_group_index_, min_nonempty_group_index()) +
                min_nonempty_group_index() + 1) /
            2;
        divide_nondominated_rectangles(
            min_nonempty_group_index(), max_group_index);
    }
 
    void iterate_globally_last() { iterate_locally_last(); }
 
    [[nodiscard]] auto min_nonempty_group_index() const -> std::size_t {
        for (std::size_t i = 0; i < groups_.size(); ++i) {
            if (!groups_[i].empty()) {
                return i;
            }
        }
        NUM_COLLECT_LOG_AND_THROW(precondition_not_satisfied,
            "adaptive_diagonal_curves::init is not called.");
    }
 
    void divide_nondominated_rectangles(
        std::size_t min_group, std::size_t max_group) {
        const auto search_rect =
            determine_nondominated_rectangles(min_group, max_group);
        for (auto iter = std::rbegin(search_rect);
            iter != std::rend(search_rect); ++iter) {
            divide_rectangle(iter->first);
        }
    }
 
    [[nodiscard]] auto determine_nondominated_rectangles(
        std::size_t min_group, std::size_t max_group) const
        -> std::vector<std::pair<std::size_t, value_type>> {
        NUM_COLLECT_DEBUG_ASSERT(min_group <= max_group);
        NUM_COLLECT_DEBUG_ASSERT(max_group < groups_.size());
        NUM_COLLECT_DEBUG_ASSERT(!groups_[min_group].empty());
 
        std::vector<std::pair<std::size_t, value_type>> search_rects;
        search_rects.emplace_back(
            min_group, std::numeric_limits<value_type>::max());
 
        // convex full scan
        for (std::size_t i = min_group + 1; i <= max_group; ++i) {
            if (groups_[i].empty()) {
                continue;
            }
            while (true) {
                const auto& [last_i, last_slope] = search_rects.back();
                const auto slope = calculate_slope(last_i, i);
                if (slope <= last_slope) {
                    search_rects.emplace_back(i, slope);
                    break;
                }
                search_rects.pop_back();
            }
        }
 
        // remove rectangles which won't update optimal value
        using std::abs;
        const auto value_bound =
            optimal_value_ - min_rate_imp_ * abs(optimal_value_);
        for (auto iter = search_rects.begin(); iter != search_rects.end();) {
            const auto& [ind, slope] = *iter;
            if (groups_[ind].min_rect()->ave_value() -
                    slope * groups_[ind].dist() <=
                value_bound) {
                ++iter;
            } else {
                iter = search_rects.erase(iter);
            }
        }
 
        return search_rects;
    }
 
    [[nodiscard]] auto calculate_slope(
        std::size_t group_ind1, std::size_t group_ind2) const -> value_type {
        return (groups_[group_ind1].min_rect()->ave_value() -
                   groups_[group_ind2].min_rect()->ave_value()) /
            (groups_[group_ind1].dist() - groups_[group_ind2].dist());
    }
 
    void divide_rectangle(std::size_t group_ind) {
        impl::ternary_vector vertex = groups_[group_ind].pop()->vertex();
        index_type divided_dim = 1;
        for (; divided_dim < vertex.dim(); ++divided_dim) {
            if (vertex.digits(divided_dim) < vertex.digits(0)) {
                break;
            }
        }
        if (divided_dim == vertex.dim()) {
            divided_dim = 0;
        }
 
        vertex.push_back(divided_dim, 0);
        const auto rect0 = create_rect(vertex, group_ind + 1);
        vertex(divided_dim, vertex.digits(divided_dim) - 1) = 1;
        const auto rect1 = create_rect(vertex, group_ind + 1);
        vertex(divided_dim, vertex.digits(divided_dim) - 1) = 2;
        const auto rect2 = create_rect(vertex, group_ind + 1);
 
        if (groups_.size() == group_ind + 1) {
            groups_.emplace_back(rect0->dist());
        }
        groups_[group_ind + 1].push(rect0);
        groups_[group_ind + 1].push(rect1);
        groups_[group_ind + 1].push(rect2);
    }
 
    [[nodiscard]] auto create_rect(const impl::ternary_vector& vertex,
        std::size_t group_ind) -> std::shared_ptr<rectangle_type> {
        const auto [vertex1, vertex2] =
            rectangle_type::determine_sample_points(vertex);
        const auto value1 = value_dict_(vertex1);
        const auto value2 = value_dict_(vertex2);
        const auto ave_value = half * (value1 + value2);
        if (value1 < optimal_value_) {
            optimal_value_ = value1;
            optimal_group_index_ = group_ind;
        } else if (value1 <= optimal_value_ &&
            group_ind > optimal_group_index_) {
            optimal_group_index_ = group_ind;
        }
        if (value2 < optimal_value_) {
            optimal_value_ = value2;
            optimal_group_index_ = group_ind;
        } else if (value2 <= optimal_value_ &&
            group_ind > optimal_group_index_) {
            optimal_group_index_ = group_ind;
        }
        return std::make_shared<rectangle_type>(vertex, ave_value);
    }
 
    static inline const auto half = static_cast<value_type>(0.5);
 
    dict_type value_dict_;
 
    std::vector<group_type> groups_{};
 
    index_type iterations_{0};
 
    state_type state_{state_type::none};
 
    value_type optimal_value_{};
 
    std::size_t optimal_group_index_{0};
 
    value_type prec_optimal_value_{};
 
    std::size_t prec_optimal_group_index_{0};
 
    index_type iterations_in_current_phase_{0};
 
    static constexpr index_type default_max_evaluations = 10000;
 
    index_type max_evaluations_{default_max_evaluations};
 
    static inline const auto default_min_rate_imp =
        static_cast<value_type>(1e-4);
 
    value_type min_rate_imp_{default_min_rate_imp};
 
    static inline const auto default_decrease_rate_bound =
        static_cast<value_type>(0.01);
 
    value_type decrease_rate_bound_{default_decrease_rate_bound};
};
 
}  // namespace num_collect::opt