dividing_rectangles.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 <iterator>
#include <limits>
#include <queue>
#include <tuple>
#include <utility>
#include <vector>
 
#include "num_collect/base/get_size.h"
#include "num_collect/base/index_type.h"
#include "num_collect/base/isfinite.h"
#include "num_collect/base/norm.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/opt/concepts/objective_function.h"
#include "num_collect/opt/optimizer_base.h"
#include "num_collect/util/assert.h"
#include "num_collect/util/is_eigen_vector.h"  // IWYU pragma: keep
#include "num_collect/util/safe_cast.h"
 
namespace num_collect::opt {

constexpr auto dividing_rectangles_tag =
    logging::log_tag_view("num_collect::opt::dividing_rectangles");

template <concepts::objective_function ObjectiveFunction>
class dividing_rectangles
    : public optimizer_base<dividing_rectangles<ObjectiveFunction>> {
public:
    using this_type = dividing_rectangles<ObjectiveFunction>;

    using objective_function_type = ObjectiveFunction;

    using variable_type = typename objective_function_type::variable_type;

    using value_type = typename objective_function_type::value_type;

    explicit dividing_rectangles(
        const objective_function_type& obj_fun = objective_function_type())
        : optimizer_base<dividing_rectangles<ObjectiveFunction>>(
              dividing_rectangles_tag),
          obj_fun_(obj_fun) {}

    void change_objective_function(const objective_function_type& obj_fun) {
        obj_fun_ = obj_fun;
    }

    void init(const variable_type& lower, const variable_type& upper) {
        if constexpr (is_eigen_vector_v<variable_type>) {
            NUM_COLLECT_PRECONDITION(lower.size() == upper.size(),
                this->logger(),
                "Lower and upper limits must have the same size.");
        }
        lower_ = lower;
        upper_ = upper;
        width_ = upper_ - lower_;
        dim_ = get_size(lower_);
 
        const auto half = static_cast<value_type>(0.5);
        const variable_type first_var = half * width_ + lower_;
        obj_fun_.evaluate_on(first_var);
        opt_variable_ = first_var;
        opt_value_ = correct_value_if_needed(obj_fun_.value());
 
        iterations_ = 0;
        evaluations_ = 1;
 
        rects_.clear();
        rects_.emplace_back();
        if constexpr (is_eigen_vector_v<variable_type>) {
            rects_[0].push(rectangle(variable_type::Zero(dim_),
                variable_type::Ones(dim_), opt_value_));
        } else {
            rects_[0].push(rectangle(static_cast<variable_type>(0),
                static_cast<variable_type>(1), opt_value_));
        }
    }

    void iterate() {
        const auto divided_rects = determine_rects();
        for (auto iter = std::rbegin(divided_rects);
            iter != std::rend(divided_rects); ++iter) {
            divide_rect(iter->first);
        }
        ++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);
    }

    [[nodiscard]] auto opt_variable() const -> const variable_type& {
        return opt_variable_;
    }

    [[nodiscard]] auto opt_value() const -> const value_type& {
        return opt_value_;
    }

    [[nodiscard]] auto iterations() const noexcept -> index_type {
        return iterations_;
    }

    [[nodiscard]] auto evaluations() const noexcept -> index_type {
        return evaluations_;
    }

    auto max_evaluations(index_type value) -> dividing_rectangles& {
        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) -> dividing_rectangles& {
        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;
    }
 
private:
    class rectangle {
    public:
        rectangle(const variable_type& lower, const variable_type& upper,
            const value_type& value)
            : lower_(lower),
              upper_(upper),
              is_divided_(util::safe_cast<std::size_t>(get_size(lower)), false),
              dist_(norm(lower - upper) * half),
              value_(value) {}

        [[nodiscard]] auto lower() const -> const variable_type& {
            return lower_;
        }

        [[nodiscard]] auto upper() const -> const variable_type& {
            return upper_;
        }

        [[nodiscard]] auto is_divided() const -> const std::vector<bool>& {
            return is_divided_;
        }

        [[nodiscard]] auto dist() const -> const value_type& { return dist_; }

        [[nodiscard]] auto value() const -> const value_type& { return value_; }

        void divide_in_place(index_type dim, value_type lower, value_type upper,
            value_type value) {
            NUM_COLLECT_DEBUG_ASSERT(lower < upper);
            NUM_COLLECT_DEBUG_ASSERT(
                !is_divided_[static_cast<std::size_t>(dim)]);
            if constexpr (is_eigen_vector_v<variable_type>) {
                lower_(dim) = lower;
                upper_(dim) = upper;
            } else {
                lower_ = lower;
                upper_ = upper;
            }
            dist_ = norm(lower_ - upper_) * half;
            value_ = value;
            is_divided_[static_cast<std::size_t>(dim)] = true;
            if (std::all_of(std::begin(is_divided_), std::end(is_divided_),
                    [](bool elem) { return elem; })) {
                std::fill(
                    std::begin(is_divided_), std::end(is_divided_), false);
            }
        }

        [[nodiscard]] auto divide(index_type dim, value_type lower,
            value_type upper, value_type value) const -> rectangle {
            auto copy = rectangle(*this);
            copy.divide_in_place(dim, lower, upper, value);
            return copy;
        }
 
    private:
        static inline const auto half = static_cast<value_type>(0.5);

        variable_type lower_;

        variable_type upper_;

        std::vector<bool> is_divided_;

        value_type dist_;

        value_type value_;
    };

    struct greater_rectangle {
        [[nodiscard]] auto operator()(
            const rectangle& left, const rectangle& right) const -> bool {
            return left.value() > right.value();
        }
    };

    [[nodiscard]] auto evaluate_on(const variable_type& variable)
        -> value_type {
        variable_type actual_var;
        if constexpr (is_eigen_vector_v<variable_type>) {
            actual_var = variable.cwiseProduct(width_) + lower_;
        } else {
            actual_var = variable * width_ + lower_;
        }
        obj_fun_.evaluate_on(actual_var);
        const value_type value = correct_value_if_needed(obj_fun_.value());
        if (value < opt_value_) {
            opt_variable_ = actual_var;
            opt_value_ = value;
        }
        ++evaluations_;
        return value;
    }

    [[nodiscard]] auto determine_rects() const
        -> std::vector<std::pair<std::size_t, value_type>> {
        std::vector<std::pair<std::size_t, value_type>> search_rects;
        std::size_t ind = 0;
        // find the largest rectangle
        for (; ind < rects_.size(); ++ind) {
            if (!rects_[ind].empty()) {
                break;
            }
        }
        search_rects.emplace_back(ind, std::numeric_limits<value_type>::max());
        ++ind;
        // convex full scan
        for (; ind < rects_.size(); ++ind) {
            if (rects_[ind].empty()) {
                continue;
            }
            while (true) {
                const auto& [last_ind, last_slope] = search_rects.back();
                const auto slope =
                    calculate_slope(rects_[last_ind].top(), rects_[ind].top());
                if (slope <= last_slope) {
                    search_rects.emplace_back(ind, slope);
                    break;
                }
                search_rects.pop_back();
            }
        }
        // remove rectangles which won't update optimal value
        using std::abs;
        const auto value_bound = opt_value_ - min_rate_imp_ * abs(opt_value_);
        for (auto iter = search_rects.begin(); iter != search_rects.end();) {
            const auto& [ind, slope] = *iter;
            if (rects_[ind].top().value() - slope * rects_[ind].top().dist() <=
                value_bound) {
                ++iter;
            } else {
                iter = search_rects.erase(iter);
            }
        }
        return search_rects;
    }

    [[nodiscard]] static auto calculate_slope(const rectangle& larger_rect,
        const rectangle& smaller_rect) -> value_type {
        return (larger_rect.value() - smaller_rect.value()) /
            (larger_rect.dist() - smaller_rect.dist());
    }

    void divide_rect(std::size_t index) {
        if (index + 1 == rects_.size()) {
            rects_.emplace_back();
        }
        rectangle origin = std::move(rects_[index].top());
        rects_[index].pop();
 
        const auto [divided_dim, lower_value, upper_value] =
            determine_divided_dimension(origin);
        value_type divided_lowest;
        value_type divided_highest;
        if constexpr (is_eigen_vector_v<variable_type>) {
            divided_lowest = origin.lower()(divided_dim);
            divided_highest = origin.upper()(divided_dim);
        } else {
            divided_lowest = origin.lower();
            divided_highest = origin.upper();
        }
        const auto [divided_lower, divided_upper] =
            separate_section(divided_lowest, divided_highest);
 
        rects_[index + 1].push(origin.divide(
            divided_dim, divided_lowest, divided_lower, lower_value));
        rects_[index + 1].push(origin.divide(
            divided_dim, divided_upper, divided_highest, upper_value));
        origin.divide_in_place(
            divided_dim, divided_lower, divided_upper, origin.value());
        rects_[index + 1].push(std::move(origin));
    }

    [[nodiscard]] auto determine_divided_dimension(const rectangle& rect)
        -> std::tuple<index_type, value_type, value_type> {
        value_type min_value = std::numeric_limits<value_type>::max();
        index_type dim = 0;
        value_type dim_lower_value{};
        value_type dim_upper_value{};
        for (index_type i = 0; i < dim_; ++i) {
            if (rect.is_divided()[util::safe_cast<std::size_t>(i)]) {
                continue;
            }
 
            value_type width;
            if constexpr (is_eigen_vector_v<variable_type>) {
                width = rect.upper()[i] - rect.lower()[i];
            } else {
                width = rect.upper() - rect.lower();
            }
            static const auto one_third = static_cast<value_type>(1.0 / 3.0);
            const value_type diff = width * one_third;
 
            static const auto half = static_cast<value_type>(0.5);
            variable_type center = (rect.lower() + rect.upper()) * half;
            value_type origin_center;
            if constexpr (is_eigen_vector_v<variable_type>) {
                origin_center = center[i];
            } else {
                origin_center = center;
            }
 
            value_type lower_value;
            value_type upper_value;
            if constexpr (is_eigen_vector_v<variable_type>) {
                center[i] = origin_center - diff;
                lower_value = evaluate_on(center);
                center[i] = origin_center + diff;
                upper_value = evaluate_on(center);
            } else {
                lower_value = evaluate_on(origin_center - diff);
                upper_value = evaluate_on(origin_center + diff);
            }
 
            if (lower_value < min_value) {
                min_value = lower_value;
                dim = i;
                dim_lower_value = lower_value;
                dim_upper_value = upper_value;
            }
            if (upper_value < min_value) {
                min_value = upper_value;
                dim = i;
                dim_lower_value = lower_value;
                dim_upper_value = upper_value;
            }
        }
 
        return {dim, dim_lower_value, dim_upper_value};
    }

    [[nodiscard]] static auto separate_section(value_type lowest,
        value_type highest) -> std::pair<value_type, value_type> {
        const auto width = highest - lowest;
        static const auto one_third = static_cast<value_type>(1.0 / 3.0);
        static const auto two_thirds = static_cast<value_type>(2.0 / 3.0);
        return {lowest + one_third * width, lowest + two_thirds * width};
    }

    [[nodiscard]] static auto correct_value_if_needed(value_type value) noexcept
        -> value_type {
        constexpr auto safe_limit =
            std::numeric_limits<value_type>::max() * 1e-2;
        if (!isfinite(value) || value > safe_limit) {
            return safe_limit;
        }
        return value;
    }

    objective_function_type obj_fun_;

    std::vector<std::priority_queue<rectangle, std::vector<rectangle>,
        greater_rectangle>>
        rects_{};

    variable_type lower_{};

    variable_type upper_{};

    variable_type width_{};

    index_type dim_{0};

    variable_type opt_variable_{};

    value_type opt_value_{};

    index_type iterations_{0};

    index_type evaluations_{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};
};
 
}  // namespace num_collect::opt