1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_ALGEBRA_SOLVER_HPP 2 #define STAN_MATH_REV_MAT_FUNCTOR_ALGEBRA_SOLVER_HPP 13 #include <unsupported/Eigen/NonLinearOptimization> 27 template <
typename Fs,
typename F,
typename T,
typename Fx>
41 const Eigen::Matrix<T, Eigen::Dynamic, 1>& y,
42 const std::vector<double>& dat,
43 const std::vector<int>& dat_int,
44 const Eigen::VectorXd& theta_dbl, Fx& fx,
55 using Eigen::MatrixXd;
56 for (
int i = 0; i < y.size(); ++i)
60 for (
int i = 1; i < x.size(); ++i)
61 theta_[i] =
new vari(theta_dbl(i),
false);
68 f_y(fs, f, theta_dbl,
value_of(y), dat, dat_int, msgs)
73 for (
int j = 0; j <
y_size_; j++)
74 for (
int i = 0; i <
x_size_; i++)
75 y_[j]->
adj_ += theta_[i]->
adj_ * Jx_y_[j * x_size_ + i];
123 template <
typename F,
typename T>
125 const F& f,
const Eigen::Matrix<T, Eigen::Dynamic, 1>& x,
126 const Eigen::VectorXd& y,
const std::vector<double>& dat,
127 const std::vector<int>& dat_int, std::ostream* msgs =
nullptr,
128 double relative_tolerance = 1
e-10,
double function_tolerance = 1
e-6,
129 long int max_num_steps = 1
e+3) {
131 for (
int i = 0; i < x.size(); i++)
133 for (
int i = 0; i < y.size(); i++)
134 check_finite(
"algebra_solver",
"parameter vector", y(i));
137 for (
int x : dat_int)
140 if (relative_tolerance < 0)
142 relative_tolerance,
"",
143 ", must be greater than or equal to 0");
144 if (function_tolerance < 0)
146 function_tolerance,
"",
147 ", must be greater than or equal to 0");
148 if (max_num_steps <= 0)
150 ", must be greater than 0");
155 Fx fx(Fs(), f,
value_of(x), y, dat, dat_int, msgs);
156 Eigen::HybridNonLinearSolver<Fx> solver(fx);
160 fx.get_value(
value_of(x)),
"the vector of unknowns, x,",
164 Eigen::VectorXd theta_dbl =
value_of(x);
165 solver.parameters.xtol = relative_tolerance;
166 solver.parameters.maxfev = max_num_steps;
167 solver.solve(theta_dbl);
170 if (solver.nfev >= max_num_steps) {
171 std::ostringstream message;
172 message <<
"algebra_solver: max number of iterations: " << max_num_steps
174 throw boost::math::evaluation_error(message.str());
178 double system_norm = fx.get_value(theta_dbl).stableNorm();
179 if (system_norm > function_tolerance) {
180 std::ostringstream message2;
181 message2 <<
"algebra_solver: the norm of the algebraic function is: " 182 << system_norm <<
" but should be lower than the function " 183 <<
"tolerance: " << function_tolerance <<
". Consider " 184 <<
"decreasing the relative tolerance and increasing the " 186 throw boost::math::evaluation_error(message2.str());
235 template <
typename F,
typename T1,
typename T2>
237 const F& f,
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& x,
238 const Eigen::Matrix<T2, Eigen::Dynamic, 1>& y,
239 const std::vector<double>& dat,
const std::vector<int>& dat_int,
240 std::ostream* msgs =
nullptr,
double relative_tolerance = 1
e-10,
241 double function_tolerance = 1
e-6,
242 long int max_num_steps = 1
e+3) {
243 Eigen::VectorXd theta_dbl
245 function_tolerance, max_num_steps);
259 dat_int, theta_dbl, fx, msgs);
260 Eigen::Matrix<var, Eigen::Dynamic, 1> theta(x.size());
262 for (
int i = 1; i < x.size(); ++i)
Eigen::Matrix< fvar< T >, R1, C2 > mdivide_left(const Eigen::Matrix< fvar< T >, R1, C1 > &A, const Eigen::Matrix< fvar< T >, R2, C2 > &b)
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
void check_nonzero_size(const char *function, const char *name, const T_y &y)
Check if the specified matrix/vector is of non-zero size.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
A functor that allows us to treat either x or y as the independent variable.
Eigen::VectorXd algebra_solver(const F &f, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &x, const Eigen::VectorXd &y, const std::vector< double > &dat, const std::vector< int > &dat_int, std::ostream *msgs=nullptr, double relative_tolerance=1e-10, double function_tolerance=1e-6, long int max_num_steps=1e+3)
Return the solution to the specified system of algebraic equations given an initial guess...
double * Jx_y_
Jacobian of the solution w.r.t parameters.
The variable implementation base class.
algebra_solver_vari(const Fs &fs, const F &f, const Eigen::VectorXd &x, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &y, const std::vector< double > &dat, const std::vector< int > &dat_int, const Eigen::VectorXd &theta_dbl, Fx &fx, std::ostream *msgs)
int x_size_
number of unknowns
vari ** y_
vector of parameters
A functor with the required operators to call Eigen's algebraic solver.
The vari class for the algebraic solver.
void invalid_argument(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw an invalid_argument exception with a consistently formatted message.
void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
vari ** theta_
vector of solution
double e()
Return the base of the natural logarithm.
vari(double x)
Construct a variable implementation from a value.
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
void check_matching_sizes(const char *function, const char *name1, const T_y1 &y1, const char *name2, const T_y2 &y2)
Check if two structures at the same size.
int y_size_
number of parameters
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
This struct always provides access to the autodiff stack using the singleton pattern.