de_semi_infinite_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>
#include <type_traits>
 
#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/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_mixin.h"
 
namespace num_collect::integration {

constexpr auto de_semi_infinite_integrator_tag = logging::log_tag_view(
    "num_collect::integration::de_semi_infinite_integrator");
 
namespace impl {

template <typename Variable>
struct de_semi_infinite_integrator_traits;

template <>
struct de_semi_infinite_integrator_traits<float> {
public:
    static constexpr float default_max_point = 3.0F;
};

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

template <typename Signature>
class de_semi_infinite_integrator;

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

    using result_type = std::decay_t<Result>;

    de_semi_infinite_integrator()
        : logging::logging_mixin(de_semi_infinite_integrator_tag) {}

    template <base::concepts::invocable_as<result_type(variable_type)> Function>
    [[nodiscard]] auto integrate(const Function& function) const
        -> result_type {
        const variable_type interval =
            max_point_ / static_cast<variable_type>(points_);
 
        constexpr variable_type diff_coeff_center =
            constants::pi<variable_type>;
        result_type sum =
            function(constants::one<variable_type>) * diff_coeff_center;
 
        for (index_type i = 1; i < points_; ++i) {
            const variable_type changed_var =
                interval * static_cast<variable_type>(i);
            const variable_type pi_sinh_value =
                constants::pi<variable_type> * std::sinh(changed_var);
 
            const variable_type var_plus = std::exp(pi_sinh_value);
            const variable_type diff_coeff_plus = diff_coeff(changed_var);
            sum += function(var_plus) * diff_coeff_plus;
 
            const variable_type var_minus =
                constants::one<variable_type> / var_plus;
            const variable_type diff_coeff_minus = diff_coeff(-changed_var);
            sum += function(var_minus) * diff_coeff_minus;
        }
 
        return sum * interval;
    }

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

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

    auto points(index_type val) -> de_semi_infinite_integrator& {
        NUM_COLLECT_PRECONDITION(val > 0, this->logger(),
            "Number of points must be a positive integer.");
        points_ = val;
        return *this;
    }
 
private:
    [[nodiscard]] static auto diff_coeff(variable_type changed_var)
        -> variable_type {
        return constants::pi<variable_type> *
            std::exp(constants::pi<variable_type> * std::sinh(changed_var)) *
            std::cosh(changed_var);
    }

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

    variable_type max_point_{default_max_point};

    static constexpr index_type default_points = 20;

    index_type points_{default_points};
};
 
}  // namespace num_collect::integration