From d65d1b1dbfc06bfe2aa038349f382b80ca4082e3 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Fri, 9 Feb 2024 13:31:14 -0500 Subject: [PATCH] Drop cppoptlib (#66) * Drop cppoptlib from Problem and Solver * Cleanup changes * Private CMake target_link_libraries * Use TVector * Remove cppoptlib line-search methods * Add custom statuses * Move extra solver logic to Criteria * Add our cases to operator<< * Update Criteria::print * Removed nan setting * Remove unused line search names * Fix not descent direction failsafe * Add codecov.yml --- CMakeLists.txt | 8 +- cmake/recipes/cppoptlib.cmake | 19 -- codecov.yml | 11 + nonlinear-solver-spec.json | 5 +- src/polysolve/nonlinear/CMakeLists.txt | 12 +- src/polysolve/nonlinear/Criteria.cpp | 112 ++++++++ src/polysolve/nonlinear/Criteria.hpp | 52 ++++ src/polysolve/nonlinear/Problem.hpp | 92 ++++++- src/polysolve/nonlinear/Solver.cpp | 241 ++++++++---------- src/polysolve/nonlinear/Solver.hpp | 85 +++--- .../nonlinear/descent_strategies/Newton.cpp | 3 +- .../nonlinear/line_search/Backtracking.cpp | 2 +- .../nonlinear/line_search/CMakeLists.txt | 4 - .../nonlinear/line_search/CppOptArmijo.cpp | 26 -- .../nonlinear/line_search/CppOptArmijo.hpp | 32 --- .../nonlinear/line_search/LineSearch.cpp | 38 +-- .../nonlinear/line_search/MoreThuente.cpp | 25 -- .../nonlinear/line_search/MoreThuente.hpp | 32 --- 18 files changed, 421 insertions(+), 378 deletions(-) delete mode 100644 cmake/recipes/cppoptlib.cmake create mode 100644 codecov.yml create mode 100644 src/polysolve/nonlinear/Criteria.cpp create mode 100644 src/polysolve/nonlinear/Criteria.hpp delete mode 100644 src/polysolve/nonlinear/line_search/CppOptArmijo.cpp delete mode 100644 src/polysolve/nonlinear/line_search/CppOptArmijo.hpp delete mode 100644 src/polysolve/nonlinear/line_search/MoreThuente.cpp delete mode 100644 src/polysolve/nonlinear/line_search/MoreThuente.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 74e16933..474ef9ba 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -311,17 +311,13 @@ target_link_libraries(polysolve_linear PRIVATE polysolve::warnings) # polysolve::linear target_link_libraries(polysolve PUBLIC polysolve::linear) -# CppNumericalSolvers -include(cppoptlib) -target_link_libraries(polysolve PUBLIC cppoptlib) - # LBFGSpp include(LBFGSpp) -target_link_libraries(polysolve PUBLIC LBFGSpp::LBFGSpp) +target_link_libraries(polysolve PRIVATE LBFGSpp::LBFGSpp) # finite-diff (include this after eigen) include(finite-diff) -target_link_libraries(polysolve PUBLIC finitediff::finitediff) +target_link_libraries(polysolve PRIVATE finitediff::finitediff) # Sanitizers if(POLYSOLVE_WITH_SANITIZERS) diff --git a/cmake/recipes/cppoptlib.cmake b/cmake/recipes/cppoptlib.cmake deleted file mode 100644 index e19709e5..00000000 --- a/cmake/recipes/cppoptlib.cmake +++ /dev/null @@ -1,19 +0,0 @@ -# CppNumericalSolvers (https://github.com/PatWie/CppNumericalSolvers.git) -# License: MIT - -if(TARGET cppoptlib) - return() -endif() - -message(STATUS "Third-party: creating target 'cppoptlib'") - -include(CPM) -CPMAddPackage( - NAME cppoptlib - GITHUB_REPOSITORY PatWie/CppNumericalSolvers - GIT_TAG 7eddf28fa5a8872a956d3c8666055cac2f5a535d - DOWNLOAD_ONLY TRUE -) - -add_library(cppoptlib INTERFACE) -target_include_directories(cppoptlib SYSTEM INTERFACE ${cppoptlib_SOURCE_DIR}/include) \ No newline at end of file diff --git a/codecov.yml b/codecov.yml new file mode 100644 index 00000000..3418f88c --- /dev/null +++ b/codecov.yml @@ -0,0 +1,11 @@ +coverage: + status: + project: + default: + target: auto + threshold: 1% + patch: + default: + target: 75% + threshold: 5% + only_pulls: true diff --git a/nonlinear-solver-spec.json b/nonlinear-solver-spec.json index e615c73e..fc1ff34b 100644 --- a/nonlinear-solver-spec.json +++ b/nonlinear-solver-spec.json @@ -1,5 +1,4 @@ -[ - { +[{ "pointer": "/", "default": null, "type": "object", @@ -562,10 +561,8 @@ "type": "string", "options": [ "Armijo", - "ArmijoAlt", "RobustArmijo", "Backtracking", - "MoreThuente", "None" ], "doc": "Line-search type" diff --git a/src/polysolve/nonlinear/CMakeLists.txt b/src/polysolve/nonlinear/CMakeLists.txt index c842dbe3..dfbb760c 100644 --- a/src/polysolve/nonlinear/CMakeLists.txt +++ b/src/polysolve/nonlinear/CMakeLists.txt @@ -1,12 +1,14 @@ set(SOURCES BoxConstraintSolver.cpp BoxConstraintSolver.hpp - Solver.hpp - Solver.cpp - Problem.hpp - Problem.cpp - PostStepData.hpp + Criteria.cpp + Criteria.hpp PostStepData.cpp + PostStepData.hpp + Problem.cpp + Problem.hpp + Solver.cpp + Solver.hpp ) source_group(TREE "${CMAKE_CURRENT_SOURCE_DIR}" PREFIX "Source Files" FILES ${SOURCES}) diff --git a/src/polysolve/nonlinear/Criteria.cpp b/src/polysolve/nonlinear/Criteria.cpp new file mode 100644 index 00000000..37fed882 --- /dev/null +++ b/src/polysolve/nonlinear/Criteria.cpp @@ -0,0 +1,112 @@ +// Source: https://github.com/PatWie/CppNumericalSolvers/blob/7eddf28fa5a8872a956d3c8666055cac2f5a535d/include/cppoptlib/meta.h +// License: MIT +#include "Criteria.hpp" + +namespace polysolve::nonlinear +{ + bool is_converged_status(const Status s) + { + return s == Status::XDeltaTolerance || s == Status::FDeltaTolerance || s == Status::GradNormTolerance; + } + + Criteria::Criteria() + { + reset(); + } + + void Criteria::reset() + { + iterations = 0; + xDelta = 0; + fDelta = 0; + gradNorm = 0; + firstGradNorm = 0; + fDeltaCount = 0; + xDeltaDotGrad = 0; + } + + void Criteria::print(std::ostream &os) const + { + os << "Iterations: " << iterations << std::endl; + os << "xDelta: " << xDelta << std::endl; + os << "fDelta: " << fDelta << std::endl; + os << "GradNorm: " << gradNorm << std::endl; + os << "xDeltaDotGrad: " << xDeltaDotGrad << std::endl; + os << "FirstGradNorm: " << firstGradNorm << std::endl; + os << "fDeltaCount: " << fDeltaCount << std::endl; + } + + Status checkConvergence(const Criteria &stop, const Criteria ¤t) + { + if (stop.iterations > 0 && current.iterations > stop.iterations) + { + return Status::IterationLimit; + } + if (stop.xDelta > 0 && current.xDelta < stop.xDelta) + { + return Status::XDeltaTolerance; + } + if (stop.fDelta > 0 && current.fDelta < stop.fDelta && current.fDeltaCount >= stop.fDeltaCount) + { + return Status::FDeltaTolerance; + } + const double stopGradNorm = current.iterations == 0 ? stop.firstGradNorm : stop.gradNorm; + if (stopGradNorm > 0 && current.gradNorm < stopGradNorm) + { + return Status::GradNormTolerance; + } + // Δx⋅∇f ≥ 0 means that the search direction is not a descent direction + if (stop.xDeltaDotGrad < 0 && current.xDeltaDotGrad > stop.xDeltaDotGrad) + { + return Status::NotDescentDirection; + } + return Status::Continue; + } + + std::ostream &operator<<(std::ostream &os, const Status &s) + { + switch (s) + { + case Status::NotStarted: + os << "Solver not started."; + break; + case Status::Continue: + os << "Convergence criteria not reached."; + break; + case Status::IterationLimit: + os << "Iteration limit reached."; + break; + case Status::XDeltaTolerance: + os << "Change in parameter vector too small."; + break; + case Status::FDeltaTolerance: + os << "Change in cost function value too small."; + break; + case Status::GradNormTolerance: + os << "Gradient vector norm too small."; + break; + case Status::ObjectiveCustomStop: + os << "Objective function specified to stop."; + break; + case Status::NanEncountered: + os << "Objective or gradient function returned NaN."; + break; + case Status::NotDescentDirection: + os << "Search direction not a descent direction."; + break; + case Status::LineSearchFailed: + os << "Line search failed."; + break; + default: + os << "Unknown status."; + break; + } + return os; + } + + std::ostream &operator<<(std::ostream &os, const Criteria &c) + { + c.print(os); + return os; + } +} // namespace polysolve::nonlinear \ No newline at end of file diff --git a/src/polysolve/nonlinear/Criteria.hpp b/src/polysolve/nonlinear/Criteria.hpp new file mode 100644 index 00000000..f318ed5c --- /dev/null +++ b/src/polysolve/nonlinear/Criteria.hpp @@ -0,0 +1,52 @@ +#pragma once + +#include +#include + +namespace polysolve::nonlinear +{ + // Source: https://github.com/PatWie/CppNumericalSolvers/blob/7eddf28fa5a8872a956d3c8666055cac2f5a535d/include/cppoptlib/meta.h + // License: MIT + + enum class Status + { + NotStarted = -1, ///< The solver has not been started + Continue = 0, ///< The solver should continue + // Success cases + IterationLimit, ///< The maximum number of iterations has been reached + XDeltaTolerance, ///< The change in the parameter vector is below the tolerance + FDeltaTolerance, ///< The change in the cost function is below the tolerance + GradNormTolerance, ///< The norm of the gradient vector is below the tolerance + ObjectiveCustomStop, ///< The objective function specified to stop + // Failure cases + NanEncountered, ///< The objective function returned NaN + NotDescentDirection, ///< The search direction is not a descent direction + LineSearchFailed, ///< The line search failed + }; + + bool is_converged_status(const Status s); + + class Criteria + { + public: + size_t iterations; ///< Maximum number of iterations + double xDelta; ///< Minimum change in parameter vector + double fDelta; ///< Minimum change in cost function + double gradNorm; ///< Minimum norm of gradient vector + double firstGradNorm; ///< Initial norm of gradient vector + double xDeltaDotGrad; ///< Dot product of parameter vector and gradient vector + unsigned fDeltaCount; ///< Number of steps where fDelta is satisfied + + Criteria(); + + void reset(); + + void print(std::ostream &os) const; + }; + + Status checkConvergence(const Criteria &stop, const Criteria ¤t); + + std::ostream &operator<<(std::ostream &os, const Status &s); + + std::ostream &operator<<(std::ostream &os, const Criteria &c); +} // namespace polysolve::nonlinear \ No newline at end of file diff --git a/src/polysolve/nonlinear/Problem.hpp b/src/polysolve/nonlinear/Problem.hpp index 93fa0db6..8c501130 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -2,47 +2,115 @@ #include +#include "Criteria.hpp" #include "PostStepData.hpp" -#include - #include #include namespace polysolve::nonlinear { - class Problem : public cppoptlib::Problem + class Problem { public: - using typename cppoptlib::Problem::Scalar; - using typename cppoptlib::Problem::TVector; - typedef polysolve::StiffnessMatrix THessian; - - // disable warning for dense hessian - using cppoptlib::Problem::hessian; + static constexpr int Dim = Eigen::Dynamic; + using Scalar = double; + using TVector = Eigen::Matrix; + using TMatrix = Eigen::Matrix; + using THessian = StiffnessMatrix; + public: Problem() {} - ~Problem() = default; + virtual ~Problem() = default; + /// @brief Initialize the problem. + /// @param x0 Initial guess. virtual void init(const TVector &x0) {} - virtual double value(const TVector &x) override = 0; - virtual void gradient(const TVector &x, TVector &gradv) override = 0; + /// @brief Compute the value of the function at x. + /// @param x Degrees of freedom. + /// @return The value of the function at x. + Scalar operator()(const TVector &x) { return value(x); } + + /// @brief Compute the value of the function at x. + /// @param x Degrees of freedom. + /// @return The value of the function at x. + virtual Scalar value(const TVector &x) = 0; + + /// @brief Compute the gradient of the function at x. + /// @param[in] x Degrees of freedom. + /// @param[out] grad Gradient of the function at x. + virtual void gradient(const TVector &x, TVector &grad) = 0; + + /// @brief Compute the Hessian of the function at x. + /// @param[in] x Degrees of freedom. + /// @param[out] hessian Hessian of the function at x. + virtual void hessian(const TVector &x, TMatrix &hessian) + { + throw std::runtime_error("Dense Hessian not implemented."); + } + + /// @brief Compute the Hessian of the function at x. + /// @param[in] x Degrees of freedom. + /// @param[out] hessian Hessian of the function at x. virtual void hessian(const TVector &x, THessian &hessian) = 0; + /// @brief Determine if the step from x0 to x1 is valid. + /// @param x0 Starting point. + /// @param x1 Ending point. + /// @return True if the step is valid, false otherwise. virtual bool is_step_valid(const TVector &x0, const TVector &x1) { return true; } + + /// @brief Determine a maximum step size from x0 to x1. + /// @param x0 Starting point. + /// @param x1 Ending point. + /// @return Maximum step size. virtual double max_step_size(const TVector &x0, const TVector &x1) { return 1; } + // --- Callbacks ------------------------------------------------------ + + /// @brief Callback function for the start of a line search. + /// @param x0 Starting point. + /// @param x1 Ending point. virtual void line_search_begin(const TVector &x0, const TVector &x1) {} + + /// @brief Callback function for the end of a line search. virtual void line_search_end() {} + + /// @brief Callback function for the end of a step. + /// @param data Post step data. virtual void post_step(const PostStepData &data) {} + /// @brief Set the project to PSD flag. + /// @param val True if the problem should be projected to PSD, false otherwise. virtual void set_project_to_psd(bool val) {} + /// @brief Callback function for when the solution changes. + /// @param new_x New solution. virtual void solution_changed(const TVector &new_x) {} + /// @brief Callback function used to determine if the solver should stop. + /// @param state Current state of the solver. + /// @param x Current solution. + /// @return True if the solver should stop, false otherwise. + virtual bool callback(const Criteria &state, const TVector &x) { return true; } + + /// @brief Callback function used Determine if the solver should stop. + /// @param x Current solution. + /// @return True if the solver should stop, false otherwise. virtual bool stop(const TVector &x) { return false; } + /// --- Misc ---------------------------------------------------------- + + /// @brief Sample the function along a direction. + /// @param[in] x Starting point. + /// @param[in] direction Direction to sample along. + /// @param[in] start Starting step size. + /// @param[in] end Ending step size. + /// @param[in] num_samples Number of samples to take. + /// @param[out] alphas Sampled step sizes. + /// @param[out] fs Sampled function values. + /// @param[out] valid If each sample is valid. void sample_along_direction( const Problem::TVector &x, const Problem::TVector &direction, diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index 61e7942e..976115e6 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -196,27 +196,25 @@ namespace polysolve::nonlinear spdlog::logger &logger) : m_logger(logger), characteristic_length(characteristic_length) { - TCriteria criteria = TCriteria::defaults(); - criteria.xDelta = solver_params["x_delta"]; - criteria.fDelta = solver_params["advanced"]["f_delta"]; - criteria.gradNorm = solver_params["grad_norm"]; - - criteria.xDelta *= characteristic_length; - criteria.fDelta *= characteristic_length; - criteria.gradNorm *= characteristic_length; - - criteria.iterations = solver_params["max_iterations"]; - // criteria.condition = solver_params["condition"]; - this->setStopCriteria(criteria); - - use_grad_norm_tol = solver_params["line_search"]["use_grad_norm_tol"]; - first_grad_norm_tol = solver_params["first_grad_norm_tol"]; + m_current.reset(); + + m_stop.xDelta = solver_params["x_delta"]; + m_stop.fDelta = solver_params["advanced"]["f_delta"]; + m_stop.gradNorm = solver_params["grad_norm"]; + m_stop.firstGradNorm = solver_params["first_grad_norm_tol"]; + m_stop.xDeltaDotGrad = -solver_params["advanced"]["derivative_along_delta_x_tol"].get(); + + // Make these relative to the characteristic length + m_stop.xDelta *= characteristic_length; + m_stop.fDelta *= characteristic_length; + m_stop.gradNorm *= characteristic_length; + m_stop.firstGradNorm *= characteristic_length; + // m_stop.xDeltaDotGrad *= characteristic_length; + + m_stop.iterations = solver_params["max_iterations"]; allow_out_of_iterations = solver_params["allow_out_of_iterations"]; - use_grad_norm_tol *= characteristic_length; - first_grad_norm_tol *= characteristic_length; - - f_delta_step_tol = solver_params["advanced"]["f_delta_step_tol"]; + m_stop.fDeltaCount = solver_params["advanced"]["f_delta_step_tol"]; m_descent_strategy = 0; @@ -224,7 +222,6 @@ namespace polysolve::nonlinear gradient_fd_strategy = solver_params["advanced"]["apply_gradient_fd"]; gradient_fd_eps = solver_params["advanced"]["gradient_fd_eps"]; - derivative_along_delta_x_tol = solver_params["advanced"]["derivative_along_delta_x_tol"]; } void Solver::set_strategies_iterations(const json &solver_params) @@ -242,7 +239,7 @@ namespace polysolve::nonlinear m_iter_per_strategy.assign(m_strategies.size() + 1, solver_params["iterations_per_strategy"].get()); } - double Solver::compute_grad_norm(const Eigen::VectorXd &x, const Eigen::VectorXd &grad) const + double Solver::compute_grad_norm(const TVector &x, const TVector &grad) const { return grad.norm(); } @@ -251,6 +248,8 @@ namespace polysolve::nonlinear { m_line_search = line_search::LineSearch::create(params, m_logger); solver_info["line_search"] = params["line_search"]["method"]; + m_line_search->use_grad_norm_tol = params["line_search"]["use_grad_norm_tol"]; + m_line_search->use_grad_norm_tol *= characteristic_length; } void Solver::minimize(Problem &objFunc, TVector &x) @@ -264,8 +263,6 @@ namespace polysolve::nonlinear // --------------------------- reset(x.size()); // place for children to initialize their fields - m_line_search->use_grad_norm_tol = use_grad_norm_tol; - TVector grad = TVector::Zero(x.rows()); TVector delta_x = TVector::Zero(x.rows()); @@ -276,56 +273,41 @@ namespace polysolve::nonlinear objFunc.solution_changed(x); } - const auto g_norm_tol = this->m_stop.gradNorm; - this->m_stop.gradNorm = first_grad_norm_tol; - - StopWatch stop_watch("nonlinear solver", this->total_time, m_logger); + StopWatch stop_watch("nonlinear solver", total_time, m_logger); stop_watch.start(); m_logger.debug( "Starting {} with {} solve f₀={:g} " "(stopping criteria: max_iters={:d} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g})", descent_strategy_name(), m_line_search->name(), - objFunc.value(x), this->m_stop.iterations, - this->m_stop.fDelta, this->m_stop.gradNorm, this->m_stop.xDelta); + objFunc(x), m_stop.iterations, + m_stop.fDelta, m_stop.gradNorm, m_stop.xDelta); - update_solver_info(objFunc.value(x)); - objFunc.post_step(PostStepData(this->m_current.iterations, solver_info, x, grad)); - - int f_delta_step_cnt = 0; - double f_delta = 0; - - // Used for logging - double xDelta = 0, gradNorm = 0; + update_solver_info(objFunc(x)); + objFunc.post_step(PostStepData(m_current.iterations, solver_info, x, grad)); do { m_line_search->set_is_final_strategy(m_descent_strategy == m_strategies.size() - 1); - this->m_current.xDelta = NaN; - this->m_current.fDelta = NaN; - this->m_current.gradNorm = NaN; + // --- Energy ------------------------------------------------------ - //////////// Energy double energy; { POLYSOLVE_SCOPED_STOPWATCH("compute objective function", obj_fun_time, m_logger); - energy = objFunc.value(x); + energy = objFunc(x); } if (!std::isfinite(energy)) { - this->m_status = cppoptlib::Status::UserDefined; - m_error_code = ErrorCode::NAN_ENCOUNTERED; + m_status = Status::NanEncountered; log_and_throw_error(m_logger, "[{}][{}] f(x) is nan or inf; stopping", descent_strategy_name(), m_line_search->name()); break; } - f_delta = std::abs(old_energy - energy); - // stop based on f_delta only if the solver has taken over f_delta_step_tol steps with small f_delta - this->m_current.fDelta = (f_delta_step_cnt >= f_delta_step_tol) ? f_delta : NaN; + m_current.fDelta = std::abs(old_energy - energy); - ///////////// gradient + // --- Gradient ---------------------------------------------------- { POLYSOLVE_SCOPED_STOPWATCH("compute gradient", grad_time, m_logger); objFunc.gradient(x, grad); @@ -333,57 +315,52 @@ namespace polysolve::nonlinear { POLYSOLVE_SCOPED_STOPWATCH("verify gradient", grad_time, m_logger); - this->verify_gradient(objFunc, x, grad); + verify_gradient(objFunc, x, grad); } - const double grad_norm = compute_grad_norm(x, grad); - if (std::isnan(grad_norm)) + m_current.gradNorm = compute_grad_norm(x, grad); + if (std::isnan(m_current.gradNorm)) { - this->m_status = cppoptlib::Status::UserDefined; - m_error_code = ErrorCode::NAN_ENCOUNTERED; + m_status = Status::NanEncountered; log_and_throw_error(m_logger, "[{}][{}] Gradient is nan; stopping", descent_strategy_name(), m_line_search->name()); - break; } - this->m_current.gradNorm = grad_norm; - gradNorm = this->m_current.gradNorm; - - this->m_status = checkConvergence(this->m_stop, this->m_current); - if (this->m_status != cppoptlib::Status::Continue) + m_status = checkConvergence(m_stop, m_current); + if (m_status != Status::Continue) break; - // ------------------------ - // Compute update direction - // ------------------------ - // Compute a Δx to update the variable - // - bool ok = compute_update_direction(objFunc, x, grad, delta_x); - const double derivative_along_delta_x = delta_x.dot(grad); + // --- Update direction -------------------------------------------- + + const bool ok = compute_update_direction(objFunc, x, grad, delta_x); + m_current.xDeltaDotGrad = delta_x.dot(grad); - if (!ok || std::isnan(grad_norm) || (m_strategies[m_descent_strategy]->is_direction_descent() && grad_norm != 0 && derivative_along_delta_x >= 0)) + if (!ok || (m_strategies[m_descent_strategy]->is_direction_descent() && m_current.gradNorm != 0 && m_current.xDeltaDotGrad >= 0)) { - const auto current_name = descent_strategy_name(); + const std::string current_name = descent_strategy_name(); if (!m_strategies[m_descent_strategy]->handle_error()) ++m_descent_strategy; + if (m_descent_strategy >= m_strategies.size()) { - this->m_status = cppoptlib::Status::UserDefined; - log_and_throw_error(m_logger, "[{}][{}] direction is not a descent direction on last strategy (‖Δx‖={:g}; ‖g‖={:g}; Δx⋅g={:g}≥0); stopping", - current_name, m_line_search->name(), - delta_x.norm(), compute_grad_norm(x, grad), derivative_along_delta_x); + m_status = Status::NotDescentDirection; + log_and_throw_error( + m_logger, "[{}][{}] direction is not a descent direction on last strategy (‖Δx‖={:g}; ‖g‖={:g}; Δx⋅g={:g}≥0); stopping", + current_name, m_line_search->name(), delta_x.norm(), compute_grad_norm(x, grad), m_current.xDeltaDotGrad); + } + else + { + m_status = Status::Continue; + m_logger.debug( + "[{}][{}] direction is not a descent direction (‖Δx‖={:g}; ‖g‖={:g}; Δx⋅g={:g}≥0); reverting to {}", + current_name, m_line_search->name(), delta_x.norm(), compute_grad_norm(x, grad), m_current.xDeltaDotGrad, + descent_strategy_name()); } - - m_logger.debug( - "[{}][{}] direction is not a descent direction (‖Δx‖={:g}; ‖g‖={:g}; Δx⋅g={:g}≥0); reverting to {}", - current_name, m_line_search->name(), - delta_x.norm(), compute_grad_norm(x, grad), derivative_along_delta_x, descent_strategy_name()); - this->m_status = cppoptlib::Status::Continue; continue; } - const double delta_x_norm = delta_x.norm(); - if (std::isnan(delta_x_norm)) + m_current.xDelta = delta_x.norm(); + if (std::isnan(m_current.xDelta)) { const auto current_name = descent_strategy_name(); if (!m_strategies[m_descent_strategy]->handle_error()) @@ -391,46 +368,42 @@ namespace polysolve::nonlinear if (m_descent_strategy >= m_strategies.size()) { - this->m_status = cppoptlib::Status::UserDefined; + m_status = Status::NanEncountered; log_and_throw_error(m_logger, "[{}][{}] Δx is nan on last strategy; stopping", current_name, m_line_search->name()); } - this->m_status = cppoptlib::Status::UserDefined; m_logger.debug("[{}][{}] Δx is nan; reverting to {}", current_name, m_line_search->name(), descent_strategy_name()); - this->m_status = cppoptlib::Status::Continue; + m_status = Status::Continue; continue; } - // Use the maximum absolute displacement value divided by the timestep, - this->m_current.xDelta = delta_x_norm; - xDelta = this->m_current.xDelta; - this->m_status = checkConvergence(this->m_stop, this->m_current); - if (this->m_status == cppoptlib::Status::Continue && derivative_along_delta_x > -derivative_along_delta_x_tol) - this->m_status = cppoptlib::Status::XDeltaTolerance; - if (this->m_status != cppoptlib::Status::Continue) + // --- Check convergence ------------------------------------------- + + m_status = checkConvergence(m_stop, m_current); + + if (m_status != Status::Continue) break; - // --------------- - // Variable update - // --------------- + // --- Variable update --------------------------------------------- m_logger.trace( "[{}][{}] pre LS iter={:d} f={:g} ‖∇f‖={:g}", descent_strategy_name(), m_line_search->name(), - this->m_current.iterations, energy, gradNorm); + m_current.iterations, energy, m_current.gradNorm); // Perform a line_search to compute step scale double rate = m_line_search->line_search(x, delta_x, objFunc); if (std::isnan(rate)) { const auto current_name = descent_strategy_name(); - assert(this->m_status == cppoptlib::Status::Continue); + assert(m_status == Status::Continue); if (!m_strategies[m_descent_strategy]->handle_error()) ++m_descent_strategy; if (m_descent_strategy >= m_strategies.size()) { - this->m_status = cppoptlib::Status::UserDefined; // Line search failed on gradient descent, so quit! + // Line search failed on gradient descent, so quit! + m_status = Status::LineSearchFailed; log_and_throw_error(m_logger, "[{}][{}] Line search failed on last strategy; stopping", current_name, m_line_search->name()); } @@ -470,34 +443,27 @@ namespace polysolve::nonlinear const double step = (rate * delta_x).norm(); update_solver_info(energy); - objFunc.post_step(PostStepData(this->m_current.iterations, solver_info, x, grad)); + objFunc.post_step(PostStepData(m_current.iterations, solver_info, x, grad)); if (objFunc.stop(x)) { - this->m_status = cppoptlib::Status::UserDefined; - m_error_code = ErrorCode::SUCCESS; + m_status = Status::ObjectiveCustomStop; m_logger.debug("[{}][{}] Objective decided to stop", descent_strategy_name(), m_line_search->name()); } - if (f_delta < this->m_stop.fDelta) - f_delta_step_cnt++; - else - f_delta_step_cnt = 0; + m_current.fDeltaCount = (m_current.fDelta < m_stop.fDelta) ? (m_current.fDeltaCount + 1) : 0; m_logger.debug( "[{}][{}] iter={:d} f={:g} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g} Δx⋅∇f(x)={:g} rate={:g} ‖step‖={:g}" " (stopping criteria: max_iters={:d} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g})", descent_strategy_name(), m_line_search->name(), - this->m_current.iterations, energy, f_delta, - gradNorm, xDelta, delta_x.dot(grad), rate, step, - this->m_stop.iterations, this->m_stop.fDelta, this->m_stop.gradNorm, this->m_stop.xDelta); - - if (++this->m_current.iterations >= this->m_stop.iterations) - this->m_status = cppoptlib::Status::IterationLimit; + m_current.iterations, energy, m_current.fDelta, + m_current.gradNorm, m_current.xDelta, delta_x.dot(grad), rate, step, + m_stop.iterations, m_stop.fDelta, m_stop.gradNorm, m_stop.xDelta); - // reset the tolerance, since in the first iter it might be smaller - this->m_stop.gradNorm = g_norm_tol; - } while (objFunc.callback(this->m_current, x) && (this->m_status == cppoptlib::Status::Continue)); + if (++m_current.iterations >= m_stop.iterations) + m_status = Status::IterationLimit; + } while (objFunc.callback(m_current, x) && (m_status == Status::Continue)); stop_watch.stop(); @@ -505,31 +471,31 @@ namespace polysolve::nonlinear // Log results // ----------- - if (!allow_out_of_iterations && this->m_status == cppoptlib::Status::IterationLimit) - log_and_throw_error(m_logger, "[{}][{}] Reached iteration limit (limit={})", descent_strategy_name(), m_line_search->name(), this->m_stop.iterations); - if (this->m_status == cppoptlib::Status::UserDefined && m_error_code != ErrorCode::SUCCESS) + if (!allow_out_of_iterations && m_status == Status::IterationLimit) + log_and_throw_error(m_logger, "[{}][{}] Reached iteration limit (limit={})", descent_strategy_name(), m_line_search->name(), m_stop.iterations); + if (m_status == Status::NanEncountered) log_and_throw_error(m_logger, "[{}][{}] Failed to find minimizer", descent_strategy_name(), m_line_search->name()); double tot_time = stop_watch.getElapsedTimeInSec(); - const bool succeeded = this->m_status == cppoptlib::Status::GradNormTolerance; + const bool succeeded = m_status == Status::GradNormTolerance; m_logger.log( succeeded ? spdlog::level::info : spdlog::level::err, "[{}][{}] Finished: {} Took {:g}s (niters={:d} f={:g} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g})" " (stopping criteria: max_iters={:d} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g})", descent_strategy_name(), m_line_search->name(), - this->m_status, tot_time, this->m_current.iterations, - old_energy, f_delta, gradNorm, xDelta, - this->m_stop.iterations, this->m_stop.fDelta, this->m_stop.gradNorm, this->m_stop.xDelta); + m_status, tot_time, m_current.iterations, + old_energy, m_current.fDelta, m_current.gradNorm, m_current.xDelta, + m_stop.iterations, m_stop.fDelta, m_stop.gradNorm, m_stop.xDelta); log_times(); - update_solver_info(objFunc.value(x)); + update_solver_info(objFunc(x)); } void Solver::reset(const int ndof) { - this->m_current.reset(); + m_current.reset(); m_descent_strategy = 0; - m_error_code = ErrorCode::SUCCESS; + m_status = Status::NotStarted; const std::string line_search_name = solver_info["line_search"]; solver_info = json(); @@ -559,17 +525,14 @@ namespace polysolve::nonlinear void Solver::update_solver_info(const double energy) { - solver_info["status"] = this->status(); - solver_info["error_code"] = m_error_code; + solver_info["status"] = status(); solver_info["energy"] = energy; - const auto &crit = this->criteria(); - solver_info["iterations"] = crit.iterations; - solver_info["xDelta"] = crit.xDelta; - solver_info["fDelta"] = crit.fDelta; - solver_info["gradNorm"] = crit.gradNorm; - solver_info["condition"] = crit.condition; + solver_info["iterations"] = m_current.iterations; + solver_info["xDelta"] = m_current.xDelta; + solver_info["fDelta"] = m_current.fDelta; + solver_info["gradNorm"] = m_current.gradNorm; - double per_iteration = crit.iterations ? crit.iterations : 1; + double per_iteration = m_current.iterations ? m_current.iterations : 1; solver_info["total_time"] = total_time; solver_info["time_grad"] = grad_time / per_iteration; @@ -625,15 +588,15 @@ namespace polysolve::nonlinear return; case FiniteDiffStrategy::DIRECTIONAL_DERIVATIVE: { - Eigen::VectorXd direc = grad.normalized(); - Eigen::VectorXd x2 = x + direc * gradient_fd_eps; - Eigen::VectorXd x1 = x - direc * gradient_fd_eps; + TVector direc = grad.normalized(); + TVector x2 = x + direc * gradient_fd_eps; + TVector x1 = x - direc * gradient_fd_eps; objFunc.solution_changed(x2); - double J2 = objFunc.value(x2); + double J2 = objFunc(x2); objFunc.solution_changed(x1); - double J1 = objFunc.value(x1); + double J1 = objFunc(x1); double fd = (J2 - J1) / 2 / gradient_fd_eps; double analytic = direc.dot(grad); @@ -649,11 +612,11 @@ namespace polysolve::nonlinear break; case FiniteDiffStrategy::FULL_FINITE_DIFF: { - Eigen::VectorXd grad_fd; + TVector grad_fd; fd::finite_gradient( - x, [&](const Eigen::VectorXd &x_) { + x, [&](const TVector &x_) { objFunc.solution_changed(x_); - return objFunc.value(x_); + return objFunc(x_); }, grad_fd, fd::AccuracyOrder::SECOND, gradient_fd_eps); diff --git a/src/polysolve/nonlinear/Solver.hpp b/src/polysolve/nonlinear/Solver.hpp index c66a8568..c00f7397 100644 --- a/src/polysolve/nonlinear/Solver.hpp +++ b/src/polysolve/nonlinear/Solver.hpp @@ -1,11 +1,10 @@ #pragma once +#include "Criteria.hpp" #include "descent_strategies/DescentStrategy.hpp" // Line search methods #include "line_search/LineSearch.hpp" -#include - namespace spdlog { class logger; @@ -13,14 +12,6 @@ namespace spdlog namespace polysolve::nonlinear { - - enum class ErrorCode - { - NAN_ENCOUNTERED = -10, - STEP_TOO_SMALL = -1, - SUCCESS = 0, - }; - enum class FiniteDiffStrategy { NONE = 0, @@ -28,9 +19,16 @@ namespace polysolve::nonlinear FULL_FINITE_DIFF = 2 }; - class Solver : public cppoptlib::ISolver + class Solver { public: + using Scalar = typename Problem::Scalar; + using TVector = typename Problem::TVector; + using THessian = typename Problem::THessian; + + public: + // --- Static methods ------------------------------------------------- + // Static constructor // // @param[in] solver Solver type @@ -46,11 +44,7 @@ namespace polysolve::nonlinear // List available solvers static std::vector available_solvers(); - using Superclass = ISolver; - using typename Superclass::Scalar; - using typename Superclass::TCriteria; - using typename Superclass::TVector; - + public: /// @brief Construct a new Nonlinear Solver object /// @param solver_params JSON of solver parameters /// @param characteristic_length used to scale tolerances @@ -59,34 +53,36 @@ namespace polysolve::nonlinear const double characteristic_length, spdlog::logger &logger); - void set_strategies_iterations(const json &solver_params); + virtual ~Solver() = default; - virtual double compute_grad_norm(const Eigen::VectorXd &x, const Eigen::VectorXd &grad) const; - - void set_line_search(const json ¶ms); + /// @brief Add a descent strategy to the solver + /// @param s Descent strategy + void add_strategy(const std::shared_ptr &s) { m_strategies.push_back(s); } - void minimize(Problem &objFunc, TVector &x) override; + /// @brief Minimize the objective function + /// @param objFunc Objective function + /// @param x Initial guess + void minimize(Problem &objFunc, TVector &x); - const json &get_info() const { return solver_info; } + virtual double compute_grad_norm(const TVector &x, const TVector &grad) const; - ErrorCode error_code() const { return m_error_code; } + // ==================================================================== + // Getters and setters + // ==================================================================== - const typename Superclass::TCriteria &getStopCriteria() { return this->m_stop; } - // setStopCriteria already in ISolver + Criteria &stop_criteria() { return m_stop; } + const Criteria &stop_criteria() const { return m_stop; } + const Criteria ¤t_criteria() const { return m_current; } + const Status &status() const { return m_status; } - bool converged() const - { - return this->m_status == cppoptlib::Status::XDeltaTolerance - || this->m_status == cppoptlib::Status::FDeltaTolerance - || this->m_status == cppoptlib::Status::GradNormTolerance; - } + void set_strategies_iterations(const json &solver_params); + void set_line_search(const json ¶ms); + const json &info() const { return solver_info; } - size_t max_iterations() const { return this->m_stop.iterations; } - size_t &max_iterations() { return this->m_stop.iterations; } + /// @brief If true the solver will not throw an error if the maximum number of iterations is reached bool allow_out_of_iterations = false; - void add_strategy(const std::shared_ptr &s) { m_strategies.push_back(s); } - + /// @brief Get the line search object const std::shared_ptr &line_search() const { return m_line_search; }; protected: @@ -99,8 +95,19 @@ namespace polysolve::nonlinear return m_strategies[m_descent_strategy]->compute_update_direction(objFunc, x, grad, direction); } + /// @brief Stopping criteria + Criteria m_stop; + + /// @brief Current criteria + Criteria m_current; + + /// @brief Current status + Status m_status = Status::NotStarted; + + /// @brief Index into m_strategies int m_descent_strategy; + /// @brief Logger to use spdlog::logger &m_logger; // ==================================================================== @@ -116,13 +123,8 @@ namespace polysolve::nonlinear // Solver parameters // ==================================================================== - double use_grad_norm_tol; - double first_grad_norm_tol; - double derivative_along_delta_x_tol; - const double characteristic_length; - int f_delta_step_tol; // ==================================================================== // Solver state // ==================================================================== @@ -147,14 +149,13 @@ namespace polysolve::nonlinear json solver_info; + // Timers double total_time; double grad_time; double line_search_time; double constraint_set_update_time; double obj_fun_time; - ErrorCode m_error_code; - // ==================================================================== // END // ==================================================================== diff --git a/src/polysolve/nonlinear/descent_strategies/Newton.cpp b/src/polysolve/nonlinear/descent_strategies/Newton.cpp index 7bb88a3b..2fea6d67 100644 --- a/src/polysolve/nonlinear/descent_strategies/Newton.cpp +++ b/src/polysolve/nonlinear/descent_strategies/Newton.cpp @@ -136,8 +136,7 @@ namespace polysolve::nonlinear TVector &direction) { const double residual = - is_sparse ? // - solve_sparse_linear_system(objFunc, x, grad, direction) + is_sparse ? solve_sparse_linear_system(objFunc, x, grad, direction) : solve_dense_linear_system(objFunc, x, grad, direction); if (std::isnan(residual) || residual > m_residual_tolerance * m_characteristic_length) diff --git a/src/polysolve/nonlinear/line_search/Backtracking.cpp b/src/polysolve/nonlinear/line_search/Backtracking.cpp index a783443d..da8ee665 100644 --- a/src/polysolve/nonlinear/line_search/Backtracking.cpp +++ b/src/polysolve/nonlinear/line_search/Backtracking.cpp @@ -45,7 +45,7 @@ namespace polysolve::nonlinear::line_search continue; } - const double new_energy = objFunc.value(new_x); + const double new_energy = objFunc(new_x); if (!std::isfinite(new_energy)) { diff --git a/src/polysolve/nonlinear/line_search/CMakeLists.txt b/src/polysolve/nonlinear/line_search/CMakeLists.txt index 4533d7ae..a1966da3 100644 --- a/src/polysolve/nonlinear/line_search/CMakeLists.txt +++ b/src/polysolve/nonlinear/line_search/CMakeLists.txt @@ -5,10 +5,6 @@ set(SOURCES Armijo.hpp Backtracking.cpp Backtracking.hpp - CppOptArmijo.cpp - CppOptArmijo.hpp - MoreThuente.cpp - MoreThuente.hpp NoLineSearch.cpp NoLineSearch.hpp RobustArmijo.cpp diff --git a/src/polysolve/nonlinear/line_search/CppOptArmijo.cpp b/src/polysolve/nonlinear/line_search/CppOptArmijo.cpp deleted file mode 100644 index f4144d1f..00000000 --- a/src/polysolve/nonlinear/line_search/CppOptArmijo.cpp +++ /dev/null @@ -1,26 +0,0 @@ -#include "CppOptArmijo.hpp" - -#include - -namespace polysolve::nonlinear::line_search -{ - CppOptArmijo::CppOptArmijo(const json ¶ms, spdlog::logger &logger) - : Superclass(params, logger), after_check(params, logger) - { - } - - double CppOptArmijo::compute_descent_step_size( - const TVector &x, - const TVector &delta_x, - Problem &objFunc, - const bool use_grad_norm, - const double old_energy, - const TVector &old_grad, - const double starting_step_size) - { - double step_size = cppoptlib::Armijo::linesearch(x, delta_x, objFunc, starting_step_size); - // this ensures no collisions and decrease in energy - return after_check.compute_descent_step_size(x, delta_x, objFunc, use_grad_norm, old_energy, old_grad, step_size); - } - -} // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/CppOptArmijo.hpp b/src/polysolve/nonlinear/line_search/CppOptArmijo.hpp deleted file mode 100644 index fc277132..00000000 --- a/src/polysolve/nonlinear/line_search/CppOptArmijo.hpp +++ /dev/null @@ -1,32 +0,0 @@ -#pragma once - -#include "LineSearch.hpp" -#include "Backtracking.hpp" - -namespace polysolve::nonlinear::line_search -{ - class CppOptArmijo : public LineSearch - { - public: - using Superclass = LineSearch; - using typename Superclass::Scalar; - using typename Superclass::TVector; - - CppOptArmijo(const json ¶ms, spdlog::logger &logger); - - virtual std::string name() override { return "ArmijoAlt"; } - - protected: - double compute_descent_step_size( - const TVector &x, - const TVector &delta_x, - Problem &objFunc, - const bool use_grad_norm, - const double old_energy, - const TVector &old_grad, - const double starting_step_size) override; - - private: - Backtracking after_check; - }; -} // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/LineSearch.cpp b/src/polysolve/nonlinear/line_search/LineSearch.cpp index e3f9cb7f..f30e5a24 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.cpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.cpp @@ -3,8 +3,6 @@ #include "Armijo.hpp" #include "Backtracking.hpp" #include "RobustArmijo.hpp" -#include "CppOptArmijo.hpp" -#include "MoreThuente.hpp" #include "NoLineSearch.hpp" #include @@ -22,32 +20,19 @@ namespace polysolve::nonlinear::line_search std::shared_ptr LineSearch::create(const json ¶ms, spdlog::logger &logger) { const std::string name = params["line_search"]["method"]; - if (name == "armijo" || name == "Armijo") + if (name == "Armijo") { return std::make_shared(params, logger); } - else if (name == "armijo_alt" || name == "ArmijoAlt") - { - return std::make_shared(params, logger); - } - else if (name == "robust_armijo" || name == "RobustArmijo") + else if (name == "RobustArmijo") { return std::make_shared(params, logger); } - else if (name == "bisection" || name == "Bisection") + else if (name == "Backtracking") { - logger.warn("{} linesearch was renamed to \"backtracking\"; using backtracking line-search", name); return std::make_shared(params, logger); } - else if (name == "backtracking" || name == "Backtracking") - { - return std::make_shared(params, logger); - } - else if (name == "more_thuente" || name == "MoreThuente") - { - return std::make_shared(params, logger); - } - else if (name == "none" || name == "None") + else if (name == "None") { return std::make_shared(params, logger); } @@ -60,12 +45,7 @@ namespace polysolve::nonlinear::line_search std::vector LineSearch::available_methods() { - return {{"Armijo", - "ArmijoAlt", - "RobustArmijo", - "Backtracking", - "MoreThuente", - "None"}}; + return {{"Armijo", "RobustArmijo", "Backtracking", "None"}}; } LineSearch::LineSearch(const json ¶ms, spdlog::logger &logger) @@ -96,7 +76,7 @@ namespace polysolve::nonlinear::line_search cur_iter = 0; - initial_energy = objFunc.value(x); + initial_energy = objFunc(x); if (std::isnan(initial_energy)) { m_logger.error("Original energy in line search is nan!"); @@ -166,7 +146,7 @@ namespace polysolve::nonlinear::line_search } } - const double cur_energy = objFunc.value(x + step_size * delta_x); + const double cur_energy = objFunc(x + step_size * delta_x); const double descent_step_size = step_size; @@ -179,7 +159,7 @@ namespace polysolve::nonlinear::line_search objFunc.solution_changed(x); // tolerance for rounding error due to multithreading - assert(abs(initial_energy - objFunc.value(x)) < 1e-15); + assert(abs(initial_energy - objFunc(x)) < 1e-15); objFunc.line_search_end(); return NaN; @@ -211,7 +191,7 @@ namespace polysolve::nonlinear::line_search while (step_size > current_min_step_size() && cur_iter < current_max_step_size_iter()) { // Compute the new energy value without contacts - const double energy = objFunc.value(new_x); + const double energy = objFunc(new_x); const bool is_step_valid = objFunc.is_step_valid(x, new_x); if (!std::isfinite(energy) || !is_step_valid) diff --git a/src/polysolve/nonlinear/line_search/MoreThuente.cpp b/src/polysolve/nonlinear/line_search/MoreThuente.cpp deleted file mode 100644 index 0c78bf93..00000000 --- a/src/polysolve/nonlinear/line_search/MoreThuente.cpp +++ /dev/null @@ -1,25 +0,0 @@ -#include "MoreThuente.hpp" - -#include - -namespace polysolve::nonlinear::line_search -{ - MoreThuente::MoreThuente(const json ¶ms, spdlog::logger &logger) - : Superclass(params, logger), after_check(params, logger) - { - } - - double MoreThuente::compute_descent_step_size( - const TVector &x, - const TVector &delta_x, - Problem &objFunc, - const bool use_grad_norm, - const double old_energy, - const TVector &old_grad, - const double starting_step_size) - { - const double tmp = cppoptlib::MoreThuente::linesearch(x, delta_x, objFunc, starting_step_size); - // this ensures no collisions and decrease in energy - return after_check.compute_descent_step_size(x, delta_x, objFunc, use_grad_norm, old_energy, old_grad, std::min(tmp, starting_step_size)); - } -} // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/MoreThuente.hpp b/src/polysolve/nonlinear/line_search/MoreThuente.hpp deleted file mode 100644 index 96f7765b..00000000 --- a/src/polysolve/nonlinear/line_search/MoreThuente.hpp +++ /dev/null @@ -1,32 +0,0 @@ -#pragma once - -#include "LineSearch.hpp" -#include "Backtracking.hpp" - -namespace polysolve::nonlinear::line_search -{ - class MoreThuente : public LineSearch - { - public: - using Superclass = LineSearch; - using typename Superclass::Scalar; - using typename Superclass::TVector; - - MoreThuente(const json ¶ms, spdlog::logger &logger); - - virtual std::string name() override { return "MoreThuente"; } - - protected: - double compute_descent_step_size( - const TVector &x, - const TVector &delta_x, - Problem &objFunc, - const bool use_grad_norm, - const double old_energy, - const TVector &old_grad, - const double starting_step_size) override; - - private: - Backtracking after_check; - }; -} // namespace polysolve::nonlinear::line_search