rbf_polynomial_interpolator.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 <cstddef>
#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/precondition.h"
#include "num_collect/logging/log_tag_view.h"
#include "num_collect/logging/logging_macros.h"
#include "num_collect/logging/logging_mixin.h"
#include "num_collect/opt/dividing_rectangles.h"
#include "num_collect/opt/function_object_wrapper.h"
#include "num_collect/rbf/compute_kernel_matrix.h"
#include "num_collect/rbf/concepts/csrbf.h"
#include "num_collect/rbf/concepts/distance_function.h"
#include "num_collect/rbf/concepts/length_parameter_calculator.h"
#include "num_collect/rbf/concepts/rbf.h"
#include "num_collect/rbf/distance_functions/euclidean_distance_function.h"
#include "num_collect/rbf/impl/general_spline_equation_solver.h"
#include "num_collect/rbf/impl/get_default_scalar_type.h"
#include "num_collect/rbf/impl/get_variable_type.h"
#include "num_collect/rbf/kernel_matrix_type.h"
#include "num_collect/rbf/length_parameter_calculators/global_length_parameter_calculator.h"
#include "num_collect/rbf/length_parameter_calculators/local_length_parameter_calculator.h"
#include "num_collect/rbf/polynomial_calculator.h"
#include "num_collect/rbf/rbfs/gaussian_rbf.h"
 
namespace num_collect::rbf {

constexpr auto rbf_polynomial_interpolator_tag =
    logging::log_tag_view("num_collect::rbf::rbf_polynomial_interpolator");

template <typename FunctionSignature,
    concepts::rbf RBF =
        rbfs::gaussian_rbf<impl::get_default_scalar_type<FunctionSignature>>,
    index_type PolynomialDegree = 1,
    kernel_matrix_type KernelMatrixType = kernel_matrix_type::dense,
    concepts::distance_function DistanceFunction =
        distance_functions::euclidean_distance_function<
            impl::get_variable_type_t<FunctionSignature>>,
    concepts::length_parameter_calculator LengthParameterCalculator =
        length_parameter_calculators::local_length_parameter_calculator<
            DistanceFunction>>
class rbf_polynomial_interpolator;

template <typename Variable, typename FunctionValue, concepts::rbf RBF,
    index_type PolynomialDegree, kernel_matrix_type KernelMatrixType,
    concepts::distance_function DistanceFunction,
    concepts::length_parameter_calculator LengthParameterCalculator>
class rbf_polynomial_interpolator<FunctionValue(Variable), RBF,
    PolynomialDegree, KernelMatrixType, DistanceFunction,
    LengthParameterCalculator> : public logging::logging_mixin {
public:
    using variable_type = Variable;

    using function_value_type = FunctionValue;

    using rbf_type = RBF;

    using distance_function_type = DistanceFunction;

    using length_parameter_calculator_type = LengthParameterCalculator;

    static constexpr bool uses_global_length_parameter =
        length_parameter_calculator_type::uses_global_length_parameter;

    using kernel_value_type = typename DistanceFunction::value_type;

    using equation_solver_type =
        impl::general_spline_equation_solver<kernel_value_type,
            function_value_type, KernelMatrixType,
            uses_global_length_parameter>;

    using kernel_matrix_type =
        typename equation_solver_type::kernel_matrix_type;

    using polynomial_matrix_type =
        typename equation_solver_type::additional_matrix_type;

    using function_value_vector_type = Eigen::VectorX<function_value_type>;

    static constexpr index_type default_max_mle_evaluations = 20;

    explicit rbf_polynomial_interpolator(
        distance_function_type distance_function = distance_function_type(),
        rbf_type rbf = rbf_type())
        : logging::logging_mixin(rbf_polynomial_interpolator_tag),
          distance_function_(std::move(distance_function)),
          rbf_(std::move(rbf)) {}

    void compute(const std::vector<variable_type>& variables,
        const function_value_vector_type& function_values) {
        const auto num_variables = static_cast<index_type>(variables.size());
        NUM_COLLECT_PRECONDITION(
            num_variables > 0, this->logger(), "Variables must be given.");
        const index_type num_dimensions = base::get_size(variables.front());
        polynomial_calculator_.prepare(num_dimensions);
 
        compute_kernel_matrix(distance_function_, rbf_,
            length_parameter_calculator_, variables, kernel_matrix_);
        polynomial_calculator_.compute_polynomial_term_matrix(
            variables, polynomial_matrix_);
        equation_solver_.compute(
            kernel_matrix_, polynomial_matrix_, function_values);
        equation_solver_.solve(kernel_coeffs_, polynomial_coeffs_, reg_param);
        variables_ = &variables;
    }

    [[nodiscard]] auto interpolate(const variable_type& variable) const
        -> function_value_type {
        auto value = static_cast<function_value_type>(0);
 
        for (std::size_t i = 0; i < variables_->size(); ++i) {
            const kernel_value_type distance_rate =
                distance_function_(variable, (*variables_)[i]) /
                length_parameter_calculator_.length_parameter_at(
                    static_cast<index_type>(i));
            if constexpr (concepts::csrbf<rbf_type>) {
                if (distance_rate < rbf_type::support_boundary()) {
                    value += kernel_coeffs_(static_cast<index_type>(i)) *
                        rbf_(distance_rate);
                }
            } else {
                value += kernel_coeffs_(static_cast<index_type>(i)) *
                    rbf_(distance_rate);
            }
        }
 
        value += polynomial_calculator_.evaluate_polynomial_for_variable(
            variable, polynomial_coeffs_);
 
        return value;
    }

    void fix_length_parameter_scale(kernel_value_type value) {
        length_parameter_calculator_.scale(value);
    }

    void optimize_length_parameter_scale(
        const std::vector<variable_type>& variables,
        const function_value_vector_type& function_values,
        index_type max_mle_evaluations = default_max_mle_evaluations)
        requires uses_global_length_parameter
    {
        const auto num_variables = static_cast<index_type>(variables.size());
        NUM_COLLECT_PRECONDITION(
            num_variables > 0, this->logger(), "Variables must be given.");
        const index_type num_dimensions = base::get_size(variables.front());
        polynomial_calculator_.prepare(num_dimensions);
 
        static constexpr auto base = static_cast<kernel_value_type>(10);
        auto objective_function =
            [this, &variables, &function_values](
                kernel_value_type log_scale) -> kernel_value_type {
            const kernel_value_type scale = std::pow(base, log_scale);
            length_parameter_calculator_.scale(scale);
            compute_kernel_matrix(distance_function_, rbf_,
                length_parameter_calculator_, variables, kernel_matrix_);
            polynomial_calculator_.compute_polynomial_term_matrix(
                variables, polynomial_matrix_);
            equation_solver_.compute(
                kernel_matrix_, polynomial_matrix_, function_values);
            return std::log10(equation_solver_.calc_mle_objective(reg_param));
        };
 
        using objective_function_object_type = decltype(objective_function);
        using objective_function_wrapper_type =
            opt::function_object_wrapper<kernel_value_type(kernel_value_type),
                objective_function_object_type>;
        using optimizer_type =
            opt::dividing_rectangles<objective_function_wrapper_type>;
 
        optimizer_type optimizer{
            objective_function_wrapper_type{objective_function}};
        configure_child_algorithm_logger_if_exists(optimizer);
        constexpr auto lower_boundary = static_cast<kernel_value_type>(-1);
        constexpr auto upper_boundary = static_cast<kernel_value_type>(2);
        optimizer.max_evaluations(max_mle_evaluations);
        optimizer.init(lower_boundary, upper_boundary);
        optimizer.solve();
 
        const kernel_value_type log_scale = optimizer.opt_variable();
        const kernel_value_type scale = std::pow(base, log_scale);
        NUM_COLLECT_LOG_DEBUG(this->logger(),
            "Selected an optimized scale of length parameters: {}", scale);
        this->length_parameter_calculator_.scale(scale);
    }

    [[nodiscard]] auto distance_function() const noexcept
        -> const distance_function_type& {
        return distance_function_;
    }

    [[nodiscard]] auto rbf() const noexcept -> const rbf_type& { return rbf_; }

    [[nodiscard]] auto length_parameter_calculator() const noexcept
        -> const length_parameter_calculator_type& {
        return length_parameter_calculator_;
    }

    [[nodiscard]] auto kernel_coeffs() const noexcept
        -> const function_value_vector_type& {
        return kernel_coeffs_;
    }

    [[nodiscard]] auto polynomial_coeffs() const noexcept
        -> const function_value_vector_type& {
        return polynomial_coeffs_;
    }
 
private:
    static constexpr auto reg_param = static_cast<kernel_value_type>(0);
 
    // TODO Implementation of regularization.

    distance_function_type distance_function_;

    rbf_type rbf_;

    length_parameter_calculator_type length_parameter_calculator_{};

    polynomial_calculator<variable_type, PolynomialDegree>
        polynomial_calculator_{};

    kernel_matrix_type kernel_matrix_{};

    polynomial_matrix_type polynomial_matrix_{};

    const std::vector<variable_type>* variables_{nullptr};

    equation_solver_type equation_solver_{};

    function_value_vector_type kernel_coeffs_{};

    function_value_vector_type polynomial_coeffs_{};
};

template <typename FunctionSignature,
    concepts::rbf RBF =
        rbfs::gaussian_rbf<impl::get_default_scalar_type<FunctionSignature>>,
    index_type PolynomialDegree = 1,
    kernel_matrix_type KernelMatrixType = kernel_matrix_type::dense,
    concepts::distance_function DistanceFunction =
        distance_functions::euclidean_distance_function<
            impl::get_variable_type_t<FunctionSignature>>>
using global_rbf_polynomial_interpolator =
    rbf_polynomial_interpolator<FunctionSignature, RBF, PolynomialDegree,
        KernelMatrixType, DistanceFunction,
        length_parameter_calculators::global_length_parameter_calculator<
            DistanceFunction>>;
 
}  // namespace num_collect::rbf