rbf_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 <type_traits>  // IWYU pragma: keep
#include <vector>
 
#include <Eigen/Core>
 
#include "num_collect/base/index_type.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/get_default_scalar_type.h"
#include "num_collect/rbf/impl/get_variable_type.h"
#include "num_collect/rbf/impl/kernel_matrix_solver.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/rbfs/gaussian_rbf.h"
#include "num_collect/rbf/rbfs/wendland_csrbf.h"
 
namespace num_collect::rbf {

constexpr auto rbf_interpolator_tag =
    logging::log_tag_view("num_collect::rbf::rbf_interpolator");

template <typename FunctionSignature,
    concepts::rbf RBF =
        rbfs::gaussian_rbf<impl::get_default_scalar_type<FunctionSignature>>,
    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_interpolator;

template <typename Variable, typename FunctionValue, concepts::rbf RBF,
    kernel_matrix_type KernelMatrixType,
    concepts::distance_function DistanceFunction,
    concepts::length_parameter_calculator LengthParameterCalculator>
    requires std::is_same_v<Variable,
                 typename DistanceFunction::variable_type> &&
    std::is_same_v<typename DistanceFunction::value_type,
        typename RBF::scalar_type> &&
    std::is_same_v<DistanceFunction,
        typename LengthParameterCalculator::distance_function_type>
class rbf_interpolator<FunctionValue(Variable), RBF, 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 kernel_matrix_solver_type =
        impl::kernel_matrix_solver<kernel_value_type, function_value_type,
            KernelMatrixType, uses_global_length_parameter>;

    using kernel_matrix_type =
        typename kernel_matrix_solver_type::kernel_matrix_type;

    using function_value_vector_type = Eigen::VectorX<function_value_type>;

    static constexpr index_type default_max_mle_evaluations = 20;

    explicit rbf_interpolator(
        distance_function_type distance_function = distance_function_type(),
        rbf_type rbf = rbf_type())
        : logging::logging_mixin(rbf_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) {
        compute_kernel_matrix(distance_function_, rbf_,
            length_parameter_calculator_, variables, kernel_matrix_);
        kernel_matrix_solver_.compute(kernel_matrix_, function_values);
        kernel_matrix_solver_.solve(coeffs_, reg_param, function_values);
        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 += coeffs_(static_cast<index_type>(i)) *
                        rbf_(distance_rate);
                }
            } else {
                value +=
                    coeffs_(static_cast<index_type>(i)) * rbf_(distance_rate);
            }
        }
        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
    {
        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);
            this->length_parameter_calculator_.scale(scale);
            compute_kernel_matrix(this->distance_function_, this->rbf_,
                this->length_parameter_calculator_, variables,
                this->kernel_matrix_);
            this->kernel_matrix_solver_.compute(
                this->kernel_matrix_, function_values);
            return std::log10(
                this->kernel_matrix_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 coeffs() const noexcept
        -> const function_value_vector_type& {
        return coeffs_;
    }
 
protected:
    [[nodiscard]] auto variables() const noexcept
        -> const std::vector<variable_type>& {
        return *variables_;
    }

    [[nodiscard]] auto kernel_matrix_solver() const noexcept
        -> const kernel_matrix_solver_type& {
        return kernel_matrix_solver_;
    }
 
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_{};

    kernel_matrix_type kernel_matrix_{};

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

    kernel_matrix_solver_type kernel_matrix_solver_{};

    function_value_vector_type coeffs_{};
};

template <typename FunctionSignature,
    concepts::rbf RBF =
        rbfs::gaussian_rbf<impl::get_default_scalar_type<FunctionSignature>>,
    kernel_matrix_type KernelMatrixType = kernel_matrix_type::dense,
    concepts::distance_function DistanceFunction =
        distance_functions::euclidean_distance_function<
            impl::get_variable_type_t<FunctionSignature>>>
using local_rbf_interpolator =
    rbf_interpolator<FunctionSignature, RBF, KernelMatrixType, DistanceFunction,
        length_parameter_calculators::local_length_parameter_calculator<
            DistanceFunction>>;

template <typename FunctionSignature,
    concepts::rbf RBF =
        rbfs::gaussian_rbf<impl::get_default_scalar_type<FunctionSignature>>,
    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_interpolator =
    rbf_interpolator<FunctionSignature, RBF, KernelMatrixType, DistanceFunction,
        length_parameter_calculators::global_length_parameter_calculator<
            DistanceFunction>>;

template <typename FunctionSignature,
    concepts::rbf RBF = rbfs::wendland_csrbf<
        impl::get_default_scalar_type<FunctionSignature>, 3, 1>,
    concepts::distance_function DistanceFunction =
        distance_functions::euclidean_distance_function<
            impl::get_variable_type_t<FunctionSignature>>>
using local_csrbf_interpolator = rbf_interpolator<FunctionSignature, RBF,
    kernel_matrix_type::sparse, DistanceFunction,
    length_parameter_calculators::local_length_parameter_calculator<
        DistanceFunction>>;
 
}  // namespace num_collect::rbf