Merge pull request #71 from maxpaik16/add-comments

Add comments to nonlinear solver code

src/polysolve/linear/Solver.hpp

@@ -51,22 +51,22 @@ namespace polysolve::linear
         /// Selects the correct solver based on params using the fallback or the list of solvers if necessary
         static void select_valid_solver(json &params, spdlog::logger &logger);
 
-        // Static constructor
-        //
-        // @param[in]  params   Parameter of the solver, including name and preconditioner
-        // @param[in]  logger   Logger used for error
-        // @param[in]  strict_validation    strict validation of the input paraams
-        // @param[out] a pointer to a linear solver
+        /// @brief Static constructor
+        ///
+        /// @param[in]  params   Parameter of the solver, including name and preconditioner
+        /// @param[in]  logger   Logger used for error
+        /// @param[in]  strict_validation    strict validation of the input paraams
+        /// @return a pointer to a linear solver
         //
         static std::unique_ptr<Solver> create(const json &params,
                                               spdlog::logger &logger,
                                               const bool strict_validation = true);
 
-        // Static constructor
-        //
-        // @param[in]  solver   Solver type
-        // @param[in]  precond  Preconditioner for iterative solvers
-        //
+        /// @brief Static constructor
+        ///
+        /// @param[in]  solver   Solver type
+        /// @param[in]  precond  Preconditioner for iterative solvers
+        ///
         static std::unique_ptr<Solver> create(const std::string &solver, const std::string &precond);
 
         // List available solvers
@@ -86,50 +86,45 @@ namespace polysolve::linear
         // Public interface //
         //////////////////////
 
-        // Set solver parameters
+        /// Set solver parameters
         virtual void set_parameters(const json &params) {}
 
-        // Get info on the last solve step
+        /// Get info on the last solve step
         virtual void get_info(json &params) const {};
 
-        // Analyze sparsity pattern
+        /// Analyze sparsity pattern
         virtual void analyze_pattern(const StiffnessMatrix &A, const int precond_num) {}
 
-        // Factorize system matrix
+        /// Factorize system matrix
         virtual void factorize(const StiffnessMatrix &A) {}
 
-        // Analyze sparsity pattern of a dense matrix
+        /// Analyze sparsity pattern of a dense matrix
         virtual void analyze_pattern_dense(const Eigen::MatrixXd &A, const int precond_num) {}
 
-        // Factorize system matrix of a dense matrix
+        /// Factorize system matrix of a dense matrix
         virtual void factorize_dense(const Eigen::MatrixXd &A) {}
 
-        // If solver uses dense matrices
+        /// If solver uses dense matrices
         virtual bool is_dense() const { return false; }
 
-        // Set block size for multigrid solvers
+        /// Set block size for multigrid solvers
         virtual void set_block_size(int block_size) {}
 
-        // If the problem is nullspace for multigrid solvers
+        /// If the problem is nullspace for multigrid solvers
         virtual void set_is_nullspace(const VectorXd &x) {}
 
-        //
-        // @brief         { Solve the linear system Ax = b }
-        //
-        // @param[in]     b     { Right-hand side. }
-        // @param[in,out] x     { Unknown to compute. When using an iterative
-        //                      solver, the input unknown vector is used as an
-        //                      initial guess, and must thus be properly allocated
-        //                      and initialized. }
-        //
+        ///
+        /// @brief         { Solve the linear system Ax = b }
+        ///
+        /// @param[in]     b     { Right-hand side. }
+        /// @param[in,out] x     { Unknown to compute. When using an iterative
+        ///                      solver, the input unknown vector is used as an
+        ///                      initial guess, and must thus be properly allocated
+        ///                      and initialized. }
+        ///
         virtual void solve(const Ref<const VectorXd> b, Ref<VectorXd> x) = 0;
 
-    public:
-        ///////////
-        // Debug //
-        ///////////
-
-        // Name of the solver type (for debugging purposes)
+        /// @brief Name of the solver type (for debugging purposes)
         virtual std::string name() const { return ""; }
     };
 

src/polysolve/nonlinear/Problem.hpp

@@ -10,6 +10,7 @@
 
 namespace polysolve::nonlinear
 {
+    /// @brief Class defining optimization problem to be solved. To be defined by user code
     class Problem
     {
     public:

src/polysolve/nonlinear/Solver.hpp

@@ -41,7 +41,7 @@ namespace polysolve::nonlinear
             spdlog::logger &logger,
             const bool strict_validation = true);
 
-        // List available solvers
+        /// @brief List available solvers
         static std::vector<std::string> available_solvers();
 
     public:
@@ -82,10 +82,17 @@ namespace polysolve::nonlinear
         /// @brief If true the solver will not throw an error if the maximum number of iterations is reached
         bool allow_out_of_iterations = false;
 
+
         /// @brief Get the line search object
         const std::shared_ptr<line_search::LineSearch> &line_search() const { return m_line_search; };
 
     protected:
+        /// @brief Compute direction in which the argument should be updated 
+        /// @param objFunc Problem to be minimized
+        /// @param x Current input (n x 1)
+        /// @param grad Gradient at current step (n x 1)
+        /// @param[out] direction Current update direction (n x 1)
+        /// @return True if update direction was found, False otherwises
         virtual bool compute_update_direction(
             Problem &objFunc,
             const TVector &x,
@@ -116,6 +123,10 @@ namespace polysolve::nonlinear
 
         FiniteDiffStrategy gradient_fd_strategy = FiniteDiffStrategy::NONE;
         double gradient_fd_eps = 1e-7;
+        /// @brief Check gradient versus finite difference results
+        /// @param objFunc Problem defining relevant objective function
+        /// @param x Current input (n x 1)
+        /// @param grad Current gradient (n x 1)
         virtual void verify_gradient(Problem &objFunc, const TVector &x, const TVector &grad) final;
 
     private:
@@ -128,8 +139,9 @@ namespace polysolve::nonlinear
         // ====================================================================
         //                           Solver state
         // ====================================================================
-
-        // Reset the solver at the start of a minimization
+
+        /// @brief Reset the solver at the start of a minimization
+        /// @param ndof number of degrees of freedom
         void reset(const int ndof);
 
         std::string descent_strategy_name() const { return m_strategies[m_descent_strategy]->name(); };
@@ -143,8 +155,14 @@ namespace polysolve::nonlinear
         //                            Solver info
         // ====================================================================
 
+        /// @brief Update solver info JSON object
+        /// @param energy 
         void update_solver_info(const double energy);
+
+        /// @brief Reset timing members to 0
         void reset_times();
+
+        /// @brief Log time taken in different phases of the solve
         void log_times() const;
 
         json solver_info;

src/polysolve/nonlinear/descent_strategies/DescentStrategy.hpp

@@ -13,6 +13,10 @@ namespace polysolve::nonlinear
         using TVector = Problem::TVector;
         using Scalar = Problem::Scalar;
 
+        /// @brief Constructor 
+        /// @param solver_params_ JSON of solver parameters
+        /// @param characteristic_length 
+        /// @param logger 
         DescentStrategy(const json &solver_params_,
                         const double characteristic_length,
                         spdlog::logger &logger)
@@ -24,12 +28,22 @@ namespace polysolve::nonlinear
 
         virtual void reset(const int ndof) {}
         virtual void reset_times() {}
+
+        /// @brief Update solver info after finding descent direction
+        /// @param solver_info JSON of solver parameters
+        /// @param per_iteration Number of iterations (used to normalize timings)
         virtual void update_solver_info(json &solver_info, const double per_iteration) {}
         virtual void log_times() const {}
 
         virtual bool is_direction_descent() { return true; }
         virtual bool handle_error() { return false; }
 
+        /// @brief Compute descent direction along which to do line search
+        /// @param objFunc Problem to be minimized
+        /// @param x Current input (n x 1)
+        /// @param grad Current gradient (n x 1)
+        /// @param direction Current descent direction (n x 1)
+        /// @return True if a descent direction was successfully found
         virtual bool compute_update_direction(
             Problem &objFunc,
             const TVector &x,

src/polysolve/nonlinear/line_search/LineSearch.cpp

@@ -119,7 +119,7 @@ namespace polysolve::nonlinear::line_search
         {
             POLYSOLVE_SCOPED_STOPWATCH("CCD narrow-phase", narrow_phase_ccd_time, m_logger);
             m_logger.trace("Performing narrow-phase CCD");
-            step_size = compute_collision_free_step_size(x, delta_x, objFunc, step_size);
+            step_size = compute_max_step_size(x, delta_x, objFunc, step_size);
             if (std::isnan(step_size))
                 return NaN;
         }
@@ -217,7 +217,8 @@ namespace polysolve::nonlinear::line_search
         return step_size;
     }
 
-    double LineSearch::compute_collision_free_step_size(
+    // change to max_step_size
+    double LineSearch::compute_max_step_size(
         const TVector &x,
         const TVector &delta_x,
         Problem &objFunc,

src/polysolve/nonlinear/line_search/LineSearch.hpp

@@ -14,20 +14,34 @@ namespace polysolve::nonlinear::line_search
     public:
         using Scalar = typename Problem::Scalar;
         using TVector = typename Problem::TVector;
-
+      
     public:
+        /// @brief Constructor for creating new LineSearch Object
+        /// @param params JSON of solver parameters
+        /// @param m_logger 
         LineSearch(const json &params, spdlog::logger &m_logger);
         virtual ~LineSearch() = default;
 
+        /// @brief 
+        /// @param x Current input vector (n x 1)
+        /// @param x_delta Current descent direction (n x 1)
+        /// @param objFunc Objective function
+        /// @return 
         double line_search(
             const TVector &x,
-            const TVector &grad,
+            const TVector &x_delta,
             Problem &objFunc);
 
+        /// @brief Dispatch function for creating appropriate subclass
+        /// @param params JSON of solver parameters
+        /// @param logger 
+        /// @return Pointer to object of the specified subclass
         static std::shared_ptr<LineSearch> create(
             const json &params,
             spdlog::logger &logger);
 
+        /// @brief  Get list of available method names
+        /// @return Vector of names of available methods
         static std::vector<std::string> available_methods();
 
         virtual std::string name() const = 0;
@@ -62,6 +76,15 @@ namespace polysolve::nonlinear::line_search
         double use_grad_norm_tol = -1;
 
     protected:
+        /// @brief Compute step size to use during line search 
+        /// @param x Current input (n x 1)
+        /// @param delta_x Current step direction (n x 1)
+        /// @param objFunc Problem to be minimized
+        /// @param use_grad_norm Whether to compare grad norm or energy norm in stopping criteria
+        /// @param old_energy Previous energy (scalar)
+        /// @param old_grad Previous gradient (n x 1)
+        /// @param starting_step_size Initial step size
+        /// @return Step size to use in line search
         virtual double compute_descent_step_size(
             const TVector &x,
             const TVector &delta_x,
@@ -76,13 +99,26 @@ namespace polysolve::nonlinear::line_search
         int cur_iter;
 
     private:
+        /// @brief Compute step size that avoids nan/infinite energy
+        /// @param x Current input (n x 1)
+        /// @param delta_x Current step direction (n x 1)
+        /// @param objFunc Problem to be minimized
+        /// @param starting_step_size Initial step size
+        /// @param rate Rate at which to decrease step size if too large
+        /// @return Nan free step size
         double compute_nan_free_step_size(
             const TVector &x,
             const TVector &delta_x,
             Problem &objFunc,
             const double starting_step_size, const double rate);
 
-        double compute_collision_free_step_size(
+        /// @brief Compute maximum valid step size
+        /// @param x Current input (n x 1)
+        /// @param delta_x Current step direction (n x 1)
+        /// @param objFunc Problem to be minimized
+        /// @param starting_step_size Initial step size
+        /// @return Maximum valid step size (NaN if it is 0)
+        double compute_max_step_size(
             const TVector &x,
             const TVector &delta_x,
             Problem &objFunc,