Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Revert "Feature/issue 2401 alg solver adjoint" #2563

Merged
merged 1 commit into from
Aug 24, 2021
+808 −1,191
Merged

Revert "Feature/issue 2401 alg solver adjoint" #2563

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 0 additions & 26 deletions stan/math/prim/functor/algebra_solver_adapter.hpp
View file

This file was deleted.

55 changes: 0 additions & 55 deletions stan/math/rev/core/chainable_object.hpp
View file
Original file line number Diff line number Diff line change
@@ -62,61 +62,6 @@ 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 <typename T>
class unsafe_chainable_object : public chainable_alloc {
private:
std::decay_t<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 <typename S,
require_same_t<plain_type_t<T>, plain_type_t<S>>* = nullptr>
explicit unsafe_chainable_object(S&& obj) : obj_(std::forward<S>(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 <typename T>
auto make_unsafe_chainable_ptr(T&& obj) {
auto ptr = new unsafe_chainable_object<T>(std::forward<T>(obj));
return &ptr->get();
}

} // namespace math
} // namespace stan
#endif
230 changes: 63 additions & 167 deletions stan/math/rev/functor/algebra_solver_newton.hpp
View file
Original file line number Diff line number Diff line change
@@ -3,10 +3,11 @@

#include <stan/math/rev/core.hpp>
#include <stan/math/rev/functor/algebra_system.hpp>
#include <stan/math/rev/functor/algebra_solver_powell.hpp>
#include <stan/math/rev/functor/kinsol_solve.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/mdivide_left.hpp>
#include <stan/math/prim/fun/value_of.hpp>
#include <stan/math/prim/functor/algebra_solver_adapter.hpp>
#include <unsupported/Eigen/NonLinearOptimization>
#include <iostream>
#include <string>
@@ -24,180 +25,53 @@ 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, 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 <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... 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) {
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] 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
* 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 <code>std::invalid_argument</code> if x has size zero.
* * @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 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.
* @throw <code>std::domain_error</code> if solver exceeds max_num_steps.
*/
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) {
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<var, Eigen::Dynamic, -1>;
arena_t<ret_type> 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);
template <typename F, typename T, require_eigen_vector_t<T>* = nullptr>
Eigen::VectorXd algebra_solver_newton(
const F& f, const T& x, const Eigen::VectorXd& y,
const std::vector<double>& dat, const std::vector<int>& 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);
}

/**
@@ -209,15 +83,19 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
* The user can also specify the scaled step size, the function
* tolerance, and the maximum number of steps.
*
* Signature to maintain backward compatibility, will be removed
* in the future.
* 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 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.
* @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.
@@ -241,16 +119,34 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
* @throw <code>std::domain_error if solver exceeds max_num_steps.
*/
template <typename F, typename T1, typename T2,
require_all_eigen_vector_t<T1, T2>* = nullptr>
require_all_eigen_vector_t<T1, T2>* = nullptr,
require_st_var<T2>* = nullptr>
Eigen::Matrix<scalar_type_t<T2>, Eigen::Dynamic, 1> algebra_solver_newton(
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 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);
const std::vector<int>& 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<F, double, double, false> Fy;
typedef system_functor<F, double, double, true> Fs;
typedef hybrj_functor_solver<Fs, F, double, double> 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<Fy, F, scalar_type_t<T2>, Fx>(
Fy(), f, value_of(x_eval), y_eval, dat, dat_int, theta_dbl, fx, msgs);
Eigen::Matrix<var, Eigen::Dynamic, 1> 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;
}

} // namespace math
399 changes: 160 additions & 239 deletions stan/math/rev/functor/algebra_solver_powell.hpp
View file

Large diffs are not rendered by default.

80 changes: 70 additions & 10 deletions stan/math/rev/functor/algebra_system.hpp
View file
Original file line number Diff line number Diff line change
@@ -12,6 +12,59 @@
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 <typename F, typename T0, typename T1, bool x_is_iv>
struct system_functor {
/** algebraic system functor */
F f_;
/** unknowns */
Eigen::Matrix<T0, Eigen::Dynamic, 1> x_;
/** auxiliary parameters */
Eigen::Matrix<T1, Eigen::Dynamic, 1> y_;
/** real data */
std::vector<double> dat_;
/** integer data */
std::vector<int> dat_int_;
/** stream message */
std::ostream* msgs_;

system_functor() {}

system_functor(const F& f, const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& dat,
const std::vector<int>& 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 <typename T>
inline Eigen::Matrix<T, Eigen::Dynamic, 1> operator()(
const Eigen::Matrix<T, Eigen::Dynamic, 1>& 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.
@@ -33,28 +86,37 @@ struct nlo_functor {
};

/**
* A functor with the required operators to call Eigen's algebraic solver.
* 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.
*
* @tparam S wrapper around the algebraic system functor. Has the
* signature required for jacobian (i.e takes only one argument).
* @tparam F algebraic system functor
* @tparam T0 scalar type for unknowns
* @tparam T1 scalar type for auxiliary parameters
*/
template <typename S>
template <typename S, typename F, typename T0, typename T1>
struct hybrj_functor_solver : nlo_functor<double> {
/** Wrapper around algebraic system */
S fs_;

/** number of unknowns */
int x_size_;
/** Jacobian of algebraic function wrt unknowns */
Eigen::MatrixXd J_;

explicit hybrj_functor_solver(const S& fs) : fs_(fs) {}
hybrj_functor_solver(const S& fs, const F& f,
const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& dat,
const std::vector<int>& dat_int, std::ostream* msgs)
: fs_(f, x, y, dat, dat_int, msgs), x_size_(x.size()) {}

/**
* Computes the value the algebraic function, f, when pluging in the
* independent variables, and the Jacobian w.r.t unknowns.
*
* independent variables, and the Jacobian w.r.t unknowns. Required
* by Eigen.
* @param [in] iv independent variables
* @param [in, out] fvec value of algebraic function when plugging in iv.
*/
@@ -64,8 +126,8 @@ struct hybrj_functor_solver : nlo_functor<double> {
}

/**
* Assign the Jacobian to fjac.
*
* Assign the Jacobian to fjac (signature required by Eigen). Required
* by Eigen.
* @param [in] iv independent variables.
* @param [in, out] fjac matrix container for jacobian
*/
@@ -95,8 +157,6 @@ struct hybrj_functor_solver : nlo_functor<double> {
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 <typename T1, typename T2>
void algebra_solver_check(const Eigen::Matrix<T1, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T2, Eigen::Dynamic, 1> y,
101 changes: 51 additions & 50 deletions stan/math/rev/functor/kinsol_data.hpp
View file
Original file line number Diff line number Diff line change
@@ -5,17 +5,40 @@
#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>
#include <nvector/nvector_serial.h>
#include <vector>
#include <tuple>

namespace stan {
namespace math {

/**
* Default Jacobian builder using reverse-mode autodiff.
*/
struct kinsol_J_f {
template <typename F>
inline int operator()(const F& f, const Eigen::VectorXd& x,
const Eigen::VectorXd& y,
const std::vector<double>& dat,
const std::vector<int>& dat_int, std::ostream* msgs,
const double x_sun[], SUNMatrix J) const {
size_t N = x.size();
const std::vector<double> x_vec(x_sun, x_sun + N);
system_functor<F, double, double, 1> 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.
@@ -24,15 +47,18 @@ namespace math {
* @tparam F2 functor type for jacobian function. Default is 0.
* If 0, use rev mode autodiff to compute the Jacobian.
*/
template <typename F1, typename... Args>
template <typename F1, typename F2>
class kinsol_system_data {
const F1& f_;
const F2& J_f_;
const Eigen::VectorXd& x_;
const Eigen::VectorXd& y_;
const size_t N_;
std::ostream* const msgs_;
const std::tuple<const Args&...> args_tuple_;
const std::vector<double>& dat_;
const std::vector<int>& dat_int_;
std::ostream* msgs_;

typedef kinsol_system_data<F1, Args...> system_data;
typedef kinsol_system_data<F1, F2> system_data;

public:
N_Vector nv_x_;
@@ -41,13 +67,17 @@ class kinsol_system_data {
void* kinsol_memory_;

/* Constructor */
kinsol_system_data(const F1& f, const Eigen::VectorXd& x,
std::ostream* const msgs, const Args&... args)
kinsol_system_data(const F1& f, const F2& J_f, const Eigen::VectorXd& x,
const Eigen::VectorXd& y, const std::vector<double>& dat,
const std::vector<int>& dat_int, std::ostream* msgs)
: 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_)),
@@ -61,24 +91,18 @@ class kinsol_system_data {
}

/* Implements the user-defined function passed to KINSOL. */
static int kinsol_f_system(const N_Vector x, const N_Vector f_eval,
void* const user_data) {
static int kinsol_f_system(N_Vector x, N_Vector f, void* user_data) {
const system_data* explicit_system
= static_cast<const system_data*>(user_data);

Eigen::VectorXd x_eigen(
Eigen::Map<Eigen::VectorXd>(NV_DATA_S(x), explicit_system->N_));

Eigen::Map<Eigen::VectorXd> 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;
Eigen::Map<Eigen::VectorXd>(N_VGetArrayPointer(f), explicit_system->N_)
= explicit_system->f_(x_eigen, explicit_system->y_,
explicit_system->dat_, explicit_system->dat_int_,
explicit_system->msgs_);

return 0;
}

@@ -94,37 +118,14 @@ class kinsol_system_data {
* https://computation.llnl.gov/sites/default/files/public/kin_guide-dev.pdf,
* page 55.
*/
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) {
static int kinsol_jacobian(N_Vector x, N_Vector f, SUNMatrix J,
void* user_data, N_Vector tmp1, N_Vector tmp2) {
const system_data* explicit_system
= static_cast<const system_data*>(user_data);

Eigen::VectorXd x_eigen(
Eigen::Map<Eigen::VectorXd>(NV_DATA_S(x), explicit_system->N_));
Eigen::Map<Eigen::VectorXd> 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;
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);
}
};

22 changes: 11 additions & 11 deletions stan/math/rev/functor/kinsol_solve.hpp
View file
Original file line number Diff line number Diff line change
@@ -55,18 +55,18 @@ namespace math {
* @throw <code>std::runtime_error</code> if Kinsol returns a
* negative flag that is not due to hitting max_num_steps.
*/
template <typename F1, typename... Args>
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) {
template <typename F1, typename F2 = kinsol_J_f>
Eigen::VectorXd kinsol_solve(
const F1& f, const Eigen::VectorXd& x, const Eigen::VectorXd& y,
const std::vector<double>& dat, const std::vector<int>& 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) {
int N = x.size();
typedef kinsol_system_data<F1, Args...> system_data;
system_data kinsol_data(f, x, msgs, args...);
typedef kinsol_system_data<F1, F2> system_data;
system_data kinsol_data(f, J_f, x, y, dat, dat_int, msgs);

check_flag_sundials(KINInit(kinsol_data.kinsol_memory_,
&system_data::kinsol_f_system, kinsol_data.nv_x_),
12 changes: 6 additions & 6 deletions test/expressions/expression_test_helpers.hpp
View file
Original file line number Diff line number Diff line change
@@ -195,13 +195,13 @@ struct test_functor {
};

struct simple_eq_functor {
template <typename T1, typename T2, typename T3, typename T4,
require_eigen_vector_t<T1>* = nullptr,
require_all_eigen_vector_t<T1, T2>* = nullptr>
inline Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1>
operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int,
template <typename T0, typename T1>
inline Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1>
operator()(const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& dat, const std::vector<int>& dat_int,
std::ostream* pstream__) const {
Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1> z(1);
Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1> z(1);
z(0) = x(0) - y(0);
return z;
}
67 changes: 0 additions & 67 deletions test/unit/math/rev/core/chainable_object_test.cpp
View file
Original file line number Diff line number Diff line change
@@ -66,70 +66,3 @@ 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());
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());
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);
}
14 changes: 8 additions & 6 deletions test/unit/math/rev/functor/algebra_solver_fp_test.cpp
View file
Original file line number Diff line number Diff line change
@@ -217,9 +217,11 @@ struct FP_direct_prod_func_test : public ::testing::Test {
* RHS functor
*/
struct FP_direct_prod_func {
template <typename T0, typename T1, typename T2, typename T3>
template <typename T0, typename T1>
inline Eigen::Matrix<stan::return_type_t<T0, T1>, -1, 1> operator()(
const T0& x, const T1& y, const T2& x_r, const T3& x_i,
const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& x_r, const std::vector<int>& x_i,
std::ostream* pstream__) const {
using scalar = stan::return_type_t<T0, T1>;
const size_t n = x.size();
@@ -246,9 +248,11 @@ struct FP_direct_prod_func_test : public ::testing::Test {
* Newton root functor
*/
struct FP_direct_prod_newton_func {
template <typename T0, typename T1, typename T2, typename T3>
template <typename T0, typename T1>
inline Eigen::Matrix<stan::return_type_t<T0, T1>, -1, 1> operator()(
const T0& x, const T1& y, const T2& x_r, const T3& x_i,
const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& x_r, const std::vector<int>& x_i,
std::ostream* pstream__) const {
using scalar = stan::return_type_t<T0, T1>;
const size_t n = x.size();
@@ -511,8 +515,6 @@ 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) {
259 changes: 0 additions & 259 deletions test/unit/math/rev/functor/algebra_solver_newton_test.cpp
View file

This file was deleted.

427 changes: 378 additions & 49 deletions test/unit/math/rev/functor/algebra_solver_test.cpp
View file

Large diffs are not rendered by default.

307 changes: 61 additions & 246 deletions test/unit/math/rev/functor/util_algebra_solver.hpp
View file
Original file line number Diff line number Diff line change
@@ -8,176 +8,6 @@
#include <limits>
#include <string>

// 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<stan::math::var, Eigen::Dynamic, 1> x_var;
std::vector<double> dat;
std::vector<int> 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<double> dat = {5, 4, 2};
Eigen::VectorXd y_dummy;
std::vector<int> 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<stan::math::var, Eigen::Dynamic, 1> 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<double> dat;
std::vector<int> 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<stan::math::var, Eigen::Dynamic, Eigen::Dynamic> A;
double y_1_dbl;
double y_2_dbl;
double y_3_dbl;
int i;
Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1> 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 <typename F, typename T1, typename T2>
Eigen::Matrix<T2, Eigen::Dynamic, 1> general_algebra_solver(
@@ -203,35 +33,41 @@ Eigen::Matrix<T2, Eigen::Dynamic, 1> general_algebra_solver(
/* define algebraic functions which get solved. */

struct simple_eq_functor {
template <typename T1, typename T2, typename T3, typename T4>
inline Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1>
operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int,
template <typename T0, typename T1>
inline Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1>
operator()(const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& dat, const std::vector<int>& dat_int,
std::ostream* pstream__) const {
Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1> z(2);
Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1> z(2);
z(0) = x(0) - y(0) * y(1);
z(1) = x(1) - y(2);
return z;
}
};

struct simple_eq_functor_nopara {
template <typename T1, typename T2, typename T3, typename T4>
inline Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1>
operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int,
template <typename T0, typename T1>
inline Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1>
operator()(const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& dat, const std::vector<int>& dat_int,
std::ostream* pstream__) const {
Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1> z(2);
Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1> z(2);
z(0) = x(0) - dat[0] * dat[1];
z(1) = x(1) - dat[2];
return z;
}
};

struct non_linear_eq_functor {
template <typename T1, typename T2, typename T3, typename T4>
inline Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1>
operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int,
template <typename T0, typename T1>
inline Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1>
operator()(const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& dat, const std::vector<int>& dat_int,
std::ostream* pstream__) const {
Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1> z(3);
Eigen::Matrix<stan::return_type_t<T0, T1>, 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);
@@ -240,23 +76,27 @@ struct non_linear_eq_functor {
};

struct non_square_eq_functor {
template <typename T1, typename T2, typename T3, typename T4>
inline Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1>
operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int,
template <typename T0, typename T1>
inline Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1>
operator()(const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& dat, const std::vector<int>& dat_int,
std::ostream* pstream__) const {
Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1> z(2);
Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1> z(2);
z(0) = x(0) - y(0);
z(1) = x(1) * x(2) - y(1);
return z;
}
};

struct unsolvable_eq_functor {
template <typename T1, typename T2, typename T3, typename T4>
inline Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1>
operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int,
template <typename T0, typename T1>
inline Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1>
operator()(const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& dat, const std::vector<int>& dat_int,
std::ostream* pstream__) const {
Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1> z(2);
Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1> z(2);
z(0) = x(0) * x(0) + y(0);
z(1) = x(1) * x(1) + y(1);
return z;
@@ -265,31 +105,19 @@ struct unsolvable_eq_functor {

// Degenerate roots: each solution can either be y(0) or y(1).
struct degenerate_eq_functor {
template <typename T1, typename T2, typename T3, typename T4>
inline Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1>
operator()(const T1& x, const T2& y, const T3& dat, const T4& dat_int,
template <typename T0, typename T1>
inline Eigen::Matrix<stan::return_type_t<T0, T1>, Eigen::Dynamic, 1>
operator()(const Eigen::Matrix<T0, Eigen::Dynamic, 1>& x,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
const std::vector<double>& dat, const std::vector<int>& dat_int,
std::ostream* pstream__) const {
Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1> z(2);
Eigen::Matrix<stan::return_type_t<T0, T1>, 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 <typename T1, typename T2, typename T3, typename T4>
inline Eigen::Matrix<stan::return_type_t<T1, T2, T3, T4>, 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<stan::return_type_t<T1, T2, T3, T4>, 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 <typename F, typename T>
@@ -352,16 +180,16 @@ inline void error_conditions_test(const F& f,
std::vector<int> dat_int;

std::stringstream err_msg;
err_msg << "size of the algebraic system's output (2) "
err_msg << "algebra_solver: 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(general_algebra_solver(is_newton, non_square_eq_functor(), x,
y, dat, dat_int),
std::invalid_argument, msg);
EXPECT_THROW_MSG(
algebra_solver_powell(non_square_eq_functor(), x, y, dat, dat_int),
std::invalid_argument, msg);

Eigen::VectorXd x_bad(static_cast<Eigen::VectorXd::Index>(0));
std::stringstream err_msg2;
err_msg2 << "initial guess has size 0, but "
err_msg2 << "algebra_solver: 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),
@@ -372,7 +200,26 @@ 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, "initial guess[1] is inf, but must be finite!");
std::domain_error,
"algebra_solver: initial guess[1] is inf, but must "
"be finite!");

typedef Eigen::Matrix<T, Eigen::Dynamic, 1> 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<double> 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!");

if (!is_newton) {
EXPECT_THROW_MSG(general_algebra_solver(is_newton, f, x, y, dat, dat_int, 0,
@@ -455,35 +302,3 @@ inline void max_num_steps_test(Eigen::Matrix<T, Eigen::Dynamic, 1>& y,
function_tolerance, max_num_steps),
std::domain_error, msg);
}

template <typename F>
Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1> variadic_eq_test(
const F& f,
const Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic> 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<var, Eigen::Dynamic, 1> 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;
}