annealing_downhill_simplex.h

Go to the documentation of this file.

/*
 * Copyright 2024 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 <random>
#include <string_view>
#include <vector>
 
#include <Eigen/Core>
 
#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/optimizer_base.h"
#include "num_collect/util/safe_cast.h"
 
namespace num_collect::opt {

constexpr auto annealing_downhill_simplex_tag =
    logging::log_tag_view("num_collect::opt::annealing_downhill_simplex");

template <concepts::objective_function ObjectiveFunction>
class annealing_downhill_simplex;

template <concepts::multi_variate_objective_function ObjectiveFunction>
class annealing_downhill_simplex<ObjectiveFunction>
    : public optimizer_base<annealing_downhill_simplex<ObjectiveFunction>> {
public:
    // TODO Refactoring.

    using this_type = annealing_downhill_simplex<ObjectiveFunction>;

    using objective_function_type = ObjectiveFunction;

    using variable_type = typename objective_function_type::variable_type;

    using variable_scalar_type = typename variable_type::Scalar;

    using value_type = typename objective_function_type::value_type;

    using random_number_generator_type = std::mt19937;

    enum class process_type : std::uint8_t {
        none,                      
        reflection,                
        reflection_and_expansion,  
        contraction,               
        multiple_contraction       
    };

    [[nodiscard]] static auto process_name(process_type process)
        -> std::string_view {
        switch (process) {
        case process_type::none:
            return "none";
        case process_type::reflection:
            return "reflection";
        case process_type::reflection_and_expansion:
            return "reflection and expansion";
        case process_type::contraction:
            return "contraction";
        case process_type::multiple_contraction:
            return "multiple contraction";
        default:
            return "invalid process";
        }
    }

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

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

    void init(const variable_type& init_var,
        const variable_scalar_type& width = default_width) {
        dim_ = init_var.size();
        iterations_ = 0;
        evaluations_ = 0;
        iterations_in_current_trial_ = 0;
        temperature_ = highest_temperature_;
        opt_value_ = std::numeric_limits<value_type>::max();
 
        variable_type var = init_var;
        points_.reserve(util::safe_cast<std::size_t>(dim_ + 1));
        values_.reserve(util::safe_cast<std::size_t>(dim_ + 1));
        points_.push_back(var);
        values_.push_back(evaluate_on(var));
        for (index_type i = 0; i < dim_; ++i) {
            const variable_scalar_type val = var[i];
            var[i] += width;
            points_.push_back(var);
            values_.push_back(evaluate_on(var));
            var[i] = val;
        }
 
        value_order_.reserve(util::safe_cast<std::size_t>(dim_ + 1));
        for (std::size_t i = 0; i < util::safe_cast<std::size_t>(dim_ + 1);
            ++i) {
            value_order_.push_back(i);
        }
        fluctuated_values_.resize(util::safe_cast<std::size_t>(dim_ + 1));
        reorder();
    }

    void iterate() {
        const variable_type face_center = calc_face_center();
        const auto min_ind = value_order_.front();
        const auto second_max_ind =
            value_order_[util::safe_cast<std::size_t>(dim_ - 1)];
        const auto max_ind = value_order_.back();
 
        reflect(face_center);
        if (fluctuated_values_[max_ind] < fluctuated_values_[min_ind]) {
            expand(face_center);
        } else if (fluctuated_values_[max_ind] >=
            fluctuated_values_[second_max_ind]) {
            contract(face_center);
            if (fluctuated_values_[max_ind] >=
                fluctuated_values_[second_max_ind]) {
                multi_contract();
            }
        } else {
            // no operation
        }
 
        reorder();
        ++iterations_;
        ++iterations_in_current_trial_;
        if (iterations_in_current_trial_ >= max_iterations_per_trial_) {
            iterations_in_current_trial_ = 0;
        }
        using std::pow;
        temperature_ = highest_temperature_ *
            pow(static_cast<value_type>(1) -
                    static_cast<value_type>(iterations_in_current_trial_) /
                        static_cast<value_type>(max_iterations_per_trial_),
                4);
    }

    [[nodiscard]] auto is_stop_criteria_satisfied() const -> bool {
        return iterations() >= max_iterations_;
    }

    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);
        iteration_logger.template append<variable_scalar_type>(
            "SimplexSize", &this_type::simplex_size);
        iteration_logger.template append<variable_scalar_type>(
            "Temperature", &this_type::temperature);
        constexpr index_type process_width = 26;
        iteration_logger
            .template append<std::string_view>(
                "Process", &this_type::last_process_name)
            ->width(process_width);
    }

    [[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_;
    }

    [[nodiscard]] auto last_process() const noexcept -> process_type {
        return process_;
    }

    [[nodiscard]] auto last_process_name() const noexcept -> std::string_view {
        return process_name(last_process());
    }

    [[nodiscard]] auto simplex_size() const -> variable_scalar_type {
        return (points_[value_order_.front()] - points_[value_order_.back()])
            .norm();
    }

    [[nodiscard]] auto temperature() const noexcept -> value_type {
        return temperature_;
    }

    auto seed(random_number_generator_type::result_type value)
        -> annealing_downhill_simplex& {
        random_number_generator_.seed(value);
        return *this;
    }

    auto highest_temperature(value_type value) -> annealing_downhill_simplex& {
        NUM_COLLECT_PRECONDITION(value > static_cast<value_type>(0),
            this->logger(), "Highest temperature must be a positive number.");
        highest_temperature_ = value;
        return *this;
    }

    auto max_iterations_per_trial(index_type value)
        -> annealing_downhill_simplex& {
        NUM_COLLECT_PRECONDITION(value > 0, this->logger(),
            "Maximum number of iterations in each trial must be a positive "
            "number.");
        max_iterations_per_trial_ = value;
        return *this;
    }

    auto max_iterations(index_type value) -> annealing_downhill_simplex& {
        NUM_COLLECT_PRECONDITION(value > 0, this->logger(),
            "Maximum number of iterations must be a positive number.");
        max_iterations_ = value;
        return *this;
    }

    static inline const auto default_width =
        static_cast<variable_scalar_type>(0.1);
 
private:
    [[nodiscard]] auto evaluate_on(const variable_type& variable) {
        obj_fun_.evaluate_on(variable);
        ++evaluations_;
        const auto& value = obj_fun_.value();
        if (value < opt_value_) {
            opt_value_ = value;
            opt_variable_ = variable;
        }
        return value;
    }

    void reorder() {
        for (std::size_t i = 0; i < values_.size(); ++i) {
            fluctuated_values_[i] = values_[i] + generate_fluctuation();
        }
        std::sort(std::begin(value_order_), std::end(value_order_),
            [this](std::size_t i, std::size_t j) {
                return fluctuated_values_[i] < fluctuated_values_[j];
            });
    }

    [[nodiscard]] auto calc_face_center() const -> variable_type {
        variable_type face_center = variable_type::Zero(dim_);
        for (std::size_t i = 0; i < util::safe_cast<std::size_t>(dim_); ++i) {
            face_center += points_[value_order_[i]];
        }
        face_center /= static_cast<variable_scalar_type>(dim_);
        return face_center;
    }

    static inline const auto twice = static_cast<variable_scalar_type>(2.0);

    static inline const auto half = static_cast<variable_scalar_type>(0.5);

    void reflect(const variable_type& face_center) {
        const std::size_t ind_move = value_order_.back();
        points_[ind_move] = twice * face_center - points_[ind_move];
        values_[ind_move] = evaluate_on(points_[ind_move]);
        fluctuated_values_[ind_move] =
            values_[ind_move] - generate_fluctuation();
        process_ = process_type::reflection;
    }

    void expand(const variable_type& face_center) {
        const std::size_t ind_move = value_order_.back();
        points_[ind_move] = twice * points_[ind_move] - face_center;
        values_[ind_move] = evaluate_on(points_[ind_move]);
        fluctuated_values_[ind_move] =
            values_[ind_move] - generate_fluctuation();
        process_ = process_type::reflection_and_expansion;
    }

    void contract(const variable_type& face_center) {
        const std::size_t ind_move = value_order_.back();
        points_[ind_move] = half * (points_[ind_move] + face_center);
        values_[ind_move] = evaluate_on(points_[ind_move]);
        fluctuated_values_[ind_move] =
            values_[ind_move] - generate_fluctuation();
        process_ = process_type::contraction;
    }

    void multi_contract() {
        const auto& min_point = points_[value_order_.front()];
        for (std::size_t i = 1; i <= util::safe_cast<std::size_t>(dim_); ++i) {
            const std::size_t ind_move = value_order_[i];
            points_[ind_move] = half * (points_[ind_move] + min_point);
            values_[ind_move] = evaluate_on(points_[ind_move]);
            fluctuated_values_[ind_move] =
                values_[ind_move] - generate_fluctuation();
        }
        process_ = process_type::multiple_contraction;
    }

    [[nodiscard]] auto generate_fluctuation() -> value_type {
        using std::log;
        return temperature_ *
            log(static_cast<value_type>(random_number_generator_type::max()) /
                (static_cast<value_type>(random_number_generator_()) +
                    static_cast<value_type>(1)));
    }

    objective_function_type obj_fun_;

    index_type dim_{0};

    std::vector<variable_type, Eigen::aligned_allocator<variable_type>>
        points_{};

    std::vector<value_type> values_{};

    std::vector<value_type> fluctuated_values_{};

    std::vector<std::size_t> value_order_{};

    process_type process_{process_type::none};

    variable_type opt_variable_{};

    value_type opt_value_{};

    index_type iterations_{0};

    index_type evaluations_{0};

    index_type iterations_in_current_trial_{0};

    random_number_generator_type random_number_generator_{
        std::random_device()()};

    value_type temperature_{};

    static inline const auto default_highest_temperature =
        static_cast<value_type>(10);

    value_type highest_temperature_{default_highest_temperature};

    static constexpr index_type default_max_iterations_per_trial = 100;

    index_type max_iterations_per_trial_{default_max_iterations_per_trial};

    static constexpr index_type default_max_iterations = 1000;

    index_type max_iterations_{default_max_iterations};
};
 
}  // namespace num_collect::opt