Drop cppoptlib

CMakeLists.txt

@@ -311,17 +311,13 @@ target_link_libraries(polysolve_linear PRIVATE polysolve::warnings)
 # polysolve::linear
 target_link_libraries(polysolve PUBLIC polysolve::linear)
 
-# CppNumericalSolvers
-include(cppoptlib)
-target_link_libraries(polysolve PUBLIC cppoptlib)
-
 # LBFGSpp
 include(LBFGSpp)
-target_link_libraries(polysolve PUBLIC LBFGSpp::LBFGSpp)
+target_link_libraries(polysolve PRIVATE LBFGSpp::LBFGSpp)
 
 # finite-diff (include this after eigen)
 include(finite-diff)
-target_link_libraries(polysolve PUBLIC finitediff::finitediff)
+target_link_libraries(polysolve PRIVATE finitediff::finitediff)
 
 # Sanitizers
 if(POLYSOLVE_WITH_SANITIZERS)

codecov.yml

@@ -0,0 +1,11 @@
+coverage:
+  status:
+    project:
+      default:
+        target: auto
+        threshold: 1%
+    patch:
+      default:
+        target: 75%
+        threshold: 5%
+        only_pulls: true

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"

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})

src/polysolve/nonlinear/Criteria.cpp

@@ -0,0 +1,112 @@
+// Source: https://github.com/PatWie/CppNumericalSolvers/blob/7eddf28fa5a8872a956d3c8666055cac2f5a535d/include/cppoptlib/meta.h
+// License: MIT
+#include "Criteria.hpp"
+
+namespace polysolve::nonlinear
+{
+    bool is_converged_status(const Status s)
+    {
+        return s == Status::XDeltaTolerance || s == Status::FDeltaTolerance || s == Status::GradNormTolerance;
+    }
+
+    Criteria::Criteria()
+    {
+        reset();
+    }
+
+    void Criteria::reset()
+    {
+        iterations = 0;
+        xDelta = 0;
+        fDelta = 0;
+        gradNorm = 0;
+        firstGradNorm = 0;
+        fDeltaCount = 0;
+        xDeltaDotGrad = 0;
+    }
+
+    void Criteria::print(std::ostream &os) const
+    {
+        os << "Iterations: " << iterations << std::endl;
+        os << "xDelta:     " << xDelta << std::endl;
+        os << "fDelta:     " << fDelta << std::endl;
+        os << "GradNorm:   " << gradNorm << std::endl;
+        os << "xDeltaDotGrad: " << xDeltaDotGrad << std::endl;
+        os << "FirstGradNorm: " << firstGradNorm << std::endl;
+        os << "fDeltaCount: " << fDeltaCount << std::endl;
+    }
+
+    Status checkConvergence(const Criteria &stop, const Criteria &current)
+    {
+        if (stop.iterations > 0 && current.iterations > stop.iterations)
+        {
+            return Status::IterationLimit;
+        }
+        if (stop.xDelta > 0 && current.xDelta < stop.xDelta)
+        {
+            return Status::XDeltaTolerance;
+        }
+        if (stop.fDelta > 0 && current.fDelta < stop.fDelta && current.fDeltaCount >= stop.fDeltaCount)
+        {
+            return Status::FDeltaTolerance;
+        }
+        const double stopGradNorm = current.iterations == 0 ? stop.firstGradNorm : stop.gradNorm;
+        if (stopGradNorm > 0 && current.gradNorm < stopGradNorm)
+        {
+            return Status::GradNormTolerance;
+        }
+        // Δx⋅∇f ≥ 0 means that the search direction is not a descent direction
+        if (stop.xDeltaDotGrad < 0 && current.xDeltaDotGrad > stop.xDeltaDotGrad)
+        {
+            return Status::NotDescentDirection;
+        }
+        return Status::Continue;
+    }
+
+    std::ostream &operator<<(std::ostream &os, const Status &s)
+    {
+        switch (s)
+        {
+        case Status::NotStarted:
+            os << "Solver not started.";
+            break;
+        case Status::Continue:
+            os << "Convergence criteria not reached.";
+            break;
+        case Status::IterationLimit:
+            os << "Iteration limit reached.";
+            break;
+        case Status::XDeltaTolerance:
+            os << "Change in parameter vector too small.";
+            break;
+        case Status::FDeltaTolerance:
+            os << "Change in cost function value too small.";
+            break;
+        case Status::GradNormTolerance:
+            os << "Gradient vector norm too small.";
+            break;
+        case Status::ObjectiveCustomStop:
+            os << "Objective function specified to stop.";
+            break;
+        case Status::NanEncountered:
+            os << "Objective or gradient function returned NaN.";
+            break;
+        case Status::NotDescentDirection:
+            os << "Search direction not a descent direction.";
+            break;
+        case Status::LineSearchFailed:
+            os << "Line search failed.";
+            break;
+        default:
+            os << "Unknown status.";
+            break;
+        }
+        return os;
+    }
+
+    std::ostream &operator<<(std::ostream &os, const Criteria &c)
+    {
+        c.print(os);
+        return os;
+    }
+} // namespace polysolve::nonlinear

src/polysolve/nonlinear/Criteria.hpp

@@ -0,0 +1,52 @@
+#pragma once
+
+#include <cstddef>
+#include <iostream>
+
+namespace polysolve::nonlinear
+{
+    // Source: https://github.com/PatWie/CppNumericalSolvers/blob/7eddf28fa5a8872a956d3c8666055cac2f5a535d/include/cppoptlib/meta.h
+    // License: MIT
+
+    enum class Status
+    {
+        NotStarted = -1, ///< The solver has not been started
+        Continue = 0,    ///< The solver should continue
+        // Success cases
+        IterationLimit,      ///< The maximum number of iterations has been reached
+        XDeltaTolerance,     ///< The change in the parameter vector is below the tolerance
+        FDeltaTolerance,     ///< The change in the cost function is below the tolerance
+        GradNormTolerance,   ///< The norm of the gradient vector is below the tolerance
+        ObjectiveCustomStop, ///< The objective function specified to stop
+        // Failure cases
+        NanEncountered,      ///< The objective function returned NaN
+        NotDescentDirection, ///< The search direction is not a descent direction
+        LineSearchFailed,    ///< The line search failed
+    };
+
+    bool is_converged_status(const Status s);
+
+    class Criteria
+    {
+    public:
+        size_t iterations;    ///< Maximum number of iterations
+        double xDelta;        ///< Minimum change in parameter vector
+        double fDelta;        ///< Minimum change in cost function
+        double gradNorm;      ///< Minimum norm of gradient vector
+        double firstGradNorm; ///< Initial norm of gradient vector
+        double xDeltaDotGrad; ///< Dot product of parameter vector and gradient vector
+        unsigned fDeltaCount; ///< Number of steps where fDelta is satisfied
+
+        Criteria();
+
+        void reset();
+
+        void print(std::ostream &os) const;
+    };
+
+    Status checkConvergence(const Criteria &stop, const Criteria &current);
+
+    std::ostream &operator<<(std::ostream &os, const Status &s);
+
+    std::ostream &operator<<(std::ostream &os, const Criteria &c);
+} // namespace polysolve::nonlinear

src/polysolve/nonlinear/Problem.hpp

@@ -2,47 +2,115 @@
 
 #include <polysolve/Types.hpp>
 
+#include "Criteria.hpp"
 #include "PostStepData.hpp"
 
-#include <cppoptlib/problem.h>
-
 #include <memory>
 #include <vector>
 
 namespace polysolve::nonlinear
 {
-    class Problem : public cppoptlib::Problem<double>
+    class Problem
     {
     public:
-        using typename cppoptlib::Problem<double>::Scalar;
-        using typename cppoptlib::Problem<double>::TVector;
-        typedef polysolve::StiffnessMatrix THessian;
-
-        // disable warning for dense hessian
-        using cppoptlib::Problem<double>::hessian;
+        static constexpr int Dim = Eigen::Dynamic;
+        using Scalar = double;
+        using TVector = Eigen::Matrix<Scalar, Dim, 1>;
+        using TMatrix = Eigen::Matrix<Scalar, Dim, Dim>;
+        using THessian = StiffnessMatrix;
 
+    public:
         Problem() {}
-        ~Problem() = default;
+        virtual ~Problem() = default;
 
+        /// @brief Initialize the problem.
+        /// @param x0 Initial guess.
         virtual void init(const TVector &x0) {}
 
-        virtual double value(const TVector &x) override = 0;
-        virtual void gradient(const TVector &x, TVector &gradv) override = 0;
+        /// @brief Compute the value of the function at x.
+        /// @param x Degrees of freedom.
+        /// @return The value of the function at x.
+        Scalar operator()(const TVector &x) { return value(x); }
+
+        /// @brief Compute the value of the function at x.
+        /// @param x Degrees of freedom.
+        /// @return The value of the function at x.
+        virtual Scalar value(const TVector &x) = 0;
+
+        /// @brief Compute the gradient of the function at x.
+        /// @param[in] x Degrees of freedom.
+        /// @param[out] grad Gradient of the function at x.
+        virtual void gradient(const TVector &x, TVector &grad) = 0;
+
+        /// @brief Compute the Hessian of the function at x.
+        /// @param[in] x Degrees of freedom.
+        /// @param[out] hessian Hessian of the function at x.
+        virtual void hessian(const TVector &x, TMatrix &hessian)
+        {
+            throw std::runtime_error("Dense Hessian not implemented.");
+        }
+
+        /// @brief Compute the Hessian of the function at x.
+        /// @param[in] x Degrees of freedom.
+        /// @param[out] hessian Hessian of the function at x.
         virtual void hessian(const TVector &x, THessian &hessian) = 0;
 
+        /// @brief Determine if the step from x0 to x1 is valid.
+        /// @param x0 Starting point.
+        /// @param x1 Ending point.
+        /// @return True if the step is valid, false otherwise.
         virtual bool is_step_valid(const TVector &x0, const TVector &x1) { return true; }
+
+        /// @brief Determine a maximum step size from x0 to x1.
+        /// @param x0 Starting point.
+        /// @param x1 Ending point.
+        /// @return Maximum step size.
         virtual double max_step_size(const TVector &x0, const TVector &x1) { return 1; }
 
+        // --- Callbacks ------------------------------------------------------
+
+        /// @brief Callback function for the start of a line search.
+        /// @param x0 Starting point.
+        /// @param x1 Ending point.
         virtual void line_search_begin(const TVector &x0, const TVector &x1) {}
+
+        /// @brief Callback function for the end of a line search.
         virtual void line_search_end() {}
+
+        /// @brief Callback function for the end of a step.
+        /// @param data Post step data.
         virtual void post_step(const PostStepData &data) {}
 
+        /// @brief Set the project to PSD flag.
+        /// @param val True if the problem should be projected to PSD, false otherwise.
         virtual void set_project_to_psd(bool val) {}
 
+        /// @brief Callback function for when the solution changes.
+        /// @param new_x New solution.
         virtual void solution_changed(const TVector &new_x) {}
 
+        /// @brief Callback function used to determine if the solver should stop.
+        /// @param state Current state of the solver.
+        /// @param x Current solution.
+        /// @return True if the solver should stop, false otherwise.
+        virtual bool callback(const Criteria &state, const TVector &x) { return true; }
+
+        /// @brief Callback function used Determine if the solver should stop.
+        /// @param x Current solution.
+        /// @return True if the solver should stop, false otherwise.
         virtual bool stop(const TVector &x) { return false; }
 
+        /// --- Misc ----------------------------------------------------------
+
+        /// @brief Sample the function along a direction.
+        /// @param[in] x Starting point.
+        /// @param[in] direction Direction to sample along.
+        /// @param[in] start Starting step size.
+        /// @param[in] end Ending step size.
+        /// @param[in] num_samples Number of samples to take.
+        /// @param[out] alphas Sampled step sizes.
+        /// @param[out] fs Sampled function values.
+        /// @param[out] valid If each sample is valid.
         void sample_along_direction(
             const Problem::TVector &x,
             const Problem::TVector &direction,