tanh_finite_integrator.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 <cmath>  // IWYU pragma: keep
#include <cstddef>
#include <type_traits>
#include <vector>
 
#include "num_collect/base/concepts/invocable_as.h"
#include "num_collect/base/concepts/real_scalar.h"
#include "num_collect/base/exception.h"
#include "num_collect/base/index_type.h"
#include "num_collect/base/isfinite.h"
#include "num_collect/base/precondition.h"
#include "num_collect/constants/half.h"
#include "num_collect/constants/one.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/util/kahan_adder.h"
 
namespace num_collect::integration {

constexpr auto tanh_finite_integrator_tag =
    logging::log_tag_view("num_collect::integration::tanh_finite_integrator");
 
namespace impl {

template <typename Variable>
struct tanh_finite_integrator_traits;

template <>
struct tanh_finite_integrator_traits<float> {
public:
    static constexpr double default_max_point = 10.0F;
};

template <>
struct tanh_finite_integrator_traits<double> {
public:
    static constexpr double default_max_point = 20.0;
};
 
}  // namespace impl

template <typename Signature>
class tanh_finite_integrator;

template <typename Result, base::concepts::real_scalar Variable>
class tanh_finite_integrator<Result(Variable)> : public logging::logging_mixin {
public:
    using variable_type = std::decay_t<Variable>;

    using result_type = std::decay_t<Result>;

    tanh_finite_integrator()
        : logging::logging_mixin(tanh_finite_integrator_tag) {
        calculate_coefficients();
    }

    template <base::concepts::invocable_as<result_type(variable_type)> Function>
    [[nodiscard]] auto integrate(const Function& function, variable_type left,
        variable_type right) const -> result_type {
        const variable_type center =
            constants::half<variable_type> * (left + right);
        const variable_type width = right - left;
 
        constexpr auto center_weight_rate = static_cast<variable_type>(0.5);
        const variable_type center_weight = width * center_weight_rate;
        util::kahan_adder<result_type> sum;
        sum += function(center) * center_weight;
 
        for (index_type i = 0; i < points_; ++i) {
            const variable_type variable_distance =
                width * variable_rate_list_[static_cast<std::size_t>(i)];
            const variable_type weight =
                width * weight_rate_list_[static_cast<std::size_t>(i)];
 
            const variable_type var_plus = right - variable_distance;
            const variable_type var_minus = left + variable_distance;
            const result_type function_values =
                function(var_plus) + function(var_minus);
            if (!base::isfinite(function_values)) [[unlikely]] {
                NUM_COLLECT_LOG_WARNING(this->logger(),
                    "A function value was not a finite value. "
                    "Stopped numerical integration.");
                break;
            }
            sum += function_values * weight;
        }
 
        return sum.sum() * interval_;
    }

    template <base::concepts::invocable_as<result_type(variable_type)>
                  LeftBoundaryFunction,
        base::concepts::invocable_as<result_type(variable_type)>
            RightBoundaryFunction>
    [[nodiscard]] auto integrate(
        const LeftBoundaryFunction& left_boundary_function,
        const RightBoundaryFunction& right_boundary_function,
        variable_type left, variable_type right) const -> result_type {
        using constants::half;
 
        const variable_type width = right - left;
        const variable_type half_width = half<variable_type> * width;
 
        constexpr auto center_weight_rate = static_cast<variable_type>(0.5);
        const variable_type center_weight = width * center_weight_rate;
        util::kahan_adder<result_type> sum;
        sum += left_boundary_function(half_width) * center_weight;
 
        for (index_type i = 0; i < points_; ++i) {
            const variable_type variable_distance =
                width * variable_rate_list_[static_cast<std::size_t>(i)];
            const variable_type weight =
                width * weight_rate_list_[static_cast<std::size_t>(i)];
 
            const result_type function_values =
                left_boundary_function(variable_distance) +
                right_boundary_function(-variable_distance);
            if (!base::isfinite(function_values)) [[unlikely]] {
                NUM_COLLECT_LOG_WARNING(this->logger(),
                    "A function value was not a finite value. "
                    "Stopped numerical integration.");
                break;
            }
            sum += function_values * weight;
        }
 
        return sum.sum() * interval_;
    }

    template <base::concepts::invocable_as<result_type(variable_type)> Function>
    [[nodiscard]] auto operator()(const Function& function, variable_type left,
        variable_type right) const -> result_type {
        return integrate(function, left, right);
    }

    template <base::concepts::invocable_as<result_type(variable_type)>
                  LeftBoundaryFunction,
        base::concepts::invocable_as<result_type(variable_type)>
            RightBoundaryFunction>
    [[nodiscard]] auto operator()(
        const LeftBoundaryFunction& left_boundary_function,
        const RightBoundaryFunction& right_boundary_function,
        variable_type left, variable_type right) const -> result_type {
        return integrate(
            left_boundary_function, right_boundary_function, left, right);
    }

    auto max_point(variable_type val) -> tanh_finite_integrator& {
        NUM_COLLECT_PRECONDITION(val > static_cast<variable_type>(0),
            this->logger(), "Maximum point must be a positive value.");
        max_point_ = val;
        calculate_coefficients();
        return *this;
    }

    auto points(index_type val) -> tanh_finite_integrator& {
        NUM_COLLECT_PRECONDITION(val > 0, this->logger(),
            "Number of points must a positive integer.");
        points_ = val;
        calculate_coefficients();
        return *this;
    }
 
private:
    void calculate_coefficients() {
        using constants::one;
 
        variable_rate_list_.clear();
        variable_rate_list_.reserve(static_cast<std::size_t>(points_));
        weight_rate_list_.clear();
        weight_rate_list_.reserve(static_cast<std::size_t>(points_));
 
        interval_ = max_point_ / static_cast<variable_type>(points_);
        for (index_type i = 1; i <= points_; ++i) {
            const variable_type changed_variable =
                interval_ * static_cast<variable_type>(i);
            const variable_type exp_value =
                std::exp(static_cast<variable_type>(-2) * changed_variable);
            const variable_type denominator = one<variable_type> + exp_value;
            variable_rate_list_.push_back(exp_value / denominator);
            weight_rate_list_.push_back(static_cast<variable_type>(2) *
                exp_value / (denominator * denominator));
        }
    }

    static constexpr auto default_max_point =
        impl::tanh_finite_integrator_traits<Variable>::default_max_point;

    variable_type max_point_{default_max_point};

    static constexpr index_type default_points = 50;

    index_type points_{default_points};

    variable_type interval_{};

    std::vector<variable_type> variable_rate_list_{};

    std::vector<variable_type> weight_rate_list_{};
};
 
}  // namespace num_collect::integration