gaussian_process_optimizer.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 <cstdlib>
#include <limits>
#include <vector>
 
#include <Eigen/Core>
 
#include "num_collect/base/exception.h"
#include "num_collect/base/get_size.h"
#include "num_collect/base/index_type.h"
#include "num_collect/base/isfinite.h"
#include "num_collect/base/precondition.h"
#include "num_collect/constants/pi.h"  // IWYU pragma: keep
#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/any_objective_function.h"
#include "num_collect/opt/concepts/objective_function.h"
#include "num_collect/opt/dividing_rectangles.h"
#include "num_collect/opt/optimizer_base.h"
#include "num_collect/rbf/gaussian_process_interpolator.h"
#include "num_collect/util/assert.h"
#include "num_collect/util/is_eigen_vector.h"  // IWYU pragma: keep
 
namespace num_collect::opt {

constexpr auto gaussian_process_optimizer_tag =
    logging::log_tag_view("num_collect::opt::gaussian_process_optimizer");

template <concepts::objective_function ObjectiveFunction>
class gaussian_process_optimizer
    : public optimizer_base<gaussian_process_optimizer<ObjectiveFunction>> {
public:
    using this_type = gaussian_process_optimizer<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 gaussian_process_optimizer(
        const objective_function_type& obj_fun = objective_function_type())
        : optimizer_base<gaussian_process_optimizer<ObjectiveFunction>>(
              gaussian_process_optimizer_tag),
          obj_fun_(obj_fun) {
        this->configure_child_algorithm_logger_if_exists(interpolator_);
        this->configure_child_algorithm_logger_if_exists(
            lower_bound_optimizer_);
        lower_bound_optimizer_.max_evaluations(
            default_max_lower_bound_evaluations);
    }

    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;
        dim_ = get_size(lower_);
 
        // First sample.
        obj_fun_.evaluate_on(upper_);
        opt_variable_ = upper_;
        opt_value_ = correct_value_if_needed(obj_fun_.value());
        variables_.resize(1);
        values_.resize(1);
        variables_[0] = opt_variable_;
        values_[0] = opt_value_;
 
        iterations_ = 0;
        evaluations_ = 1;
 
        // Second sample.
        evaluate_on(lower_);
    }

    void iterate() {
        NUM_COLLECT_ASSERT(!variables_.empty());
        NUM_COLLECT_ASSERT(
            variables_.size() == static_cast<std::size_t>(values_.size()));
 
        interpolator_.optimize_length_parameter_scale(variables_, values_);
        interpolator_.compute(variables_, values_);
 
        (void)try_find_and_evaluate_using_lower_bound();
 
        ++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) -> gaussian_process_optimizer& {
        NUM_COLLECT_PRECONDITION(value > 0, this->logger(),
            "Maximum number of function evaluations must be a positive "
            "integer.");
        max_evaluations_ = value;
        return *this;
    }

    auto max_lower_bound_evaluations(index_type value)
        -> gaussian_process_optimizer& {
        NUM_COLLECT_PRECONDITION(value > 0, this->logger(),
            "Maximum number of evaluations of lower bounds must be a positive "
            "integer.");
        lower_bound_optimizer_.max_evaluations(value);
        return *this;
    }
 
private:
    using interpolator_type =
        rbf::gaussian_process_interpolator<value_type(variable_type)>;

    void evaluate_on(const variable_type& variable) {
        obj_fun_.evaluate_on(variable);
        const value_type value = correct_value_if_needed(obj_fun_.value());
        if (value < opt_value_) {
            opt_variable_ = variable;
            opt_value_ = value;
        }
        variables_.push_back(variable);
        values_.conservativeResize(values_.size() + 1);
        values_(values_.size() - 1) = value;
        ++evaluations_;
    }

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

    [[nodiscard]] auto try_find_and_evaluate_using_lower_bound() -> bool {
        const auto lower_bound_function =
            [this](const variable_type& variable) -> value_type {
            using std::log;
            using std::sqrt;
            using std::max;
            using std::pow;
 
            const auto [mean, variance] =
                interpolator_.evaluate_mean_and_variance_on(variable);
            const auto variance_coeff = static_cast<value_type>(2) *
                log(pow(static_cast<value_type>(evaluations_),
                        static_cast<value_type>(0.5) *
                                static_cast<value_type>(dim_) +
                            static_cast<value_type>(2)) *
                    constants::pi<value_type> * constants::pi<value_type> /
                    static_cast<value_type>(0.3));
            return mean -
                sqrt(
                    max(variance_coeff * variance, static_cast<value_type>(0)));
        };
        lower_bound_optimizer_.change_objective_function(lower_bound_function);
 
        lower_bound_optimizer_.init(lower_, upper_);
        lower_bound_optimizer_.solve();
 
        evaluate_on(lower_bound_optimizer_.opt_variable());
 
        return true;
    }

    dividing_rectangles<any_objective_function<value_type(variable_type)>>
        lower_bound_optimizer_{};

    objective_function_type obj_fun_;

    interpolator_type interpolator_{};

    std::vector<variable_type> variables_{};

    Eigen::VectorX<value_type> values_{};

    variable_type lower_{};

    variable_type upper_{};

    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 = 20;

    index_type max_evaluations_{default_max_evaluations};

    static constexpr index_type default_max_lower_bound_evaluations = 100;
};
 
}  // namespace num_collect::opt