Automatic Differentiation
 
Loading...
Searching...
No Matches
algebra_solver_fp.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_REV_FUNCTOR_FP_SOLVER_HPP
2#define STAN_MATH_REV_FUNCTOR_FP_SOLVER_HPP
3
14#include <sundials/sundials_context.h>
15#include <kinsol/kinsol.h>
16#include <sunmatrix/sunmatrix_dense.h>
17#include <sunlinsol/sunlinsol_dense.h>
18#include <nvector/nvector_serial.h>
19#include <algorithm>
20#include <iostream>
21#include <string>
22#include <vector>
23
24namespace stan {
25namespace math {
26
34template <typename F>
37 sundials::Context sundials_context_;
39 const F& f_;
42 const Eigen::VectorXd y_dummy;
44 const Eigen::VectorXd& y_;
46 const size_t N_;
48 const size_t M_;
50 const std::vector<double>& x_r_;
52 const std::vector<int>& x_i_;
54 std::ostream* msgs_;
56 void* mem_;
58 N_Vector nv_x_;
60 N_Vector nv_u_scal_;
62 N_Vector nv_f_scal_;
63
65 template <typename T, typename T_u, typename T_f>
66 KinsolFixedPointEnv(const F& f, const Eigen::Matrix<T, -1, 1>& x,
67 const Eigen::VectorXd& y, const std::vector<double>& x_r,
68 const std::vector<int>& x_i, std::ostream* msgs,
69 const std::vector<T_u>& u_scale,
70 const std::vector<T_f>& f_scale)
72 f_(f),
73 y_dummy(),
74 y_(y),
75 N_(x.size()),
76 M_(y.size()),
77 x_r_(x_r),
78 x_i_(x_i),
79 msgs_(msgs),
80 mem_(KINCreate(sundials_context_)),
81 nv_x_(N_VNew_Serial(N_, sundials_context_)),
82 nv_u_scal_(N_VNew_Serial(N_, sundials_context_)),
83 nv_f_scal_(N_VNew_Serial(N_, sundials_context_)) {
84 for (int i = 0; i < N_; ++i) {
85 NV_Ith_S(nv_x_, i) = stan::math::value_of(x(i));
86 NV_Ith_S(nv_u_scal_, i) = stan::math::value_of(u_scale[i]);
87 NV_Ith_S(nv_f_scal_, i) = stan::math::value_of(f_scale[i]);
88 }
89 }
90
92 template <typename T, typename T_u, typename T_f>
93 KinsolFixedPointEnv(const F& f, const Eigen::Matrix<T, -1, 1>& x,
94 const Eigen::Matrix<stan::math::var, -1, 1>& y,
95 const std::vector<double>& x_r,
96 const std::vector<int>& x_i, std::ostream* msgs,
97 const std::vector<T_u>& u_scale,
98 const std::vector<T_f>& f_scale)
100 f_(f),
101 y_dummy(stan::math::value_of(y)),
102 y_(y_dummy),
103 N_(x.size()),
104 M_(y.size()),
105 x_r_(x_r),
106 x_i_(x_i),
107 msgs_(msgs),
108 mem_(KINCreate(sundials_context_)),
109 nv_x_(N_VNew_Serial(N_, sundials_context_)),
110 nv_u_scal_(N_VNew_Serial(N_, sundials_context_)),
111 nv_f_scal_(N_VNew_Serial(N_, sundials_context_)) {
112 for (int i = 0; i < N_; ++i) {
113 NV_Ith_S(nv_x_, i) = stan::math::value_of(x(i));
114 NV_Ith_S(nv_u_scal_, i) = stan::math::value_of(u_scale[i]);
115 NV_Ith_S(nv_f_scal_, i) = stan::math::value_of(f_scale[i]);
116 }
117 }
118
120 N_VDestroy_Serial(nv_x_);
121 N_VDestroy_Serial(nv_u_scal_);
122 N_VDestroy_Serial(nv_f_scal_);
123 KINFree(&mem_);
124 }
125
127 static int kinsol_f_system(N_Vector x, N_Vector f, void* user_data) {
128 auto g = static_cast<const KinsolFixedPointEnv<F>*>(user_data);
129 Eigen::VectorXd x_eigen(Eigen::Map<Eigen::VectorXd>(NV_DATA_S(x), g->N_));
130 Eigen::Map<Eigen::VectorXd>(N_VGetArrayPointer(f), g->N_)
131 = g->f_(x_eigen, g->y_, g->x_r_, g->x_i_, g->msgs_);
132 return 0;
133 }
134};
135
162 template <typename F>
163 inline Eigen::Matrix<stan::math::var, -1, 1> operator()(
164 const Eigen::VectorXd& x, const Eigen::Matrix<stan::math::var, -1, 1>& y,
168 using stan::math::var;
169
170 auto g = [&env](const Eigen::Matrix<var, -1, 1>& par_) {
171 Eigen::Matrix<var, -1, 1> x_(par_.head(env.N_));
172 Eigen::Matrix<var, -1, 1> y_(par_.tail(env.M_));
173 return env.f_(x_, y_, env.x_r_, env.x_i_, env.msgs_);
174 };
175
176 Eigen::VectorXd theta(x.size() + y.size());
177 for (int i = 0; i < env.N_; ++i) {
178 theta(i) = x(i);
179 }
180 for (int i = 0; i < env.M_; ++i) {
181 theta(i + env.N_) = env.y_(i);
182 }
183 Eigen::Matrix<double, -1, 1> fx;
184 Eigen::Matrix<double, -1, -1> J_theta;
185 stan::math::jacobian(g, theta, fx, J_theta);
186 Eigen::MatrixXd A(J_theta.block(0, 0, env.N_, env.N_));
187 Eigen::MatrixXd b(J_theta.block(0, env.N_, env.N_, env.M_));
188 A = Eigen::MatrixXd::Identity(env.N_, env.N_) - A;
189 Eigen::MatrixXd Jxy = A.colPivHouseholderQr().solve(b);
190 std::vector<double> gradients(env.M_);
191 Eigen::Matrix<var, -1, 1> x_sol(env.N_);
192 std::vector<stan::math::var> yv(to_array_1d(y));
193 for (int i = 0; i < env.N_; ++i) {
194 gradients = to_array_1d(Eigen::VectorXd(Jxy.row(i)));
195 x_sol[i] = precomputed_gradients(x(i), yv, gradients);
196 }
197 return x_sol;
198 }
199};
200
220template <typename fp_env_type, typename fp_jac_type>
222
229template <typename F, typename fp_jac_type>
230struct FixedPointSolver<KinsolFixedPointEnv<F>, fp_jac_type> {
239 void kinsol_solve_fp(Eigen::VectorXd& x, KinsolFixedPointEnv<F>& env,
240 double f_tol, int max_num_steps) {
241 int N = env.N_;
242 void* mem = env.mem_;
243
244 const int default_anderson_depth = 4;
245 int anderson_depth = std::min(N, default_anderson_depth);
246
247 CHECK_KINSOL_CALL(KINSetNumMaxIters(mem, max_num_steps));
248 CHECK_KINSOL_CALL(KINSetMAA(mem, anderson_depth));
249 CHECK_KINSOL_CALL(KINInit(mem, &env.kinsol_f_system, env.nv_x_));
250 CHECK_KINSOL_CALL(KINSetFuncNormTol(mem, f_tol));
251 CHECK_KINSOL_CALL(KINSetUserData(mem, static_cast<void*>(&env)));
252 kinsol_check(KINSol(mem, env.nv_x_, KIN_FP, env.nv_u_scal_, env.nv_f_scal_),
253 "KINSol", max_num_steps);
254
255 for (int i = 0; i < N; ++i) {
256 x(i) = NV_Ith_S(env.nv_x_, i);
257 }
258 }
259
271 template <typename T1>
272 Eigen::Matrix<double, -1, 1> solve(const Eigen::Matrix<T1, -1, 1>& x,
273 const Eigen::Matrix<double, -1, 1>& y,
274 KinsolFixedPointEnv<F>& env, double f_tol,
275 int max_num_steps) {
276 Eigen::VectorXd xd(stan::math::value_of(x));
277 kinsol_solve_fp(xd, env, f_tol, max_num_steps);
278 return xd;
279 }
280
292 template <typename T1>
293 Eigen::Matrix<stan::math::var, -1, 1> solve(
294 const Eigen::Matrix<T1, -1, 1>& x,
295 const Eigen::Matrix<stan::math::var, -1, 1>& y,
296 KinsolFixedPointEnv<F>& env, double f_tol, int max_num_steps) {
298 using stan::math::var;
299
300 // FP solution
301 Eigen::VectorXd xd(solve(x, Eigen::VectorXd(), env, f_tol, max_num_steps));
302
303 fp_jac_type jac_sol;
304 return jac_sol(xd, y, env);
305 }
306};
307
364template <typename F, typename T1, typename T2, typename T_u, typename T_f>
365Eigen::Matrix<T2, -1, 1> algebra_solver_fp(
366 const F& f, const Eigen::Matrix<T1, -1, 1>& x,
367 const Eigen::Matrix<T2, -1, 1>& y, const std::vector<double>& x_r,
368 const std::vector<int>& x_i, const std::vector<T_u>& u_scale,
369 const std::vector<T_f>& f_scale, std::ostream* msgs = nullptr,
370 double f_tol = 1e-8,
371 int max_num_steps = 200) { // NOLINT(runtime/int)
372 algebra_solver_check(x, y, x_r, x_i, f_tol, max_num_steps);
373 check_nonnegative("algebra_solver", "u_scale", u_scale);
374 check_nonnegative("algebra_solver", "f_scale", f_scale);
375 check_matching_sizes("algebra_solver", "the algebraic system's output",
376 value_of(f(x, y, x_r, x_i, msgs)),
377 "the vector of unknowns, x,", x);
378
379 KinsolFixedPointEnv<F> env(f, x, y, x_r, x_i, msgs, u_scale,
380 f_scale); // NOLINT
382 return fp.solve(x, y, env, f_tol, max_num_steps);
383}
384
385} // namespace math
386} // namespace stan
387
388#endif
#define CHECK_KINSOL_CALL(call)
auto to_array_1d(T_x &&x)
Returns input matrix reshaped into a vector.
int64_t size(const T &m)
Returns the size (number of the elements) of a matrix_cl or var_value<matrix_cl<T>>.
Definition size.hpp:19
void jacobian(const F &f, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &fx, Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &J)
Definition jacobian.hpp:11
void check_nonnegative(const char *function, const char *name, const T_y &y)
Check if y is non-negative.
void algebra_solver_check(const Eigen::Matrix< T1, Eigen::Dynamic, 1 > &x, const Eigen::Matrix< T2, Eigen::Dynamic, 1 > y, const std::vector< double > &dat, const std::vector< int > &dat_int, double function_tolerance, long int max_num_steps)
static constexpr double e()
Return the base of the natural logarithm.
Definition constants.hpp:20
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
Eigen::Matrix< T2, -1, 1 > algebra_solver_fp(const F &f, const Eigen::Matrix< T1, -1, 1 > &x, const Eigen::Matrix< T2, -1, 1 > &y, const std::vector< double > &x_r, const std::vector< int > &x_i, const std::vector< T_u > &u_scale, const std::vector< T_f > &f_scale, std::ostream *msgs=nullptr, double f_tol=1e-8, int max_num_steps=200)
Return a fixed pointer to the specified system of algebraic equations of form.
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.
var_value< double > var
Definition var.hpp:1187
void kinsol_check(int flag, const char *func_name)
Throws an exception message when the functions in KINSOL fails.
var precomputed_gradients(Arith value, const VecVar &operands, const VecArith &gradients, const std::tuple< ContainerOperands... > &container_operands=std::tuple<>(), const std::tuple< ContainerGradients... > &container_gradients=std::tuple<>())
This function returns a var for an expression that has the specified value, vector of operands,...
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
Eigen::Matrix< stan::math::var, -1, 1 > operator()(const Eigen::VectorXd &x, const Eigen::Matrix< stan::math::var, -1, 1 > &y, KinsolFixedPointEnv< F > &env)
Calculate Jacobian Jxy.
Calculate Jacobian Jxy(Jacobian of unknown x w.r.t.
Eigen::Matrix< stan::math::var, -1, 1 > solve(const Eigen::Matrix< T1, -1, 1 > &x, const Eigen::Matrix< stan::math::var, -1, 1 > &y, KinsolFixedPointEnv< F > &env, double f_tol, int max_num_steps)
Solve FP problem and calculate jacobian.
Eigen::Matrix< double, -1, 1 > solve(const Eigen::Matrix< T1, -1, 1 > &x, const Eigen::Matrix< double, -1, 1 > &y, KinsolFixedPointEnv< F > &env, double f_tol, int max_num_steps)
Solve data-only FP problem so no need to calculate jacobian.
void kinsol_solve_fp(Eigen::VectorXd &x, KinsolFixedPointEnv< F > &env, double f_tol, int max_num_steps)
Solve FP using KINSOL.
Fixed point solver for problem of form.
const std::vector< int > & x_i_
integer data
void * mem_
KINSOL memory block.
const std::vector< double > & x_r_
real data
std::ostream * msgs_
message stream
N_Vector nv_u_scal_
NVECTOR for scaling u.
N_Vector nv_x_
NVECTOR for unknowns.
KinsolFixedPointEnv(const F &f, const Eigen::Matrix< T, -1, 1 > &x, const Eigen::VectorXd &y, const std::vector< double > &x_r, const std::vector< int > &x_i, std::ostream *msgs, const std::vector< T_u > &u_scale, const std::vector< T_f > &f_scale)
Constructor when y is data.
static int kinsol_f_system(N_Vector x, N_Vector f, void *user_data)
Implements the user-defined function passed to KINSOL.
sundials::Context sundials_context_
Sundials context.
N_Vector nv_f_scal_
NVECTOR for scaling f.
const Eigen::VectorXd & y_
ref to val of params
const Eigen::VectorXd y_dummy
val of params for y_ to refer to when params are var type
KinsolFixedPointEnv(const F &f, const Eigen::Matrix< T, -1, 1 > &x, const Eigen::Matrix< stan::math::var, -1, 1 > &y, const std::vector< double > &x_r, const std::vector< int > &x_i, std::ostream *msgs, const std::vector< T_u > &u_scale, const std::vector< T_f > &f_scale)
Constructor when y is param.
KINSOL algebraic system data holder that handles construction & destruction of SUNDIALS data,...