From 00474f435ed2e4a7de62da2d72d51d4c9d64327b Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Fri, 5 Jan 2024 13:17:39 -0500 Subject: [PATCH 01/13] Drop cppoptlib from Problem and Solver --- src/polysolve/nonlinear/Criteria.hpp | 130 ++++++++++++++++++ src/polysolve/nonlinear/Problem.hpp | 100 ++++++++++++-- src/polysolve/nonlinear/Solver.cpp | 126 +++++++++-------- src/polysolve/nonlinear/Solver.hpp | 44 +++--- .../nonlinear/descent_strategies/Newton.cpp | 3 +- 5 files changed, 310 insertions(+), 93 deletions(-) create mode 100644 src/polysolve/nonlinear/Criteria.hpp diff --git a/src/polysolve/nonlinear/Criteria.hpp b/src/polysolve/nonlinear/Criteria.hpp new file mode 100644 index 0000000..aaa3082 --- /dev/null +++ b/src/polysolve/nonlinear/Criteria.hpp @@ -0,0 +1,130 @@ +#pragma once + +#include +#include + +namespace polysolve::nonlinear +{ + enum class Status + { + NotStarted = -1, + Continue = 0, + IterationLimit, + XDeltaTolerance, + FDeltaTolerance, + GradNormTolerance, + Condition, + UserDefined + }; + + template + class Criteria + { + public: + size_t iterations; //!< Maximum number of iterations + T xDelta; //!< Minimum change in parameter vector + T fDelta; //!< Minimum change in cost function + T gradNorm; //!< Minimum norm of gradient vector + T condition; //!< Maximum condition number of Hessian + + Criteria() + { + reset(); + } + + static Criteria defaults() + { + Criteria d; + d.iterations = 10000; + d.xDelta = 0; + d.fDelta = 0; + d.gradNorm = 1e-4; + d.condition = 0; + return d; + } + + void reset() + { + iterations = 0; + xDelta = 0; + fDelta = 0; + gradNorm = 0; + condition = 0; + } + + void 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 << "Condition: " << condition << std::endl; + } + }; + + template + 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)) + { + return Status::FDeltaTolerance; + } + if ((stop.gradNorm > 0) && (current.gradNorm < stop.gradNorm)) + { + return Status::GradNormTolerance; + } + if ((stop.condition > 0) && (current.condition > stop.condition)) + { + return Status::Condition; + } + return Status::Continue; + } + + inline 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::Condition: + os << "Condition of Hessian/Covariance matrix too large."; + break; + case Status::UserDefined: + os << "Stop condition defined in the callback."; + break; + } + return os; + } + + template + 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/Problem.hpp b/src/polysolve/nonlinear/Problem.hpp index 93fa0db..024d1a6 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -2,47 +2,119 @@ #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 const 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; + + // TODO: Add dense Hessian + + /// @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) + { + StiffnessMatrix sparse_hessian; + hessian(x, sparse_hessian); + hessian = sparse_hessian; + } + + /// @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; - virtual bool is_step_valid(const TVector &x0, const TVector &x1) { return true; } - virtual double max_step_size(const TVector &x0, const TVector &x1) { return 1; } + /// @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) const { 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) const { 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 fa55cbd..e751969 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -196,18 +196,19 @@ 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"]; + m_current.reset(); - criteria.xDelta *= characteristic_length; - criteria.fDelta *= characteristic_length; - criteria.gradNorm *= characteristic_length; + m_stop = TCriteria::defaults(); + m_stop.xDelta = solver_params["x_delta"]; + m_stop.fDelta = solver_params["advanced"]["f_delta"]; + m_stop.gradNorm = solver_params["grad_norm"]; - criteria.iterations = solver_params["max_iterations"]; - // criteria.condition = solver_params["condition"]; - this->setStopCriteria(criteria); + m_stop.xDelta *= characteristic_length; + m_stop.fDelta *= characteristic_length; + m_stop.gradNorm *= characteristic_length; + + m_stop.iterations = solver_params["max_iterations"]; + // m_stop.condition = solver_params["condition"]; use_grad_norm_tol = solver_params["line_search"]["use_grad_norm_tol"]; first_grad_norm_tol = solver_params["first_grad_norm_tol"]; @@ -275,21 +276,21 @@ namespace polysolve::nonlinear objFunc.solution_changed(x); } - const auto g_norm_tol = this->m_stop.gradNorm; - this->m_stop.gradNorm = first_grad_norm_tol; + const auto g_norm_tol = m_stop.gradNorm; + 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.value(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)); + objFunc.post_step(PostStepData(m_current.iterations, solver_info, x, grad)); int f_delta_step_cnt = 0; double f_delta = 0; @@ -301,9 +302,9 @@ namespace polysolve::nonlinear { 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; + m_current.xDelta = NaN; + m_current.fDelta = NaN; + m_current.gradNorm = NaN; //////////// Energy double energy; @@ -314,7 +315,7 @@ namespace polysolve::nonlinear if (!std::isfinite(energy)) { - this->m_status = cppoptlib::Status::UserDefined; + m_status = Status::UserDefined; m_error_code = ErrorCode::NAN_ENCOUNTERED; log_and_throw_error(m_logger, "[{}][{}] f(x) is nan or inf; stopping", descent_strategy_name(), m_line_search->name()); break; @@ -322,7 +323,7 @@ namespace polysolve::nonlinear 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 = (f_delta_step_cnt >= f_delta_step_tol) ? f_delta : NaN; ///////////// gradient { @@ -332,23 +333,23 @@ 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)) { - this->m_status = cppoptlib::Status::UserDefined; + m_status = Status::UserDefined; m_error_code = ErrorCode::NAN_ENCOUNTERED; 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; + m_current.gradNorm = grad_norm; + gradNorm = 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; // ------------------------ @@ -366,17 +367,20 @@ namespace polysolve::nonlinear ++m_descent_strategy; if (m_descent_strategy >= m_strategies.size()) { - this->m_status = cppoptlib::Status::UserDefined; + m_status = 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), delta_x.dot(grad)); } - 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), delta_x.dot(grad), descent_strategy_name()); - this->m_status = cppoptlib::Status::Continue; + if (ok) // if ok, then the direction is not a descent direction + { + 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), delta_x.dot(grad), descent_strategy_name()); + } + m_status = Status::Continue; continue; } @@ -389,22 +393,22 @@ namespace polysolve::nonlinear if (m_descent_strategy >= m_strategies.size()) { - this->m_status = cppoptlib::Status::UserDefined; + m_status = Status::UserDefined; 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_status = 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) + m_current.xDelta = delta_x_norm; + xDelta = m_current.xDelta; + m_status = checkConvergence(m_stop, m_current); + if (m_status != Status::Continue) break; // --------------- @@ -414,19 +418,19 @@ namespace polysolve::nonlinear 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, 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! + m_status = Status::UserDefined; // Line search failed on gradient descent, so quit! log_and_throw_error(m_logger, "[{}][{}] Line search failed on last strategy; stopping", current_name, m_line_search->name()); } @@ -466,16 +470,16 @@ 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_status = Status::UserDefined; m_error_code = ErrorCode::SUCCESS; m_logger.debug("[{}][{}] Objective decided to stop", descent_strategy_name(), m_line_search->name()); } - if (f_delta < this->m_stop.fDelta) + if (f_delta < m_stop.fDelta) f_delta_step_cnt++; else f_delta_step_cnt = 0; @@ -484,16 +488,16 @@ namespace polysolve::nonlinear "[{}][{}] 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, + 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); + m_stop.iterations, m_stop.fDelta, m_stop.gradNorm, m_stop.xDelta); - if (++this->m_current.iterations >= this->m_stop.iterations) - this->m_status = cppoptlib::Status::IterationLimit; + if (++m_current.iterations >= m_stop.iterations) + m_status = Status::IterationLimit; // 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)); + m_stop.gradNorm = g_norm_tol; + } while (objFunc.callback(m_current, x) && (m_status == Status::Continue)); stop_watch.stop(); @@ -501,21 +505,21 @@ 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::UserDefined && m_error_code != ErrorCode::SUCCESS) 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, + m_status, tot_time, 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_stop.iterations, m_stop.fDelta, m_stop.gradNorm, m_stop.xDelta); log_times(); update_solver_info(objFunc.value(x)); @@ -523,7 +527,7 @@ namespace polysolve::nonlinear void Solver::reset(const int ndof) { - this->m_current.reset(); + m_current.reset(); m_descent_strategy = 0; m_error_code = ErrorCode::SUCCESS; @@ -555,10 +559,10 @@ namespace polysolve::nonlinear void Solver::update_solver_info(const double energy) { - solver_info["status"] = this->status(); + solver_info["status"] = status(); solver_info["error_code"] = m_error_code; solver_info["energy"] = energy; - const auto &crit = this->criteria(); + const auto &crit = criteria(); solver_info["iterations"] = crit.iterations; solver_info["xDelta"] = crit.xDelta; solver_info["fDelta"] = crit.fDelta; diff --git a/src/polysolve/nonlinear/Solver.hpp b/src/polysolve/nonlinear/Solver.hpp index 45ee0b9..f1f0203 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; @@ -28,9 +27,17 @@ 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; + using TCriteria = Criteria; + + public: + // --- Static methods ------------------------------------------------- + // Static constructor // // @param[in] solver Solver type @@ -46,11 +53,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,26 +62,31 @@ namespace polysolve::nonlinear const double characteristic_length, spdlog::logger &logger); + virtual ~Solver() = default; + + const TCriteria &get_stop_criteria() { return this->m_stop; } + void set_stop_criteria(const TCriteria &s) { m_stop = s; } + + const TCriteria &criteria() { return m_current; } + const Status &status() { return m_status; } + void set_strategies_iterations(const json &solver_params); - virtual double compute_grad_norm(const Eigen::VectorXd &x, const Eigen::VectorXd &grad) const; + virtual double compute_grad_norm(const TVector &x, const TVector &grad) const; void set_line_search(const json ¶ms); - void minimize(Problem &objFunc, TVector &x) override; + void minimize(Problem &objFunc, TVector &x); const json &get_info() const { return solver_info; } ErrorCode error_code() const { return m_error_code; } - const typename Superclass::TCriteria &getStopCriteria() { return this->m_stop; } - // setStopCriteria already in ISolver - bool converged() const { - return this->m_status == cppoptlib::Status::XDeltaTolerance - || this->m_status == cppoptlib::Status::FDeltaTolerance - || this->m_status == cppoptlib::Status::GradNormTolerance; + return this->m_status == Status::XDeltaTolerance + || this->m_status == Status::FDeltaTolerance + || this->m_status == Status::GradNormTolerance; } size_t max_iterations() const { return this->m_stop.iterations; } @@ -99,6 +107,10 @@ namespace polysolve::nonlinear return m_strategies[m_descent_strategy]->compute_update_direction(objFunc, x, grad, direction); } + TCriteria m_stop; + TCriteria m_current; + Status m_status = Status::NotStarted; + int m_descent_strategy; spdlog::logger &m_logger; diff --git a/src/polysolve/nonlinear/descent_strategies/Newton.cpp b/src/polysolve/nonlinear/descent_strategies/Newton.cpp index 7bb88a3..2fea6d6 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) From a094f637f1e26d3bfdcae813266eec9faf94937d Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Fri, 5 Jan 2024 13:45:14 -0500 Subject: [PATCH 02/13] Cleanup changes --- src/polysolve/nonlinear/Criteria.hpp | 5 ++++- src/polysolve/nonlinear/Problem.hpp | 10 +++------- src/polysolve/nonlinear/Solver.cpp | 14 +++++++------- src/polysolve/nonlinear/Solver.hpp | 12 ++++++------ .../nonlinear/line_search/Backtracking.cpp | 2 +- src/polysolve/nonlinear/line_search/LineSearch.cpp | 8 ++++---- 6 files changed, 25 insertions(+), 26 deletions(-) diff --git a/src/polysolve/nonlinear/Criteria.hpp b/src/polysolve/nonlinear/Criteria.hpp index aaa3082..7aad171 100644 --- a/src/polysolve/nonlinear/Criteria.hpp +++ b/src/polysolve/nonlinear/Criteria.hpp @@ -1,10 +1,13 @@ #pragma once +#include #include -#include namespace polysolve::nonlinear { + // Source: https://github.com/PatWie/CppNumericalSolvers/blob/7eddf28fa5a8872a956d3c8666055cac2f5a535d/include/cppoptlib/meta.h + // License: MIT + enum class Status { NotStarted = -1, diff --git a/src/polysolve/nonlinear/Problem.hpp b/src/polysolve/nonlinear/Problem.hpp index 024d1a6..3a9ecbe 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -42,16 +42,12 @@ namespace polysolve::nonlinear /// @param[out] grad Gradient of the function at x. virtual void gradient(const TVector &x, TVector &grad) = 0; - // TODO: Add dense Hessian - /// @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) { - StiffnessMatrix sparse_hessian; - hessian(x, sparse_hessian); - hessian = sparse_hessian; + throw std::runtime_error("Dense Hessian not implemented."); } /// @brief Compute the Hessian of the function at x. @@ -63,13 +59,13 @@ namespace polysolve::nonlinear /// @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) const { return true; } + 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) const { return 1; } + virtual double max_step_size(const TVector &x0, const TVector &x1) { return 1; } // --- Callbacks ------------------------------------------------------ diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index e751969..c89eb05 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -286,10 +286,10 @@ namespace polysolve::nonlinear "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), m_stop.iterations, + objFunc(x), m_stop.iterations, m_stop.fDelta, m_stop.gradNorm, m_stop.xDelta); - update_solver_info(objFunc.value(x)); + update_solver_info(objFunc(x)); objFunc.post_step(PostStepData(m_current.iterations, solver_info, x, grad)); int f_delta_step_cnt = 0; @@ -310,7 +310,7 @@ namespace polysolve::nonlinear double energy; { POLYSOLVE_SCOPED_STOPWATCH("compute objective function", obj_fun_time, m_logger); - energy = objFunc.value(x); + energy = objFunc(x); } if (!std::isfinite(energy)) @@ -522,7 +522,7 @@ namespace polysolve::nonlinear 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) @@ -630,10 +630,10 @@ namespace polysolve::nonlinear Eigen::VectorXd 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); @@ -653,7 +653,7 @@ namespace polysolve::nonlinear fd::finite_gradient( x, [&](const Eigen::VectorXd &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 f1f0203..1096716 100644 --- a/src/polysolve/nonlinear/Solver.hpp +++ b/src/polysolve/nonlinear/Solver.hpp @@ -64,7 +64,7 @@ namespace polysolve::nonlinear virtual ~Solver() = default; - const TCriteria &get_stop_criteria() { return this->m_stop; } + const TCriteria &get_stop_criteria() { return m_stop; } void set_stop_criteria(const TCriteria &s) { m_stop = s; } const TCriteria &criteria() { return m_current; } @@ -84,13 +84,13 @@ namespace polysolve::nonlinear bool converged() const { - return this->m_status == Status::XDeltaTolerance - || this->m_status == Status::FDeltaTolerance - || this->m_status == Status::GradNormTolerance; + return m_status == Status::XDeltaTolerance + || m_status == Status::FDeltaTolerance + || m_status == Status::GradNormTolerance; } - size_t max_iterations() const { return this->m_stop.iterations; } - size_t &max_iterations() { return this->m_stop.iterations; } + size_t max_iterations() const { return m_stop.iterations; } + size_t &max_iterations() { return m_stop.iterations; } bool allow_out_of_iterations = false; void add_strategy(const std::shared_ptr &s) { m_strategies.push_back(s); } diff --git a/src/polysolve/nonlinear/line_search/Backtracking.cpp b/src/polysolve/nonlinear/line_search/Backtracking.cpp index a783443..da8ee66 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/LineSearch.cpp b/src/polysolve/nonlinear/line_search/LineSearch.cpp index e3f9cb7..b3c2f98 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.cpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.cpp @@ -96,7 +96,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 +166,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 +179,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 +211,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) From ce5b3a2dc0d3fcd28eb8d22f7521a5dbeaeac788 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Fri, 5 Jan 2024 13:50:18 -0500 Subject: [PATCH 03/13] Private CMake target_link_libraries --- CMakeLists.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 587f37e..b868c2a 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -313,15 +313,15 @@ target_link_libraries(polysolve PUBLIC polysolve::linear) # CppNumericalSolvers include(cppoptlib) -target_link_libraries(polysolve PUBLIC cppoptlib) +target_link_libraries(polysolve PRIVATE 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) From 272681c11108787c57c614568c5a201197d05690 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Fri, 5 Jan 2024 13:53:32 -0500 Subject: [PATCH 04/13] Use TVector --- src/polysolve/nonlinear/Solver.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index c89eb05..172f3ae 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -242,7 +242,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(); } @@ -625,9 +625,9 @@ 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(x2); @@ -649,9 +649,9 @@ 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(x_); }, From 7830b156796dedfc867da7b8ad61da207fd95edc Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Thu, 8 Feb 2024 16:37:41 -0500 Subject: [PATCH 05/13] Remove cppoptlib line-search methods --- CMakeLists.txt | 4 --- cmake/recipes/cppoptlib.cmake | 19 ----------- nonlinear-solver-spec.json | 5 +-- .../nonlinear/line_search/CMakeLists.txt | 4 --- .../nonlinear/line_search/CppOptArmijo.cpp | 26 --------------- .../nonlinear/line_search/CppOptArmijo.hpp | 32 ------------------- .../nonlinear/line_search/LineSearch.cpp | 10 ------ .../nonlinear/line_search/MoreThuente.cpp | 25 --------------- .../nonlinear/line_search/MoreThuente.hpp | 32 ------------------- 9 files changed, 1 insertion(+), 156 deletions(-) delete mode 100644 cmake/recipes/cppoptlib.cmake 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 b868c2a..a3ccfe6 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -311,10 +311,6 @@ target_link_libraries(polysolve_linear PRIVATE polysolve::warnings) # polysolve::linear target_link_libraries(polysolve PUBLIC polysolve::linear) -# CppNumericalSolvers -include(cppoptlib) -target_link_libraries(polysolve PRIVATE cppoptlib) - # LBFGSpp include(LBFGSpp) target_link_libraries(polysolve PRIVATE LBFGSpp::LBFGSpp) diff --git a/cmake/recipes/cppoptlib.cmake b/cmake/recipes/cppoptlib.cmake deleted file mode 100644 index e19709e..0000000 --- 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/nonlinear-solver-spec.json b/nonlinear-solver-spec.json index 6e9baa4..7887519 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/line_search/CMakeLists.txt b/src/polysolve/nonlinear/line_search/CMakeLists.txt index 4533d7a..a1966da 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 f4144d1..0000000 --- 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 fc27713..0000000 --- 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 b3c2f98..83dcb9d 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 @@ -26,10 +24,6 @@ namespace polysolve::nonlinear::line_search { 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") { return std::make_shared(params, logger); @@ -43,10 +37,6 @@ namespace polysolve::nonlinear::line_search { 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") { return std::make_shared(params, logger); diff --git a/src/polysolve/nonlinear/line_search/MoreThuente.cpp b/src/polysolve/nonlinear/line_search/MoreThuente.cpp deleted file mode 100644 index 0c78bf9..0000000 --- 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 96f7765..0000000 --- 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 From 85df8457ac74a29dc41f4832634c10ece4856f04 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Thu, 8 Feb 2024 17:31:38 -0500 Subject: [PATCH 06/13] Add custom statuses --- src/polysolve/nonlinear/CMakeLists.txt | 12 ++- src/polysolve/nonlinear/Criteria.cpp | 84 ++++++++++++++++ src/polysolve/nonlinear/Criteria.hpp | 130 +++++-------------------- src/polysolve/nonlinear/Problem.hpp | 2 +- src/polysolve/nonlinear/Solver.cpp | 24 ++--- src/polysolve/nonlinear/Solver.hpp | 26 ++--- 6 files changed, 129 insertions(+), 149 deletions(-) create mode 100644 src/polysolve/nonlinear/Criteria.cpp diff --git a/src/polysolve/nonlinear/CMakeLists.txt b/src/polysolve/nonlinear/CMakeLists.txt index c842dbe..dfbb760 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 0000000..aaa48e8 --- /dev/null +++ b/src/polysolve/nonlinear/Criteria.cpp @@ -0,0 +1,84 @@ +// Source: https://github.com/PatWie/CppNumericalSolvers/blob/7eddf28fa5a8872a956d3c8666055cac2f5a535d/include/cppoptlib/meta.h +// License: MIT +#include "Criteria.hpp" + +namespace polysolve::nonlinear +{ + Criteria::Criteria() + { + reset(); + } + + void Criteria::reset() + { + iterations = 0; + xDelta = 0; + fDelta = 0; + gradNorm = 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; + } + + 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)) + { + return Status::FDeltaTolerance; + } + if ((stop.gradNorm > 0) && (current.gradNorm < stop.gradNorm)) + { + return Status::GradNormTolerance; + } + 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; + 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 index 7aad171..63bb856 100644 --- a/src/polysolve/nonlinear/Criteria.hpp +++ b/src/polysolve/nonlinear/Criteria.hpp @@ -10,124 +10,38 @@ namespace polysolve::nonlinear enum class Status { - NotStarted = -1, - Continue = 0, - IterationLimit, - XDeltaTolerance, - FDeltaTolerance, - GradNormTolerance, - Condition, - UserDefined + 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 }; - template class Criteria { public: - size_t iterations; //!< Maximum number of iterations - T xDelta; //!< Minimum change in parameter vector - T fDelta; //!< Minimum change in cost function - T gradNorm; //!< Minimum norm of gradient vector - T condition; //!< Maximum condition number of Hessian + 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 - Criteria() - { - reset(); - } + Criteria(); - static Criteria defaults() - { - Criteria d; - d.iterations = 10000; - d.xDelta = 0; - d.fDelta = 0; - d.gradNorm = 1e-4; - d.condition = 0; - return d; - } + void reset(); - void reset() - { - iterations = 0; - xDelta = 0; - fDelta = 0; - gradNorm = 0; - condition = 0; - } - - void 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 << "Condition: " << condition << std::endl; - } + void print(std::ostream &os) const; }; - template - Status checkConvergence(const Criteria &stop, const Criteria ¤t) - { + 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)) - { - return Status::FDeltaTolerance; - } - if ((stop.gradNorm > 0) && (current.gradNorm < stop.gradNorm)) - { - return Status::GradNormTolerance; - } - if ((stop.condition > 0) && (current.condition > stop.condition)) - { - return Status::Condition; - } - return Status::Continue; - } + std::ostream &operator<<(std::ostream &os, const Status &s); - inline 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::Condition: - os << "Condition of Hessian/Covariance matrix too large."; - break; - case Status::UserDefined: - os << "Stop condition defined in the callback."; - break; - } - return os; - } - - template - std::ostream &operator<<(std::ostream &os, const Criteria &c) - { - c.print(os); - return os; - } + 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 3a9ecbe..8df1454 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -93,7 +93,7 @@ namespace polysolve::nonlinear /// @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; } + 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. diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index 172f3ae..90d4cdd 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -198,7 +198,6 @@ namespace polysolve::nonlinear { m_current.reset(); - m_stop = TCriteria::defaults(); m_stop.xDelta = solver_params["x_delta"]; m_stop.fDelta = solver_params["advanced"]["f_delta"]; m_stop.gradNorm = solver_params["grad_norm"]; @@ -315,8 +314,7 @@ namespace polysolve::nonlinear if (!std::isfinite(energy)) { - m_status = 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; } @@ -339,8 +337,7 @@ namespace polysolve::nonlinear const double grad_norm = compute_grad_norm(x, grad); if (std::isnan(grad_norm)) { - m_status = 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; } @@ -367,7 +364,7 @@ namespace polysolve::nonlinear ++m_descent_strategy; if (m_descent_strategy >= m_strategies.size()) { - m_status = Status::UserDefined; + 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), delta_x.dot(grad)); @@ -393,12 +390,11 @@ namespace polysolve::nonlinear if (m_descent_strategy >= m_strategies.size()) { - m_status = 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()); } - m_status = Status::UserDefined; m_logger.debug("[{}][{}] Δx is nan; reverting to {}", current_name, m_line_search->name(), descent_strategy_name()); m_status = Status::Continue; continue; @@ -430,7 +426,8 @@ namespace polysolve::nonlinear ++m_descent_strategy; if (m_descent_strategy >= m_strategies.size()) { - m_status = 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()); } @@ -474,8 +471,7 @@ namespace polysolve::nonlinear if (objFunc.stop(x)) { - m_status = 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()); } @@ -507,7 +503,7 @@ namespace polysolve::nonlinear 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::UserDefined && m_error_code != ErrorCode::SUCCESS) + 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(); @@ -529,7 +525,7 @@ namespace polysolve::nonlinear { 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(); @@ -560,14 +556,12 @@ namespace polysolve::nonlinear void Solver::update_solver_info(const double energy) { solver_info["status"] = status(); - solver_info["error_code"] = m_error_code; solver_info["energy"] = energy; const auto &crit = 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; double per_iteration = crit.iterations ? crit.iterations : 1; diff --git a/src/polysolve/nonlinear/Solver.hpp b/src/polysolve/nonlinear/Solver.hpp index 1096716..eb7ea05 100644 --- a/src/polysolve/nonlinear/Solver.hpp +++ b/src/polysolve/nonlinear/Solver.hpp @@ -12,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, @@ -33,7 +25,6 @@ namespace polysolve::nonlinear using Scalar = typename Problem::Scalar; using TVector = typename Problem::TVector; using THessian = typename Problem::THessian; - using TCriteria = Criteria; public: // --- Static methods ------------------------------------------------- @@ -64,11 +55,10 @@ namespace polysolve::nonlinear virtual ~Solver() = default; - const TCriteria &get_stop_criteria() { return m_stop; } - void set_stop_criteria(const TCriteria &s) { m_stop = s; } - - const TCriteria &criteria() { return m_current; } - const Status &status() { return m_status; } + Criteria &stop_criteria() { return m_stop; } + const Criteria &stop_criteria() const { return m_stop; } + const Criteria &criteria() const { return m_current; } + const Status &status() const { return m_status; } void set_strategies_iterations(const json &solver_params); @@ -80,8 +70,6 @@ namespace polysolve::nonlinear const json &get_info() const { return solver_info; } - ErrorCode error_code() const { return m_error_code; } - bool converged() const { return m_status == Status::XDeltaTolerance @@ -107,8 +95,8 @@ namespace polysolve::nonlinear return m_strategies[m_descent_strategy]->compute_update_direction(objFunc, x, grad, direction); } - TCriteria m_stop; - TCriteria m_current; + Criteria m_stop; + Criteria m_current; Status m_status = Status::NotStarted; int m_descent_strategy; @@ -164,8 +152,6 @@ namespace polysolve::nonlinear double constraint_set_update_time; double obj_fun_time; - ErrorCode m_error_code; - // ==================================================================== // END // ==================================================================== From 7059e03e51698034c11b0ba66fcaf4819da8d7e3 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Thu, 8 Feb 2024 18:13:08 -0500 Subject: [PATCH 07/13] Move extra solver logic to Criteria --- src/polysolve/nonlinear/Criteria.cpp | 17 +++-- src/polysolve/nonlinear/Criteria.hpp | 12 ++-- src/polysolve/nonlinear/Solver.cpp | 93 +++++++++++----------------- src/polysolve/nonlinear/Solver.hpp | 50 ++++++++------- 4 files changed, 82 insertions(+), 90 deletions(-) diff --git a/src/polysolve/nonlinear/Criteria.cpp b/src/polysolve/nonlinear/Criteria.cpp index aaa48e8..3b30b19 100644 --- a/src/polysolve/nonlinear/Criteria.cpp +++ b/src/polysolve/nonlinear/Criteria.cpp @@ -4,6 +4,11 @@ namespace polysolve::nonlinear { + bool is_converged_status(const Status s) + { + return s == Status::XDeltaTolerance || s == Status::FDeltaTolerance || s == Status::GradNormTolerance; + } + Criteria::Criteria() { reset(); @@ -15,6 +20,8 @@ namespace polysolve::nonlinear xDelta = 0; fDelta = 0; gradNorm = 0; + firstGradNorm = 0; + fDeltaCount = 0; } void Criteria::print(std::ostream &os) const @@ -27,20 +34,20 @@ namespace polysolve::nonlinear Status checkConvergence(const Criteria &stop, const Criteria ¤t) { - - if ((stop.iterations > 0) && (current.iterations > stop.iterations)) + if (stop.iterations > 0 && current.iterations > stop.iterations) { return Status::IterationLimit; } - if ((stop.xDelta > 0) && (current.xDelta < stop.xDelta)) + if (stop.xDelta > 0 && current.xDelta < stop.xDelta) { return Status::XDeltaTolerance; } - if ((stop.fDelta > 0) && (current.fDelta < stop.fDelta)) + if (stop.fDelta > 0 && current.fDelta < stop.fDelta && current.fDeltaCount >= stop.fDeltaCount) { return Status::FDeltaTolerance; } - if ((stop.gradNorm > 0) && (current.gradNorm < stop.gradNorm)) + const double stopGradNorm = current.iterations == 0 ? stop.firstGradNorm : stop.gradNorm; + if (stopGradNorm > 0 && current.gradNorm < stopGradNorm) { return Status::GradNormTolerance; } diff --git a/src/polysolve/nonlinear/Criteria.hpp b/src/polysolve/nonlinear/Criteria.hpp index 63bb856..e3551c7 100644 --- a/src/polysolve/nonlinear/Criteria.hpp +++ b/src/polysolve/nonlinear/Criteria.hpp @@ -24,13 +24,17 @@ namespace polysolve::nonlinear 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 + 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 + unsigned fDeltaCount; ///< Number of steps where fDelta is satisfied Criteria(); diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index 90d4cdd..6a64924 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -201,22 +201,18 @@ namespace polysolve::nonlinear 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"]; + // 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.iterations = solver_params["max_iterations"]; - // m_stop.condition = solver_params["condition"]; - - use_grad_norm_tol = solver_params["line_search"]["use_grad_norm_tol"]; - first_grad_norm_tol = solver_params["first_grad_norm_tol"]; 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; @@ -250,6 +246,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) @@ -263,8 +261,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()); @@ -275,9 +271,6 @@ namespace polysolve::nonlinear objFunc.solution_changed(x); } - const auto g_norm_tol = m_stop.gradNorm; - m_stop.gradNorm = first_grad_norm_tol; - StopWatch stop_watch("nonlinear solver", total_time, m_logger); stop_watch.start(); @@ -291,11 +284,8 @@ namespace polysolve::nonlinear update_solver_info(objFunc(x)); objFunc.post_step(PostStepData(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; + // double xDelta = 0, gradNorm = 0; do { @@ -305,7 +295,8 @@ namespace polysolve::nonlinear m_current.fDelta = NaN; m_current.gradNorm = NaN; - //////////// Energy + // --- Energy ------------------------------------------------------ + double energy; { POLYSOLVE_SCOPED_STOPWATCH("compute objective function", obj_fun_time, m_logger); @@ -319,11 +310,9 @@ namespace polysolve::nonlinear 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 - 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); @@ -334,34 +323,28 @@ namespace polysolve::nonlinear 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)) { m_status = Status::NanEncountered; log_and_throw_error(m_logger, "[{}][{}] Gradient is nan; stopping", descent_strategy_name(), m_line_search->name()); - break; } - m_current.gradNorm = grad_norm; - gradNorm = m_current.gradNorm; - 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); + // --- Update direction -------------------------------------------- - if (!ok || std::isnan(grad_norm) || (m_strategies[m_descent_strategy]->is_direction_descent() && grad_norm != 0 && delta_x.dot(grad) >= 0)) + const bool ok = compute_update_direction(objFunc, x, grad, delta_x); + + if (!ok || (m_strategies[m_descent_strategy]->is_direction_descent() && m_current.gradNorm != 0 && delta_x.dot(grad) >= 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()) { m_status = Status::NotDescentDirection; @@ -377,12 +360,13 @@ namespace polysolve::nonlinear current_name, m_line_search->name(), delta_x.norm(), compute_grad_norm(x, grad), delta_x.dot(grad), descent_strategy_name()); } + m_status = 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()) @@ -400,10 +384,10 @@ namespace polysolve::nonlinear continue; } - // Use the maximum absolute displacement value divided by the timestep, - m_current.xDelta = delta_x_norm; - xDelta = m_current.xDelta; + // --- Check convergence ------------------------------------------- + m_status = checkConvergence(m_stop, m_current); + if (m_status != Status::Continue) break; @@ -414,7 +398,7 @@ namespace polysolve::nonlinear m_logger.trace( "[{}][{}] pre LS iter={:d} f={:g} ‖∇f‖={:g}", descent_strategy_name(), m_line_search->name(), - 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); @@ -475,24 +459,18 @@ namespace polysolve::nonlinear m_logger.debug("[{}][{}] Objective decided to stop", descent_strategy_name(), m_line_search->name()); } - if (f_delta < 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(), - m_current.iterations, energy, f_delta, - gradNorm, xDelta, delta_x.dot(grad), rate, step, + 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); if (++m_current.iterations >= m_stop.iterations) m_status = Status::IterationLimit; - - // reset the tolerance, since in the first iter it might be smaller - m_stop.gradNorm = g_norm_tol; } while (objFunc.callback(m_current, x) && (m_status == Status::Continue)); stop_watch.stop(); @@ -514,7 +492,7 @@ namespace polysolve::nonlinear " (stopping criteria: max_iters={:d} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g})", descent_strategy_name(), m_line_search->name(), m_status, tot_time, m_current.iterations, - old_energy, f_delta, gradNorm, xDelta, + 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(); @@ -557,13 +535,12 @@ namespace polysolve::nonlinear { solver_info["status"] = status(); solver_info["energy"] = energy; - const auto &crit = criteria(); - solver_info["iterations"] = crit.iterations; - solver_info["xDelta"] = crit.xDelta; - solver_info["fDelta"] = crit.fDelta; - solver_info["gradNorm"] = crit.gradNorm; + 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; diff --git a/src/polysolve/nonlinear/Solver.hpp b/src/polysolve/nonlinear/Solver.hpp index eb7ea05..c00f739 100644 --- a/src/polysolve/nonlinear/Solver.hpp +++ b/src/polysolve/nonlinear/Solver.hpp @@ -55,34 +55,34 @@ namespace polysolve::nonlinear virtual ~Solver() = default; - Criteria &stop_criteria() { return m_stop; } - const Criteria &stop_criteria() const { return m_stop; } - const Criteria &criteria() const { return m_current; } - const Status &status() const { return m_status; } + /// @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 set_strategies_iterations(const json &solver_params); + /// @brief Minimize the objective function + /// @param objFunc Objective function + /// @param x Initial guess + void minimize(Problem &objFunc, TVector &x); virtual double compute_grad_norm(const TVector &x, const TVector &grad) const; - void set_line_search(const json ¶ms); - - void minimize(Problem &objFunc, TVector &x); + // ==================================================================== + // Getters and setters + // ==================================================================== - const json &get_info() const { return solver_info; } + 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 m_status == Status::XDeltaTolerance - || m_status == Status::FDeltaTolerance - || m_status == 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 m_stop.iterations; } - size_t &max_iterations() { return 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: @@ -95,12 +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,12 +123,8 @@ namespace polysolve::nonlinear // Solver parameters // ==================================================================== - double use_grad_norm_tol; - double first_grad_norm_tol; - const double characteristic_length; - int f_delta_step_tol; // ==================================================================== // Solver state // ==================================================================== @@ -146,6 +149,7 @@ namespace polysolve::nonlinear json solver_info; + // Timers double total_time; double grad_time; double line_search_time; From 6c9181409e1c12fd67ea8c22684c4ab6caf71fb5 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Thu, 8 Feb 2024 18:14:27 -0500 Subject: [PATCH 08/13] Add our cases to operator<< --- src/polysolve/nonlinear/Criteria.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/polysolve/nonlinear/Criteria.cpp b/src/polysolve/nonlinear/Criteria.cpp index 3b30b19..de502fb 100644 --- a/src/polysolve/nonlinear/Criteria.cpp +++ b/src/polysolve/nonlinear/Criteria.cpp @@ -76,6 +76,18 @@ namespace polysolve::nonlinear 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; From 4454147aeeacded41a5d6f18d88daee90eeb3d93 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Thu, 8 Feb 2024 18:15:32 -0500 Subject: [PATCH 09/13] Update Criteria::print --- src/polysolve/nonlinear/Criteria.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/polysolve/nonlinear/Criteria.cpp b/src/polysolve/nonlinear/Criteria.cpp index de502fb..e8d56ac 100644 --- a/src/polysolve/nonlinear/Criteria.cpp +++ b/src/polysolve/nonlinear/Criteria.cpp @@ -30,6 +30,8 @@ namespace polysolve::nonlinear os << "xDelta: " << xDelta << std::endl; os << "fDelta: " << fDelta << std::endl; os << "GradNorm: " << gradNorm << std::endl; + os << "FirstGradNorm: " << firstGradNorm << std::endl; + os << "fDeltaCount: " << fDeltaCount << std::endl; } Status checkConvergence(const Criteria &stop, const Criteria ¤t) From 1d5566bc9c028f356b7fde5696989afe4a7d2e55 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Thu, 8 Feb 2024 18:41:05 -0500 Subject: [PATCH 10/13] Removed nan setting --- src/polysolve/nonlinear/Problem.hpp | 2 +- src/polysolve/nonlinear/Solver.cpp | 11 +---------- 2 files changed, 2 insertions(+), 11 deletions(-) diff --git a/src/polysolve/nonlinear/Problem.hpp b/src/polysolve/nonlinear/Problem.hpp index 8df1454..8c50113 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -13,7 +13,7 @@ namespace polysolve::nonlinear class Problem { public: - static const int Dim = Eigen::Dynamic; + static constexpr int Dim = Eigen::Dynamic; using Scalar = double; using TVector = Eigen::Matrix; using TMatrix = Eigen::Matrix; diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index 7b8b9a2..b12a788 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -286,17 +286,10 @@ namespace polysolve::nonlinear update_solver_info(objFunc(x)); objFunc.post_step(PostStepData(m_current.iterations, solver_info, x, grad)); - // Used for logging - // double xDelta = 0, gradNorm = 0; - do { m_line_search->set_is_final_strategy(m_descent_strategy == m_strategies.size() - 1); - m_current.xDelta = NaN; - m_current.fDelta = NaN; - m_current.gradNorm = NaN; - // --- Energy ------------------------------------------------------ double energy; @@ -390,9 +383,7 @@ namespace polysolve::nonlinear if (m_status != Status::Continue) break; - // --------------- - // Variable update - // --------------- + // --- Variable update --------------------------------------------- m_logger.trace( "[{}][{}] pre LS iter={:d} f={:g} ‖∇f‖={:g}", From 82038dd84f60cce8c1ef40ca03ac2aa8291a49b8 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Thu, 8 Feb 2024 18:47:44 -0500 Subject: [PATCH 11/13] Remove unused line search names --- .../nonlinear/line_search/LineSearch.cpp | 20 +++++-------------- 1 file changed, 5 insertions(+), 15 deletions(-) diff --git a/src/polysolve/nonlinear/line_search/LineSearch.cpp b/src/polysolve/nonlinear/line_search/LineSearch.cpp index 83dcb9d..f30e5a2 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.cpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.cpp @@ -20,24 +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 == "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 == "none" || name == "None") + else if (name == "None") { return std::make_shared(params, logger); } @@ -50,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) From 144a3559864e2733fd4dab4535dc43000c866664 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Thu, 8 Feb 2024 18:54:24 -0500 Subject: [PATCH 12/13] Fix not descent direction failsafe --- src/polysolve/nonlinear/Solver.cpp | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index b12a788..976115e 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -344,16 +344,18 @@ namespace polysolve::nonlinear if (m_descent_strategy >= m_strategies.size()) { 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); + 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), m_current.xDeltaDotGrad, descent_strategy_name()); - m_status = Status::NotDescentDirection; continue; } From ca0611ece291678a582c565dece2c357f179d159 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Thu, 8 Feb 2024 18:58:52 -0500 Subject: [PATCH 13/13] Add codecov.yml --- codecov.yml | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 codecov.yml diff --git a/codecov.yml b/codecov.yml new file mode 100644 index 0000000..3418f88 --- /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