From 61b1588d3b1301a1e128ee009f2fa0c937ed4ad7 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Wed, 2 Mar 2022 17:50:06 -0500 Subject: [PATCH 01/18] prim version of fft, inv_fft, with tests --- stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/fft.hpp | 66 ++++++++++++++++++++++++++++ test/unit/math/prim/fun/fft_test.cpp | 64 +++++++++++++++++++++++++++ 3 files changed, 131 insertions(+) create mode 100644 stan/math/prim/fun/fft.hpp create mode 100644 test/unit/math/prim/fun/fft_test.cpp diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 79b5bd2c230..d55c7a286cb 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -99,6 +99,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp new file mode 100644 index 00000000000..8f3bc4855d0 --- /dev/null +++ b/stan/math/prim/fun/fft.hpp @@ -0,0 +1,66 @@ +#ifndef STAN_MATH_PRIM_FUN_FFT_HPP +#define STAN_MATH_PRIM_FUN_FFT_HPP + +#include +#include +#include +#include + +namespace stan { +namespace math { + + /** + * Return the forward discrete Fourier transform of the specified + * complex vector. If the input is size zero, the output will be + * size zero. + * + * Given an input complex vector `x[0:N-1]` of size `N`, the forward + * FFT computes entries of the resulting complex vector `y[0:N-1]` by + * + * ``` + * y[i] = SUM_{j < N} x[j] exp(-2 * pi * i * j * sqrt(-1) / n) + * ``` + * + * @tparam T scalar type of real and imaginary components + * @param x complex time vector to transform + * @return discrete Fourier transform of `x` + */ + template + Eigen::Matrix, -1, 1> + fft(const Eigen::Matrix, -1, 1>& x) { + if (x.size() <= 1) return x; + Eigen::FFT fft; + auto y = fft.fwd(x); + return y; + } + + + /** + * Return the inverse discrete Fourier transform of the specified + * complex vector. If the input is size zero, the output will be + * size zero. + * + * Given an input complex vector `y[0:N-1]` of size `N`, the reverse + * FFT computes entries of the resulting complex vector `x[0:N-1]` by + * + * ``` + * x[i] = SUM_{j < N} y[j] exp(2 * pi * i * j * sqrt(-1) / n) + * ``` + * + * @tparam T scalar type of real and imaginary components + * @param y complex frequency domain vector to inverse transform + * @return inverse discrete Fourier transform of `x` + */ + template + Eigen::Matrix, -1, 1> + inv_fft(const Eigen::Matrix, -1, 1>& y) { + if (y.size() <= 1) return y; + Eigen::FFT fft; + auto x = fft.inv(y); + return x; + } + +} +} + +#endif diff --git a/test/unit/math/prim/fun/fft_test.cpp b/test/unit/math/prim/fun/fft_test.cpp new file mode 100644 index 00000000000..78e78538d25 --- /dev/null +++ b/test/unit/math/prim/fun/fft_test.cpp @@ -0,0 +1,64 @@ +#include +#include +#include + +TEST(primFun, fft) { + typedef std::complex c_t; + typedef Eigen::Matrix, -1, 1> cv_t; + using stan::math::fft; + + // reference answers calculated using Scipy.fft with double precision + + cv_t x0(0); + cv_t y0 = fft(x0); + EXPECT_EQ(0, y0.size()); + + cv_t x1(1); + x1 << c_t(-3.247, 1.98555); + cv_t y1 = fft(x1); + EXPECT_EQ(1, y1.size()); + EXPECT_EQ(real(y1[0]), -3.247); + EXPECT_EQ(imag(y1[0]), 1.98555); + + cv_t x(3); + x << c_t(1, -2), c_t(-3, 5), c_t(-7, 11); + cv_t y = fft(x); + EXPECT_EQ(3, y.size()); + EXPECT_NEAR(real(y[0]), -9, 1e-6); + EXPECT_NEAR(imag(y[0]), 14, 1e-6); + EXPECT_NEAR(real(y[1]), 0.80384758, 1e-6); + EXPECT_NEAR(imag(y[1]), -13.46410162, 1e-6); + EXPECT_NEAR(real(y[2]), 11.19615242, 1e-6); + EXPECT_NEAR(imag(y[2]), -6.53589838, 1e-6); +} + + +TEST(primFun, inv_fft) { + typedef std::complex c_t; + typedef Eigen::Matrix, -1, 1> cv_t; + using stan::math::inv_fft; + + // reference answers calculated using Scipy.fft with double precision + + cv_t y0(0); + cv_t x0 = inv_fft(y0); + EXPECT_EQ(0, x0.size()); + + cv_t y1(1); + y1 << c_t(-3.247, 1.98555); + cv_t x1 = inv_fft(y1); + EXPECT_EQ(1, x1.size()); + EXPECT_EQ(real(x1[0]), -3.247); + EXPECT_EQ(imag(x1[0]), 1.98555); + + cv_t y(3); + y << c_t(-9, 14), c_t(0.80384758, -13.46410162), c_t(11.19615242, -6.53589838); + cv_t x = inv_fft(y); + EXPECT_EQ(3, y.size()); + EXPECT_NEAR(real(x[0]), 1, 1e-6); + EXPECT_NEAR(imag(x[0]), -2, 1e-6); + EXPECT_NEAR(real(x[1]), -3, 1e-6); + EXPECT_NEAR(imag(x[1]), 5, 1e-6); + EXPECT_NEAR(real(x[2]), -7, 1e-6); + EXPECT_NEAR(imag(x[2]), 11, 1e-6); +} From e35b7628317cf771bed0d8c1d9ce3bddf56ae542 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Thu, 3 Mar 2022 13:00:22 -0500 Subject: [PATCH 02/18] fft expression template args --- stan/math/prim/fun/fft.hpp | 26 ++++++++++++++------------ test/unit/math/prim/fun/fft_test.cpp | 16 ++++++++++++++-- 2 files changed, 28 insertions(+), 14 deletions(-) diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index 8f3bc4855d0..976838471d4 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -1,9 +1,11 @@ #ifndef STAN_MATH_PRIM_FUN_FFT_HPP #define STAN_MATH_PRIM_FUN_FFT_HPP +#include #include #include #include +#include #include namespace stan { @@ -25,12 +27,12 @@ namespace math { * @param x complex time vector to transform * @return discrete Fourier transform of `x` */ - template - Eigen::Matrix, -1, 1> - fft(const Eigen::Matrix, -1, 1>& x) { - if (x.size() <= 1) return x; - Eigen::FFT fft; - auto y = fft.fwd(x); + template * = nullptr> + inline Eigen::Matrix, -1, 1> fft(const V& x) { + Eigen::Matrix, -1, 1> xv = x; + if (xv.size() <= 1) return xv; + Eigen::FFT::value_type> fft; + auto y = fft.fwd(xv); return y; } @@ -51,12 +53,12 @@ namespace math { * @param y complex frequency domain vector to inverse transform * @return inverse discrete Fourier transform of `x` */ - template - Eigen::Matrix, -1, 1> - inv_fft(const Eigen::Matrix, -1, 1>& y) { - if (y.size() <= 1) return y; - Eigen::FFT fft; - auto x = fft.inv(y); + template * = nullptr> + inline Eigen::Matrix, -1, 1> inv_fft(const V& y) { + Eigen::Matrix, -1, 1> yv = y; + if (y.size() <= 1) return yv; + Eigen::FFT::value_type> fft; + auto x = fft.inv(yv); return x; } diff --git a/test/unit/math/prim/fun/fft_test.cpp b/test/unit/math/prim/fun/fft_test.cpp index 78e78538d25..24f93479c51 100644 --- a/test/unit/math/prim/fun/fft_test.cpp +++ b/test/unit/math/prim/fun/fft_test.cpp @@ -19,7 +19,11 @@ TEST(primFun, fft) { EXPECT_EQ(1, y1.size()); EXPECT_EQ(real(y1[0]), -3.247); EXPECT_EQ(imag(y1[0]), 1.98555); - + + cv_t y1b = fft(x1 + x1); + EXPECT_NEAR(real(y1b[0]), -3.247 * 2, 1e-6); + EXPECT_NEAR(imag(y1b[0]), 1.98555 * 2, 1e-6); + cv_t x(3); x << c_t(1, -2), c_t(-3, 5), c_t(-7, 11); cv_t y = fft(x); @@ -30,8 +34,16 @@ TEST(primFun, fft) { EXPECT_NEAR(imag(y[1]), -13.46410162, 1e-6); EXPECT_NEAR(real(y[2]), 11.19615242, 1e-6); EXPECT_NEAR(imag(y[2]), -6.53589838, 1e-6); -} + cv_t yb = fft(x + x); + EXPECT_EQ(3, yb.size()); + EXPECT_NEAR(real(yb[0]), 2 * -9, 1e-6); + EXPECT_NEAR(imag(yb[0]), 2 * 14, 1e-6); + EXPECT_NEAR(real(yb[1]), 2 * 0.80384758, 1e-6); + EXPECT_NEAR(imag(yb[1]), 2 * -13.46410162, 1e-6); + EXPECT_NEAR(real(yb[2]), 2 * 11.19615242, 1e-6); + EXPECT_NEAR(imag(yb[2]), 2 * -6.53589838, 1e-6); +} TEST(primFun, inv_fft) { typedef std::complex c_t; From ff0248c5a38194eee0acf90ccb5117326bd2a754 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Thu, 3 Mar 2022 17:51:17 -0500 Subject: [PATCH 03/18] generic type guards for fft and inv_fft --- stan/math/prim/fun/fft.hpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index 976838471d4..b07a9511bfa 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -27,13 +27,13 @@ namespace math { * @param x complex time vector to transform * @return discrete Fourier transform of `x` */ - template * = nullptr> + template * = nullptr> + // version w/o is_complex restriction: template * = nullptr> inline Eigen::Matrix, -1, 1> fft(const V& x) { Eigen::Matrix, -1, 1> xv = x; if (xv.size() <= 1) return xv; Eigen::FFT::value_type> fft; - auto y = fft.fwd(xv); - return y; + return fft.fwd(xv); } @@ -53,16 +53,15 @@ namespace math { * @param y complex frequency domain vector to inverse transform * @return inverse discrete Fourier transform of `x` */ - template * = nullptr> + template * = nullptr> inline Eigen::Matrix, -1, 1> inv_fft(const V& y) { Eigen::Matrix, -1, 1> yv = y; if (y.size() <= 1) return yv; Eigen::FFT::value_type> fft; - auto x = fft.inv(yv); - return x; + return fft.inv(yv); } -} -} +} // namespace math +} // namespace stan #endif From 9b0028dde34c9714864a47ac8a66fc244ec8e0c0 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Mon, 7 Mar 2022 17:27:51 -0500 Subject: [PATCH 04/18] autodiff tests for fft (failing due to finite diffs) --- stan/math/prim/fun/fft.hpp | 27 ++++++----- test/unit/math/mix/fun/fft_test.cpp | 75 +++++++++++++++++++++++++++++ 2 files changed, 91 insertions(+), 11 deletions(-) create mode 100644 test/unit/math/mix/fun/fft_test.cpp diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index b07a9511bfa..85b2b60bdae 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -12,23 +12,25 @@ namespace stan { namespace math { /** - * Return the forward discrete Fourier transform of the specified - * complex vector. If the input is size zero, the output will be - * size zero. + * Return the discrete Fourier transform of the specified complex + * vector. The input vector may be considered to be in the time + * domain and the output will be in the frequency domain. * - * Given an input complex vector `x[0:N-1]` of size `N`, the forward - * FFT computes entries of the resulting complex vector `y[0:N-1]` by + * Given an input complex vector `x[0:N-1]` of size `N`, the discrete + * Fourier transform computes entries of the resulting complex + * vector `y[0:N-1]` by * * ``` * y[i] = SUM_{j < N} x[j] exp(-2 * pi * i * j * sqrt(-1) / n) * ``` + * + * If `y` is of size zero, the result is a size zero vector. * * @tparam T scalar type of real and imaginary components - * @param x complex time vector to transform + * @param x complex time domain vector to transform * @return discrete Fourier transform of `x` */ template * = nullptr> - // version w/o is_complex restriction: template * = nullptr> inline Eigen::Matrix, -1, 1> fft(const V& x) { Eigen::Matrix, -1, 1> xv = x; if (xv.size() <= 1) return xv; @@ -39,16 +41,19 @@ namespace math { /** * Return the inverse discrete Fourier transform of the specified - * complex vector. If the input is size zero, the output will be - * size zero. + * complex vector. The input may be considered to be in the + * frequency domain and the output will be in the time domain. * - * Given an input complex vector `y[0:N-1]` of size `N`, the reverse - * FFT computes entries of the resulting complex vector `x[0:N-1]` by + * Given an input complex vector `y[0:N-1]` of size `N`, the inverse + * discrete Fourier transform computes entries of the resulting + * complex vector `x[0:N-1]` by * * ``` * x[i] = SUM_{j < N} y[j] exp(2 * pi * i * j * sqrt(-1) / n) * ``` * + * If the input is size zero, the output will be size zero. + * * @tparam T scalar type of real and imaginary components * @param y complex frequency domain vector to inverse transform * @return inverse discrete Fourier transform of `x` diff --git a/test/unit/math/mix/fun/fft_test.cpp b/test/unit/math/mix/fun/fft_test.cpp new file mode 100644 index 00000000000..f8b40efc109 --- /dev/null +++ b/test/unit/math/mix/fun/fft_test.cpp @@ -0,0 +1,75 @@ +#include + +template +Eigen::Matrix +flatten(const Eigen::Matrix, -1, 1>& x) { + Eigen::Matrix xf(x.size() / 2); + for (int i = 0; i < x.size(); ++i) { + xf[i / 2] = std::real(xf[i]); + xf[i / 2 + 1] = std::imag(xf[i]); + } + return xf; +} + +template +Eigen::Matrix to_real_vec(const Eigen::Matrix, -1, 1>& x) { + Eigen::Matrix y(x.size() * 2); + for (int i = 0; i < y.size(); i += 2) { + y[i] = real(x[i / 2]); + y[i + 1] = imag(x[i / 2]); + } + return y; +} + +template +Eigen::Matrix, -1, 1> to_complex_vec(const Eigen::Matrix& x) { + Eigen::Matrix, -1, 1> y(x.size() / 2); + for (int i = 0; i < y.size(); ++i) + y[i] = { x[2 * i], x[2 * i + 1] }; + return y; +} + +struct foo { + template + Eigen::Matrix operator()(const Eigen::Matrix& x) const { + auto xc = to_complex_vec(x); + auto y = stan::math::fft(xc); + return to_real_vec(y); + } +}; + + + + +TEST(mathMixFun, fft) { + typedef Eigen::Matrix, -1, 1> cvec_t; + typedef std::complex c_t; + auto f = [](const auto& x) { + using stan::math::fft; + return fft(x); + }; + + cvec_t x0(0); + stan::test::expect_ad(f, x0); + + cvec_t x1(1); + x1[0] = { 1, 2 }; + stan::test::expect_ad(f, x1); + + cvec_t x2(2); + x2[0] = { 1, 2}; + x2[1] = {-1.3, 2.9}; + + // stan::test::expect_ad(f, x2); + + Eigen::VectorXd x(6); + x << 1, -1.3, 2.9, 14.7, -12.9, -4.8; + + Eigen::VectorXd gx; + Eigen::MatrixXd J; + foo g; + stan::math::jacobian(g, x, gx, J); + std::cout << "J = " << std::endl << J << std::endl; + +} + From cf71d6264d30982483ca3139a5c41149d22b54d2 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Mon, 7 Mar 2022 17:36:11 -0500 Subject: [PATCH 05/18] failing fft test example --- test/unit/math/mix/fun/fft_test.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/unit/math/mix/fun/fft_test.cpp b/test/unit/math/mix/fun/fft_test.cpp index f8b40efc109..4b832e6d195 100644 --- a/test/unit/math/mix/fun/fft_test.cpp +++ b/test/unit/math/mix/fun/fft_test.cpp @@ -60,11 +60,12 @@ TEST(mathMixFun, fft) { x2[0] = { 1, 2}; x2[1] = {-1.3, 2.9}; - // stan::test::expect_ad(f, x2); - Eigen::VectorXd x(6); x << 1, -1.3, 2.9, 14.7, -12.9, -4.8; + auto xcv = to_complex_vec(x); + stan::test::expect_ad(f, xcv); + Eigen::VectorXd gx; Eigen::MatrixXd J; foo g; From c843aacb41a4f6b43246606d29ed725c02ac705c Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Tue, 8 Mar 2022 17:57:41 -0500 Subject: [PATCH 06/18] forward 2D FFT with tests, but low precision match to numpy --- stan/math/prim/fun/fft.hpp | 92 ++++++++++++++++++++++++++-- test/unit/math/prim/fun/fft_test.cpp | 80 ++++++++++++++++++++++++ 2 files changed, 168 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index 85b2b60bdae..5b4a6fdda81 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -21,12 +21,12 @@ namespace math { * vector `y[0:N-1]` by * * ``` - * y[i] = SUM_{j < N} x[j] exp(-2 * pi * i * j * sqrt(-1) / n) + * y[n] = SUM_{i < N} x[i] * exp(n * i * -2 * pi * sqrt(-1) / N) * ``` * * If `y` is of size zero, the result is a size zero vector. * - * @tparam T scalar type of real and imaginary components + * @tparam V type of complex vector argument * @param x complex time domain vector to transform * @return discrete Fourier transform of `x` */ @@ -49,12 +49,15 @@ namespace math { * complex vector `x[0:N-1]` by * * ``` - * x[i] = SUM_{j < N} y[j] exp(2 * pi * i * j * sqrt(-1) / n) + * x[n] = SUM_{i < N} y[i] * exp(n * i * 2 * pi * sqrt(-1) / N) * ``` * + * The only difference between the discrete DFT and its inverse is + * the sign of the exponent. + * * If the input is size zero, the output will be size zero. * - * @tparam T scalar type of real and imaginary components + * @tparam V type of complex vector argument * @param y complex frequency domain vector to inverse transform * @return inverse discrete Fourier transform of `x` */ @@ -65,6 +68,87 @@ namespace math { Eigen::FFT::value_type> fft; return fft.inv(yv); } + + /** + * Return the two-dimensional discrete Fourier transform to the + * specified complex matrix. Given a complex matrix (M x N) matrix + * x, the definition of the 2D discrete Fourier transform is + * + * ``` + * y[m, n] = SUM_{i < M, j < N} x[i, j] + * * exp(i * m * -2 * pi * sqrt(-1) / N) + * * exp(j * n * -2 * pi * sqrt(-1) / M) + * ``` + * + * Another way to view the 2D FFT is as first running a 1D FFT on + * each row, then running a 1D FFT on each resulting column. + * + * @tparam M type of complex matrix argument + * @param x complex time-domain matrix + * @return discrete 2D Fourier transform of `x` + */ + template // , require_eigen_matrix_dynamic_t* = nullptr> + inline Eigen::Matrix, -1, -1> fft2(const M& x) { + if (x.size() <= 1) return x; + int rows = x.rows(); + int cols = x.cols(); + Eigen::FFT::value_type> fft; + Eigen::Matrix, -1, -1> y(rows, cols); + Eigen::Matrix, -1, 1> v_temp; + if (cols > 1) { + for (int i = 0; i < rows; ++i) { + // implicit transposition by assignment + Eigen::Matrix, -1, 1> row_i_t = x.row(i); + fft.fwd(v_temp, row_i_t); + y.row(i) = v_temp; + } + } + if (rows > 1) { + for (int j = 0; j < cols; ++j) { + Eigen::Matrix, -1, 1> col_j = y.col(j); + fft.fwd(v_temp, col_j); + y.col(j) = v_temp; + } + } + return y; + } + + // /** + // * Return the two-dimensional discrete Fourier transform to the + // * specified complex matrix. Given a complex matrix (M x N) matrix + // * x, the definition of the 2D discrete Fourier transform is + // * + // * ``` + // * x[m, n] = SUM_{i < M, j < N} y[i, j] + // * * exp(i * m * 2 * pi * sqrt(-1) / N) + // * * exp(j * n * 2 * pi * sqrt(-1) / M) + // * ``` + // * + // * Another way to view the 2D inverse FFT is as first running a 1D + // * inverse FFT on each row, then running a 1D invese FFT on each + // * resulting column. + // * + // * @tparam M type of complex matrix argument + // * @param y complex frequency-domain matrix + // * @return inverse discrete 2D Fourier transform of `x` + // */ + // template * = nullptr> + // inline Eigen::Matrix, -1, -1> inv_fft2(const M& y) { + // int rows = y.rows(); + // int cols = y.cols(); + // Eigen::FFT::value_type> fft; + // Eigen::Matrix, -1, -1> x; + // Eigen::Matrix, -1, 1> v_temp; + // for (i = 0; i < rows; ++i) { + // fft.inv(v_temp, y.row(i)); + // x.row(i) = v_temp; + // } + // for (j = 0; j < cols; ++j) { + // fft.inv(v_temp, x.col(j)); + // x.col(j) = v_temp; + // } + // return x; + // } } // namespace math } // namespace stan diff --git a/test/unit/math/prim/fun/fft_test.cpp b/test/unit/math/prim/fun/fft_test.cpp index 24f93479c51..78d8f90f303 100644 --- a/test/unit/math/prim/fun/fft_test.cpp +++ b/test/unit/math/prim/fun/fft_test.cpp @@ -74,3 +74,83 @@ TEST(primFun, inv_fft) { EXPECT_NEAR(real(x[2]), -7, 1e-6); EXPECT_NEAR(imag(x[2]), 11, 1e-6); } + +TEST(primFun, fft2) { + typedef std::complex c_t; + typedef Eigen::Matrix, -1, -1> cm_t; + using stan::math::fft2; + + // reference answers calculated using Scipy.fft with double precision + + cm_t x(0, 0); + cm_t y = fft2(x); + EXPECT_EQ(0, y.rows()); + EXPECT_EQ(0, y.cols()); + + cm_t x11(1, 1); + x11 << c_t(1.0, -3.9); + cm_t y11 = fft2(x11); + EXPECT_NEAR(1.0, std::real(y11(0, 0)), 1e-6); + EXPECT_NEAR(-3.9, std::imag(y11(0, 0)), 1e-6); + + cm_t x12(1, 2); + x12 << c_t(1.0, -3.9), c_t(-8.6, 0.2); + cm_t y12 = fft2(x12); + EXPECT_EQ(1, y12.rows()); + EXPECT_EQ(2, y12.cols()); + EXPECT_NEAR(-7.6, real(y12(0, 0)), 1e-6); + EXPECT_NEAR(-3.7, imag(y12(0, 0)), 1e-6); + EXPECT_NEAR(9.6, real(y12(0, 1)), 1e-6); + EXPECT_NEAR(-4.1, imag(y12(0, 1)), 1e-6); + + cm_t x33(3, 3); + x33 << + c_t(1, 2), c_t(3, -1.4), c_t(2, 1), + c_t(3, -9), c_t(2, 1.3), c_t(3.9, -1.8), + c_t(13, -1.8), c_t(1.3, 1.9), c_t(-2.2, -2.2); + cm_t y33 = fft2(x33); + EXPECT_EQ(3, y33.rows()); + EXPECT_EQ(3, y33.cols()); + + // no idea why this is only two decimal places of precision + // relative to numpy.fft.fft2() + EXPECT_NEAR(27, real(y33(0, 0)), 1e-1); + EXPECT_NEAR(-10.1, imag(y33(0, 0)), 1e-1); + + EXPECT_NEAR(16.0703194, real(y33(0, 1)), 1e-1); + EXPECT_NEAR(-10.40166605, imag(y33(0, 1)), 1e-1); + + EXPECT_NEAR(7.9296806, real(y33(0, 2)), 1e-1); + EXPECT_NEAR(-5.89833395, imag(y33(0, 2)), 1e-1); + + EXPECT_NEAR(-10.90858799, real(y33(1, 0)), 1e-1); + EXPECT_NEAR(10.07128129, imag(y33(1, 0)), 2e-1); + + EXPECT_NEAR(-15.63153533, real(y33(1, 1)), 1e-1); + EXPECT_NEAR(19.63153533, imag(y33(1, 1)), 1e-1); + + EXPECT_NEAR(-13.1660254, real(y33(1, 2)), 1e-1); + EXPECT_NEAR(+18.47794549, imag(y33(1, 2)), 1e-1); + + EXPECT_NEAR(1.90858799, real(y33(2, 0)), 1e-1); + EXPECT_NEAR(4.52871871, imag(y33(2, 0)), 1e-1); + + EXPECT_NEAR(-11.4339746, real(y33(2, 1)), 1e-1); + EXPECT_NEAR(-5.07794549, imag(y33(2, 1)), 1e-1); + + EXPECT_NEAR(7.23153533, real(y33(2, 2)), 1e-1); + EXPECT_NEAR(-3.23153533, imag(y33(2, 2)), 1e-1); + + // this is return from numpy + // >>> a + // array([[ 1. +2.j , 3. -1.5j, 2. +1.j ], + // [ 3. -9.j , 2. +1.3j, 3.9-1.8j], + // [13. -1.8j, 1.3+1.9j, -2.2-2.2j]]) + // >>> np.fft.fft2(a) + // array([[ 27. -10.1j , 16.0703194 -10.40166605j, + // 7.9296806 -5.89833395j], + // [-10.90858799+10.07128129j, -15.63153533+19.63153533j, + // -13.1660254 +18.47794549j], + // [ 1.90858799 +4.52871871j, -11.4339746 -5.07794549j, + // 7.23153533 -3.23153533j]]) +} From 81d190e5e5e15ac9488a5c2c64595f05fe223230 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Mon, 14 Mar 2022 18:12:07 -0400 Subject: [PATCH 07/18] 2D FFT and inverse passing tests --- stan/math/prim/fun/fft.hpp | 102 +++++++++------------------ test/unit/math/prim/fun/fft_test.cpp | 86 +++++++++++++--------- 2 files changed, 83 insertions(+), 105 deletions(-) diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index 5b4a6fdda81..0cdde64ed9a 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -70,85 +71,46 @@ namespace math { } /** - * Return the two-dimensional discrete Fourier transform to the - * specified complex matrix. Given a complex matrix (M x N) matrix - * x, the definition of the 2D discrete Fourier transform is - * - * ``` - * y[m, n] = SUM_{i < M, j < N} x[i, j] - * * exp(i * m * -2 * pi * sqrt(-1) / N) - * * exp(j * n * -2 * pi * sqrt(-1) / M) - * ``` - * - * Another way to view the 2D FFT is as first running a 1D FFT on - * each row, then running a 1D FFT on each resulting column. + * Return the two-dimensional discrete Fourier transform of the + * specified complex matrix. The 2D discrete Fourier transform + * first runs the DFT on the each row, then on each column of the + * result. * * @tparam M type of complex matrix argument * @param x complex time-domain matrix * @return discrete 2D Fourier transform of `x` */ - template // , require_eigen_matrix_dynamic_t* = nullptr> + template * = nullptr> inline Eigen::Matrix, -1, -1> fft2(const M& x) { - if (x.size() <= 1) return x; - int rows = x.rows(); - int cols = x.cols(); - Eigen::FFT::value_type> fft; - Eigen::Matrix, -1, -1> y(rows, cols); - Eigen::Matrix, -1, 1> v_temp; - if (cols > 1) { - for (int i = 0; i < rows; ++i) { - // implicit transposition by assignment - Eigen::Matrix, -1, 1> row_i_t = x.row(i); - fft.fwd(v_temp, row_i_t); - y.row(i) = v_temp; - } - } - if (rows > 1) { - for (int j = 0; j < cols; ++j) { - Eigen::Matrix, -1, 1> col_j = y.col(j); - fft.fwd(v_temp, col_j); - y.col(j) = v_temp; - } - } + Eigen::Matrix, -1, -1> y(x.rows(), x.cols()); + for (int i = 0; i < y.rows(); ++i) + y.row(i) = fft(x.row(i)); + for (int j = 0; j < y.cols(); ++j) + y.col(j) = fft(y.col(j).eval()); + return y; + } + + /** + * Return the two-dimensional inverse discrete Fourier transform of + * the specified complex matrix. The 2D inverse discrete Fourier + * transform first runs the 1D inverse Fourier transform on the + * columns, and then on the resulting rows. The composition of the + * FFT and inverse FFT (or vice-versa) is the identity. + * + * @tparam M type of complex matrix argument + * @param x complex frequency-domain matrix + * @return inverse discrete 2D Fourier transform of `x` + */ + template * = nullptr> + inline Eigen::Matrix, -1, -1> inv_fft2(const M& y) { + Eigen::Matrix, -1, -1> x(y.rows(), y.cols()); + for (int j = 0; j < x.cols(); ++j) + x.col(j) = inv_fft(y.col(j)); + for (int i = 0; i < x.rows(); ++i) + x.row(i) = inv_fft(x.row(i)); return y; } - // /** - // * Return the two-dimensional discrete Fourier transform to the - // * specified complex matrix. Given a complex matrix (M x N) matrix - // * x, the definition of the 2D discrete Fourier transform is - // * - // * ``` - // * x[m, n] = SUM_{i < M, j < N} y[i, j] - // * * exp(i * m * 2 * pi * sqrt(-1) / N) - // * * exp(j * n * 2 * pi * sqrt(-1) / M) - // * ``` - // * - // * Another way to view the 2D inverse FFT is as first running a 1D - // * inverse FFT on each row, then running a 1D invese FFT on each - // * resulting column. - // * - // * @tparam M type of complex matrix argument - // * @param y complex frequency-domain matrix - // * @return inverse discrete 2D Fourier transform of `x` - // */ - // template * = nullptr> - // inline Eigen::Matrix, -1, -1> inv_fft2(const M& y) { - // int rows = y.rows(); - // int cols = y.cols(); - // Eigen::FFT::value_type> fft; - // Eigen::Matrix, -1, -1> x; - // Eigen::Matrix, -1, 1> v_temp; - // for (i = 0; i < rows; ++i) { - // fft.inv(v_temp, y.row(i)); - // x.row(i) = v_temp; - // } - // for (j = 0; j < cols; ++j) { - // fft.inv(v_temp, x.col(j)); - // x.col(j) = v_temp; - // } - // return x; - // } } // namespace math } // namespace stan diff --git a/test/unit/math/prim/fun/fft_test.cpp b/test/unit/math/prim/fun/fft_test.cpp index 78d8f90f303..26046a45f3d 100644 --- a/test/unit/math/prim/fun/fft_test.cpp +++ b/test/unit/math/prim/fun/fft_test.cpp @@ -79,6 +79,7 @@ TEST(primFun, fft2) { typedef std::complex c_t; typedef Eigen::Matrix, -1, -1> cm_t; using stan::math::fft2; + using stan::math::inv_fft2; // reference answers calculated using Scipy.fft with double precision @@ -87,61 +88,64 @@ TEST(primFun, fft2) { EXPECT_EQ(0, y.rows()); EXPECT_EQ(0, y.cols()); - cm_t x11(1, 1); - x11 << c_t(1.0, -3.9); - cm_t y11 = fft2(x11); - EXPECT_NEAR(1.0, std::real(y11(0, 0)), 1e-6); - EXPECT_NEAR(-3.9, std::imag(y11(0, 0)), 1e-6); - - cm_t x12(1, 2); - x12 << c_t(1.0, -3.9), c_t(-8.6, 0.2); - cm_t y12 = fft2(x12); - EXPECT_EQ(1, y12.rows()); - EXPECT_EQ(2, y12.cols()); - EXPECT_NEAR(-7.6, real(y12(0, 0)), 1e-6); - EXPECT_NEAR(-3.7, imag(y12(0, 0)), 1e-6); - EXPECT_NEAR(9.6, real(y12(0, 1)), 1e-6); - EXPECT_NEAR(-4.1, imag(y12(0, 1)), 1e-6); + // cm_t x11(1, 1); + // x11 << c_t(1.0, -3.9); + // cm_t y11 = fft2(x11); + // EXPECT_NEAR(1.0, std::real(y11(0, 0)), 1e-6); + // EXPECT_NEAR(-3.9, std::imag(y11(0, 0)), 1e-6); + + // cm_t x12(1, 2); + // x12 << c_t(1.0, -3.9), c_t(-8.6, 0.2); + // cm_t y12 = fft2(x12); + // EXPECT_EQ(1, y12.rows()); + // EXPECT_EQ(2, y12.cols()); + // EXPECT_NEAR(-7.6, real(y12(0, 0)), 1e-6); + // EXPECT_NEAR(-3.7, imag(y12(0, 0)), 1e-6); + // EXPECT_NEAR(9.6, real(y12(0, 1)), 1e-6); + // EXPECT_NEAR(-4.1, imag(y12(0, 1)), 1e-6); cm_t x33(3, 3); x33 << c_t(1, 2), c_t(3, -1.4), c_t(2, 1), - c_t(3, -9), c_t(2, 1.3), c_t(3.9, -1.8), + c_t(3, -9), c_t(2, -1.3), c_t(3.9, -1.8), c_t(13, -1.8), c_t(1.3, 1.9), c_t(-2.2, -2.2); cm_t y33 = fft2(x33); EXPECT_EQ(3, y33.rows()); EXPECT_EQ(3, y33.cols()); + // cm_t xcopy = inv_fft2(y33); + // std::cout << "round trip=" << std::endl << xcopy << std::endl; + // no idea why this is only two decimal places of precision // relative to numpy.fft.fft2() - EXPECT_NEAR(27, real(y33(0, 0)), 1e-1); - EXPECT_NEAR(-10.1, imag(y33(0, 0)), 1e-1); + EXPECT_NEAR(27, real(y33(0, 0)), 1e-6); + EXPECT_NEAR(-12.6, imag(y33(0, 0)), 1e-6); - EXPECT_NEAR(16.0703194, real(y33(0, 1)), 1e-1); - EXPECT_NEAR(-10.40166605, imag(y33(0, 1)), 1e-1); + EXPECT_NEAR(13.90525589, real(y33(0, 1)), 1e-6); + EXPECT_NEAR(-9.15166605, imag(y33(0, 1)), 1e-6); - EXPECT_NEAR(7.9296806, real(y33(0, 2)), 1e-1); - EXPECT_NEAR(-5.89833395, imag(y33(0, 2)), 1e-1); + EXPECT_NEAR(10.09474411, real(y33(0, 2)), 1e-1); + EXPECT_NEAR(-4.64833395, imag(y33(0, 2)), 2e-1); - EXPECT_NEAR(-10.90858799, real(y33(1, 0)), 1e-1); - EXPECT_NEAR(10.07128129, imag(y33(1, 0)), 2e-1); + EXPECT_NEAR(-13.160254038, real(y33(1, 0)), 1e-1); + EXPECT_NEAR(11.471281292, imag(y33(1, 0)), 2e-1); - EXPECT_NEAR(-15.63153533, real(y33(1, 1)), 1e-1); - EXPECT_NEAR(19.63153533, imag(y33(1, 1)), 1e-1); + EXPECT_NEAR(-13.29326674, real(y33(1, 1)), 1e-1); + EXPECT_NEAR(20.88153533, imag(y33(1, 1)), 1e-1); - EXPECT_NEAR(-13.1660254, real(y33(1, 2)), 1e-1); - EXPECT_NEAR(+18.47794549, imag(y33(1, 2)), 1e-1); + EXPECT_NEAR(-13.25262794, real(y33(1, 2)), 1e-1); + EXPECT_NEAR(15.82794549, imag(y33(1, 2)), 1e-1); - EXPECT_NEAR(1.90858799, real(y33(2, 0)), 1e-1); - EXPECT_NEAR(4.52871871, imag(y33(2, 0)), 1e-1); + EXPECT_NEAR(4.160254038, real(y33(2, 0)), 1e-1); + EXPECT_NEAR(5.928718708, imag(y33(2, 0)), 1e-1); - EXPECT_NEAR(-11.4339746, real(y33(2, 1)), 1e-1); - EXPECT_NEAR(-5.07794549, imag(y33(2, 1)), 1e-1); + EXPECT_NEAR(-11.34737206, real(y33(2, 1)), 1e-1); + EXPECT_NEAR(-7.72794549, imag(y33(2, 1)), 1e-1); - EXPECT_NEAR(7.23153533, real(y33(2, 2)), 1e-1); - EXPECT_NEAR(-3.23153533, imag(y33(2, 2)), 1e-1); + EXPECT_NEAR(4.89326674, real(y33(2, 2)), 1e-1); + EXPECT_NEAR(-1.98153533, imag(y33(2, 2)), 1e-1); - // this is return from numpy + // Python (Scipy) // >>> a // array([[ 1. +2.j , 3. -1.5j, 2. +1.j ], // [ 3. -9.j , 2. +1.3j, 3.9-1.8j], @@ -153,4 +157,16 @@ TEST(primFun, fft2) { // -13.1660254 +18.47794549j], // [ 1.90858799 +4.52871871j, -11.4339746 -5.07794549j, // 7.23153533 -3.23153533j]]) + + // R (stats) + // > b + // [,1] [,2] [,3] + // [1,] 1+2.0i 3.0-1.5i 2.0+1.0i + // [2,] 3-9.0i 2.0+1.3i 3.9-1.8i + // [3,] 13-1.8i 1.3+1.9i -2.2-2.2i + // > fft(b) + // [,1] [,2] [,3] + // [1,] 27.000000-10.100000i 16.07032-10.40167i 7.929681- 5.898334i + // [2,] -10.908588+10.071281i -15.63154+19.63154i -13.166025+18.477945i + // [3,] 1.908588+ 4.528719i -11.43397- 5.07795i 7.231535- 3.231535i } From e5e04f73ad2e8ca496fd929e7099161aa66cd998 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Tue, 15 Mar 2022 12:21:43 -0400 Subject: [PATCH 08/18] fwd fft tested --- test/unit/math/mix/fun/fft_test.cpp | 43 ++++++++++++++++++----------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/test/unit/math/mix/fun/fft_test.cpp b/test/unit/math/mix/fun/fft_test.cpp index 4b832e6d195..8cea181f961 100644 --- a/test/unit/math/mix/fun/fft_test.cpp +++ b/test/unit/math/mix/fun/fft_test.cpp @@ -44,33 +44,44 @@ struct foo { TEST(mathMixFun, fft) { typedef Eigen::Matrix, -1, 1> cvec_t; typedef std::complex c_t; - auto f = [](const auto& x) { + int n = 0; int i = 0; + auto f = [&](const auto& x) { using stan::math::fft; - return fft(x); + return i ? real(fft(x)(n)) : imag(fft(x)(n)); }; - + cvec_t x0(0); - stan::test::expect_ad(f, x0); + EXPECT_EQ(0, stan::math::fft(x0).size()); // no AD to test cvec_t x1(1); x1[0] = { 1, 2 }; - stan::test::expect_ad(f, x1); + for (i = 0; i <= 1; ++i) + for (n = 0; n <= 0; ++n) + stan::test::expect_ad(f, x1); cvec_t x2(2); - x2[0] = { 1, 2}; + x2[0] = {1, 2}; x2[1] = {-1.3, 2.9}; - - Eigen::VectorXd x(6); - x << 1, -1.3, 2.9, 14.7, -12.9, -4.8; + for (i = 0; i <= 1; ++i) + for (n = 0; n <= 1; ++n) + stan::test::expect_ad(f, x2); - auto xcv = to_complex_vec(x); - stan::test::expect_ad(f, xcv); + Eigen::VectorXcd x3(3); + x3[0] = {1, -1.3}; + x3[1] = {2.9, 14.7}; + x3[2] = {-12.9, -4.8}; + for (i = 0; i <= 1; ++i) + for (n = 0; n <= 2; ++n) + stan::test::expect_ad(f, x3); - Eigen::VectorXd gx; - Eigen::MatrixXd J; - foo g; - stan::math::jacobian(g, x, gx, J); - std::cout << "J = " << std::endl << J << std::endl; + Eigen::VectorXcd x4(4); + x4[0] = {1, -1.3}; + x4[1] = {-2.9, 14.7}; + x4[2] = {-12.9, -4.8}; + x4[3] = {8.398, 9.387}; + for (i = 0; i <= 1; ++i) + for (n = 0; n <= 3; ++n) + stan::test::expect_ad(f, x4); } From b345c5141ca598ac41acc2a32604e577b9e73650 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Tue, 15 Mar 2022 15:02:43 -0400 Subject: [PATCH 09/18] fft funs passing tests, fixes #20 --- stan/math/prim/fun/fft.hpp | 12 +- test/unit/math/mix/fun/fft_test.cpp | 210 +++++++++++++++++++-------- test/unit/math/prim/fun/fft_test.cpp | 145 +++++++++++------- 3 files changed, 250 insertions(+), 117 deletions(-) diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index 0cdde64ed9a..030035961b4 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -22,10 +22,10 @@ namespace math { * vector `y[0:N-1]` by * * ``` - * y[n] = SUM_{i < N} x[i] * exp(n * i * -2 * pi * sqrt(-1) / N) + * y[n] = SUM_{i < N} x[i] * exp(-n * i * 2 * pi * sqrt(-1) / N) * ``` * - * If `y` is of size zero, the result is a size zero vector. + * If the input is of size zero, the result is a size zero vector. * * @tparam V type of complex vector argument * @param x complex time domain vector to transform @@ -53,11 +53,10 @@ namespace math { * x[n] = SUM_{i < N} y[i] * exp(n * i * 2 * pi * sqrt(-1) / N) * ``` * + * If the input is of size zero, the result is a size zero vector. * The only difference between the discrete DFT and its inverse is * the sign of the exponent. * - * If the input is size zero, the output will be size zero. - * * @tparam V type of complex vector argument * @param y complex frequency domain vector to inverse transform * @return inverse discrete Fourier transform of `x` @@ -86,7 +85,7 @@ namespace math { for (int i = 0; i < y.rows(); ++i) y.row(i) = fft(x.row(i)); for (int j = 0; j < y.cols(); ++j) - y.col(j) = fft(y.col(j).eval()); + y.col(j) = fft(y.col(j)); return y; } @@ -108,9 +107,8 @@ namespace math { x.col(j) = inv_fft(y.col(j)); for (int i = 0; i < x.rows(); ++i) x.row(i) = inv_fft(x.row(i)); - return y; + return x; } - } // namespace math } // namespace stan diff --git a/test/unit/math/mix/fun/fft_test.cpp b/test/unit/math/mix/fun/fft_test.cpp index 8cea181f961..b5aeddeee34 100644 --- a/test/unit/math/mix/fun/fft_test.cpp +++ b/test/unit/math/mix/fun/fft_test.cpp @@ -1,87 +1,183 @@ #include -template -Eigen::Matrix -flatten(const Eigen::Matrix, -1, 1>& x) { - Eigen::Matrix xf(x.size() / 2); - for (int i = 0; i < x.size(); ++i) { - xf[i / 2] = std::real(xf[i]); - xf[i / 2 + 1] = std::imag(xf[i]); - } - return xf; +void expect_fft(const Eigen::VectorXcd& x) { + int i = 0; + int m = 0; + auto g = [&](const auto& x) { + using stan::math::fft; + return i ? real(fft(x)(m)) : imag(fft(x)(m)); + }; + for (int i = 0; i < 2; ++i) + for (int m = 0; m < x.rows(); ++m) + stan::test::expect_ad(g, x); } -template -Eigen::Matrix to_real_vec(const Eigen::Matrix, -1, 1>& x) { - Eigen::Matrix y(x.size() * 2); - for (int i = 0; i < y.size(); i += 2) { - y[i] = real(x[i / 2]); - y[i + 1] = imag(x[i / 2]); - } - return y; -} +TEST(mathMixFun, fft) { + typedef Eigen::VectorXcd cvec_t; -template -Eigen::Matrix, -1, 1> to_complex_vec(const Eigen::Matrix& x) { - Eigen::Matrix, -1, 1> y(x.size() / 2); - for (int i = 0; i < y.size(); ++i) - y[i] = { x[2 * i], x[2 * i + 1] }; - return y; -} + cvec_t x0(0); + expect_fft(x0); -struct foo { - template - Eigen::Matrix operator()(const Eigen::Matrix& x) const { - auto xc = to_complex_vec(x); - auto y = stan::math::fft(xc); - return to_real_vec(y); - } -}; + cvec_t x1(1); + x1[0] = {1, 2}; + expect_fft(x1); + cvec_t x2(2); + x2[0] = {1, 2}; + x2[1] = {-1.3, 2.9}; + expect_fft(x2); + Eigen::VectorXcd x3(3); + x3[0] = {1, -1.3}; + x3[1] = {2.9, 14.7}; + x3[2] = {-12.9, -4.8}; + expect_fft(x3); + Eigen::VectorXcd x4(4); + x4[0] = {1, -1.3}; + x4[1] = {-2.9, 14.7}; + x4[2] = {-12.9, -4.8}; + x4[3] = {8.398, 9.387}; + expect_fft(x4); +} -TEST(mathMixFun, fft) { - typedef Eigen::Matrix, -1, 1> cvec_t; - typedef std::complex c_t; - int n = 0; int i = 0; - auto f = [&](const auto& x) { - using stan::math::fft; - return i ? real(fft(x)(n)) : imag(fft(x)(n)); +void expect_inv_fft(const Eigen::VectorXcd& x) { + int i = 0; + int m = 0; + auto g = [&](const auto& x) { + using stan::math::inv_fft; + return i ? real(inv_fft(x)(m)) : imag(inv_fft(x)(m)); }; + for (int i = 0; i < 2; ++i) + for (int m = 0; m < x.rows(); ++m) + stan::test::expect_ad(g, x); +} + + +TEST(mathMixFun, invFft) { + typedef Eigen::VectorXcd cvec_t; cvec_t x0(0); - EXPECT_EQ(0, stan::math::fft(x0).size()); // no AD to test + expect_inv_fft(x0); cvec_t x1(1); - x1[0] = { 1, 2 }; - for (i = 0; i <= 1; ++i) - for (n = 0; n <= 0; ++n) - stan::test::expect_ad(f, x1); + x1[0] = {1, 2}; + expect_inv_fft(x1); cvec_t x2(2); x2[0] = {1, 2}; x2[1] = {-1.3, 2.9}; - for (i = 0; i <= 1; ++i) - for (n = 0; n <= 1; ++n) - stan::test::expect_ad(f, x2); - + expect_inv_fft(x2); + Eigen::VectorXcd x3(3); x3[0] = {1, -1.3}; x3[1] = {2.9, 14.7}; x3[2] = {-12.9, -4.8}; - for (i = 0; i <= 1; ++i) - for (n = 0; n <= 2; ++n) - stan::test::expect_ad(f, x3); - + expect_inv_fft(x3); + Eigen::VectorXcd x4(4); x4[0] = {1, -1.3}; x4[1] = {-2.9, 14.7}; x4[2] = {-12.9, -4.8}; x4[3] = {8.398, 9.387}; - for (i = 0; i <= 1; ++i) - for (n = 0; n <= 3; ++n) - stan::test::expect_ad(f, x4); - + expect_inv_fft(x4); +} + +void expect_fft2(const Eigen::MatrixXcd& x) { + int i = 0; + int m = 0; + int n = 0; + auto g = [&](const auto& x) { + using stan::math::fft2; + return i ? real(fft2(x)(m, n)) : imag(fft2(x)(m, n)); + }; + for (int i = 0; i < 2; ++i) + for (int n = 0; n < x.cols(); ++n) + for (int m = 0; m < x.rows(); ++m) + stan::test::expect_ad(g, x); +} + +TEST(mathMixFun, fft2) { + typedef Eigen::MatrixXcd cmat_t; + + cmat_t x00(0, 0); + expect_fft2(x00); + + cmat_t x11(1, 1); + x11(0, 0) = { 1, 2 }; + expect_fft2(x11); + + cmat_t x21(2, 1); + x21(0, 0) = {1, 2}; + x21(1, 0) = {-1.3, 2.9}; + expect_fft2(x21); + + cmat_t x22(2, 2); + x22(0, 0) = {3, 9}; + x22(0, 1) = {-13.2, 8.345}; + x22(1, 0) = {-4.23, 7.566}; + x22(1, 1) = {1, -12.9}; + expect_fft2(x22); + + cmat_t x33(3, 3); + x33(0, 0) = {3, 9}; + x33(0, 1) = {-13.2, 8.345}; + x33(0, 2) = {4.333, -1.9}; + x33(1, 0) = {-4.23, 7.566}; + x33(1, 1) = {1, -12.9}; + x33(1, 2) = {-1.01, -4.01}; + x33(2, 0) = {3.87, 5.89}; + x33(2, 1) = {-2.875, -2.999}; + x33(2, 2) = {12.98, 14.5555}; + expect_fft2(x33); +} + +void expect_inv_fft2(const Eigen::MatrixXcd& x) { + int i = 0; + int m = 0; + int n = 0; + auto g = [&](const auto& x) { + using stan::math::inv_fft2; + return i ? real(inv_fft2(x)(m, n)) : imag(inv_fft2(x)(m, n)); + }; + for (int i = 0; i < 2; ++i) + for (int n = 0; n < x.cols(); ++n) + for (int m = 0; m < x.rows(); ++m) + stan::test::expect_ad(g, x); } +TEST(mathMixFun, inv_fft2) { + typedef Eigen::MatrixXcd cmat_t; + + cmat_t x00(0, 0); + expect_inv_fft2(x00); + + cmat_t x11(1, 1); + x11(0, 0) = { 1, 2 }; + expect_inv_fft2(x11); + + cmat_t x21(2, 1); + x21(0, 0) = {1, 2}; + x21(1, 0) = {-1.3, 2.9}; + expect_inv_fft2(x21); + + cmat_t x22(2, 2); + x22(0, 0) = {3, 9}; + x22(0, 1) = {-13.2, 8.345}; + x22(1, 0) = {-4.23, 7.566}; + x22(1, 1) = {1, -12.9}; + expect_inv_fft2(x22); + + cmat_t x33(3, 3); + x33(0, 0) = {3, 9}; + x33(0, 1) = {-13.2, 8.345}; + x33(0, 2) = {4.333, -1.9}; + x33(1, 0) = {-4.23, 7.566}; + x33(1, 1) = {1, -12.9}; + x33(1, 2) = {-1.01, -4.01}; + x33(2, 0) = {3.87, 5.89}; + x33(2, 1) = {-2.875, -2.999}; + x33(2, 2) = {12.98, 14.5555}; + expect_inv_fft2(x33); +} + diff --git a/test/unit/math/prim/fun/fft_test.cpp b/test/unit/math/prim/fun/fft_test.cpp index 26046a45f3d..d6a015f8f0b 100644 --- a/test/unit/math/prim/fun/fft_test.cpp +++ b/test/unit/math/prim/fun/fft_test.cpp @@ -81,28 +81,29 @@ TEST(primFun, fft2) { using stan::math::fft2; using stan::math::inv_fft2; - // reference answers calculated using Scipy.fft with double precision + // reference answers calculated using Scipy.fft with double + // precision, and verified with R cm_t x(0, 0); cm_t y = fft2(x); EXPECT_EQ(0, y.rows()); EXPECT_EQ(0, y.cols()); - // cm_t x11(1, 1); - // x11 << c_t(1.0, -3.9); - // cm_t y11 = fft2(x11); - // EXPECT_NEAR(1.0, std::real(y11(0, 0)), 1e-6); - // EXPECT_NEAR(-3.9, std::imag(y11(0, 0)), 1e-6); - - // cm_t x12(1, 2); - // x12 << c_t(1.0, -3.9), c_t(-8.6, 0.2); - // cm_t y12 = fft2(x12); - // EXPECT_EQ(1, y12.rows()); - // EXPECT_EQ(2, y12.cols()); - // EXPECT_NEAR(-7.6, real(y12(0, 0)), 1e-6); - // EXPECT_NEAR(-3.7, imag(y12(0, 0)), 1e-6); - // EXPECT_NEAR(9.6, real(y12(0, 1)), 1e-6); - // EXPECT_NEAR(-4.1, imag(y12(0, 1)), 1e-6); + cm_t x11(1, 1); + x11 << c_t(1.0, -3.9); + cm_t y11 = fft2(x11); + EXPECT_NEAR(1.0, std::real(y11(0, 0)), 1e-6); + EXPECT_NEAR(-3.9, std::imag(y11(0, 0)), 1e-6); + + cm_t x12(1, 2); + x12 << c_t(1.0, -3.9), c_t(-8.6, 0.2); + cm_t y12 = fft2(x12); + EXPECT_EQ(1, y12.rows()); + EXPECT_EQ(2, y12.cols()); + EXPECT_NEAR(-7.6, real(y12(0, 0)), 1e-6); + EXPECT_NEAR(-3.7, imag(y12(0, 0)), 1e-6); + EXPECT_NEAR(9.6, real(y12(0, 1)), 1e-6); + EXPECT_NEAR(-4.1, imag(y12(0, 1)), 1e-6); cm_t x33(3, 3); x33 << @@ -113,60 +114,98 @@ TEST(primFun, fft2) { EXPECT_EQ(3, y33.rows()); EXPECT_EQ(3, y33.cols()); - // cm_t xcopy = inv_fft2(y33); - // std::cout << "round trip=" << std::endl << xcopy << std::endl; - - // no idea why this is only two decimal places of precision - // relative to numpy.fft.fft2() + // check vs. results from R (matches SciPy) EXPECT_NEAR(27, real(y33(0, 0)), 1e-6); EXPECT_NEAR(-12.6, imag(y33(0, 0)), 1e-6); - EXPECT_NEAR(13.90525589, real(y33(0, 1)), 1e-6); EXPECT_NEAR(-9.15166605, imag(y33(0, 1)), 1e-6); - EXPECT_NEAR(10.09474411, real(y33(0, 2)), 1e-1); EXPECT_NEAR(-4.64833395, imag(y33(0, 2)), 2e-1); - EXPECT_NEAR(-13.160254038, real(y33(1, 0)), 1e-1); EXPECT_NEAR(11.471281292, imag(y33(1, 0)), 2e-1); - EXPECT_NEAR(-13.29326674, real(y33(1, 1)), 1e-1); EXPECT_NEAR(20.88153533, imag(y33(1, 1)), 1e-1); - EXPECT_NEAR(-13.25262794, real(y33(1, 2)), 1e-1); EXPECT_NEAR(15.82794549, imag(y33(1, 2)), 1e-1); - EXPECT_NEAR(4.160254038, real(y33(2, 0)), 1e-1); EXPECT_NEAR(5.928718708, imag(y33(2, 0)), 1e-1); - EXPECT_NEAR(-11.34737206, real(y33(2, 1)), 1e-1); EXPECT_NEAR(-7.72794549, imag(y33(2, 1)), 1e-1); - EXPECT_NEAR(4.89326674, real(y33(2, 2)), 1e-1); EXPECT_NEAR(-1.98153533, imag(y33(2, 2)), 1e-1); +} + +void expect_complex_matrix_float_eq( + const Eigen::Matrix, -1, -1>& x, + const Eigen::Matrix, -1, -1>& y) { + EXPECT_EQ(x.rows(), y.rows()); + EXPECT_EQ(x.cols(), y.cols()); + for (int j = 0; j < y.cols(); ++j) { + for (int i = 0; i < y.rows(); ++i) { + EXPECT_FLOAT_EQ(real(x(i, j)), real(y(i, j))); + EXPECT_FLOAT_EQ(imag(x(i, j)), imag(y(i, j))); + } + } +} + +TEST(primFunFFT, invfft2) { + typedef std::complex c_t; + typedef Eigen::Matrix, -1, -1> cm_t; + using stan::math::fft2; + using stan::math::inv_fft2; + using stan::math::inv_fft; + + // reference answers calculated using R + cm_t x(0, 0); + cm_t y = inv_fft2(x); + EXPECT_EQ(0, y.rows()); + EXPECT_EQ(0, y.cols()); + + cm_t x11(1, 1); + x11 << c_t(1.0, -3.9); + cm_t y11 = inv_fft2(x11); + EXPECT_NEAR(1.0, std::real(y11(0, 0)), 1e-6); + EXPECT_NEAR(-3.9, std::imag(y11(0, 0)), 1e-6); + + cm_t x13(1, 3); + x13 << c_t(-2.3, 1.82), c_t(1.18, 9.32), c_t(1.15, -14.1); + cm_t y13 = inv_fft2(x13); + cm_t y13copy = inv_fft(x13.row(0)); + expect_complex_matrix_float_eq(y13, y13copy.transpose()); + + cm_t x33(3, 3); + x33 << + c_t(1, 2), c_t(3, -1.4), c_t(2, 1), + c_t(3, -9), c_t(2, -1.3), c_t(3.9, -1.8), + c_t(13, -1.8), c_t(1.3, 1.9), c_t(-2.2, -2.2); + cm_t y33 = fft2(x33); - // Python (Scipy) - // >>> a - // array([[ 1. +2.j , 3. -1.5j, 2. +1.j ], - // [ 3. -9.j , 2. +1.3j, 3.9-1.8j], - // [13. -1.8j, 1.3+1.9j, -2.2-2.2j]]) - // >>> np.fft.fft2(a) - // array([[ 27. -10.1j , 16.0703194 -10.40166605j, - // 7.9296806 -5.89833395j], - // [-10.90858799+10.07128129j, -15.63153533+19.63153533j, - // -13.1660254 +18.47794549j], - // [ 1.90858799 +4.52871871j, -11.4339746 -5.07794549j, - // 7.23153533 -3.23153533j]]) - - // R (stats) - // > b - // [,1] [,2] [,3] - // [1,] 1+2.0i 3.0-1.5i 2.0+1.0i - // [2,] 3-9.0i 2.0+1.3i 3.9-1.8i - // [3,] 13-1.8i 1.3+1.9i -2.2-2.2i - // > fft(b) - // [,1] [,2] [,3] - // [1,] 27.000000-10.100000i 16.07032-10.40167i 7.929681- 5.898334i - // [2,] -10.908588+10.071281i -15.63154+19.63154i -13.166025+18.477945i - // [3,] 1.908588+ 4.528719i -11.43397- 5.07795i 7.231535- 3.231535i + // check versus results from R + EXPECT_NEAR(27.000000, real(y33(0, 0)), 1e-5); + EXPECT_NEAR(-12.600000, imag(y33(0, 0)), 1e-5); + EXPECT_NEAR(13.90526, real(y33(0, 1)), 1e-5); + EXPECT_NEAR(-9.15167, imag(y33(0, 1)), 1e-5); + EXPECT_NEAR(10.094744, real(y33(0, 2)), 1e-5); + EXPECT_NEAR(-4.648334, imag(y33(0, 2)), 1e-5); + EXPECT_NEAR(-13.160254, real(y33(1, 0)), 1e-5); + EXPECT_NEAR(11.471281, imag(y33(1, 0)), 1e-5); + EXPECT_NEAR(-13.29327, real(y33(1, 1)), 1e-5); + EXPECT_NEAR(20.88154, imag(y33(1, 1)), 1e-5); + EXPECT_NEAR(-13.252628, real(y33(1, 2)), 1e-5); + EXPECT_NEAR(15.827945, imag(y33(1, 2)), 1e-5); + EXPECT_NEAR(4.160254, real(y33(2, 0)), 1e-5); + EXPECT_NEAR(5.928719, imag(y33(2, 0)), 1e-5); + EXPECT_NEAR(-11.34737, real(y33(2, 1)), 1e-5); + EXPECT_NEAR(-7.72795, imag(y33(2, 1)), 1e-5); + EXPECT_NEAR(4.893267, real(y33(2, 2)), 1e-5); + EXPECT_NEAR(-1.981535, imag(y33(2, 2)), 1e-5); + + // check round trips inv_fft(fft(x)) + cm_t x33copy = inv_fft2(y33); + expect_complex_matrix_float_eq(x33, x33copy); + + // check round trip fft(inv_fft(x)) + cm_t z33 = inv_fft2(x33); + cm_t x33copy2 = fft2(z33); + expect_complex_matrix_float_eq(x33, x33copy2); } From f023e83cfd5d763415aa8e8b91a71fbc261fca6e Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Tue, 15 Mar 2022 15:43:04 -0400 Subject: [PATCH 10/18] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/fft.hpp | 189 ++++++++++++++------------- test/unit/math/mix/fun/fft_test.cpp | 18 ++- test/unit/math/prim/fun/fft_test.cpp | 27 ++-- 3 files changed, 115 insertions(+), 119 deletions(-) diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index 030035961b4..032dcf4782f 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -12,104 +12,105 @@ namespace stan { namespace math { - /** - * Return the discrete Fourier transform of the specified complex - * vector. The input vector may be considered to be in the time - * domain and the output will be in the frequency domain. - * - * Given an input complex vector `x[0:N-1]` of size `N`, the discrete - * Fourier transform computes entries of the resulting complex - * vector `y[0:N-1]` by - * - * ``` - * y[n] = SUM_{i < N} x[i] * exp(-n * i * 2 * pi * sqrt(-1) / N) - * ``` - * - * If the input is of size zero, the result is a size zero vector. - * - * @tparam V type of complex vector argument - * @param x complex time domain vector to transform - * @return discrete Fourier transform of `x` - */ - template * = nullptr> - inline Eigen::Matrix, -1, 1> fft(const V& x) { - Eigen::Matrix, -1, 1> xv = x; - if (xv.size() <= 1) return xv; - Eigen::FFT::value_type> fft; - return fft.fwd(xv); - } +/** + * Return the discrete Fourier transform of the specified complex + * vector. The input vector may be considered to be in the time + * domain and the output will be in the frequency domain. + * + * Given an input complex vector `x[0:N-1]` of size `N`, the discrete + * Fourier transform computes entries of the resulting complex + * vector `y[0:N-1]` by + * + * ``` + * y[n] = SUM_{i < N} x[i] * exp(-n * i * 2 * pi * sqrt(-1) / N) + * ``` + * + * If the input is of size zero, the result is a size zero vector. + * + * @tparam V type of complex vector argument + * @param x complex time domain vector to transform + * @return discrete Fourier transform of `x` + */ +template * = nullptr> +inline Eigen::Matrix, -1, 1> fft(const V& x) { + Eigen::Matrix, -1, 1> xv = x; + if (xv.size() <= 1) + return xv; + Eigen::FFT::value_type> fft; + return fft.fwd(xv); +} +/** + * Return the inverse discrete Fourier transform of the specified + * complex vector. The input may be considered to be in the + * frequency domain and the output will be in the time domain. + * + * Given an input complex vector `y[0:N-1]` of size `N`, the inverse + * discrete Fourier transform computes entries of the resulting + * complex vector `x[0:N-1]` by + * + * ``` + * x[n] = SUM_{i < N} y[i] * exp(n * i * 2 * pi * sqrt(-1) / N) + * ``` + * + * If the input is of size zero, the result is a size zero vector. + * The only difference between the discrete DFT and its inverse is + * the sign of the exponent. + * + * @tparam V type of complex vector argument + * @param y complex frequency domain vector to inverse transform + * @return inverse discrete Fourier transform of `x` + */ +template * = nullptr> +inline Eigen::Matrix, -1, 1> inv_fft(const V& y) { + Eigen::Matrix, -1, 1> yv = y; + if (y.size() <= 1) + return yv; + Eigen::FFT::value_type> fft; + return fft.inv(yv); +} - /** - * Return the inverse discrete Fourier transform of the specified - * complex vector. The input may be considered to be in the - * frequency domain and the output will be in the time domain. - * - * Given an input complex vector `y[0:N-1]` of size `N`, the inverse - * discrete Fourier transform computes entries of the resulting - * complex vector `x[0:N-1]` by - * - * ``` - * x[n] = SUM_{i < N} y[i] * exp(n * i * 2 * pi * sqrt(-1) / N) - * ``` - * - * If the input is of size zero, the result is a size zero vector. - * The only difference between the discrete DFT and its inverse is - * the sign of the exponent. - * - * @tparam V type of complex vector argument - * @param y complex frequency domain vector to inverse transform - * @return inverse discrete Fourier transform of `x` - */ - template * = nullptr> - inline Eigen::Matrix, -1, 1> inv_fft(const V& y) { - Eigen::Matrix, -1, 1> yv = y; - if (y.size() <= 1) return yv; - Eigen::FFT::value_type> fft; - return fft.inv(yv); - } +/** + * Return the two-dimensional discrete Fourier transform of the + * specified complex matrix. The 2D discrete Fourier transform + * first runs the DFT on the each row, then on each column of the + * result. + * + * @tparam M type of complex matrix argument + * @param x complex time-domain matrix + * @return discrete 2D Fourier transform of `x` + */ +template * = nullptr> +inline Eigen::Matrix, -1, -1> fft2(const M& x) { + Eigen::Matrix, -1, -1> y(x.rows(), x.cols()); + for (int i = 0; i < y.rows(); ++i) + y.row(i) = fft(x.row(i)); + for (int j = 0; j < y.cols(); ++j) + y.col(j) = fft(y.col(j)); + return y; +} - /** - * Return the two-dimensional discrete Fourier transform of the - * specified complex matrix. The 2D discrete Fourier transform - * first runs the DFT on the each row, then on each column of the - * result. - * - * @tparam M type of complex matrix argument - * @param x complex time-domain matrix - * @return discrete 2D Fourier transform of `x` - */ - template * = nullptr> - inline Eigen::Matrix, -1, -1> fft2(const M& x) { - Eigen::Matrix, -1, -1> y(x.rows(), x.cols()); - for (int i = 0; i < y.rows(); ++i) - y.row(i) = fft(x.row(i)); - for (int j = 0; j < y.cols(); ++j) - y.col(j) = fft(y.col(j)); - return y; - } +/** + * Return the two-dimensional inverse discrete Fourier transform of + * the specified complex matrix. The 2D inverse discrete Fourier + * transform first runs the 1D inverse Fourier transform on the + * columns, and then on the resulting rows. The composition of the + * FFT and inverse FFT (or vice-versa) is the identity. + * + * @tparam M type of complex matrix argument + * @param x complex frequency-domain matrix + * @return inverse discrete 2D Fourier transform of `x` + */ +template * = nullptr> +inline Eigen::Matrix, -1, -1> inv_fft2(const M& y) { + Eigen::Matrix, -1, -1> x(y.rows(), y.cols()); + for (int j = 0; j < x.cols(); ++j) + x.col(j) = inv_fft(y.col(j)); + for (int i = 0; i < x.rows(); ++i) + x.row(i) = inv_fft(x.row(i)); + return x; +} - /** - * Return the two-dimensional inverse discrete Fourier transform of - * the specified complex matrix. The 2D inverse discrete Fourier - * transform first runs the 1D inverse Fourier transform on the - * columns, and then on the resulting rows. The composition of the - * FFT and inverse FFT (or vice-versa) is the identity. - * - * @tparam M type of complex matrix argument - * @param x complex frequency-domain matrix - * @return inverse discrete 2D Fourier transform of `x` - */ - template * = nullptr> - inline Eigen::Matrix, -1, -1> inv_fft2(const M& y) { - Eigen::Matrix, -1, -1> x(y.rows(), y.cols()); - for (int j = 0; j < x.cols(); ++j) - x.col(j) = inv_fft(y.col(j)); - for (int i = 0; i < x.rows(); ++i) - x.row(i) = inv_fft(x.row(i)); - return x; - } - } // namespace math } // namespace stan diff --git a/test/unit/math/mix/fun/fft_test.cpp b/test/unit/math/mix/fun/fft_test.cpp index b5aeddeee34..6279ea4cbca 100644 --- a/test/unit/math/mix/fun/fft_test.cpp +++ b/test/unit/math/mix/fun/fft_test.cpp @@ -53,7 +53,6 @@ void expect_inv_fft(const Eigen::VectorXcd& x) { stan::test::expect_ad(g, x); } - TEST(mathMixFun, invFft) { typedef Eigen::VectorXcd cvec_t; @@ -68,13 +67,13 @@ TEST(mathMixFun, invFft) { x2[0] = {1, 2}; x2[1] = {-1.3, 2.9}; expect_inv_fft(x2); - + Eigen::VectorXcd x3(3); x3[0] = {1, -1.3}; x3[1] = {2.9, 14.7}; x3[2] = {-12.9, -4.8}; expect_inv_fft(x3); - + Eigen::VectorXcd x4(4); x4[0] = {1, -1.3}; x4[1] = {-2.9, 14.7}; @@ -94,7 +93,7 @@ void expect_fft2(const Eigen::MatrixXcd& x) { for (int i = 0; i < 2; ++i) for (int n = 0; n < x.cols(); ++n) for (int m = 0; m < x.rows(); ++m) - stan::test::expect_ad(g, x); + stan::test::expect_ad(g, x); } TEST(mathMixFun, fft2) { @@ -104,7 +103,7 @@ TEST(mathMixFun, fft2) { expect_fft2(x00); cmat_t x11(1, 1); - x11(0, 0) = { 1, 2 }; + x11(0, 0) = {1, 2}; expect_fft2(x11); cmat_t x21(2, 1); @@ -130,7 +129,7 @@ TEST(mathMixFun, fft2) { x33(2, 1) = {-2.875, -2.999}; x33(2, 2) = {12.98, 14.5555}; expect_fft2(x33); -} +} void expect_inv_fft2(const Eigen::MatrixXcd& x) { int i = 0; @@ -143,7 +142,7 @@ void expect_inv_fft2(const Eigen::MatrixXcd& x) { for (int i = 0; i < 2; ++i) for (int n = 0; n < x.cols(); ++n) for (int m = 0; m < x.rows(); ++m) - stan::test::expect_ad(g, x); + stan::test::expect_ad(g, x); } TEST(mathMixFun, inv_fft2) { @@ -153,7 +152,7 @@ TEST(mathMixFun, inv_fft2) { expect_inv_fft2(x00); cmat_t x11(1, 1); - x11(0, 0) = { 1, 2 }; + x11(0, 0) = {1, 2}; expect_inv_fft2(x11); cmat_t x21(2, 1); @@ -179,5 +178,4 @@ TEST(mathMixFun, inv_fft2) { x33(2, 1) = {-2.875, -2.999}; x33(2, 2) = {12.98, 14.5555}; expect_inv_fft2(x33); -} - +} diff --git a/test/unit/math/prim/fun/fft_test.cpp b/test/unit/math/prim/fun/fft_test.cpp index d6a015f8f0b..988d5ce6da7 100644 --- a/test/unit/math/prim/fun/fft_test.cpp +++ b/test/unit/math/prim/fun/fft_test.cpp @@ -25,7 +25,7 @@ TEST(primFun, fft) { EXPECT_NEAR(imag(y1b[0]), 1.98555 * 2, 1e-6); cv_t x(3); - x << c_t(1, -2), c_t(-3, 5), c_t(-7, 11); + x << c_t(1, -2), c_t(-3, 5), c_t(-7, 11); cv_t y = fft(x); EXPECT_EQ(3, y.size()); EXPECT_NEAR(real(y[0]), -9, 1e-6); @@ -62,9 +62,10 @@ TEST(primFun, inv_fft) { EXPECT_EQ(1, x1.size()); EXPECT_EQ(real(x1[0]), -3.247); EXPECT_EQ(imag(x1[0]), 1.98555); - + cv_t y(3); - y << c_t(-9, 14), c_t(0.80384758, -13.46410162), c_t(11.19615242, -6.53589838); + y << c_t(-9, 14), c_t(0.80384758, -13.46410162), + c_t(11.19615242, -6.53589838); cv_t x = inv_fft(y); EXPECT_EQ(3, y.size()); EXPECT_NEAR(real(x[0]), 1, 1e-6); @@ -106,10 +107,8 @@ TEST(primFun, fft2) { EXPECT_NEAR(-4.1, imag(y12(0, 1)), 1e-6); cm_t x33(3, 3); - x33 << - c_t(1, 2), c_t(3, -1.4), c_t(2, 1), - c_t(3, -9), c_t(2, -1.3), c_t(3.9, -1.8), - c_t(13, -1.8), c_t(1.3, 1.9), c_t(-2.2, -2.2); + x33 << c_t(1, 2), c_t(3, -1.4), c_t(2, 1), c_t(3, -9), c_t(2, -1.3), + c_t(3.9, -1.8), c_t(13, -1.8), c_t(1.3, 1.9), c_t(-2.2, -2.2); cm_t y33 = fft2(x33); EXPECT_EQ(3, y33.rows()); EXPECT_EQ(3, y33.cols()); @@ -136,8 +135,8 @@ TEST(primFun, fft2) { } void expect_complex_matrix_float_eq( - const Eigen::Matrix, -1, -1>& x, - const Eigen::Matrix, -1, -1>& y) { + const Eigen::Matrix, -1, -1>& x, + const Eigen::Matrix, -1, -1>& y) { EXPECT_EQ(x.rows(), y.rows()); EXPECT_EQ(x.cols(), y.cols()); for (int j = 0; j < y.cols(); ++j) { @@ -152,8 +151,8 @@ TEST(primFunFFT, invfft2) { typedef std::complex c_t; typedef Eigen::Matrix, -1, -1> cm_t; using stan::math::fft2; - using stan::math::inv_fft2; using stan::math::inv_fft; + using stan::math::inv_fft2; // reference answers calculated using R cm_t x(0, 0); @@ -172,12 +171,10 @@ TEST(primFunFFT, invfft2) { cm_t y13 = inv_fft2(x13); cm_t y13copy = inv_fft(x13.row(0)); expect_complex_matrix_float_eq(y13, y13copy.transpose()); - + cm_t x33(3, 3); - x33 << - c_t(1, 2), c_t(3, -1.4), c_t(2, 1), - c_t(3, -9), c_t(2, -1.3), c_t(3.9, -1.8), - c_t(13, -1.8), c_t(1.3, 1.9), c_t(-2.2, -2.2); + x33 << c_t(1, 2), c_t(3, -1.4), c_t(2, 1), c_t(3, -9), c_t(2, -1.3), + c_t(3.9, -1.8), c_t(13, -1.8), c_t(1.3, 1.9), c_t(-2.2, -2.2); cm_t y33 = fft2(x33); // check versus results from R From 4b45aba7e1bc0a06be0ad0f4c9ca6f57c023a5e2 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Wed, 16 Mar 2022 09:55:15 -0400 Subject: [PATCH 11/18] fft doc --- stan/math/prim/fun/fft.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index 030035961b4..814fdd55cd4 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -97,7 +97,7 @@ namespace math { * FFT and inverse FFT (or vice-versa) is the identity. * * @tparam M type of complex matrix argument - * @param x complex frequency-domain matrix + * @param y complex frequency-domain matrix * @return inverse discrete 2D Fourier transform of `x` */ template * = nullptr> From 1954e9677dcdbbf39efe3ba92723d1748856b991 Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Wed, 16 Mar 2022 11:21:05 -0400 Subject: [PATCH 12/18] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/fft.hpp | 42 +++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index 40e7ccfa7ed..76f4f8ad8d8 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -90,27 +90,27 @@ inline Eigen::Matrix, -1, -1> fft2(const M& x) { return y; } - /** - * Return the two-dimensional inverse discrete Fourier transform of - * the specified complex matrix. The 2D inverse discrete Fourier - * transform first runs the 1D inverse Fourier transform on the - * columns, and then on the resulting rows. The composition of the - * FFT and inverse FFT (or vice-versa) is the identity. - * - * @tparam M type of complex matrix argument - * @param y complex frequency-domain matrix - * @return inverse discrete 2D Fourier transform of `x` - */ - template * = nullptr> - inline Eigen::Matrix, -1, -1> inv_fft2(const M& y) { - Eigen::Matrix, -1, -1> x(y.rows(), y.cols()); - for (int j = 0; j < x.cols(); ++j) - x.col(j) = inv_fft(y.col(j)); - for (int i = 0; i < x.rows(); ++i) - x.row(i) = inv_fft(x.row(i)); - return x; - } - +/** + * Return the two-dimensional inverse discrete Fourier transform of + * the specified complex matrix. The 2D inverse discrete Fourier + * transform first runs the 1D inverse Fourier transform on the + * columns, and then on the resulting rows. The composition of the + * FFT and inverse FFT (or vice-versa) is the identity. + * + * @tparam M type of complex matrix argument + * @param y complex frequency-domain matrix + * @return inverse discrete 2D Fourier transform of `x` + */ +template * = nullptr> +inline Eigen::Matrix, -1, -1> inv_fft2(const M& y) { + Eigen::Matrix, -1, -1> x(y.rows(), y.cols()); + for (int j = 0; j < x.cols(); ++j) + x.col(j) = inv_fft(y.col(j)); + for (int i = 0; i < x.rows(); ++i) + x.row(i) = inv_fft(x.row(i)); + return x; +} + } // namespace math } // namespace stan From 6b3e712a5e59d0f67dd0f78e6bba068b33c7660b Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Tue, 12 Apr 2022 16:41:02 -0400 Subject: [PATCH 13/18] elide test if tolerance is inf; fft test rev only --- test/unit/math/ad_tolerances.hpp | 31 +++++++++++++++++++++++++++ test/unit/math/mix/fun/fft_test.cpp | 12 +++++++---- test/unit/math/relative_tolerance.hpp | 10 +++++++++ test/unit/math/test_ad.hpp | 29 ++++++++++++++++++++----- 4 files changed, 73 insertions(+), 9 deletions(-) diff --git a/test/unit/math/ad_tolerances.hpp b/test/unit/math/ad_tolerances.hpp index 9d8e84e4b35..a67bc5803dd 100644 --- a/test/unit/math/ad_tolerances.hpp +++ b/test/unit/math/ad_tolerances.hpp @@ -67,8 +67,39 @@ struct ad_tolerances { grad_hessian_val_(1e-8), grad_hessian_hessian_(1e-3), grad_hessian_grad_hessian_(1e-2) {} + }; +/** + * Return tolerances that are infinite other than for the value + * tolerance and gradient tolerance for reverse mode. + * + * @param val_tol value relative tolerance (default `1e-8`) + * @param grad_tol gradient relative tolerance (default `1e-4`) + */ +ad_tolerances reverse_only_ad_tolerances(double val_tol = 1e-8, + double grad_tol = 1e-4) { + ad_tolerances tols; + tols.gradient_val_ = val_tol; + tols.gradient_grad_ = grad_tol; + double inf = std::numeric_limits::infinity(); + tols.gradient_grad_varmat_matvar_ = inf; + tols.gradient_fvar_val_ = inf; + tols.gradient_fvar_grad_ = inf; + tols.hessian_val_ = inf; + tols.hessian_grad_ = inf; + tols.hessian_hessian_ = inf; + tols.hessian_fvar_val_ = inf; + tols.hessian_fvar_grad_ = inf; + tols.hessian_fvar_hessian_ = inf; + tols.grad_hessian_val_ = inf; + tols.grad_hessian_hessian_ = inf; + tols.grad_hessian_grad_hessian_ = inf; + return tols; +} + + + } // namespace test } // namespace stan #endif diff --git a/test/unit/math/mix/fun/fft_test.cpp b/test/unit/math/mix/fun/fft_test.cpp index 6279ea4cbca..dd145fa5908 100644 --- a/test/unit/math/mix/fun/fft_test.cpp +++ b/test/unit/math/mix/fun/fft_test.cpp @@ -1,6 +1,7 @@ #include void expect_fft(const Eigen::VectorXcd& x) { + stan::test::ad_tolerances tols = stan::test::reverse_only_ad_tolerances(); int i = 0; int m = 0; auto g = [&](const auto& x) { @@ -9,7 +10,7 @@ void expect_fft(const Eigen::VectorXcd& x) { }; for (int i = 0; i < 2; ++i) for (int m = 0; m < x.rows(); ++m) - stan::test::expect_ad(g, x); + stan::test::expect_ad(tols, g, x); } TEST(mathMixFun, fft) { @@ -42,6 +43,7 @@ TEST(mathMixFun, fft) { } void expect_inv_fft(const Eigen::VectorXcd& x) { + stan::test::ad_tolerances tols = stan::test::reverse_only_ad_tolerances(); int i = 0; int m = 0; auto g = [&](const auto& x) { @@ -50,7 +52,7 @@ void expect_inv_fft(const Eigen::VectorXcd& x) { }; for (int i = 0; i < 2; ++i) for (int m = 0; m < x.rows(); ++m) - stan::test::expect_ad(g, x); + stan::test::expect_ad(tols, g, x); } TEST(mathMixFun, invFft) { @@ -83,6 +85,7 @@ TEST(mathMixFun, invFft) { } void expect_fft2(const Eigen::MatrixXcd& x) { + stan::test::ad_tolerances tols = stan::test::reverse_only_ad_tolerances(); int i = 0; int m = 0; int n = 0; @@ -93,7 +96,7 @@ void expect_fft2(const Eigen::MatrixXcd& x) { for (int i = 0; i < 2; ++i) for (int n = 0; n < x.cols(); ++n) for (int m = 0; m < x.rows(); ++m) - stan::test::expect_ad(g, x); + stan::test::expect_ad(tols, g, x); } TEST(mathMixFun, fft2) { @@ -132,6 +135,7 @@ TEST(mathMixFun, fft2) { } void expect_inv_fft2(const Eigen::MatrixXcd& x) { + stan::test::ad_tolerances tols = stan::test::reverse_only_ad_tolerances(); int i = 0; int m = 0; int n = 0; @@ -142,7 +146,7 @@ void expect_inv_fft2(const Eigen::MatrixXcd& x) { for (int i = 0; i < 2; ++i) for (int n = 0; n < x.cols(); ++n) for (int m = 0; m < x.rows(); ++m) - stan::test::expect_ad(g, x); + stan::test::expect_ad(tols, g, x); } TEST(mathMixFun, inv_fft2) { diff --git a/test/unit/math/relative_tolerance.hpp b/test/unit/math/relative_tolerance.hpp index d6ea28b2293..e9384836cc2 100644 --- a/test/unit/math/relative_tolerance.hpp +++ b/test/unit/math/relative_tolerance.hpp @@ -3,6 +3,7 @@ #include #include +#include namespace stan { namespace test { @@ -79,6 +80,15 @@ class relative_tolerance { return std::max(tol_ * 0.5 * (fabs(x) + fabs(y)), tol_min_); } + /** + * Returns `true` if the tolerance is infinite. + * + * @return `true` if the tolernace is infinite. + */ + bool is_inf() const { + return std::isinf(tol_); + } + private: /** * The relative tolerance diff --git a/test/unit/math/test_ad.hpp b/test/unit/math/test_ad.hpp index c5bf2732a6b..23ae8f3a6c3 100644 --- a/test/unit/math/test_ad.hpp +++ b/test/unit/math/test_ad.hpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -333,12 +334,30 @@ template void expect_ad_derivatives(const ad_tolerances& tols, const G& g, const Eigen::VectorXd& x) { double gx = g(x); - test_gradient(tols, g, x, gx); + if (!tols.gradient_val_.is_inf() + || !tols.gradient_grad_.is_inf()) { + test_gradient(tols, g, x, gx); + } #ifndef STAN_MATH_TESTS_REV_ONLY - test_gradient_fvar(tols, g, x, gx); - test_hessian(tols, g, x, gx); - test_hessian_fvar(tols, g, x, gx); - test_grad_hessian(tols, g, x, gx); + if (!tols.gradient_fvar_val_.is_inf() + || !tols.gradient_fvar_grad_.is_inf()) { + test_gradient_fvar(tols, g, x, gx); + } + if (!tols.hessian_val_.is_inf() + || !tols.hessian_grad_.is_inf() + || !tols.hessian_hessian_.is_inf()) { + test_hessian(tols, g, x, gx); + } + if (!tols.hessian_fvar_val_.is_inf() + || !tols.hessian_fvar_grad_.is_inf() + || !tols.hessian_hessian_.is_inf()) { + test_hessian_fvar(tols, g, x, gx); + } + if (!tols.grad_hessian_val_.is_inf() + || !tols.grad_hessian_hessian_.is_inf() + || !tols.grad_hessian_grad_hessian_.is_inf()) { + test_grad_hessian(tols, g, x, gx); + } #endif } From cb269222d1a107af74b22b7f6cb7d7b70dc0b069 Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Tue, 12 Apr 2022 16:42:01 -0400 Subject: [PATCH 14/18] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- test/unit/math/ad_tolerances.hpp | 7 ++----- test/unit/math/relative_tolerance.hpp | 4 +--- test/unit/math/test_ad.hpp | 15 +++++---------- 3 files changed, 8 insertions(+), 18 deletions(-) diff --git a/test/unit/math/ad_tolerances.hpp b/test/unit/math/ad_tolerances.hpp index a67bc5803dd..7465bc44b71 100644 --- a/test/unit/math/ad_tolerances.hpp +++ b/test/unit/math/ad_tolerances.hpp @@ -67,18 +67,17 @@ struct ad_tolerances { grad_hessian_val_(1e-8), grad_hessian_hessian_(1e-3), grad_hessian_grad_hessian_(1e-2) {} - }; /** * Return tolerances that are infinite other than for the value - * tolerance and gradient tolerance for reverse mode. + * tolerance and gradient tolerance for reverse mode. * * @param val_tol value relative tolerance (default `1e-8`) * @param grad_tol gradient relative tolerance (default `1e-4`) */ ad_tolerances reverse_only_ad_tolerances(double val_tol = 1e-8, - double grad_tol = 1e-4) { + double grad_tol = 1e-4) { ad_tolerances tols; tols.gradient_val_ = val_tol; tols.gradient_grad_ = grad_tol; @@ -98,8 +97,6 @@ ad_tolerances reverse_only_ad_tolerances(double val_tol = 1e-8, return tols; } - - } // namespace test } // namespace stan #endif diff --git a/test/unit/math/relative_tolerance.hpp b/test/unit/math/relative_tolerance.hpp index e9384836cc2..16eaca8cc4d 100644 --- a/test/unit/math/relative_tolerance.hpp +++ b/test/unit/math/relative_tolerance.hpp @@ -85,9 +85,7 @@ class relative_tolerance { * * @return `true` if the tolernace is infinite. */ - bool is_inf() const { - return std::isinf(tol_); - } + bool is_inf() const { return std::isinf(tol_); } private: /** diff --git a/test/unit/math/test_ad.hpp b/test/unit/math/test_ad.hpp index 23ae8f3a6c3..5a8fff9ec1d 100644 --- a/test/unit/math/test_ad.hpp +++ b/test/unit/math/test_ad.hpp @@ -334,27 +334,22 @@ template void expect_ad_derivatives(const ad_tolerances& tols, const G& g, const Eigen::VectorXd& x) { double gx = g(x); - if (!tols.gradient_val_.is_inf() - || !tols.gradient_grad_.is_inf()) { + if (!tols.gradient_val_.is_inf() || !tols.gradient_grad_.is_inf()) { test_gradient(tols, g, x, gx); } #ifndef STAN_MATH_TESTS_REV_ONLY - if (!tols.gradient_fvar_val_.is_inf() - || !tols.gradient_fvar_grad_.is_inf()) { + if (!tols.gradient_fvar_val_.is_inf() || !tols.gradient_fvar_grad_.is_inf()) { test_gradient_fvar(tols, g, x, gx); } - if (!tols.hessian_val_.is_inf() - || !tols.hessian_grad_.is_inf() + if (!tols.hessian_val_.is_inf() || !tols.hessian_grad_.is_inf() || !tols.hessian_hessian_.is_inf()) { test_hessian(tols, g, x, gx); } - if (!tols.hessian_fvar_val_.is_inf() - || !tols.hessian_fvar_grad_.is_inf() + if (!tols.hessian_fvar_val_.is_inf() || !tols.hessian_fvar_grad_.is_inf() || !tols.hessian_hessian_.is_inf()) { test_hessian_fvar(tols, g, x, gx); } - if (!tols.grad_hessian_val_.is_inf() - || !tols.grad_hessian_hessian_.is_inf() + if (!tols.grad_hessian_val_.is_inf() || !tols.grad_hessian_hessian_.is_inf() || !tols.grad_hessian_grad_hessian_.is_inf()) { test_grad_hessian(tols, g, x, gx); } From f2e4dba181967106cd94d74d638d6e4845f0034c Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Mon, 18 Apr 2022 16:57:54 -0400 Subject: [PATCH 15/18] fft done w/o analytic rev mode --- stan/math/prim/fun/fft.hpp | 36 +++++----- test/unit/math/ad_tolerances.hpp | 2 +- test/unit/math/prim/fun/fft_test.cpp | 101 +++++++++++---------------- 3 files changed, 60 insertions(+), 79 deletions(-) diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index 76f4f8ad8d8..19241771367 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -14,8 +14,7 @@ namespace math { /** * Return the discrete Fourier transform of the specified complex - * vector. The input vector may be considered to be in the time - * domain and the output will be in the frequency domain. + * vector. * * Given an input complex vector `x[0:N-1]` of size `N`, the discrete * Fourier transform computes entries of the resulting complex @@ -28,22 +27,22 @@ namespace math { * If the input is of size zero, the result is a size zero vector. * * @tparam V type of complex vector argument - * @param x complex time domain vector to transform + * @param[in] x vector to transform * @return discrete Fourier transform of `x` */ template * = nullptr> inline Eigen::Matrix, -1, 1> fft(const V& x) { + // copy because fft() requires Eigen::Matrix type Eigen::Matrix, -1, 1> xv = x; if (xv.size() <= 1) return xv; - Eigen::FFT::value_type> fft; + Eigen::FFT> fft; return fft.fwd(xv); } /** * Return the inverse discrete Fourier transform of the specified - * complex vector. The input may be considered to be in the - * frequency domain and the output will be in the time domain. + * complex vector. * * Given an input complex vector `y[0:N-1]` of size `N`, the inverse * discrete Fourier transform computes entries of the resulting @@ -58,29 +57,30 @@ inline Eigen::Matrix, -1, 1> fft(const V& x) { * the sign of the exponent. * * @tparam V type of complex vector argument - * @param y complex frequency domain vector to inverse transform - * @return inverse discrete Fourier transform of `x` + * @param[in] y vector to inverse transform + * @return inverse discrete Fourier transform of `y` */ template * = nullptr> inline Eigen::Matrix, -1, 1> inv_fft(const V& y) { - Eigen::Matrix, -1, 1> yv = y; + // copy because fft() requires Eigen::Matrix type + Eigen::Matrix, -1, 1> yv = y; if (y.size() <= 1) return yv; - Eigen::FFT::value_type> fft; + Eigen::FFT> fft; return fft.inv(yv); } /** * Return the two-dimensional discrete Fourier transform of the - * specified complex matrix. The 2D discrete Fourier transform - * first runs the DFT on the each row, then on each column of the - * result. + * specified complex matrix. The 2D discrete Fourier transform first + * runs the discrete Fourier transform on the each row, then on each + * column of the result. * * @tparam M type of complex matrix argument - * @param x complex time-domain matrix + * @param[in] x matrix to transform * @return discrete 2D Fourier transform of `x` */ -template * = nullptr> +template * = nullptr> inline Eigen::Matrix, -1, -1> fft2(const M& x) { Eigen::Matrix, -1, -1> y(x.rows(), x.cols()); for (int i = 0; i < y.rows(); ++i) @@ -98,10 +98,10 @@ inline Eigen::Matrix, -1, -1> fft2(const M& x) { * FFT and inverse FFT (or vice-versa) is the identity. * * @tparam M type of complex matrix argument - * @param y complex frequency-domain matrix - * @return inverse discrete 2D Fourier transform of `x` + * @param[in] y matrix to inverse trnasform + * @return inverse discrete 2D Fourier transform of `y` */ -template * = nullptr> +template * = nullptr> inline Eigen::Matrix, -1, -1> inv_fft2(const M& y) { Eigen::Matrix, -1, -1> x(y.rows(), y.cols()); for (int j = 0; j < x.cols(); ++j) diff --git a/test/unit/math/ad_tolerances.hpp b/test/unit/math/ad_tolerances.hpp index 7465bc44b71..a9a2087ab64 100644 --- a/test/unit/math/ad_tolerances.hpp +++ b/test/unit/math/ad_tolerances.hpp @@ -81,7 +81,7 @@ ad_tolerances reverse_only_ad_tolerances(double val_tol = 1e-8, ad_tolerances tols; tols.gradient_val_ = val_tol; tols.gradient_grad_ = grad_tol; - double inf = std::numeric_limits::infinity(); + constexpr double inf = std::numeric_limits::infinity(); tols.gradient_grad_varmat_matvar_ = inf; tols.gradient_fvar_val_ = inf; tols.gradient_fvar_grad_ = inf; diff --git a/test/unit/math/prim/fun/fft_test.cpp b/test/unit/math/prim/fun/fft_test.cpp index 988d5ce6da7..bcd2830be5e 100644 --- a/test/unit/math/prim/fun/fft_test.cpp +++ b/test/unit/math/prim/fun/fft_test.cpp @@ -3,8 +3,8 @@ #include TEST(primFun, fft) { - typedef std::complex c_t; - typedef Eigen::Matrix, -1, 1> cv_t; + using c_t = std::complex; + using cv_t = Eigen::Matrix, -1, 1>; using stan::math::fft; // reference answers calculated using Scipy.fft with double precision @@ -45,20 +45,37 @@ TEST(primFun, fft) { EXPECT_NEAR(imag(yb[2]), 2 * -6.53589838, 1e-6); } +template +void expect_complex_mat_eq(const T& x, const U& y, double tol = 1e-8) { + EXPECT_EQ(x.rows(), y.rows()); + EXPECT_EQ(x.cols(), y.cols()); + for (int j = 0; j < x.cols(); ++j) { + for (int i = 0; i < x.rows(); ++i) { + EXPECT_FLOAT_EQ(real(x(i, j)), real(y(i, j))); + EXPECT_FLOAT_EQ(imag(x(i, j)), imag(y(i, j))); + } + } +} + TEST(primFun, inv_fft) { - typedef std::complex c_t; - typedef Eigen::Matrix, -1, 1> cv_t; + using c_t = std::complex; + using cv_t = Eigen::Matrix, -1, 1>; using stan::math::inv_fft; // reference answers calculated using Scipy.fft with double precision cv_t y0(0); cv_t x0 = inv_fft(y0); - EXPECT_EQ(0, x0.size()); + cv_t x0_expected(0); + EXPECT_EQ(x0_expected, x0); cv_t y1(1); y1 << c_t(-3.247, 1.98555); cv_t x1 = inv_fft(y1); + cv_t x1_expected(1); + x1_expected << c_t(-3.247, 1.98555); + expect_complex_mat_eq(x1_expected, x1); + EXPECT_EQ(1, x1.size()); EXPECT_EQ(real(x1[0]), -3.247); EXPECT_EQ(imag(x1[0]), 1.98555); @@ -68,17 +85,14 @@ TEST(primFun, inv_fft) { c_t(11.19615242, -6.53589838); cv_t x = inv_fft(y); EXPECT_EQ(3, y.size()); - EXPECT_NEAR(real(x[0]), 1, 1e-6); - EXPECT_NEAR(imag(x[0]), -2, 1e-6); - EXPECT_NEAR(real(x[1]), -3, 1e-6); - EXPECT_NEAR(imag(x[1]), 5, 1e-6); - EXPECT_NEAR(real(x[2]), -7, 1e-6); - EXPECT_NEAR(imag(x[2]), 11, 1e-6); + Eigen::VectorXcd x_expected(3); + x_expected << c_t(1, -2), c_t(-3, 5), c_t(-7, 11); + expect_complex_mat_eq(x_expected, x); } TEST(primFun, fft2) { - typedef std::complex c_t; - typedef Eigen::Matrix, -1, -1> cm_t; + using c_t = std::complex; + using cm_t = Eigen::Matrix, -1, -1>; using stan::math::fft2; using stan::math::inv_fft2; @@ -99,57 +113,24 @@ TEST(primFun, fft2) { cm_t x12(1, 2); x12 << c_t(1.0, -3.9), c_t(-8.6, 0.2); cm_t y12 = fft2(x12); - EXPECT_EQ(1, y12.rows()); - EXPECT_EQ(2, y12.cols()); - EXPECT_NEAR(-7.6, real(y12(0, 0)), 1e-6); - EXPECT_NEAR(-3.7, imag(y12(0, 0)), 1e-6); - EXPECT_NEAR(9.6, real(y12(0, 1)), 1e-6); - EXPECT_NEAR(-4.1, imag(y12(0, 1)), 1e-6); + cm_t y12_expected(1, 2); + y12_expected << c_t(-7.6, -3.7), c_t(9.6, -4.1); + expect_complex_mat_eq(y12_expected, y12); cm_t x33(3, 3); x33 << c_t(1, 2), c_t(3, -1.4), c_t(2, 1), c_t(3, -9), c_t(2, -1.3), c_t(3.9, -1.8), c_t(13, -1.8), c_t(1.3, 1.9), c_t(-2.2, -2.2); cm_t y33 = fft2(x33); - EXPECT_EQ(3, y33.rows()); - EXPECT_EQ(3, y33.cols()); - - // check vs. results from R (matches SciPy) - EXPECT_NEAR(27, real(y33(0, 0)), 1e-6); - EXPECT_NEAR(-12.6, imag(y33(0, 0)), 1e-6); - EXPECT_NEAR(13.90525589, real(y33(0, 1)), 1e-6); - EXPECT_NEAR(-9.15166605, imag(y33(0, 1)), 1e-6); - EXPECT_NEAR(10.09474411, real(y33(0, 2)), 1e-1); - EXPECT_NEAR(-4.64833395, imag(y33(0, 2)), 2e-1); - EXPECT_NEAR(-13.160254038, real(y33(1, 0)), 1e-1); - EXPECT_NEAR(11.471281292, imag(y33(1, 0)), 2e-1); - EXPECT_NEAR(-13.29326674, real(y33(1, 1)), 1e-1); - EXPECT_NEAR(20.88153533, imag(y33(1, 1)), 1e-1); - EXPECT_NEAR(-13.25262794, real(y33(1, 2)), 1e-1); - EXPECT_NEAR(15.82794549, imag(y33(1, 2)), 1e-1); - EXPECT_NEAR(4.160254038, real(y33(2, 0)), 1e-1); - EXPECT_NEAR(5.928718708, imag(y33(2, 0)), 1e-1); - EXPECT_NEAR(-11.34737206, real(y33(2, 1)), 1e-1); - EXPECT_NEAR(-7.72794549, imag(y33(2, 1)), 1e-1); - EXPECT_NEAR(4.89326674, real(y33(2, 2)), 1e-1); - EXPECT_NEAR(-1.98153533, imag(y33(2, 2)), 1e-1); + cm_t y33_expected(3, 3); + y33_expected << c_t(27, -12.6), c_t(13.90525589, -9.15166605), c_t(10.09474411, -4.64833395), + c_t(-13.160254038, 11.471281292), c_t(-13.29326674, 20.88153533), c_t(-13.25262794, 15.82794549), + c_t(4.160254038, 5.928718708), c_t(-11.34737206, -7.72794549), c_t(4.89326674, -1.98153533); + expect_complex_mat_eq(y33_expected, y33); } - -void expect_complex_matrix_float_eq( - const Eigen::Matrix, -1, -1>& x, - const Eigen::Matrix, -1, -1>& y) { - EXPECT_EQ(x.rows(), y.rows()); - EXPECT_EQ(x.cols(), y.cols()); - for (int j = 0; j < y.cols(); ++j) { - for (int i = 0; i < y.rows(); ++i) { - EXPECT_FLOAT_EQ(real(x(i, j)), real(y(i, j))); - EXPECT_FLOAT_EQ(imag(x(i, j)), imag(y(i, j))); - } - } -} - + TEST(primFunFFT, invfft2) { - typedef std::complex c_t; - typedef Eigen::Matrix, -1, -1> cm_t; + using c_t = std::complex; + using cm_t = Eigen::Matrix, -1, -1>; using stan::math::fft2; using stan::math::inv_fft; using stan::math::inv_fft2; @@ -170,7 +151,7 @@ TEST(primFunFFT, invfft2) { x13 << c_t(-2.3, 1.82), c_t(1.18, 9.32), c_t(1.15, -14.1); cm_t y13 = inv_fft2(x13); cm_t y13copy = inv_fft(x13.row(0)); - expect_complex_matrix_float_eq(y13, y13copy.transpose()); + expect_complex_mat_eq(y13, y13copy.transpose()); cm_t x33(3, 3); x33 << c_t(1, 2), c_t(3, -1.4), c_t(2, 1), c_t(3, -9), c_t(2, -1.3), @@ -199,10 +180,10 @@ TEST(primFunFFT, invfft2) { // check round trips inv_fft(fft(x)) cm_t x33copy = inv_fft2(y33); - expect_complex_matrix_float_eq(x33, x33copy); + expect_complex_mat_eq(x33, x33copy); // check round trip fft(inv_fft(x)) cm_t z33 = inv_fft2(x33); cm_t x33copy2 = fft2(z33); - expect_complex_matrix_float_eq(x33, x33copy2); + expect_complex_mat_eq(x33, x33copy2); } From a1a4cd2a79f15cc933339a251202b6859d78aadb Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Mon, 18 Apr 2022 16:59:04 -0400 Subject: [PATCH 16/18] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/fft.hpp | 2 +- test/unit/math/prim/fun/fft_test.cpp | 12 +++++++----- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/fun/fft.hpp b/stan/math/prim/fun/fft.hpp index 19241771367..ffca15a0108 100644 --- a/stan/math/prim/fun/fft.hpp +++ b/stan/math/prim/fun/fft.hpp @@ -63,7 +63,7 @@ inline Eigen::Matrix, -1, 1> fft(const V& x) { template * = nullptr> inline Eigen::Matrix, -1, 1> inv_fft(const V& y) { // copy because fft() requires Eigen::Matrix type - Eigen::Matrix, -1, 1> yv = y; + Eigen::Matrix, -1, 1> yv = y; if (y.size() <= 1) return yv; Eigen::FFT> fft; diff --git a/test/unit/math/prim/fun/fft_test.cpp b/test/unit/math/prim/fun/fft_test.cpp index bcd2830be5e..78a97bee665 100644 --- a/test/unit/math/prim/fun/fft_test.cpp +++ b/test/unit/math/prim/fun/fft_test.cpp @@ -75,7 +75,7 @@ TEST(primFun, inv_fft) { cv_t x1_expected(1); x1_expected << c_t(-3.247, 1.98555); expect_complex_mat_eq(x1_expected, x1); - + EXPECT_EQ(1, x1.size()); EXPECT_EQ(real(x1[0]), -3.247); EXPECT_EQ(imag(x1[0]), 1.98555); @@ -122,12 +122,14 @@ TEST(primFun, fft2) { c_t(3.9, -1.8), c_t(13, -1.8), c_t(1.3, 1.9), c_t(-2.2, -2.2); cm_t y33 = fft2(x33); cm_t y33_expected(3, 3); - y33_expected << c_t(27, -12.6), c_t(13.90525589, -9.15166605), c_t(10.09474411, -4.64833395), - c_t(-13.160254038, 11.471281292), c_t(-13.29326674, 20.88153533), c_t(-13.25262794, 15.82794549), - c_t(4.160254038, 5.928718708), c_t(-11.34737206, -7.72794549), c_t(4.89326674, -1.98153533); + y33_expected << c_t(27, -12.6), c_t(13.90525589, -9.15166605), + c_t(10.09474411, -4.64833395), c_t(-13.160254038, 11.471281292), + c_t(-13.29326674, 20.88153533), c_t(-13.25262794, 15.82794549), + c_t(4.160254038, 5.928718708), c_t(-11.34737206, -7.72794549), + c_t(4.89326674, -1.98153533); expect_complex_mat_eq(y33_expected, y33); } - + TEST(primFunFFT, invfft2) { using c_t = std::complex; using cm_t = Eigen::Matrix, -1, -1>; From 593e5ec5ebaa33db2f5b98293bed90dac93fc035 Mon Sep 17 00:00:00 2001 From: Bob Carpenter Date: Wed, 20 Apr 2022 16:45:38 -0400 Subject: [PATCH 17/18] requested fixes to testing; fixes #20 --- test/unit/math/mix/fun/fft_test.cpp | 80 +++++++++++++---------------- 1 file changed, 36 insertions(+), 44 deletions(-) diff --git a/test/unit/math/mix/fun/fft_test.cpp b/test/unit/math/mix/fun/fft_test.cpp index dd145fa5908..e0b9a644999 100644 --- a/test/unit/math/mix/fun/fft_test.cpp +++ b/test/unit/math/mix/fun/fft_test.cpp @@ -2,19 +2,17 @@ void expect_fft(const Eigen::VectorXcd& x) { stan::test::ad_tolerances tols = stan::test::reverse_only_ad_tolerances(); - int i = 0; - int m = 0; - auto g = [&](const auto& x) { - using stan::math::fft; - return i ? real(fft(x)(m)) : imag(fft(x)(m)); - }; - for (int i = 0; i < 2; ++i) - for (int m = 0; m < x.rows(); ++m) - stan::test::expect_ad(tols, g, x); + for (int m = 0; m < x.rows(); ++m) { + auto g = [m](const auto& x) { + using stan::math::fft; + return fft(x)(m); + }; + stan::test::expect_ad(tols, g, x); + } } TEST(mathMixFun, fft) { - typedef Eigen::VectorXcd cvec_t; + using cvec_t = Eigen::VectorXcd; cvec_t x0(0); expect_fft(x0); @@ -44,19 +42,17 @@ TEST(mathMixFun, fft) { void expect_inv_fft(const Eigen::VectorXcd& x) { stan::test::ad_tolerances tols = stan::test::reverse_only_ad_tolerances(); - int i = 0; - int m = 0; - auto g = [&](const auto& x) { - using stan::math::inv_fft; - return i ? real(inv_fft(x)(m)) : imag(inv_fft(x)(m)); - }; - for (int i = 0; i < 2; ++i) - for (int m = 0; m < x.rows(); ++m) - stan::test::expect_ad(tols, g, x); + for (int m = 0; m < x.rows(); ++m) { + auto g = [m](const auto& x) { + using stan::math::inv_fft; + return inv_fft(x)(m); + }; + stan::test::expect_ad(tols, g, x); + } } TEST(mathMixFun, invFft) { - typedef Eigen::VectorXcd cvec_t; + using cvec_t = Eigen::VectorXcd; cvec_t x0(0); expect_inv_fft(x0); @@ -86,21 +82,19 @@ TEST(mathMixFun, invFft) { void expect_fft2(const Eigen::MatrixXcd& x) { stan::test::ad_tolerances tols = stan::test::reverse_only_ad_tolerances(); - int i = 0; - int m = 0; - int n = 0; - auto g = [&](const auto& x) { - using stan::math::fft2; - return i ? real(fft2(x)(m, n)) : imag(fft2(x)(m, n)); - }; - for (int i = 0; i < 2; ++i) - for (int n = 0; n < x.cols(); ++n) - for (int m = 0; m < x.rows(); ++m) - stan::test::expect_ad(tols, g, x); + for (int n = 0; n < x.cols(); ++n) { + for (int m = 0; m < x.rows(); ++m) { + auto g = [m, n](const auto& x) { + using stan::math::fft2; + return fft2(x)(m, n); + }; + stan::test::expect_ad(tols, g, x); + } + } } TEST(mathMixFun, fft2) { - typedef Eigen::MatrixXcd cmat_t; + using cmat_t = Eigen::MatrixXcd; cmat_t x00(0, 0); expect_fft2(x00); @@ -136,21 +130,19 @@ TEST(mathMixFun, fft2) { void expect_inv_fft2(const Eigen::MatrixXcd& x) { stan::test::ad_tolerances tols = stan::test::reverse_only_ad_tolerances(); - int i = 0; - int m = 0; - int n = 0; - auto g = [&](const auto& x) { - using stan::math::inv_fft2; - return i ? real(inv_fft2(x)(m, n)) : imag(inv_fft2(x)(m, n)); - }; - for (int i = 0; i < 2; ++i) - for (int n = 0; n < x.cols(); ++n) - for (int m = 0; m < x.rows(); ++m) - stan::test::expect_ad(tols, g, x); + for (int n = 0; n < x.cols(); ++n) { + for (int m = 0; m < x.rows(); ++m) { + auto g = [m, n](const auto& x) { + using stan::math::inv_fft2; + return inv_fft2(x)(m, n); + }; + stan::test::expect_ad(tols, g, x); + } + } } TEST(mathMixFun, inv_fft2) { - typedef Eigen::MatrixXcd cmat_t; + using cmat_t = Eigen::MatrixXcd; cmat_t x00(0, 0); expect_inv_fft2(x00); From cc813b454098cba4fc4cccb51c091438685a91ff Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Wed, 20 Apr 2022 16:46:38 -0400 Subject: [PATCH 18/18] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- test/unit/math/mix/fun/fft_test.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/unit/math/mix/fun/fft_test.cpp b/test/unit/math/mix/fun/fft_test.cpp index e0b9a644999..78ec560693a 100644 --- a/test/unit/math/mix/fun/fft_test.cpp +++ b/test/unit/math/mix/fun/fft_test.cpp @@ -85,8 +85,8 @@ void expect_fft2(const Eigen::MatrixXcd& x) { for (int n = 0; n < x.cols(); ++n) { for (int m = 0; m < x.rows(); ++m) { auto g = [m, n](const auto& x) { - using stan::math::fft2; - return fft2(x)(m, n); + using stan::math::fft2; + return fft2(x)(m, n); }; stan::test::expect_ad(tols, g, x); } @@ -133,8 +133,8 @@ void expect_inv_fft2(const Eigen::MatrixXcd& x) { for (int n = 0; n < x.cols(); ++n) { for (int m = 0; m < x.rows(); ++m) { auto g = [m, n](const auto& x) { - using stan::math::inv_fft2; - return inv_fft2(x)(m, n); + using stan::math::inv_fft2; + return inv_fft2(x)(m, n); }; stan::test::expect_ad(tols, g, x); }