Skip to content

Commit

Permalink
Updated newton test
Browse files Browse the repository at this point in the history
  • Loading branch information
bbbales2 committed Jul 11, 2021
1 parent d7d9356 commit 8d864cb
Show file tree
Hide file tree
Showing 10 changed files with 414 additions and 431 deletions.
4 changes: 1 addition & 3 deletions stan/math/prim/functor/algebra_solver_fp.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -99,9 +99,7 @@ struct KinsolFixedPointEnv {
};

/**
* Private interface for solving fixed point problems using KINSOL. Users should
* call the KINSOL fixed point solver through `algebra_solver_fp` or
* `algebra_solver_fp_impl`.
* Solve algebraic equations using KINSOL fixed point solver
*
* @tparam F type of the equation system functor f
* @tparam T_u type of scaling vector for unknowns. We allow
Expand Down
16 changes: 4 additions & 12 deletions stan/math/prim/functor/algebra_solver_powell.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,15 +18,13 @@ namespace stan {
namespace math {

/**
* Private interface for calling the Powell solver. Users should call the Powell
* solver through `algebra_solver_powell` or `algebra_solver_powell_impl`.
* Solve algebraic equations using Powell solver
*
* @tparam F type of equation system function, curried with respect to inputs
* @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, curried with
* respect to its input (i.e., "y") values.
* @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
Expand Down Expand Up @@ -95,9 +93,7 @@ T& algebra_solver_powell_call_solver(
* (xtol in Eigen's code), the function tolerance,
* and the maximum number of steps (maxfev in Eigen's code).
*
* This function is overloaded to handle both constant and var-type parameters.
* This overload handles non-var parameters, and checks the input and calls the
* algebraic solver only.
* This overload handles non-autodiff parameters.
*
* @tparam F type of equation system function
* @tparam T type of elements in the x vector
Expand Down Expand Up @@ -137,7 +133,6 @@ Eigen::VectorXd algebra_solver_powell_impl(const F& f, const T& x,
const auto& x_ref = to_ref(x);
auto x_val = to_ref(value_of(x_ref));

// Curry the input function w.r.t. y
auto f_wrt_x = [&f, msgs, &args...](const auto& x) { return f(x, msgs, args...); };

check_nonzero_size("algebra_solver_powell", "initial guess", x_val);
Expand Down Expand Up @@ -166,9 +161,6 @@ Eigen::VectorXd algebra_solver_powell_impl(const F& f, const T& x,
* (xtol in Eigen's code), the function tolerance,
* and the maximum number of steps (maxfev in Eigen's code).
*
* Signature to maintain backward compatibility, will be removed
* in the future.
*
* @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
Expand Down
7 changes: 2 additions & 5 deletions stan/math/rev/functor/algebra_solver_fp.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -36,10 +36,7 @@ namespace math {
* The user can also specify the scaling controls, the function
* tolerance, and the maximum number of steps.
*
* This function is overloaded to handle both constant and var-type parameters.
* This overload handles var parameters, and checks the input, calls the
* algebraic solver, and appropriately handles derivative propagation through
* the `reverse_pass_callback`.
* 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
Expand Down Expand Up @@ -139,7 +136,7 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_fp_impl(
arena_t<ret_type> ret = theta_dbl;

auto Jf_xT_lu_ptr
= make_unsafe_chainable_ptr((Eigen::MatrixXd::Identity(x.size(), x.size()) - Jf_x).transpose().eval().partialPivLu()); // Lu
= make_unsafe_chainable_ptr((Eigen::MatrixXd::Identity(x.size(), x.size()) - Jf_x).transpose().eval().partialPivLu());

reverse_pass_callback([f, ret, arena_args_tuple, Jf_xT_lu_ptr, msgs]() mutable {
// Contract specificities with inverse Jacobian of f with respect to x.
Expand Down
10 changes: 2 additions & 8 deletions stan/math/rev/functor/algebra_solver_newton.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -24,9 +24,7 @@ namespace math {
* The user can also specify the scaled step size, the function
* tolerance, and the maximum number of steps.
*
* This function is overloaded to handle both constant and var-type parameters.
* This overload handles non-var parameters, and checks the input and calls the
* algebraic solver only.
* This overload handles non-autodiff parameters.
*
* @tparam F type of equation system function.
* @tparam T type of initial guess vector.
Expand Down Expand Up @@ -85,10 +83,7 @@ Eigen::VectorXd algebra_solver_newton_impl(const F& f, const T& x,
* The user can also specify the scaled step size, the function
* tolerance, and the maximum number of steps.
*
* This function is overloaded to handle both constant and var-type parameters.
* This overload handles var parameters, and checks the input, calls the
* algebraic solver, and appropriately handles derivative propagation through
* the `reverse_pass_callback`.
* 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
Expand Down Expand Up @@ -169,7 +164,6 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
},
args_vals_tuple);

// Evaluate and store the Jacobian.
auto f_wrt_x = [&](const auto& x) {
return apply([&](const auto&... args) { return f(x, msgs, args...); },
args_vals_tuple);
Expand Down
6 changes: 1 addition & 5 deletions stan/math/rev/functor/algebra_solver_powell.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -28,10 +28,7 @@ namespace math {
* (xtol in Eigen's code), the function tolerance,
* and the maximum number of steps (maxfev in Eigen's code).
*
* This function is overloaded to handle both constant and var-type parameters.
* This overload handles var parameters, and checks the input, calls the
* algebraic solver, and appropriately handles derivative propagation through
* the `reverse_pass_callback`.
* 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
Expand Down Expand Up @@ -97,7 +94,6 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_powell_impl(
},
arena_args_tuple);

// Curry the input function w.r.t. y
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);
Expand Down
13 changes: 5 additions & 8 deletions stan/math/rev/functor/algebra_system.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -33,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).
Expand All @@ -56,8 +53,8 @@ struct hybrj_functor_solver : nlo_functor<double> {

/**
* 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.
*/
Expand All @@ -67,8 +64,8 @@ struct hybrj_functor_solver : nlo_functor<double> {
}

/**
* 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
*/
Expand Down
1 change: 1 addition & 0 deletions stan/math/rev/functor/kinsol_data.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include <stan/math/rev/functor/jacobian.hpp>
#include <stan/math/prim/fun/to_array_1d.hpp>
#include <stan/math/prim/fun/to_vector.hpp>
#include <stan/math/prim/functor/apply.hpp>
#include <kinsol/kinsol.h>
#include <sunmatrix/sunmatrix_dense.h>
#include <sunlinsol/sunlinsol_dense.h>
Expand Down
Loading

0 comments on commit 8d864cb

Please sign in to comment.