downhill_simplex.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 <cstddef>
#include <cstdint>
#include <iterator>
#include <string_view>
#include <vector>
 
#include <Eigen/Core>
 
#include "num_collect/base/index_type.h"
#include "num_collect/logging/iterations/iteration_logger.h"
#include "num_collect/logging/log_tag_view.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 downhill_simplex_tag =
    logging::log_tag_view("num_collect::opt::downhill_simplex");

template <concepts::objective_function ObjectiveFunction>
class downhill_simplex;

template <concepts::multi_variate_objective_function ObjectiveFunction>
class downhill_simplex<ObjectiveFunction>
    : public optimizer_base<downhill_simplex<ObjectiveFunction>> {
public:
    using this_type = 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;

    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 downhill_simplex(
        const objective_function_type& obj_fun = objective_function_type())
        : optimizer_base<downhill_simplex<ObjectiveFunction>>(
              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;
 
        variable_type var = init_var;
        points_.clear();
        values_.clear();
        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_.clear();
        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);
        }
        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 (values_[max_ind] < values_[min_ind]) {
            expand(face_center);
        } else if (values_[max_ind] >= values_[second_max_ind]) {
            contract(face_center);
            if (values_[max_ind] >= values_[second_max_ind]) {
                multi_contract();
            }
        } else {
            // no operation
        }
 
        reorder();
        ++iterations_;
    }

    [[nodiscard]] auto is_stop_criteria_satisfied() const -> bool {
        return (simplex_size() < tol_simplex_size_) ||
            (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);
        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 points_[value_order_.front()];
    }

    [[nodiscard]] auto opt_value() const -> const value_type& {
        return values_[value_order_.front()];
    }

    [[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();
    }

    static inline const auto default_width =
        static_cast<variable_scalar_type>(0.1);

    auto tol_simplex_size(const variable_scalar_type& value)
        -> downhill_simplex& {
        tol_simplex_size_ = value;
        return *this;
    }
 
private:
    [[nodiscard]] auto evaluate_on(const variable_type& variable) {
        obj_fun_.evaluate_on(variable);
        ++evaluations_;
        return obj_fun_.value();
    }

    void reorder() {
        std::sort(std::begin(value_order_), std::end(value_order_),
            [this](std::size_t i, std::size_t j) {
                return values_[i] < 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]);
        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]);
        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]);
        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]);
        }
        process_ = process_type::multiple_contraction;
    }

    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<std::size_t> value_order_{};

    process_type process_{process_type::none};

    index_type iterations_{0};

    index_type evaluations_{0};

    static inline const auto default_tol_simplex_size =
        static_cast<variable_scalar_type>(1e-4);

    variable_scalar_type tol_simplex_size_{default_tol_simplex_size};

    static constexpr index_type default_max_iterations = 1000;

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