diff --git a/stan/math/prim/functor/algebra_solver_adapter.hpp b/stan/math/prim/functor/algebra_solver_adapter.hpp new file mode 100644 index 00000000000..528912c660e --- /dev/null +++ b/stan/math/prim/functor/algebra_solver_adapter.hpp @@ -0,0 +1,26 @@ +#ifndef STAN_MATH_PRIM_FUNCTOR_ALGEBRA_SOLVER_ADAPTER_HPP +#define STAN_MATH_PRIM_FUNCTOR_ALGEBRA_SOLVER_ADAPTER_HPP + +#include +#include + +/** + * Adapt the non-variadic algebra_solver_newton and algebra_solver_powell + * arguemts to the variadic algebra_solver_newton_impl and + * algebra_solver_powell_impl interfaces. + * + * @tparam F type of function to adapt. + */ +template +struct algebra_solver_adapter { + const F& f_; + + explicit algebra_solver_adapter(const F& f) : f_(f) {} + + template + auto operator()(const T1& x, std::ostream* msgs, T2&&... args) const { + return f_(x, args..., msgs); + } +}; + +#endif diff --git a/stan/math/rev/core/chainable_object.hpp b/stan/math/rev/core/chainable_object.hpp index 50232476354..0631d88354c 100644 --- a/stan/math/rev/core/chainable_object.hpp +++ b/stan/math/rev/core/chainable_object.hpp @@ -62,6 +62,61 @@ auto make_chainable_ptr(T&& obj) { return &ptr->get(); } +/** + * `unsafe_chainable_object` hold another object and is useful for connecting + * the lifetime of a specific object to the chainable stack. This class + * differs from `chainable_object` in that this class does not evaluate + * expressions. + * + * `unsafe_chainable_object` objects should only be allocated with `new`. + * `unsafe_chainable_object` objects allocated on the stack will result + * in a double free (`obj_` will get destructed once when the + * `unsafe_chainable_object` leaves scope and once when the chainable + * stack memory is recovered). + * + * @tparam T type of object to hold + */ +template +class unsafe_chainable_object : public chainable_alloc { + private: + std::decay_t obj_; + + public: + /** + * Construct chainable object from another object + * + * @tparam S type of object to hold (must have the same plain type as `T`) + */ + template , plain_type_t>* = nullptr> + explicit unsafe_chainable_object(S&& obj) : obj_(std::forward(obj)) {} + + /** + * Return a reference to the underlying object + * + * @return reference to underlying object + */ + inline auto& get() noexcept { return obj_; } + inline const auto& get() const noexcept { return obj_; } +}; + +/** + * Store the given object in a `chainable_object` so it is destructed + * only when the chainable stack memory is recovered and return + * a pointer to the underlying object This function + * differs from `make_chainable_object` in that this class does not evaluate + * expressions. + * + * @tparam T type of object to hold + * @param obj object to hold + * @return pointer to object held in `chainable_object` + */ +template +auto make_unsafe_chainable_ptr(T&& obj) { + auto ptr = new unsafe_chainable_object(std::forward(obj)); + return &ptr->get(); +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/functor/algebra_solver_newton.hpp b/stan/math/rev/functor/algebra_solver_newton.hpp index f3a04386789..51809bed762 100644 --- a/stan/math/rev/functor/algebra_solver_newton.hpp +++ b/stan/math/rev/functor/algebra_solver_newton.hpp @@ -3,11 +3,10 @@ #include #include -#include #include #include -#include #include +#include #include #include #include @@ -25,16 +24,14 @@ namespace math { * The user can also specify the scaled step size, the function * tolerance, and the maximum number of steps. * + * This overload handles non-autodiff parameters. + * * @tparam F type of equation system function. * @tparam T type of initial guess vector. + * @tparam Args types of additional parameters to the equation system functor * * @param[in] f Functor that evaluated the system of equations. * @param[in] x Vector of starting values. - * @param[in] y Parameter vector for the equation system. The function - * is overloaded to treat y as a vector of doubles or of a - * a template type T. - * @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] scaling_step_size Scaled-step stopping tolerance. If * a Newton step is smaller than the scaling step @@ -42,36 +39,165 @@ namespace math { * 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. - * * @throw std::invalid_argument if x has size zero. + * @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. + * @throw std::invalid_argument if x has size zero. + * @throw std::invalid_argument if x has non-finite elements. + * @throw std::invalid_argument if scaled_step_size is strictly + * negative. + * @throw std::invalid_argument if function_tolerance is strictly + * negative. + * @throw std::invalid_argument if max_num_steps is not positive. + * @throw std::domain_error if solver exceeds max_num_steps. + */ +template * = nullptr, + require_all_st_arithmetic* = 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) { + 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); + + return kinsol_solve(f, x_ref, scaling_step_size, function_tolerance, + max_num_steps, 1, 10, KIN_LINESEARCH, msgs, args...); +} + +/** + * 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. + * + * The user can also specify the scaled step size, the function + * tolerance, and the maximum number of steps. + * + * This overload handles var parameters. + * + * The Jacobian \(J_{xy}\) (i.e., Jacobian of unknown \(x\) w.r.t. the parameter + * \(y\)) is calculated given the solution as follows. Since + * \[ + * f(x, y) = 0, + * \] + * we have (\(J_{pq}\) being the Jacobian matrix \(\tfrac {dq} {dq}\)) + * \[ + * - J_{fx} J_{xy} = J_{fy}, + * \] + * and therefore \(J_{xy}\) can be solved from system + * \[ + * - J_{fx} J_{xy} = J_{fy}. + * \] + * Let \(eta\) be the adjoint with respect to \(x\); then to calculate + * \[ + * \eta J_{xy}, + * \] + * we solve + * \[ + * - (\eta J_{fx}^{-1}) J_{fy}. + * \] + * + * @tparam F type of equation system function. + * @tparam T type of initial guess vector. + * @tparam Args types of additional parameters to the equation system functor + * + * @param[in] f Functor that evaluated the system of equations. + * @param[in] x Vector of starting values. + * @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] 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. + * @throw std::invalid_argument if x has size zero. * @throw std::invalid_argument if x has non-finite elements. - * @throw std::invalid_argument if y has non-finite elements. - * @throw std::invalid_argument if dat has non-finite elements. - * @throw std::invalid_argument if dat_int has non-finite elements. * @throw std::invalid_argument if scaled_step_size is strictly * negative. * @throw std::invalid_argument if function_tolerance is strictly * negative. * @throw std::invalid_argument if max_num_steps is not positive. - * @throw std::domain_error if solver exceeds max_num_steps. + * @throw std::domain_error if solver exceeds max_num_steps. */ -template * = nullptr> -Eigen::VectorXd algebra_solver_newton( - const F& f, const T& x, const Eigen::VectorXd& y, - const std::vector& dat, const std::vector& dat_int, - std::ostream* msgs = nullptr, double scaling_step_size = 1e-3, - double function_tolerance = 1e-6, - long int max_num_steps = 200) { // NOLINT(runtime/int) - const auto& x_eval = x.eval(); - algebra_solver_check(x_eval, y, dat, dat_int, function_tolerance, - max_num_steps); - check_nonnegative("algebra_solver", "scaling_step_size", scaling_step_size); - - check_matching_sizes("algebra_solver", "the algebraic system's output", - value_of(f(x_eval, y, dat, dat_int, msgs)), - "the vector of unknowns, x,", x); - - return kinsol_solve(f, value_of(x_eval), y, dat, dat_int, 0, - scaling_step_size, function_tolerance, max_num_steps); +template * = nullptr, + require_any_st_var* = nullptr> +Eigen::Matrix 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) { + const auto& x_ref = to_ref(value_of(x)); + auto arena_args_tuple = std::make_tuple(to_arena(args)...); + auto args_vals_tuple = apply( + [&](const auto&... args) { + return std::make_tuple(to_ref(value_of(args))...); + }, + 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); + + // Solve the system + Eigen::VectorXd theta_dbl = apply( + [&](const auto&... vals) { + return kinsol_solve(f, x_ref, scaling_step_size, function_tolerance, + max_num_steps, 1, 10, KIN_LINESEARCH, msgs, + vals...); + }, + args_vals_tuple); + + auto f_wrt_x = [&](const auto& x) { + return apply([&](const auto&... args) { return f(x, msgs, args...); }, + args_vals_tuple); + }; + + Eigen::MatrixXd Jf_x; + Eigen::VectorXd f_x; + + jacobian(f_wrt_x, theta_dbl, f_x, Jf_x); + + using ret_type = Eigen::Matrix; + arena_t ret = theta_dbl; + auto Jf_x_T_lu_ptr + = make_unsafe_chainable_ptr(Jf_x.transpose().partialPivLu()); // Lu + + reverse_pass_callback([f, ret, arena_args_tuple, Jf_x_T_lu_ptr, + msgs]() mutable { + Eigen::VectorXd eta = -Jf_x_T_lu_ptr->solve(ret.adj().eval()); + + // Contract with Jacobian of f with respect to y using a nested reverse + // autodiff pass. + { + nested_rev_autodiff rev; + + Eigen::VectorXd ret_val = ret.val(); + auto x_nrad_ = apply( + [&](const auto&... args) { return eval(f(ret_val, msgs, args...)); }, + arena_args_tuple); + x_nrad_.adj() = eta; + grad(); + } + }); + + return ret_type(ret); } /** @@ -83,19 +209,15 @@ Eigen::VectorXd algebra_solver_newton( * The user can also specify the scaled step size, the function * tolerance, and the maximum number of steps. * - * Overload the previous definition to handle the case where y - * is a vector of parameters (var). The overload calls the - * algebraic solver defined above and builds a vari object on - * top, using the algebra_solver_vari class. + * Signature to maintain backward compatibility, will be removed + * in the future. * * @tparam F type of equation system function. * @tparam T type of initial guess vector. * * @param[in] f Functor that evaluated the system of equations. * @param[in] x Vector of starting values. - * @param[in] y Parameter vector for the equation system. The function - * is overloaded to treat y as a vector of doubles or of a - * a template type T. + * @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. @@ -119,34 +241,16 @@ Eigen::VectorXd algebra_solver_newton( * @throw std::domain_error if solver exceeds max_num_steps. */ template * = nullptr, - require_st_var* = nullptr> + require_all_eigen_vector_t* = nullptr> Eigen::Matrix, Eigen::Dynamic, 1> algebra_solver_newton( const F& f, const T1& x, const T2& y, const std::vector& dat, - const std::vector& dat_int, std::ostream* msgs = nullptr, - double scaling_step_size = 1e-3, double function_tolerance = 1e-6, - long int max_num_steps = 200) { // NOLINT(runtime/int) - - const auto& x_eval = x.eval(); - const auto& y_eval = y.eval(); - Eigen::VectorXd theta_dbl = algebra_solver_newton( - f, x_eval, value_of(y_eval), dat, dat_int, msgs, scaling_step_size, - function_tolerance, max_num_steps); - - typedef system_functor Fy; - typedef system_functor Fs; - typedef hybrj_functor_solver Fx; - Fx fx(Fs(), f, value_of(x_eval), value_of(y_eval), dat, dat_int, msgs); - - // Construct vari - auto* vi0 = new algebra_solver_vari, Fx>( - Fy(), f, value_of(x_eval), y_eval, dat, dat_int, theta_dbl, fx, msgs); - Eigen::Matrix theta(x.size()); - theta(0) = var(vi0->theta_[0]); - for (int i = 1; i < x.size(); ++i) - theta(i) = var(vi0->theta_[i]); - - return theta; + const std::vector& dat_int, std::ostream* const msgs = nullptr, + 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), x, msgs, + scaling_step_size, function_tolerance, + max_num_steps, y, dat, dat_int); } } // namespace math diff --git a/stan/math/rev/functor/algebra_solver_powell.hpp b/stan/math/rev/functor/algebra_solver_powell.hpp index a36705c70ac..5415f361cad 100644 --- a/stan/math/rev/functor/algebra_solver_powell.hpp +++ b/stan/math/rev/functor/algebra_solver_powell.hpp @@ -5,8 +5,10 @@ #include #include #include -#include #include +#include +#include +#include #include #include #include @@ -16,66 +18,72 @@ namespace stan { namespace math { /** - * The vari class for the algebraic solver. We compute the Jacobian of - * the solutions with respect to the parameters using the implicit - * function theorem. The call to Jacobian() occurs outside the call to - * chain() -- this prevents malloc issues. + * Solve algebraic equations using Powell solver + * + * @tparam F type of equation system function + * @tparam T type of elements in the x vector + * @tparam Args types of additional parameters 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] 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_dbl Double vector of solutions to the system of equations. + * @pre x has size greater than zero. + * @pre x has only finite elements. + * @pre f returns finite values when passed any value of x and the given args. + * @pre relative_tolerance is non-negative. + * @pre function_tolerance is non-negative. + * @pre max_num_steps is positive. + * @throw std::domain_error solver exceeds max_num_steps. + * @throw std::domain_error if the norm of the solution exceeds + * the function tolerance. */ -template -struct algebra_solver_vari : public vari { - /** vector of parameters */ - vari** y_; - /** number of parameters */ - int y_size_; - /** number of unknowns */ - int x_size_; - /** vector of solution */ - vari** theta_; - /** Jacobian of the solution w.r.t parameters */ - double* Jx_y_; - - algebra_solver_vari(const Fs& fs, const F& f, const Eigen::VectorXd& x, - const Eigen::Matrix& y, - const std::vector& dat, - const std::vector& dat_int, - const Eigen::VectorXd& theta_dbl, Fx& fx, - std::ostream* msgs) - : vari(theta_dbl(0)), - y_(ChainableStack::instance_->memalloc_.alloc_array(y.size())), - y_size_(y.size()), - x_size_(x.size()), - theta_( - ChainableStack::instance_->memalloc_.alloc_array(x_size_)), - Jx_y_(ChainableStack::instance_->memalloc_.alloc_array( - x_size_ * y_size_)) { - using Eigen::Map; - using Eigen::MatrixXd; - for (int i = 0; i < y.size(); ++i) { - y_[i] = y(i).vi_; - } +template * = 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) { + // Construct the solver + hybrj_functor_solver hfs(f); + Eigen::HybridNonLinearSolver> solver(hfs); - theta_[0] = this; - for (int i = 1; i < x.size(); ++i) { - theta_[i] = new vari(theta_dbl(i), false); - } + // Compute theta_dbl + solver.parameters.xtol = relative_tolerance; + solver.parameters.maxfev = max_num_steps; + solver.solve(x); - // Compute the Jacobian and store in array, using the - // implicit function theorem, i.e. Jx_y = Jf_y / Jf_x - using f_y = hybrj_functor_solver; - Map(&Jx_y_[0], x_size_, y_size_) - = -mdivide_left(fx.get_jacobian(theta_dbl), - f_y(fs, f, theta_dbl, value_of(y), dat, dat_int, msgs) - .get_jacobian(value_of(y))); + // Check if the max number of steps has been exceeded + if (solver.nfev >= max_num_steps) { + [max_num_steps]() STAN_COLD_PATH { + throw_domain_error("algebra_solver", "maximum number of iterations", + max_num_steps, "(", ") was exceeded in the solve."); + }(); } - void chain() { - for (int j = 0; j < y_size_; j++) { - for (int i = 0; i < x_size_; i++) { - y_[j]->adj_ += theta_[i]->adj_ * Jx_y_[j * x_size_ + i]; - } - } + // Check solution is a root + double system_norm = f(x).stableNorm(); + if (system_norm > function_tolerance) { + [function_tolerance, system_norm]() STAN_COLD_PATH { + std::ostringstream message; + message << "the norm of the algebraic function is " << system_norm + << " but should be lower than the function " + << "tolerance:"; + throw_domain_error("algebra_solver", message.str().c_str(), + function_tolerance, "", + ". Consider decreasing the relative tolerance and " + "increasing max_num_steps."); + }(); } -}; + + return x; +} /** * Return the solution to the specified system of algebraic @@ -87,32 +95,25 @@ struct algebra_solver_vari : public vari { * (xtol in Eigen's code), the function tolerance, * and the maximum number of steps (maxfev in Eigen's code). * - * Throw an exception if the norm of f(x), where f is the - * output of the algebraic system and x the proposed solution, - * is greater than the function tolerance. We here use the - * norm as a metric to measure how far we are from the origin (0). + * This overload handles non-autodiff parameters. * - * @tparam F type of equation system function. - * @tparam T type of initial guess vector. + * @tparam F type of equation system function + * @tparam T type of elements in the x vector + * @tparam Args types of additional parameters to the equation system functor * * @param[in] f Functor that evaluates the system of equations. - * @param[in] x Vector of starting values. - * @param[in] y parameter vector for the equation system. The function - * is overloaded to treat y as a vector of doubles or of a - * a template type T. - * @param[in] dat continuous data vector for the equation system. - * @param[in] dat_int integer data vector for the equation system. + * @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] 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. + * @pre f returns finite values when passed any value of x and the given args. * @throw std::invalid_argument if x has size zero. * @throw std::invalid_argument if x has non-finite elements. - * @throw std::invalid_argument if y has non-finite elements. - * @throw std::invalid_argument if dat has non-finite elements. - * @throw std::invalid_argument if dat_int has non-finite elements. + * elements. * @throw std::invalid_argument if relative_tolerance is strictly * negative. * @throw std::invalid_argument if function_tolerance is strictly @@ -122,59 +123,38 @@ struct algebra_solver_vari : public vari { * @throw std::domain_error if the norm of the solution exceeds * the function tolerance. */ -template * = nullptr> -Eigen::VectorXd algebra_solver_powell( - const F& f, const T& x, const Eigen::VectorXd& y, - const std::vector& dat, const std::vector& dat_int, - std::ostream* msgs = nullptr, double relative_tolerance = 1e-10, - double function_tolerance = 1e-6, - long int max_num_steps = 1e+3) { // NOLINT(runtime/int) - const auto& x_eval = x.eval(); - const auto& x_val = (value_of(x_eval)).eval(); - algebra_solver_check(x_val, y, dat, dat_int, function_tolerance, - max_num_steps); - check_nonnegative("alegbra_solver", "relative_tolerance", relative_tolerance); - // if (relative_tolerance < 0) - // invalid_argument("algebra_solver", "relative_tolerance,", - // function_tolerance, "", - // ", must be greater than or equal to 0"); +template * = nullptr, + require_all_st_arithmetic* = 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) { + auto x_ref = eval(value_of(x)); + auto args_vals_tuple = std::make_tuple(to_ref(args)...); - // Create functor for algebraic system - using Fs = system_functor; - using Fx = hybrj_functor_solver; - Fx fx(Fs(), f, x_val, y, dat, dat_int, msgs); - Eigen::HybridNonLinearSolver solver(fx); + auto f_wrt_x = [&args_vals_tuple, &f, msgs](const auto& x) { + return apply( + [&x, &f, msgs](const auto&... args) { return f(x, msgs, args...); }, + args_vals_tuple); + }; - // Check dimension unknowns equals dimension of system output + check_nonzero_size("algebra_solver_powell", "initial guess", x_ref); + check_finite("algebra_solver_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", - fx.get_value(x_val), "the vector of unknowns, x,", x); + f_wrt_x(x_ref), "the vector of unknowns, x,", x_ref); - // Compute theta_dbl - Eigen::VectorXd theta_dbl = x_val; - solver.parameters.xtol = relative_tolerance; - solver.parameters.maxfev = max_num_steps; - solver.solve(theta_dbl); - - // Check if the max number of steps has been exceeded - if (solver.nfev >= max_num_steps) { - throw_domain_error("algebra_solver", "maximum number of iterations", - max_num_steps, "(", ") was exceeded in the solve."); - } - - // Check solution is a root - double system_norm = fx.get_value(theta_dbl).stableNorm(); - if (system_norm > function_tolerance) { - std::ostringstream message; - message << "the norm of the algebraic function is " << system_norm - << " but should be lower than the function " - << "tolerance:"; - throw_domain_error("algebra_solver", message.str().c_str(), - function_tolerance, "", - ". Consider decreasing the relative tolerance and " - "increasing max_num_steps."); - } - - return theta_dbl; + // Solve the system + return algebra_solver_powell_call_solver(f_wrt_x, x_ref, msgs, + relative_tolerance, + function_tolerance, max_num_steps); } /** @@ -187,11 +167,6 @@ Eigen::VectorXd algebra_solver_powell( * (xtol in Eigen's code), the function tolerance, * and the maximum number of steps (maxfev in Eigen's code). * - * Overload the previous definition to handle the case where y - * is a vector of parameters (var). The overload calls the - * algebraic solver defined above and builds a vari object on - * top, using the algebra_solver_vari class. - * * @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 @@ -223,41 +198,16 @@ Eigen::VectorXd algebra_solver_powell( * the function tolerance. */ template * = nullptr, - require_st_var* = nullptr> + require_all_eigen_vector_t* = nullptr> Eigen::Matrix, Eigen::Dynamic, 1> algebra_solver_powell( const F& f, const T1& x, const T2& y, const std::vector& dat, - const std::vector& dat_int, std::ostream* msgs = nullptr, - double relative_tolerance = 1e-10, double function_tolerance = 1e-6, - long int max_num_steps = 1e+3) { // NOLINT(runtime/int) - const auto& x_eval = x.eval(); - const auto& y_eval = y.eval(); - const auto& x_val = (value_of(x_eval)).eval(); - const auto& y_val = (value_of(y_eval)).eval(); - Eigen::VectorXd theta_dbl = algebra_solver_powell( - f, x_eval, y_val, dat, dat_int, 0, relative_tolerance, function_tolerance, - max_num_steps); - - using Fy = system_functor; - - // TODO(charlesm93): a similar object gets constructed inside - // the call to algebra_solver. Cache the previous result - // and use it here (if possible). - using Fs = system_functor; - using Fx = hybrj_functor_solver; - Fx fx(Fs(), f, x_val, y_val, dat, dat_int, msgs); - - // Construct vari - algebra_solver_vari, Fx>* vi0 - = new algebra_solver_vari, Fx>( - Fy(), f, x_val, y_eval, dat, dat_int, theta_dbl, fx, msgs); - Eigen::Matrix theta(x.size()); - theta(0) = var(vi0->theta_[0]); - for (int i = 1; i < x.size(); ++i) { - theta(i) = var(vi0->theta_[i]); - } - - return theta; + const std::vector& 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), x, msgs, + relative_tolerance, function_tolerance, + max_num_steps, y, dat, dat_int); } /** @@ -288,6 +238,8 @@ Eigen::Matrix, Eigen::Dynamic, 1> algebra_solver_powell( * @param[in] function_tolerance determines whether roots are acceptable. * @param[in] max_num_steps maximum number of function evaluations. * @return theta Vector of solutions to the system of equations. + * @pre f returns finite values when passed any value of x and the given y, dat, + * and dat_int. * @throw std::invalid_argument if x has size zero. * @throw std::invalid_argument if x has non-finite elements. * @throw std::invalid_argument if y has non-finite elements. @@ -310,12 +262,139 @@ template , Eigen::Dynamic, 1> algebra_solver( const F& f, const T1& x, const T2& y, const std::vector& dat, const std::vector& dat_int, std::ostream* msgs = nullptr, - double relative_tolerance = 1e-10, double function_tolerance = 1e-6, - long int max_num_steps = 1e+3) { // NOLINT(runtime/int) + 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 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 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 overload handles var parameters. + * + * The Jacobian \(J_{xy}\) (i.e., Jacobian of unknown \(x\) w.r.t. the parameter + * \(y\)) is calculated given the solution as follows. Since + * \[ + * f(x, y) = 0, + * \] + * we have (\(J_{pq}\) being the Jacobian matrix \(\tfrac {dq} {dq}\)) + * \[ + * - J_{fx} J_{xy} = J_{fy}, + * \] + * and therefore \(J_{xy}\) can be solved from system + * \[ + * - J_{fx} J_{xy} = J_{fy}. + * \] + * Let \(eta\) be the adjoint with respect to \(x\); then to calculate + * \[ + * \eta J_{xy}, + * \] + * we solve + * \[ + * - (\eta J_{fx}^{-1}) J_{fy}. + * \] + * + * @tparam F type of equation system function + * @tparam T type of elements in the x vector + * @tparam Args types of additional parameters 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] 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. + * @pre f returns finite values when passed any value of x and the given args. + * @throw std::invalid_argument if x has size zero. + * @throw std::invalid_argument if x has non-finite elements. + * elements. + * @throw std::invalid_argument if relative_tolerance is strictly + * negative. + * @throw std::invalid_argument if function_tolerance is strictly + * negative. + * @throw std::invalid_argument if max_num_steps is not positive. + * @throw std::domain_error solver exceeds max_num_steps. + * @throw std::domain_error if the norm of the solution exceeds + * the function tolerance. + */ +template * = nullptr, + require_any_st_var* = nullptr> +Eigen::Matrix 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) { + auto x_ref = eval(value_of(x)); + auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(args)...)); + auto args_vals_tuple = apply( + [&](const auto&... args) { + return std::make_tuple(to_ref(value_of(args))...); + }, + *arena_args_tuple); + + auto f_wrt_x = [&args_vals_tuple, &f, msgs](const auto& x) { + return apply( + [&x, &f, msgs](const auto&... args) { return f(x, msgs, args...); }, + args_vals_tuple); + }; + + check_nonzero_size("algebra_solver_powell", "initial guess", x_ref); + check_finite("algebra_solver_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", + 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); + + Eigen::MatrixXd Jf_x; + Eigen::VectorXd f_x; + + jacobian(f_wrt_x, x_ref, f_x, Jf_x); + + using ret_type = Eigen::Matrix; + auto Jf_x_T_lu_ptr + = make_unsafe_chainable_ptr(Jf_x.transpose().partialPivLu()); // Lu + + arena_t ret = x_ref; + + reverse_pass_callback([f, ret, arena_args_tuple, Jf_x_T_lu_ptr, + msgs]() mutable { + Eigen::VectorXd eta = -Jf_x_T_lu_ptr->solve(ret.adj().eval()); + + // Contract with Jacobian of f with respect to y using a nested reverse + // autodiff pass. + { + nested_rev_autodiff rev; + Eigen::VectorXd ret_val = ret.val(); + auto x_nrad_ = apply( + [&](const auto&... args) { return eval(f(ret_val, msgs, args...)); }, + *arena_args_tuple); + x_nrad_.adj() = eta; + grad(); + } + }); + + return ret_type(ret); +} + } // namespace math } // namespace stan diff --git a/stan/math/rev/functor/algebra_system.hpp b/stan/math/rev/functor/algebra_system.hpp index f165ae56624..76b32ffe77f 100644 --- a/stan/math/rev/functor/algebra_system.hpp +++ b/stan/math/rev/functor/algebra_system.hpp @@ -12,59 +12,6 @@ namespace stan { namespace math { -/** - * A functor that allows us to treat either x or y as - * the independent variable. If x_is_iv = true, than the - * Jacobian is computed w.r.t x, else it is computed - * w.r.t y. - * - * @tparam F type for algebraic system functor - * @tparam T0 type for unknowns - * @tparam T1 type for auxiliary parameters - * @tparam x_is_iv true if x is the independent variable - */ -template -struct system_functor { - /** algebraic system functor */ - F f_; - /** unknowns */ - Eigen::Matrix x_; - /** auxiliary parameters */ - Eigen::Matrix y_; - /** real data */ - std::vector dat_; - /** integer data */ - std::vector dat_int_; - /** stream message */ - std::ostream* msgs_; - - system_functor() {} - - system_functor(const F& f, const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& dat, - const std::vector& dat_int, std::ostream* msgs) - : f_(f), x_(x), y_(y), dat_(dat), dat_int_(dat_int), msgs_(msgs) {} - - /** - * An operator that takes in an independent variable. The - * independent variable is either passed as the unknown x, - * or the auxiliary parameter y. The x_is_iv template parameter - * allows us to determine whether the jacobian is computed - * with respect to x or y. - * @tparam T the scalar type of the independent variable - */ - template - inline Eigen::Matrix operator()( - const Eigen::Matrix& iv) const { - if (x_is_iv) { - return f_(iv, y_, dat_, dat_int_, msgs_); - } else { - return f_(x_, iv, dat_, dat_int_, msgs_); - } - } -}; - /** * A structure which gets passed to Eigen's dogleg * algebraic solver. @@ -86,10 +33,7 @@ struct nlo_functor { }; /** - * A functor with the required operators to call Eigen's - * algebraic solver. - * It is also used in the vari classes of the algebraic solvers - * to compute the requisite sensitivities. + * A functor with the required operators to call Eigen's algebraic solver. * * @tparam S wrapper around the algebraic system functor. Has the * signature required for jacobian (i.e takes only one argument). @@ -97,26 +41,20 @@ struct nlo_functor { * @tparam T0 scalar type for unknowns * @tparam T1 scalar type for auxiliary parameters */ -template +template struct hybrj_functor_solver : nlo_functor { /** Wrapper around algebraic system */ S fs_; - /** number of unknowns */ - int x_size_; + /** Jacobian of algebraic function wrt unknowns */ Eigen::MatrixXd J_; - hybrj_functor_solver(const S& fs, const F& f, - const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& dat, - const std::vector& dat_int, std::ostream* msgs) - : fs_(f, x, y, dat, dat_int, msgs), x_size_(x.size()) {} + explicit hybrj_functor_solver(const S& fs) : fs_(fs) {} /** * Computes the value the algebraic function, f, when pluging in the - * independent variables, and the Jacobian w.r.t unknowns. Required - * by Eigen. + * independent variables, and the Jacobian w.r.t unknowns. + * * @param [in] iv independent variables * @param [in, out] fvec value of algebraic function when plugging in iv. */ @@ -126,8 +64,8 @@ struct hybrj_functor_solver : nlo_functor { } /** - * Assign the Jacobian to fjac (signature required by Eigen). Required - * by Eigen. + * Assign the Jacobian to fjac. + * * @param [in] iv independent variables. * @param [in, out] fjac matrix container for jacobian */ @@ -157,6 +95,8 @@ struct hybrj_functor_solver : nlo_functor { Eigen::VectorXd get_value(const Eigen::VectorXd& iv) const { return fs_(iv); } }; +// TODO(jgaeb): Remove this when the chain method of the fixed point solver is +// updated. template void algebra_solver_check(const Eigen::Matrix& x, const Eigen::Matrix y, diff --git a/stan/math/rev/functor/kinsol_data.hpp b/stan/math/rev/functor/kinsol_data.hpp index 5dd9b288338..8557c29dace 100644 --- a/stan/math/rev/functor/kinsol_data.hpp +++ b/stan/math/rev/functor/kinsol_data.hpp @@ -5,40 +5,17 @@ #include #include #include +#include #include #include #include #include #include +#include namespace stan { namespace math { -/** - * Default Jacobian builder using reverse-mode autodiff. - */ -struct kinsol_J_f { - template - inline int operator()(const F& f, const Eigen::VectorXd& x, - const Eigen::VectorXd& y, - const std::vector& dat, - const std::vector& dat_int, std::ostream* msgs, - const double x_sun[], SUNMatrix J) const { - size_t N = x.size(); - const std::vector x_vec(x_sun, x_sun + N); - system_functor system(f, x, y, dat, dat_int, msgs); - Eigen::VectorXd fx; - Eigen::MatrixXd Jac; - jacobian(system, to_vector(x_vec), fx, Jac); - - for (int j = 0; j < Jac.cols(); j++) - for (int i = 0; i < Jac.rows(); i++) - SM_ELEMENT_D(J, i, j) = Jac(i, j); - - return 0; - } -}; - /** * KINSOL algebraic system data holder. * Based on cvodes_ode_data. @@ -47,18 +24,15 @@ struct kinsol_J_f { * @tparam F2 functor type for jacobian function. Default is 0. * If 0, use rev mode autodiff to compute the Jacobian. */ -template +template class kinsol_system_data { const F1& f_; - const F2& J_f_; const Eigen::VectorXd& x_; - const Eigen::VectorXd& y_; const size_t N_; - const std::vector& dat_; - const std::vector& dat_int_; - std::ostream* msgs_; + std::ostream* const msgs_; + const std::tuple args_tuple_; - typedef kinsol_system_data system_data; + typedef kinsol_system_data system_data; public: N_Vector nv_x_; @@ -67,17 +41,13 @@ class kinsol_system_data { void* kinsol_memory_; /* Constructor */ - kinsol_system_data(const F1& f, const F2& J_f, const Eigen::VectorXd& x, - const Eigen::VectorXd& y, const std::vector& dat, - const std::vector& dat_int, std::ostream* msgs) + kinsol_system_data(const F1& f, const Eigen::VectorXd& x, + std::ostream* const msgs, const Args&... args) : f_(f), - J_f_(J_f), x_(x), - y_(y), N_(x.size()), - dat_(dat), - dat_int_(dat_int), msgs_(msgs), + args_tuple_(args...), nv_x_(N_VMake_Serial(N_, &to_array_1d(x_)[0])), J_(SUNDenseMatrix(N_, N_)), LS_(SUNLinSol_Dense(nv_x_, J_)), @@ -91,18 +61,24 @@ class kinsol_system_data { } /* Implements the user-defined function passed to KINSOL. */ - static int kinsol_f_system(N_Vector x, N_Vector f, void* user_data) { + static int kinsol_f_system(const N_Vector x, const N_Vector f_eval, + void* const user_data) { const system_data* explicit_system = static_cast(user_data); Eigen::VectorXd x_eigen( Eigen::Map(NV_DATA_S(x), explicit_system->N_)); - Eigen::Map(N_VGetArrayPointer(f), explicit_system->N_) - = explicit_system->f_(x_eigen, explicit_system->y_, - explicit_system->dat_, explicit_system->dat_int_, - explicit_system->msgs_); - + Eigen::Map f_eval_map(N_VGetArrayPointer(f_eval), + explicit_system->N_); + auto result = apply( + [&](const auto&... args) { + return explicit_system->f_(x_eigen, explicit_system->msgs_, args...); + }, + explicit_system->args_tuple_); + check_matching_sizes("", "the algebraic system's output", result, + "the vector of unknowns, x,", f_eval_map); + f_eval_map = result; return 0; } @@ -118,14 +94,37 @@ class kinsol_system_data { * https://computation.llnl.gov/sites/default/files/public/kin_guide-dev.pdf, * page 55. */ - static int kinsol_jacobian(N_Vector x, N_Vector f, SUNMatrix J, - void* user_data, N_Vector tmp1, N_Vector tmp2) { + static int kinsol_jacobian(const N_Vector x, const N_Vector f_eval, + const SUNMatrix J, void* const user_data, + const N_Vector tmp1, const N_Vector tmp2) { const system_data* explicit_system = static_cast(user_data); - return explicit_system->J_f_(explicit_system->f_, explicit_system->x_, - explicit_system->y_, explicit_system->dat_, - explicit_system->dat_int_, - explicit_system->msgs_, NV_DATA_S(x), J); + + Eigen::VectorXd x_eigen( + Eigen::Map(NV_DATA_S(x), explicit_system->N_)); + Eigen::Map f_eval_map(N_VGetArrayPointer(f_eval), + explicit_system->N_); + + auto f_wrt_x = [&](const auto& x) { + return apply( + [&](const auto&... args) { + return explicit_system->f_(x, explicit_system->msgs_, args...); + }, + explicit_system->args_tuple_); + }; + + Eigen::MatrixXd Jf_x; + Eigen::VectorXd f_x; + + jacobian(f_wrt_x, x_eigen, f_x, Jf_x); + + f_eval_map = f_x; + + for (int j = 0; j < Jf_x.cols(); j++) + for (int i = 0; i < Jf_x.rows(); i++) + SM_ELEMENT_D(J, i, j) = Jf_x(i, j); + + return 0; } }; diff --git a/stan/math/rev/functor/kinsol_solve.hpp b/stan/math/rev/functor/kinsol_solve.hpp index 3da9c825938..8d03912adf6 100644 --- a/stan/math/rev/functor/kinsol_solve.hpp +++ b/stan/math/rev/functor/kinsol_solve.hpp @@ -55,18 +55,18 @@ namespace math { * @throw std::runtime_error if Kinsol returns a * negative flag that is not due to hitting max_num_steps. */ -template -Eigen::VectorXd kinsol_solve( - const F1& f, const Eigen::VectorXd& x, const Eigen::VectorXd& y, - const std::vector& dat, const std::vector& dat_int, - std::ostream* msgs = nullptr, double scaling_step_tol = 1e-3, - double function_tolerance = 1e-6, - long int max_num_steps = 200, // NOLINT(runtime/int) - bool custom_jacobian = 1, const F2& J_f = kinsol_J_f(), - int steps_eval_jacobian = 10, int global_line_search = KIN_LINESEARCH) { +template +Eigen::VectorXd kinsol_solve(const F1& f, const Eigen::VectorXd& x, + const double scaling_step_tol, // = 1e-3 + const double function_tolerance, // = 1e-6 + const int64_t max_num_steps, // = 200 + const bool custom_jacobian, // = 1 + const int steps_eval_jacobian, // = 10 + const int global_line_search, // = KIN_LINESEARCH + std::ostream* const msgs, const Args&... args) { int N = x.size(); - typedef kinsol_system_data system_data; - system_data kinsol_data(f, J_f, x, y, dat, dat_int, msgs); + typedef kinsol_system_data system_data; + system_data kinsol_data(f, x, msgs, args...); check_flag_sundials(KINInit(kinsol_data.kinsol_memory_, &system_data::kinsol_f_system, kinsol_data.nv_x_), diff --git a/test/expressions/expression_test_helpers.hpp b/test/expressions/expression_test_helpers.hpp index 7528f2ef833..0c6dc72ad83 100644 --- a/test/expressions/expression_test_helpers.hpp +++ b/test/expressions/expression_test_helpers.hpp @@ -195,13 +195,13 @@ struct test_functor { }; struct simple_eq_functor { - template - inline Eigen::Matrix, Eigen::Dynamic, 1> - operator()(const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& dat, const std::vector& dat_int, + template * = nullptr, + require_all_eigen_vector_t* = nullptr> + inline Eigen::Matrix, Eigen::Dynamic, 1> + operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int, std::ostream* pstream__) const { - Eigen::Matrix, Eigen::Dynamic, 1> z(1); + Eigen::Matrix, Eigen::Dynamic, 1> z(1); z(0) = x(0) - y(0); return z; } diff --git a/test/unit/math/rev/core/chainable_object_test.cpp b/test/unit/math/rev/core/chainable_object_test.cpp index 5523b05ed63..fddd02fd3cf 100644 --- a/test/unit/math/rev/core/chainable_object_test.cpp +++ b/test/unit/math/rev/core/chainable_object_test.cpp @@ -66,3 +66,70 @@ TEST(AgradRev, make_chainable_ptr_nested_test) { EXPECT_EQ((ChainableObjectTest::counter), 1); } + +class UnsafeChainableObjectTest { + public: + static int counter; + + ~UnsafeChainableObjectTest() { counter++; } +}; + +int UnsafeChainableObjectTest::counter = 0; + +TEST(AgradRev, unsafe_chainable_object_test) { + { + auto ptr + = new stan::math::unsafe_chainable_object( + UnsafeChainableObjectTest()); + UnsafeChainableObjectTest::counter = 0; + } + + EXPECT_EQ((UnsafeChainableObjectTest::counter), 0); + stan::math::recover_memory(); + EXPECT_EQ((UnsafeChainableObjectTest::counter), 1); +} + +TEST(AgradRev, unsafe_chainable_object_nested_test) { + stan::math::start_nested(); + + { + auto ptr + = new stan::math::unsafe_chainable_object( + UnsafeChainableObjectTest()); + UnsafeChainableObjectTest::counter = 0; + } + + EXPECT_EQ((UnsafeChainableObjectTest::counter), 0); + + stan::math::recover_memory_nested(); + + EXPECT_EQ((UnsafeChainableObjectTest::counter), 1); +} + +TEST(AgradRev, make_unsafe_chainable_ptr_test) { + { + UnsafeChainableObjectTest* ptr + = stan::math::make_unsafe_chainable_ptr(UnsafeChainableObjectTest()); + UnsafeChainableObjectTest::counter = 0; + } + + EXPECT_EQ((UnsafeChainableObjectTest::counter), 0); + stan::math::recover_memory(); + EXPECT_EQ((UnsafeChainableObjectTest::counter), 1); +} + +TEST(AgradRev, make_unsafe_chainable_ptr_nested_test) { + stan::math::start_nested(); + + { + UnsafeChainableObjectTest* ptr + = stan::math::make_unsafe_chainable_ptr(UnsafeChainableObjectTest()); + UnsafeChainableObjectTest::counter = 0; + } + + EXPECT_EQ((UnsafeChainableObjectTest::counter), 0); + + stan::math::recover_memory_nested(); + + EXPECT_EQ((UnsafeChainableObjectTest::counter), 1); +} diff --git a/test/unit/math/rev/functor/algebra_solver_fp_test.cpp b/test/unit/math/rev/functor/algebra_solver_fp_test.cpp index 4e30a9b92fe..67b742968bf 100644 --- a/test/unit/math/rev/functor/algebra_solver_fp_test.cpp +++ b/test/unit/math/rev/functor/algebra_solver_fp_test.cpp @@ -217,11 +217,9 @@ struct FP_direct_prod_func_test : public ::testing::Test { * RHS functor */ struct FP_direct_prod_func { - template + template inline Eigen::Matrix, -1, 1> operator()( - const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& x_r, const std::vector& x_i, + const T0& x, const T1& y, const T2& x_r, const T3& x_i, std::ostream* pstream__) const { using scalar = stan::return_type_t; const size_t n = x.size(); @@ -248,11 +246,9 @@ struct FP_direct_prod_func_test : public ::testing::Test { * Newton root functor */ struct FP_direct_prod_newton_func { - template + template inline Eigen::Matrix, -1, 1> operator()( - const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& x_r, const std::vector& x_i, + const T0& x, const T1& y, const T2& x_r, const T3& x_i, std::ostream* pstream__) const { using scalar = stan::return_type_t; const size_t n = x.size(); @@ -515,6 +511,8 @@ TEST_F(FP_direct_prod_func_test, algebra_solver_fp) { EXPECT_NEAR(g_newton[j], g_fp[j], 1.e-8); } } + + stan::math::recover_memory(); } TEST_F(FP_2d_func_test, exception_handling) { diff --git a/test/unit/math/rev/functor/algebra_solver_newton_test.cpp b/test/unit/math/rev/functor/algebra_solver_newton_test.cpp new file mode 100644 index 00000000000..c88db8a8914 --- /dev/null +++ b/test/unit/math/rev/functor/algebra_solver_newton_test.cpp @@ -0,0 +1,259 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +////////////////////////////////////////////////////////////////////////// +// Tests for newton solver. + +TEST_F(algebra_solver_simple_eq_test, newton_dbl) { + bool is_newton = true; + Eigen::VectorXd theta = simple_eq_test(simple_eq_functor(), y_dbl, is_newton); +} + +TEST_F(algebra_solver_simple_eq_test, newton_tuned_dbl) { + bool is_newton = true; + Eigen::VectorXd theta = simple_eq_test(simple_eq_functor(), y_dbl, is_newton, + true, scale_step, xtol, ftol, maxfev); +} + +TEST_F(algebra_solver_simple_eq_nopara_test, newton_dbl) { + using stan::math::algebra_solver_newton; + Eigen::VectorXd theta = algebra_solver_newton(simple_eq_functor_nopara(), x, + y_dummy, dat, dummy_dat_int); + EXPECT_EQ(20, theta(0)); + EXPECT_EQ(2, theta(1)); +} + +TEST_F(algebra_solver_non_linear_eq_test, newton_dbl) { + bool is_newton = true; + Eigen::VectorXd theta + = non_linear_eq_test(non_linear_eq_functor(), y_dbl, is_newton); + EXPECT_FLOAT_EQ(-y_dbl(0), theta(0)); + EXPECT_FLOAT_EQ(-y_dbl(1), theta(1)); + EXPECT_FLOAT_EQ(y_dbl(2), theta(2)); +} + +TEST_F(error_message_test, newton_dbl) { + bool is_newton = true; + error_conditions_test(non_linear_eq_functor(), y_3, is_newton); +} + +TEST_F(max_steps_test, newton_dbl) { + bool is_newton = true; + max_num_steps_test(y, is_newton); +} + +TEST(MathMatrixRevMat, unsolvable_flag_newton_dbl) { + Eigen::VectorXd y(2); + y << 1, 1; + + unsolvable_flag_test(y); +} + +TEST_F(degenerate_eq_test, newton_guess1_dbl) { + using stan::math::algebra_solver_newton; + + // This first initial guess produces the + // solution x = {8, 8} + + Eigen::VectorXd theta = algebra_solver_newton(degenerate_eq_functor(), + x_guess_1, y_dbl, dat, dat_int); + EXPECT_FLOAT_EQ(8, theta(0)); + EXPECT_FLOAT_EQ(8, theta(1)); +} + +TEST_F(degenerate_eq_test, newton_guess2_dbl) { + using stan::math::algebra_solver_newton; + // This next initial guess produces the + // solution x = {5, 5} + + Eigen::VectorXd theta = algebra_solver_newton(degenerate_eq_functor(), + x_guess_2, y_dbl, dat, dat_int); + EXPECT_FLOAT_EQ(5, theta(0)); + EXPECT_FLOAT_EQ(5, theta(1)); +} + +// For the next two unit tests, see if the initial +// guess determines neighborhood of the +// solution, when solutions have different scales, +// using y_scale. + +TEST_F(degenerate_eq_test, newton_guess2_scale_dbl) { + using stan::math::algebra_solver_newton; + + Eigen::VectorXd theta = algebra_solver_newton( + degenerate_eq_functor(), x_guess_2, y_scale, dat, dat_int); + EXPECT_FLOAT_EQ(5, theta(0)); + EXPECT_FLOAT_EQ(5, theta(1)); +} + +TEST_F(degenerate_eq_test, newton_guess_saddle_point_dbl) { + // Newton solver fails this test because the initial point is + // a saddle point. + using stan::math::algebra_solver_newton; + std::stringstream err_msg; + err_msg << "algebra_solver failed with error flag -11."; + std::string msg = err_msg.str(); + + EXPECT_THROW_MSG(algebra_solver_newton(degenerate_eq_functor(), x_guess_3, + y_scale, dat, dat_int), + std::runtime_error, msg); +} + +TEST_F(algebra_solver_simple_eq_test, newton) { + using stan::math::var; + bool is_newton = true; + for (int k = 0; k < n_x; k++) { + Eigen::Matrix y = y_dbl; + + Eigen::Matrix theta + = simple_eq_test(simple_eq_functor(), y, is_newton); + + std::vector y_vec{y(0), y(1), y(2)}; + std::vector g; + theta(k).grad(y_vec, g); + + for (int i = 0; i < n_y; i++) + EXPECT_EQ(J(k, i), g[i]); + } +} + +TEST_F(algebra_solver_simple_eq_test, newton_tuned) { + using stan::math::var; + bool is_newton = true; + for (int k = 0; k < n_x; k++) { + Eigen::Matrix y = y_dbl; + + Eigen::Matrix theta + = simple_eq_test(simple_eq_functor(), y, is_newton, true, scale_step, + xtol, ftol, maxfev); + + std::vector y_vec{y(0), y(1), y(2)}; + std::vector g; + theta(k).grad(y_vec, g); + + for (int i = 0; i < n_y; i++) + EXPECT_EQ(J(k, i), g[i]); + } +} + +TEST_F(algebra_solver_simple_eq_test, newton_init_is_para) { + using stan::math::algebra_solver_newton; + Eigen::VectorXd theta + = algebra_solver_newton(simple_eq_functor(), x_var, y_dbl, dat, dat_int); + EXPECT_EQ(20, theta(0)); + EXPECT_EQ(2, theta(1)); +} + +TEST_F(algebra_solver_non_linear_eq_test, newton) { + using stan::math::var; + bool is_newton = true; + for (int k = 0; k < n_x; k++) { + Eigen::Matrix y = y_dbl; + Eigen::Matrix theta + = non_linear_eq_test(non_linear_eq_functor(), y, is_newton); + + EXPECT_FLOAT_EQ(-y(0).val(), theta(0).val()); + EXPECT_FLOAT_EQ(-y(1).val(), theta(1).val()); + EXPECT_FLOAT_EQ(y(2).val(), theta(2).val()); + + std::vector y_vec{y(0), y(1), y(2)}; + std::vector g; + theta(k).grad(y_vec, g); + + for (int i = 0; i < n_y; i++) + EXPECT_NEAR(J(k, i), g[i], err); + } +} + +TEST_F(error_message_test, newton) { + using stan::math::var; + bool is_newton = true; + Eigen::Matrix y = y_2; + error_conditions_test(non_linear_eq_functor(), y, is_newton); +} + +TEST_F(max_steps_test, newton) { + bool is_newton = true; + max_num_steps_test(y_var, is_newton); +} + +TEST(MathMatrixRevMat, unsolvable_flag_newton) { + Eigen::Matrix y(2); + y << 1, 1; + + unsolvable_flag_test(y); +} + +TEST_F(degenerate_eq_test, newton_guess1) { + using stan::math::algebra_solver_newton; + // using stan::math::sum; + using stan::math::var; + + // This first initial guess produces the + // solution x = {8, 8} + for (int k = 0; k < n_x; k++) { + Eigen::Matrix y = y_dbl; + Eigen::Matrix theta = algebra_solver_newton( + degenerate_eq_functor(), x_guess_1, y, dat, dat_int); + EXPECT_FLOAT_EQ(8, theta(0).val()); + EXPECT_FLOAT_EQ(8, theta(1).val()); + + std::vector y_vec{y(0), y(1)}; + std::vector g; + theta(k).grad(y_vec, g); + + for (int l = 0; l < n_y; l++) + EXPECT_NEAR(J1(k, l), g[l], tolerance); + } +} + +TEST_F(degenerate_eq_test, newton_guess2) { + using stan::math::algebra_solver_newton; + using stan::math::var; + // This next initial guess produces the + // solution x = {5, 5} + for (int k = 0; k < 1; k++) { + Eigen::Matrix y = y_dbl; + Eigen::Matrix theta = algebra_solver_newton( + degenerate_eq_functor(), x_guess_2, y, dat, dat_int); + EXPECT_FLOAT_EQ(5, theta(0).val()); + EXPECT_FLOAT_EQ(5, theta(0).val()); + + std::vector y_vec{y(0), y(1)}; + std::vector g; + theta(k).grad(y_vec, g); + + for (int l = 0; l < n_y; l++) + EXPECT_NEAR(J2(k, l), g[l], tolerance); + } +} + +TEST_F(variadic_test, newton) { + using stan::math::var; + bool is_newton = true; + for (int k = 0; k < n_x; k++) { + var y_1 = y_1_dbl; + var y_2 = y_2_dbl; + var y_3 = y_3_dbl; + + Eigen::Matrix theta + = variadic_eq_test(variadic_eq_functor(), A, y_1, y_2, y_3, i, + is_newton, scaling_step_size, relative_tolerance, + function_tolerance, max_num_steps); + + std::vector y_vec{y_1, y_2, y_3}; + std::vector g; + theta(k).grad(y_vec, g); + + for (int i = 0; i < n_y; i++) + EXPECT_NEAR(J(k, i), g[i], 1e-6); + } +} diff --git a/test/unit/math/rev/functor/algebra_solver_test.cpp b/test/unit/math/rev/functor/algebra_solver_test.cpp index 3bda620c702..481aa94bd4d 100644 --- a/test/unit/math/rev/functor/algebra_solver_test.cpp +++ b/test/unit/math/rev/functor/algebra_solver_test.cpp @@ -1,6 +1,5 @@ #include #include -#include #include #include #include @@ -10,160 +9,9 @@ #include #include -// Every test exists in four versions for the cases -// where y (the auxiliary parameters) are passed as -// data (double type) or parameters (var types), -// and the cases where the solver is based on Powell's -// or Newton's method. - -class algebra_solver_simple_eq_test : public ::testing::Test { - protected: - void SetUp() override { - n_x = 2; - n_y = 3; - - y_dbl = stan::math::to_vector({5, 4, 2}); - Eigen::MatrixXd J_(n_x, n_y); - J_ << 4, 5, 0, 0, 0, 1; - J = J_; - - x_var = stan::math::to_vector({1, 1}); - } - - int n_x; - int n_y; - Eigen::VectorXd y_dbl; - Eigen::Matrix x_var; - std::vector dat; - std::vector dat_int; - double scale_step = 1e-3; - double xtol = 1e-6; - double ftol = 1e-3; - int maxfev = 1e+2; - - Eigen::MatrixXd J; -}; - -class algebra_solver_simple_eq_nopara_test : public ::testing::Test { - protected: - void SetUp() override { x = stan::math::to_vector({1, 1}); } - - int n_x = 2; - Eigen::VectorXd x; - std::vector dat = {5, 4, 2}; - Eigen::VectorXd y_dummy; - std::vector dummy_dat_int; -}; - -class algebra_solver_non_linear_eq_test : public ::testing::Test { - protected: - void SetUp() override { - y_dbl = stan::math::to_vector({4, 6, 3}); - Eigen::MatrixXd J_(n_x, n_y); - J_ << -1, 0, 0, 0, -1, 0, 0, 0, 1; - J = J_; - } - int n_x = 3; - int n_y = 3; - double err = 1e-11; - - Eigen::VectorXd y_dbl; - Eigen::MatrixXd J; -}; - -class error_message_test : public ::testing::Test { - protected: - void SetUp() override { - y_2 = stan::math::to_vector({4, 6}); - y_3 = stan::math::to_vector({4, 6, 3}); - } - - Eigen::VectorXd y_2; - Eigen::VectorXd y_3; -}; - -class max_steps_test : public ::testing::Test { - protected: - void SetUp() override { - y = stan::math::to_vector({1, 1, 1}); - y_var = stan::math::to_vector({1, 1, 1}); - } - - Eigen::VectorXd y; - Eigen::Matrix y_var; -}; - -class degenerate_eq_test : public ::testing::Test { - protected: - void SetUp() override { - using stan::math::to_vector; - y_dbl = to_vector({5, 8}); - y_scale = to_vector({5, 100}); - x_guess_1 = to_vector({10, 1}); - x_guess_2 = to_vector({1, 1}); - x_guess_3 = to_vector({5, 100}); // 80, 80 - - Eigen::MatrixXd J_(n_x, n_y); - J_ << 0, 1, 0, 1; - J1 = J_; - J_ << 1, 0, 1, 0; - J2 = J_; - } - - int n_x = 2; - int n_y = 2; - double tolerance = 1e-10; - Eigen::VectorXd y_dbl; - Eigen::VectorXd y_scale; - Eigen::VectorXd x_guess_1; - Eigen::VectorXd x_guess_2; - Eigen::VectorXd x_guess_3; - Eigen::MatrixXd J1; - Eigen::MatrixXd J2; - std::vector dat; - std::vector dat_int; -}; - ////////////////////////////////////////////////////////////////////////// // Tests for powell solver. -TEST_F(algebra_solver_simple_eq_test, powell) { - using stan::math::var; - bool is_newton = false; - for (int k = 0; k < n_x; k++) { - Eigen::Matrix y = y_dbl; - - Eigen::Matrix theta - = simple_eq_test(simple_eq_functor(), y, is_newton); - - std::vector y_vec{y(0), y(1), y(2)}; - std::vector g; - theta(k).grad(y_vec, g); - - for (int i = 0; i < n_y; i++) - EXPECT_EQ(J(k, i), g[i]); - } -} - -TEST_F(algebra_solver_simple_eq_test, powell_tuned) { - using stan::math::var; - bool is_newton = false; - for (int k = 0; k < n_x; k++) { - Eigen::Matrix y = y_dbl; - - Eigen::Matrix theta - = simple_eq_test(simple_eq_functor(), y, is_newton, true, scale_step, - xtol, ftol, maxfev); - - std::vector y_vec{y(0), y(1), y(2)}; - std::vector g; - theta(k).grad(y_vec, g); - - for (int i = 0; i < n_y; i++) - EXPECT_EQ(J(k, i), g[i]); - } -} - TEST_F(algebra_solver_simple_eq_test, powell_dbl) { bool is_newton = false; Eigen::VectorXd theta = simple_eq_test(simple_eq_functor(), y_dbl, is_newton); @@ -183,35 +31,6 @@ TEST_F(algebra_solver_simple_eq_nopara_test, powell) { EXPECT_EQ(2, theta(1)); } -TEST_F(algebra_solver_simple_eq_test, powell_init_is_para) { - using stan::math::algebra_solver_powell; - Eigen::VectorXd theta - = algebra_solver_powell(simple_eq_functor(), x_var, y_dbl, dat, dat_int); - EXPECT_EQ(20, theta(0)); - EXPECT_EQ(2, theta(1)); -} - -TEST_F(algebra_solver_non_linear_eq_test, powell) { - using stan::math::var; - bool is_newton = false; - for (int k = 0; k < n_x; k++) { - Eigen::Matrix y = y_dbl; - Eigen::Matrix theta - = non_linear_eq_test(non_linear_eq_functor(), y, is_newton); - - EXPECT_FLOAT_EQ(-y(0).val(), theta(0).val()); - EXPECT_FLOAT_EQ(-y(1).val(), theta(1).val()); - EXPECT_FLOAT_EQ(y(2).val(), theta(2).val()); - - std::vector y_vec{y(0), y(1), y(2)}; - std::vector g; - theta(k).grad(y_vec, g); - - for (int i = 0; i < n_y; i++) - EXPECT_NEAR(J(k, i), g[i], err); - } -} - TEST_F(algebra_solver_non_linear_eq_test, powell_dbl) { bool is_newton = false; Eigen::VectorXd theta @@ -221,11 +40,12 @@ TEST_F(algebra_solver_non_linear_eq_test, powell_dbl) { EXPECT_FLOAT_EQ(y_dbl(2), theta(2)); } -TEST_F(error_message_test, powell) { - using stan::math::var; - bool is_newton = false; - Eigen::Matrix y = y_2; - error_conditions_test(non_linear_eq_functor(), y, is_newton); +TEST_F(algebra_solver_simple_eq_nopara_test, powell_double) { + using stan::math::algebra_solver_powell; + Eigen::VectorXd theta = algebra_solver_powell(simple_eq_functor_nopara(), x, + y_dummy, dat, dummy_dat_int); + EXPECT_EQ(20, theta(0)); + EXPECT_EQ(2, theta(1)); } TEST_F(error_message_test, powell_dbl) { @@ -233,72 +53,17 @@ TEST_F(error_message_test, powell_dbl) { error_conditions_test(non_linear_eq_functor(), y_3, is_newton); } -TEST(unsolvable_test, powell) { - using stan::math::var; - Eigen::Matrix y(2); - y << 1, 1; - unsolvable_test(y); -} - TEST(unsolvable_test, powell_dbl) { Eigen::VectorXd y(2); y << 1, 1; unsolvable_test(y); } -TEST_F(max_steps_test, powell) { - bool is_newton = false; - max_num_steps_test(y_var, is_newton); -} - TEST_F(max_steps_test, powell_dbl) { bool is_newton = false; max_num_steps_test(y, is_newton); } -TEST_F(degenerate_eq_test, powell_guess1) { - using stan::math::algebra_solver_powell; - using stan::math::var; - - // This first initial guess produces the - // solution x = {8, 8} - for (int k = 0; k < n_x; k++) { - Eigen::Matrix y = y_dbl; - Eigen::Matrix theta = algebra_solver_powell( - degenerate_eq_functor(), x_guess_1, y, dat, dat_int); - EXPECT_FLOAT_EQ(8, theta(0).val()); - EXPECT_FLOAT_EQ(8, theta(1).val()); - - std::vector y_vec{y(0), y(1)}; - std::vector g; - theta(k).grad(y_vec, g); - - for (int l = 0; l < n_y; l++) - EXPECT_NEAR(J1(k, l), g[l], tolerance); - } -} - -TEST_F(degenerate_eq_test, powell_guess2) { - using stan::math::algebra_solver_powell; - using stan::math::var; - // This next initial guess produces the - // solution x = {5, 5} - for (int k = 0; k < 1; k++) { - Eigen::Matrix y = y_dbl; - Eigen::Matrix theta = algebra_solver_powell( - degenerate_eq_functor(), x_guess_2, y, dat, dat_int); - EXPECT_FLOAT_EQ(5, theta(0).val()); - EXPECT_FLOAT_EQ(5, theta(0).val()); - - std::vector y_vec{y(0), y(1)}; - std::vector g; - theta(k).grad(y_vec, g); - - for (int l = 0; l < n_y; l++) - EXPECT_NEAR(J2(k, l), g[l], tolerance); - } -} - TEST_F(degenerate_eq_test, powell_guess1_dbl) { using stan::math::algebra_solver_powell; @@ -345,30 +110,9 @@ TEST_F(degenerate_eq_test, powell_guess_saddle_point_dbl) { EXPECT_FLOAT_EQ(100, theta(1)); } -// unit test to demo issue #696 -// system functor init bug issue #696 -TEST(MathMatrixRevMat, system_functor_constructor) { - using stan::math::system_functor; - - Eigen::VectorXd y(2); - y << 5, 8; - Eigen::VectorXd x(2); - x << 10, 1; - std::vector dat{0.0, 0.0}; - std::vector dat_int{0, 0}; - std::ostream* msgs = 0; - int f = 99; - - system_functor fs(f, x, y, dat, dat_int, msgs); - - EXPECT_EQ(fs.f_, f); -} -////////////////////////////////////////////////////////////////////////// -// Tests for newton solver. - -TEST_F(algebra_solver_simple_eq_test, newton) { +TEST_F(algebra_solver_simple_eq_test, powell) { using stan::math::var; - bool is_newton = true; + bool is_newton = false; for (int k = 0; k < n_x; k++) { Eigen::Matrix y = y_dbl; @@ -379,14 +123,15 @@ TEST_F(algebra_solver_simple_eq_test, newton) { std::vector g; theta(k).grad(y_vec, g); - for (int i = 0; i < n_y; i++) + for (int i = 0; i < n_y; i++) { EXPECT_EQ(J(k, i), g[i]); + } } } -TEST_F(algebra_solver_simple_eq_test, newton_tuned) { +TEST_F(algebra_solver_simple_eq_test, powell_tuned) { using stan::math::var; - bool is_newton = true; + bool is_newton = false; for (int k = 0; k < n_x; k++) { Eigen::Matrix y = y_dbl; @@ -403,36 +148,17 @@ TEST_F(algebra_solver_simple_eq_test, newton_tuned) { } } -TEST_F(algebra_solver_simple_eq_test, newton_dbl) { - bool is_newton = true; - Eigen::VectorXd theta = simple_eq_test(simple_eq_functor(), y_dbl, is_newton); -} - -TEST_F(algebra_solver_simple_eq_test, newton_tuned_dbl) { - bool is_newton = true; - Eigen::VectorXd theta = simple_eq_test(simple_eq_functor(), y_dbl, is_newton, - true, scale_step, xtol, ftol, maxfev); -} - -TEST_F(algebra_solver_simple_eq_nopara_test, newton) { - using stan::math::algebra_solver_newton; - Eigen::VectorXd theta = algebra_solver_newton(simple_eq_functor_nopara(), x, - y_dummy, dat, dummy_dat_int); - EXPECT_EQ(20, theta(0)); - EXPECT_EQ(2, theta(1)); -} - -TEST_F(algebra_solver_simple_eq_test, newton_init_is_para) { - using stan::math::algebra_solver_newton; +TEST_F(algebra_solver_simple_eq_test, powell_init_is_para) { + using stan::math::algebra_solver_powell; Eigen::VectorXd theta - = algebra_solver_newton(simple_eq_functor(), x_var, y_dbl, dat, dat_int); + = algebra_solver_powell(simple_eq_functor(), x_var, y_dbl, dat, dat_int); EXPECT_EQ(20, theta(0)); EXPECT_EQ(2, theta(1)); } -TEST_F(algebra_solver_non_linear_eq_test, newton) { +TEST_F(algebra_solver_non_linear_eq_test, powell) { using stan::math::var; - bool is_newton = true; + bool is_newton = false; for (int k = 0; k < n_x; k++) { Eigen::Matrix y = y_dbl; Eigen::Matrix theta @@ -451,61 +177,34 @@ TEST_F(algebra_solver_non_linear_eq_test, newton) { } } -TEST_F(algebra_solver_non_linear_eq_test, newton_dbl) { - bool is_newton = true; - Eigen::VectorXd theta - = non_linear_eq_test(non_linear_eq_functor(), y_dbl, is_newton); - EXPECT_FLOAT_EQ(-y_dbl(0), theta(0)); - EXPECT_FLOAT_EQ(-y_dbl(1), theta(1)); - EXPECT_FLOAT_EQ(y_dbl(2), theta(2)); -} - -TEST_F(error_message_test, newton) { +TEST_F(error_message_test, powell) { using stan::math::var; - bool is_newton = true; + bool is_newton = false; Eigen::Matrix y = y_2; error_conditions_test(non_linear_eq_functor(), y, is_newton); } -TEST_F(error_message_test, newton_dbl) { - bool is_newton = true; - error_conditions_test(non_linear_eq_functor(), y_3, is_newton); -} - -TEST_F(max_steps_test, newton) { - bool is_newton = true; - max_num_steps_test(y_var, is_newton); -} - -TEST_F(max_steps_test, newton_dbl) { - bool is_newton = true; - max_num_steps_test(y, is_newton); -} - -TEST(MathMatrixRevMat, unsolvable_flag_newton) { - Eigen::Matrix y(2); +TEST(unsolvable_test, powell) { + using stan::math::var; + Eigen::Matrix y(2); y << 1, 1; - - unsolvable_flag_test(y); + unsolvable_test(y); } -TEST(MathMatrixRevMat, unsolvable_flag_newton_dbl) { - Eigen::VectorXd y(2); - y << 1, 1; - - unsolvable_flag_test(y); +TEST_F(max_steps_test, powell) { + bool is_newton = false; + max_num_steps_test(y_var, is_newton); } -TEST_F(degenerate_eq_test, newton_guess1) { - using stan::math::algebra_solver_newton; - // using stan::math::sum; +TEST_F(degenerate_eq_test, powell_guess1) { + using stan::math::algebra_solver_powell; using stan::math::var; // This first initial guess produces the // solution x = {8, 8} for (int k = 0; k < n_x; k++) { Eigen::Matrix y = y_dbl; - Eigen::Matrix theta = algebra_solver_newton( + Eigen::Matrix theta = algebra_solver_powell( degenerate_eq_functor(), x_guess_1, y, dat, dat_int); EXPECT_FLOAT_EQ(8, theta(0).val()); EXPECT_FLOAT_EQ(8, theta(1).val()); @@ -519,14 +218,14 @@ TEST_F(degenerate_eq_test, newton_guess1) { } } -TEST_F(degenerate_eq_test, newton_guess2) { - using stan::math::algebra_solver_newton; +TEST_F(degenerate_eq_test, powell_guess2) { + using stan::math::algebra_solver_powell; using stan::math::var; // This next initial guess produces the // solution x = {5, 5} for (int k = 0; k < 1; k++) { Eigen::Matrix y = y_dbl; - Eigen::Matrix theta = algebra_solver_newton( + Eigen::Matrix theta = algebra_solver_powell( degenerate_eq_functor(), x_guess_2, y, dat, dat_int); EXPECT_FLOAT_EQ(5, theta(0).val()); EXPECT_FLOAT_EQ(5, theta(0).val()); @@ -540,52 +239,24 @@ TEST_F(degenerate_eq_test, newton_guess2) { } } -TEST_F(degenerate_eq_test, newton_guess1_dbl) { - using stan::math::algebra_solver_newton; - - // This first initial guess produces the - // solution x = {8, 8} - - Eigen::VectorXd theta = algebra_solver_newton(degenerate_eq_functor(), - x_guess_1, y_dbl, dat, dat_int); - EXPECT_FLOAT_EQ(8, theta(0)); - EXPECT_FLOAT_EQ(8, theta(1)); -} - -TEST_F(degenerate_eq_test, newton_guess2_dbl) { - using stan::math::algebra_solver_newton; - // This next initial guess produces the - // solution x = {5, 5} - - Eigen::VectorXd theta = algebra_solver_newton(degenerate_eq_functor(), - x_guess_2, y_dbl, dat, dat_int); - EXPECT_FLOAT_EQ(5, theta(0)); - EXPECT_FLOAT_EQ(5, theta(1)); -} - -// For the next two unit tests, see if the initial -// guess determines neighborhood of the -// solution, when solutions have different scales, -// using y_scale. +TEST_F(variadic_test, powell) { + using stan::math::var; + bool is_newton = true; + for (int k = 0; k < n_x; k++) { + var y_1 = y_1_dbl; + var y_2 = y_2_dbl; + var y_3 = y_3_dbl; -TEST_F(degenerate_eq_test, newton_guess2_scale_dbl) { - using stan::math::algebra_solver_newton; + Eigen::Matrix theta + = variadic_eq_test(variadic_eq_functor(), A, y_1, y_2, y_3, i, + is_newton, scaling_step_size, relative_tolerance, + function_tolerance, max_num_steps); - Eigen::VectorXd theta = algebra_solver_newton( - degenerate_eq_functor(), x_guess_2, y_scale, dat, dat_int); - EXPECT_FLOAT_EQ(5, theta(0)); - EXPECT_FLOAT_EQ(5, theta(1)); -} + std::vector y_vec{y_1, y_2, y_3}; + std::vector g; + theta(k).grad(y_vec, g); -TEST_F(degenerate_eq_test, newton_guess_saddle_point_dbl) { - // Newton solver fails this test because the initial point is - // a saddle point. - using stan::math::algebra_solver_newton; - std::stringstream err_msg; - err_msg << "algebra_solver failed with error flag -11."; - std::string msg = err_msg.str(); - - EXPECT_THROW_MSG(algebra_solver_newton(degenerate_eq_functor(), x_guess_3, - y_scale, dat, dat_int), - std::runtime_error, msg); + for (int i = 0; i < n_y; i++) + EXPECT_NEAR(J(k, i), g[i], 1e-6); + } } diff --git a/test/unit/math/rev/functor/util_algebra_solver.hpp b/test/unit/math/rev/functor/util_algebra_solver.hpp index 1ae2c426ce1..63833170bfb 100644 --- a/test/unit/math/rev/functor/util_algebra_solver.hpp +++ b/test/unit/math/rev/functor/util_algebra_solver.hpp @@ -8,6 +8,176 @@ #include #include +// Every test exists in four versions for the cases +// where y (the auxiliary parameters) are passed as +// data (double type) or parameters (var types), +// and the cases where the solver is based on Powell's +// or Newton's method. + +class algebra_solver_simple_eq_test : public ::testing::Test { + protected: + void SetUp() override { + n_x = 2; + n_y = 3; + + y_dbl = stan::math::to_vector({5, 4, 2}); + Eigen::MatrixXd J_(n_x, n_y); + J_ << 4, 5, 0, 0, 0, 1; + J = J_; + + x_var = stan::math::to_vector({1, 1}); + } + + void TearDown() override { stan::math::recover_memory(); } + + int n_x; + int n_y; + Eigen::VectorXd y_dbl; + Eigen::Matrix x_var; + std::vector dat; + std::vector dat_int; + double scale_step = 1e-3; + double xtol = 1e-6; + double ftol = 1e-3; + int maxfev = 1e+2; + + Eigen::MatrixXd J; +}; + +class algebra_solver_simple_eq_nopara_test : public ::testing::Test { + protected: + void SetUp() override { x = stan::math::to_vector({1, 1}); } + + void TearDown() override { stan::math::recover_memory(); } + + int n_x = 2; + Eigen::VectorXd x; + std::vector dat = {5, 4, 2}; + Eigen::VectorXd y_dummy; + std::vector dummy_dat_int; +}; + +class algebra_solver_non_linear_eq_test : public ::testing::Test { + protected: + void SetUp() override { + y_dbl = stan::math::to_vector({4, 6, 3}); + Eigen::MatrixXd J_(n_x, n_y); + J_ << -1, 0, 0, 0, -1, 0, 0, 0, 1; + J = J_; + } + + void TearDown() override { stan::math::recover_memory(); } + + int n_x = 3; + int n_y = 3; + double err = 1e-11; + + Eigen::VectorXd y_dbl; + Eigen::MatrixXd J; +}; + +class error_message_test : public ::testing::Test { + protected: + void SetUp() override { + y_2 = stan::math::to_vector({4, 6}); + y_3 = stan::math::to_vector({4, 6, 3}); + } + + void TearDown() override { stan::math::recover_memory(); } + + Eigen::VectorXd y_2; + Eigen::VectorXd y_3; +}; + +class max_steps_test : public ::testing::Test { + protected: + void SetUp() override { + y = stan::math::to_vector({1, 1, 1}); + y_var = stan::math::to_vector({1, 1, 1}); + } + + void TearDown() override { stan::math::recover_memory(); } + + Eigen::VectorXd y; + Eigen::Matrix y_var; +}; + +class degenerate_eq_test : public ::testing::Test { + protected: + void SetUp() override { + using stan::math::to_vector; + y_dbl = to_vector({5, 8}); + y_scale = to_vector({5, 100}); + x_guess_1 = to_vector({10, 1}); + x_guess_2 = to_vector({1, 1}); + x_guess_3 = to_vector({5, 100}); // 80, 80 + + Eigen::MatrixXd J_(n_x, n_y); + J_ << 0, 1, 0, 1; + J1 = J_; + J_ << 1, 0, 1, 0; + J2 = J_; + } + + void TearDown() override { stan::math::recover_memory(); } + + int n_x = 2; + int n_y = 2; + double tolerance = 1e-10; + Eigen::VectorXd y_dbl; + Eigen::VectorXd y_scale; + Eigen::VectorXd x_guess_1; + Eigen::VectorXd x_guess_2; + Eigen::VectorXd x_guess_3; + Eigen::MatrixXd J1; + Eigen::MatrixXd J2; + std::vector dat; + std::vector dat_int; +}; + +class variadic_test : public ::testing::Test { + protected: + void SetUp() override { + n_x = 2; + n_y = 3; + + y_1_dbl = 5; + y_2_dbl = 4; + y_3_dbl = 2; + + Eigen::MatrixXd A_(n_x, n_x); + A_ << 1, 2, 2, 1; + A = A_; + + x_var = stan::math::to_vector({1, 3}); + + Eigen::MatrixXd J_(n_x, n_y); + J_ << 4, 5, 0, 0, 0, 1; + J = J_; + + i = 3; + } + + void TearDown() override { stan::math::recover_memory(); } + + int n_x; + int n_y; + Eigen::MatrixXd A_; + Eigen::Matrix A; + double y_1_dbl; + double y_2_dbl; + double y_3_dbl; + int i; + Eigen::Matrix x_var; + + double scaling_step_size = 1e-3; + double relative_tolerance = 1e-6; + double function_tolerance = 1e-3; + int max_num_steps = 1e+2; + + Eigen::MatrixXd J; +}; + /* wrapper function that either calls the dogleg or the newton solver. */ template Eigen::Matrix general_algebra_solver( @@ -33,13 +203,11 @@ Eigen::Matrix general_algebra_solver( /* define algebraic functions which get solved. */ struct simple_eq_functor { - template - inline Eigen::Matrix, Eigen::Dynamic, 1> - operator()(const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& dat, const std::vector& dat_int, + template + inline Eigen::Matrix, Eigen::Dynamic, 1> + operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int, std::ostream* pstream__) const { - Eigen::Matrix, Eigen::Dynamic, 1> z(2); + Eigen::Matrix, Eigen::Dynamic, 1> z(2); z(0) = x(0) - y(0) * y(1); z(1) = x(1) - y(2); return z; @@ -47,13 +215,11 @@ struct simple_eq_functor { }; struct simple_eq_functor_nopara { - template - inline Eigen::Matrix, Eigen::Dynamic, 1> - operator()(const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& dat, const std::vector& dat_int, + template + inline Eigen::Matrix, Eigen::Dynamic, 1> + operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int, std::ostream* pstream__) const { - Eigen::Matrix, Eigen::Dynamic, 1> z(2); + Eigen::Matrix, Eigen::Dynamic, 1> z(2); z(0) = x(0) - dat[0] * dat[1]; z(1) = x(1) - dat[2]; return z; @@ -61,13 +227,11 @@ struct simple_eq_functor_nopara { }; struct non_linear_eq_functor { - template - inline Eigen::Matrix, Eigen::Dynamic, 1> - operator()(const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& dat, const std::vector& dat_int, + template + inline Eigen::Matrix, Eigen::Dynamic, 1> + operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int, std::ostream* pstream__) const { - Eigen::Matrix, Eigen::Dynamic, 1> z(3); + Eigen::Matrix, Eigen::Dynamic, 1> z(3); z(0) = x(2) - y(2); z(1) = x(0) * x(1) - y(0) * y(1); z(2) = x(2) / x(0) + y(2) / y(0); @@ -76,13 +240,11 @@ struct non_linear_eq_functor { }; struct non_square_eq_functor { - template - inline Eigen::Matrix, Eigen::Dynamic, 1> - operator()(const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& dat, const std::vector& dat_int, + template + inline Eigen::Matrix, Eigen::Dynamic, 1> + operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int, std::ostream* pstream__) const { - Eigen::Matrix, Eigen::Dynamic, 1> z(2); + Eigen::Matrix, Eigen::Dynamic, 1> z(2); z(0) = x(0) - y(0); z(1) = x(1) * x(2) - y(1); return z; @@ -90,13 +252,11 @@ struct non_square_eq_functor { }; struct unsolvable_eq_functor { - template - inline Eigen::Matrix, Eigen::Dynamic, 1> - operator()(const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& dat, const std::vector& dat_int, + template + inline Eigen::Matrix, Eigen::Dynamic, 1> + operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int, std::ostream* pstream__) const { - Eigen::Matrix, Eigen::Dynamic, 1> z(2); + Eigen::Matrix, Eigen::Dynamic, 1> z(2); z(0) = x(0) * x(0) + y(0); z(1) = x(1) * x(1) + y(1); return z; @@ -105,19 +265,31 @@ struct unsolvable_eq_functor { // Degenerate roots: each solution can either be y(0) or y(1). struct degenerate_eq_functor { - template - inline Eigen::Matrix, Eigen::Dynamic, 1> - operator()(const Eigen::Matrix& x, - const Eigen::Matrix& y, - const std::vector& dat, const std::vector& dat_int, + template + inline Eigen::Matrix, Eigen::Dynamic, 1> + operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int, std::ostream* pstream__) const { - Eigen::Matrix, Eigen::Dynamic, 1> z(2); + Eigen::Matrix, Eigen::Dynamic, 1> z(2); z(0) = (x(0) - y(0)) * (x(1) - y(1)); z(1) = (x(0) - y(1)) * (x(1) - y(0)); return z; } }; +struct variadic_eq_functor { + template + inline Eigen::Matrix, Eigen::Dynamic, 1> + operator()(const T1& x, std::ostream* pstream__, const T2& A, const T3& y_1, + const T3& y_2, const T3& y_3, const T4& i) const { + Eigen::Matrix, Eigen::Dynamic, 1> z(2); + z(0) = x(0) - y_1 * y_2; + z(1) = x(1) - y_3; + for (int j = 0; j < i; ++j) + z = A * z; + return z; + } +}; + /* template code for running tests in the prim and rev regime */ template @@ -180,16 +352,16 @@ inline void error_conditions_test(const F& f, std::vector dat_int; std::stringstream err_msg; - err_msg << "algebra_solver: size of the algebraic system's output (2) " + err_msg << "size of the algebraic system's output (2) " << "and size of the vector of unknowns, x, (3) must match in size"; std::string msg = err_msg.str(); - EXPECT_THROW_MSG( - algebra_solver_powell(non_square_eq_functor(), x, y, dat, dat_int), - std::invalid_argument, msg); + EXPECT_THROW_MSG(general_algebra_solver(is_newton, non_square_eq_functor(), x, + y, dat, dat_int), + std::invalid_argument, msg); Eigen::VectorXd x_bad(static_cast(0)); std::stringstream err_msg2; - err_msg2 << "algebra_solver: initial guess has size 0, but " + err_msg2 << "initial guess has size 0, but " << "must have a non-zero size"; std::string msg2 = err_msg2.str(); EXPECT_THROW_MSG(general_algebra_solver(is_newton, f, x_bad, y, dat, dat_int), @@ -200,26 +372,7 @@ inline void error_conditions_test(const F& f, x_bad_inf << inf, 1, 1; EXPECT_THROW_MSG( general_algebra_solver(is_newton, f, x_bad_inf, y, dat, dat_int), - std::domain_error, - "algebra_solver: initial guess[1] is inf, but must " - "be finite!"); - - typedef Eigen::Matrix matrix; - matrix y_bad_inf(3); - y_bad_inf << inf, 1, 1; - EXPECT_THROW_MSG( - general_algebra_solver(is_newton, f, x, y_bad_inf, dat, dat_int), - std::domain_error, - "algebra_solver: parameter vector[1] is inf, but must " - "be finite!"); - - std::vector dat_bad_inf(1); - dat_bad_inf[0] = inf; - EXPECT_THROW_MSG( - general_algebra_solver(is_newton, f, x, y, dat_bad_inf, dat_int), - std::domain_error, - "algebra_solver: continuous data[1] is inf, but must " - "be finite!"); + std::domain_error, "initial guess[1] is inf, but must be finite!"); if (!is_newton) { EXPECT_THROW_MSG(general_algebra_solver(is_newton, f, x, y, dat, dat_int, 0, @@ -302,3 +455,35 @@ inline void max_num_steps_test(Eigen::Matrix& y, function_tolerance, max_num_steps), std::domain_error, msg); } + +template +Eigen::Matrix variadic_eq_test( + const F& f, + const Eigen::Matrix A, + const stan::math::var& y_1, const stan::math::var& y_2, + const stan::math::var& y_3, const int i, bool is_newton = false, + double scaling_step_size = 1e-3, double relative_tolerance = 1e-10, + double function_tolerance = 1e-6, int32_t max_num_steps = 1e+3) { + using stan::math::algebra_solver_newton; + using stan::math::algebra_solver_powell; + using stan::math::var; + + int n_x = 2; + Eigen::VectorXd x(2); + x << 1, 1; // initial guess + + Eigen::Matrix theta; + + theta = is_newton + ? algebra_solver_newton_impl( + variadic_eq_functor(), x, &std::cout, scaling_step_size, + function_tolerance, max_num_steps, A, y_1, y_2, y_3, i) + : algebra_solver_powell_impl( + variadic_eq_functor(), x, &std::cout, relative_tolerance, + function_tolerance, max_num_steps, A, y_1, y_2, y_3, i); + + EXPECT_NEAR(20, value_of(theta(0)), 1e-6); + EXPECT_NEAR(2, value_of(theta(1)), 1e-6); + + return theta; +}