error_tolerances.h

Go to the documentation of this file.

/*
 * Copyright 2022 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 <variant>
 
#include "num_collect/base/concepts/real_scalar.h"
#include "num_collect/base/concepts/real_scalar_dense_vector.h"
#include "num_collect/base/exception.h"
#include "num_collect/base/precondition.h"
#include "num_collect/logging/logging_macros.h"
 
namespace num_collect::ode {
 
namespace impl {

template <base::concepts::real_scalar Scalar>
constexpr Scalar default_tol_rel_error = static_cast<Scalar>(1e-6);

template <base::concepts::real_scalar Scalar>
constexpr Scalar default_tol_abs_error = static_cast<Scalar>(1e-6);
 
}  // namespace impl

template <typename Variable>
class error_tolerances;

template <base::concepts::real_scalar_dense_vector Variable>
class error_tolerances<Variable> {
public:
    using variable_type = Variable;

    using scalar_type = typename variable_type::Scalar;

    error_tolerances()
        : tol_rel_error_(impl::default_tol_rel_error<scalar_type>),
          tol_abs_error_(impl::default_tol_abs_error<scalar_type>) {}

    [[nodiscard]] auto check(const variable_type& variable,
        const variable_type& error) const -> bool {
        class evaluator {
        public:
            evaluator(const variable_type& variable, const variable_type& error)
                : variable_(variable), error_(error) {}
 
            [[nodiscard]] auto operator()(const scalar_type& tol_rel_error,
                const scalar_type& tol_abs_error) const -> bool {
                return (error_.array().abs() <=
                    (tol_rel_error * variable_.array().abs() + tol_abs_error))
                    .all();
            }
 
            [[nodiscard]] auto operator()(const variable_type& tol_rel_error,
                const scalar_type& tol_abs_error) const -> bool {
                return (error_.array().abs() <=
                    (tol_rel_error.array() * variable_.array().abs() +
                        tol_abs_error))
                    .all();
            }
 
            [[nodiscard]] auto operator()(const scalar_type& tol_rel_error,
                const variable_type& tol_abs_error) const -> bool {
                return (error_.array().abs() <=
                    (tol_rel_error * variable_.array().abs() +
                        tol_abs_error.array()))
                    .all();
            }
 
            [[nodiscard]] auto operator()(const variable_type& tol_rel_error,
                const variable_type& tol_abs_error) const -> bool {
                return (error_.array().abs() <=
                    (tol_rel_error.array() * variable_.array().abs() +
                        tol_abs_error.array()))
                    .all();
            }
 
        private:
            const variable_type& variable_;
            const variable_type& error_;
        };
 
        evaluator eval{variable, error};
        return std::visit(eval, tol_rel_error_, tol_abs_error_);
    }

    [[nodiscard]] auto calc_norm(const variable_type& variable,
        const variable_type& error) const -> scalar_type {
        class evaluator {
        public:
            evaluator(const variable_type& variable, const variable_type& error)
                : variable_(variable), error_(error) {}
 
            [[nodiscard]] auto operator()(const scalar_type& tol_rel_error,
                const scalar_type& tol_abs_error) const -> scalar_type {
                using std::sqrt;
                return sqrt((error_.array() /
                                (tol_rel_error * variable_.array().abs() +
                                    tol_abs_error))
                                .abs2()
                                .sum() /
                    static_cast<scalar_type>(variable_.size()));
            }
 
            [[nodiscard]] auto operator()(const variable_type& tol_rel_error,
                const scalar_type& tol_abs_error) const -> scalar_type {
                using std::sqrt;
                return sqrt(
                    (error_.array() /
                        (tol_rel_error.array() * variable_.array().abs() +
                            tol_abs_error))
                        .abs2()
                        .sum() /
                    static_cast<scalar_type>(variable_.size()));
            }
 
            [[nodiscard]] auto operator()(const scalar_type& tol_rel_error,
                const variable_type& tol_abs_error) const -> scalar_type {
                using std::sqrt;
                return sqrt((error_.array() /
                                (tol_rel_error * variable_.array().abs() +
                                    tol_abs_error.array()))
                                .abs2()
                                .sum() /
                    static_cast<scalar_type>(variable_.size()));
            }
 
            [[nodiscard]] auto operator()(const variable_type& tol_rel_error,
                const variable_type& tol_abs_error) const -> scalar_type {
                using std::sqrt;
                return sqrt(
                    (error_.array() /
                        (tol_rel_error.array() * variable_.array().abs() +
                            tol_abs_error.array()))
                        .abs2()
                        .sum() /
                    static_cast<scalar_type>(variable_.size()));
            }
 
        private:
            const variable_type& variable_;
            const variable_type& error_;
        };
 
        evaluator eval{variable, error};
        return std::visit(eval, tol_rel_error_, tol_abs_error_);
    }

    auto tol_rel_error(const scalar_type& val) -> error_tolerances& {
        NUM_COLLECT_PRECONDITION(val >= static_cast<scalar_type>(0),
            "Tolerance of relative error must be a non-negative value.");
        tol_rel_error_ = val;
        return *this;
    }

    auto tol_rel_error(const variable_type& val) -> error_tolerances& {
        NUM_COLLECT_PRECONDITION(
            (val.array() >= static_cast<scalar_type>(0)).all(),
            "Tolerance of relative error must be non-negative values.");
        tol_rel_error_ = val;
        return *this;
    }

    auto tol_abs_error(const scalar_type& val) -> error_tolerances& {
        NUM_COLLECT_PRECONDITION(val >= static_cast<scalar_type>(0),
            "Tolerance of absolute error must be a non-negative value.");
        tol_abs_error_ = val;
        return *this;
    }

    auto tol_abs_error(const variable_type& val) -> error_tolerances& {
        NUM_COLLECT_PRECONDITION(
            (val.array() >= static_cast<scalar_type>(0)).all(),
            "Tolerance of absolute error must be non-negative values.");
        tol_abs_error_ = val;
        return *this;
    }
 
private:
    using variant_type = std::variant<scalar_type, variable_type>;

    variant_type tol_rel_error_;

    variant_type tol_abs_error_;
};

template <base::concepts::real_scalar Variable>
class error_tolerances<Variable> {
public:
    using variable_type = Variable;

    using scalar_type = variable_type;

    explicit error_tolerances()
        : tol_rel_error_(impl::default_tol_rel_error<scalar_type>),
          tol_abs_error_(impl::default_tol_abs_error<scalar_type>) {}

    [[nodiscard]] auto check(const variable_type& variable,
        const variable_type& error) const -> bool {
        using std::abs;
        return abs(error) <= tol_rel_error_ * abs(variable) + tol_abs_error_;
    }

    [[nodiscard]] auto calc_norm(const variable_type& variable,
        const variable_type& error) const -> scalar_type {
        using std::abs;
        return abs(error / (tol_rel_error_ * abs(variable) + tol_abs_error_));
    }

    auto tol_rel_error(const variable_type& val) -> error_tolerances& {
        NUM_COLLECT_PRECONDITION(val >= static_cast<scalar_type>(0),
            "Tolerance of relative error must be a non-negative value.");
        tol_rel_error_ = val;
        return *this;
    }

    auto tol_abs_error(const variable_type& val) -> error_tolerances& {
        NUM_COLLECT_PRECONDITION(val >= static_cast<scalar_type>(0),
            "Tolerance of absolute error must be a non-negative value.");
        tol_abs_error_ = val;
        return *this;
    }
 
private:
    variable_type tol_rel_error_;

    variable_type tol_abs_error_;
};
 
}  // namespace num_collect::ode