put back f_delta and x_delta (#49)

* put back f_delta and x_delta

* revert change in line search

* more info in log

* address comments

* move f_delta to advanced

* fix compile

non-linear-solver-spec.json

@@ -15,7 +15,8 @@
             "LBFGS",
             "LBFGSB",
             "Newton",
-            "box_constraints"
+            "box_constraints",
+            "advanced"
         ],
         "doc": "Settings for nonlinear solver. Interior-loop linear solver settings are defined in the solver/linear section."
     },
@@ -295,5 +296,28 @@
         "pointer": "/box_constraints/max_change/*",
         "type": "float",
         "doc": "Maximum change of every optimization variable in one iteration, only for solvers with box constraints."
+    },
+    {
+        "pointer": "/advanced",
+        "default": null,
+        "type": "object",
+        "optional": [
+            "f_delta",
+            "f_delta_step_tol"
+        ],
+        "doc": "Nonlinear solver advanced options"
+    },
+    {
+        "pointer": "/advanced/f_delta",
+        "default": 0,
+        "min": 0,
+        "type": "float",
+        "doc": "Dangerous Option: Quit the optimization if the solver reduces the energy by less than f_delta for consecutive f_delta_step_tol steps."
+    },
+    {
+        "pointer": "/advanced/f_delta_step_tol",
+        "default": 100,
+        "type": "int",
+        "doc": "Dangerous Option: Quit the optimization if the solver reduces the energy by less than f_delta for consecutive f_delta_step_tol steps."
     }
 ]

src/polysolve/nonlinear/Solver.cpp

@@ -107,7 +107,7 @@ namespace polysolve::nonlinear
     {
         TCriteria criteria = TCriteria::defaults();
         criteria.xDelta = solver_params["x_delta"];
-        criteria.fDelta = -1;
+        criteria.fDelta = solver_params["advanced"]["f_delta"];
         criteria.gradNorm = solver_params["grad_norm"];
 
         criteria.xDelta *= characteristic_length;
@@ -125,6 +125,8 @@ namespace polysolve::nonlinear
         use_grad_norm_tol *= characteristic_length;
         first_grad_norm_tol *= characteristic_length;
 
+        f_delta_step_tol = solver_params["advanced"]["f_delta_step_tol"];
+
         set_line_search(solver_params);
     }
 
@@ -194,6 +196,9 @@ namespace polysolve::nonlinear
 
         update_solver_info(objFunc.value(x));
 
+        int f_delta_step_cnt = 0;
+        double f_delta = 0;
+
         do
         {
             m_line_search->set_is_final_strategy(m_descent_strategy == m_strategies.size() - 1);
@@ -217,7 +222,9 @@ namespace polysolve::nonlinear
                 break;
             }
 
-            this->m_current.fDelta = std::abs(old_energy - energy);
+            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;
 
             ///////////// gradient
             {
@@ -285,7 +292,7 @@ namespace polysolve::nonlinear
             }
 
             // Use the maximum absolute displacement value divided by the timestep,
-            this->m_current.xDelta = (m_descent_strategy == m_strategies.size() - 1) ? NaN : delta_x_norm;
+            this->m_current.xDelta = delta_x_norm;
             this->m_status = checkConvergence(this->m_stop, this->m_current);
             if (this->m_status != cppoptlib::Status::Continue)
                 break;
@@ -318,12 +325,18 @@ namespace polysolve::nonlinear
             // if the strategy got changed, we start counting
             if (m_descent_strategy != previous_strategy)
                 current_strategy_iter = 0;
-            // if we did enoug lower strategy, we revert back to normal
+            // if we did enough lower strategy, we revert back to normal
             if (m_descent_strategy != 0 && current_strategy_iter >= m_iter_per_strategy[m_descent_strategy])
             {
+                const std::string prev_strategy_name = descent_strategy_name();
+
                 m_descent_strategy = 0;
                 for (auto &s : m_strategies)
                     s->reset(x.size());
+
+                m_logger.debug(
+                    "[{}][{}] {} was successful for {} iterations; resetting to {}",
+                    name(), m_line_search->name(), prev_strategy_name, current_strategy_iter, descent_strategy_name());
             }
 
             previous_strategy = m_descent_strategy;
@@ -343,11 +356,18 @@ namespace polysolve::nonlinear
 
             objFunc.post_step(this->m_current.iterations, x);
 
+            if (f_delta < this->m_stop.fDelta)
+                f_delta_step_cnt++;
+            else
+                f_delta_step_cnt = 0;
+
             m_logger.debug(
-                "[{}][{}] iter={:d} f={:g} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g} Δx⋅∇f(x)={:g} rate={:g} ‖step‖={:g}",
+                "[{}][{}] 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})",
                 name(), m_line_search->name(),
-                this->m_current.iterations, energy, this->m_current.fDelta,
-                this->m_current.gradNorm, this->m_current.xDelta, delta_x.dot(grad), rate, step);
+                this->m_current.iterations, energy, f_delta,
+                this->m_current.gradNorm, this->m_current.xDelta, delta_x.dot(grad), rate, step,
+                this->m_stop.iterations, this->m_stop.fDelta, this->m_stop.gradNorm, this->m_stop.xDelta);
 
             if (++this->m_current.iterations >= this->m_stop.iterations)
                 this->m_status = cppoptlib::Status::IterationLimit;
@@ -370,11 +390,14 @@ namespace polysolve::nonlinear
             log_and_throw_error(m_logger, "[{}][{}] Failed to find minimizer", name(), m_line_search->name());
 
         double tot_time = stop_watch.getElapsedTimeInSec();
-        m_logger.info(
-            "[{}][{}] Finished: {} Took {:g}s (niters={:d} f={:g} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g} ftol={})",
+        const bool succeeded = this->m_status == cppoptlib::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})",
             name(), m_line_search->name(),
             this->m_status, tot_time, this->m_current.iterations,
-            old_energy, this->m_current.fDelta, this->m_current.gradNorm, this->m_current.xDelta, this->m_stop.fDelta);
+            old_energy, f_delta, this->m_current.gradNorm, this->m_current.xDelta,
+            this->m_stop.iterations, this->m_stop.fDelta, this->m_stop.gradNorm, this->m_stop.xDelta);
 
         log_times();
         update_solver_info(objFunc.value(x));

src/polysolve/nonlinear/Solver.hpp

@@ -111,6 +111,7 @@ namespace polysolve::nonlinear
 
         const double characteristic_length;
 
+        int f_delta_step_tol;
         // ====================================================================
         //                           Solver state
         // ====================================================================