From 381ad687a68cbb661df2610fe2e176ea2626576e Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Mon, 22 Nov 2021 15:53:44 -0800 Subject: [PATCH 01/28] revamp matrix_exp_multiply & more unit test matrices --- .../prim/fun/matrix_exp_action_handler.hpp | 208 +++++++++++++----- .../prim/fun/matrix_exp_multiply_test.cpp | 82 +++++-- 2 files changed, 221 insertions(+), 69 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index 9c7bcfc1230..49a8d0b588f 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -3,6 +3,8 @@ #include #include +#include +#include #include #include @@ -14,6 +16,7 @@ namespace math { * "Computing the Action of the Matrix Exponential, * with an Application to Exponential Integrators" * Read More: https://epubs.siam.org/doi/abs/10.1137/100788860 + * See also: https://www.mathworks.com/matlabcentral/fileexchange/29576-matrix-exponential-times-a-vector * * Calculates exp(mat*t)*b, where mat & b are matrices, * and t is double. @@ -21,19 +24,51 @@ namespace math { class matrix_exp_action_handler { static constexpr int p_max = 8; static constexpr int m_max = 55; - static constexpr double tol = 1.1e-16; - - // table 3.1 in the reference - const std::vector theta_m_single_precision{ - 1.3e-1, 1.0e0, 2.2e0, 3.6e0, 4.9e0, 6.3e0, - 7.7e0, 9.1e0, 1.1e1, 1.2e1, 1.3e1}; - const std::vector theta_m_double_precision{ - 2.4e-3, 1.4e-1, 6.4e-1, 1.4e0, 2.4e0, 3.5e0, - 4.7e0, 6.0e0, 7.2e0, 8.5e0, 9.9e0}; - - double l1norm(const Eigen::MatrixXd& m) { - return m.colwise().lpNorm<1>().maxCoeff(); - } + static constexpr double tol = 1.0e-16; + + // + const std::vector theta_m_single{ + 1.19209280e-07, 5.97885889e-04, 1.12338647e-02, 5.11661936e-02, + 1.30848716e-01, 2.49528932e-01, 4.01458242e-01, 5.80052463e-01, + 7.79511337e-01, 9.95184079e-01, 1.22347954e+00, 1.46166151e+00, + 1.70764853e+00, 1.95985059e+00, 2.21704439e+00, 2.47828088e+00, + 2.74281711e+00, 3.01006636e+00, 3.27956121e+00, 3.55092621e+00, + 3.82385743e+00, 4.09810697e+00, 4.37347131e+00, 4.64978222e+00, + 4.92689984e+00, 5.20470723e+00, 5.48310609e+00, 5.76201341e+00, + 6.04135876e+00, 6.32108213e+00, 6.60113218e+00, 6.88146485e+00, + 7.16204215e+00, 7.44283129e+00, 7.72380381e+00, 8.00493499e+00, + 8.28620327e+00, 8.56758983e+00, 8.84907819e+00, 9.13065388e+00, + 9.41230420e+00, 9.69401796e+00, 9.97578527e+00, 1.02575974e+01, + 1.05394467e+01, 1.08213263e+01, 1.11032302e+01, 1.13851530e+01, + 1.16670899e+01, 1.19490366e+01, 1.22309896e+01, 1.25129453e+01, + 1.27949008e+01, 1.30768536e+01, 1.33588011e+01, 1.36407415e+01, + 1.39226727e+01, 1.42045932e+01, 1.44865015e+01, 1.47683963e+01}; + const std::vector theta_m_double{ + 2.22044605e-16, 2.58095680e-08, 1.38634787e-05, 3.39716884e-04, + 2.40087636e-03, 9.06565641e-03, 2.38445553e-02, 4.99122887e-02, + 8.95776020e-02, 1.44182976e-01, 2.14235807e-01, 2.99615891e-01, + 3.99777534e-01, 5.13914694e-01, 6.41083523e-01, 7.80287426e-01, + 9.30532846e-01, 1.09086372e+00, 1.26038106e+00, 1.43825260e+00, + 1.62371595e+00, 1.81607782e+00, 2.01471078e+00, 2.21904887e+00, + 2.42858252e+00, 2.64285346e+00, 2.86144963e+00, 3.08400054e+00, + 3.31017284e+00, 3.53966635e+00, 3.77221050e+00, 4.00756109e+00, + 4.24549744e+00, 4.48581986e+00, 4.72834735e+00, 4.97291563e+00, + 5.21937537e+00, 5.46759063e+00, 5.71743745e+00, 5.96880263e+00, + 6.22158266e+00, 6.47568274e+00, 6.73101590e+00, 6.98750228e+00, + 7.24506843e+00, 7.50364669e+00, 7.76317466e+00, 8.02359473e+00, + 8.28485363e+00, 8.54690205e+00, 8.80969427e+00, 9.07318789e+00, + 9.33734351e+00, 9.60212447e+00, 9.86749668e+00, 1.01334283e+01, + 1.03998897e+01, 1.06668532e+01, 1.09342929e+01, 1.12021845e+01, + 1.14705053e+01, 1.17392341e+01, 1.20083509e+01, 1.22778370e+01, + 1.25476748e+01, 1.28178476e+01, 1.30883399e+01, 1.33591369e+01, + 1.36302250e+01, 1.39015909e+01, 1.41732223e+01, 1.44451076e+01, + 1.47172357e+01, 1.49895963e+01, 1.52621795e+01, 1.55349758e+01, + 1.58079765e+01, 1.60811732e+01, 1.63545578e+01, 1.66281227e+01, + 1.69018609e+01, 1.71757655e+01, 1.74498298e+01, 1.77240478e+01, + 1.79984136e+01, 1.82729215e+01, 1.85475662e+01, 1.88223426e+01, + 1.90972458e+01, 1.93722711e+01, 1.96474142e+01, 1.99226707e+01, + 2.01980367e+01, 2.04735082e+01, 2.07490816e+01, 2.10247533e+01, + 2.13005199e+01, 2.15763782e+01, 2.18523250e+01, 2.21283574e+01}; public: /** @@ -54,39 +89,65 @@ class matrix_exp_action_handler { A(i, i) -= mu; } - int m{0}, s{0}; - set_approximation_parameter(A, t, m, s); - - Eigen::MatrixXd res(A.rows(), b.cols()); - - for (int col = 0; col < b.cols(); ++col) { - bool conv = false; - Eigen::VectorXd B = b.col(col); - Eigen::VectorXd F = B; - const auto eta = std::exp(t * mu / s); - for (int i = 1; i < s + 1; ++i) { - auto c1 = B.template lpNorm(); - if (m > 0) { - for (int j = 1; j < m + 1; ++j) { - B = t * A * B / (s * j); - auto c2 = B.template lpNorm(); - F += B; - if (c1 + c2 < tol * F.template lpNorm()) { - conv = true; - break; - } - c1 = c2; - } - } - F *= eta; - B = F; - if (conv) { - break; - } + int m, s; + set_approx_order(A, b, t, m, s); + + double eta = exp(t * mu/s); + + Eigen::MatrixXd f = b; + Eigen::MatrixXd bi = b; + +#define MAT_OP_INF_NORM(X) X.cwiseAbs().rowwise().sum().maxCoeff() + + for (auto i = 0; i < s; ++i) { + // maximum absolute row sum, aka inf-norm of matrix operator + double c1 = MAT_OP_INF_NORM(bi); + for (auto j = 0; j < m; ++j) { + bi = (t/(s * (j + 1))) * (A * bi); + f += bi; + double c2 = MAT_OP_INF_NORM(bi); + if (c1 + c2 < tol * MAT_OP_INF_NORM(f)) break; + c1 = c2; + } + f *= eta; + bi = f; + } + +#undef MAT_OP_INF_NORM + + return f; + } + + /** + * Estimate the 1-norm of mat^m + * + * See A. H. Al-Mohy and N. J. Higham, A New Scaling and Squaring + * Algorithm for the Matrix Exponential, SIAM J. Matrix Anal. Appl. 31(3): + * 970-989, 2009. + * + * For positive matrices the results is exact. Otherwise it falls + * back to Eigen's norm, which is only + * efficient for small & medium-size matrices (n < 100). Large size + * matrices require a more efficient 1-norm approximation + * algorithm such as normest1. See, e.g., + * https://hg.savannah.gnu.org/hgweb/octave/file/e35866e6a2e0/scripts/linear-algebra/normest1.m + * + * @param mat matrix + * @param m power + */ + template * = nullptr, + require_all_st_same* = nullptr> + double mat_power_1_norm(const EigMat1& mat, int m) { + if ((mat.array() > 0.0).all()) { + Eigen::VectorXd e = Eigen::VectorXd::Constant(mat.rows(), 1.0); + for (auto j = 0; j < m; ++j) { + e = mat.transpose() * e; } - res.col(col) = F; - } // loop b columns - return res; + return e.lpNorm(); + } else { + return mat.pow(m).cwiseAbs().colwise().sum().maxCoeff(); + } } /** @@ -104,22 +165,57 @@ class matrix_exp_action_handler { * @param [out] m int parameter m * @param [out] s int parameter s */ - inline void set_approximation_parameter(const Eigen::MatrixXd& mat, - const double& t, int& m, int& s) { - if (l1norm(mat) < tol || t < tol) { + template * = nullptr, + require_all_st_same* = nullptr> + inline void set_approx_order(const EigMat1& mat, const EigMat2& b, + const double& t, int& m, int& s) { + // L1 norm + double normA = mat.colwise(). template lpNorm<1>().maxCoeff(); + + if (normA < tol || t < tol) { m = 0; s = 1; } else { - m = m_max; + const std::vector& theta = theta_m_double; + Eigen::VectorXd alpha(p_max - 1); - Eigen::MatrixXd a = mat * t; - const double theta_m = theta_m_double_precision.back(); - for (auto i = 0; i < std::ceil(std::log2(p_max)); ++i) { - a *= a; + if (normA < 4.0*theta[m_max]*p_max*(p_max + 3)/(m_max*b.cols())) { + alpha = Eigen::VectorXd::Constant(p_max-1, normA); + } else { + Eigen::VectorXd eta(p_max); + for (auto p = 0; p < p_max; ++p) { + double c = mat_power_1_norm(mat, p+2); + c = std::pow(c, 1.0/(p+2.0)); + eta[p] = c; + } + for (auto p = 0; p < p_max - 1; ++p) { + alpha[p] = std::max(eta[p],eta(p+1)); + } + } + + Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(p_max-1, m_max); + for (auto p = 1; p < p_max; ++p) { + for (auto i = p*(p+1)-2; i < m_max; ++i) { + mt(p-1, i) = alpha[p-1]/theta[i]; + } + } + + Eigen::Matrix u(m_max); + for (auto i = 0; i < m_max; ++i) { + u(i) = i + 1; + } + + Eigen::MatrixXd c = stan::math::ceil(mt) * u.asDiagonal(); + for (auto i = 0; i < c.size(); ++i) { + if (c(i) <= 1.e-16) { + c(i) = std::numeric_limits::infinity(); + } } - double ap = std::pow(l1norm(a), 1.0 / p_max); - int c = std::ceil(ap / theta_m); - s = (c < 1 ? 1 : c); + + int cost = c.colwise().minCoeff().minCoeff(&m); + if (std::isinf(cost)) cost = 0; + s = std::max(cost/m, 1); } } }; diff --git a/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp b/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp index 27772f4e348..558980371a5 100644 --- a/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp +++ b/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -21,23 +22,37 @@ inline void test_matrix_exp_multiply() { // matrix_exp_multiply Eigen::Matrix res = matrix_exp_multiply(A, B); - EXPECT_EQ(res.size(), expAB.size()); - for (int l = 0; l < res.size(); ++l) { - EXPECT_FLOAT_EQ(res(l), expAB(l)); - } + + + EXPECT_MATRIX_FLOAT_EQ(expAB, res); } -TEST(MathMatrixPrimMat, matrix_exp_multiply) { - // the helper above doesn't handle 0 size inputs - Eigen::MatrixXd A(0, 0); - Eigen::MatrixXd B(0, 0); - EXPECT_EQ(stan::math::matrix_exp_multiply(A, B).size(), 0); +TEST(MathMatrixPrimMat, matrix_exp_multiply_power_1_norm_fun) { + stan::math::matrix_exp_action_handler maexp; + Eigen::MatrixXd x(2,2); + x << 1, -2, -3, 4; + EXPECT_FLOAT_EQ(maexp.mat_power_1_norm(x, 2), 32.0); + EXPECT_FLOAT_EQ(maexp.mat_power_1_norm(x, 3), 172.0); - Eigen::MatrixXd C(0, 2); - Eigen::MatrixXd M = stan::math::matrix_exp_multiply(A, C); - EXPECT_EQ(A.rows(), M.rows()); - EXPECT_EQ(C.cols(), M.cols()); + Eigen::MatrixXd y(3,3); + y << 1, 2, 3, 4, 5, 6, 7, 1, 2; + EXPECT_FLOAT_EQ(maexp.mat_power_1_norm(y, 3), 1163.0); + EXPECT_FLOAT_EQ(maexp.mat_power_1_norm(y, 12), 8.3805595e+11); +} + +TEST(MathMatrixPrimMat, matrix_exp_multiply_approx_order) { + stan::math::matrix_exp_action_handler maexp; + Eigen::MatrixXd x(2,2); + x << 1, -2, -3, 4; + Eigen::VectorXd b(2); + b << 1.5, 2.5; + int m, s; + maexp.set_approx_order(x, b, 1.0, m, s); + EXPECT_EQ(m, 40); + EXPECT_EQ(s, 1); +} +TEST(MathMatrixPrimMat, matrix_exp_multiply) { test_matrix_exp_multiply<1, 1>(); test_matrix_exp_multiply<1, 5>(); test_matrix_exp_multiply<5, 1>(); @@ -45,6 +60,47 @@ TEST(MathMatrixPrimMat, matrix_exp_multiply) { test_matrix_exp_multiply<20, 2>(); } +TEST(MathMatrixPrimMat, matrix_exp_multiply_issue_2529) { + // issue #2529 https://github.com/stan-dev/math/issues/2529 + Eigen::MatrixXd a(3, 3); + a << -1.2, 1.2, 0, + 0, -0.5, 0.5, + 0, 0, -0.3; + Eigen::VectorXd b(3); + b << 1.0, 1.0, 1.0; + for (auto i = 15; i < 40; ++i) { + double t = double(i); + Eigen::MatrixXd m1 = stan::math::matrix_exp_multiply(a * t, b); + Eigen::MatrixXd m2 = stan::math::matrix_exp(a * t) * b; + EXPECT_MATRIX_FLOAT_EQ(m1, m2); + } +} + +TEST(MathMatrixPrimMat, matrix_exp_multiply_poisson_5) { +// Block tridiag mat by discretizing Poisson's Eq. on a 5x5 grid + int n = 5; + int nd = n * n; + Eigen::MatrixXd m = Eigen::MatrixXd::Zero(nd, nd); + Eigen::VectorXd b(nd); + for (auto i = 0; i < nd; ++i) { + b(i) = -1.0 + i * 2.0/(nd - 1); + m(i, i) = 4.0; + if (i + 1 < nd) { + m(i, i+1) = -1.0; + m(i+1, i) = -1.0; + } + if (i + 5 < nd) { + m(i, std::min(i + 5, nd - 1)) = -1.0; + m(std::min(i + 5, nd - 1), i) = -1.0; + } + } + + double t = 1.0; + Eigen::MatrixXd p1 = stan::math::matrix_exp_multiply(m * t, b); + Eigen::MatrixXd p2 = stan::math::matrix_exp(m * t) * b; + EXPECT_MATRIX_NEAR(p1, p2, 1.e-12); +} + TEST(MathMatrixPrimMat, matrix_exp_multiply_exception) { using stan::math::matrix_exp_multiply; { // multiplicable From 25a7d8be94de7a5d07c95a0bd3c0118bb724806a Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Tue, 23 Nov 2021 08:41:18 -0800 Subject: [PATCH 02/28] fix header --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index 49a8d0b588f..af21d92b6ac 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -3,7 +3,7 @@ #include #include -#include +#include #include #include #include From 97fc786d1841d3b710099b4a22a79479b3fff93d Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 23 Nov 2021 16:43:13 +0000 Subject: [PATCH 03/28] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- .../prim/fun/matrix_exp_action_handler.hpp | 127 +++++++++--------- .../prim/fun/matrix_exp_multiply_test.cpp | 19 ++- 2 files changed, 73 insertions(+), 73 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index af21d92b6ac..e8f26664363 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -16,7 +16,8 @@ namespace math { * "Computing the Action of the Matrix Exponential, * with an Application to Exponential Integrators" * Read More: https://epubs.siam.org/doi/abs/10.1137/100788860 - * See also: https://www.mathworks.com/matlabcentral/fileexchange/29576-matrix-exponential-times-a-vector + * See also: + * https://www.mathworks.com/matlabcentral/fileexchange/29576-matrix-exponential-times-a-vector * * Calculates exp(mat*t)*b, where mat & b are matrices, * and t is double. @@ -26,49 +27,49 @@ class matrix_exp_action_handler { static constexpr int m_max = 55; static constexpr double tol = 1.0e-16; - // + // const std::vector theta_m_single{ - 1.19209280e-07, 5.97885889e-04, 1.12338647e-02, 5.11661936e-02, - 1.30848716e-01, 2.49528932e-01, 4.01458242e-01, 5.80052463e-01, - 7.79511337e-01, 9.95184079e-01, 1.22347954e+00, 1.46166151e+00, - 1.70764853e+00, 1.95985059e+00, 2.21704439e+00, 2.47828088e+00, - 2.74281711e+00, 3.01006636e+00, 3.27956121e+00, 3.55092621e+00, - 3.82385743e+00, 4.09810697e+00, 4.37347131e+00, 4.64978222e+00, - 4.92689984e+00, 5.20470723e+00, 5.48310609e+00, 5.76201341e+00, - 6.04135876e+00, 6.32108213e+00, 6.60113218e+00, 6.88146485e+00, - 7.16204215e+00, 7.44283129e+00, 7.72380381e+00, 8.00493499e+00, - 8.28620327e+00, 8.56758983e+00, 8.84907819e+00, 9.13065388e+00, - 9.41230420e+00, 9.69401796e+00, 9.97578527e+00, 1.02575974e+01, - 1.05394467e+01, 1.08213263e+01, 1.11032302e+01, 1.13851530e+01, - 1.16670899e+01, 1.19490366e+01, 1.22309896e+01, 1.25129453e+01, - 1.27949008e+01, 1.30768536e+01, 1.33588011e+01, 1.36407415e+01, - 1.39226727e+01, 1.42045932e+01, 1.44865015e+01, 1.47683963e+01}; + 1.19209280e-07, 5.97885889e-04, 1.12338647e-02, 5.11661936e-02, + 1.30848716e-01, 2.49528932e-01, 4.01458242e-01, 5.80052463e-01, + 7.79511337e-01, 9.95184079e-01, 1.22347954e+00, 1.46166151e+00, + 1.70764853e+00, 1.95985059e+00, 2.21704439e+00, 2.47828088e+00, + 2.74281711e+00, 3.01006636e+00, 3.27956121e+00, 3.55092621e+00, + 3.82385743e+00, 4.09810697e+00, 4.37347131e+00, 4.64978222e+00, + 4.92689984e+00, 5.20470723e+00, 5.48310609e+00, 5.76201341e+00, + 6.04135876e+00, 6.32108213e+00, 6.60113218e+00, 6.88146485e+00, + 7.16204215e+00, 7.44283129e+00, 7.72380381e+00, 8.00493499e+00, + 8.28620327e+00, 8.56758983e+00, 8.84907819e+00, 9.13065388e+00, + 9.41230420e+00, 9.69401796e+00, 9.97578527e+00, 1.02575974e+01, + 1.05394467e+01, 1.08213263e+01, 1.11032302e+01, 1.13851530e+01, + 1.16670899e+01, 1.19490366e+01, 1.22309896e+01, 1.25129453e+01, + 1.27949008e+01, 1.30768536e+01, 1.33588011e+01, 1.36407415e+01, + 1.39226727e+01, 1.42045932e+01, 1.44865015e+01, 1.47683963e+01}; const std::vector theta_m_double{ - 2.22044605e-16, 2.58095680e-08, 1.38634787e-05, 3.39716884e-04, - 2.40087636e-03, 9.06565641e-03, 2.38445553e-02, 4.99122887e-02, - 8.95776020e-02, 1.44182976e-01, 2.14235807e-01, 2.99615891e-01, - 3.99777534e-01, 5.13914694e-01, 6.41083523e-01, 7.80287426e-01, - 9.30532846e-01, 1.09086372e+00, 1.26038106e+00, 1.43825260e+00, - 1.62371595e+00, 1.81607782e+00, 2.01471078e+00, 2.21904887e+00, - 2.42858252e+00, 2.64285346e+00, 2.86144963e+00, 3.08400054e+00, - 3.31017284e+00, 3.53966635e+00, 3.77221050e+00, 4.00756109e+00, - 4.24549744e+00, 4.48581986e+00, 4.72834735e+00, 4.97291563e+00, - 5.21937537e+00, 5.46759063e+00, 5.71743745e+00, 5.96880263e+00, - 6.22158266e+00, 6.47568274e+00, 6.73101590e+00, 6.98750228e+00, - 7.24506843e+00, 7.50364669e+00, 7.76317466e+00, 8.02359473e+00, - 8.28485363e+00, 8.54690205e+00, 8.80969427e+00, 9.07318789e+00, - 9.33734351e+00, 9.60212447e+00, 9.86749668e+00, 1.01334283e+01, - 1.03998897e+01, 1.06668532e+01, 1.09342929e+01, 1.12021845e+01, - 1.14705053e+01, 1.17392341e+01, 1.20083509e+01, 1.22778370e+01, - 1.25476748e+01, 1.28178476e+01, 1.30883399e+01, 1.33591369e+01, - 1.36302250e+01, 1.39015909e+01, 1.41732223e+01, 1.44451076e+01, - 1.47172357e+01, 1.49895963e+01, 1.52621795e+01, 1.55349758e+01, - 1.58079765e+01, 1.60811732e+01, 1.63545578e+01, 1.66281227e+01, - 1.69018609e+01, 1.71757655e+01, 1.74498298e+01, 1.77240478e+01, - 1.79984136e+01, 1.82729215e+01, 1.85475662e+01, 1.88223426e+01, - 1.90972458e+01, 1.93722711e+01, 1.96474142e+01, 1.99226707e+01, - 2.01980367e+01, 2.04735082e+01, 2.07490816e+01, 2.10247533e+01, - 2.13005199e+01, 2.15763782e+01, 2.18523250e+01, 2.21283574e+01}; + 2.22044605e-16, 2.58095680e-08, 1.38634787e-05, 3.39716884e-04, + 2.40087636e-03, 9.06565641e-03, 2.38445553e-02, 4.99122887e-02, + 8.95776020e-02, 1.44182976e-01, 2.14235807e-01, 2.99615891e-01, + 3.99777534e-01, 5.13914694e-01, 6.41083523e-01, 7.80287426e-01, + 9.30532846e-01, 1.09086372e+00, 1.26038106e+00, 1.43825260e+00, + 1.62371595e+00, 1.81607782e+00, 2.01471078e+00, 2.21904887e+00, + 2.42858252e+00, 2.64285346e+00, 2.86144963e+00, 3.08400054e+00, + 3.31017284e+00, 3.53966635e+00, 3.77221050e+00, 4.00756109e+00, + 4.24549744e+00, 4.48581986e+00, 4.72834735e+00, 4.97291563e+00, + 5.21937537e+00, 5.46759063e+00, 5.71743745e+00, 5.96880263e+00, + 6.22158266e+00, 6.47568274e+00, 6.73101590e+00, 6.98750228e+00, + 7.24506843e+00, 7.50364669e+00, 7.76317466e+00, 8.02359473e+00, + 8.28485363e+00, 8.54690205e+00, 8.80969427e+00, 9.07318789e+00, + 9.33734351e+00, 9.60212447e+00, 9.86749668e+00, 1.01334283e+01, + 1.03998897e+01, 1.06668532e+01, 1.09342929e+01, 1.12021845e+01, + 1.14705053e+01, 1.17392341e+01, 1.20083509e+01, 1.22778370e+01, + 1.25476748e+01, 1.28178476e+01, 1.30883399e+01, 1.33591369e+01, + 1.36302250e+01, 1.39015909e+01, 1.41732223e+01, 1.44451076e+01, + 1.47172357e+01, 1.49895963e+01, 1.52621795e+01, 1.55349758e+01, + 1.58079765e+01, 1.60811732e+01, 1.63545578e+01, 1.66281227e+01, + 1.69018609e+01, 1.71757655e+01, 1.74498298e+01, 1.77240478e+01, + 1.79984136e+01, 1.82729215e+01, 1.85475662e+01, 1.88223426e+01, + 1.90972458e+01, 1.93722711e+01, 1.96474142e+01, 1.99226707e+01, + 2.01980367e+01, 2.04735082e+01, 2.07490816e+01, 2.10247533e+01, + 2.13005199e+01, 2.15763782e+01, 2.18523250e+01, 2.21283574e+01}; public: /** @@ -92,7 +93,7 @@ class matrix_exp_action_handler { int m, s; set_approx_order(A, b, t, m, s); - double eta = exp(t * mu/s); + double eta = exp(t * mu / s); Eigen::MatrixXd f = b; Eigen::MatrixXd bi = b; @@ -103,10 +104,11 @@ class matrix_exp_action_handler { // maximum absolute row sum, aka inf-norm of matrix operator double c1 = MAT_OP_INF_NORM(bi); for (auto j = 0; j < m; ++j) { - bi = (t/(s * (j + 1))) * (A * bi); + bi = (t / (s * (j + 1))) * (A * bi); f += bi; double c2 = MAT_OP_INF_NORM(bi); - if (c1 + c2 < tol * MAT_OP_INF_NORM(f)) break; + if (c1 + c2 < tol * MAT_OP_INF_NORM(f)) + break; c1 = c2; } f *= eta; @@ -118,9 +120,9 @@ class matrix_exp_action_handler { return f; } - /** + /** * Estimate the 1-norm of mat^m - * + * * See A. H. Al-Mohy and N. J. Higham, A New Scaling and Squaring * Algorithm for the Matrix Exponential, SIAM J. Matrix Anal. Appl. 31(3): * 970-989, 2009. @@ -135,8 +137,7 @@ class matrix_exp_action_handler { * @param mat matrix * @param m power */ - template * = nullptr, + template * = nullptr, require_all_st_same* = nullptr> double mat_power_1_norm(const EigMat1& mat, int m) { if ((mat.array() > 0.0).all()) { @@ -171,7 +172,7 @@ class matrix_exp_action_handler { inline void set_approx_order(const EigMat1& mat, const EigMat2& b, const double& t, int& m, int& s) { // L1 norm - double normA = mat.colwise(). template lpNorm<1>().maxCoeff(); + double normA = mat.colwise().template lpNorm<1>().maxCoeff(); if (normA < tol || t < tol) { m = 0; @@ -180,24 +181,25 @@ class matrix_exp_action_handler { const std::vector& theta = theta_m_double; Eigen::VectorXd alpha(p_max - 1); - if (normA < 4.0*theta[m_max]*p_max*(p_max + 3)/(m_max*b.cols())) { - alpha = Eigen::VectorXd::Constant(p_max-1, normA); + if (normA + < 4.0 * theta[m_max] * p_max * (p_max + 3) / (m_max * b.cols())) { + alpha = Eigen::VectorXd::Constant(p_max - 1, normA); } else { Eigen::VectorXd eta(p_max); for (auto p = 0; p < p_max; ++p) { - double c = mat_power_1_norm(mat, p+2); - c = std::pow(c, 1.0/(p+2.0)); + double c = mat_power_1_norm(mat, p + 2); + c = std::pow(c, 1.0 / (p + 2.0)); eta[p] = c; } for (auto p = 0; p < p_max - 1; ++p) { - alpha[p] = std::max(eta[p],eta(p+1)); + alpha[p] = std::max(eta[p], eta(p + 1)); } } - Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(p_max-1, m_max); + Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(p_max - 1, m_max); for (auto p = 1; p < p_max; ++p) { - for (auto i = p*(p+1)-2; i < m_max; ++i) { - mt(p-1, i) = alpha[p-1]/theta[i]; + for (auto i = p * (p + 1) - 2; i < m_max; ++i) { + mt(p - 1, i) = alpha[p - 1] / theta[i]; } } @@ -206,16 +208,17 @@ class matrix_exp_action_handler { u(i) = i + 1; } - Eigen::MatrixXd c = stan::math::ceil(mt) * u.asDiagonal(); + Eigen::MatrixXd c = stan::math::ceil(mt) * u.asDiagonal(); for (auto i = 0; i < c.size(); ++i) { if (c(i) <= 1.e-16) { c(i) = std::numeric_limits::infinity(); } } - + int cost = c.colwise().minCoeff().minCoeff(&m); - if (std::isinf(cost)) cost = 0; - s = std::max(cost/m, 1); + if (std::isinf(cost)) + cost = 0; + s = std::max(cost / m, 1); } } }; diff --git a/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp b/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp index 558980371a5..c8d6c6a1b1d 100644 --- a/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp +++ b/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp @@ -23,18 +23,17 @@ inline void test_matrix_exp_multiply() { // matrix_exp_multiply Eigen::Matrix res = matrix_exp_multiply(A, B); - EXPECT_MATRIX_FLOAT_EQ(expAB, res); } TEST(MathMatrixPrimMat, matrix_exp_multiply_power_1_norm_fun) { stan::math::matrix_exp_action_handler maexp; - Eigen::MatrixXd x(2,2); + Eigen::MatrixXd x(2, 2); x << 1, -2, -3, 4; EXPECT_FLOAT_EQ(maexp.mat_power_1_norm(x, 2), 32.0); EXPECT_FLOAT_EQ(maexp.mat_power_1_norm(x, 3), 172.0); - Eigen::MatrixXd y(3,3); + Eigen::MatrixXd y(3, 3); y << 1, 2, 3, 4, 5, 6, 7, 1, 2; EXPECT_FLOAT_EQ(maexp.mat_power_1_norm(y, 3), 1163.0); EXPECT_FLOAT_EQ(maexp.mat_power_1_norm(y, 12), 8.3805595e+11); @@ -42,7 +41,7 @@ TEST(MathMatrixPrimMat, matrix_exp_multiply_power_1_norm_fun) { TEST(MathMatrixPrimMat, matrix_exp_multiply_approx_order) { stan::math::matrix_exp_action_handler maexp; - Eigen::MatrixXd x(2,2); + Eigen::MatrixXd x(2, 2); x << 1, -2, -3, 4; Eigen::VectorXd b(2); b << 1.5, 2.5; @@ -63,9 +62,7 @@ TEST(MathMatrixPrimMat, matrix_exp_multiply) { TEST(MathMatrixPrimMat, matrix_exp_multiply_issue_2529) { // issue #2529 https://github.com/stan-dev/math/issues/2529 Eigen::MatrixXd a(3, 3); - a << -1.2, 1.2, 0, - 0, -0.5, 0.5, - 0, 0, -0.3; + a << -1.2, 1.2, 0, 0, -0.5, 0.5, 0, 0, -0.3; Eigen::VectorXd b(3); b << 1.0, 1.0, 1.0; for (auto i = 15; i < 40; ++i) { @@ -77,17 +74,17 @@ TEST(MathMatrixPrimMat, matrix_exp_multiply_issue_2529) { } TEST(MathMatrixPrimMat, matrix_exp_multiply_poisson_5) { -// Block tridiag mat by discretizing Poisson's Eq. on a 5x5 grid + // Block tridiag mat by discretizing Poisson's Eq. on a 5x5 grid int n = 5; int nd = n * n; Eigen::MatrixXd m = Eigen::MatrixXd::Zero(nd, nd); Eigen::VectorXd b(nd); for (auto i = 0; i < nd; ++i) { - b(i) = -1.0 + i * 2.0/(nd - 1); + b(i) = -1.0 + i * 2.0 / (nd - 1); m(i, i) = 4.0; if (i + 1 < nd) { - m(i, i+1) = -1.0; - m(i+1, i) = -1.0; + m(i, i + 1) = -1.0; + m(i + 1, i) = -1.0; } if (i + 5 < nd) { m(i, std::min(i + 5, nd - 1)) = -1.0; From c15984aa96b1384a14413ca3c57eefc0e9bc5c42 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Tue, 23 Nov 2021 09:53:27 -0800 Subject: [PATCH 04/28] add matrix model from #1146 to unit test. --- .../fun/scale_matrix_exp_multiply_test.cpp | 47 +++++++++++++++++-- 1 file changed, 43 insertions(+), 4 deletions(-) diff --git a/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp b/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp index 1dd1bf30a70..834a79f2671 100644 --- a/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp +++ b/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp @@ -23,10 +23,7 @@ inline void test_scale_matrix_exp_multiply() { // matrix_exp_multiply const double t = 1.0; Eigen::Matrix res = scale_matrix_exp_multiply(t, A, B); - EXPECT_EQ(res.size(), expAB.size()); - for (int l = 0; l < res.size(); ++l) { - EXPECT_FLOAT_EQ(res(l), expAB(l)); - } + EXPECT_MATRIX_FLOAT_EQ(expAB, res); } TEST(MathMatrixPrimMat, scale_matrix_exp_multiply) { @@ -48,6 +45,48 @@ TEST(MathMatrixPrimMat, scale_matrix_exp_multiply) { test_scale_matrix_exp_multiply<20, 2>(); } +TEST(MathMatrixPrimMat, scale_matrix_exp_multiply_issue_1146) { + // matrix from bug report + // https://github.com/stan-dev/math/issues/1146 + int m = 17; + Eigen::MatrixXd Ac(m, m); + Ac << + -4.4758400, 0.0732692, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0183173, -0.0366346, 0.0183173, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0146538, -0.1794070, 0.000000, 0.000000, 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, 0.00000000, 0.0000000, 0.0426379, 0.0249012, 0.00748033, 0.00498689, 0.00415574, + 0.0000000, 0.0000000, 0.0000000, -5.868960, 0.447212, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.223606, -0.447212, 0.2236060, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.3287580, 0.000000, 0.447212, -2.3328900, 0.0000, 0.04931380, 0.000000, 0.4931380, 0.00000000, 0.0000000, 0.3287580, 0.6857140, 0.05767690, 0.03845130, 0.03204270, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, -42.1500, 11.43980000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 2.3972000, 0.000000, 0.000000, 1.5102300, 11.4398, -26.34030000, 0.000000, 3.5957900, 0.00000000, 0.0000000, 2.3972000, 5.0000000, 0.42056100, 0.28037400, 0.23364500, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, -0.122394, 0.1223940, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 1.0958600, 0.000000, 0.000000, 0.6903930, 0.0000, 0.16437900, 1.101550, -6.4337600, 0.00000000, 0.0000000, 1.0958600, 2.2857100, 0.19225600, 0.12817100, 0.10680900, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, -0.00471794, 0.0014969, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.01901160, -0.0402984, 0.0212869, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0426379, 0.000000, 0.000000, 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, 0.00000000, 0.0170295, -0.3347470, 0.1778660, 0.00748033, 0.00498689, 0.00415574, + 0.0000000, 0.0000000, 0.2695880, 0.000000, 0.000000, 0.6065730, 0.0000, 0.14442200, 0.000000, 1.4442200, 0.00000000, 0.0000000, 1.9256300, -4.3904300, 0.03851260, 0.15405000, 0.19256300, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, -1.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, -1.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, -1.00000000; + + Eigen::MatrixXd phi = Eigen::MatrixXd::Zero(2 * m, 2 * m); + Eigen::MatrixXd OI = Eigen::MatrixXd::Zero(2 * m, m); + + OI.bottomRightCorner(m, m) = Eigen::VectorXd::Constant(m, 1.0).asDiagonal(); + phi.topLeftCorner(m, m) = Ac; + phi.topRightCorner(m, m) = Eigen::VectorXd::Constant(m, 0.1).asDiagonal(); + phi.bottomRightCorner(m, m) = -Ac.transpose(); + + Eigen::MatrixXd expAB(2 * m, m); + Eigen::MatrixXd res(2 * m, m); + for (auto i = 0; i < 20; ++i) { + double dt = 0.1 + i * 0.05; + expAB = stan::math::matrix_exp(phi * dt) * OI; + res = stan::math::scale_matrix_exp_multiply(dt, phi, OI); + EXPECT_MATRIX_FLOAT_EQ(expAB, res); + } +} + TEST(MathMatrixPrimMat, scale_matrix_exp_multiply_exception) { using stan::math::scale_matrix_exp_multiply; const double t = 1.0; From 564757e506a30f749cc66eefdfa4ab87b2a437f2 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 23 Nov 2021 18:00:22 +0000 Subject: [PATCH 05/28] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- .../fun/scale_matrix_exp_multiply_test.cpp | 66 +++++++++++++------ 1 file changed, 47 insertions(+), 19 deletions(-) diff --git a/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp b/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp index 834a79f2671..cfd75abb990 100644 --- a/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp +++ b/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp @@ -50,28 +50,56 @@ TEST(MathMatrixPrimMat, scale_matrix_exp_multiply_issue_1146) { // https://github.com/stan-dev/math/issues/1146 int m = 17; Eigen::MatrixXd Ac(m, m); - Ac << - -4.4758400, 0.0732692, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0183173, -0.0366346, 0.0183173, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0146538, -0.1794070, 0.000000, 0.000000, 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, 0.00000000, 0.0000000, 0.0426379, 0.0249012, 0.00748033, 0.00498689, 0.00415574, - 0.0000000, 0.0000000, 0.0000000, -5.868960, 0.447212, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.223606, -0.447212, 0.2236060, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.3287580, 0.000000, 0.447212, -2.3328900, 0.0000, 0.04931380, 0.000000, 0.4931380, 0.00000000, 0.0000000, 0.3287580, 0.6857140, 0.05767690, 0.03845130, 0.03204270, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, -42.1500, 11.43980000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 2.3972000, 0.000000, 0.000000, 1.5102300, 11.4398, -26.34030000, 0.000000, 3.5957900, 0.00000000, 0.0000000, 2.3972000, 5.0000000, 0.42056100, 0.28037400, 0.23364500, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, -0.122394, 0.1223940, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 1.0958600, 0.000000, 0.000000, 0.6903930, 0.0000, 0.16437900, 1.101550, -6.4337600, 0.00000000, 0.0000000, 1.0958600, 2.2857100, 0.19225600, 0.12817100, 0.10680900, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, -0.00471794, 0.0014969, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.01901160, -0.0402984, 0.0212869, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0426379, 0.000000, 0.000000, 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, 0.00000000, 0.0170295, -0.3347470, 0.1778660, 0.00748033, 0.00498689, 0.00415574, - 0.0000000, 0.0000000, 0.2695880, 0.000000, 0.000000, 0.6065730, 0.0000, 0.14442200, 0.000000, 1.4442200, 0.00000000, 0.0000000, 1.9256300, -4.3904300, 0.03851260, 0.15405000, 0.19256300, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, -1.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, -1.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, -1.00000000; + Ac << -4.4758400, 0.0732692, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, + 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.00000000, 0.00000000, 0.00000000, 0.0183173, -0.0366346, + 0.0183173, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, + 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, 0.0000000, 0.0146538, -0.1794070, 0.000000, + 0.000000, 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, 0.00000000, + 0.0000000, 0.0426379, 0.0249012, 0.00748033, 0.00498689, 0.00415574, + 0.0000000, 0.0000000, 0.0000000, -5.868960, 0.447212, 0.0000000, 0.0000, + 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.00000000, 0.00000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.223606, -0.447212, 0.2236060, 0.0000, 0.00000000, 0.000000, + 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, 0.0000000, 0.0000000, 0.3287580, 0.000000, + 0.447212, -2.3328900, 0.0000, 0.04931380, 0.000000, 0.4931380, 0.00000000, + 0.0000000, 0.3287580, 0.6857140, 0.05767690, 0.03845130, 0.03204270, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, -42.1500, + 11.43980000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.00000000, 0.00000000, 0.00000000, 0.0000000, 0.0000000, + 2.3972000, 0.000000, 0.000000, 1.5102300, 11.4398, -26.34030000, 0.000000, + 3.5957900, 0.00000000, 0.0000000, 2.3972000, 5.0000000, 0.42056100, + 0.28037400, 0.23364500, 0.0000000, 0.0000000, 0.0000000, 0.000000, + 0.000000, 0.0000000, 0.0000, 0.00000000, -0.122394, 0.1223940, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 1.0958600, 0.000000, 0.000000, 0.6903930, 0.0000, + 0.16437900, 1.101550, -6.4337600, 0.00000000, 0.0000000, 1.0958600, + 2.2857100, 0.19225600, 0.12817100, 0.10680900, 0.0000000, 0.0000000, + 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, + 0.0000000, -0.00471794, 0.0014969, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.000000, + 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.01901160, + -0.0402984, 0.0212869, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0426379, 0.000000, 0.000000, 0.0268619, 0.0000, + 0.00639568, 0.000000, 0.0639568, 0.00000000, 0.0170295, -0.3347470, + 0.1778660, 0.00748033, 0.00498689, 0.00415574, 0.0000000, 0.0000000, + 0.2695880, 0.000000, 0.000000, 0.6065730, 0.0000, 0.14442200, 0.000000, + 1.4442200, 0.00000000, 0.0000000, 1.9256300, -4.3904300, 0.03851260, + 0.15405000, 0.19256300, 0.0000000, 0.0000000, 0.0000000, 0.000000, + 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, -1.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, + 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.00000000, -1.00000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, + 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, -1.00000000; Eigen::MatrixXd phi = Eigen::MatrixXd::Zero(2 * m, 2 * m); Eigen::MatrixXd OI = Eigen::MatrixXd::Zero(2 * m, m); - + OI.bottomRightCorner(m, m) = Eigen::VectorXd::Constant(m, 1.0).asDiagonal(); phi.topLeftCorner(m, m) = Ac; phi.topRightCorner(m, m) = Eigen::VectorXd::Constant(m, 0.1).asDiagonal(); From 4b78bf6deca538ee2213958b0ae00dacf377cfaa Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Tue, 23 Nov 2021 10:18:49 -0800 Subject: [PATCH 06/28] linting --- .../fun/scale_matrix_exp_multiply_test.cpp | 85 +++++++++++++++---- 1 file changed, 68 insertions(+), 17 deletions(-) diff --git a/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp b/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp index 834a79f2671..49f57172cf6 100644 --- a/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp +++ b/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp @@ -51,23 +51,74 @@ TEST(MathMatrixPrimMat, scale_matrix_exp_multiply_issue_1146) { int m = 17; Eigen::MatrixXd Ac(m, m); Ac << - -4.4758400, 0.0732692, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0183173, -0.0366346, 0.0183173, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0146538, -0.1794070, 0.000000, 0.000000, 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, 0.00000000, 0.0000000, 0.0426379, 0.0249012, 0.00748033, 0.00498689, 0.00415574, - 0.0000000, 0.0000000, 0.0000000, -5.868960, 0.447212, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.223606, -0.447212, 0.2236060, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.3287580, 0.000000, 0.447212, -2.3328900, 0.0000, 0.04931380, 0.000000, 0.4931380, 0.00000000, 0.0000000, 0.3287580, 0.6857140, 0.05767690, 0.03845130, 0.03204270, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, -42.1500, 11.43980000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 2.3972000, 0.000000, 0.000000, 1.5102300, 11.4398, -26.34030000, 0.000000, 3.5957900, 0.00000000, 0.0000000, 2.3972000, 5.0000000, 0.42056100, 0.28037400, 0.23364500, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, -0.122394, 0.1223940, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 1.0958600, 0.000000, 0.000000, 0.6903930, 0.0000, 0.16437900, 1.101550, -6.4337600, 0.00000000, 0.0000000, 1.0958600, 2.2857100, 0.19225600, 0.12817100, 0.10680900, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, -0.00471794, 0.0014969, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.01901160, -0.0402984, 0.0212869, 0.0000000, 0.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0426379, 0.000000, 0.000000, 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, 0.00000000, 0.0170295, -0.3347470, 0.1778660, 0.00748033, 0.00498689, 0.00415574, - 0.0000000, 0.0000000, 0.2695880, 0.000000, 0.000000, 0.6065730, 0.0000, 0.14442200, 0.000000, 1.4442200, 0.00000000, 0.0000000, 1.9256300, -4.3904300, 0.03851260, 0.15405000, 0.19256300, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, -1.00000000, 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, -1.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, -1.00000000; + -4.4758400, 0.0732692, 0.0000000, 0.000000, 0.000000, + 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, + 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, + 0.0183173, -0.0366346, 0.0183173, 0.000000, 0.000000, + 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, + 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, + 0.0000000, 0.0146538, -0.1794070, 0.000000, 0.000000, + 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, + 0.00000000, 0.0000000, 0.0426379, 0.0249012, 0.00748033, + 0.00498689, 0.00415574, + 0.0000000, 0.0000000, 0.0000000, -5.868960, 0.447212, + 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, + 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.223606, -0.447212, + 0.2236060, 0.0000, 0.00000000, 0.000000, 0.0000000, + 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.3287580, 0.000000, 0.447212, + -2.3328900, 0.0000, 0.04931380, 0.000000, 0.4931380, + 0.00000000, 0.0000000, 0.3287580, 0.6857140, 0.05767690, + 0.03845130, 0.03204270, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, + 0.0000000, -42.1500, 11.43980000, 0.000000, 0.0000000, + 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 2.3972000, 0.000000, 0.000000, + 1.5102300, 11.4398, -26.34030000, 0.000000, 3.5957900, + 0.00000000, 0.0000000, 2.3972000, 5.0000000, 0.42056100, + 0.28037400, 0.23364500, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, + 0.0000000, 0.0000, 0.00000000, -0.122394, 0.1223940, + 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 1.0958600, 0.000000, 0.000000, + 0.6903930, 0.0000, 0.16437900, 1.101550, -6.4337600, + 0.00000000, 0.0000000, 1.0958600, 2.2857100, 0.19225600, + 0.12817100, 0.10680900, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, + 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, + -0.00471794, 0.0014969, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, + 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, + 0.01901160, -0.0402984, 0.0212869, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0426379, 0.000000, 0.000000, + 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, + 0.00000000, 0.0170295, -0.3347470, 0.1778660, 0.00748033, + 0.00498689, 0.00415574, + 0.0000000, 0.0000000, 0.2695880, 0.000000, 0.000000, + 0.6065730, 0.0000, 0.14442200, 0.000000, 1.4442200, + 0.00000000, 0.0000000, 1.9256300, -4.3904300, 0.03851260, + 0.15405000, 0.19256300, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, + 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, + 0.00000000, 0.0000000, 0.0000000, 0.0000000, -1.00000000, + 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, + 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, + 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + -1.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, + 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, + 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, -1.00000000; Eigen::MatrixXd phi = Eigen::MatrixXd::Zero(2 * m, 2 * m); Eigen::MatrixXd OI = Eigen::MatrixXd::Zero(2 * m, m); From 49dc3e5595dd53f26de2e071f19bcde2c40a698d Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 23 Nov 2021 18:23:55 +0000 Subject: [PATCH 07/28] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- .../fun/scale_matrix_exp_multiply_test.cpp | 115 +++++++----------- 1 file changed, 46 insertions(+), 69 deletions(-) diff --git a/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp b/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp index 935b6ada8fb..cfd75abb990 100644 --- a/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp +++ b/test/unit/math/prim/fun/scale_matrix_exp_multiply_test.cpp @@ -50,75 +50,52 @@ TEST(MathMatrixPrimMat, scale_matrix_exp_multiply_issue_1146) { // https://github.com/stan-dev/math/issues/1146 int m = 17; Eigen::MatrixXd Ac(m, m); - Ac << - -4.4758400, 0.0732692, 0.0000000, 0.000000, 0.000000, - 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, - 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, - 0.00000000, 0.00000000, - 0.0183173, -0.0366346, 0.0183173, 0.000000, 0.000000, - 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, - 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, - 0.00000000, 0.00000000, - 0.0000000, 0.0146538, -0.1794070, 0.000000, 0.000000, - 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, - 0.00000000, 0.0000000, 0.0426379, 0.0249012, 0.00748033, - 0.00498689, 0.00415574, - 0.0000000, 0.0000000, 0.0000000, -5.868960, 0.447212, - 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, - 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, - 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.223606, -0.447212, - 0.2236060, 0.0000, 0.00000000, 0.000000, 0.0000000, - 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, - 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.3287580, 0.000000, 0.447212, - -2.3328900, 0.0000, 0.04931380, 0.000000, 0.4931380, - 0.00000000, 0.0000000, 0.3287580, 0.6857140, 0.05767690, - 0.03845130, 0.03204270, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, - 0.0000000, -42.1500, 11.43980000, 0.000000, 0.0000000, - 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, - 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 2.3972000, 0.000000, 0.000000, - 1.5102300, 11.4398, -26.34030000, 0.000000, 3.5957900, - 0.00000000, 0.0000000, 2.3972000, 5.0000000, 0.42056100, - 0.28037400, 0.23364500, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, - 0.0000000, 0.0000, 0.00000000, -0.122394, 0.1223940, - 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, - 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 1.0958600, 0.000000, 0.000000, - 0.6903930, 0.0000, 0.16437900, 1.101550, -6.4337600, - 0.00000000, 0.0000000, 1.0958600, 2.2857100, 0.19225600, - 0.12817100, 0.10680900, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, - 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, - -0.00471794, 0.0014969, 0.0000000, 0.0000000, 0.00000000, - 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, - 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, - 0.01901160, -0.0402984, 0.0212869, 0.0000000, 0.00000000, - 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0426379, 0.000000, 0.000000, - 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, - 0.00000000, 0.0170295, -0.3347470, 0.1778660, 0.00748033, - 0.00498689, 0.00415574, - 0.0000000, 0.0000000, 0.2695880, 0.000000, 0.000000, - 0.6065730, 0.0000, 0.14442200, 0.000000, 1.4442200, - 0.00000000, 0.0000000, 1.9256300, -4.3904300, 0.03851260, - 0.15405000, 0.19256300, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, - 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, - 0.00000000, 0.0000000, 0.0000000, 0.0000000, -1.00000000, - 0.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, - 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, - 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, - -1.00000000, 0.00000000, - 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, - 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, - 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, - 0.00000000, -1.00000000; + Ac << -4.4758400, 0.0732692, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, + 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.00000000, 0.00000000, 0.00000000, 0.0183173, -0.0366346, + 0.0183173, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, + 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, 0.0000000, 0.0146538, -0.1794070, 0.000000, + 0.000000, 0.0268619, 0.0000, 0.00639568, 0.000000, 0.0639568, 0.00000000, + 0.0000000, 0.0426379, 0.0249012, 0.00748033, 0.00498689, 0.00415574, + 0.0000000, 0.0000000, 0.0000000, -5.868960, 0.447212, 0.0000000, 0.0000, + 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.00000000, 0.00000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.223606, -0.447212, 0.2236060, 0.0000, 0.00000000, 0.000000, + 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, 0.0000000, 0.0000000, 0.3287580, 0.000000, + 0.447212, -2.3328900, 0.0000, 0.04931380, 0.000000, 0.4931380, 0.00000000, + 0.0000000, 0.3287580, 0.6857140, 0.05767690, 0.03845130, 0.03204270, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, -42.1500, + 11.43980000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.00000000, 0.00000000, 0.00000000, 0.0000000, 0.0000000, + 2.3972000, 0.000000, 0.000000, 1.5102300, 11.4398, -26.34030000, 0.000000, + 3.5957900, 0.00000000, 0.0000000, 2.3972000, 5.0000000, 0.42056100, + 0.28037400, 0.23364500, 0.0000000, 0.0000000, 0.0000000, 0.000000, + 0.000000, 0.0000000, 0.0000, 0.00000000, -0.122394, 0.1223940, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 1.0958600, 0.000000, 0.000000, 0.6903930, 0.0000, + 0.16437900, 1.101550, -6.4337600, 0.00000000, 0.0000000, 1.0958600, + 2.2857100, 0.19225600, 0.12817100, 0.10680900, 0.0000000, 0.0000000, + 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, + 0.0000000, -0.00471794, 0.0014969, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.000000, + 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.01901160, + -0.0402984, 0.0212869, 0.0000000, 0.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0426379, 0.000000, 0.000000, 0.0268619, 0.0000, + 0.00639568, 0.000000, 0.0639568, 0.00000000, 0.0170295, -0.3347470, + 0.1778660, 0.00748033, 0.00498689, 0.00415574, 0.0000000, 0.0000000, + 0.2695880, 0.000000, 0.000000, 0.6065730, 0.0000, 0.14442200, 0.000000, + 1.4442200, 0.00000000, 0.0000000, 1.9256300, -4.3904300, 0.03851260, + 0.15405000, 0.19256300, 0.0000000, 0.0000000, 0.0000000, 0.000000, + 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, 0.0000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, -1.00000000, 0.00000000, 0.00000000, + 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, + 0.00000000, 0.000000, 0.0000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.00000000, -1.00000000, 0.00000000, 0.0000000, 0.0000000, + 0.0000000, 0.000000, 0.000000, 0.0000000, 0.0000, 0.00000000, 0.000000, + 0.0000000, 0.00000000, 0.0000000, 0.0000000, 0.0000000, 0.00000000, + 0.00000000, -1.00000000; Eigen::MatrixXd phi = Eigen::MatrixXd::Zero(2 * m, 2 * m); Eigen::MatrixXd OI = Eigen::MatrixXd::Zero(2 * m, m); From bf2e7452fcbbded0f459fb214fb8054b160cc6db Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Tue, 23 Nov 2021 10:32:34 -0800 Subject: [PATCH 08/28] fix linting --- test/unit/math/prim/fun/matrix_exp_multiply_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp b/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp index c8d6c6a1b1d..2ab0da7f7e0 100644 --- a/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp +++ b/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp @@ -66,7 +66,7 @@ TEST(MathMatrixPrimMat, matrix_exp_multiply_issue_2529) { Eigen::VectorXd b(3); b << 1.0, 1.0, 1.0; for (auto i = 15; i < 40; ++i) { - double t = double(i); + double t = i; Eigen::MatrixXd m1 = stan::math::matrix_exp_multiply(a * t, b); Eigen::MatrixXd m2 = stan::math::matrix_exp(a * t) * b; EXPECT_MATRIX_FLOAT_EQ(m1, m2); From f708734d34ac33275079bef8bdb2432e4eea5062 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Tue, 23 Nov 2021 11:58:30 -0800 Subject: [PATCH 09/28] fix unit test --- .../unit/math/prim/fun/matrix_exp_action_handler_test.cpp | 8 -------- 1 file changed, 8 deletions(-) diff --git a/test/unit/math/prim/fun/matrix_exp_action_handler_test.cpp b/test/unit/math/prim/fun/matrix_exp_action_handler_test.cpp index c240df6bc54..78e0ed48f2f 100644 --- a/test/unit/math/prim/fun/matrix_exp_action_handler_test.cpp +++ b/test/unit/math/prim/fun/matrix_exp_action_handler_test.cpp @@ -68,14 +68,6 @@ TEST(MathMatrixRevMat, matrix_exp_action_vector) { for (size_t i = 0; i < n; ++i) { EXPECT_NEAR(res(i), expb(i), 1.e-6); } - - int m1, s1, m2, s2; - const double t1 = 9.9, t2 = 1.0; - handler.set_approximation_parameter(A, t1, m1, s1); - A *= t1; - handler.set_approximation_parameter(A, t2, m2, s2); - EXPECT_EQ(m1, m2); - EXPECT_EQ(s1, s2); } } From 1081352b5f5f7c053d69570bc5c67319a83d0bc7 Mon Sep 17 00:00:00 2001 From: rok-cesnovar Date: Mon, 29 Nov 2021 19:17:33 +0100 Subject: [PATCH 10/28] fix multiple evals of b --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index e8f26664363..ca20c57377c 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -95,8 +95,9 @@ class matrix_exp_action_handler { double eta = exp(t * mu / s); - Eigen::MatrixXd f = b; - Eigen::MatrixXd bi = b; + const auto& b_eval = b.eval(); + Eigen::MatrixXd f = b_eval; + Eigen::MatrixXd bi = b_eval; #define MAT_OP_INF_NORM(X) X.cwiseAbs().rowwise().sum().maxCoeff() From ea55ec7d4177d2862ca2a9211cd87158de947891 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Mon, 29 Nov 2021 10:54:33 -0800 Subject: [PATCH 11/28] misc review comments --- .../prim/fun/matrix_exp_action_handler.hpp | 26 ++++++++++--------- 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index ca20c57377c..67a13ad03d7 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -27,7 +27,6 @@ class matrix_exp_action_handler { static constexpr int m_max = 55; static constexpr double tol = 1.0e-16; - // const std::vector theta_m_single{ 1.19209280e-07, 5.97885889e-04, 1.12338647e-02, 5.11661936e-02, 1.30848716e-01, 2.49528932e-01, 4.01458242e-01, 5.80052463e-01, @@ -86,9 +85,7 @@ class matrix_exp_action_handler { const double& t = 1.0) { Eigen::MatrixXd A = mat; double mu = A.trace() / A.rows(); - for (int i = 0; i < A.rows(); ++i) { - A(i, i) -= mu; - } + A.diagonal().array() -= mu; int m, s; set_approx_order(A, b, t, m, s); @@ -99,16 +96,14 @@ class matrix_exp_action_handler { Eigen::MatrixXd f = b_eval; Eigen::MatrixXd bi = b_eval; -#define MAT_OP_INF_NORM(X) X.cwiseAbs().rowwise().sum().maxCoeff() - for (auto i = 0; i < s; ++i) { // maximum absolute row sum, aka inf-norm of matrix operator - double c1 = MAT_OP_INF_NORM(bi); + double c1 = matrix_operator_inf_norm(bi); for (auto j = 0; j < m; ++j) { bi = (t / (s * (j + 1))) * (A * bi); f += bi; - double c2 = MAT_OP_INF_NORM(bi); - if (c1 + c2 < tol * MAT_OP_INF_NORM(f)) + double c2 = matrix_operator_inf_norm(bi); + if (c1 + c2 < tol * matrix_operator_inf_norm(f)) break; c1 = c2; } @@ -116,11 +111,18 @@ class matrix_exp_action_handler { bi = f; } -#undef MAT_OP_INF_NORM - return f; } + /** + * Eigen expression for matrix operator infinity norm + * + * @param mat matrix + */ + double matrix_operator_inf_norm(Eigen::MatrixXd const& x) { + return x.cwiseAbs().rowwise().sum().maxCoeff(); + } + /** * Estimate the 1-norm of mat^m * @@ -143,7 +145,7 @@ class matrix_exp_action_handler { double mat_power_1_norm(const EigMat1& mat, int m) { if ((mat.array() > 0.0).all()) { Eigen::VectorXd e = Eigen::VectorXd::Constant(mat.rows(), 1.0); - for (auto j = 0; j < m; ++j) { + for (int j = 0; j < m; ++j) { e = mat.transpose() * e; } return e.lpNorm(); From fac3d0c1b9f5b36998f76f6a788ef7bf25dcb814 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Fri, 3 Dec 2021 09:53:46 -0800 Subject: [PATCH 12/28] additional unit tests for expm of large condition number --- .../prim/fun/matrix_exp_multiply_test.cpp | 46 +++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp b/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp index 2ab0da7f7e0..79036c329db 100644 --- a/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp +++ b/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp @@ -98,6 +98,52 @@ TEST(MathMatrixPrimMat, matrix_exp_multiply_poisson_5) { EXPECT_MATRIX_NEAR(p1, p2, 1.e-12); } +TEST(MathMatrixPrimMat, matrix_exp_multiply_matlab) { + // https://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/#7401fe8a-5a7d-40df-92d7-1ae34f45adf2 // NOLINT + { + Eigen::MatrixXd m(3, 3); + m << + 0, 1e-08, 0, + -20066666666.6667, -3, 20000000000, + 66.6666666666667, 0, -66.6666666666667; + Eigen::MatrixXd expm(3, 3); + expm << + 0.446849468283175, 1.54044157383952e-09, 0.462811453558774, + -5743067.77947947, -0.0152830038686819, -4526542.71278401, + 0.447722977849494, 1.54270484519591e-09, 0.463480648837651; + std::vector r{-1, 0, 1}; + Eigen::VectorXd p1(3), p2(3); + for (auto i : r) { for (auto j : r) { for (auto k : r) { + Eigen::VectorXd b(3); + b << i, j, k; + p1 = stan::math::matrix_exp_multiply(m, b); + p2 = expm * b; + EXPECT_MATRIX_NEAR(p1, p2, 5.e-6); + } + } + } + } + + // https://discourse.mc-stan.org/t/documentation-for-the-action-of-the-matrix-exponential/25498/8?u=yizhang // NOLINT + { + Eigen::MatrixXd m(2, 2); + m << -49, 24, -64, 31; + Eigen::MatrixXd expm(2, 2); + expm << -0.735758758144755, 0.551819099658099, + -1.47151759908826, 1.10363824071558; + std::vector r{0, 1, 10, 100, 1000, 10000}; + Eigen::VectorXd p1(2), p2(2); + for (auto i : r) { for (auto j : r) { + Eigen::VectorXd b(2); + b << i, j; + p1 = stan::math::matrix_exp_multiply(m, b); + p2 = expm * b; + EXPECT_MATRIX_NEAR(p1, p2, 1.e-8); + } + } + } +} + TEST(MathMatrixPrimMat, matrix_exp_multiply_exception) { using stan::math::matrix_exp_multiply; { // multiplicable From f0f462fc73322c6377f4a99ebd7907e294b225c8 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 3 Dec 2021 17:58:09 +0000 Subject: [PATCH 13/28] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- .../prim/fun/matrix_exp_multiply_test.cpp | 30 ++++++++++--------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp b/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp index 79036c329db..b81df9f10a2 100644 --- a/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp +++ b/test/unit/math/prim/fun/matrix_exp_multiply_test.cpp @@ -99,21 +99,21 @@ TEST(MathMatrixPrimMat, matrix_exp_multiply_poisson_5) { } TEST(MathMatrixPrimMat, matrix_exp_multiply_matlab) { - // https://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/#7401fe8a-5a7d-40df-92d7-1ae34f45adf2 // NOLINT + // https://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/#7401fe8a-5a7d-40df-92d7-1ae34f45adf2 + // // NOLINT { Eigen::MatrixXd m(3, 3); - m << - 0, 1e-08, 0, - -20066666666.6667, -3, 20000000000, - 66.6666666666667, 0, -66.6666666666667; + m << 0, 1e-08, 0, -20066666666.6667, -3, 20000000000, 66.6666666666667, 0, + -66.6666666666667; Eigen::MatrixXd expm(3, 3); - expm << - 0.446849468283175, 1.54044157383952e-09, 0.462811453558774, - -5743067.77947947, -0.0152830038686819, -4526542.71278401, - 0.447722977849494, 1.54270484519591e-09, 0.463480648837651; + expm << 0.446849468283175, 1.54044157383952e-09, 0.462811453558774, + -5743067.77947947, -0.0152830038686819, -4526542.71278401, + 0.447722977849494, 1.54270484519591e-09, 0.463480648837651; std::vector r{-1, 0, 1}; Eigen::VectorXd p1(3), p2(3); - for (auto i : r) { for (auto j : r) { for (auto k : r) { + for (auto i : r) { + for (auto j : r) { + for (auto k : r) { Eigen::VectorXd b(3); b << i, j, k; p1 = stan::math::matrix_exp_multiply(m, b); @@ -124,16 +124,18 @@ TEST(MathMatrixPrimMat, matrix_exp_multiply_matlab) { } } - // https://discourse.mc-stan.org/t/documentation-for-the-action-of-the-matrix-exponential/25498/8?u=yizhang // NOLINT + // https://discourse.mc-stan.org/t/documentation-for-the-action-of-the-matrix-exponential/25498/8?u=yizhang + // // NOLINT { Eigen::MatrixXd m(2, 2); m << -49, 24, -64, 31; Eigen::MatrixXd expm(2, 2); - expm << -0.735758758144755, 0.551819099658099, - -1.47151759908826, 1.10363824071558; + expm << -0.735758758144755, 0.551819099658099, -1.47151759908826, + 1.10363824071558; std::vector r{0, 1, 10, 100, 1000, 10000}; Eigen::VectorXd p1(2), p2(2); - for (auto i : r) { for (auto j : r) { + for (auto i : r) { + for (auto j : r) { Eigen::VectorXd b(2); b << i, j; p1 = stan::math::matrix_exp_multiply(m, b); From 8d9bbad075ae21623f5d351fec78d7af3af071ee Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 10:37:30 -0400 Subject: [PATCH 14/28] matrix_exp_action_handler: adding definition of constexpr variables --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index 67a13ad03d7..72ef95daa76 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -226,6 +226,11 @@ class matrix_exp_action_handler { } }; +// C++-11 requires definition of the static member in addition +// to the declaration. This is fixed in C++-17. +constexpr int matrix_exp_action_handler::p_max; +constexpr int matrix_exp_action_handler::m_max; + } // namespace math } // namespace stan From d6777b02ba9aefbee58a889cec9dbf7c1b618451 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 19 May 2023 10:38:35 -0400 Subject: [PATCH 15/28] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index 72ef95daa76..b4b7c53f852 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -230,7 +230,7 @@ class matrix_exp_action_handler { // to the declaration. This is fixed in C++-17. constexpr int matrix_exp_action_handler::p_max; constexpr int matrix_exp_action_handler::m_max; - + } // namespace math } // namespace stan From d018c4c5c3dd9de5aa401a3b86471633bbab714b Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 11:04:10 -0400 Subject: [PATCH 16/28] matrix_exp_action_handler: set_approx_order -- eager return --- .../prim/fun/matrix_exp_action_handler.hpp | 82 ++++++++++--------- 1 file changed, 45 insertions(+), 37 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index b4b7c53f852..a508bc9b6f7 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -1,6 +1,8 @@ #ifndef STAN_MATH_PRIM_FUN_MATRIX_EXP_ACTION_HANDLER_HPP #define STAN_MATH_PRIM_FUN_MATRIX_EXP_ACTION_HANDLER_HPP +#include + #include #include #include @@ -174,55 +176,61 @@ class matrix_exp_action_handler { require_all_st_same* = nullptr> inline void set_approx_order(const EigMat1& mat, const EigMat2& b, const double& t, int& m, int& s) { + if (t < tol) { + m = 0; + s = 1; + return; + } + // L1 norm double normA = mat.colwise().template lpNorm<1>().maxCoeff(); - - if (normA < tol || t < tol) { + if (normA < tol) { m = 0; s = 1; + return; + } + + const std::vector& theta = theta_m_double; + + Eigen::VectorXd alpha(p_max - 1); + if (normA + < 4.0 * theta[m_max] * p_max * (p_max + 3) / (m_max * b.cols())) { + alpha = Eigen::VectorXd::Constant(p_max - 1, normA); } else { - const std::vector& theta = theta_m_double; - Eigen::VectorXd alpha(p_max - 1); - - if (normA - < 4.0 * theta[m_max] * p_max * (p_max + 3) / (m_max * b.cols())) { - alpha = Eigen::VectorXd::Constant(p_max - 1, normA); - } else { - Eigen::VectorXd eta(p_max); - for (auto p = 0; p < p_max; ++p) { - double c = mat_power_1_norm(mat, p + 2); - c = std::pow(c, 1.0 / (p + 2.0)); - eta[p] = c; - } - for (auto p = 0; p < p_max - 1; ++p) { - alpha[p] = std::max(eta[p], eta(p + 1)); - } + Eigen::VectorXd eta(p_max); + for (auto p = 0; p < p_max; ++p) { + double c = mat_power_1_norm(mat, p + 2); + c = std::pow(c, 1.0 / (p + 2.0)); + eta[p] = c; } - - Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(p_max - 1, m_max); - for (auto p = 1; p < p_max; ++p) { - for (auto i = p * (p + 1) - 2; i < m_max; ++i) { - mt(p - 1, i) = alpha[p - 1] / theta[i]; - } + for (auto p = 0; p < p_max - 1; ++p) { + alpha[p] = std::max(eta[p], eta(p + 1)); } + } - Eigen::Matrix u(m_max); - for (auto i = 0; i < m_max; ++i) { - u(i) = i + 1; + Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(p_max - 1, m_max); + for (auto p = 1; p < p_max; ++p) { + for (auto i = p * (p + 1) - 2; i < m_max; ++i) { + mt(p - 1, i) = alpha[p - 1] / theta[i]; } + } - Eigen::MatrixXd c = stan::math::ceil(mt) * u.asDiagonal(); - for (auto i = 0; i < c.size(); ++i) { - if (c(i) <= 1.e-16) { - c(i) = std::numeric_limits::infinity(); - } - } + Eigen::Matrix u(m_max); + for (auto i = 0; i < m_max; ++i) { + u(i) = i + 1; + } - int cost = c.colwise().minCoeff().minCoeff(&m); - if (std::isinf(cost)) - cost = 0; - s = std::max(cost / m, 1); + Eigen::MatrixXd c = stan::math::ceil(mt) * u.asDiagonal(); + for (auto i = 0; i < c.size(); ++i) { + if (c(i) <= 1.e-16) { + c(i) = std::numeric_limits::infinity(); + } } + + int cost = c.colwise().minCoeff().minCoeff(&m); + if (std::isinf(cost)) + cost = 0; + s = std::max(cost / m, 1); } }; From 9ea04de66c57bbd1fbfa41a258dbde479ea06e71 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 11:04:40 -0400 Subject: [PATCH 17/28] matrix_exp_action_handler: set_approx_order -- adjust condition --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index a508bc9b6f7..2e3f6af3c57 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -193,8 +193,8 @@ class matrix_exp_action_handler { const std::vector& theta = theta_m_double; Eigen::VectorXd alpha(p_max - 1); - if (normA - < 4.0 * theta[m_max] * p_max * (p_max + 3) / (m_max * b.cols())) { + if (normA * (m_max * b.cols()) + < 4.0 * theta[m_max] * p_max * (p_max + 3)) { alpha = Eigen::VectorXd::Constant(p_max - 1, normA); } else { Eigen::VectorXd eta(p_max); From 079c66e7e253a62c53a5610a081d7d5fbf48a5f5 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 11:06:05 -0400 Subject: [PATCH 18/28] matrix_exp_action_handler: set_approx_order -- replace theta --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index 2e3f6af3c57..f53769ee921 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -190,11 +190,9 @@ class matrix_exp_action_handler { return; } - const std::vector& theta = theta_m_double; - Eigen::VectorXd alpha(p_max - 1); if (normA * (m_max * b.cols()) - < 4.0 * theta[m_max] * p_max * (p_max + 3)) { + < 4.0 * theta_m_double[m_max] * p_max * (p_max + 3)) { alpha = Eigen::VectorXd::Constant(p_max - 1, normA); } else { Eigen::VectorXd eta(p_max); @@ -211,7 +209,7 @@ class matrix_exp_action_handler { Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(p_max - 1, m_max); for (auto p = 1; p < p_max; ++p) { for (auto i = p * (p + 1) - 2; i < m_max; ++i) { - mt(p - 1, i) = alpha[p - 1] / theta[i]; + mt(p - 1, i) = alpha[p - 1] / theta_m_double[i]; } } From 74df8a54cb6a452b2c2077cd20b2c60c02383f51 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 11:07:23 -0400 Subject: [PATCH 19/28] matrix_exp_action_handler: set_approx_order -- consistent indexing --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index f53769ee921..334110d429e 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -202,7 +202,7 @@ class matrix_exp_action_handler { eta[p] = c; } for (auto p = 0; p < p_max - 1; ++p) { - alpha[p] = std::max(eta[p], eta(p + 1)); + alpha[p] = std::max(eta[p], eta[p + 1]); } } From e40da15d75c96ad4b7ef5f496e57aead52f0fccc Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 11:09:19 -0400 Subject: [PATCH 20/28] matrix_exp_action_handler: set_approx_order -- simplifying code --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index 334110d429e..fa51b805fe6 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -196,10 +196,9 @@ class matrix_exp_action_handler { alpha = Eigen::VectorXd::Constant(p_max - 1, normA); } else { Eigen::VectorXd eta(p_max); - for (auto p = 0; p < p_max; ++p) { - double c = mat_power_1_norm(mat, p + 2); - c = std::pow(c, 1.0 / (p + 2.0)); - eta[p] = c; + for (int p = 0; p < p_max; ++p) { + eta[p] = std::pow(mat_power_1_norm(mat, p + 2), + 1.0 / (p + 2)); } for (auto p = 0; p < p_max - 1; ++p) { alpha[p] = std::max(eta[p], eta[p + 1]); From 35c259b54321b9924ee76bb09d2e7890d166f6cb Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 11:09:28 -0400 Subject: [PATCH 21/28] matrix_exp_action_handler: set_approx_order -- simplifying code --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index fa51b805fe6..48650679742 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -200,7 +200,7 @@ class matrix_exp_action_handler { eta[p] = std::pow(mat_power_1_norm(mat, p + 2), 1.0 / (p + 2)); } - for (auto p = 0; p < p_max - 1; ++p) { + for (int p = 0; p < p_max - 1; ++p) { alpha[p] = std::max(eta[p], eta[p + 1]); } } From f9476a12f3d100096f4c261ad5a296681cc15691 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 11:10:30 -0400 Subject: [PATCH 22/28] matrix_exp_action_handler: set_approx_order -- replacing auto with int --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index 48650679742..693a45d9e2b 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -206,14 +206,14 @@ class matrix_exp_action_handler { } Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(p_max - 1, m_max); - for (auto p = 1; p < p_max; ++p) { - for (auto i = p * (p + 1) - 2; i < m_max; ++i) { + for (int p = 1; p < p_max; ++p) { + for (int i = p * (p + 1) - 2; i < m_max; ++i) { mt(p - 1, i) = alpha[p - 1] / theta_m_double[i]; } } Eigen::Matrix u(m_max); - for (auto i = 0; i < m_max; ++i) { + for (int i = 0; i < m_max; ++i) { u(i) = i + 1; } From df936684a44bd6c75748fbeb965362fa6db5332e Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 11:21:23 -0400 Subject: [PATCH 23/28] matrix_exp_action_handler: set_approx_order -- using LinSpaced to set u --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index 693a45d9e2b..f976a158361 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -212,10 +212,7 @@ class matrix_exp_action_handler { } } - Eigen::Matrix u(m_max); - for (int i = 0; i < m_max; ++i) { - u(i) = i + 1; - } + Eigen::VectorXd u = Eigen::VectorXd::LinSpaced(Eigen::Sequential, m_max, 1, m_max); Eigen::MatrixXd c = stan::math::ceil(mt) * u.asDiagonal(); for (auto i = 0; i < c.size(); ++i) { From 8feefd430aeb0a3fef989337323131f4a11741d1 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 16:51:17 -0400 Subject: [PATCH 24/28] matrix_exp_action_handler: set_approx_order -- returning special case --- .../prim/fun/matrix_exp_action_handler.hpp | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index f976a158361..d2b1541f7b3 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -1,8 +1,6 @@ #ifndef STAN_MATH_PRIM_FUN_MATRIX_EXP_ACTION_HANDLER_HPP #define STAN_MATH_PRIM_FUN_MATRIX_EXP_ACTION_HANDLER_HPP -#include - #include #include #include @@ -167,9 +165,10 @@ class matrix_exp_action_handler { * (paragraph between eq. 3.11 & eq. 3.12). * * @param [in] mat matrix A + * @param [in] b matrix B * @param [in] t double t in exp(A*t)*B - * @param [out] m int parameter m - * @param [out] s int parameter s + * @param [out] m int parameter m; m >= 0, m < m_max + * @param [out] s int parameter s; s >= 1 */ template * = nullptr, @@ -212,18 +211,18 @@ class matrix_exp_action_handler { } } - Eigen::VectorXd u = Eigen::VectorXd::LinSpaced(Eigen::Sequential, m_max, 1, m_max); - + Eigen::VectorXd u = Eigen::VectorXd::LinSpaced(m_max, 1, m_max); Eigen::MatrixXd c = stan::math::ceil(mt) * u.asDiagonal(); - for (auto i = 0; i < c.size(); ++i) { + for (Eigen::MatrixXd::Index i = 0; i < c.size(); ++i) { if (c(i) <= 1.e-16) { c(i) = std::numeric_limits::infinity(); } } - int cost = c.colwise().minCoeff().minCoeff(&m); - if (std::isinf(cost)) - cost = 0; + if (std::isinf(cost)) { + s = 1; + return; + } s = std::max(cost / m, 1); } }; From 42228ef515512b1e1c15c417f087e6f983e263bb Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 19 May 2023 16:52:44 -0400 Subject: [PATCH 25/28] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index d2b1541f7b3..7daaf8da3d4 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -180,7 +180,7 @@ class matrix_exp_action_handler { s = 1; return; } - + // L1 norm double normA = mat.colwise().template lpNorm<1>().maxCoeff(); if (normA < tol) { @@ -191,23 +191,22 @@ class matrix_exp_action_handler { Eigen::VectorXd alpha(p_max - 1); if (normA * (m_max * b.cols()) - < 4.0 * theta_m_double[m_max] * p_max * (p_max + 3)) { + < 4.0 * theta_m_double[m_max] * p_max * (p_max + 3)) { alpha = Eigen::VectorXd::Constant(p_max - 1, normA); } else { Eigen::VectorXd eta(p_max); for (int p = 0; p < p_max; ++p) { - eta[p] = std::pow(mat_power_1_norm(mat, p + 2), - 1.0 / (p + 2)); + eta[p] = std::pow(mat_power_1_norm(mat, p + 2), 1.0 / (p + 2)); } for (int p = 0; p < p_max - 1; ++p) { - alpha[p] = std::max(eta[p], eta[p + 1]); + alpha[p] = std::max(eta[p], eta[p + 1]); } } Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(p_max - 1, m_max); for (int p = 1; p < p_max; ++p) { for (int i = p * (p + 1) - 2; i < m_max; ++i) { - mt(p - 1, i) = alpha[p - 1] / theta_m_double[i]; + mt(p - 1, i) = alpha[p - 1] / theta_m_double[i]; } } @@ -215,14 +214,14 @@ class matrix_exp_action_handler { Eigen::MatrixXd c = stan::math::ceil(mt) * u.asDiagonal(); for (Eigen::MatrixXd::Index i = 0; i < c.size(); ++i) { if (c(i) <= 1.e-16) { - c(i) = std::numeric_limits::infinity(); + c(i) = std::numeric_limits::infinity(); } } int cost = c.colwise().minCoeff().minCoeff(&m); if (std::isinf(cost)) { s = 1; return; - } + } s = std::max(cost / m, 1); } }; From f4b6c42ad14634b7b0b912be96e2a7cd42a10f4e Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 16:54:26 -0400 Subject: [PATCH 26/28] matrix_exp_action_handler: adding one more definition --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index 7daaf8da3d4..5b57b4aa673 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -230,6 +230,7 @@ class matrix_exp_action_handler { // to the declaration. This is fixed in C++-17. constexpr int matrix_exp_action_handler::p_max; constexpr int matrix_exp_action_handler::m_max; +constexpr double matrix_exp_action_handler::tol; } // namespace math } // namespace stan From 8f1c923e1c247eafc3cfc021a6adcb174c69d603 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 19 May 2023 17:26:01 -0400 Subject: [PATCH 27/28] matrix_exp_action_handler: removing auto --- stan/math/prim/fun/matrix_exp_action_handler.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index 5b57b4aa673..e038ba22e2a 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -96,10 +96,10 @@ class matrix_exp_action_handler { Eigen::MatrixXd f = b_eval; Eigen::MatrixXd bi = b_eval; - for (auto i = 0; i < s; ++i) { + for (int i = 0; i < s; ++i) { // maximum absolute row sum, aka inf-norm of matrix operator double c1 = matrix_operator_inf_norm(bi); - for (auto j = 0; j < m; ++j) { + for (int j = 0; j < m; ++j) { bi = (t / (s * (j + 1))) * (A * bi); f += bi; double c2 = matrix_operator_inf_norm(bi); From 3e12a84917043904be63c0a2f38e9ef1d89309aa Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Mon, 22 May 2023 09:45:24 -0400 Subject: [PATCH 28/28] matrix_exp_action_handler: converting static constexpr to const int, renaming for consistency --- .../prim/fun/matrix_exp_action_handler.hpp | 63 ++++++------------- 1 file changed, 20 insertions(+), 43 deletions(-) diff --git a/stan/math/prim/fun/matrix_exp_action_handler.hpp b/stan/math/prim/fun/matrix_exp_action_handler.hpp index e038ba22e2a..92bb8117671 100644 --- a/stan/math/prim/fun/matrix_exp_action_handler.hpp +++ b/stan/math/prim/fun/matrix_exp_action_handler.hpp @@ -23,27 +23,10 @@ namespace math { * and t is double. */ class matrix_exp_action_handler { - static constexpr int p_max = 8; - static constexpr int m_max = 55; - static constexpr double tol = 1.0e-16; - - const std::vector theta_m_single{ - 1.19209280e-07, 5.97885889e-04, 1.12338647e-02, 5.11661936e-02, - 1.30848716e-01, 2.49528932e-01, 4.01458242e-01, 5.80052463e-01, - 7.79511337e-01, 9.95184079e-01, 1.22347954e+00, 1.46166151e+00, - 1.70764853e+00, 1.95985059e+00, 2.21704439e+00, 2.47828088e+00, - 2.74281711e+00, 3.01006636e+00, 3.27956121e+00, 3.55092621e+00, - 3.82385743e+00, 4.09810697e+00, 4.37347131e+00, 4.64978222e+00, - 4.92689984e+00, 5.20470723e+00, 5.48310609e+00, 5.76201341e+00, - 6.04135876e+00, 6.32108213e+00, 6.60113218e+00, 6.88146485e+00, - 7.16204215e+00, 7.44283129e+00, 7.72380381e+00, 8.00493499e+00, - 8.28620327e+00, 8.56758983e+00, 8.84907819e+00, 9.13065388e+00, - 9.41230420e+00, 9.69401796e+00, 9.97578527e+00, 1.02575974e+01, - 1.05394467e+01, 1.08213263e+01, 1.11032302e+01, 1.13851530e+01, - 1.16670899e+01, 1.19490366e+01, 1.22309896e+01, 1.25129453e+01, - 1.27949008e+01, 1.30768536e+01, 1.33588011e+01, 1.36407415e+01, - 1.39226727e+01, 1.42045932e+01, 1.44865015e+01, 1.47683963e+01}; - const std::vector theta_m_double{ + const int _p_max = 8; + const int _m_max = 55; + const double _tol = 1.1e-16; // from the paper, double precision: 2^-53 + const std::vector _theta_m{ 2.22044605e-16, 2.58095680e-08, 1.38634787e-05, 3.39716884e-04, 2.40087636e-03, 9.06565641e-03, 2.38445553e-02, 4.99122887e-02, 8.95776020e-02, 1.44182976e-01, 2.14235807e-01, 2.99615891e-01, @@ -103,7 +86,7 @@ class matrix_exp_action_handler { bi = (t / (s * (j + 1))) * (A * bi); f += bi; double c2 = matrix_operator_inf_norm(bi); - if (c1 + c2 < tol * matrix_operator_inf_norm(f)) + if (c1 + c2 < _tol * matrix_operator_inf_norm(f)) break; c1 = c2; } @@ -167,7 +150,7 @@ class matrix_exp_action_handler { * @param [in] mat matrix A * @param [in] b matrix B * @param [in] t double t in exp(A*t)*B - * @param [out] m int parameter m; m >= 0, m < m_max + * @param [out] m int parameter m; m >= 0, m < _m_max * @param [out] s int parameter s; s >= 1 */ template * = nullptr> inline void set_approx_order(const EigMat1& mat, const EigMat2& b, const double& t, int& m, int& s) { - if (t < tol) { + if (t < _tol) { m = 0; s = 1; return; @@ -183,34 +166,34 @@ class matrix_exp_action_handler { // L1 norm double normA = mat.colwise().template lpNorm<1>().maxCoeff(); - if (normA < tol) { + if (normA < _tol) { m = 0; s = 1; return; } - Eigen::VectorXd alpha(p_max - 1); - if (normA * (m_max * b.cols()) - < 4.0 * theta_m_double[m_max] * p_max * (p_max + 3)) { - alpha = Eigen::VectorXd::Constant(p_max - 1, normA); + Eigen::VectorXd alpha(_p_max - 1); + if (normA * (_m_max * b.cols()) + < 4.0 * _theta_m[_m_max] * _p_max * (_p_max + 3)) { + alpha = Eigen::VectorXd::Constant(_p_max - 1, normA); } else { - Eigen::VectorXd eta(p_max); - for (int p = 0; p < p_max; ++p) { + Eigen::VectorXd eta(_p_max); + for (int p = 0; p < _p_max; ++p) { eta[p] = std::pow(mat_power_1_norm(mat, p + 2), 1.0 / (p + 2)); } - for (int p = 0; p < p_max - 1; ++p) { + for (int p = 0; p < _p_max - 1; ++p) { alpha[p] = std::max(eta[p], eta[p + 1]); } } - Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(p_max - 1, m_max); - for (int p = 1; p < p_max; ++p) { - for (int i = p * (p + 1) - 2; i < m_max; ++i) { - mt(p - 1, i) = alpha[p - 1] / theta_m_double[i]; + Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(_p_max - 1, _m_max); + for (int p = 1; p < _p_max; ++p) { + for (int i = p * (p + 1) - 2; i < _m_max; ++i) { + mt(p - 1, i) = alpha[p - 1] / _theta_m[i]; } } - Eigen::VectorXd u = Eigen::VectorXd::LinSpaced(m_max, 1, m_max); + Eigen::VectorXd u = Eigen::VectorXd::LinSpaced(_m_max, 1, _m_max); Eigen::MatrixXd c = stan::math::ceil(mt) * u.asDiagonal(); for (Eigen::MatrixXd::Index i = 0; i < c.size(); ++i) { if (c(i) <= 1.e-16) { @@ -226,12 +209,6 @@ class matrix_exp_action_handler { } }; -// C++-11 requires definition of the static member in addition -// to the declaration. This is fixed in C++-17. -constexpr int matrix_exp_action_handler::p_max; -constexpr int matrix_exp_action_handler::m_max; -constexpr double matrix_exp_action_handler::tol; - } // namespace math } // namespace stan