Merge pull request #2820 from stan-dev/feature/issue-2828-algebra_solver_variadic

Feature/issue 2828 algebra solver variadic

Author: WardBrian

make/tests

@@ -82,7 +82,7 @@ $(IDAS_TESTS) : $(LIBSUNDIALS)
 # KINSOL tests
 ##
 
-ALGEBRA_SOLVER_TESTS := $(subst .cpp,$(EXE),$(call findfiles,test,*algebra_solver*_test.cpp))
+ALGEBRA_SOLVER_TESTS := $(subst .cpp,$(EXE),$(call findfiles,test,*solve*_test.cpp))
 $(ALGEBRA_SOLVER_TESTS) : $(LIBSUNDIALS)
 
 ### These can be generated by the jumbo tests and include the above kinds

stan/math/rev/functor.hpp

@@ -2,8 +2,8 @@
 #define STAN_MATH_REV_FUNCTOR_HPP
 
 #include <stan/math/rev/functor/algebra_solver_fp.hpp>
-#include <stan/math/rev/functor/algebra_solver_powell.hpp>
-#include <stan/math/rev/functor/algebra_solver_newton.hpp>
+#include <stan/math/rev/functor/solve_powell.hpp>
+#include <stan/math/rev/functor/solve_newton.hpp>
 #include <stan/math/rev/functor/algebra_system.hpp>
 #include <stan/math/rev/functor/apply_scalar_unary.hpp>
 #include <stan/math/rev/functor/apply_scalar_binary.hpp>

…th/rev/functor/algebra_solver_newton.hpp stan/math/rev/functor/solve_newton.hppstan/math/rev/functor/algebra_solver_newton.hpp renamed to stan/math/rev/functor/solve_newton.hpp stan/math/rev/functor/algebra_solver_newton.hpp renamed to stan/math/rev/functor/solve_newton.hpp

@@ -1,5 +1,5 @@
-#ifndef STAN_MATH_REV_FUNCTOR_ALGEBRA_SOLVER_NEWTON_HPP
-#define STAN_MATH_REV_FUNCTOR_ALGEBRA_SOLVER_NEWTON_HPP
+#ifndef STAN_MATH_REV_FUNCTOR_SOLVE_NEWTON_HPP
+#define STAN_MATH_REV_FUNCTOR_SOLVE_NEWTON_HPP
 
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/functor/algebra_system.hpp>
@@ -56,21 +56,19 @@ namespace math {
 template <typename F, typename T, typename... Args,
           require_eigen_vector_t<T>* = nullptr,
           require_all_st_arithmetic<Args...>* = nullptr>
-Eigen::VectorXd algebra_solver_newton_impl(const F& f, const T& x,
-                                           std::ostream* const msgs,
-                                           const double scaling_step_size,
-                                           const double function_tolerance,
-                                           const int64_t max_num_steps,
-                                           const Args&... args) {
+Eigen::VectorXd solve_newton_tol(const F& f, const T& x,
+                                 const double scaling_step_size,
+                                 const double function_tolerance,
+                                 const int64_t max_num_steps,
+                                 std::ostream* const msgs,
+                                 const Args&... args) {
   const auto& x_ref = to_ref(value_of(x));
 
-  check_nonzero_size("algebra_solver_newton", "initial guess", x_ref);
-  check_finite("algebra_solver_newton", "initial guess", x_ref);
-  check_nonnegative("algebra_solver_newton", "scaling_step_size",
-                    scaling_step_size);
-  check_nonnegative("algebra_solver_newton", "function_tolerance",
-                    function_tolerance);
-  check_positive("algebra_solver_newton", "max_num_steps", max_num_steps);
+  check_nonzero_size("solve_newton", "initial guess", x_ref);
+  check_finite("solve_newton", "initial guess", x_ref);
+  check_nonnegative("solve_newton", "scaling_step_size", scaling_step_size);
+  check_nonnegative("solve_newton", "function_tolerance", function_tolerance);
+  check_positive("solve_newton", "max_num_steps", max_num_steps);
 
   return kinsol_solve(f, x_ref, scaling_step_size, function_tolerance,
                       max_num_steps, 1, 10, KIN_LINESEARCH, msgs, args...);
@@ -137,10 +135,10 @@ Eigen::VectorXd algebra_solver_newton_impl(const F& f, const T& x,
 template <typename F, typename T, typename... T_Args,
           require_eigen_vector_t<T>* = nullptr,
           require_any_st_var<T_Args...>* = nullptr>
-Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
-    const F& f, const T& x, std::ostream* const msgs,
-    const double scaling_step_size, const double function_tolerance,
-    const int64_t max_num_steps, const T_Args&... args) {
+Eigen::Matrix<var, Eigen::Dynamic, 1> solve_newton_tol(
+    const F& f, const T& x, const double scaling_step_size,
+    const double function_tolerance, const int64_t max_num_steps,
+    std::ostream* const msgs, const T_Args&... args) {
   const auto& x_ref = to_ref(value_of(x));
   auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(args)...));
   auto args_vals_tuple = math::apply(
@@ -149,13 +147,11 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
       },
       *arena_args_tuple);
 
-  check_nonzero_size("algebra_solver_newton", "initial guess", x_ref);
-  check_finite("algebra_solver_newton", "initial guess", x_ref);
-  check_nonnegative("algebra_solver_newton", "scaling_step_size",
-                    scaling_step_size);
-  check_nonnegative("algebra_solver_newton", "function_tolerance",
-                    function_tolerance);
-  check_positive("algebra_solver_newton", "max_num_steps", max_num_steps);
+  check_nonzero_size("solve_newton", "initial guess", x_ref);
+  check_finite("solve_newton", "initial guess", x_ref);
+  check_nonnegative("solve_newton", "scaling_step_size", scaling_step_size);
+  check_nonnegative("solve_newton", "function_tolerance", function_tolerance);
+  check_positive("solve_newton", "max_num_steps", max_num_steps);
 
   // Solve the system
   Eigen::VectorXd theta_dbl = math::apply(
@@ -204,6 +200,56 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
   return ret_type(ret);
 }
 
+/**
+ * Return the solution to the specified system of algebraic
+ * equations given an initial guess, and parameters and data,
+ * which get passed into the algebraic system. Use the
+ * KINSOL solver from the SUNDIALS suite.
+ *
+ * This signature does not give users control over the tuning parameters
+ * and instead relies on default values.
+ *
+ * @tparam F type of equation system function
+ * @tparam T type of elements in the x vector
+ * @tparam Args types of additional input to the equation system functor
+ *
+ * @param[in] f Functor that evaluates the system of equations.
+ * @param[in] x Vector of starting values (initial guess).
+ * @param[in, out] msgs The print stream for warning messages.
+ * @param[in] scaling_step_size Scaled-step stopping tolerance. If
+ *            a Newton step is smaller than the scaling step
+ *            tolerance, the code breaks, assuming the solver is no
+ *            longer making significant progress (i.e. is stuck)
+ * @param[in] function_tolerance determines whether roots are acceptable.
+ * @param[in] max_num_steps  maximum number of function evaluations.
+ * @param[in, out] msgs the print stream for warning messages.
+ * @param[in] args Additional parameters to the equation system functor.
+ * @return theta Vector of solutions to the system of equations.
+ * @throw <code>std::invalid_argument</code> if x has size zero.
+ * @throw <code>std::invalid_argument</code> if x has non-finite elements.
+ * @throw <code>std::invalid_argument</code> if scaled_step_size is strictly
+ * negative.
+ * @throw <code>std::invalid_argument</code> if function_tolerance is strictly
+ * negative.
+ * @throw <code>std::invalid_argument</code> if max_num_steps is not positive.
+ * @throw <code>std::domain_error if solver exceeds max_num_steps.
+ */
+template <typename F, typename T, typename... T_Args,
+          require_eigen_vector_t<T>* = nullptr>
+Eigen::Matrix<stan::return_type_t<T_Args...>, Eigen::Dynamic, 1> solve_newton(
+    const F& f, const T& x, std::ostream* const msgs, const T_Args&... args) {
+  double scaling_step_size = 1e-3;
+  double function_tolerance = 1e-6;
+  int64_t max_num_steps = 200;
+  const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
+  return math::apply(
+      [&](const auto&... args_refs) {
+        return solve_newton_tol(f, x, scaling_step_size, function_tolerance,
+                                max_num_steps, msgs, args_refs...);
+      },
+      args_ref_tuple);
+}
+
 /**
  * Return the solution to the specified system of algebraic
  * equations given an initial guess, and parameters and data,
@@ -252,9 +298,9 @@ Eigen::Matrix<scalar_type_t<T2>, Eigen::Dynamic, 1> algebra_solver_newton(
     const double scaling_step_size = 1e-3,
     const double function_tolerance = 1e-6,
     const long int max_num_steps = 200) {  // NOLINT(runtime/int)
-  return algebra_solver_newton_impl(algebra_solver_adapter<F>(f), x, msgs,
-                                    scaling_step_size, function_tolerance,
-                                    max_num_steps, y, dat, dat_int);
+  return solve_newton_tol(algebra_solver_adapter<F>(f), x, scaling_step_size,
+                          function_tolerance, max_num_steps, msgs, y, dat,
+                          dat_int);
 }
 
 }  // namespace math

…th/rev/functor/algebra_solver_powell.hpp stan/math/rev/functor/solve_powell.hppstan/math/rev/functor/algebra_solver_powell.hpp renamed to stan/math/rev/functor/solve_powell.hpp stan/math/rev/functor/algebra_solver_powell.hpp renamed to stan/math/rev/functor/solve_powell.hpp

@@ -1,5 +1,5 @@
-#ifndef STAN_MATH_REV_FUNCTOR_ALGEBRA_SOLVER_POWELL_HPP
-#define STAN_MATH_REV_FUNCTOR_ALGEBRA_SOLVER_POWELL_HPP
+#ifndef STAN_MATH_REV_FUNCTOR_SOLVE_POWELL_HPP
+#define STAN_MATH_REV_FUNCTOR_SOLVE_POWELL_HPP
 
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
@@ -45,11 +45,10 @@ namespace math {
  */
 template <typename F, typename T, typename... Args,
           require_eigen_vector_t<T>* = nullptr>
-T& algebra_solver_powell_call_solver(const F& f, T& x, std::ostream* const msgs,
-                                     const double relative_tolerance,
-                                     const double function_tolerance,
-                                     const int64_t max_num_steps,
-                                     const Args&... args) {
+T& solve_powell_call_solver(const F& f, T& x, std::ostream* const msgs,
+                            const double relative_tolerance,
+                            const double function_tolerance,
+                            const int64_t max_num_steps, const Args&... args) {
   // Construct the solver
   hybrj_functor_solver<F> hfs(f);
   Eigen::HybridNonLinearSolver<hybrj_functor_solver<F>> solver(hfs);
@@ -103,11 +102,11 @@ T& algebra_solver_powell_call_solver(const F& f, T& x, std::ostream* const msgs,
  *
  * @param[in] f Functor that evaluates the system of equations.
  * @param[in] x Vector of starting values (initial guess).
- * @param[in, out] msgs the print stream for warning messages.
  * @param[in] relative_tolerance determines the convergence criteria
  *            for the solution.
  * @param[in] function_tolerance determines whether roots are acceptable.
  * @param[in] max_num_steps maximum number of function evaluations.
+ * @param[in, out] msgs the print stream for warning messages.
  * @param[in] args additional parameters to the equation system functor.
  * @return theta Vector of solutions to the system of equations.
  * @pre f returns finite values when passed any value of x and the given args.
@@ -126,12 +125,12 @@ T& algebra_solver_powell_call_solver(const F& f, T& x, std::ostream* const msgs,
 template <typename F, typename T, typename... Args,
           require_eigen_vector_t<T>* = nullptr,
           require_all_st_arithmetic<Args...>* = nullptr>
-Eigen::VectorXd algebra_solver_powell_impl(const F& f, const T& x,
-                                           std::ostream* const msgs,
-                                           const double relative_tolerance,
-                                           const double function_tolerance,
-                                           const int64_t max_num_steps,
-                                           const Args&... args) {
+Eigen::VectorXd solve_powell_tol(const F& f, const T& x,
+                                 const double relative_tolerance,
+                                 const double function_tolerance,
+                                 const int64_t max_num_steps,
+                                 std::ostream* const msgs,
+                                 const Args&... args) {
   auto x_ref = eval(value_of(x));
   auto args_vals_tuple = std::make_tuple(to_ref(args)...);
 
@@ -141,20 +140,18 @@ Eigen::VectorXd algebra_solver_powell_impl(const F& f, const T& x,
         args_vals_tuple);
   };
 
-  check_nonzero_size("algebra_solver_powell", "initial guess", x_ref);
-  check_finite("algebra_solver_powell", "initial guess", x_ref);
+  check_nonzero_size("solve_powell", "initial guess", x_ref);
+  check_finite("solve_powell", "initial guess", x_ref);
   check_nonnegative("alegbra_solver_powell", "relative_tolerance",
                     relative_tolerance);
-  check_nonnegative("algebra_solver_powell", "function_tolerance",
-                    function_tolerance);
-  check_positive("algebra_solver_powell", "max_num_steps", max_num_steps);
-  check_matching_sizes("algebra_solver", "the algebraic system's output",
+  check_nonnegative("solve_powell", "function_tolerance", function_tolerance);
+  check_positive("solve_powell", "max_num_steps", max_num_steps);
+  check_matching_sizes("solve_powell", "the algebraic system's output",
                        f_wrt_x(x_ref), "the vector of unknowns, x,", x_ref);
 
   // Solve the system
-  return algebra_solver_powell_call_solver(f_wrt_x, x_ref, msgs,
-                                           relative_tolerance,
-                                           function_tolerance, max_num_steps);
+  return solve_powell_call_solver(f_wrt_x, x_ref, msgs, relative_tolerance,
+                                  function_tolerance, max_num_steps);
 }
 
 /**
@@ -163,31 +160,20 @@ Eigen::VectorXd algebra_solver_powell_impl(const F& f, const T& x,
  * which get passed into the algebraic system.
  * Use Powell's dogleg solver.
  *
- * The user can also specify the relative tolerance
- * (xtol in Eigen's code), the function tolerance,
- * and the maximum number of steps (maxfev in Eigen's code).
+ * This signature does not let the user specify the tuning parameters of the
+ * solver (instead default values are used).
  *
  * @tparam F type of equation system function
- * @tparam T1 type of elements in the x vector
- * @tparam T2 type of elements in the y vector
+ * @tparam T type of elements in the x vector
+ * @tparam Args types of additional input to the equation system functor
  *
  * @param[in] f Functor that evaluates the system of equations.
  * @param[in] x Vector of starting values (initial guess).
- * @param[in] y parameter vector for the equation system.
- * @param[in] dat continuous data vector for the equation system.
- * @param[in] dat_int integer data vector for the equation system.
  * @param[in, out] msgs the print stream for warning messages.
- * @param[in] relative_tolerance determines the convergence criteria
- *            for the solution.
- * @param[in] function_tolerance determines whether roots are acceptable.
- * @param[in] max_num_steps  maximum number of function evaluations.
+ * @param[in] args Additional parameters to the equation system functor.
  * @return theta Vector of solutions to the system of equations.
  * @throw <code>std::invalid_argument</code> if x has size zero.
  * @throw <code>std::invalid_argument</code> if x has non-finite elements.
- * @throw <code>std::invalid_argument</code> if y has non-finite elements.
- * @throw <code>std::invalid_argument</code> if dat has non-finite elements.
- * @throw <code>std::invalid_argument</code> if dat_int has non-finite
- * elements.
  * @throw <code>std::invalid_argument</code> if relative_tolerance is strictly
  * negative.
  * @throw <code>std::invalid_argument</code> if function_tolerance is strictly
@@ -197,17 +183,20 @@ Eigen::VectorXd algebra_solver_powell_impl(const F& f, const T& x,
  * @throw <code>std::domain_error</code> if the norm of the solution exceeds
  * the function tolerance.
  */
-template <typename F, typename T1, typename T2,
-          require_all_eigen_vector_t<T1, T2>* = nullptr>
-Eigen::Matrix<value_type_t<T2>, Eigen::Dynamic, 1> algebra_solver_powell(
-    const F& f, const T1& x, const T2& y, const std::vector<double>& dat,
-    const std::vector<int>& dat_int, std::ostream* const msgs = nullptr,
-    const double relative_tolerance = 1e-10,
-    const double function_tolerance = 1e-6,
-    const int64_t max_num_steps = 1e+3) {
-  return algebra_solver_powell_impl(algebra_solver_adapter<F>(f), x, msgs,
-                                    relative_tolerance, function_tolerance,
-                                    max_num_steps, y, dat, dat_int);
+template <typename F, typename T, typename... T_Args,
+          require_eigen_vector_t<T>* = nullptr>
+Eigen::Matrix<stan::return_type_t<T_Args...>, Eigen::Dynamic, 1> solve_powell(
+    const F& f, const T& x, std::ostream* const msgs, const T_Args&... args) {
+  double relative_tolerance = 1e-10;
+  double function_tolerance = 1e-6;
+  int64_t max_num_steps = 200;
+  const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
+  return math::apply(
+      [&](const auto&... args_refs) {
+        return solve_powell_tol(f, x, relative_tolerance, function_tolerance,
+                                max_num_steps, msgs, args_refs...);
+      },
+      args_ref_tuple);
 }
 
 /**
@@ -265,8 +254,9 @@ Eigen::Matrix<value_type_t<T2>, Eigen::Dynamic, 1> algebra_solver(
     const double relative_tolerance = 1e-10,
     const double function_tolerance = 1e-6,
     const int64_t max_num_steps = 1e+3) {
-  return algebra_solver_powell(f, x, y, dat, dat_int, msgs, relative_tolerance,
-                               function_tolerance, max_num_steps);
+  return solve_powell_tol(algebra_solver_adapter<F>(f), x, relative_tolerance,
+                          function_tolerance, max_num_steps, msgs, y, dat,
+                          dat_int);
 }
 
 /**
@@ -309,11 +299,11 @@ Eigen::Matrix<value_type_t<T2>, Eigen::Dynamic, 1> algebra_solver(
  *
  * @param[in] f Functor that evaluates the system of equations.
  * @param[in] x Vector of starting values (initial guess).
- * @param[in, out] msgs the print stream for warning messages.
  * @param[in] relative_tolerance determines the convergence criteria
  *            for the solution.
  * @param[in] function_tolerance determines whether roots are acceptable.
  * @param[in] max_num_steps  maximum number of function evaluations.
+ * @param[in, out] msgs the print stream for warning messages.
  * @param[in] args Additional parameters to the equation system functor.
  * @return theta Vector of solutions to the system of equations.
  * @pre f returns finite values when passed any value of x and the given args.
@@ -332,10 +322,10 @@ Eigen::Matrix<value_type_t<T2>, Eigen::Dynamic, 1> algebra_solver(
 template <typename F, typename T, typename... T_Args,
           require_eigen_vector_t<T>* = nullptr,
           require_any_st_var<T_Args...>* = nullptr>
-Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_powell_impl(
-    const F& f, const T& x, std::ostream* const msgs,
-    const double relative_tolerance, const double function_tolerance,
-    const int64_t max_num_steps, const T_Args&... args) {
+Eigen::Matrix<var, Eigen::Dynamic, 1> solve_powell_tol(
+    const F& f, const T& x, const double relative_tolerance,
+    const double function_tolerance, const int64_t max_num_steps,
+    std::ostream* const msgs, const T_Args&... args) {
   auto x_ref = eval(value_of(x));
   auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(args)...));
   auto args_vals_tuple = math::apply(
@@ -350,19 +340,18 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_powell_impl(
         args_vals_tuple);
   };
 
-  check_nonzero_size("algebra_solver_powell", "initial guess", x_ref);
-  check_finite("algebra_solver_powell", "initial guess", x_ref);
+  check_nonzero_size("solve_powell", "initial guess", x_ref);
+  check_finite("solve_powell", "initial guess", x_ref);
   check_nonnegative("alegbra_solver_powell", "relative_tolerance",
                     relative_tolerance);
-  check_nonnegative("algebra_solver_powell", "function_tolerance",
-                    function_tolerance);
-  check_positive("algebra_solver_powell", "max_num_steps", max_num_steps);
-  check_matching_sizes("algebra_solver", "the algebraic system's output",
+  check_nonnegative("solve_powell", "function_tolerance", function_tolerance);
+  check_positive("solve_powell", "max_num_steps", max_num_steps);
+  check_matching_sizes("solve_powell", "the algebraic system's output",
                        f_wrt_x(x_ref), "the vector of unknowns, x,", x_ref);
 
   // Solve the system
-  algebra_solver_powell_call_solver(f_wrt_x, x_ref, msgs, relative_tolerance,
-                                    function_tolerance, max_num_steps);
+  solve_powell_call_solver(f_wrt_x, x_ref, msgs, relative_tolerance,
+                           function_tolerance, max_num_steps);
 
   Eigen::MatrixXd Jf_x;
   Eigen::VectorXd f_x;