inexact_newton_update_equation_solver.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 <optional>
 
#include <Eigen/Core>
#include <Eigen/LU>
 
#include "num_collect/base/exception.h"
#include "num_collect/base/index_type.h"
#include "num_collect/base/iterative_solver_base.h"
#include "num_collect/base/precondition.h"
#include "num_collect/logging/iterations/iteration_logger.h"
#include "num_collect/logging/log_tag_view.h"
#include "num_collect/logging/logging_macros.h"
#include "num_collect/ode/concepts/differentiable_problem.h"
#include "num_collect/ode/concepts/multi_variate_differentiable_problem.h"
#include "num_collect/ode/concepts/single_variate_differentiable_problem.h"
#include "num_collect/ode/error_tolerances.h"
#include "num_collect/ode/evaluation_type.h"
 
namespace num_collect::ode {
 
constexpr auto inexact_newton_update_equation_solver_tag =
    logging::log_tag_view(
        "num_collect::ode::inexact_newton_update_equation_solver");
 
template <concepts::differentiable_problem Problem>
class inexact_newton_update_equation_solver;
 
template <concepts::single_variate_differentiable_problem Problem>
class inexact_newton_update_equation_solver<Problem>
    : public iterative_solver_base<
          inexact_newton_update_equation_solver<Problem>> {
public:
    using this_type = inexact_newton_update_equation_solver<Problem>;
 
    using problem_type = Problem;
 
    using variable_type = typename problem_type::variable_type;
 
    using scalar_type = typename problem_type::scalar_type;
 
    using jacobian_type = typename problem_type::jacobian_type;
 
    static_assert(!problem_type::allowed_evaluations.mass,
        "Mass matrix is not supported.");
    // TODO: Actually, support is not difficult, but I want to test after
    // implementation, so I postpone the implementation.
 
    inexact_newton_update_equation_solver()
        : iterative_solver_base<inexact_newton_update_equation_solver<Problem>>(
              inexact_newton_update_equation_solver_tag) {}
 
    void update_jacobian(problem_type& problem, scalar_type time,
        scalar_type step_size, const variable_type& variable,
        scalar_type slope_coeff) {
        problem_ = &problem;
        time_ = time;
        step_size_ = step_size;
        variable_ = variable;
        slope_coeff_ = slope_coeff;
 
        problem_->evaluate_on(time_, variable_,
            evaluation_type{.diff_coeff = true, .jacobian = true});
 
        coeff_inverse_ = static_cast<jacobian_type>(1) /
            (static_cast<jacobian_type>(1) -
                step_size_ * slope_coeff_ * problem_->jacobian());
        using std::isfinite;
        if (!isfinite(coeff_inverse_)) {
            NUM_COLLECT_LOG_AND_THROW(
                algorithm_failure, "Failed to calculate inverse.");
        }
    }
 
    void init(const variable_type& solution_offset, variable_type& solution) {
        solution_offset_ = solution_offset;
        solution_ = &solution;
        update_norm_.reset();
        if (update_reduction_rate_) {
            constexpr auto exponent = static_cast<scalar_type>(0.8);
            constexpr auto min_rate = static_cast<scalar_type>(0.5);
            using std::pow;
            *update_reduction_rate_ = pow(*update_reduction_rate_, exponent);
            if (*update_reduction_rate_ < min_rate) {
                *update_reduction_rate_ = min_rate;
            }
        }
        iterations_ = 0;
    }
 
    void init(const scalar_type& time, const variable_type& solution_offset,
        variable_type& solution) {
        time_ = time;
        init(solution_offset, solution);
    }
 
    void iterate() {
        NUM_COLLECT_PRECONDITION(problem_ != nullptr && solution_ != nullptr,
            this->logger(), "Initialization must be done before iterations.");
 
        temp_variable_ = variable_ + (*solution_);
 
        problem_->evaluate_on(
            time_, temp_variable_, evaluation_type{.diff_coeff = true});
        residual_ = (*solution_) -
            step_size_ * slope_coeff_ * problem_->diff_coeff() -
            solution_offset_;
        update_ = -coeff_inverse_ * residual_;
        *solution_ += update_;
 
        const scalar_type update_norm =
            tolerances().calc_norm(variable_, update_);
        if (update_norm_) {
            update_reduction_rate_ = update_norm / (*update_norm_);
        }
        update_norm_ = update_norm;
 
        ++iterations_;
    }
 
    [[nodiscard]] auto is_stop_criteria_satisfied() const -> bool {
        bool converged = false;
        if (update_norm_ && update_reduction_rate_ &&
            *update_reduction_rate_ < static_cast<scalar_type>(1)) {
            converged =
                (*update_reduction_rate_ /
                    (static_cast<scalar_type>(1) - *update_reduction_rate_)) *
                    (*update_norm_) <=
                tolerance_rate_;
        }
 
        constexpr index_type max_iterations = 1000;  // safe guard
        return converged || (iterations_ > max_iterations);
    }
 
    void configure_iteration_logger(
        logging::iterations::iteration_logger<this_type>& iteration_logger)
        const {
        iteration_logger.template append<index_type>(
            "Iter.", &this_type::iterations);
        iteration_logger.template append<scalar_type>(
            "Update", &this_type::update_norm);
    }
 
    [[nodiscard]] auto solution_offset() const -> const variable_type& {
        return solution_offset_;
    }
 
    [[nodiscard]] auto update_norm() const -> scalar_type {
        if (!update_norm_) {
            return static_cast<scalar_type>(0);
        }
        return *update_norm_;
    }
 
    [[nodiscard]] auto iterations() const -> index_type { return iterations_; }
 
    auto tolerances(const error_tolerances<variable_type>& val)
        -> inexact_newton_update_equation_solver& {
        tolerances_ = val;
        return *this;
    }
 
    [[nodiscard]] auto tolerances() const
        -> const error_tolerances<variable_type>& {
        return tolerances_;
    }
 
private:
    problem_type* problem_{nullptr};
 
    scalar_type time_{};
 
    scalar_type step_size_{};
 
    scalar_type slope_coeff_{};
 
    variable_type variable_{};
 
    variable_type solution_offset_{};
 
    variable_type* solution_{nullptr};
 
    variable_type coeff_inverse_{};
 
    variable_type temp_variable_{};
 
    variable_type residual_{};
 
    variable_type update_{};
 
    std::optional<scalar_type> update_norm_{};
 
    std::optional<scalar_type> update_reduction_rate_{};
 
    static constexpr auto default_tolerance_rate =
        static_cast<scalar_type>(1e-2);
 
    scalar_type tolerance_rate_{default_tolerance_rate};
 
    index_type iterations_{0};
 
    error_tolerances<variable_type> tolerances_{};
};
 
template <concepts::multi_variate_differentiable_problem Problem>
class inexact_newton_update_equation_solver<Problem>
    : public iterative_solver_base<
          inexact_newton_update_equation_solver<Problem>> {
public:
    using this_type = inexact_newton_update_equation_solver<Problem>;
 
    using problem_type = Problem;
 
    using variable_type = typename problem_type::variable_type;
 
    using scalar_type = typename problem_type::scalar_type;
 
    using jacobian_type = typename problem_type::jacobian_type;
 
    static_assert(!problem_type::allowed_evaluations.mass,
        "Mass matrix is not supported.");
    // TODO: Actually, support is not difficult, but I want to test after
    // implementation, so I postpone the implementation.
 
    inexact_newton_update_equation_solver()
        : iterative_solver_base<inexact_newton_update_equation_solver<Problem>>(
              inexact_newton_update_equation_solver_tag) {}
 
    template <typename VariableExpression>
    void update_jacobian(problem_type& problem, scalar_type time,
        scalar_type step_size,
        const Eigen::MatrixBase<VariableExpression>& variable,
        scalar_type slope_coeff) {
        problem_ = &problem;
        time_ = time;
        step_size_ = step_size;
        variable_ = variable;
        slope_coeff_ = slope_coeff;
 
        problem_->evaluate_on(time_, variable_,
            evaluation_type{.diff_coeff = true, .jacobian = true});
 
        const index_type dim = variable_.size();
        lu_.compute(jacobian_type::Identity(dim, dim) -
            step_size_ * slope_coeff_ * problem_->jacobian());
    }
 
    template <typename OffsetExpression>
    void init(const Eigen::MatrixBase<OffsetExpression>& solution_offset,
        variable_type& solution) {
        solution_offset_ = solution_offset;
        solution_ = &solution;
        update_norm_.reset();
        if (update_reduction_rate_) {
            constexpr auto exponent = static_cast<scalar_type>(0.8);
            constexpr auto min_rate = static_cast<scalar_type>(0.5);
            using std::pow;
            *update_reduction_rate_ = pow(*update_reduction_rate_, exponent);
            if (*update_reduction_rate_ < min_rate) {
                *update_reduction_rate_ = min_rate;
            }
        }
        iterations_ = 0;
    }
 
    template <typename OffsetExpression>
    void init(const scalar_type& time,
        const Eigen::MatrixBase<OffsetExpression>& solution_offset,
        variable_type& solution) {
        time_ = time;
        init(solution_offset, solution);
    }
 
    void iterate() {
        NUM_COLLECT_PRECONDITION(problem_ != nullptr && solution_ != nullptr,
            this->logger(), "Initialization must be done before iterations.");
 
        temp_variable_ = variable_ + (*solution_);
 
        problem_->evaluate_on(
            time_, temp_variable_, evaluation_type{.diff_coeff = true});
        residual_ = (*solution_) -
            step_size_ * slope_coeff_ * problem_->diff_coeff() -
            solution_offset_;
        update_ = -lu_.solve(residual_);
        if (!update_.array().isFinite().all()) {
            NUM_COLLECT_LOG_AND_THROW(algorithm_failure,
                "Failed to solve an equation. step_size={}, cond={}.",
                step_size_, lu_.rcond());
        }
        *solution_ += update_;
 
        const scalar_type update_norm =
            tolerances().calc_norm(variable_, update_);
        if (update_norm_) {
            update_reduction_rate_ = update_norm / (*update_norm_);
        }
        update_norm_ = update_norm;
 
        ++iterations_;
    }
 
    [[nodiscard]] auto is_stop_criteria_satisfied() const -> bool {
        bool converged = false;
        if (update_norm_ && update_reduction_rate_ &&
            *update_reduction_rate_ < static_cast<scalar_type>(1)) {
            converged =
                (*update_reduction_rate_ /
                    (static_cast<scalar_type>(1) - *update_reduction_rate_)) *
                    (*update_norm_) <=
                tolerance_rate_;
        }
 
        constexpr index_type max_iterations = 100;  // safe guard
        return converged || (iterations_ > max_iterations);
    }
 
    void configure_iteration_logger(logging::iterations::iteration_logger<
        inexact_newton_update_equation_solver<Problem>>& iteration_logger)
        const {
        iteration_logger.template append<index_type>(
            "Iter.", &this_type::iterations);
        iteration_logger.template append<scalar_type>(
            "Update", &this_type::update_norm);
    }
 
    [[nodiscard]] auto solution_offset() const -> const variable_type& {
        return solution_offset_;
    }
 
    [[nodiscard]] auto update_norm() const -> scalar_type {
        if (!update_norm_) {
            return static_cast<scalar_type>(0);
        }
        return *update_norm_;
    }
 
    [[nodiscard]] auto iterations() const -> index_type { return iterations_; }
 
    auto tolerances(const error_tolerances<variable_type>& val)
        -> inexact_newton_update_equation_solver& {
        tolerances_ = val;
        return *this;
    }
 
    [[nodiscard]] auto tolerances() const
        -> const error_tolerances<variable_type>& {
        return tolerances_;
    }
 
private:
    problem_type* problem_{nullptr};
 
    scalar_type time_{};
 
    scalar_type step_size_{};
 
    scalar_type slope_coeff_{};
 
    variable_type variable_{};
 
    variable_type solution_offset_{};
 
    variable_type* solution_{nullptr};
 
    Eigen::PartialPivLU<jacobian_type> lu_{};
 
    variable_type temp_variable_{};
 
    variable_type residual_{};
 
    variable_type update_{};
 
    std::optional<scalar_type> update_norm_{};
 
    std::optional<scalar_type> update_reduction_rate_{};
 
    static constexpr auto default_tolerance_rate =
        static_cast<scalar_type>(1e-2);
 
    scalar_type tolerance_rate_{default_tolerance_rate};
 
    index_type iterations_{0};
 
    error_tolerances<variable_type> tolerances_{};
};
 
}  // namespace num_collect::ode