From cfc409292ea9c73c28a25e0d2a9c80ed8f670858 Mon Sep 17 00:00:00 2001 From: Teseo Schneider Date: Fri, 13 Oct 2023 10:20:42 -0700 Subject: [PATCH] major cleanup with Zach and Zizhou --- src/polysolve/nonlinear/BFGS.cpp | 34 +---- src/polysolve/nonlinear/BFGS.hpp | 4 +- .../nonlinear/BoxConstraintSolver.cpp | 8 +- .../nonlinear/BoxConstraintSolver.hpp | 2 - src/polysolve/nonlinear/DenseNewton.cpp | 3 +- src/polysolve/nonlinear/DenseNewton.hpp | 2 +- src/polysolve/nonlinear/GradientDescent.cpp | 7 +- src/polysolve/nonlinear/GradientDescent.hpp | 3 +- src/polysolve/nonlinear/LBFGS.cpp | 6 +- src/polysolve/nonlinear/LBFGS.hpp | 3 +- src/polysolve/nonlinear/LBFGSB.cpp | 8 +- src/polysolve/nonlinear/LBFGSB.hpp | 3 +- src/polysolve/nonlinear/MMA.cpp | 8 +- src/polysolve/nonlinear/MMA.hpp | 3 +- src/polysolve/nonlinear/Newton.cpp | 20 +-- src/polysolve/nonlinear/Newton.hpp | 3 +- src/polysolve/nonlinear/Solver.cpp | 144 ++++++------------ src/polysolve/nonlinear/Solver.hpp | 8 +- src/polysolve/nonlinear/SparseNewton.cpp | 2 - src/polysolve/nonlinear/SparseNewton.hpp | 2 +- tests/test_nonlinear_solver.cpp | 2 - 21 files changed, 83 insertions(+), 192 deletions(-) diff --git a/src/polysolve/nonlinear/BFGS.cpp b/src/polysolve/nonlinear/BFGS.cpp index 52790c1..ad1d4bb 100644 --- a/src/polysolve/nonlinear/BFGS.cpp +++ b/src/polysolve/nonlinear/BFGS.cpp @@ -7,9 +7,9 @@ namespace polysolve::nonlinear BFGS::BFGS(const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, + const double characteristic_length, spdlog::logger &logger) - : Superclass(solver_params, dt, characteristic_length, logger) + : Superclass(solver_params, characteristic_length, logger) { linear_solver = polysolve::linear::Solver::create( linear_solver_params["solver"], linear_solver_params["precond"]); @@ -34,6 +34,8 @@ namespace polysolve::nonlinear if (this->descent_strategy == 1) this->descent_strategy++; + assert(m_prev_x.size() > 0); + reset_history(m_prev_x.size()); assert(this->descent_strategy <= 2); } @@ -54,7 +56,7 @@ namespace polysolve::nonlinear this->descent_strategy = 2; } - bool BFGS::compute_update_direction( + void BFGS::compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, @@ -79,7 +81,8 @@ namespace polysolve::nonlinear // warn if using gradient descent m_logger.warn("Unable to factorize Hessian: \"{}\"; reverting to {}", err.what(), this->descent_strategy_name()); - return compute_update_direction(objFunc, x, grad, direction); + compute_update_direction(objFunc, x, grad, direction); + return; } TVector y = grad - m_prev_grad; @@ -94,28 +97,5 @@ namespace polysolve::nonlinear m_prev_x = x; m_prev_grad = grad; - - if (std::isnan(direction.squaredNorm())) - { - reset_history(x.size()); - increase_descent_strategy(); - m_logger.log( - this->descent_strategy == 2 ? spdlog::level::warn : spdlog::level::debug, - "nan in direction {} (||∇f||={}); reverting to {}", - direction.dot(grad), this->descent_strategy_name()); - return compute_update_direction(objFunc, x, grad, direction); - } - else if (grad.squaredNorm() != 0 && direction.dot(grad) >= 0) - { - reset_history(x.size()); - increase_descent_strategy(); - m_logger.log( - this->descent_strategy == 2 ? spdlog::level::warn : spdlog::level::debug, - "BFGS direction is not a descent direction (Δx⋅g={}≥0); reverting to {}", - direction.dot(grad), this->descent_strategy_name()); - return compute_update_direction(objFunc, x, grad, direction); - } - - return true; } } // namespace polysolve::nonlinear diff --git a/src/polysolve/nonlinear/BFGS.hpp b/src/polysolve/nonlinear/BFGS.hpp index 0d0f0a7..1761a5e 100644 --- a/src/polysolve/nonlinear/BFGS.hpp +++ b/src/polysolve/nonlinear/BFGS.hpp @@ -20,7 +20,7 @@ namespace polysolve::nonlinear BFGS(const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, + const double characteristic_length, spdlog::logger &logger); std::string name() const override { return "BFGS"; } @@ -43,7 +43,7 @@ namespace polysolve::nonlinear void reset_history(const int ndof); - virtual bool compute_update_direction( + virtual void compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, diff --git a/src/polysolve/nonlinear/BoxConstraintSolver.cpp b/src/polysolve/nonlinear/BoxConstraintSolver.cpp index 3801194..feb8ba5 100644 --- a/src/polysolve/nonlinear/BoxConstraintSolver.cpp +++ b/src/polysolve/nonlinear/BoxConstraintSolver.cpp @@ -10,17 +10,16 @@ namespace polysolve::nonlinear std::unique_ptr BoxConstraintSolver::create(const std::string &solver, const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger) { if (solver == "LBFGSB" || solver == "L-BFGS-B") { - return std::make_unique(solver_params, dt, characteristic_length, logger); + return std::make_unique(solver_params, characteristic_length, logger); } else if (solver == "MMA") { - return std::make_unique(solver_params, dt, characteristic_length, logger); + return std::make_unique(solver_params, characteristic_length, logger); } throw std::runtime_error("Unrecognized solver type: " + solver); } @@ -32,10 +31,9 @@ namespace polysolve::nonlinear } BoxConstraintSolver::BoxConstraintSolver(const json &solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger) - : Superclass(solver_params, dt, characteristic_length, logger) + : Superclass(solver_params, characteristic_length, logger) { if (solver_params["max_change"].is_number()) max_change_val_ = solver_params["max_change"]; diff --git a/src/polysolve/nonlinear/BoxConstraintSolver.hpp b/src/polysolve/nonlinear/BoxConstraintSolver.hpp index e36d866..9152c88 100644 --- a/src/polysolve/nonlinear/BoxConstraintSolver.hpp +++ b/src/polysolve/nonlinear/BoxConstraintSolver.hpp @@ -16,7 +16,6 @@ namespace polysolve::nonlinear static std::unique_ptr create(const std::string &solver, const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger); // List available solvers @@ -27,7 +26,6 @@ namespace polysolve::nonlinear using typename Superclass::TVector; BoxConstraintSolver(const json &solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger); diff --git a/src/polysolve/nonlinear/DenseNewton.cpp b/src/polysolve/nonlinear/DenseNewton.cpp index a189073..bbe5e7d 100644 --- a/src/polysolve/nonlinear/DenseNewton.cpp +++ b/src/polysolve/nonlinear/DenseNewton.cpp @@ -6,11 +6,10 @@ namespace polysolve::nonlinear DenseNewton::DenseNewton( const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, + const double characteristic_length, spdlog::logger &logger) : Superclass(solver_params, linear_solver_params, - dt, characteristic_length, logger) { diff --git a/src/polysolve/nonlinear/DenseNewton.hpp b/src/polysolve/nonlinear/DenseNewton.hpp index c55ca3a..36c5791 100644 --- a/src/polysolve/nonlinear/DenseNewton.hpp +++ b/src/polysolve/nonlinear/DenseNewton.hpp @@ -16,7 +16,7 @@ namespace polysolve::nonlinear DenseNewton(const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, + const double characteristic_length, spdlog::logger &logger); std::string name() const override { return "DenseNewton"; } diff --git a/src/polysolve/nonlinear/GradientDescent.cpp b/src/polysolve/nonlinear/GradientDescent.cpp index 2983536..398ec1b 100644 --- a/src/polysolve/nonlinear/GradientDescent.cpp +++ b/src/polysolve/nonlinear/GradientDescent.cpp @@ -4,10 +4,9 @@ namespace polysolve::nonlinear { GradientDescent::GradientDescent(const json &solver_params_, - const double dt, const double characteristic_length, spdlog::logger &logger) - : Superclass(solver_params_, dt, characteristic_length, logger) + : Superclass(solver_params_, characteristic_length, logger) { } @@ -33,15 +32,13 @@ namespace polysolve::nonlinear this->descent_strategy = 1; } - bool GradientDescent::compute_update_direction( + void GradientDescent::compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, TVector &direction) { direction = -grad; - - return true; } } // namespace polysolve::nonlinear diff --git a/src/polysolve/nonlinear/GradientDescent.hpp b/src/polysolve/nonlinear/GradientDescent.hpp index 9326961..4af0845 100644 --- a/src/polysolve/nonlinear/GradientDescent.hpp +++ b/src/polysolve/nonlinear/GradientDescent.hpp @@ -15,7 +15,6 @@ namespace polysolve::nonlinear using typename Superclass::TVector; GradientDescent(const json &solver_params_, - const double dt, const double characteristic_length, spdlog::logger &logger); @@ -31,7 +30,7 @@ namespace polysolve::nonlinear protected: void reset(const int ndof) override; - virtual bool compute_update_direction( + virtual void compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, diff --git a/src/polysolve/nonlinear/LBFGS.cpp b/src/polysolve/nonlinear/LBFGS.cpp index 924aab3..71970e9 100644 --- a/src/polysolve/nonlinear/LBFGS.cpp +++ b/src/polysolve/nonlinear/LBFGS.cpp @@ -5,11 +5,9 @@ namespace polysolve::nonlinear { LBFGS::LBFGS(const json &solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger) : Superclass(solver_params, - dt, characteristic_length, logger) { @@ -49,7 +47,7 @@ namespace polysolve::nonlinear this->descent_strategy = 2; } - bool LBFGS::compute_update_direction( + void LBFGS::compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, @@ -75,7 +73,5 @@ namespace polysolve::nonlinear m_prev_x = x; m_prev_grad = grad; - - return true; } } // namespace polysolve::nonlinear diff --git a/src/polysolve/nonlinear/LBFGS.hpp b/src/polysolve/nonlinear/LBFGS.hpp index 6c87726..e1c390d 100644 --- a/src/polysolve/nonlinear/LBFGS.hpp +++ b/src/polysolve/nonlinear/LBFGS.hpp @@ -17,7 +17,6 @@ namespace polysolve::nonlinear using typename Superclass::TVector; LBFGS(const json &solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger); @@ -33,7 +32,7 @@ namespace polysolve::nonlinear void reset(const int ndof) override; - bool compute_update_direction( + void compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, diff --git a/src/polysolve/nonlinear/LBFGSB.cpp b/src/polysolve/nonlinear/LBFGSB.cpp index 243967b..a6cf46c 100644 --- a/src/polysolve/nonlinear/LBFGSB.cpp +++ b/src/polysolve/nonlinear/LBFGSB.cpp @@ -6,10 +6,9 @@ namespace polysolve::nonlinear { LBFGSB::LBFGSB(const json &solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger) - : Superclass(solver_params, dt, characteristic_length, logger) + : Superclass(solver_params, characteristic_length, logger) { m_history_size = solver_params.value("history_size", 6); } @@ -54,7 +53,7 @@ namespace polysolve::nonlinear this->descent_strategy = 1; } - bool LBFGSB::compute_update_direction( + void LBFGSB::compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, @@ -105,6 +104,7 @@ namespace polysolve::nonlinear // m_logger.debug("direc: {}", direction.transpose()); // } + // maybe remove me? if (std::isnan(direction.squaredNorm())) { reset_history(x.size()); @@ -128,7 +128,5 @@ namespace polysolve::nonlinear m_prev_x = x; m_prev_grad = grad; - - return true; } } // namespace polysolve::nonlinear diff --git a/src/polysolve/nonlinear/LBFGSB.hpp b/src/polysolve/nonlinear/LBFGSB.hpp index a450235..8127099 100644 --- a/src/polysolve/nonlinear/LBFGSB.hpp +++ b/src/polysolve/nonlinear/LBFGSB.hpp @@ -17,7 +17,6 @@ namespace polysolve::nonlinear using typename Superclass::TVector; LBFGSB(const json &solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger); @@ -50,7 +49,7 @@ namespace polysolve::nonlinear void reset_history(const int ndof); - virtual bool compute_update_direction( + virtual void compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, diff --git a/src/polysolve/nonlinear/MMA.cpp b/src/polysolve/nonlinear/MMA.cpp index 96079f5..aecdb4a 100644 --- a/src/polysolve/nonlinear/MMA.cpp +++ b/src/polysolve/nonlinear/MMA.cpp @@ -5,10 +5,9 @@ namespace polysolve::nonlinear { MMA::MMA(const json &solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger) - : Superclass(solver_params, dt, characteristic_length, logger) + : Superclass(solver_params, characteristic_length, logger) { this->m_line_search = nullptr; } @@ -36,7 +35,7 @@ namespace polysolve::nonlinear this->descent_strategy++; } - bool MMA::compute_update_direction( + void MMA::compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, @@ -67,6 +66,7 @@ namespace polysolve::nonlinear mma->Update(y.data(), grad.data(), g.data(), dg.data(), lower_bound.data(), upper_bound.data()); direction = y - x; + // maybe remove me if (std::isnan(direction.squaredNorm())) { log_and_throw_error(m_logger, "nan in direction."); @@ -76,7 +76,5 @@ namespace polysolve::nonlinear // polyfem::logger().error("Direction is not a descent direction, stop."); // throw std::runtime_error("Direction is not a descent direction, stop."); // } - - return true; } } // namespace polysolve::nonlinear \ No newline at end of file diff --git a/src/polysolve/nonlinear/MMA.hpp b/src/polysolve/nonlinear/MMA.hpp index 7822282..f97aad2 100644 --- a/src/polysolve/nonlinear/MMA.hpp +++ b/src/polysolve/nonlinear/MMA.hpp @@ -13,7 +13,6 @@ namespace polysolve::nonlinear using typename Superclass::TVector; MMA(const json &solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger); @@ -37,7 +36,7 @@ namespace polysolve::nonlinear void reset(const int ndof) override; - virtual bool compute_update_direction( + virtual void compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, diff --git a/src/polysolve/nonlinear/Newton.cpp b/src/polysolve/nonlinear/Newton.cpp index e69d62a..13740d7 100644 --- a/src/polysolve/nonlinear/Newton.cpp +++ b/src/polysolve/nonlinear/Newton.cpp @@ -6,10 +6,9 @@ namespace polysolve::nonlinear Newton::Newton( const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger) - : Superclass(solver_params, dt, characteristic_length, logger) + : Superclass(solver_params, characteristic_length, logger) { force_psd_projection = solver_params["force_psd_projection"]; } @@ -56,7 +55,7 @@ namespace polysolve::nonlinear // ======================================================================= - bool Newton::compute_update_direction( + void Newton::compute_update_direction( Problem &objFunc, const TVector &x, const TVector &grad, @@ -65,7 +64,7 @@ namespace polysolve::nonlinear if (this->descent_strategy == 2) { direction = -grad; - return true; + return; } if (this->descent_strategy == 1) @@ -75,21 +74,12 @@ namespace polysolve::nonlinear else assert(false); + // solve_linear_system will increase descent_strategy if needed const double residual = solve_linear_system(objFunc, x, grad, direction); - if (std::isnan(residual)) - // solve_linear_system will increase descent_strategy if needed - return compute_update_direction(objFunc, x, grad, direction); - if (!check_direction(residual, grad, direction)) // check_direction will increase descent_strategy if needed - return compute_update_direction(objFunc, x, grad, direction); - - // reg_weight /= reg_weight_dec; - // if (reg_weight < reg_weight_min) - // reg_weight = 0; - - return true; + compute_update_direction(objFunc, x, grad, direction); } // ======================================================================= diff --git a/src/polysolve/nonlinear/Newton.hpp b/src/polysolve/nonlinear/Newton.hpp index af547f9..6fcd03e 100644 --- a/src/polysolve/nonlinear/Newton.hpp +++ b/src/polysolve/nonlinear/Newton.hpp @@ -14,7 +14,6 @@ namespace polysolve::nonlinear Newton(const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger); @@ -25,7 +24,7 @@ namespace polysolve::nonlinear const TVector &x, const TVector &grad, TVector &direction) = 0; - bool compute_update_direction(Problem &objFunc, const TVector &x, const TVector &grad, TVector &direction) override; + void compute_update_direction(Problem &objFunc, const TVector &x, const TVector &grad, TVector &direction) override; bool check_direction(const double residual, const TVector &grad, const TVector &direction); // ==================================================================== diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index a2e791d..3e2ac2b 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -22,29 +22,28 @@ namespace polysolve::nonlinear std::unique_ptr Solver::create(const std::string &solver, const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger) { if (solver == "BFGS") { - return std::make_unique(solver_params, linear_solver_params, dt, characteristic_length, logger); + return std::make_unique(solver_params, linear_solver_params, characteristic_length, logger); } else if (solver == "DenseNewton" || solver == "dense_newton") { - return std::make_unique(solver_params, linear_solver_params, dt, characteristic_length, logger); + return std::make_unique(solver_params, linear_solver_params, characteristic_length, logger); } else if (solver == "SparseNewton" || solver == "sparse_newton") { - return std::make_unique(solver_params, linear_solver_params, dt, characteristic_length, logger); + return std::make_unique(solver_params, linear_solver_params, characteristic_length, logger); } else if (solver == "GradientDescent" || solver == "gradient_descent") { - return std::make_unique(solver_params, dt, characteristic_length, logger); + return std::make_unique(solver_params, characteristic_length, logger); } else if (solver == "LBFGS" || solver == "L-BFGS") { - return std::make_unique(solver_params, dt, characteristic_length, logger); + return std::make_unique(solver_params, characteristic_length, logger); } throw std::runtime_error("Unrecognized solver type: " + solver); } @@ -59,10 +58,9 @@ namespace polysolve::nonlinear } Solver::Solver(const json &solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger) - : dt(dt), m_logger(logger), characteristic_length(characteristic_length) + : m_logger(logger), characteristic_length(characteristic_length) { TCriteria criteria = TCriteria::defaults(); criteria.xDelta = solver_params["x_delta"]; @@ -77,7 +75,6 @@ namespace polysolve::nonlinear // criteria.condition = solver_params["condition"]; this->setStopCriteria(criteria); - normalize_gradient = solver_params["relative_gradient"]; use_grad_norm_tol = solver_params["line_search"]["use_grad_norm_tol"]; first_grad_norm_tol = solver_params["first_grad_norm_tol"]; @@ -106,59 +103,26 @@ namespace polysolve::nonlinear // --------------------------- // Initialize the minimization // --------------------------- - 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()); - // double factor = 1e-5; - // Set these to nan to indicate they have not been computed yet double old_energy = NaN; - { POLYSOLVE_SCOPED_STOPWATCH("constraint set update", constraint_set_update_time, m_logger); objFunc.solution_changed(x); } - { - POLYSOLVE_SCOPED_STOPWATCH("compute gradient", grad_time, m_logger); - objFunc.gradient(x, grad); - } - double first_grad_norm = compute_grad_norm(x, grad); - if (std::isnan(first_grad_norm)) - { - this->m_status = cppoptlib::Status::UserDefined; - m_error_code = ErrorCode::NAN_ENCOUNTERED; - log_and_throw_error(m_logger, "[{}] Initial gradient is nan; stopping", name()); - return; - } - this->m_current.xDelta = NaN; // we don't know the initial step size - this->m_current.fDelta = old_energy; - this->m_current.gradNorm = first_grad_norm / (normalize_gradient ? first_grad_norm : 1); - - const auto current_g_norm = this->m_stop.gradNorm; + const auto g_norm_tol = this->m_stop.gradNorm; this->m_stop.gradNorm = first_grad_norm_tol; - this->m_status = checkConvergence(this->m_stop, this->m_current); - if (this->m_status != cppoptlib::Status::Continue) - { - POLYSOLVE_SCOPED_STOPWATCH("compute objective function", obj_fun_time, m_logger); - this->m_current.fDelta = objFunc.value(x); - m_logger.info( - "[{}] Not even starting, {} (f={:g} ‖∇f‖={:g} g={:g} tol={:g})", - name(), this->m_status, this->m_current.fDelta, first_grad_norm, this->m_current.gradNorm, this->m_stop.gradNorm); - update_solver_info(this->m_current.fDelta); - return; - } - this->m_stop.gradNorm = current_g_norm; StopWatch stop_watch("non-linear solver", this->total_time, m_logger); stop_watch.start(); - if (m_line_search) - m_line_search->use_grad_norm_tol = use_grad_norm_tol; - objFunc.save_to_file(x); m_logger.debug( @@ -171,11 +135,17 @@ namespace polysolve::nonlinear do { + this->m_current.xDelta = NaN; + this->m_current.fDelta = NaN; + this->m_current.gradNorm = NaN; + + //////////// Energy double energy; { POLYSOLVE_SCOPED_STOPWATCH("compute objective function", obj_fun_time, m_logger); energy = objFunc.value(x); } + if (!std::isfinite(energy)) { this->m_status = cppoptlib::Status::UserDefined; @@ -184,6 +154,13 @@ namespace polysolve::nonlinear break; } + this->m_current.fDelta = std::abs(old_energy - energy); // / std::abs(old_energy); + old_energy = energy; + this->m_status = checkConvergence(this->m_stop, this->m_current); + if (this->m_status != cppoptlib::Status::Continue) + break; + + ///////////// gradient { POLYSOLVE_SCOPED_STOPWATCH("compute gradient", grad_time, m_logger); objFunc.gradient(x, grad); @@ -197,17 +174,17 @@ namespace polysolve::nonlinear log_and_throw_error(m_logger, "[{}] Gradient is nan; stopping", name()); break; } + this->m_current.gradNorm = grad_norm; + this->m_status = checkConvergence(this->m_stop, this->m_current); + if (this->m_status != cppoptlib::Status::Continue) + break; // ------------------------ // Compute update direction // ------------------------ - // Compute a Δx to update the variable - if (!compute_update_direction(objFunc, x, grad, delta_x)) - { - this->m_status = cppoptlib::Status::Continue; - continue; - } + // + compute_update_direction(objFunc, x, grad, delta_x); if (is_direction_descent() && grad_norm != 0 && delta_x.dot(grad) >= 0) { @@ -230,30 +207,34 @@ namespace polysolve::nonlinear } // Use the maximum absolute displacement value divided by the timestep, - // so the units are in velocity units. - // TODO: Also divide by the world scale to make this criteria scale invariant. - this->m_current.xDelta = descent_strategy == 2 ? NaN : (delta_x_norm / dt); - this->m_current.fDelta = std::abs(old_energy - energy); // / std::abs(old_energy); - // if normalize_gradient, use relative to first norm - this->m_current.gradNorm = grad_norm / (normalize_gradient ? first_grad_norm : 1); - + this->m_current.xDelta = descent_strategy == 2 ? NaN : delta_x_norm; this->m_status = checkConvergence(this->m_stop, this->m_current); - - old_energy = energy; + if (this->m_status != cppoptlib::Status::Continue) + break; // --------------- // Variable update // --------------- // Perform a line_search to compute step scale - double rate = line_search(x, delta_x, objFunc); + double rate = m_line_search->line_search(x, delta_x, objFunc); if (std::isnan(rate)) { - // descent_strategy set by line_search upon failure - if (this->m_status == cppoptlib::Status::Continue) + assert(this->m_status == cppoptlib::Status::Continue); + + if (descent_strategy < 2) // 2 is the max, grad descent + { + increase_descent_strategy(); + m_logger.warn( + "[{}] Line search failed; reverting to {}", name(), descent_strategy_name()); continue; + } else - break; + { + assert(descent_strategy == 2); // failed on gradient descent + this->m_status = cppoptlib::Status::UserDefined; // Line search failed on gradient descent, so quit! + log_and_throw_error(m_logger, "[{}] Line search failed on gradient descent; stopping", name()); + } } x += rate * delta_x; @@ -287,6 +268,9 @@ namespace polysolve::nonlinear objFunc.save_to_file(x); + // 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)); stop_watch.stop(); @@ -297,8 +281,6 @@ namespace polysolve::nonlinear if (!allow_out_of_iterations && this->m_status == cppoptlib::Status::IterationLimit) log_and_throw_error(m_logger, "[{}] Reached iteration limit (limit={})", name(), this->m_stop.iterations); - if (this->m_current.iterations == 0) - log_and_throw_error(m_logger, "[{}] Unable to take a step", name()); if (this->m_status == cppoptlib::Status::UserDefined && m_error_code != ErrorCode::SUCCESS) log_and_throw_error(m_logger, "[{}] Failed to find minimizer", name()); @@ -312,35 +294,6 @@ namespace polysolve::nonlinear update_solver_info(objFunc.value(x)); } - double Solver::line_search(const TVector &x, const TVector &delta_x, Problem &objFunc) - { - POLYSOLVE_SCOPED_STOPWATCH("line search", line_search_time, m_logger); - - if (!m_line_search) - { - objFunc.solution_changed(x + delta_x); - return 1; // no linesearch - } - - double rate = m_line_search->line_search(x, delta_x, objFunc); - - if (std::isnan(rate) && descent_strategy < 2) // 2 is the max, grad descent - { - increase_descent_strategy(); - m_logger.warn( - "[{}] Line search failed; reverting to {}", name(), descent_strategy_name()); - this->m_status = cppoptlib::Status::Continue; // Try the step again with gradient descent - } - else if (std::isnan(rate)) - { - assert(descent_strategy == 2); // failed on gradient descent - this->m_status = cppoptlib::Status::UserDefined; // Line search failed on gradient descent, so quit! - log_and_throw_error(m_logger, "[{}] Line search failed on gradient descent; stopping", name()); - } - - return rate; - } - void Solver::reset(const int ndof) { this->m_current.reset(); @@ -381,7 +334,6 @@ namespace polysolve::nonlinear solver_info["fDelta"] = crit.fDelta; solver_info["gradNorm"] = crit.gradNorm; solver_info["condition"] = crit.condition; - solver_info["relative_gradient"] = normalize_gradient; double per_iteration = crit.iterations ? crit.iterations : 1; diff --git a/src/polysolve/nonlinear/Solver.hpp b/src/polysolve/nonlinear/Solver.hpp index 67f8d6b..14ee1f8 100644 --- a/src/polysolve/nonlinear/Solver.hpp +++ b/src/polysolve/nonlinear/Solver.hpp @@ -31,7 +31,6 @@ namespace polysolve::nonlinear static std::unique_ptr create(const std::string &solver, const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger); // List available solvers @@ -47,7 +46,6 @@ namespace polysolve::nonlinear /// @param dt time step size (use 1 if not time-dependent) TODO /// @param logger Solver(const json &solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger); @@ -59,8 +57,6 @@ namespace polysolve::nonlinear void minimize(Problem &objFunc, TVector &x) override; - double line_search(const TVector &x, const TVector &delta_x, Problem &objFunc); - const json &get_info() const { return solver_info; } ErrorCode error_code() const { return m_error_code; } @@ -84,10 +80,8 @@ namespace polysolve::nonlinear // Solver parameters // ==================================================================== - bool normalize_gradient; double use_grad_norm_tol; double first_grad_norm_tol; - double dt; const double characteristic_length; @@ -99,7 +93,7 @@ namespace polysolve::nonlinear virtual void reset(const int ndof); // Compute the search/update direction - virtual bool compute_update_direction(Problem &objFunc, const TVector &x_vec, const TVector &grad, TVector &direction) = 0; + virtual void compute_update_direction(Problem &objFunc, const TVector &x_vec, const TVector &grad, TVector &direction) = 0; virtual void set_default_descent_strategy() = 0; virtual void increase_descent_strategy() = 0; diff --git a/src/polysolve/nonlinear/SparseNewton.cpp b/src/polysolve/nonlinear/SparseNewton.cpp index 8a978ce..5e6fc9a 100644 --- a/src/polysolve/nonlinear/SparseNewton.cpp +++ b/src/polysolve/nonlinear/SparseNewton.cpp @@ -8,12 +8,10 @@ namespace polysolve::nonlinear SparseNewton::SparseNewton( const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, spdlog::logger &logger) : Superclass(solver_params, linear_solver_params, - dt, characteristic_length, logger) { diff --git a/src/polysolve/nonlinear/SparseNewton.hpp b/src/polysolve/nonlinear/SparseNewton.hpp index 034201f..c004ef1 100644 --- a/src/polysolve/nonlinear/SparseNewton.hpp +++ b/src/polysolve/nonlinear/SparseNewton.hpp @@ -16,7 +16,7 @@ namespace polysolve::nonlinear SparseNewton(const json &solver_params, const json &linear_solver_params, - const double dt, const double characteristic_length, + const double characteristic_length, spdlog::logger &logger); std::string name() const override { return "SparseNewton"; } diff --git a/tests/test_nonlinear_solver.cpp b/tests/test_nonlinear_solver.cpp index 1904e9f..2d272de 100644 --- a/tests/test_nonlinear_solver.cpp +++ b/tests/test_nonlinear_solver.cpp @@ -258,7 +258,6 @@ TEST_CASE("non-linear", "[solver]") linear_solver_params["solver"] = "Eigen::SimplicialLDLT"; linear_solver_params["precond"] = polysolve::linear::Solver::default_precond(); - const double dt = 0.1; const double characteristic_length = 1; static std::shared_ptr logger = spdlog::stdout_color_mt("test_logger"); @@ -286,7 +285,6 @@ TEST_CASE("non-linear", "[solver]") auto solver = Solver::create(solver_name, solver_params, linear_solver_params, - dt, characteristic_length, *logger); try