1 #ifndef STAN_MATH_TEST_FIXTURE_ODE_LORENZ_HPP 2 #define STAN_MATH_TEST_FIXTURE_ODE_LORENZ_HPP 3 4 #include <stan/math/rev.hpp> 5 #include <boost/numeric/odeint.hpp> 6 #include <gtest/gtest.h> 7 #include <test/unit/math/prim/functor/lorenz.hpp> 8 #include <test/unit/math/rev/functor/test_fixture_ode.hpp> 9 #include <test/unit/util.hpp> 10 #include <iostream> 11 #include <sstream> 12 #include <vector> 13 #include <limits> 14 #include <string> 15 16 struct lorenz_ode_base { 17 struct lorenz_rhs { 18 template <typename T0, typename T1, typename T2> operator ()lorenz_ode_base::lorenz_rhs19 inline Eigen::Matrix<stan::return_type_t<T1, T2>, -1, 1> operator()( 20 const T0& t_in, const T1& y_in, std::ostream* msgs, 21 const T2& theta) const { 22 Eigen::Matrix<stan::return_type_t<T1, T2>, -1, 1> res(3); 23 res << theta.at(0) * (y_in(1) - y_in(0)), 24 theta.at(1) * y_in(0) - y_in(1) - y_in(0) * y_in(2), 25 -theta.at(2) * y_in(2) + y_in(0) * y_in(1); 26 return res; 27 } 28 }; 29 30 lorenz_rhs f; 31 Eigen::VectorXd y0; 32 std::vector<double> theta; 33 double t0; 34 std::vector<double> ts; 35 double rtol; 36 double atol; 37 int max_num_step; 38 lorenz_ode_baselorenz_ode_base39 lorenz_ode_base() 40 : f(), 41 y0(3), 42 theta{10.0, 28.0, 8.0 / 3.0}, 43 t0(0.0), 44 ts(100), 45 rtol(1.e-8), 46 atol(1.e-10), 47 max_num_step(100000) { 48 y0 << 10.0, 1.0, 1.0; 49 for (int i = 0; i < 100; i++) 50 ts[i] = t0 + (0.1 * (i + 1)); 51 } 52 timeslorenz_ode_base53 std::vector<double>& times() { return ts; } initlorenz_ode_base54 Eigen::VectorXd init() { return y0; } paramlorenz_ode_base55 std::vector<double> param() { return theta; } 56 }; 57 58 /** 59 * Inheriting base type, various fixtures differs by the type of ODE 60 * functor used in <code>apply_solver</code> calls, intended for 61 * different kind of tests. 62 * 63 */ 64 template <typename T> 65 struct lorenz_test : public lorenz_ode_base, 66 public ODETestFixture<lorenz_test<T>> { lorenz_testlorenz_test67 lorenz_test() : lorenz_ode_base() {} 68 apply_solverlorenz_test69 auto apply_solver() { 70 std::tuple_element_t<0, T> sol; 71 return sol(this->f, this->y0, this->t0, this->ts, nullptr, this->theta); 72 } 73 74 template <typename T1, typename T2> apply_solverlorenz_test75 auto apply_solver(T1&& init, T2&& theta_in) { 76 std::tuple_element_t<0, T> sol; 77 return sol(this->f, init, this->t0, this->ts, nullptr, theta_in); 78 } 79 apply_solver_tollorenz_test80 auto apply_solver_tol() { 81 std::tuple_element_t<1, T> sol; 82 return sol(this->f, this->y0, this->t0, this->ts, this->rtol, this->atol, 83 this->max_num_step, nullptr, this->theta); 84 } 85 }; 86 87 #endif 88