1 #include <stan/math/rev.hpp>
2 #include <test/unit/math/rev/util.hpp>
3 #include <test/unit/math/rev/prob/expect_eq_diffs.hpp>
4 #include <test/unit/math/rev/prob/test_gradients.hpp>
5 #include <test/unit/math/prim/prob/agrad_distributions_multi_gp_cholesky.hpp>
6 #include <gtest/gtest.h>
7 #include <string>
8 #include <vector>
9
10 template <typename T_y, typename T_scale, typename T_w>
expect_propto_multi_gp_cholesky_log(T_y y1,T_scale L1,T_w w1,T_y y2,T_scale L2,T_w w2,std::string message="")11 void expect_propto_multi_gp_cholesky_log(T_y y1, T_scale L1, T_w w1, T_y y2,
12 T_scale L2, T_w w2,
13 std::string message = "") {
14 expect_eq_diffs(stan::math::multi_gp_cholesky_log<false>(y1, L1, w1),
15 stan::math::multi_gp_cholesky_log<false>(y2, L2, w2),
16 stan::math::multi_gp_cholesky_log<true>(y1, L1, w1),
17 stan::math::multi_gp_cholesky_log<true>(y2, L2, w2), message);
18 }
19
TEST_F(agrad_distributions_multi_gp_cholesky,Propto)20 TEST_F(agrad_distributions_multi_gp_cholesky, Propto) {
21 using stan::math::to_var;
22 expect_propto_multi_gp_cholesky_log(to_var(y), to_var(L), to_var(w),
23 to_var(y2), to_var(L2), to_var(w2),
24 "All vars: y, w, sigma");
25 }
TEST_F(agrad_distributions_multi_gp_cholesky,ProptoY)26 TEST_F(agrad_distributions_multi_gp_cholesky, ProptoY) {
27 using stan::math::to_var;
28 expect_propto_multi_gp_cholesky_log(to_var(y), L, w, to_var(y2), L, w,
29 "var: y");
30 }
TEST_F(agrad_distributions_multi_gp_cholesky,ProptoYMu)31 TEST_F(agrad_distributions_multi_gp_cholesky, ProptoYMu) {
32 using stan::math::to_var;
33 expect_propto_multi_gp_cholesky_log(to_var(y), L, to_var(w), to_var(y2), L,
34 to_var(w2), "var: y and w");
35 }
TEST_F(agrad_distributions_multi_gp_cholesky,ProptoYSigma)36 TEST_F(agrad_distributions_multi_gp_cholesky, ProptoYSigma) {
37 using stan::math::to_var;
38 expect_propto_multi_gp_cholesky_log(to_var(y), to_var(L), w, to_var(y2),
39 to_var(L2), w, "var: y and sigma");
40 }
TEST_F(agrad_distributions_multi_gp_cholesky,ProptoMu)41 TEST_F(agrad_distributions_multi_gp_cholesky, ProptoMu) {
42 using stan::math::to_var;
43 expect_propto_multi_gp_cholesky_log(y, L, to_var(w), y, L, to_var(w2),
44 "var: w");
45 }
TEST_F(agrad_distributions_multi_gp_cholesky,ProptoMuSigma)46 TEST_F(agrad_distributions_multi_gp_cholesky, ProptoMuSigma) {
47 using stan::math::to_var;
48 expect_propto_multi_gp_cholesky_log(y, to_var(L), to_var(w), y, to_var(L2),
49 to_var(w2), "var: w and sigma");
50 }
TEST_F(agrad_distributions_multi_gp_cholesky,ProptoSigma)51 TEST_F(agrad_distributions_multi_gp_cholesky, ProptoSigma) {
52 using stan::math::to_var;
53 expect_propto_multi_gp_cholesky_log(y, to_var(L), w, y, to_var(L2), w,
54 "var: sigma");
55 }
56
TEST(ProbDistributionsMultiGPCholesky,MultiGPCholeskyVar)57 TEST(ProbDistributionsMultiGPCholesky, MultiGPCholeskyVar) {
58 using Eigen::Dynamic;
59 using Eigen::Matrix;
60 using stan::math::var;
61 Matrix<var, Dynamic, Dynamic> y(3, 3);
62 y << 2.0, -2.0, 11.0, -4.0, 0.0, 2.0, 1.0, 5.0, 3.3;
63 Matrix<var, Dynamic, 1> w(3, 1);
64 w << 1.0, 0.5, 3.0;
65 Matrix<var, Dynamic, Dynamic> Sigma(3, 3);
66 Sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 0.0, 0.0, 0.0, 5.0;
67 Matrix<var, Dynamic, Dynamic> L = Sigma.llt().matrixL();
68 EXPECT_FLOAT_EQ(-46.087162, stan::math::multi_gp_cholesky_log(y, L, w).val());
69 }
70
TEST(ProbDistributionsMultiGPCholesky,MultiGPCholeskyGradientUnivariate)71 TEST(ProbDistributionsMultiGPCholesky, MultiGPCholeskyGradientUnivariate) {
72 using Eigen::Dynamic;
73 using Eigen::Matrix;
74 using Eigen::VectorXd;
75 using stan::math::multi_gp_cholesky_log;
76 using stan::math::var;
77 using std::vector;
78
79 Matrix<var, Dynamic, Dynamic> y_var(1, 1);
80 y_var << 2.0;
81
82 Matrix<var, Dynamic, 1> w_var(1, 1);
83 w_var << 1.0;
84
85 Matrix<var, Dynamic, Dynamic> L_var(1, 1);
86 L_var(0, 0) = 3.0;
87
88 std::vector<var> x;
89 x.push_back(y_var(0));
90 x.push_back(w_var(0));
91 x.push_back(L_var(0, 0));
92
93 var lp = stan::math::multi_gp_cholesky_log(y_var, L_var, w_var);
94 vector<double> grad;
95 lp.grad(x, grad);
96
97 // ===================================
98
99 Matrix<double, Dynamic, Dynamic> y(1, 1);
100 y << 2.0;
101
102 Matrix<double, Dynamic, 1> w(1, 1);
103 w << 1.0;
104
105 Matrix<double, Dynamic, Dynamic> L(1, 1);
106 L << 3.0;
107
108 double epsilon = 1e-6;
109
110 Matrix<double, Dynamic, Dynamic> y_m(1, 1);
111 Matrix<double, Dynamic, Dynamic> y_p(1, 1);
112 y_p(0) = y(0) + epsilon;
113 y_m(0) = y(0) - epsilon;
114 double grad_diff
115 = (multi_gp_cholesky_log(y_p, L, w) - multi_gp_cholesky_log(y_m, L, w))
116 / (2 * epsilon);
117 EXPECT_FLOAT_EQ(grad_diff, grad[0]);
118
119 Matrix<double, Dynamic, 1> w_m(1, 1);
120 Matrix<double, Dynamic, 1> w_p(1, 1);
121 w_p[0] = w[0] + epsilon;
122 w_m[0] = w[0] - epsilon;
123 grad_diff
124 = (multi_gp_cholesky_log(y, L, w_p) - multi_gp_cholesky_log(y, L, w_m))
125 / (2 * epsilon);
126 EXPECT_FLOAT_EQ(grad_diff, grad[1]);
127
128 Matrix<double, Dynamic, Dynamic> L_m(1, 1);
129 Matrix<double, Dynamic, Dynamic> L_p(1, 1);
130 L_p(0) = L(0) + epsilon;
131 L_m(0) = L(0) - epsilon;
132 grad_diff
133 = (multi_gp_cholesky_log(y, L_p, w) - multi_gp_cholesky_log(y, L_m, w))
134 / (2 * epsilon);
135 EXPECT_FLOAT_EQ(grad_diff, grad[2]);
136 }
137
138 struct multi_gp_cholesky_fun {
139 const int K_, N_;
140
multi_gp_cholesky_funmulti_gp_cholesky_fun141 multi_gp_cholesky_fun(int K, int N) : K_(K), N_(N) {}
142
143 template <typename T>
operator ()multi_gp_cholesky_fun144 T operator()(const std::vector<T>& x) const {
145 using Eigen::Dynamic;
146 using Eigen::Matrix;
147 using stan::math::var;
148 Matrix<T, Dynamic, Dynamic> y(K_, N_);
149 Matrix<T, Dynamic, Dynamic> Sigma(N_, N_);
150 Matrix<T, Dynamic, 1> w(K_);
151
152 int pos = 0;
153 for (int j = 0; j < N_; ++j)
154 for (int i = 0; i < K_; ++i)
155 y(i, j) = x[pos++];
156 for (int j = 0; j < N_; ++j) {
157 for (int i = 0; i <= j; ++i) {
158 Sigma(i, j) = x[pos++];
159 Sigma(j, i) = Sigma(i, j);
160 }
161 }
162 Matrix<T, Dynamic, Dynamic> L = Sigma.llt().matrixL();
163 for (int i = 0; i < K_; ++i)
164 w(i) = x[pos++];
165 return stan::math::multi_gp_cholesky_log<false>(y, L, w);
166 }
167 };
168
TEST(MultiGPCholesky,TestGradFunctional)169 TEST(MultiGPCholesky, TestGradFunctional) {
170 std::vector<double> x(3 * 2 + 3 + 3);
171 // y
172 x[0] = 1.0;
173 x[1] = 2.0;
174 x[2] = -3.0;
175
176 x[3] = 0.0;
177 x[4] = -2.0;
178 x[5] = -3.0;
179
180 // Sigma
181 x[6] = 1;
182 x[7] = -1;
183 x[8] = 10;
184
185 // w
186 x[9] = 1;
187 x[10] = 10;
188 x[11] = 5;
189
190 test_grad(multi_gp_cholesky_fun(3, 2), x);
191
192 std::vector<double> u(3);
193 u[0] = 1.9;
194 u[1] = 0.48;
195 u[2] = 2.7;
196
197 test_grad(multi_gp_cholesky_fun(1, 1), u);
198 }
199
TEST(ProbDistributionsMultiGPCholesky,check_varis_on_stack)200 TEST(ProbDistributionsMultiGPCholesky, check_varis_on_stack) {
201 using Eigen::Dynamic;
202 using Eigen::Matrix;
203 using stan::math::to_var;
204 Matrix<double, Dynamic, Dynamic> y(3, 3);
205 y << 2.0, -2.0, 11.0, -4.0, 0.0, 2.0, 1.0, 5.0, 3.3;
206 Matrix<double, Dynamic, 1> w(3, 1);
207 w << 1.0, 0.5, 3.0;
208 Matrix<double, Dynamic, Dynamic> Sigma(3, 3);
209 Sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 0.0, 0.0, 0.0, 5.0;
210 Matrix<double, Dynamic, Dynamic> L = Sigma.llt().matrixL();
211
212 test::check_varis_on_stack(
213 stan::math::multi_gp_cholesky_log<true>(to_var(y), to_var(L), to_var(w)));
214 test::check_varis_on_stack(
215 stan::math::multi_gp_cholesky_log<true>(to_var(y), to_var(L), w));
216 test::check_varis_on_stack(
217 stan::math::multi_gp_cholesky_log<true>(to_var(y), L, to_var(w)));
218 test::check_varis_on_stack(
219 stan::math::multi_gp_cholesky_log<true>(to_var(y), L, w));
220 test::check_varis_on_stack(
221 stan::math::multi_gp_cholesky_log<true>(y, to_var(L), to_var(w)));
222 test::check_varis_on_stack(
223 stan::math::multi_gp_cholesky_log<true>(y, to_var(L), w));
224 test::check_varis_on_stack(
225 stan::math::multi_gp_cholesky_log<true>(y, L, to_var(w)));
226
227 test::check_varis_on_stack(stan::math::multi_gp_cholesky_log<false>(
228 to_var(y), to_var(L), to_var(w)));
229 test::check_varis_on_stack(
230 stan::math::multi_gp_cholesky_log<false>(to_var(y), to_var(L), w));
231 test::check_varis_on_stack(
232 stan::math::multi_gp_cholesky_log<false>(to_var(y), L, to_var(w)));
233 test::check_varis_on_stack(
234 stan::math::multi_gp_cholesky_log<false>(to_var(y), L, w));
235 test::check_varis_on_stack(
236 stan::math::multi_gp_cholesky_log<false>(y, to_var(L), to_var(w)));
237 test::check_varis_on_stack(
238 stan::math::multi_gp_cholesky_log<false>(y, to_var(L), w));
239 test::check_varis_on_stack(
240 stan::math::multi_gp_cholesky_log<false>(y, L, to_var(w)));
241 }
242