de_finite_integrator.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 <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/index_type.h"
#include "num_collect/base/isfinite.h"
#include "num_collect/base/precondition.h"
#include "num_collect/constants/half.h"  // IWYU pragma: keep
#include "num_collect/constants/one.h"   // IWYU pragma: keep
#include "num_collect/constants/pi.h"    // IWYU pragma: keep
#include "num_collect/constants/two.h"   // IWYU pragma: keep
#include "num_collect/constants/zero.h"  // IWYU pragma: keep
#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 de_finite_integrator_tag =
    logging::log_tag_view("num_collect::integration::de_finite_integrator");

template <typename Signature>
class de_finite_integrator;

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

    using result_type = std::decay_t<Result>;

    de_finite_integrator() : logging::logging_mixin(de_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 {
        using constants::pi;
 
        const variable_type center =
            constants::half<variable_type> * (left + right);
        const variable_type width = right - left;
 
        constexpr variable_type center_weight_rate =
            pi<variable_type> / static_cast<variable_type>(4);
        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;
        using constants::pi;
 
        const variable_type width = right - left;
        const variable_type half_width = half<variable_type> * width;
 
        constexpr variable_type center_weight_rate =
            pi<variable_type> / static_cast<variable_type>(4);
        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) -> de_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) -> de_finite_integrator& {
        NUM_COLLECT_PRECONDITION(val > 0, this->logger(),
            "Number of points must be a positive integer.");
        points_ = val;
        calculate_coefficients();
        return *this;
    }
 
private:
    [[nodiscard]] static auto diff_coeff(
        variable_type changed_var, variable_type half_width) -> variable_type {
        const variable_type exp_value =
            std::exp(-constants::pi<variable_type> * std::sinh(changed_var));
        const variable_type exp_value_p1 =
            constants::one<variable_type> + exp_value;
        return constants::two<variable_type> * constants::pi<variable_type> *
            half_width * std::cosh(changed_var) * exp_value /
            (exp_value_p1 * exp_value_p1);
    }

    void calculate_coefficients() {
        using constants::one;
        using constants::pi;
 
        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(-pi<variable_type> * std::sinh(changed_variable));
            const variable_type denominator = one<variable_type> + exp_value;
            variable_rate_list_.push_back(exp_value / denominator);
            weight_rate_list_.push_back(pi<variable_type> *
                std::cosh(changed_variable) * exp_value /
                (denominator * denominator));
        }
    }

    static constexpr auto default_max_point = static_cast<variable_type>(4);

    variable_type max_point_{default_max_point};

    static constexpr index_type default_points = 20;

    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