From f340d6fcfd7c84389a61571ef129d0f668116851 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Tue, 31 Mar 2020 15:56:39 -0400 Subject: [PATCH 01/57] added interpolation functions --- stan/math/prim/fun.hpp | 2 + stan/math/prim/fun/conv_gaus_line.hpp | 65 +++++++ stan/math/prim/fun/gaus_interp.hpp | 174 ++++++++++++++++++ stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/conv_gaus_line.hpp | 44 +++++ .../unit/math/mix/fun/conv_gaus_line_test.cpp | 85 +++++++++ test/unit/math/mix/fun/gaus_interp_test.cpp | 144 +++++++++++++++ 7 files changed, 515 insertions(+) create mode 100644 stan/math/prim/fun/conv_gaus_line.hpp create mode 100644 stan/math/prim/fun/gaus_interp.hpp create mode 100644 stan/math/rev/fun/conv_gaus_line.hpp create mode 100644 test/unit/math/mix/fun/conv_gaus_line_test.cpp create mode 100644 test/unit/math/mix/fun/gaus_interp_test.cpp diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 72d6a26d557..1d41bcfb4e3 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include #include @@ -95,6 +96,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp new file mode 100644 index 00000000000..8fe6c6a3855 --- /dev/null +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -0,0 +1,65 @@ +#ifndef STAN_MATH_PRIM_FUN_CONV_GAUS_LINE +#define STAN_MATH_PRIM_FUN_CONV_GAUS_LINE + +#include +#include +#include +#include + +namespace stan { +namespace math { + +/* +evaluate the derivative of conv_gaus_line with respect to x +*/ +double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, + double sig2) { + using stan::math::normal_cdf; + double pi = stan::math::pi(); + using std::exp; + using std::pow; + using std::sqrt; + double sig = sqrt(sig2); + double y; + + double alpha = sqrt(2 * pi * sig2); + y = (a * x0 + b) / alpha * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) + + exp(-pow(t0 - x0, 2) / (2 * sig2))); + y += a * (normal_cdf(t1, x0, sig) - normal_cdf(t0, x0, sig)); + y -= a / alpha * ((t1-x0)*exp(-pow(t1 - x0, 2) / (2 * sig2)) + -(t0-x0)*exp(-pow(t0 - x0, 2) / (2 * sig2))); + return y; +} + + +/* +evaluate the integral + + | t1 +(2*pi*sig2)**-(1/2) | (a*t + b) * exp(-(t - x0)^2 / (2*sig2)) dt + | t0 + +*/ +template +double conv_gaus_line(double t0, double t1, double a, double b, Tx const& x, + double sig2) { + using stan::math::normal_cdf; + double pi = stan::math::pi(); + using std::exp; + using std::pow; + using std::sqrt; + double sig = sqrt(sig2); + + double y = (a * value_of(x) + b) * (normal_cdf(t1, value_of(x), sig) + - normal_cdf(t0, value_of(x), sig)); + y += - a * sig2 / sqrt(2 * pi * sig2) * + (exp(-pow(t1 - value_of(x), 2)/ (2 * sig2)) + - exp(-pow(t0 - value_of(x), 2) / (2 * sig2))); + + return y; +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp new file mode 100644 index 00000000000..b77b58b35c8 --- /dev/null +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -0,0 +1,174 @@ +#ifndef STAN_MATH_PRIM_FUN_GAUS_INTERP +#define STAN_MATH_PRIM_FUN_GAUS_INTERP + +#include +#include +#include + +using namespace std; + +namespace stan { +namespace math { + +// set tolerance for how close interpolated values need to be to +// the specified function values +const double INTERP_TOL = 1e-8; +const double SIG2_SCALE = 0.1; + +/* +given a set of points (x_i, y_i) with x_i's in increasing order, find +a_i, b_i such that the line between (x_i, y_i) and (x_{i+1}, y_{i+1}) is +f(t) = a_i*t + b_i, store the a_i's and b_i's in as and bs. +*/ +void lin_interp_coefs(int n, std::vector xs,std::vector ys, + std::vector& as, std::vector& bs) { + // initialize size of vectors of linear coefficients + as.resize(n - 1); + bs.resize(n - 1); + + // find slope and intercept between each point + for (int i=0; i xs, + vector ys, + vector as, + vector bs, + double x) { + // find interval where x lives + if (x <= xs[0]) return ys[0]; + if (x >= xs[n-1]) return ys[n-1]; + + auto lb = upper_bound(xs.begin(), xs.end(), x); + int ind = distance(xs.begin(), lb) - 1; + ind = std::max(0, ind); + + return as[ind]*x + bs[ind]; +} + +/* +linear interpolation at a vector of points +*/ +vector lin_interp(int n, + std::vector xs, + std::vector ys, + int n_new, + std::vector xs_new) { + // compute coefficients of linear interpolation + vector as, bs; + lin_interp_coefs(n, xs, ys, as, bs); + + // evaluate at new points + vector ys_new; + for (int i = 0; i < n_new; i++) { + ys_new.push_back(lin_interp_pt(n, xs, ys, as, bs, xs_new[i])); + } + return ys_new; +} + +/* +given a set of points and a width for the gaussian kernel, do a convolution +and evaluate at one point, x +*/ +template +inline return_type_t interp_gaus_pt(int n, std::vector xs, + std::vector ys, + std::vector as, + std::vector bs, + Tx const& x, double sig2) { + + // extend out first and last lines for convolution + double sig = std::sqrt(sig2); + xs[0] += - 10 * sig; + xs[n-1] += 10 * sig; + + // no need to convolve far from center of gaussian, so + // get lower and upper indexes for integration bounds + auto lb = upper_bound(xs.begin(), xs.end(), x - 10 * sig); + int ind_start = distance(xs.begin(), lb) - 1; + ind_start = std::max(0, ind_start); + + auto ub = upper_bound(xs.begin(), xs.end(), x + 10 * sig); + int ind_end = distance(xs.begin(), ub); + ind_end = std::min(n - 1, ind_end); + + // sum convolutions over intervals + using return_t = return_type_t; + return_t y = 0; + for (int i = ind_start; i < ind_end; i++) { + y += conv_gaus_line(xs[i], xs[i + 1], as[i], bs[i], x, sig2); + } + + return y; +} + +/* +find the smallest difference between successive elements in a sorted vector +*/ +template +double min_diff(int n, std::vector xs) { + double dmin = value_of(xs[1]) - value_of(xs[0]); + for (int i = 1; i < n - 1; i++) { + if (value_of(xs[i + 1]) - value_of(xs[i]) < dmin) { + dmin = value_of(xs[i + 1]) - value_of(xs[i]); + } + } + return dmin; +} + +/* +given a set of pairs (x_i, y_i), do a gaussian interpolation through those +points and evaluate the interpolation at the points xs_new +*/ +template +inline vector gaus_interp(int n, + std::vector xs, + std::vector ys, + int n_new, + std::vector xs_new) { + + // find minimum distance between points for std of gaussian kernel + double sig2 = square(min_diff(n, xs) * SIG2_SCALE); + + // copy ys into a new vector + std::vector y2s; + y2s.resize(n); + for (int i = 0; i < n; i++) { + y2s[i] = ys[i]; + } + + // interatively find interpolation that coincides with ys at xs + int max_iters = 50; + double dmax, dd; + std::vector as, bs; + for (int j=0; j < max_iters; j++) { + // linear interpolation for new ys + lin_interp_coefs(n, xs, y2s, as, bs); + dmax = 0; + for (int i = 0; i < n; i++) { + dd = ys[i] - interp_gaus_pt(n, xs, y2s, as, bs, xs[i], sig2); + y2s[i] += dd; + dmax = std::max(std::abs(dd), dmax); + } + if (dmax < INTERP_TOL) break; + } + + // fill ys_new with interpolated values at new_xs + std::vector ys_new; + for (int i = 0; i < n_new; i++) { + ys_new.push_back(interp_gaus_pt(n, xs, y2s, as, bs, xs_new[i], sig2)); + } + return ys_new; +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index b26b0c057fe..236efe26086 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp new file mode 100644 index 00000000000..fa116f01ffd --- /dev/null +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -0,0 +1,44 @@ +#ifndef STAN_MATH_REV_FUN_CONV_GAUS_LINE +#define STAN_MATH_REV_FUN_CONV_GAUS_LINE + +#include +#include +#include + +namespace stan { +namespace math { + +namespace internal { +class conv_gaus_line_vari : public op_v_vari { + double t0_; + double t1_; + double a_; + double b_; + double sig2_; + + public: + explicit conv_gaus_line_vari(double t0, double t1, double a, double b, + vari* avi, double sig2) : + t0_(t0), + t1_(t1), + a_(a), + b_(b), + sig2_(sig2), + op_v_vari(conv_gaus_line(t0, t1, a, b, avi->val_, sig2), avi) {} + + void chain() { avi_->adj_ += adj_ * der_conv_gaus_line(t0_, t1_, a_, b_, + avi_->val_, sig2_); } +}; + +} // namespace internal + + +inline var conv_gaus_line(double t0, double t1, double a, double b, + const var& x, double sig2) { + return var(new internal::conv_gaus_line_vari(t0, t1, a, b, x.vi_, sig2)); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/mix/fun/conv_gaus_line_test.cpp b/test/unit/math/mix/fun/conv_gaus_line_test.cpp new file mode 100644 index 00000000000..990e3e71fef --- /dev/null +++ b/test/unit/math/mix/fun/conv_gaus_line_test.cpp @@ -0,0 +1,85 @@ +#include +#include +#include + +/* +evaluate the convolution using trapezoid rule +*/ +double check_int_dumb(double t0, double t1, double a, double b, double x0, + double sig2) { + double x, y, t, tot; + vector ys, ts; + double pi = stan::math::pi(); + int n; + + // evaluate function at equispaced nodes + n = 100000; + for (int i = 0; i < n; i++) { + t = t0 + i * (t1 - t0) / (n - 1); + ts.push_back(t); + + y = a * t + b; + y = y * exp(-pow(t - x0, 2) / (2 * sig2)); + ys.push_back(y); + } + + // sum up function evaluations + tot = 0; + for (int i = 0; i < n; i++) { + tot += ys[i]; + } + tot *= (t1 - t0) / n; + tot /= sqrt(2 * pi * sig2); + + return tot; +} + +/* +check that integral using formula is close to trapezoid approximation +*/ +TEST(mathMixConvGausLin, trap_test) { + double t0, t1, a, b, x0, sig2, y, y2; + using stan::math::conv_gaus_line; + t0 = 0; + t1 = 1; + a = 1; + b = 0; + x0 = 1; + sig2 = 3; + y = check_int_dumb(t0, t1, a, b, x0, sig2); + y2 = conv_gaus_line(t0, t1, a, b, x0, sig2); + + ASSERT_NEAR(y, y2, 1e-5); +} + +/* +check derivative using finite differences and autodiffing +*/ +TEST(mathMixGausInterp, gaus_conv_der) { + double t0, t1, a, b, x0, sig2, y, h, yn, yp, x0p, x0n, dder, dder2; + double alpha, sig; + using stan::math::conv_gaus_line; + + t0 = 0; + t1 = 5; + a = -1; + b = 2; + x0 = 0.99; + sig2 = 2; + + // derivative with finite difference + h = 1e-4; + x0n = x0 + h; + x0p = x0 - h; + yn = conv_gaus_line(t0, t1, a, b, x0n, sig2); + yp = conv_gaus_line(t0, t1, a, b, x0p, sig2); + dder = (yn - yp) / (2*h); + + // derivative with autodiff + using stan::math::var; + var v = x0; + var vy = conv_gaus_line(t0, t1, a, b, v, sig2); + vy.grad(); + + ASSERT_NEAR(dder, v.adj(), 1e-5); +} diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp new file mode 100644 index 00000000000..803e879a653 --- /dev/null +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -0,0 +1,144 @@ +#include +#include +#include + +double ABS_TOL = 1e-12; + +TEST(mathMixGausInterp, interp_line) { + using stan::math::gaus_interp; + + // check that interpolation of line returns the same function + // generate function tabulation + int n = 2; + double xmin = 0; + double xmax = 1; + vector xs = {0, 1}; + vector ys = {0, 2}; + + // vector of pts at which to compute interpolation + vector xs_new; + int n_interp = 100; + double t0 = xmin; + double t1 = xmax; + double t; + for (int i=0; i < n_interp; i++) { + t = t0 + i * (t1 - t0) / (n_interp - 1); + xs_new.push_back(t); + } + + // create interpolation + vector ys_new; + ys_new = gaus_interp(n, xs, ys, n_interp, xs_new); + + // test points + double tmp, y; + for (int i=0; i < n_interp; i++) { + tmp = (ys[1] - ys[0]) / (xs[1] - xs[0]); + y = tmp * xs_new[i] + ys[0] - tmp * xs[0]; + ASSERT_NEAR(ys_new[i], y, ABS_TOL); + } +} + +TEST(mathMixGausInterp, gaus_and_lin_interp) { + using stan::math::gaus_interp; + using stan::math::lin_interp; + + // check that interpolation of line returns the same function + // generate function tabulation + int n = 30; + double xmin = 0; + double xmax = 1; + double x; + vector xs, ys; + for (int i=0; i < n; i++) { + x = xmin + i * (xmax - xmin) / (n - 1); + xs.push_back(x); + ys.push_back(x*x); + } + + // vector of pts at which to compute interpolation + vector xs_new; + int n_interp = 100; + double t0 = xmin; + double t1 = xmax; + double t; + for (int i=0; i < n_interp; i++) { + t = t0 + i * (t1 - t0) / (n_interp - 1); + xs_new.push_back(t); + } + + // create interpolation + vector ys_new_gaus, ys_new_lin; + ys_new_gaus = gaus_interp(n, xs, ys, n_interp, xs_new); + ys_new_lin = lin_interp(n, xs, ys, n_interp, xs_new); + + // test points + double tmp, y; + for (int i=0; i < n_interp; i++) { + ASSERT_NEAR(ys_new_lin[i], ys_new_gaus[i], 1e-4); + } +} + + +TEST(mathMixGausInterp, derivs) { + using stan::math::var; + using stan::math::gaus_interp; + std::vector xs, ys, x2s, y2s, ts, as, bs, y3s, t3s; + double xmin, xmax, x, y, x2, y2, t, t0, t1, dd, dder, dder2; + double x0, x1, y0, y1, tmp; + int n; + + // generate function tabulation + n = 10; + xmin = 0; + xmax = 1; + for (int i = 0; i < n; i++) { + x = xmin + i * (xmax - xmin) / (n - 1); + xs.push_back(x); + y = rand() % 100; + ys.push_back(y); + } + + // create vector of interpolation pts + std::vector xs_new_v; + std::vector xs_new; + int n_interp = 100; + t0 = xs[1] - 1e-4; + t0 = xs[0]; + t1 = xmax; + for (int i=0; i < n_interp; i++) { + t = t0 + i * (t1 - t0) / (n_interp - 1); + xs_new.push_back(t); + stan::math::var tvar = t; + xs_new_v.push_back(tvar); + } + + std::vector ys_new, ys_new_p, ys_new_n, xs_new_p, xs_new_n; + ys_new = gaus_interp(n, xs, ys, n_interp, xs_new); + + // autodiff at each interpolation pt + std::vector ys_new_v; + ys_new_v = gaus_interp(n, xs, ys, n_interp, xs_new_v); + + std::vector ys_new_dder; + std::vector ys_new_v2; + for (int i=0; i < n_interp; i++) { + ys_new_v[i].grad(); + dder = xs_new_v[i].adj(); + ys_new_dder.push_back(dder); + } + + // take derivative of interpolation using finite differencing + double h = 1e-4; + xs_new_p = xs_new; + xs_new_n = xs_new; + for (int i=0; i < n_interp; i++) { + xs_new_p[i] += -h; + xs_new_n[i] += h; + ys_new_p = gaus_interp(n, xs, ys, n_interp, xs_new_p); + ys_new_n = gaus_interp(n, xs, ys, n_interp, xs_new_n); + dder2 = (ys_new_n[i] - ys_new_p[i]) / (2 * h); + } + + +} From 63172f13b315514eab966ebe9c854b6ac7dd0b24 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Tue, 31 Mar 2020 17:57:10 -0400 Subject: [PATCH 02/57] moved prim tests to new files --- .../unit/math/mix/fun/conv_gaus_line_test.cpp | 50 ----------- test/unit/math/mix/fun/gaus_interp_test.cpp | 82 +------------------ .../math/prim/fun/conv_gaus_line_test.cpp | 53 ++++++++++++ test/unit/math/prim/fun/gaus_interp_test.cpp | 80 ++++++++++++++++++ 4 files changed, 135 insertions(+), 130 deletions(-) create mode 100644 test/unit/math/prim/fun/conv_gaus_line_test.cpp create mode 100644 test/unit/math/prim/fun/gaus_interp_test.cpp diff --git a/test/unit/math/mix/fun/conv_gaus_line_test.cpp b/test/unit/math/mix/fun/conv_gaus_line_test.cpp index 990e3e71fef..58edebf9799 100644 --- a/test/unit/math/mix/fun/conv_gaus_line_test.cpp +++ b/test/unit/math/mix/fun/conv_gaus_line_test.cpp @@ -2,56 +2,6 @@ #include #include -/* -evaluate the convolution using trapezoid rule -*/ -double check_int_dumb(double t0, double t1, double a, double b, double x0, - double sig2) { - double x, y, t, tot; - vector ys, ts; - double pi = stan::math::pi(); - int n; - - // evaluate function at equispaced nodes - n = 100000; - for (int i = 0; i < n; i++) { - t = t0 + i * (t1 - t0) / (n - 1); - ts.push_back(t); - - y = a * t + b; - y = y * exp(-pow(t - x0, 2) / (2 * sig2)); - ys.push_back(y); - } - - // sum up function evaluations - tot = 0; - for (int i = 0; i < n; i++) { - tot += ys[i]; - } - tot *= (t1 - t0) / n; - tot /= sqrt(2 * pi * sig2); - - return tot; -} - -/* -check that integral using formula is close to trapezoid approximation -*/ -TEST(mathMixConvGausLin, trap_test) { - double t0, t1, a, b, x0, sig2, y, y2; - using stan::math::conv_gaus_line; - t0 = 0; - t1 = 1; - a = 1; - b = 0; - x0 = 1; - sig2 = 3; - y = check_int_dumb(t0, t1, a, b, x0, sig2); - y2 = conv_gaus_line(t0, t1, a, b, x0, sig2); - - ASSERT_NEAR(y, y2, 1e-5); -} - /* check derivative using finite differences and autodiffing */ diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index 803e879a653..82dda3902d7 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -2,83 +2,6 @@ #include #include -double ABS_TOL = 1e-12; - -TEST(mathMixGausInterp, interp_line) { - using stan::math::gaus_interp; - - // check that interpolation of line returns the same function - // generate function tabulation - int n = 2; - double xmin = 0; - double xmax = 1; - vector xs = {0, 1}; - vector ys = {0, 2}; - - // vector of pts at which to compute interpolation - vector xs_new; - int n_interp = 100; - double t0 = xmin; - double t1 = xmax; - double t; - for (int i=0; i < n_interp; i++) { - t = t0 + i * (t1 - t0) / (n_interp - 1); - xs_new.push_back(t); - } - - // create interpolation - vector ys_new; - ys_new = gaus_interp(n, xs, ys, n_interp, xs_new); - - // test points - double tmp, y; - for (int i=0; i < n_interp; i++) { - tmp = (ys[1] - ys[0]) / (xs[1] - xs[0]); - y = tmp * xs_new[i] + ys[0] - tmp * xs[0]; - ASSERT_NEAR(ys_new[i], y, ABS_TOL); - } -} - -TEST(mathMixGausInterp, gaus_and_lin_interp) { - using stan::math::gaus_interp; - using stan::math::lin_interp; - - // check that interpolation of line returns the same function - // generate function tabulation - int n = 30; - double xmin = 0; - double xmax = 1; - double x; - vector xs, ys; - for (int i=0; i < n; i++) { - x = xmin + i * (xmax - xmin) / (n - 1); - xs.push_back(x); - ys.push_back(x*x); - } - - // vector of pts at which to compute interpolation - vector xs_new; - int n_interp = 100; - double t0 = xmin; - double t1 = xmax; - double t; - for (int i=0; i < n_interp; i++) { - t = t0 + i * (t1 - t0) / (n_interp - 1); - xs_new.push_back(t); - } - - // create interpolation - vector ys_new_gaus, ys_new_lin; - ys_new_gaus = gaus_interp(n, xs, ys, n_interp, xs_new); - ys_new_lin = lin_interp(n, xs, ys, n_interp, xs_new); - - // test points - double tmp, y; - for (int i=0; i < n_interp; i++) { - ASSERT_NEAR(ys_new_lin[i], ys_new_gaus[i], 1e-4); - } -} - TEST(mathMixGausInterp, derivs) { using stan::math::var; @@ -129,7 +52,7 @@ TEST(mathMixGausInterp, derivs) { } // take derivative of interpolation using finite differencing - double h = 1e-4; + double h = 1e-8; xs_new_p = xs_new; xs_new_n = xs_new; for (int i=0; i < n_interp; i++) { @@ -138,7 +61,6 @@ TEST(mathMixGausInterp, derivs) { ys_new_p = gaus_interp(n, xs, ys, n_interp, xs_new_p); ys_new_n = gaus_interp(n, xs, ys, n_interp, xs_new_n); dder2 = (ys_new_n[i] - ys_new_p[i]) / (2 * h); + ASSERT_NEAR(ys_new_dder[i], dder2, 1e-5); } - - } diff --git a/test/unit/math/prim/fun/conv_gaus_line_test.cpp b/test/unit/math/prim/fun/conv_gaus_line_test.cpp new file mode 100644 index 00000000000..5ee590ddba7 --- /dev/null +++ b/test/unit/math/prim/fun/conv_gaus_line_test.cpp @@ -0,0 +1,53 @@ +#include +#include +#include + +/* +evaluate the convolution using trapezoid rule +*/ +double check_int_dumb(double t0, double t1, double a, double b, double x0, + double sig2) { + double x, y, t, tot; + vector ys, ts; + double pi = stan::math::pi(); + int n; + + // evaluate function at equispaced nodes + n = 100000; + for (int i = 0; i < n; i++) { + t = t0 + i * (t1 - t0) / (n - 1); + ts.push_back(t); + + y = a * t + b; + y = y * exp(-pow(t - x0, 2) / (2 * sig2)); + ys.push_back(y); + } + + // sum up function evaluations + tot = 0; + for (int i = 0; i < n; i++) { + tot += ys[i]; + } + tot *= (t1 - t0) / n; + tot /= sqrt(2 * pi * sig2); + + return tot; +} + +/* +check that integral using formula is close to trapezoid approximation +*/ +TEST(mathMixConvGausLin, trap_test) { + double t0, t1, a, b, x0, sig2, y, y2; + using stan::math::conv_gaus_line; + t0 = 0; + t1 = 1; + a = 1; + b = 0; + x0 = 1; + sig2 = 3; + y = check_int_dumb(t0, t1, a, b, x0, sig2); + y2 = conv_gaus_line(t0, t1, a, b, x0, sig2); + + ASSERT_NEAR(y, y2, 1e-5); +} diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp new file mode 100644 index 00000000000..47ef627f9e4 --- /dev/null +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -0,0 +1,80 @@ +#include +#include +#include + +double ABS_TOL = 1e-12; + +TEST(mathPrimGausInterp, interp_line) { + using stan::math::gaus_interp; + + // check that interpolation of line returns the same function + // generate function tabulation + int n = 2; + double xmin = 0; + double xmax = 1; + vector xs = {0, 1}; + vector ys = {0, 2}; + + // vector of pts at which to compute interpolation + vector xs_new; + int n_interp = 100; + double t0 = xmin; + double t1 = xmax; + double t; + for (int i=0; i < n_interp; i++) { + t = t0 + i * (t1 - t0) / (n_interp - 1); + xs_new.push_back(t); + } + + // create interpolation + vector ys_new; + ys_new = gaus_interp(n, xs, ys, n_interp, xs_new); + + // test points + double tmp, y; + for (int i=0; i < n_interp; i++) { + tmp = (ys[1] - ys[0]) / (xs[1] - xs[0]); + y = tmp * xs_new[i] + ys[0] - tmp * xs[0]; + ASSERT_NEAR(ys_new[i], y, ABS_TOL); + } +} + +TEST(mathPrimGausInterp, gaus_and_lin_interp) { + using stan::math::gaus_interp; + using stan::math::lin_interp; + + // check that interpolation of line returns the same function + // generate function tabulation + int n = 30; + double xmin = 0; + double xmax = 1; + double x; + vector xs, ys; + for (int i=0; i < n; i++) { + x = xmin + i * (xmax - xmin) / (n - 1); + xs.push_back(x); + ys.push_back(x*x); + } + + // vector of pts at which to compute interpolation + vector xs_new; + int n_interp = 100; + double t0 = xmin; + double t1 = xmax; + double t; + for (int i=0; i < n_interp; i++) { + t = t0 + i * (t1 - t0) / (n_interp - 1); + xs_new.push_back(t); + } + + // create interpolation + vector ys_new_gaus, ys_new_lin; + ys_new_gaus = gaus_interp(n, xs, ys, n_interp, xs_new); + ys_new_lin = lin_interp(n, xs, ys, n_interp, xs_new); + + // test points + double tmp, y; + for (int i=0; i < n_interp; i++) { + ASSERT_NEAR(ys_new_lin[i], ys_new_gaus[i], 1e-4); + } +} From 98fdf8bbec8eb8a07df7a99ac5b94dde87bb9654 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 31 Mar 2020 20:43:41 -0400 Subject: [PATCH 03/57] [Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/stable/2017-11-14) --- stan/math/opencl/kernel_generator/load.hpp | 2 +- stan/math/prim/fun/conv_gaus_line.hpp | 32 ++++---- stan/math/prim/fun/gaus_interp.hpp | 76 +++++++++---------- stan/math/rev/fun/conv_gaus_line.hpp | 37 ++++----- .../unit/math/mix/fun/conv_gaus_line_test.cpp | 2 +- test/unit/math/mix/fun/gaus_interp_test.cpp | 9 +-- .../math/prim/fun/conv_gaus_line_test.cpp | 2 +- test/unit/math/prim/fun/gaus_interp_test.cpp | 12 +-- 8 files changed, 83 insertions(+), 89 deletions(-) diff --git a/stan/math/opencl/kernel_generator/load.hpp b/stan/math/opencl/kernel_generator/load.hpp index 799cdb551a0..ae0f7d1eb6d 100644 --- a/stan/math/opencl/kernel_generator/load.hpp +++ b/stan/math/opencl/kernel_generator/load.hpp @@ -49,7 +49,7 @@ class load_ * Creates a deep copy of this expression. * @return copy of \c *this */ - inline load_ deep_copy() const & { return load_(a_); } + inline load_ deep_copy() const& { return load_(a_); } inline load_ deep_copy() && { return load_(std::forward(a_)); } /** diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 8fe6c6a3855..2f1f7d42220 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -12,8 +12,8 @@ namespace math { /* evaluate the derivative of conv_gaus_line with respect to x */ -double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, - double sig2) { +double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, + double sig2) { using stan::math::normal_cdf; double pi = stan::math::pi(); using std::exp; @@ -23,15 +23,16 @@ double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, double y; double alpha = sqrt(2 * pi * sig2); - y = (a * x0 + b) / alpha * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) - + exp(-pow(t0 - x0, 2) / (2 * sig2))); - y += a * (normal_cdf(t1, x0, sig) - normal_cdf(t0, x0, sig)); - y -= a / alpha * ((t1-x0)*exp(-pow(t1 - x0, 2) / (2 * sig2)) - -(t0-x0)*exp(-pow(t0 - x0, 2) / (2 * sig2))); + y = (a * x0 + b) / alpha + * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) + + exp(-pow(t0 - x0, 2) / (2 * sig2))); + y += a * (normal_cdf(t1, x0, sig) - normal_cdf(t0, x0, sig)); + y -= a / alpha + * ((t1 - x0) * exp(-pow(t1 - x0, 2) / (2 * sig2)) + - (t0 - x0) * exp(-pow(t0 - x0, 2) / (2 * sig2))); return y; } - /* evaluate the integral @@ -41,8 +42,8 @@ evaluate the integral */ template -double conv_gaus_line(double t0, double t1, double a, double b, Tx const& x, - double sig2) { +double conv_gaus_line(double t0, double t1, double a, double b, Tx const& x, + double sig2) { using stan::math::normal_cdf; double pi = stan::math::pi(); using std::exp; @@ -50,11 +51,12 @@ double conv_gaus_line(double t0, double t1, double a, double b, Tx const& x, using std::sqrt; double sig = sqrt(sig2); - double y = (a * value_of(x) + b) * (normal_cdf(t1, value_of(x), sig) - - normal_cdf(t0, value_of(x), sig)); - y += - a * sig2 / sqrt(2 * pi * sig2) * - (exp(-pow(t1 - value_of(x), 2)/ (2 * sig2)) - - exp(-pow(t0 - value_of(x), 2) / (2 * sig2))); + double y + = (a * value_of(x) + b) + * (normal_cdf(t1, value_of(x), sig) - normal_cdf(t0, value_of(x), sig)); + y += -a * sig2 / sqrt(2 * pi * sig2) + * (exp(-pow(t1 - value_of(x), 2) / (2 * sig2)) + - exp(-pow(t0 - value_of(x), 2) / (2 * sig2))); return y; } diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index b77b58b35c8..e29ce96b69e 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -11,23 +11,23 @@ namespace stan { namespace math { // set tolerance for how close interpolated values need to be to -// the specified function values +// the specified function values const double INTERP_TOL = 1e-8; const double SIG2_SCALE = 0.1; -/* -given a set of points (x_i, y_i) with x_i's in increasing order, find -a_i, b_i such that the line between (x_i, y_i) and (x_{i+1}, y_{i+1}) is +/* +given a set of points (x_i, y_i) with x_i's in increasing order, find +a_i, b_i such that the line between (x_i, y_i) and (x_{i+1}, y_{i+1}) is f(t) = a_i*t + b_i, store the a_i's and b_i's in as and bs. */ -void lin_interp_coefs(int n, std::vector xs,std::vector ys, - std::vector& as, std::vector& bs) { +void lin_interp_coefs(int n, std::vector xs, std::vector ys, + std::vector& as, std::vector& bs) { // initialize size of vectors of linear coefficients as.resize(n - 1); bs.resize(n - 1); // find slope and intercept between each point - for (int i=0; i xs,std::vector ys, /* linear interpolation at one point, x, given the coefficients */ -double lin_interp_pt(int n, - vector xs, - vector ys, - vector as, - vector bs, - double x) { +double lin_interp_pt(int n, vector xs, vector ys, + vector as, vector bs, double x) { // find interval where x lives - if (x <= xs[0]) return ys[0]; - if (x >= xs[n-1]) return ys[n-1]; + if (x <= xs[0]) + return ys[0]; + if (x >= xs[n - 1]) + return ys[n - 1]; auto lb = upper_bound(xs.begin(), xs.end(), x); int ind = distance(xs.begin(), lb) - 1; ind = std::max(0, ind); - return as[ind]*x + bs[ind]; + return as[ind] * x + bs[ind]; } /* linear interpolation at a vector of points */ -vector lin_interp(int n, - std::vector xs, - std::vector ys, - int n_new, - std::vector xs_new) { +vector lin_interp(int n, std::vector xs, std::vector ys, + int n_new, std::vector xs_new) { // compute coefficients of linear interpolation vector as, bs; lin_interp_coefs(n, xs, ys, as, bs); @@ -78,18 +73,17 @@ given a set of points and a width for the gaussian kernel, do a convolution and evaluate at one point, x */ template -inline return_type_t interp_gaus_pt(int n, std::vector xs, - std::vector ys, - std::vector as, - std::vector bs, - Tx const& x, double sig2) { - +inline return_type_t interp_gaus_pt(int n, std::vector xs, + std::vector ys, + std::vector as, + std::vector bs, Tx const& x, + double sig2) { // extend out first and last lines for convolution double sig = std::sqrt(sig2); - xs[0] += - 10 * sig; - xs[n-1] += 10 * sig; + xs[0] += -10 * sig; + xs[n - 1] += 10 * sig; - // no need to convolve far from center of gaussian, so + // no need to convolve far from center of gaussian, so // get lower and upper indexes for integration bounds auto lb = upper_bound(xs.begin(), xs.end(), x - 10 * sig); int ind_start = distance(xs.begin(), lb) - 1; @@ -123,32 +117,29 @@ double min_diff(int n, std::vector xs) { return dmin; } -/* +/* given a set of pairs (x_i, y_i), do a gaussian interpolation through those points and evaluate the interpolation at the points xs_new */ template -inline vector gaus_interp(int n, - std::vector xs, - std::vector ys, - int n_new, - std::vector xs_new) { - +inline vector gaus_interp(int n, std::vector xs, + std::vector ys, int n_new, + std::vector xs_new) { // find minimum distance between points for std of gaussian kernel double sig2 = square(min_diff(n, xs) * SIG2_SCALE); - // copy ys into a new vector + // copy ys into a new vector std::vector y2s; y2s.resize(n); for (int i = 0; i < n; i++) { y2s[i] = ys[i]; } - // interatively find interpolation that coincides with ys at xs + // interatively find interpolation that coincides with ys at xs int max_iters = 50; double dmax, dd; std::vector as, bs; - for (int j=0; j < max_iters; j++) { + for (int j = 0; j < max_iters; j++) { // linear interpolation for new ys lin_interp_coefs(n, xs, y2s, as, bs); dmax = 0; @@ -157,9 +148,10 @@ inline vector gaus_interp(int n, y2s[i] += dd; dmax = std::max(std::abs(dd), dmax); } - if (dmax < INTERP_TOL) break; + if (dmax < INTERP_TOL) + break; } - + // fill ys_new with interpolated values at new_xs std::vector ys_new; for (int i = 0; i < n_new; i++) { diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp index fa116f01ffd..fc130268ff4 100644 --- a/stan/math/rev/fun/conv_gaus_line.hpp +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -10,31 +10,32 @@ namespace math { namespace internal { class conv_gaus_line_vari : public op_v_vari { - double t0_; - double t1_; - double a_; - double b_; - double sig2_; + double t0_; + double t1_; + double a_; + double b_; + double sig2_; public: explicit conv_gaus_line_vari(double t0, double t1, double a, double b, - vari* avi, double sig2) : - t0_(t0), - t1_(t1), - a_(a), - b_(b), - sig2_(sig2), - op_v_vari(conv_gaus_line(t0, t1, a, b, avi->val_, sig2), avi) {} - - void chain() { avi_->adj_ += adj_ * der_conv_gaus_line(t0_, t1_, a_, b_, - avi_->val_, sig2_); } + vari* avi, double sig2) + : t0_(t0), + t1_(t1), + a_(a), + b_(b), + sig2_(sig2), + op_v_vari(conv_gaus_line(t0, t1, a, b, avi->val_, sig2), avi) {} + + void chain() { + avi_->adj_ + += adj_ * der_conv_gaus_line(t0_, t1_, a_, b_, avi_->val_, sig2_); + } }; } // namespace internal - -inline var conv_gaus_line(double t0, double t1, double a, double b, - const var& x, double sig2) { +inline var conv_gaus_line(double t0, double t1, double a, double b, + const var& x, double sig2) { return var(new internal::conv_gaus_line_vari(t0, t1, a, b, x.vi_, sig2)); } diff --git a/test/unit/math/mix/fun/conv_gaus_line_test.cpp b/test/unit/math/mix/fun/conv_gaus_line_test.cpp index 58edebf9799..7124bd2cf2c 100644 --- a/test/unit/math/mix/fun/conv_gaus_line_test.cpp +++ b/test/unit/math/mix/fun/conv_gaus_line_test.cpp @@ -23,7 +23,7 @@ TEST(mathMixGausInterp, gaus_conv_der) { x0p = x0 - h; yn = conv_gaus_line(t0, t1, a, b, x0n, sig2); yp = conv_gaus_line(t0, t1, a, b, x0p, sig2); - dder = (yn - yp) / (2*h); + dder = (yn - yp) / (2 * h); // derivative with autodiff using stan::math::var; diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index 82dda3902d7..bb3e3fc8748 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -2,10 +2,9 @@ #include #include - TEST(mathMixGausInterp, derivs) { - using stan::math::var; using stan::math::gaus_interp; + using stan::math::var; std::vector xs, ys, x2s, y2s, ts, as, bs, y3s, t3s; double xmin, xmax, x, y, x2, y2, t, t0, t1, dd, dder, dder2; double x0, x1, y0, y1, tmp; @@ -29,7 +28,7 @@ TEST(mathMixGausInterp, derivs) { t0 = xs[1] - 1e-4; t0 = xs[0]; t1 = xmax; - for (int i=0; i < n_interp; i++) { + for (int i = 0; i < n_interp; i++) { t = t0 + i * (t1 - t0) / (n_interp - 1); xs_new.push_back(t); stan::math::var tvar = t; @@ -45,7 +44,7 @@ TEST(mathMixGausInterp, derivs) { std::vector ys_new_dder; std::vector ys_new_v2; - for (int i=0; i < n_interp; i++) { + for (int i = 0; i < n_interp; i++) { ys_new_v[i].grad(); dder = xs_new_v[i].adj(); ys_new_dder.push_back(dder); @@ -55,7 +54,7 @@ TEST(mathMixGausInterp, derivs) { double h = 1e-8; xs_new_p = xs_new; xs_new_n = xs_new; - for (int i=0; i < n_interp; i++) { + for (int i = 0; i < n_interp; i++) { xs_new_p[i] += -h; xs_new_n[i] += h; ys_new_p = gaus_interp(n, xs, ys, n_interp, xs_new_p); diff --git a/test/unit/math/prim/fun/conv_gaus_line_test.cpp b/test/unit/math/prim/fun/conv_gaus_line_test.cpp index 5ee590ddba7..2847308cee8 100644 --- a/test/unit/math/prim/fun/conv_gaus_line_test.cpp +++ b/test/unit/math/prim/fun/conv_gaus_line_test.cpp @@ -6,7 +6,7 @@ evaluate the convolution using trapezoid rule */ double check_int_dumb(double t0, double t1, double a, double b, double x0, - double sig2) { + double sig2) { double x, y, t, tot; vector ys, ts; double pi = stan::math::pi(); diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index 47ef627f9e4..231a862d97c 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -21,7 +21,7 @@ TEST(mathPrimGausInterp, interp_line) { double t0 = xmin; double t1 = xmax; double t; - for (int i=0; i < n_interp; i++) { + for (int i = 0; i < n_interp; i++) { t = t0 + i * (t1 - t0) / (n_interp - 1); xs_new.push_back(t); } @@ -32,7 +32,7 @@ TEST(mathPrimGausInterp, interp_line) { // test points double tmp, y; - for (int i=0; i < n_interp; i++) { + for (int i = 0; i < n_interp; i++) { tmp = (ys[1] - ys[0]) / (xs[1] - xs[0]); y = tmp * xs_new[i] + ys[0] - tmp * xs[0]; ASSERT_NEAR(ys_new[i], y, ABS_TOL); @@ -50,10 +50,10 @@ TEST(mathPrimGausInterp, gaus_and_lin_interp) { double xmax = 1; double x; vector xs, ys; - for (int i=0; i < n; i++) { + for (int i = 0; i < n; i++) { x = xmin + i * (xmax - xmin) / (n - 1); xs.push_back(x); - ys.push_back(x*x); + ys.push_back(x * x); } // vector of pts at which to compute interpolation @@ -62,7 +62,7 @@ TEST(mathPrimGausInterp, gaus_and_lin_interp) { double t0 = xmin; double t1 = xmax; double t; - for (int i=0; i < n_interp; i++) { + for (int i = 0; i < n_interp; i++) { t = t0 + i * (t1 - t0) / (n_interp - 1); xs_new.push_back(t); } @@ -74,7 +74,7 @@ TEST(mathPrimGausInterp, gaus_and_lin_interp) { // test points double tmp, y; - for (int i=0; i < n_interp; i++) { + for (int i = 0; i < n_interp; i++) { ASSERT_NEAR(ys_new_lin[i], ys_new_gaus[i], 1e-4); } } From efde0255fd9e07336b466901bd815776ac572071 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 1 Apr 2020 11:29:05 -0400 Subject: [PATCH 04/57] fixing messed up branch --- stan/math/prim/fun.hpp | 2 ++ stan/math/rev/fun.hpp | 1 + 2 files changed, 3 insertions(+) diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index f214acd63ea..40e68e5c9de 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -42,6 +42,7 @@ #include #include #include +#include #include #include #include @@ -97,6 +98,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index a6025090f21..98c555d0e4f 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include From 2169b76db152297d0f9f3ae5c1370f49b4ddd36b Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 1 Apr 2020 11:30:54 -0400 Subject: [PATCH 05/57] fixing messed up branch again --- stan/math/prim/fun/conv_gaus_line.hpp | 67 ++++++ stan/math/prim/fun/gaus_interp.hpp | 191 ++++++++++++++++++ stan/math/rev/fun/conv_gaus_line.hpp | 45 +++++ .../unit/math/mix/fun/conv_gaus_line_test.cpp | 35 ++++ test/unit/math/mix/fun/gaus_interp_test.cpp | 65 ++++++ .../math/prim/fun/conv_gaus_line_test.cpp | 53 +++++ test/unit/math/prim/fun/gaus_interp_test.cpp | 80 ++++++++ 7 files changed, 536 insertions(+) create mode 100644 stan/math/prim/fun/conv_gaus_line.hpp create mode 100644 stan/math/prim/fun/gaus_interp.hpp create mode 100644 stan/math/rev/fun/conv_gaus_line.hpp create mode 100644 test/unit/math/mix/fun/conv_gaus_line_test.cpp create mode 100644 test/unit/math/mix/fun/gaus_interp_test.cpp create mode 100644 test/unit/math/prim/fun/conv_gaus_line_test.cpp create mode 100644 test/unit/math/prim/fun/gaus_interp_test.cpp diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp new file mode 100644 index 00000000000..2f1f7d42220 --- /dev/null +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -0,0 +1,67 @@ +#ifndef STAN_MATH_PRIM_FUN_CONV_GAUS_LINE +#define STAN_MATH_PRIM_FUN_CONV_GAUS_LINE + +#include +#include +#include +#include + +namespace stan { +namespace math { + +/* +evaluate the derivative of conv_gaus_line with respect to x +*/ +double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, + double sig2) { + using stan::math::normal_cdf; + double pi = stan::math::pi(); + using std::exp; + using std::pow; + using std::sqrt; + double sig = sqrt(sig2); + double y; + + double alpha = sqrt(2 * pi * sig2); + y = (a * x0 + b) / alpha + * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) + + exp(-pow(t0 - x0, 2) / (2 * sig2))); + y += a * (normal_cdf(t1, x0, sig) - normal_cdf(t0, x0, sig)); + y -= a / alpha + * ((t1 - x0) * exp(-pow(t1 - x0, 2) / (2 * sig2)) + - (t0 - x0) * exp(-pow(t0 - x0, 2) / (2 * sig2))); + return y; +} + +/* +evaluate the integral + + | t1 +(2*pi*sig2)**-(1/2) | (a*t + b) * exp(-(t - x0)^2 / (2*sig2)) dt + | t0 + +*/ +template +double conv_gaus_line(double t0, double t1, double a, double b, Tx const& x, + double sig2) { + using stan::math::normal_cdf; + double pi = stan::math::pi(); + using std::exp; + using std::pow; + using std::sqrt; + double sig = sqrt(sig2); + + double y + = (a * value_of(x) + b) + * (normal_cdf(t1, value_of(x), sig) - normal_cdf(t0, value_of(x), sig)); + y += -a * sig2 / sqrt(2 * pi * sig2) + * (exp(-pow(t1 - value_of(x), 2) / (2 * sig2)) + - exp(-pow(t0 - value_of(x), 2) / (2 * sig2))); + + return y; +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp new file mode 100644 index 00000000000..16973e4d7f6 --- /dev/null +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -0,0 +1,191 @@ +#ifndef STAN_MATH_PRIM_FUN_GAUS_INTERP +#define STAN_MATH_PRIM_FUN_GAUS_INTERP + +#include +#include +#include +#include + +namespace stan { +namespace math { + +// set tolerance for how close interpolated values need to be to +// the specified function values +const double INTERP_TOL = 1e-8; +const double SIG2_SCALE = 0.1; + +/* +given a set of points (x_i, y_i) with x_i's in increasing order, find +a_i, b_i such that the line between (x_i, y_i) and (x_{i+1}, y_{i+1}) is +f(t) = a_i*t + b_i, store the a_i's and b_i's in as and bs. +*/ +void lin_interp_coefs(int n, std::vector xs, std::vector ys, + std::vector& as, std::vector& bs) { + // initialize size of vectors of linear coefficients + as.resize(n - 1); + bs.resize(n - 1); + + // find slope and intercept between each point + for (int i = 0; i < n - 1; i++) { + as[i] = (ys[i + 1] - ys[i]) / (xs[i + 1] - xs[i]); + bs[i] = -xs[i] * as[i] + ys[i]; + } +} + +/* +linear interpolation at one point, x, given the coefficients +*/ +<<<<<<< HEAD +double lin_interp_pt(int n, + std::vector xs, + std::vector ys, + std::vector as, + std::vector bs, + double x) { +======= +double lin_interp_pt(int n, vector xs, vector ys, + vector as, vector bs, double x) { +>>>>>>> 98fdf8bbec8eb8a07df7a99ac5b94dde87bb9654 + // find interval where x lives + if (x <= xs[0]) + return ys[0]; + if (x >= xs[n - 1]) + return ys[n - 1]; + + auto lb = upper_bound(xs.begin(), xs.end(), x); + int ind = distance(xs.begin(), lb) - 1; + ind = std::max(0, ind); + + return as[ind] * x + bs[ind]; +} + +/* +linear interpolation at a vector of points +*/ +<<<<<<< HEAD +std::vector lin_interp(int n, + std::vector xs, + std::vector ys, + int n_new, + std::vector xs_new) { +======= +vector lin_interp(int n, std::vector xs, std::vector ys, + int n_new, std::vector xs_new) { +>>>>>>> 98fdf8bbec8eb8a07df7a99ac5b94dde87bb9654 + // compute coefficients of linear interpolation + std::vector as, bs; + lin_interp_coefs(n, xs, ys, as, bs); + + // evaluate at new points + std::vector ys_new; + for (int i = 0; i < n_new; i++) { + ys_new.push_back(lin_interp_pt(n, xs, ys, as, bs, xs_new[i])); + } + return ys_new; +} + +/* +given a set of points and a width for the gaussian kernel, do a convolution +and evaluate at one point, x +*/ +template +inline return_type_t interp_gaus_pt(int n, std::vector xs, + std::vector ys, + std::vector as, + std::vector bs, Tx const& x, + double sig2) { + // extend out first and last lines for convolution + double sig = std::sqrt(sig2); + xs[0] += -10 * sig; + xs[n - 1] += 10 * sig; + + // no need to convolve far from center of gaussian, so + // get lower and upper indexes for integration bounds + auto lb = upper_bound(xs.begin(), xs.end(), x - 10 * sig); + int ind_start = distance(xs.begin(), lb) - 1; + ind_start = std::max(0, ind_start); + + auto ub = upper_bound(xs.begin(), xs.end(), x + 10 * sig); + int ind_end = distance(xs.begin(), ub); + ind_end = std::min(n - 1, ind_end); + + // sum convolutions over intervals + using return_t = return_type_t; + return_t y = 0; + for (int i = ind_start; i < ind_end; i++) { + y += conv_gaus_line(xs[i], xs[i + 1], as[i], bs[i], x, sig2); + } + + return y; +} + +/* +find the smallest difference between successive elements in a sorted vector +*/ +template +double min_diff(int n, std::vector xs) { + double dmin = value_of(xs[1]) - value_of(xs[0]); + for (int i = 1; i < n - 1; i++) { + if (value_of(xs[i + 1]) - value_of(xs[i]) < dmin) { + dmin = value_of(xs[i + 1]) - value_of(xs[i]); + } + } + return dmin; +} + +/* +given a set of pairs (x_i, y_i), do a gaussian interpolation through those +points and evaluate the interpolation at the points xs_new +*/ +template +<<<<<<< HEAD +inline std::vector gaus_interp(int n, + std::vector xs, + std::vector ys, + int n_new, + std::vector xs_new) { + +======= +inline vector gaus_interp(int n, std::vector xs, + std::vector ys, int n_new, + std::vector xs_new) { +>>>>>>> 98fdf8bbec8eb8a07df7a99ac5b94dde87bb9654 + // find minimum distance between points for std of gaussian kernel + double sig2 = square(min_diff(n, xs) * SIG2_SCALE); + + // copy ys into a new vector + std::vector y2s; + y2s.resize(n); + for (int i = 0; i < n; i++) { + y2s[i] = ys[i]; + } + + // interatively find interpolation that coincides with ys at xs + int max_iters = 50; + double dmax, dd; + std::vector as, bs; + for (int j = 0; j < max_iters; j++) { + // linear interpolation for new ys + lin_interp_coefs(n, xs, y2s, as, bs); + dmax = 0; + for (int i = 0; i < n; i++) { + dd = ys[i] - interp_gaus_pt(n, xs, y2s, as, bs, xs[i], sig2); + y2s[i] += dd; + dmax = std::max(std::abs(dd), dmax); + } + if (dmax < INTERP_TOL) + break; + } + + // fill ys_new with interpolated values at new_xs + std::vector ys_new; + for (int i = 0; i < n_new; i++) { + ys_new.push_back(interp_gaus_pt(n, xs, y2s, as, bs, xs_new[i], sig2)); + } + return ys_new; +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp new file mode 100644 index 00000000000..fc130268ff4 --- /dev/null +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -0,0 +1,45 @@ +#ifndef STAN_MATH_REV_FUN_CONV_GAUS_LINE +#define STAN_MATH_REV_FUN_CONV_GAUS_LINE + +#include +#include +#include + +namespace stan { +namespace math { + +namespace internal { +class conv_gaus_line_vari : public op_v_vari { + double t0_; + double t1_; + double a_; + double b_; + double sig2_; + + public: + explicit conv_gaus_line_vari(double t0, double t1, double a, double b, + vari* avi, double sig2) + : t0_(t0), + t1_(t1), + a_(a), + b_(b), + sig2_(sig2), + op_v_vari(conv_gaus_line(t0, t1, a, b, avi->val_, sig2), avi) {} + + void chain() { + avi_->adj_ + += adj_ * der_conv_gaus_line(t0_, t1_, a_, b_, avi_->val_, sig2_); + } +}; + +} // namespace internal + +inline var conv_gaus_line(double t0, double t1, double a, double b, + const var& x, double sig2) { + return var(new internal::conv_gaus_line_vari(t0, t1, a, b, x.vi_, sig2)); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/mix/fun/conv_gaus_line_test.cpp b/test/unit/math/mix/fun/conv_gaus_line_test.cpp new file mode 100644 index 00000000000..7124bd2cf2c --- /dev/null +++ b/test/unit/math/mix/fun/conv_gaus_line_test.cpp @@ -0,0 +1,35 @@ +#include +#include +#include + +/* +check derivative using finite differences and autodiffing +*/ +TEST(mathMixGausInterp, gaus_conv_der) { + double t0, t1, a, b, x0, sig2, y, h, yn, yp, x0p, x0n, dder, dder2; + double alpha, sig; + using stan::math::conv_gaus_line; + + t0 = 0; + t1 = 5; + a = -1; + b = 2; + x0 = 0.99; + sig2 = 2; + + // derivative with finite difference + h = 1e-4; + x0n = x0 + h; + x0p = x0 - h; + yn = conv_gaus_line(t0, t1, a, b, x0n, sig2); + yp = conv_gaus_line(t0, t1, a, b, x0p, sig2); + dder = (yn - yp) / (2 * h); + + // derivative with autodiff + using stan::math::var; + var v = x0; + var vy = conv_gaus_line(t0, t1, a, b, v, sig2); + vy.grad(); + + ASSERT_NEAR(dder, v.adj(), 1e-5); +} diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp new file mode 100644 index 00000000000..86b7bf53bfc --- /dev/null +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -0,0 +1,65 @@ +#include +#include +#include + +TEST(mathMixGausInterp, derivs) { + using stan::math::gaus_interp; + using stan::math::var; + std::vector xs, ys, x2s, y2s, ts, as, bs, y3s, t3s; + double xmin, xmax, x, y, x2, y2, t, t0, t1, dd, dder, dder2; + double x0, x1, y0, y1, tmp; + int n; + + // generate function tabulation + n = 10; + xmin = 0; + xmax = 1; + for (int i = 0; i < n; i++) { + x = xmin + i * (xmax - xmin) / (n - 1); + xs.push_back(x); + y = rand_r(0) % 100; + ys.push_back(y); + } + + // create vector of interpolation pts + std::vector xs_new_v; + std::vector xs_new; + int n_interp = 100; + t0 = xs[1] - 1e-4; + t0 = xs[0]; + t1 = xmax; + for (int i = 0; i < n_interp; i++) { + t = t0 + i * (t1 - t0) / (n_interp - 1); + xs_new.push_back(t); + stan::math::var tvar = t; + xs_new_v.push_back(tvar); + } + + std::vector ys_new, ys_new_p, ys_new_n, xs_new_p, xs_new_n; + ys_new = gaus_interp(n, xs, ys, n_interp, xs_new); + + // autodiff at each interpolation pt + std::vector ys_new_v; + ys_new_v = gaus_interp(n, xs, ys, n_interp, xs_new_v); + + std::vector ys_new_dder; + std::vector ys_new_v2; + for (int i = 0; i < n_interp; i++) { + ys_new_v[i].grad(); + dder = xs_new_v[i].adj(); + ys_new_dder.push_back(dder); + } + + // take derivative of interpolation using finite differencing + double h = 1e-8; + xs_new_p = xs_new; + xs_new_n = xs_new; + for (int i = 0; i < n_interp; i++) { + xs_new_p[i] += -h; + xs_new_n[i] += h; + ys_new_p = gaus_interp(n, xs, ys, n_interp, xs_new_p); + ys_new_n = gaus_interp(n, xs, ys, n_interp, xs_new_n); + dder2 = (ys_new_n[i] - ys_new_p[i]) / (2 * h); + ASSERT_NEAR(ys_new_dder[i], dder2, 1e-5); + } +} diff --git a/test/unit/math/prim/fun/conv_gaus_line_test.cpp b/test/unit/math/prim/fun/conv_gaus_line_test.cpp new file mode 100644 index 00000000000..2847308cee8 --- /dev/null +++ b/test/unit/math/prim/fun/conv_gaus_line_test.cpp @@ -0,0 +1,53 @@ +#include +#include +#include + +/* +evaluate the convolution using trapezoid rule +*/ +double check_int_dumb(double t0, double t1, double a, double b, double x0, + double sig2) { + double x, y, t, tot; + vector ys, ts; + double pi = stan::math::pi(); + int n; + + // evaluate function at equispaced nodes + n = 100000; + for (int i = 0; i < n; i++) { + t = t0 + i * (t1 - t0) / (n - 1); + ts.push_back(t); + + y = a * t + b; + y = y * exp(-pow(t - x0, 2) / (2 * sig2)); + ys.push_back(y); + } + + // sum up function evaluations + tot = 0; + for (int i = 0; i < n; i++) { + tot += ys[i]; + } + tot *= (t1 - t0) / n; + tot /= sqrt(2 * pi * sig2); + + return tot; +} + +/* +check that integral using formula is close to trapezoid approximation +*/ +TEST(mathMixConvGausLin, trap_test) { + double t0, t1, a, b, x0, sig2, y, y2; + using stan::math::conv_gaus_line; + t0 = 0; + t1 = 1; + a = 1; + b = 0; + x0 = 1; + sig2 = 3; + y = check_int_dumb(t0, t1, a, b, x0, sig2); + y2 = conv_gaus_line(t0, t1, a, b, x0, sig2); + + ASSERT_NEAR(y, y2, 1e-5); +} diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp new file mode 100644 index 00000000000..231a862d97c --- /dev/null +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -0,0 +1,80 @@ +#include +#include +#include + +double ABS_TOL = 1e-12; + +TEST(mathPrimGausInterp, interp_line) { + using stan::math::gaus_interp; + + // check that interpolation of line returns the same function + // generate function tabulation + int n = 2; + double xmin = 0; + double xmax = 1; + vector xs = {0, 1}; + vector ys = {0, 2}; + + // vector of pts at which to compute interpolation + vector xs_new; + int n_interp = 100; + double t0 = xmin; + double t1 = xmax; + double t; + for (int i = 0; i < n_interp; i++) { + t = t0 + i * (t1 - t0) / (n_interp - 1); + xs_new.push_back(t); + } + + // create interpolation + vector ys_new; + ys_new = gaus_interp(n, xs, ys, n_interp, xs_new); + + // test points + double tmp, y; + for (int i = 0; i < n_interp; i++) { + tmp = (ys[1] - ys[0]) / (xs[1] - xs[0]); + y = tmp * xs_new[i] + ys[0] - tmp * xs[0]; + ASSERT_NEAR(ys_new[i], y, ABS_TOL); + } +} + +TEST(mathPrimGausInterp, gaus_and_lin_interp) { + using stan::math::gaus_interp; + using stan::math::lin_interp; + + // check that interpolation of line returns the same function + // generate function tabulation + int n = 30; + double xmin = 0; + double xmax = 1; + double x; + vector xs, ys; + for (int i = 0; i < n; i++) { + x = xmin + i * (xmax - xmin) / (n - 1); + xs.push_back(x); + ys.push_back(x * x); + } + + // vector of pts at which to compute interpolation + vector xs_new; + int n_interp = 100; + double t0 = xmin; + double t1 = xmax; + double t; + for (int i = 0; i < n_interp; i++) { + t = t0 + i * (t1 - t0) / (n_interp - 1); + xs_new.push_back(t); + } + + // create interpolation + vector ys_new_gaus, ys_new_lin; + ys_new_gaus = gaus_interp(n, xs, ys, n_interp, xs_new); + ys_new_lin = lin_interp(n, xs, ys, n_interp, xs_new); + + // test points + double tmp, y; + for (int i = 0; i < n_interp; i++) { + ASSERT_NEAR(ys_new_lin[i], ys_new_gaus[i], 1e-4); + } +} From 6c8726de277c04ab80b05d5067c9457a55a65666 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 1 Apr 2020 16:42:09 +0000 Subject: [PATCH 06/57] [Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (tags/RELEASE_500/final) --- stan/math/opencl/kernel_generator/load.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/opencl/kernel_generator/load.hpp b/stan/math/opencl/kernel_generator/load.hpp index ae0f7d1eb6d..799cdb551a0 100644 --- a/stan/math/opencl/kernel_generator/load.hpp +++ b/stan/math/opencl/kernel_generator/load.hpp @@ -49,7 +49,7 @@ class load_ * Creates a deep copy of this expression. * @return copy of \c *this */ - inline load_ deep_copy() const& { return load_(a_); } + inline load_ deep_copy() const & { return load_(a_); } inline load_ deep_copy() && { return load_(std::forward(a_)); } /** From ee28a652ee8d42548ac6320c2fa861a2ac254fdf Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Tue, 4 Aug 2020 17:20:20 -0400 Subject: [PATCH 07/57] changes to reflect design doc discussion --- stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/conv_gaus_line.hpp | 45 ++++-- stan/math/prim/fun/gaus_interp.hpp | 151 ++++++++++--------- test/unit/math/prim/fun/gaus_interp_test.cpp | 31 +++- 4 files changed, 135 insertions(+), 93 deletions(-) diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 40e68e5c9de..b05ca82421d 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -149,6 +149,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 2f1f7d42220..e283be3cb7b 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -1,17 +1,30 @@ #ifndef STAN_MATH_PRIM_FUN_CONV_GAUS_LINE #define STAN_MATH_PRIM_FUN_CONV_GAUS_LINE -#include -#include #include #include +#include +#include namespace stan { namespace math { -/* -evaluate the derivative of conv_gaus_line with respect to x -*/ + +/** \ingroup prob_dists + * Evaluates the derivative of the convolution of a line with a Gaussian + * kernel. + * + * \f$\frac{\partial}{\partial x} \int_{t_0}^{t_1} (at + b) \f$ + * \f$e^{\frac{-(t-x)^2}{2\sigma^2} dt \f$ + * + * @param t_0 lower integration bound + * @param t_1 upper integration bound + * @param a coefficient of t in line + * @param b constant in line + * @param x0 point at which convolution is evaluated + * @param sig2 variance of the Gaussian kernel + * @return The value of the derivative + */ double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, double sig2) { using stan::math::normal_cdf; @@ -33,15 +46,20 @@ double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, return y; } -/* -evaluate the integral - - | t1 -(2*pi*sig2)**-(1/2) | (a*t + b) * exp(-(t - x0)^2 / (2*sig2)) dt - | t0 -*/ -template +/** \ingroup prob_dists + * Evaluate the convolution of a line with a Gaussian kernel. + * + * \f$\int_{t_0}^{t_1} (at + b) e^{\frac{-(t-x)^2}{2\sigma^2} dt \f$ + * + * @param t_0 lower integration bound + * @param t_1 upper integration bound + * @param a coefficient of t in line + * @param b constant in line + * @param x0 point at which convolution is evaluated + * @param sig2 variance of the Gaussian kernel + * @return The value of the derivative + */ double conv_gaus_line(double t0, double t1, double a, double b, Tx const& x, double sig2) { using stan::math::normal_cdf; @@ -57,7 +75,6 @@ double conv_gaus_line(double t0, double t1, double a, double b, Tx const& x, y += -a * sig2 / sqrt(2 * pi * sig2) * (exp(-pow(t1 - value_of(x), 2) / (2 * sig2)) - exp(-pow(t0 - value_of(x), 2) / (2 * sig2))); - return y; } diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 10d16f841af..f402f2d2bf2 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -1,96 +1,109 @@ #ifndef STAN_MATH_PRIM_FUN_GAUS_INTERP #define STAN_MATH_PRIM_FUN_GAUS_INTERP -#include #include #include #include +#include using namespace std; namespace stan { namespace math { +// struct for passing precomputation data +struct gaus_interp_params { + std::vector as; + std::vector bs; + double sig2; +}; + +namespace internal { + +// struct for passing precomputation data +struct asbs { + std::vector as; + std::vector bs; +}; + // set tolerance for how close interpolated values need to be to // the specified function values const double INTERP_TOL = 1e-8; const double SIG2_SCALE = 0.1; +const double NSTDS = 10; /* given a set of points (x_i, y_i) with x_i's in increasing order, find a_i, b_i such that the line between (x_i, y_i) and (x_{i+1}, y_{i+1}) is f(t) = a_i*t + b_i, store the a_i's and b_i's in as and bs. */ -void lin_interp_coefs(int n, std::vector xs, std::vector ys, - std::vector& as, std::vector& bs) { +asbs lin_interp_coefs(std::vector xs, std::vector ys) { + + int n = xs.size(); + // initialize size of vectors of linear coefficients - as.resize(n - 1); - bs.resize(n - 1); + std::vector as(n - 1); + std::vector bs(n - 1); // find slope and intercept between each point for (int i = 0; i < n - 1; i++) { as[i] = (ys[i + 1] - ys[i]) / (xs[i + 1] - xs[i]); bs[i] = -xs[i] * as[i] + ys[i]; } + // return struct with as and bs + asbs coefs; + coefs.as = as; + coefs.bs = bs; + return coefs; } /* -linear interpolation at one point, x, given the coefficients -*/ -double lin_interp_pt(int n, vector xs, vector ys, - vector as, vector bs, double x) { - // find interval where x lives - if (x <= xs[0]) - return ys[0]; - if (x >= xs[n - 1]) - return ys[n - 1]; - - auto lb = upper_bound(xs.begin(), xs.end(), x); - int ind = distance(xs.begin(), lb) - 1; - ind = std::max(0, ind); - - return as[ind] * x + bs[ind]; -} - -/* -linear interpolation at a vector of points +find the smallest difference between successive elements in a sorted vector */ -vector lin_interp(int n, std::vector xs, std::vector ys, - int n_new, std::vector xs_new) { - // compute coefficients of linear interpolation - vector as, bs; - lin_interp_coefs(n, xs, ys, as, bs); - - // evaluate at new points - vector ys_new; - for (int i = 0; i < n_new; i++) { - ys_new.push_back(lin_interp_pt(n, xs, ys, as, bs, xs_new[i])); +template +double min_diff(int n, std::vector xs) { + double dmin = value_of(xs[1]) - value_of(xs[0]); + for (int i = 1; i < n - 1; i++) { + if (value_of(xs[i + 1]) - value_of(xs[i]) < dmin) { + dmin = value_of(xs[i + 1]) - value_of(xs[i]); + } } - return ys_new; + return dmin; } +} // namespace internal + /* given a set of points and a width for the gaussian kernel, do a convolution and evaluate at one point, x */ template -inline return_type_t interp_gaus_pt(int n, std::vector xs, - std::vector ys, - std::vector as, - std::vector bs, Tx const& x, - double sig2) { +inline return_type_t gaus_interp(std::vector xs, + std::vector ys, + gaus_interp_params params, + Tx const& x) { + using internal::NSTDS; + + // enforce that interpolation point is between smallest and largest + // reference point + static char const* function = "interp_gaus_pt"; + check_less_or_equal(function, "Interpolation point", x, xs.back()); + check_greater_or_equal(function, "Interpolation point", x, xs.front()); + + int n = xs.size(); + // extend out first and last lines for convolution - double sig = std::sqrt(sig2); - xs[0] += -10 * sig; - xs[n - 1] += 10 * sig; + double sig = std::sqrt(params.sig2); + xs[0] += -NSTDS * sig; + xs[n - 1] += NSTDS * sig; // no need to convolve far from center of gaussian, so // get lower and upper indexes for integration bounds - auto lb = upper_bound(xs.begin(), xs.end(), x - 10 * sig); + auto lb = upper_bound(xs.begin(), xs.end(), x - NSTDS * sig); int ind_start = distance(xs.begin(), lb) - 1; ind_start = std::max(0, ind_start); - auto ub = upper_bound(xs.begin(), xs.end(), x + 10 * sig); + auto ub = upper_bound(xs.begin(), xs.end(), x + NSTDS * sig); int ind_end = distance(xs.begin(), ub); ind_end = std::min(n - 1, ind_end); @@ -98,36 +111,32 @@ inline return_type_t interp_gaus_pt(int n, std::vector xs, using return_t = return_type_t; return_t y = 0; for (int i = ind_start; i < ind_end; i++) { - y += conv_gaus_line(xs[i], xs[i + 1], as[i], bs[i], x, sig2); + y += conv_gaus_line(xs[i], xs[i+1], params.as[i], params.bs[i], x, + params.sig2); } return y; } -/* -find the smallest difference between successive elements in a sorted vector -*/ -template -double min_diff(int n, std::vector xs) { - double dmin = value_of(xs[1]) - value_of(xs[0]); - for (int i = 1; i < n - 1; i++) { - if (value_of(xs[i + 1]) - value_of(xs[i]) < dmin) { - dmin = value_of(xs[i + 1]) - value_of(xs[i]); - } - } - return dmin; -} - /* given a set of pairs (x_i, y_i), do a gaussian interpolation through those points and evaluate the interpolation at the points xs_new */ template -inline vector gaus_interp(int n, std::vector xs, - std::vector ys, int n_new, - std::vector xs_new) { +gaus_interp_params gaus_interp_precomp(std::vector xs, + std::vector ys, + std::vector xs_new) { + using internal::min_diff; + using internal::lin_interp_coefs; + using internal::SIG2_SCALE; + using internal::INTERP_TOL; + gaus_interp_params params; + internal::asbs coefs; + int n = xs.size(); + int n_new = xs_new.size(); + // find minimum distance between points for std of gaussian kernel - double sig2 = square(min_diff(n, xs) * SIG2_SCALE); + params.sig2 = square(min_diff(n, xs) * SIG2_SCALE); // copy ys into a new vector std::vector y2s; @@ -139,13 +148,14 @@ inline vector gaus_interp(int n, std::vector xs, // interatively find interpolation that coincides with ys at xs int max_iters = 50; double dmax, dd; - std::vector as, bs; for (int j = 0; j < max_iters; j++) { // linear interpolation for new ys - lin_interp_coefs(n, xs, y2s, as, bs); + coefs = lin_interp_coefs(xs, y2s); + params.as = coefs.as; + params.bs = coefs.bs; dmax = 0; for (int i = 0; i < n; i++) { - dd = ys[i] - interp_gaus_pt(n, xs, y2s, as, bs, xs[i], sig2); + dd = ys[i] - gaus_interp(xs, y2s, params, xs[i]); y2s[i] += dd; dmax = std::max(std::abs(dd), dmax); } @@ -153,12 +163,7 @@ inline vector gaus_interp(int n, std::vector xs, break; } - // fill ys_new with interpolated values at new_xs - std::vector ys_new; - for (int i = 0; i < n_new; i++) { - ys_new.push_back(interp_gaus_pt(n, xs, y2s, as, bs, xs_new[i], sig2)); - } - return ys_new; + return params; } } // namespace math diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index 231a862d97c..44497799281 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -1,11 +1,14 @@ #include #include #include +////#include double ABS_TOL = 1e-12; TEST(mathPrimGausInterp, interp_line) { using stan::math::gaus_interp; + using stan::math::gaus_interp_params; + using stan::math::gaus_interp_precomp; // check that interpolation of line returns the same function // generate function tabulation @@ -27,20 +30,25 @@ TEST(mathPrimGausInterp, interp_line) { } // create interpolation - vector ys_new; - ys_new = gaus_interp(n, xs, ys, n_interp, xs_new); + vector ys_new_gaus(n_interp); + gaus_interp_params params = gaus_interp_precomp(xs, ys, xs_new); + for (int i=0; i ys_new_gaus, ys_new_lin; - ys_new_gaus = gaus_interp(n, xs, ys, n_interp, xs_new); - ys_new_lin = lin_interp(n, xs, ys, n_interp, xs_new); + vector ys_new_gaus(n_interp); + vector ys_new_lin(n_interp); + + // linear interpolation + ys_new_lin.resize(n_interp); + for (int i=0; i Date: Tue, 4 Aug 2020 17:52:53 -0400 Subject: [PATCH 08/57] fixing tests --- stan/math/prim/fun/conv_gaus_line.hpp | 22 ++++++++++---------- stan/math/prim/fun/gaus_interp.hpp | 19 +++++++++++++---- test/unit/math/mix/fun/gaus_interp_test.cpp | 16 +++++++------- test/unit/math/prim/fun/gaus_interp_test.cpp | 11 +++++++--- 4 files changed, 43 insertions(+), 25 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index e283be3cb7b..3bbcc976276 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -12,10 +12,10 @@ namespace math { /** \ingroup prob_dists * Evaluates the derivative of the convolution of a line with a Gaussian - * kernel. + * kernel on an interval. * * \f$\frac{\partial}{\partial x} \int_{t_0}^{t_1} (at + b) \f$ - * \f$e^{\frac{-(t-x)^2}{2\sigma^2} dt \f$ + * \f$ e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ * * @param t_0 lower integration bound * @param t_1 upper integration bound @@ -28,15 +28,14 @@ namespace math { double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, double sig2) { using stan::math::normal_cdf; - double pi = stan::math::pi(); using std::exp; using std::pow; using std::sqrt; - double sig = sqrt(sig2); - double y; + const double pi = stan::math::pi(); + const double sig = sqrt(sig2); + const double alpha = sqrt(2 * pi * sig2); - double alpha = sqrt(2 * pi * sig2); - y = (a * x0 + b) / alpha + double y = (a * x0 + b) / alpha * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) + exp(-pow(t0 - x0, 2) / (2 * sig2))); y += a * (normal_cdf(t1, x0, sig) - normal_cdf(t0, x0, sig)); @@ -48,9 +47,9 @@ double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, /** \ingroup prob_dists - * Evaluate the convolution of a line with a Gaussian kernel. + * Evaluate the convolution of a line with a Gaussian kernel on an interval. * - * \f$\int_{t_0}^{t_1} (at + b) e^{\frac{-(t-x)^2}{2\sigma^2} dt \f$ + * \f$\int_{t_0}^{t_1} (at + b) e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ * * @param t_0 lower integration bound * @param t_1 upper integration bound @@ -60,14 +59,15 @@ double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, * @param sig2 variance of the Gaussian kernel * @return The value of the derivative */ +template double conv_gaus_line(double t0, double t1, double a, double b, Tx const& x, double sig2) { using stan::math::normal_cdf; - double pi = stan::math::pi(); using std::exp; using std::pow; using std::sqrt; - double sig = sqrt(sig2); + const double pi = stan::math::pi(); + const double sig = sqrt(sig2); double y = (a * value_of(x) + b) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index f402f2d2bf2..e6b4f419c8c 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -122,10 +122,8 @@ inline return_type_t gaus_interp(std::vector xs, given a set of pairs (x_i, y_i), do a gaussian interpolation through those points and evaluate the interpolation at the points xs_new */ -template gaus_interp_params gaus_interp_precomp(std::vector xs, - std::vector ys, - std::vector xs_new) { + std::vector ys) { using internal::min_diff; using internal::lin_interp_coefs; using internal::SIG2_SCALE; @@ -133,7 +131,6 @@ gaus_interp_params gaus_interp_precomp(std::vector xs, gaus_interp_params params; internal::asbs coefs; int n = xs.size(); - int n_new = xs_new.size(); // find minimum distance between points for std of gaussian kernel params.sig2 = square(min_diff(n, xs) * SIG2_SCALE); @@ -165,6 +162,20 @@ gaus_interp_params gaus_interp_precomp(std::vector xs, return params; } +template +inline vector gaus_interp_vect(std::vector xs, + std::vector ys, + std::vector xs_new) { + int n_interp = xs_new.size(); + std::vector ys_new(n_interp); + + // create interpolation + gaus_interp_params params = gaus_interp_precomp(xs, ys); + for (int i=0; i TEST(mathMixGausInterp, derivs) { - using stan::math::gaus_interp; + using stan::math::gaus_interp_vect; using stan::math::var; std::vector xs, ys, x2s, y2s, ts, as, bs, y3s, t3s; double xmin, xmax, x, y, x2, y2, t, t0, t1, dd, dder, dder2; @@ -26,8 +26,10 @@ TEST(mathMixGausInterp, derivs) { std::vector xs_new_v; std::vector xs_new; int n_interp = 100; - t0 = xs[0]; - t1 = xmax; + + // add a cushion to each endpoint so that we don't leave interval + t0 = xs[0] + 0.01; + t1 = xmax - 0.01; for (int i = 0; i < n_interp; i++) { t = t0 + i * (t1 - t0) / (n_interp - 1); xs_new.push_back(t); @@ -36,11 +38,11 @@ TEST(mathMixGausInterp, derivs) { } std::vector ys_new, ys_new_p, ys_new_n, xs_new_p, xs_new_n; - ys_new = gaus_interp(n, xs, ys, n_interp, xs_new); + ys_new = gaus_interp_vect(xs, ys, xs_new); // autodiff at each interpolation pt std::vector ys_new_v; - ys_new_v = gaus_interp(n, xs, ys, n_interp, xs_new_v); + ys_new_v = gaus_interp_vect(xs, ys, xs_new_v); std::vector ys_new_dder; std::vector ys_new_v2; @@ -57,8 +59,8 @@ TEST(mathMixGausInterp, derivs) { for (int i = 0; i < n_interp; i++) { xs_new_p[i] += -h; xs_new_n[i] += h; - ys_new_p = gaus_interp(n, xs, ys, n_interp, xs_new_p); - ys_new_n = gaus_interp(n, xs, ys, n_interp, xs_new_n); + ys_new_p = gaus_interp_vect(xs, ys, xs_new_p); + ys_new_n = gaus_interp_vect(xs, ys, xs_new_n); dder2 = (ys_new_n[i] - ys_new_p[i]) / (2 * h); ASSERT_NEAR(ys_new_dder[i], dder2, 1e-5); } diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index 44497799281..27871bd32d9 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -9,6 +9,7 @@ TEST(mathPrimGausInterp, interp_line) { using stan::math::gaus_interp; using stan::math::gaus_interp_params; using stan::math::gaus_interp_precomp; + using stan::math::gaus_interp_vect; // check that interpolation of line returns the same function // generate function tabulation @@ -29,19 +30,23 @@ TEST(mathPrimGausInterp, interp_line) { xs_new.push_back(t); } - // create interpolation + // create interpolation using precomp vector ys_new_gaus(n_interp); - gaus_interp_params params = gaus_interp_precomp(xs, ys, xs_new); + gaus_interp_params params = gaus_interp_precomp(xs, ys); for (int i=0; i ys_new_gaus2 = gaus_interp_vect(xs, ys, xs_new); + // test points double tmp, y; for (int i = 0; i < n_interp; i++) { tmp = (ys[1] - ys[0]) / (xs[1] - xs[0]); y = tmp * xs_new[i] + ys[0] - tmp * xs[0]; ASSERT_NEAR(ys_new_gaus[i], y, ABS_TOL); + ASSERT_NEAR(ys_new_gaus2[i], y, ABS_TOL); } } @@ -86,7 +91,7 @@ TEST(mathPrimGausInterp, gaus_and_lin_interp) { } // gaus interpolation - gaus_interp_params params = gaus_interp_precomp(xs, ys, xs_new); + gaus_interp_params params = gaus_interp_precomp(xs, ys); for (int i=0; i Date: Tue, 4 Aug 2020 18:26:40 -0400 Subject: [PATCH 09/57] reflects steve's changes --- stan/math/prim/fun/gaus_interp.hpp | 45 +++++++++++++++----------- stan/math/prim/fun/lin_interp.hpp | 47 ++++++++++++++++++++++++++++ stan/math/rev/fun/conv_gaus_line.hpp | 38 ++++++++++++++++++++++ 3 files changed, 111 insertions(+), 19 deletions(-) create mode 100644 stan/math/prim/fun/lin_interp.hpp diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index e6b4f419c8c..69294579fd2 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -26,18 +26,13 @@ struct asbs { std::vector bs; }; -// set tolerance for how close interpolated values need to be to -// the specified function values -const double INTERP_TOL = 1e-8; -const double SIG2_SCALE = 0.1; -const double NSTDS = 10; - /* given a set of points (x_i, y_i) with x_i's in increasing order, find a_i, b_i such that the line between (x_i, y_i) and (x_{i+1}, y_{i+1}) is f(t) = a_i*t + b_i, store the a_i's and b_i's in as and bs. */ -asbs lin_interp_coefs(std::vector xs, std::vector ys) { +asbs lin_interp_coefs(std::vector const& xs, + std::vector const& ys) { int n = xs.size(); @@ -61,7 +56,7 @@ asbs lin_interp_coefs(std::vector xs, std::vector ys) { find the smallest difference between successive elements in a sorted vector */ template -double min_diff(int n, std::vector xs) { +double min_diff(int n, std::vector const& xs) { double dmin = value_of(xs[1]) - value_of(xs[0]); for (int i = 1; i < n - 1; i++) { if (value_of(xs[i + 1]) - value_of(xs[i]) < dmin) { @@ -78,11 +73,11 @@ given a set of points and a width for the gaussian kernel, do a convolution and evaluate at one point, x */ template -inline return_type_t gaus_interp(std::vector xs, - std::vector ys, +inline return_type_t gaus_interp(std::vector const& xs, + std::vector const& ys, gaus_interp_params params, Tx const& x) { - using internal::NSTDS; + const double NSTDS = 10; // enforce that interpolation point is between smallest and largest // reference point @@ -92,10 +87,13 @@ inline return_type_t gaus_interp(std::vector xs, int n = xs.size(); + // create copy of xs so that endpoints can be extended + std::vector xs2 = xs; + // extend out first and last lines for convolution double sig = std::sqrt(params.sig2); - xs[0] += -NSTDS * sig; - xs[n - 1] += NSTDS * sig; + xs2[0] = xs[0] - NSTDS * sig; + xs2[n - 1] = xs[n-1] + NSTDS * sig; // no need to convolve far from center of gaussian, so // get lower and upper indexes for integration bounds @@ -111,7 +109,7 @@ inline return_type_t gaus_interp(std::vector xs, using return_t = return_type_t; return_t y = 0; for (int i = ind_start; i < ind_end; i++) { - y += conv_gaus_line(xs[i], xs[i+1], params.as[i], params.bs[i], x, + y += conv_gaus_line(xs2[i], xs2[i+1], params.as[i], params.bs[i], x, params.sig2); } @@ -126,14 +124,16 @@ gaus_interp_params gaus_interp_precomp(std::vector xs, std::vector ys) { using internal::min_diff; using internal::lin_interp_coefs; - using internal::SIG2_SCALE; - using internal::INTERP_TOL; gaus_interp_params params; internal::asbs coefs; + const double INTERP_TOL = 1e-8; + const double SIG2_SCALE = 0.1; int n = xs.size(); // find minimum distance between points for std of gaussian kernel params.sig2 = square(min_diff(n, xs) * SIG2_SCALE); + params.as.resize(n-1); + params.bs.resize(n-1); // copy ys into a new vector std::vector y2s; @@ -147,9 +147,16 @@ gaus_interp_params gaus_interp_precomp(std::vector xs, double dmax, dd; for (int j = 0; j < max_iters; j++) { // linear interpolation for new ys - coefs = lin_interp_coefs(xs, y2s); - params.as = coefs.as; - params.bs = coefs.bs; + //coefs = lin_interp_coefs(xs, y2s); + //params.as = coefs.as; + //params.bs = coefs.bs; + + // find slope and intercept of line between each point + for (int i = 0; i < n - 1; i++) { + params.as[i] = (y2s[i + 1] - y2s[i]) / (xs[i + 1] - xs[i]); + params.bs[i] = -xs[i] * params.as[i] + y2s[i]; + } + dmax = 0; for (int i = 0; i < n; i++) { dd = ys[i] - gaus_interp(xs, y2s, params, xs[i]); diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/lin_interp.hpp new file mode 100644 index 00000000000..e0f1b22b784 --- /dev/null +++ b/stan/math/prim/fun/lin_interp.hpp @@ -0,0 +1,47 @@ +#ifndef STAN_MATH_PRIM_FUN_LIN_INTERP +#define STAN_MATH_PRIM_FUN_LIN_INTERP + +#include +#include +#include + +using namespace std; + +namespace stan { +namespace math { + +double lin_interp(std::vector xs, std::vector ys, double x) { + double x1, x2, y1, y2; + int n; + + n = xs.size(); + + // if x is less than left endpoint or greater than right, return endpoint + if (x <= xs[0]) { + return ys[0]; + } + if (x >= xs[n - 1]) { + return ys[n - 1]; + } + + // find in between which points the input, x, lives + auto ub = upper_bound(xs.begin(), xs.end(), x); + int ind = distance(xs.begin(), ub); + + // check if the interpolation point falls on a reference point + if (x == *(ub-1)) return ys[ind-1]; + + // do linear interpolation + x1 = *(ub-1); + x2 = *ub; + y1 = ys[ind-1]; + y2 = ys[ind]; + + return y1 + (x - x1) * (y2 - y1) / (x2-x1); +} + + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp index fc130268ff4..3e8afc20c56 100644 --- a/stan/math/rev/fun/conv_gaus_line.hpp +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -9,6 +9,44 @@ namespace stan { namespace math { namespace internal { + +/** \ingroup prob_dists + * Evaluates the derivative of the convolution of a line with a Gaussian + * kernel on an interval. + * + * \f$\frac{\partial}{\partial x} \int_{t_0}^{t_1} (at + b) \f$ + * \f$ e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ + * + * @param t_0 lower integration bound + * @param t_1 upper integration bound + * @param a coefficient of t in line + * @param b constant in line + * @param x0 point at which convolution is evaluated + * @param sig2 variance of the Gaussian kernel + * @return The value of the derivative + */ +double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, + double sig2) { + using stan::math::normal_cdf; + using std::exp; + using std::pow; + using std::sqrt; + const double pi = stan::math::pi(); + const double sig = sqrt(sig2); + const double alpha = sqrt(2 * pi * sig2); + + double y = (a * x0 + b) / alpha + * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) + + exp(-pow(t0 - x0, 2) / (2 * sig2))); + y += a * (normal_cdf(t1, x0, sig) - normal_cdf(t0, x0, sig)); + y -= a / alpha + * ((t1 - x0) * exp(-pow(t1 - x0, 2) / (2 * sig2)) + - (t0 - x0) * exp(-pow(t0 - x0, 2) / (2 * sig2))); + return y; +} + + + class conv_gaus_line_vari : public op_v_vari { double t0_; double t1_; From eb0a0d9e71354be1ac2b901b4d1dfa2454464767 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 5 Aug 2020 12:20:24 -0400 Subject: [PATCH 10/57] adding docs --- stan/math/prim/fun/conv_gaus_line.hpp | 39 +-------- stan/math/prim/fun/gaus_interp.hpp | 113 ++++++++++++++------------ stan/math/prim/fun/lin_interp.hpp | 30 +++++-- 3 files changed, 84 insertions(+), 98 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 3bbcc976276..27c5d14c01d 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -9,44 +9,7 @@ namespace stan { namespace math { - -/** \ingroup prob_dists - * Evaluates the derivative of the convolution of a line with a Gaussian - * kernel on an interval. - * - * \f$\frac{\partial}{\partial x} \int_{t_0}^{t_1} (at + b) \f$ - * \f$ e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ - * - * @param t_0 lower integration bound - * @param t_1 upper integration bound - * @param a coefficient of t in line - * @param b constant in line - * @param x0 point at which convolution is evaluated - * @param sig2 variance of the Gaussian kernel - * @return The value of the derivative - */ -double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, - double sig2) { - using stan::math::normal_cdf; - using std::exp; - using std::pow; - using std::sqrt; - const double pi = stan::math::pi(); - const double sig = sqrt(sig2); - const double alpha = sqrt(2 * pi * sig2); - - double y = (a * x0 + b) / alpha - * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) - + exp(-pow(t0 - x0, 2) / (2 * sig2))); - y += a * (normal_cdf(t1, x0, sig) - normal_cdf(t0, x0, sig)); - y -= a / alpha - * ((t1 - x0) * exp(-pow(t1 - x0, 2) / (2 * sig2)) - - (t0 - x0) * exp(-pow(t0 - x0, 2) / (2 * sig2))); - return y; -} - - -/** \ingroup prob_dists +/** * Evaluate the convolution of a line with a Gaussian kernel on an interval. * * \f$\int_{t_0}^{t_1} (at + b) e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 69294579fd2..eecc777bf64 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -26,32 +26,6 @@ struct asbs { std::vector bs; }; -/* -given a set of points (x_i, y_i) with x_i's in increasing order, find -a_i, b_i such that the line between (x_i, y_i) and (x_{i+1}, y_{i+1}) is -f(t) = a_i*t + b_i, store the a_i's and b_i's in as and bs. -*/ -asbs lin_interp_coefs(std::vector const& xs, - std::vector const& ys) { - - int n = xs.size(); - - // initialize size of vectors of linear coefficients - std::vector as(n - 1); - std::vector bs(n - 1); - - // find slope and intercept between each point - for (int i = 0; i < n - 1; i++) { - as[i] = (ys[i + 1] - ys[i]) / (xs[i + 1] - xs[i]); - bs[i] = -xs[i] * as[i] + ys[i]; - } - // return struct with as and bs - asbs coefs; - coefs.as = as; - coefs.bs = bs; - return coefs; -} - /* find the smallest difference between successive elements in a sorted vector */ @@ -68,14 +42,30 @@ double min_diff(int n, std::vector const& xs) { } // namespace internal -/* -given a set of points and a width for the gaussian kernel, do a convolution -and evaluate at one point, x -*/ +/** + * Given a set of reference points \f$(xs_i, ys_i)\f$, create a mollifier + * that intersects the reference points. This function requires as input + * a struct created by the function gaus_interp_precomp. The algorithm + * used to create the mollifier is an iterative algorithm that works + * as follows. First a linear + * interpolation is created through the reference points. Then, the + * linear interpolation is convolved with a Gaussian whose width is + * proportional the smallest distance between successive points + * \f$xs_i\f$ and \f$xs_{i+1}\f$. Since the convolution is unlikely to + * intersect the reference points, the y-values of the reference points + * are shifted and the process repeats itself until the interpolation + * intersects all reference points. + * + * @param xs vector of independent variable of reference points + * @param ys vector of dependent variable of reference points + * @param params a struct created by gaus_interp_precomp that + * @param x the point at which to evaluate the interpolation + * @return value of the interpolation at x + */ template inline return_type_t gaus_interp(std::vector const& xs, std::vector const& ys, - gaus_interp_params params, + gaus_interp_params const& params, Tx const& x) { const double NSTDS = 10; @@ -116,14 +106,28 @@ inline return_type_t gaus_interp(std::vector const& xs, return y; } -/* -given a set of pairs (x_i, y_i), do a gaussian interpolation through those -points and evaluate the interpolation at the points xs_new -*/ -gaus_interp_params gaus_interp_precomp(std::vector xs, - std::vector ys) { +/** + * This function was written to be used with gaus_interp. This function + * computes the shifted y-values of the reference points of an interpolation + * in such a way that when that piecewise linear function is convolved + * with a Gaussian kernel, the resulting function coincides with the + * points \f$(xs_i, ys_i)\f$ inputted into this function. The output of this + * function depends heavily on the choice of width of the Gaussian + * kernel, which at the time of writing, is set to one tenth the + * minimum distance between successive elements of the vector xs. + * A tolerance for the maximum distance between the interpolation and + * all reference points is also set manually and is not an input. + * + * + * @param xs vector of independent variable of reference points + * @param ys vector of dependent variable of reference points + * @param params a struct created by gaus_interp_precomp that + * @param x the point at which to evaluate the interpolation + * @return struct containing slopes, intercepts, and width of kernel + */ +gaus_interp_params gaus_interp_precomp(std::vector const& xs, + std::vector const& ys) { using internal::min_diff; - using internal::lin_interp_coefs; gaus_interp_params params; internal::asbs coefs; const double INTERP_TOL = 1e-8; @@ -135,22 +139,13 @@ gaus_interp_params gaus_interp_precomp(std::vector xs, params.as.resize(n-1); params.bs.resize(n-1); - // copy ys into a new vector - std::vector y2s; - y2s.resize(n); - for (int i = 0; i < n; i++) { - y2s[i] = ys[i]; - } + // copy ys into a new vector that will be changed + std::vector y2s = ys; // interatively find interpolation that coincides with ys at xs int max_iters = 50; double dmax, dd; for (int j = 0; j < max_iters; j++) { - // linear interpolation for new ys - //coefs = lin_interp_coefs(xs, y2s); - //params.as = coefs.as; - //params.bs = coefs.bs; - // find slope and intercept of line between each point for (int i = 0; i < n - 1; i++) { params.as[i] = (y2s[i + 1] - y2s[i]) / (xs[i + 1] - xs[i]); @@ -169,10 +164,24 @@ gaus_interp_params gaus_interp_precomp(std::vector xs, return params; } + +/** + * This function combines gaus_interp_precomp and gaus_interp. + * It takes as input two vectors of reference points (xs and ys) + * in addition to a vector, xs_new, of points at which the + * function will evaluate the interpolation through those reference + * points. + * + * @param xs vector of independent variable of reference points + * @param ys vector of dependent variable of reference points + * @param xs_new vector of point at which to evaluate interpolation + * @return vector of interpolation values + */ + template -inline vector gaus_interp_vect(std::vector xs, - std::vector ys, - std::vector xs_new) { +inline vector gaus_interp_vect(std::vector const& xs, + std::vector const& ys, + std::vector const& xs_new) { int n_interp = xs_new.size(); std::vector ys_new(n_interp); diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/lin_interp.hpp index e0f1b22b784..daef8135818 100644 --- a/stan/math/prim/fun/lin_interp.hpp +++ b/stan/math/prim/fun/lin_interp.hpp @@ -10,11 +10,25 @@ using namespace std; namespace stan { namespace math { -double lin_interp(std::vector xs, std::vector ys, double x) { - double x1, x2, y1, y2; - int n; - n = xs.size(); +/** + * This function performs linear interpolation. The function takes as + * input two vectors of reference points \f$(xs_i, ys_i)\f$ in addition + * to one value, \f$x\f$, at which the value of the linear interpolation + * is to be returned. The vector xs is required to be in increasing order + * For values of \f$x\f$ less than xs[0], the function returns ys[0]. + * For values of \f$x\f$ greater than xs[n-1], the function returns + * ys[n-1], where xs is of length n. + * + * @param xs vector of independent variable of reference points + * @param ys vector of dependent variable of reference points + * @param x the point at which to evaluate the interpolation + * @return value of linear interpolation at x + */ +double lin_interp(std::vector const& xs, + std::vector const& ys, + double x) { + int n = xs.size(); // if x is less than left endpoint or greater than right, return endpoint if (x <= xs[0]) { @@ -32,10 +46,10 @@ double lin_interp(std::vector xs, std::vector ys, double x) { if (x == *(ub-1)) return ys[ind-1]; // do linear interpolation - x1 = *(ub-1); - x2 = *ub; - y1 = ys[ind-1]; - y2 = ys[ind]; + double x1 = *(ub-1); + double x2 = *ub; + double y1 = ys[ind-1]; + double y2 = ys[ind]; return y1 + (x - x1) * (y2 - y1) / (x2-x1); } From c4582be5522d564679ff9bf5d093b598dc5d8b59 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 5 Aug 2020 16:46:28 +0000 Subject: [PATCH 11/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/conv_gaus_line.hpp | 4 +- stan/math/prim/fun/gaus_interp.hpp | 68 ++++++++++---------- stan/math/prim/fun/lin_interp.hpp | 22 +++---- stan/math/rev/fun/conv_gaus_line.hpp | 8 +-- test/unit/math/prim/fun/gaus_interp_test.cpp | 8 +-- 5 files changed, 53 insertions(+), 57 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 27c5d14c01d..66e3e8a0c45 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -9,7 +9,7 @@ namespace stan { namespace math { -/** +/** * Evaluate the convolution of a line with a Gaussian kernel on an interval. * * \f$\int_{t_0}^{t_1} (at + b) e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ @@ -18,7 +18,7 @@ namespace math { * @param t_1 upper integration bound * @param a coefficient of t in line * @param b constant in line - * @param x0 point at which convolution is evaluated + * @param x0 point at which convolution is evaluated * @param sig2 variance of the Gaussian kernel * @return The value of the derivative */ diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index eecc777bf64..195760cabaa 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -40,36 +40,36 @@ double min_diff(int n, std::vector const& xs) { return dmin; } -} // namespace internal +} // namespace internal /** * Given a set of reference points \f$(xs_i, ys_i)\f$, create a mollifier * that intersects the reference points. This function requires as input * a struct created by the function gaus_interp_precomp. The algorithm - * used to create the mollifier is an iterative algorithm that works - * as follows. First a linear - * interpolation is created through the reference points. Then, the - * linear interpolation is convolved with a Gaussian whose width is - * proportional the smallest distance between successive points - * \f$xs_i\f$ and \f$xs_{i+1}\f$. Since the convolution is unlikely to + * used to create the mollifier is an iterative algorithm that works + * as follows. First a linear + * interpolation is created through the reference points. Then, the + * linear interpolation is convolved with a Gaussian whose width is + * proportional the smallest distance between successive points + * \f$xs_i\f$ and \f$xs_{i+1}\f$. Since the convolution is unlikely to * intersect the reference points, the y-values of the reference points - * are shifted and the process repeats itself until the interpolation - * intersects all reference points. + * are shifted and the process repeats itself until the interpolation + * intersects all reference points. * * @param xs vector of independent variable of reference points * @param ys vector of dependent variable of reference points - * @param params a struct created by gaus_interp_precomp that + * @param params a struct created by gaus_interp_precomp that * @param x the point at which to evaluate the interpolation * @return value of the interpolation at x */ template inline return_type_t gaus_interp(std::vector const& xs, - std::vector const& ys, - gaus_interp_params const& params, - Tx const& x) { + std::vector const& ys, + gaus_interp_params const& params, + Tx const& x) { const double NSTDS = 10; - // enforce that interpolation point is between smallest and largest + // enforce that interpolation point is between smallest and largest // reference point static char const* function = "interp_gaus_pt"; check_less_or_equal(function, "Interpolation point", x, xs.back()); @@ -83,7 +83,7 @@ inline return_type_t gaus_interp(std::vector const& xs, // extend out first and last lines for convolution double sig = std::sqrt(params.sig2); xs2[0] = xs[0] - NSTDS * sig; - xs2[n - 1] = xs[n-1] + NSTDS * sig; + xs2[n - 1] = xs[n - 1] + NSTDS * sig; // no need to convolve far from center of gaussian, so // get lower and upper indexes for integration bounds @@ -99,8 +99,8 @@ inline return_type_t gaus_interp(std::vector const& xs, using return_t = return_type_t; return_t y = 0; for (int i = ind_start; i < ind_end; i++) { - y += conv_gaus_line(xs2[i], xs2[i+1], params.as[i], params.bs[i], x, - params.sig2); + y += conv_gaus_line(xs2[i], xs2[i + 1], params.as[i], params.bs[i], x, + params.sig2); } return y; @@ -109,24 +109,24 @@ inline return_type_t gaus_interp(std::vector const& xs, /** * This function was written to be used with gaus_interp. This function * computes the shifted y-values of the reference points of an interpolation - * in such a way that when that piecewise linear function is convolved - * with a Gaussian kernel, the resulting function coincides with the + * in such a way that when that piecewise linear function is convolved + * with a Gaussian kernel, the resulting function coincides with the * points \f$(xs_i, ys_i)\f$ inputted into this function. The output of this - * function depends heavily on the choice of width of the Gaussian - * kernel, which at the time of writing, is set to one tenth the - * minimum distance between successive elements of the vector xs. - * A tolerance for the maximum distance between the interpolation and - * all reference points is also set manually and is not an input. + * function depends heavily on the choice of width of the Gaussian + * kernel, which at the time of writing, is set to one tenth the + * minimum distance between successive elements of the vector xs. + * A tolerance for the maximum distance between the interpolation and + * all reference points is also set manually and is not an input. * * * @param xs vector of independent variable of reference points * @param ys vector of dependent variable of reference points - * @param params a struct created by gaus_interp_precomp that + * @param params a struct created by gaus_interp_precomp that * @param x the point at which to evaluate the interpolation - * @return struct containing slopes, intercepts, and width of kernel + * @return struct containing slopes, intercepts, and width of kernel */ gaus_interp_params gaus_interp_precomp(std::vector const& xs, - std::vector const& ys) { + std::vector const& ys) { using internal::min_diff; gaus_interp_params params; internal::asbs coefs; @@ -136,8 +136,8 @@ gaus_interp_params gaus_interp_precomp(std::vector const& xs, // find minimum distance between points for std of gaussian kernel params.sig2 = square(min_diff(n, xs) * SIG2_SCALE); - params.as.resize(n-1); - params.bs.resize(n-1); + params.as.resize(n - 1); + params.bs.resize(n - 1); // copy ys into a new vector that will be changed std::vector y2s = ys; @@ -168,8 +168,8 @@ gaus_interp_params gaus_interp_precomp(std::vector const& xs, /** * This function combines gaus_interp_precomp and gaus_interp. * It takes as input two vectors of reference points (xs and ys) - * in addition to a vector, xs_new, of points at which the - * function will evaluate the interpolation through those reference + * in addition to a vector, xs_new, of points at which the + * function will evaluate the interpolation through those reference * points. * * @param xs vector of independent variable of reference points @@ -180,14 +180,14 @@ gaus_interp_params gaus_interp_precomp(std::vector const& xs, template inline vector gaus_interp_vect(std::vector const& xs, - std::vector const& ys, - std::vector const& xs_new) { + std::vector const& ys, + std::vector const& xs_new) { int n_interp = xs_new.size(); std::vector ys_new(n_interp); // create interpolation gaus_interp_params params = gaus_interp_precomp(xs, ys); - for (int i=0; i const& xs, - std::vector const& ys, - double x) { +double lin_interp(std::vector const& xs, std::vector const& ys, + double x) { int n = xs.size(); // if x is less than left endpoint or greater than right, return endpoint @@ -43,17 +41,17 @@ double lin_interp(std::vector const& xs, int ind = distance(xs.begin(), ub); // check if the interpolation point falls on a reference point - if (x == *(ub-1)) return ys[ind-1]; + if (x == *(ub - 1)) + return ys[ind - 1]; - // do linear interpolation - double x1 = *(ub-1); + // do linear interpolation + double x1 = *(ub - 1); double x2 = *ub; - double y1 = ys[ind-1]; + double y1 = ys[ind - 1]; double y2 = ys[ind]; - - return y1 + (x - x1) * (y2 - y1) / (x2-x1); -} + return y1 + (x - x1) * (y2 - y1) / (x2 - x1); +} } // namespace math } // namespace stan diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp index 3e8afc20c56..c661aa51f79 100644 --- a/stan/math/rev/fun/conv_gaus_line.hpp +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -21,7 +21,7 @@ namespace internal { * @param t_1 upper integration bound * @param a coefficient of t in line * @param b constant in line - * @param x0 point at which convolution is evaluated + * @param x0 point at which convolution is evaluated * @param sig2 variance of the Gaussian kernel * @return The value of the derivative */ @@ -36,8 +36,8 @@ double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, const double alpha = sqrt(2 * pi * sig2); double y = (a * x0 + b) / alpha - * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) - + exp(-pow(t0 - x0, 2) / (2 * sig2))); + * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) + + exp(-pow(t0 - x0, 2) / (2 * sig2))); y += a * (normal_cdf(t1, x0, sig) - normal_cdf(t0, x0, sig)); y -= a / alpha * ((t1 - x0) * exp(-pow(t1 - x0, 2) / (2 * sig2)) @@ -45,8 +45,6 @@ double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, return y; } - - class conv_gaus_line_vari : public op_v_vari { double t0_; double t1_; diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index 27871bd32d9..5642ed99daa 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -33,7 +33,7 @@ TEST(mathPrimGausInterp, interp_line) { // create interpolation using precomp vector ys_new_gaus(n_interp); gaus_interp_params params = gaus_interp_precomp(xs, ys); - for (int i=0; i Date: Wed, 5 Aug 2020 15:06:53 -0400 Subject: [PATCH 12/57] fixing formatting --- test/unit/math/prim/fun/gaus_interp_test.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index 5642ed99daa..330b742dffb 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -1,7 +1,6 @@ #include #include #include -////#include double ABS_TOL = 1e-12; From 0ce6663569011e8f827e753b7230e5f018433f77 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 5 Aug 2020 15:08:03 -0400 Subject: [PATCH 13/57] doc typo fix --- stan/math/prim/fun/conv_gaus_line.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 66e3e8a0c45..9600bba20d0 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -14,8 +14,8 @@ namespace math { * * \f$\int_{t_0}^{t_1} (at + b) e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ * - * @param t_0 lower integration bound - * @param t_1 upper integration bound + * @param t0 lower integration bound + * @param t1 upper integration bound * @param a coefficient of t in line * @param b constant in line * @param x0 point at which convolution is evaluated From c4dd75d10c41ee4c53d010dd308b7e4ca62eb69b Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 5 Aug 2020 16:10:56 -0400 Subject: [PATCH 14/57] another doc typo fix --- stan/math/rev/fun/conv_gaus_line.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp index c661aa51f79..bfb88be7b4b 100644 --- a/stan/math/rev/fun/conv_gaus_line.hpp +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -17,8 +17,8 @@ namespace internal { * \f$\frac{\partial}{\partial x} \int_{t_0}^{t_1} (at + b) \f$ * \f$ e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ * - * @param t_0 lower integration bound - * @param t_1 upper integration bound + * @param t0 lower integration bound + * @param t1 upper integration bound * @param a coefficient of t in line * @param b constant in line * @param x0 point at which convolution is evaluated From be423daa810ad838e67ac3affefa7e759a10882b Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 5 Aug 2020 16:24:53 -0400 Subject: [PATCH 15/57] another doc typo fix --- stan/math/prim/fun/conv_gaus_line.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 9600bba20d0..1221fd9ddc3 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -18,7 +18,7 @@ namespace math { * @param t1 upper integration bound * @param a coefficient of t in line * @param b constant in line - * @param x0 point at which convolution is evaluated + * @param x point at which convolution is evaluated * @param sig2 variance of the Gaussian kernel * @return The value of the derivative */ From 05737f4ff26271c3de23dd600a9e2bb9a865bf84 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 5 Aug 2020 19:20:30 -0400 Subject: [PATCH 16/57] another doc typo fix --- stan/math/prim/fun/gaus_interp.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 195760cabaa..8c9b2a779b5 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -121,8 +121,6 @@ inline return_type_t gaus_interp(std::vector const& xs, * * @param xs vector of independent variable of reference points * @param ys vector of dependent variable of reference points - * @param params a struct created by gaus_interp_precomp that - * @param x the point at which to evaluate the interpolation * @return struct containing slopes, intercepts, and width of kernel */ gaus_interp_params gaus_interp_precomp(std::vector const& xs, From 149f0d7dfa3a18a0a76122f2a5e1b3f37d6fb653 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 5 Aug 2020 19:59:42 -0400 Subject: [PATCH 17/57] adding checks and tests --- stan/math/prim/fun/gaus_interp.hpp | 21 +++++---- test/unit/math/prim/fun/gaus_interp_test.cpp | 48 ++++++++++++++++++++ 2 files changed, 59 insertions(+), 10 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 8c9b2a779b5..39970755a65 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -20,12 +20,6 @@ struct gaus_interp_params { namespace internal { -// struct for passing precomputation data -struct asbs { - std::vector as; - std::vector bs; -}; - /* find the smallest difference between successive elements in a sorted vector */ @@ -67,14 +61,17 @@ inline return_type_t gaus_interp(std::vector const& xs, std::vector const& ys, gaus_interp_params const& params, Tx const& x) { - const double NSTDS = 10; - // enforce that interpolation point is between smallest and largest // reference point - static char const* function = "interp_gaus_pt"; + static char const* function = "gaus_interp"; check_less_or_equal(function, "Interpolation point", x, xs.back()); check_greater_or_equal(function, "Interpolation point", x, xs.front()); + check_ordered(function, "xs", xs); + check_not_nan(function, "xs", xs); + check_not_nan(function, "ys", ys); + check_not_nan(function, "x", x); + const double NSTDS = 10; int n = xs.size(); // create copy of xs so that endpoints can be extended @@ -125,9 +122,13 @@ inline return_type_t gaus_interp(std::vector const& xs, */ gaus_interp_params gaus_interp_precomp(std::vector const& xs, std::vector const& ys) { + static char const* function = "gaus_interp_precomp"; + check_not_nan(function, "xs", xs); + check_not_nan(function, "ys", ys); + check_ordered(function, "xs", xs); + using internal::min_diff; gaus_interp_params params; - internal::asbs coefs; const double INTERP_TOL = 1e-8; const double SIG2_SCALE = 0.1; int n = xs.size(); diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index 330b742dffb..8e5adaadcc6 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -4,6 +4,54 @@ double ABS_TOL = 1e-12; +TEST(mathPrimGausInterp, throwing) { + using stan::math::gaus_interp; + using stan::math::gaus_interp_precomp; + using stan::math::gaus_interp_params; + double nan = std::numeric_limits::quiet_NaN(); + double x; + vector xs, ys; + + // check that when xs are not increasing, an error throws + int n = 2; + xs = {1, 1}; + ys = {0, 2}; + EXPECT_THROW(gaus_interp_precomp(xs, ys), std::domain_error); + + // check that when xs contain a nan + xs = {nan, 1}; + ys = {0, 2}; + EXPECT_THROW(gaus_interp_precomp(xs, ys), std::domain_error); + + // check that error throws when trying to interpolate out of range or nan + xs = {0, 1}; + ys = {0, 2}; + gaus_interp_params params = gaus_interp_precomp(xs, ys); + x = 1.1; + EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + x = -0.1; + EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + x = nan; + EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + + // ys can't contain nan + xs = {0, 1}; + ys = {0, nan}; + x = 0.5; + EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + + // xs can't contain nan + xs = {0, nan}; + ys = {0, 2}; + x = 0.5; + EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + + // xs must be increasing + xs = {1, 1}; + ys = {0, 2}; + EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); +} + TEST(mathPrimGausInterp, interp_line) { using stan::math::gaus_interp; using stan::math::gaus_interp_params; From 8594a4743ced9fb3700ff855e7c6c41a2133c979 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 6 Aug 2020 00:15:04 +0000 Subject: [PATCH 18/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- test/unit/math/prim/fun/gaus_interp_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index 8e5adaadcc6..67214c06650 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -6,8 +6,8 @@ double ABS_TOL = 1e-12; TEST(mathPrimGausInterp, throwing) { using stan::math::gaus_interp; - using stan::math::gaus_interp_precomp; using stan::math::gaus_interp_params; + using stan::math::gaus_interp_precomp; double nan = std::numeric_limits::quiet_NaN(); double x; vector xs, ys; From ee5b696228e2479fc91ba4fe83653f9b8c5af2e7 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Thu, 6 Aug 2020 08:00:42 -0400 Subject: [PATCH 19/57] fix namespace --- stan/math/prim/fun/gaus_interp.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 39970755a65..37bb9e75975 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -128,6 +128,7 @@ gaus_interp_params gaus_interp_precomp(std::vector const& xs, check_ordered(function, "xs", xs); using internal::min_diff; + using stan::math::square; gaus_interp_params params; const double INTERP_TOL = 1e-8; const double SIG2_SCALE = 0.1; From adbfc60fc5a6ba062f3f5a88ad1338584fc5928e Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Thu, 6 Aug 2020 08:28:29 -0400 Subject: [PATCH 20/57] fix namespace --- stan/math/prim/fun/gaus_interp.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 37bb9e75975..3af8bb8146e 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -128,14 +128,13 @@ gaus_interp_params gaus_interp_precomp(std::vector const& xs, check_ordered(function, "xs", xs); using internal::min_diff; - using stan::math::square; gaus_interp_params params; const double INTERP_TOL = 1e-8; const double SIG2_SCALE = 0.1; int n = xs.size(); // find minimum distance between points for std of gaussian kernel - params.sig2 = square(min_diff(n, xs) * SIG2_SCALE); + params.sig2 = stan::math::square(min_diff(n, xs) * SIG2_SCALE); params.as.resize(n - 1); params.bs.resize(n - 1); From 3fde18a7d01fe2b61ac13fc46d12b4e994a6fbab Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Thu, 6 Aug 2020 17:20:52 -0400 Subject: [PATCH 21/57] added test for lin_interp --- stan/math/prim/fun/conv_gaus_line.hpp | 2 +- stan/math/prim/fun/gaus_interp.hpp | 26 +++++++++++--------- stan/math/prim/fun/lin_interp.hpp | 10 +++++--- test/unit/math/prim/fun/gaus_interp_test.cpp | 10 ++++++++ 4 files changed, 32 insertions(+), 16 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 1221fd9ddc3..7f3d8bdb1ea 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -23,7 +23,7 @@ namespace math { * @return The value of the derivative */ template -double conv_gaus_line(double t0, double t1, double a, double b, Tx const& x, +double conv_gaus_line(double t0, double t1, double a, double b, const Tx& x, double sig2) { using stan::math::normal_cdf; using std::exp; diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 3af8bb8146e..72e8d88637c 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -5,6 +5,7 @@ #include #include #include +#include using namespace std; @@ -24,7 +25,7 @@ namespace internal { find the smallest difference between successive elements in a sorted vector */ template -double min_diff(int n, std::vector const& xs) { +double min_diff(int n, const std::vector& xs) { double dmin = value_of(xs[1]) - value_of(xs[0]); for (int i = 1; i < n - 1; i++) { if (value_of(xs[i + 1]) - value_of(xs[i]) < dmin) { @@ -57,10 +58,10 @@ double min_diff(int n, std::vector const& xs) { * @return value of the interpolation at x */ template -inline return_type_t gaus_interp(std::vector const& xs, - std::vector const& ys, - gaus_interp_params const& params, - Tx const& x) { +inline return_type_t gaus_interp(const std::vector& xs, + const std::vector& ys, + const gaus_interp_params& params, + const Tx& x) { // enforce that interpolation point is between smallest and largest // reference point static char const* function = "gaus_interp"; @@ -70,6 +71,7 @@ inline return_type_t gaus_interp(std::vector const& xs, check_not_nan(function, "xs", xs); check_not_nan(function, "ys", ys); check_not_nan(function, "x", x); + check_greater(function, "xs", xs.size(), 1); const double NSTDS = 10; int n = xs.size(); @@ -120,12 +122,13 @@ inline return_type_t gaus_interp(std::vector const& xs, * @param ys vector of dependent variable of reference points * @return struct containing slopes, intercepts, and width of kernel */ -gaus_interp_params gaus_interp_precomp(std::vector const& xs, - std::vector const& ys) { +gaus_interp_params gaus_interp_precomp(const std::vector& xs, + const std::vector& ys) { static char const* function = "gaus_interp_precomp"; check_not_nan(function, "xs", xs); check_not_nan(function, "ys", ys); check_ordered(function, "xs", xs); + check_greater(function, "xs", xs.size(), 1); using internal::min_diff; gaus_interp_params params; @@ -134,7 +137,7 @@ gaus_interp_params gaus_interp_precomp(std::vector const& xs, int n = xs.size(); // find minimum distance between points for std of gaussian kernel - params.sig2 = stan::math::square(min_diff(n, xs) * SIG2_SCALE); + params.sig2 = square(min_diff(n, xs) * SIG2_SCALE); params.as.resize(n - 1); params.bs.resize(n - 1); @@ -160,7 +163,6 @@ gaus_interp_params gaus_interp_precomp(std::vector const& xs, if (dmax < INTERP_TOL) break; } - return params; } @@ -178,9 +180,9 @@ gaus_interp_params gaus_interp_precomp(std::vector const& xs, */ template -inline vector gaus_interp_vect(std::vector const& xs, - std::vector const& ys, - std::vector const& xs_new) { +inline vector gaus_interp_vect(const std::vector& xs, + const std::vector& ys, + const std::vector& xs_new) { int n_interp = xs_new.size(); std::vector ys_new(n_interp); diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/lin_interp.hpp index 4622743e090..22a6c4503a4 100644 --- a/stan/math/prim/fun/lin_interp.hpp +++ b/stan/math/prim/fun/lin_interp.hpp @@ -5,8 +5,6 @@ #include #include -using namespace std; - namespace stan { namespace math { @@ -24,8 +22,14 @@ namespace math { * @param x the point at which to evaluate the interpolation * @return value of linear interpolation at x */ -double lin_interp(std::vector const& xs, std::vector const& ys, +double lin_interp(const std::vector& xs, const std::vector& ys, double x) { + static char const* function = "lin_interp"; + check_ordered(function, "xs", xs); + check_not_nan(function, "xs", xs); + check_not_nan(function, "ys", ys); + check_not_nan(function, "x", x); + check_greater(function, "xs", xs.size(), 1); int n = xs.size(); // if x is less than left endpoint or greater than right, return endpoint diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index 67214c06650..3f7b1427e7e 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -23,6 +23,11 @@ TEST(mathPrimGausInterp, throwing) { ys = {0, 2}; EXPECT_THROW(gaus_interp_precomp(xs, ys), std::domain_error); + // xs must contain at least two elements + xs = {1}; + ys = {0, 2}; + EXPECT_THROW(gaus_interp_precomp(xs, ys), std::domain_error); + // check that error throws when trying to interpolate out of range or nan xs = {0, 1}; ys = {0, 2}; @@ -50,6 +55,11 @@ TEST(mathPrimGausInterp, throwing) { xs = {1, 1}; ys = {0, 2}; EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + + // xs must contain at least two elements + xs = {1}; + ys = {0, 2}; + EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); } TEST(mathPrimGausInterp, interp_line) { From c624bd01fb6d040b2f5b2d009e335e4a728f9f5f Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Thu, 6 Aug 2020 17:21:16 -0400 Subject: [PATCH 22/57] added test for lin_interp --- test/unit/math/prim/fun/lin_interp_test.cpp | 87 +++++++++++++++++++++ 1 file changed, 87 insertions(+) create mode 100644 test/unit/math/prim/fun/lin_interp_test.cpp diff --git a/test/unit/math/prim/fun/lin_interp_test.cpp b/test/unit/math/prim/fun/lin_interp_test.cpp new file mode 100644 index 00000000000..d4576deabd8 --- /dev/null +++ b/test/unit/math/prim/fun/lin_interp_test.cpp @@ -0,0 +1,87 @@ +#include +#include +#include + +double ABS_TOL = 1e-12; + +TEST(mathPrimLinInterp, throwing) { + using stan::math::lin_interp; + double nan = std::numeric_limits::quiet_NaN(); + double x = 0.5; + vector xs, ys; + + // check that when xs are not increasing, an error throws + int n = 2; + xs = {1, 1}; + ys = {0, 2}; + EXPECT_THROW(lin_interp(xs, ys, x), std::domain_error); + + // check when xs contain a nan + xs = {nan, 1}; + ys = {0, 2}; + EXPECT_THROW(lin_interp(xs, ys, x), std::domain_error); + + // xs must contain at least two elements + xs = {1}; + ys = {0, 2}; + EXPECT_THROW(lin_interp(xs, ys, x), std::domain_error); + + // ys can't contain nan + xs = {0, 1}; + ys = {0, nan}; + x = 0.5; + EXPECT_THROW(lin_interp(xs, ys, x), std::domain_error); + + // x can't be nan + xs = {0, 1}; + ys = {0, 2}; + x = nan; + EXPECT_THROW(lin_interp(xs, ys, x), std::domain_error); +} + +TEST(mathPrimLinInterp, interp_line) { + using stan::math::lin_interp; + + // check that interpolation of line returns the same function + // generate function tabulation + int n = 2; + double xmin = 0; + double xmax = 1; + vector xs = {0, 1}; + vector ys = {0, 2}; + + // vector of pts at which to compute interpolation + vector xs_new; + int n_interp = 100; + double t0 = xmin; + double t1 = xmax; + double t; + for (int i = 0; i < n_interp; i++) { + t = t0 + i * (t1 - t0) / (n_interp - 1); + xs_new.push_back(t); + } + + // interpolate + vector ys_new(n_interp); + for (int i = 0; i < n_interp; i++) { + ys_new[i] = lin_interp(xs, ys, xs_new[i]); + } + + // test points + double tmp, y; + for (int i = 0; i < n_interp; i++) { + tmp = (ys[1] - ys[0]) / (xs[1] - xs[0]); + y = tmp * xs_new[i] + ys[0] - tmp * xs[0]; + ASSERT_NEAR(ys_new[i], y, ABS_TOL); + } + + // check values outside of range of reference points + ASSERT_NEAR(lin_interp(xs, ys, -1), ys[0], ABS_TOL); + ASSERT_NEAR(lin_interp(xs, ys, 100), ys[1], ABS_TOL); + + // xs with more than 2 elements + vector xs2 = {0, 1, 2, 3, 4, 5, 6}; + vector ys2 = {0, 2, 5, 2, 3, 2, 2}; + double x = 0.5; + ASSERT_NEAR(lin_interp(xs, ys, x), lin_interp(xs2, ys2, x), ABS_TOL); +} From 77c0d255e03fb511073434521bfb8de408453a09 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Thu, 6 Aug 2020 17:53:29 -0400 Subject: [PATCH 23/57] adding includes --- stan/math/prim/fun/conv_gaus_line.hpp | 4 ++-- stan/math/prim/fun/gaus_interp.hpp | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 7f3d8bdb1ea..2d252a3e025 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -1,10 +1,10 @@ #ifndef STAN_MATH_PRIM_FUN_CONV_GAUS_LINE #define STAN_MATH_PRIM_FUN_CONV_GAUS_LINE -#include -#include #include #include +#include +#include namespace stan { namespace math { diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 72e8d88637c..d6d086cde99 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -1,11 +1,12 @@ #ifndef STAN_MATH_PRIM_FUN_GAUS_INTERP #define STAN_MATH_PRIM_FUN_GAUS_INTERP +#include +#include +#include #include #include #include -#include -#include using namespace std; From f4dde7b37d58343b404ffddf9388cc407d9842f7 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Thu, 6 Aug 2020 18:25:14 -0400 Subject: [PATCH 24/57] adding includes --- stan/math/prim/fun/lin_interp.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/lin_interp.hpp index 22a6c4503a4..76c76d46b8a 100644 --- a/stan/math/prim/fun/lin_interp.hpp +++ b/stan/math/prim/fun/lin_interp.hpp @@ -1,6 +1,7 @@ #ifndef STAN_MATH_PRIM_FUN_LIN_INTERP #define STAN_MATH_PRIM_FUN_LIN_INTERP +#include #include #include #include From d37c84415f0299a4cb991d80aea48492bf65c4a7 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Fri, 7 Aug 2020 14:11:54 -0400 Subject: [PATCH 25/57] adding inline statements --- stan/math/prim/fun/conv_gaus_line.hpp | 4 ++-- stan/math/prim/fun/lin_interp.hpp | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 2d252a3e025..62e860249e7 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -23,8 +23,8 @@ namespace math { * @return The value of the derivative */ template -double conv_gaus_line(double t0, double t1, double a, double b, const Tx& x, - double sig2) { +inline double conv_gaus_line(double t0, double t1, double a, double b, + const Tx& x, double sig2) { using stan::math::normal_cdf; using std::exp; using std::pow; diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/lin_interp.hpp index 76c76d46b8a..e4409bb8b5e 100644 --- a/stan/math/prim/fun/lin_interp.hpp +++ b/stan/math/prim/fun/lin_interp.hpp @@ -23,8 +23,9 @@ namespace math { * @param x the point at which to evaluate the interpolation * @return value of linear interpolation at x */ -double lin_interp(const std::vector& xs, const std::vector& ys, - double x) { +inline double lin_interp(const std::vector& xs, + const std::vector& ys, + double x) { static char const* function = "lin_interp"; check_ordered(function, "xs", xs); check_not_nan(function, "xs", xs); From dea43062d007ecf31aaece7e4fffad5716597272 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 8 Aug 2020 14:02:03 +0000 Subject: [PATCH 26/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/conv_gaus_line.hpp | 4 ++-- stan/math/prim/fun/lin_interp.hpp | 5 ++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 62e860249e7..2988945b467 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -23,8 +23,8 @@ namespace math { * @return The value of the derivative */ template -inline double conv_gaus_line(double t0, double t1, double a, double b, - const Tx& x, double sig2) { +inline double conv_gaus_line(double t0, double t1, double a, double b, + const Tx& x, double sig2) { using stan::math::normal_cdf; using std::exp; using std::pow; diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/lin_interp.hpp index e4409bb8b5e..f992d95db71 100644 --- a/stan/math/prim/fun/lin_interp.hpp +++ b/stan/math/prim/fun/lin_interp.hpp @@ -23,9 +23,8 @@ namespace math { * @param x the point at which to evaluate the interpolation * @return value of linear interpolation at x */ -inline double lin_interp(const std::vector& xs, - const std::vector& ys, - double x) { +inline double lin_interp(const std::vector& xs, + const std::vector& ys, double x) { static char const* function = "lin_interp"; check_ordered(function, "xs", xs); check_not_nan(function, "xs", xs); From 2a3020bb0a72682238b191bb53634898142a5a68 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Sat, 8 Aug 2020 11:52:11 -0400 Subject: [PATCH 27/57] adding inline statements --- stan/math/prim/fun/gaus_interp.hpp | 6 +++--- stan/math/rev/fun/conv_gaus_line.hpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index d6d086cde99..4bdef248371 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -26,7 +26,7 @@ namespace internal { find the smallest difference between successive elements in a sorted vector */ template -double min_diff(int n, const std::vector& xs) { +inline double min_diff(int n, const std::vector& xs) { double dmin = value_of(xs[1]) - value_of(xs[0]); for (int i = 1; i < n - 1; i++) { if (value_of(xs[i + 1]) - value_of(xs[i]) < dmin) { @@ -123,8 +123,8 @@ inline return_type_t gaus_interp(const std::vector& xs, * @param ys vector of dependent variable of reference points * @return struct containing slopes, intercepts, and width of kernel */ -gaus_interp_params gaus_interp_precomp(const std::vector& xs, - const std::vector& ys) { +inline gaus_interp_params gaus_interp_precomp(const std::vector& xs, + const std::vector& ys) { static char const* function = "gaus_interp_precomp"; check_not_nan(function, "xs", xs); check_not_nan(function, "ys", ys); diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp index bfb88be7b4b..857ac45a8b7 100644 --- a/stan/math/rev/fun/conv_gaus_line.hpp +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -25,8 +25,8 @@ namespace internal { * @param sig2 variance of the Gaussian kernel * @return The value of the derivative */ -double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, - double sig2) { +inline double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, + double sig2) { using stan::math::normal_cdf; using std::exp; using std::pow; From cce952ba5290dd61b4f3988dbd441d9ac0e6e060 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 8 Aug 2020 15:55:24 +0000 Subject: [PATCH 28/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/gaus_interp.hpp | 2 +- stan/math/rev/fun/conv_gaus_line.hpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 4bdef248371..3a91814e20f 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -124,7 +124,7 @@ inline return_type_t gaus_interp(const std::vector& xs, * @return struct containing slopes, intercepts, and width of kernel */ inline gaus_interp_params gaus_interp_precomp(const std::vector& xs, - const std::vector& ys) { + const std::vector& ys) { static char const* function = "gaus_interp_precomp"; check_not_nan(function, "xs", xs); check_not_nan(function, "ys", ys); diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp index 857ac45a8b7..b13b130d273 100644 --- a/stan/math/rev/fun/conv_gaus_line.hpp +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -25,8 +25,8 @@ namespace internal { * @param sig2 variance of the Gaussian kernel * @return The value of the derivative */ -inline double der_conv_gaus_line(double t0, double t1, double a, double b, double x0, - double sig2) { +inline double der_conv_gaus_line(double t0, double t1, double a, double b, + double x0, double sig2) { using stan::math::normal_cdf; using std::exp; using std::pow; From e9e9786daecb876b2780dfb0286b87a130ad90d2 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 12 Aug 2020 10:48:55 -0400 Subject: [PATCH 29/57] redo order of initializations --- stan/math/rev/fun/conv_gaus_line.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp index 857ac45a8b7..0e2a5277f70 100644 --- a/stan/math/rev/fun/conv_gaus_line.hpp +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -55,12 +55,13 @@ class conv_gaus_line_vari : public op_v_vari { public: explicit conv_gaus_line_vari(double t0, double t1, double a, double b, vari* avi, double sig2) - : t0_(t0), + : op_v_vari(conv_gaus_line(t0, t1, a, b, avi->val_, sig2), avi), + t0_(t0), t1_(t1), a_(a), b_(b), - sig2_(sig2), - op_v_vari(conv_gaus_line(t0, t1, a, b, avi->val_, sig2), avi) {} + sig2_(sig2) {} + void chain() { avi_->adj_ From 4544f5e61d2c30e6f2745503c7c7240e5cabf49a Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 12 Aug 2020 15:06:27 +0000 Subject: [PATCH 30/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/conv_gaus_line.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp index 9b44c8a2f05..4fd1fb151b7 100644 --- a/stan/math/rev/fun/conv_gaus_line.hpp +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -55,14 +55,13 @@ class conv_gaus_line_vari : public op_v_vari { public: explicit conv_gaus_line_vari(double t0, double t1, double a, double b, vari* avi, double sig2) - : op_v_vari(conv_gaus_line(t0, t1, a, b, avi->val_, sig2), avi), - t0_(t0), + : op_v_vari(conv_gaus_line(t0, t1, a, b, avi->val_, sig2), avi), + t0_(t0), t1_(t1), a_(a), b_(b), sig2_(sig2) {} - void chain() { avi_->adj_ += adj_ * der_conv_gaus_line(t0_, t1_, a_, b_, avi_->val_, sig2_); From 5b9d65e74b1184aa5744fd61694712e54bf8d4a9 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 12 Aug 2020 11:34:19 -0400 Subject: [PATCH 31/57] update docs --- stan/math/prim/fun/gaus_interp.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 3a91814e20f..769e2ffb34d 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -52,6 +52,9 @@ inline double min_diff(int n, const std::vector& xs) { * are shifted and the process repeats itself until the interpolation * intersects all reference points. * + * Note: This interpolation scheme should be used when the function + * to be interpolated is well-resolved by the reference points. + * * @param xs vector of independent variable of reference points * @param ys vector of dependent variable of reference points * @param params a struct created by gaus_interp_precomp that From a5d00b6b2d839994d3c05532b27999d24aa91e1f Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 12 Aug 2020 15:35:24 +0000 Subject: [PATCH 32/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/gaus_interp.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 769e2ffb34d..c6193328990 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -52,8 +52,8 @@ inline double min_diff(int n, const std::vector& xs) { * are shifted and the process repeats itself until the interpolation * intersects all reference points. * - * Note: This interpolation scheme should be used when the function - * to be interpolated is well-resolved by the reference points. + * Note: This interpolation scheme should be used when the function + * to be interpolated is well-resolved by the reference points. * * @param xs vector of independent variable of reference points * @param ys vector of dependent variable of reference points From b401fafa01e5656cc3802cd702e7e67ac6cc6c26 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 12 Aug 2020 14:25:41 -0400 Subject: [PATCH 33/57] adding include --- test/unit/math/mix/fun/gaus_interp_test.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index da2ccf15029..ef651150d6c 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -1,6 +1,7 @@ #include #include #include +#include TEST(mathMixGausInterp, derivs) { using stan::math::gaus_interp_vect; From 2401a630aaea9649310e252cacdcb0aefd25a0f4 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 12 Aug 2020 16:50:14 -0400 Subject: [PATCH 34/57] rand_r fix --- test/unit/math/mix/fun/gaus_interp_test.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index ef651150d6c..13bc3c8a795 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -1,7 +1,6 @@ #include #include #include -#include TEST(mathMixGausInterp, derivs) { using stan::math::gaus_interp_vect; @@ -19,7 +18,7 @@ TEST(mathMixGausInterp, derivs) { for (int i = 0; i < n; i++) { x = xmin + i * (xmax - xmin) / (n - 1); xs.push_back(x); - y = rand_r(&seed) % 100; + y = rand(&seed) % 100; ys.push_back(y); } From 26b452923f1a918ee721543a3a44893325d8b1fc Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 12 Aug 2020 18:05:34 -0400 Subject: [PATCH 35/57] rand fix --- test/unit/math/mix/fun/gaus_interp_test.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index 13bc3c8a795..3b4b3595504 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -1,4 +1,5 @@ #include +#include #include #include From 561c909b7b67d6bc57bdf43c109b475ddf500671 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Thu, 13 Aug 2020 12:58:36 -0400 Subject: [PATCH 36/57] rand fix --- test/unit/math/mix/fun/gaus_interp_test.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index 3b4b3595504..fa920784b63 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -1,5 +1,4 @@ #include -#include #include #include @@ -19,7 +18,7 @@ TEST(mathMixGausInterp, derivs) { for (int i = 0; i < n; i++) { x = xmin + i * (xmax - xmin) / (n - 1); xs.push_back(x); - y = rand(&seed) % 100; + y = rand() % 100; ys.push_back(y); } From 2b98654888eec0cae1bcf56fe80735c00ed66990 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Fri, 14 Aug 2020 09:29:29 -0400 Subject: [PATCH 37/57] change rand --- test/unit/math/mix/fun/gaus_interp_test.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index fa920784b63..a4193548273 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -1,10 +1,13 @@ #include #include #include +#include TEST(mathMixGausInterp, derivs) { using stan::math::gaus_interp_vect; using stan::math::var; + std::default_random_engine generator; + std::uniform_real_distribution unif(0.0, 100.0); std::vector xs, ys, x2s, y2s, ts, as, bs, y3s, t3s; double xmin, xmax, x, y, x2, y2, t, t0, t1, dd, dder, dder2; double x0, x1, y0, y1, tmp; @@ -18,7 +21,7 @@ TEST(mathMixGausInterp, derivs) { for (int i = 0; i < n; i++) { x = xmin + i * (xmax - xmin) / (n - 1); xs.push_back(x); - y = rand() % 100; + y = unif(generator); ys.push_back(y); } From 0cc8ae90473cbeec1682bc3512c34c27f28af8f5 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Tue, 18 Aug 2020 18:06:22 -0400 Subject: [PATCH 38/57] modified mix test --- test/unit/math/mix/fun/gaus_interp_test.cpp | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index a4193548273..ac88f16206a 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -1,13 +1,10 @@ #include #include #include -#include TEST(mathMixGausInterp, derivs) { using stan::math::gaus_interp_vect; using stan::math::var; - std::default_random_engine generator; - std::uniform_real_distribution unif(0.0, 100.0); std::vector xs, ys, x2s, y2s, ts, as, bs, y3s, t3s; double xmin, xmax, x, y, x2, y2, t, t0, t1, dd, dder, dder2; double x0, x1, y0, y1, tmp; @@ -15,15 +12,20 @@ TEST(mathMixGausInterp, derivs) { unsigned int seed = 1; // generate function tabulation - n = 10; + n = 5; xmin = 0; xmax = 1; + xs.resize(n); + ys.resize(n); for (int i = 0; i < n; i++) { x = xmin + i * (xmax - xmin) / (n - 1); - xs.push_back(x); - y = unif(generator); - ys.push_back(y); + xs[i] = x; } + ys[0] = 0; + ys[1] = 0; + ys[2] = 1; + ys[3] = 0; + ys[4] = 0.01; // create vector of interpolation pts std::vector xs_new_v; From a3f988ba418f733cd174e041bbb6a3d2720abd21 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Thu, 20 Aug 2020 17:02:27 -0400 Subject: [PATCH 39/57] remove extra initializations --- test/unit/math/mix/fun/gaus_interp_test.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index ac88f16206a..ec8a7f9cc60 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -5,11 +5,9 @@ TEST(mathMixGausInterp, derivs) { using stan::math::gaus_interp_vect; using stan::math::var; - std::vector xs, ys, x2s, y2s, ts, as, bs, y3s, t3s; - double xmin, xmax, x, y, x2, y2, t, t0, t1, dd, dder, dder2; - double x0, x1, y0, y1, tmp; + std::vector xs, ys, ts; + double xmin, xmax, x, y, t, t0, t1, dder, dder2, x0, x1; int n; - unsigned int seed = 1; // generate function tabulation n = 5; From 85da4f662843b4ee1a88121a3d00341a33b78bdb Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Fri, 21 Aug 2020 08:37:34 -0400 Subject: [PATCH 40/57] remove using std --- stan/math/prim/fun/gaus_interp.hpp | 8 +++----- test/unit/math/mix/fun/gaus_interp_test.cpp | 3 --- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index c6193328990..216f160e03b 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -8,8 +8,6 @@ #include #include -using namespace std; - namespace stan { namespace math { @@ -184,9 +182,9 @@ inline gaus_interp_params gaus_interp_precomp(const std::vector& xs, */ template -inline vector gaus_interp_vect(const std::vector& xs, - const std::vector& ys, - const std::vector& xs_new) { +inline std::vector gaus_interp_vect(const std::vector& xs, + const std::vector& ys, + const std::vector& xs_new) { int n_interp = xs_new.size(); std::vector ys_new(n_interp); diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index ec8a7f9cc60..c678fb80764 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -1,6 +1,4 @@ #include -#include -#include TEST(mathMixGausInterp, derivs) { using stan::math::gaus_interp_vect; @@ -48,7 +46,6 @@ TEST(mathMixGausInterp, derivs) { ys_new_v = gaus_interp_vect(xs, ys, xs_new_v); std::vector ys_new_dder; - std::vector ys_new_v2; for (int i = 0; i < n_interp; i++) { ys_new_v[i].grad(); dder = xs_new_v[i].adj(); From 490cc03049eb2c67cfddfb642f4c8ff44eb99db4 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 21 Aug 2020 12:38:33 +0000 Subject: [PATCH 41/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/gaus_interp.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 216f160e03b..fbb34829d23 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -183,8 +183,8 @@ inline gaus_interp_params gaus_interp_precomp(const std::vector& xs, template inline std::vector gaus_interp_vect(const std::vector& xs, - const std::vector& ys, - const std::vector& xs_new) { + const std::vector& ys, + const std::vector& xs_new) { int n_interp = xs_new.size(); std::vector ys_new(n_interp); From 22280c91fbfde43db6e12e50ee7abfdcfe206173 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Fri, 21 Aug 2020 10:24:53 -0400 Subject: [PATCH 42/57] namespace issues --- .../math/prim/fun/conv_gaus_line_test.cpp | 2 +- test/unit/math/prim/fun/gaus_interp_test.cpp | 20 +++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/test/unit/math/prim/fun/conv_gaus_line_test.cpp b/test/unit/math/prim/fun/conv_gaus_line_test.cpp index 2847308cee8..67a80cfcc3d 100644 --- a/test/unit/math/prim/fun/conv_gaus_line_test.cpp +++ b/test/unit/math/prim/fun/conv_gaus_line_test.cpp @@ -8,7 +8,7 @@ evaluate the convolution using trapezoid rule double check_int_dumb(double t0, double t1, double a, double b, double x0, double sig2) { double x, y, t, tot; - vector ys, ts; + std::vector ys, ts; double pi = stan::math::pi(); int n; diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index 3f7b1427e7e..bab21ba44c4 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -10,7 +10,7 @@ TEST(mathPrimGausInterp, throwing) { using stan::math::gaus_interp_precomp; double nan = std::numeric_limits::quiet_NaN(); double x; - vector xs, ys; + std::vector xs, ys; // check that when xs are not increasing, an error throws int n = 2; @@ -73,11 +73,11 @@ TEST(mathPrimGausInterp, interp_line) { int n = 2; double xmin = 0; double xmax = 1; - vector xs = {0, 1}; - vector ys = {0, 2}; + std::vector xs = {0, 1}; + std::vector ys = {0, 2}; // vector of pts at which to compute interpolation - vector xs_new; + std::vector xs_new; int n_interp = 100; double t0 = xmin; double t1 = xmax; @@ -88,14 +88,14 @@ TEST(mathPrimGausInterp, interp_line) { } // create interpolation using precomp - vector ys_new_gaus(n_interp); + std::vector ys_new_gaus(n_interp); gaus_interp_params params = gaus_interp_precomp(xs, ys); for (int i = 0; i < n_interp; i++) { ys_new_gaus[i] = gaus_interp(xs, ys, params, xs_new[i]); } // create interpolation without precomp - vector ys_new_gaus2 = gaus_interp_vect(xs, ys, xs_new); + std::vector ys_new_gaus2 = gaus_interp_vect(xs, ys, xs_new); // test points double tmp, y; @@ -119,7 +119,7 @@ TEST(mathPrimGausInterp, gaus_and_lin_interp) { double xmin = 0; double xmax = 1; double x; - vector xs, ys; + std::vector xs, ys; for (int i = 0; i < n; i++) { x = xmin + i * (xmax - xmin) / (n - 1); xs.push_back(x); @@ -127,7 +127,7 @@ TEST(mathPrimGausInterp, gaus_and_lin_interp) { } // vector of pts at which to compute interpolation - vector xs_new; + std::vector xs_new; int n_interp = 100; double t0 = xmin; double t1 = xmax; @@ -138,8 +138,8 @@ TEST(mathPrimGausInterp, gaus_and_lin_interp) { } // create interpolation - vector ys_new_gaus(n_interp); - vector ys_new_lin(n_interp); + std::vector ys_new_gaus(n_interp); + std::vector ys_new_lin(n_interp); // linear interpolation ys_new_lin.resize(n_interp); From f5bbc96da71e3b63022b2ed5f9459a9e6da4c7c1 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Sat, 22 Aug 2020 07:40:55 -0400 Subject: [PATCH 43/57] namespace issues --- test/unit/math/prim/fun/lin_interp_test.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/unit/math/prim/fun/lin_interp_test.cpp b/test/unit/math/prim/fun/lin_interp_test.cpp index d4576deabd8..1342be5aa90 100644 --- a/test/unit/math/prim/fun/lin_interp_test.cpp +++ b/test/unit/math/prim/fun/lin_interp_test.cpp @@ -8,7 +8,7 @@ TEST(mathPrimLinInterp, throwing) { using stan::math::lin_interp; double nan = std::numeric_limits::quiet_NaN(); double x = 0.5; - vector xs, ys; + std::vector xs, ys; // check that when xs are not increasing, an error throws int n = 2; @@ -47,11 +47,11 @@ TEST(mathPrimLinInterp, interp_line) { int n = 2; double xmin = 0; double xmax = 1; - vector xs = {0, 1}; - vector ys = {0, 2}; + std::vector xs = {0, 1}; + std::vector ys = {0, 2}; // vector of pts at which to compute interpolation - vector xs_new; + std::vector xs_new; int n_interp = 100; double t0 = xmin; double t1 = xmax; @@ -62,7 +62,7 @@ TEST(mathPrimLinInterp, interp_line) { } // interpolate - vector ys_new(n_interp); + std::vector ys_new(n_interp); for (int i = 0; i < n_interp; i++) { ys_new[i] = lin_interp(xs, ys, xs_new[i]); } @@ -80,8 +80,8 @@ TEST(mathPrimLinInterp, interp_line) { ASSERT_NEAR(lin_interp(xs, ys, 100), ys[1], ABS_TOL); // xs with more than 2 elements - vector xs2 = {0, 1, 2, 3, 4, 5, 6}; - vector ys2 = {0, 2, 5, 2, 3, 2, 2}; + std::vector xs2 = {0, 1, 2, 3, 4, 5, 6}; + std::vector ys2 = {0, 2, 5, 2, 3, 2, 2}; double x = 0.5; ASSERT_NEAR(lin_interp(xs, ys, x), lin_interp(xs2, ys2, x), ABS_TOL); } From 0e32b0ce87c002e102bc6e07fe3ac3e962838064 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Mon, 31 Aug 2020 16:52:47 -0400 Subject: [PATCH 44/57] params changed to vector --- stan/math/prim/fun/gaus_interp.hpp | 61 ++++++++++++-------- test/unit/math/prim/fun/gaus_interp_test.cpp | 9 ++- 2 files changed, 42 insertions(+), 28 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index fbb34829d23..46d596be735 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -11,18 +11,23 @@ namespace stan { namespace math { + // struct for passing precomputation data struct gaus_interp_params { - std::vector as; - std::vector bs; - double sig2; + std::vector as_; + std::vector bs_; + double sig2_; + gaus_interp_params(size_t as_size, size_t bs_size, double sig2) : + as_(as_size), bs_(bs_size), sig2_(sig2) {} +// std::vector as; +// std::vector bs; +// double sig2; }; namespace internal { - /* -find the smallest difference between successive elements in a sorted vector -*/ + * find the smallest difference between successive elements in a sorted vector + */ template inline double min_diff(int n, const std::vector& xs) { double dmin = value_of(xs[1]) - value_of(xs[0]); @@ -62,7 +67,7 @@ inline double min_diff(int n, const std::vector& xs) { template inline return_type_t gaus_interp(const std::vector& xs, const std::vector& ys, - const gaus_interp_params& params, + const std::vector& params2, const Tx& x) { // enforce that interpolation point is between smallest and largest // reference point @@ -78,11 +83,14 @@ inline return_type_t gaus_interp(const std::vector& xs, const double NSTDS = 10; int n = xs.size(); + // params2 is of the form (as, bs, sig2) + // create copy of xs so that endpoints can be extended std::vector xs2 = xs; // extend out first and last lines for convolution - double sig = std::sqrt(params.sig2); + double sig = std::sqrt(params2.back()); + std::cout << "sig: " << sig << std::endl; xs2[0] = xs[0] - NSTDS * sig; xs2[n - 1] = xs[n - 1] + NSTDS * sig; @@ -100,8 +108,8 @@ inline return_type_t gaus_interp(const std::vector& xs, using return_t = return_type_t; return_t y = 0; for (int i = ind_start; i < ind_end; i++) { - y += conv_gaus_line(xs2[i], xs2[i + 1], params.as[i], params.bs[i], x, - params.sig2); + y += conv_gaus_line(xs2[i], xs2[i + 1], params2[i], params2[(n-1) + i], x, + params2.back()); } return y; @@ -124,8 +132,8 @@ inline return_type_t gaus_interp(const std::vector& xs, * @param ys vector of dependent variable of reference points * @return struct containing slopes, intercepts, and width of kernel */ -inline gaus_interp_params gaus_interp_precomp(const std::vector& xs, - const std::vector& ys) { +inline std::vector gaus_interp_precomp(const std::vector& xs, + const std::vector& ys) { static char const* function = "gaus_interp_precomp"; check_not_nan(function, "xs", xs); check_not_nan(function, "ys", ys); @@ -133,16 +141,19 @@ inline gaus_interp_params gaus_interp_precomp(const std::vector& xs, check_greater(function, "xs", xs.size(), 1); using internal::min_diff; - gaus_interp_params params; const double INTERP_TOL = 1e-8; const double SIG2_SCALE = 0.1; int n = xs.size(); - - // find minimum distance between points for std of gaussian kernel - params.sig2 = square(min_diff(n, xs) * SIG2_SCALE); - params.as.resize(n - 1); - params.bs.resize(n - 1); - + gaus_interp_params params(n-1, n-1, square(min_diff(n, xs) * SIG2_SCALE)); + std::vector params2; + params2.resize(2*n-1); + params2[2*n-2] = square(min_diff(n, xs) * SIG2_SCALE); + std::cout << "sig2: " << square(min_diff(n, xs) * SIG2_SCALE) << std::endl; + std::cout << "n: " << n << std::endl; + std::cout << "params2.size(): " << params2.size() << std::endl; + std::cout << "params2.back(): " << params2.back() << std::endl; + std::cout << "params2[2*n-2]: " << params2[2*n-2] << std::endl; + std::cout << "min_diff(n, xs): " << min_diff(n, xs) << std::endl; // copy ys into a new vector that will be changed std::vector y2s = ys; @@ -152,20 +163,20 @@ inline gaus_interp_params gaus_interp_precomp(const std::vector& xs, for (int j = 0; j < max_iters; j++) { // find slope and intercept of line between each point for (int i = 0; i < n - 1; i++) { - params.as[i] = (y2s[i + 1] - y2s[i]) / (xs[i + 1] - xs[i]); - params.bs[i] = -xs[i] * params.as[i] + y2s[i]; + params2[i] = (y2s[i + 1] - y2s[i]) / (xs[i + 1] - xs[i]); + params2[(n - 1) + i] = -xs[i] * params2[i] + y2s[i]; } dmax = 0; for (int i = 0; i < n; i++) { - dd = ys[i] - gaus_interp(xs, y2s, params, xs[i]); + dd = ys[i] - gaus_interp(xs, y2s, params2, xs[i]); y2s[i] += dd; dmax = std::max(std::abs(dd), dmax); } if (dmax < INTERP_TOL) break; } - return params; + return params2; } /** @@ -189,9 +200,9 @@ inline std::vector gaus_interp_vect(const std::vector& xs, std::vector ys_new(n_interp); // create interpolation - gaus_interp_params params = gaus_interp_precomp(xs, ys); + std::vector params2 = gaus_interp_precomp(xs, ys); for (int i = 0; i < n_interp; i++) { - ys_new[i] = gaus_interp(xs, ys, params, xs_new[i]); + ys_new[i] = gaus_interp(xs, ys, params2, xs_new[i]); } return ys_new; } diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index bab21ba44c4..5566e382a28 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -31,7 +31,8 @@ TEST(mathPrimGausInterp, throwing) { // check that error throws when trying to interpolate out of range or nan xs = {0, 1}; ys = {0, 2}; - gaus_interp_params params = gaus_interp_precomp(xs, ys); + //gaus_interp_params params = gaus_interp_precomp(xs, ys); + std::vector params = gaus_interp_precomp(xs, ys); x = 1.1; EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); x = -0.1; @@ -89,7 +90,8 @@ TEST(mathPrimGausInterp, interp_line) { // create interpolation using precomp std::vector ys_new_gaus(n_interp); - gaus_interp_params params = gaus_interp_precomp(xs, ys); + //gaus_interp_params params = gaus_interp_precomp(xs, ys); + std::vector params = gaus_interp_precomp(xs, ys); for (int i = 0; i < n_interp; i++) { ys_new_gaus[i] = gaus_interp(xs, ys, params, xs_new[i]); } @@ -148,7 +150,8 @@ TEST(mathPrimGausInterp, gaus_and_lin_interp) { } // gaus interpolation - gaus_interp_params params = gaus_interp_precomp(xs, ys); + //gaus_interp_params params = gaus_interp_precomp(xs, ys); + std::vector params = gaus_interp_precomp(xs, ys); for (int i = 0; i < n_interp; i++) { ys_new_gaus[i] = gaus_interp(xs, ys, params, xs_new[i]); } From dc49bc1b52ce35738d125022faa44de080a0b21a Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Tue, 1 Sep 2020 11:28:00 -0400 Subject: [PATCH 45/57] added new conv_gaus_line for sum of convolutions --- stan/math/prim/fun/conv_gaus_line.hpp | 41 +++++++++ stan/math/prim/fun/gaus_interp.hpp | 89 ++++++++------------ test/unit/math/prim/fun/gaus_interp_test.cpp | 32 +++++-- 3 files changed, 102 insertions(+), 60 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 2988945b467..5fac34e016b 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -41,6 +41,47 @@ inline double conv_gaus_line(double t0, double t1, double a, double b, return y; } + +/** + * Evaluate the convolution of a piecewise linear function with a + * Gaussian kernel. + * + * \f$\int_{x_0}^{x_1} (a_1t + b_1) e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ + * \f$ + ... + $\int_{x_{n-1}}^{x_{n}} (a_{n-1}t + b_{n-1}) \f$ + * \f$ e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ + * + * @param x point at which convolution is evaluated + * @param xs increasing vector with endpoints of each piecewise linear function + * @param params vector containing slopes and intercepts of peicewise linear function and width of Gaussian kernel + * @param ind_start + * @param ind_end + * @return The value of the derivative + */ +template +inline double conv_gaus_line_sum(const Tx& x, std::vector xs, + std::vector params, + int ind_start, + int ind_end) { + using stan::math::normal_cdf; + using std::exp; + using std::pow; + using std::sqrt; + const double pi = stan::math::pi(); + int n = xs.size(); + double sig2 = params[2*n-2]; + double sig = std::sqrt(sig2); + + double y=0; + for (int i = ind_start; i < ind_end; i++) { + y += (params[i] * value_of(x) + params[n-1 + i]) + * (normal_cdf(xs[i+1], value_of(x), sig) - normal_cdf(xs[i], value_of(x), sig)); + y += -params[i] * sig2 / sqrt(2 * pi * sig2) + * (exp(-pow(xs[i+1] - value_of(x), 2) / (2 * sig2)) + - exp(-pow(xs[i] - value_of(x), 2) / (2 * sig2))); + } + return y; +} + } // namespace math } // namespace stan diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 46d596be735..6ed05b2ead3 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -11,19 +11,6 @@ namespace stan { namespace math { - -// struct for passing precomputation data -struct gaus_interp_params { - std::vector as_; - std::vector bs_; - double sig2_; - gaus_interp_params(size_t as_size, size_t bs_size, double sig2) : - as_(as_size), bs_(bs_size), sig2_(sig2) {} -// std::vector as; -// std::vector bs; -// double sig2; -}; - namespace internal { /* * find the smallest difference between successive elements in a sorted vector @@ -60,14 +47,14 @@ inline double min_diff(int n, const std::vector& xs) { * * @param xs vector of independent variable of reference points * @param ys vector of dependent variable of reference points - * @param params a struct created by gaus_interp_precomp that + * @param params a vector created by gaus_interp_precomp * @param x the point at which to evaluate the interpolation * @return value of the interpolation at x */ template inline return_type_t gaus_interp(const std::vector& xs, const std::vector& ys, - const std::vector& params2, + const std::vector& params, const Tx& x) { // enforce that interpolation point is between smallest and largest // reference point @@ -80,37 +67,39 @@ inline return_type_t gaus_interp(const std::vector& xs, check_not_nan(function, "x", x); check_greater(function, "xs", xs.size(), 1); + // number of standard deviations to extend endpoints for convolution const double NSTDS = 10; int n = xs.size(); - // params2 is of the form (as, bs, sig2) + // params is vector of the form (as, bs, sig2) // create copy of xs so that endpoints can be extended std::vector xs2 = xs; // extend out first and last lines for convolution - double sig = std::sqrt(params2.back()); - std::cout << "sig: " << sig << std::endl; + double sig = std::sqrt(params[2*n-2]); xs2[0] = xs[0] - NSTDS * sig; xs2[n - 1] = xs[n - 1] + NSTDS * sig; // no need to convolve far from center of gaussian, so // get lower and upper indexes for integration bounds - auto lb = upper_bound(xs.begin(), xs.end(), x - NSTDS * sig); - int ind_start = distance(xs.begin(), lb) - 1; + auto lb = std::lower_bound(xs.begin(), xs.end(), x - NSTDS * sig); + int ind_start = &(*lb) - &xs[0] - 1; ind_start = std::max(0, ind_start); - auto ub = upper_bound(xs.begin(), xs.end(), x + NSTDS * sig); - int ind_end = distance(xs.begin(), ub); + auto ub = std::upper_bound(xs.begin(), xs.end(), x + NSTDS * sig); + int ind_end = &(*ub) - &xs[0]; ind_end = std::min(n - 1, ind_end); // sum convolutions over intervals - using return_t = return_type_t; - return_t y = 0; + /* + return_type_t y = 0; for (int i = ind_start; i < ind_end; i++) { - y += conv_gaus_line(xs2[i], xs2[i + 1], params2[i], params2[(n-1) + i], x, - params2.back()); + y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(n-1) + i], x, + params[2*n-2]); } + */ + return_type_t y = conv_gaus_line_sum(x, xs2, params, ind_start, ind_end); return y; } @@ -141,42 +130,38 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, check_greater(function, "xs", xs.size(), 1); using internal::min_diff; - const double INTERP_TOL = 1e-8; - const double SIG2_SCALE = 0.1; + static const double max_diff = 1e-8; + static const double sig2_scale = 0.1; int n = xs.size(); - gaus_interp_params params(n-1, n-1, square(min_diff(n, xs) * SIG2_SCALE)); - std::vector params2; - params2.resize(2*n-1); - params2[2*n-2] = square(min_diff(n, xs) * SIG2_SCALE); - std::cout << "sig2: " << square(min_diff(n, xs) * SIG2_SCALE) << std::endl; - std::cout << "n: " << n << std::endl; - std::cout << "params2.size(): " << params2.size() << std::endl; - std::cout << "params2.back(): " << params2.back() << std::endl; - std::cout << "params2[2*n-2]: " << params2[2*n-2] << std::endl; - std::cout << "min_diff(n, xs): " << min_diff(n, xs) << std::endl; + + // create the vector to be returned that consists of as, bs, sig2 + std::vector params; + params.resize(2*n-1); + params[2*n-2] = square(min_diff(n, xs) * sig2_scale); + // copy ys into a new vector that will be changed std::vector y2s = ys; // interatively find interpolation that coincides with ys at xs int max_iters = 50; - double dmax, dd; + double dd; for (int j = 0; j < max_iters; j++) { // find slope and intercept of line between each point - for (int i = 0; i < n - 1; i++) { - params2[i] = (y2s[i + 1] - y2s[i]) / (xs[i + 1] - xs[i]); - params2[(n - 1) + i] = -xs[i] * params2[i] + y2s[i]; + for (size_t i = 0; i < n - 1; i++) { + params[i] = (y2s[i + 1] - y2s[i]) / (xs[i + 1] - xs[i]); + params[(n - 1) + i] = -xs[i] * params[i] + y2s[i]; } - dmax = 0; - for (int i = 0; i < n; i++) { - dd = ys[i] - gaus_interp(xs, y2s, params2, xs[i]); + double dmax = 0; + for (size_t i = 0; i < n; i++) { + dd = ys[i] - gaus_interp(xs, y2s, params, xs[i]); y2s[i] += dd; dmax = std::max(std::abs(dd), dmax); } - if (dmax < INTERP_TOL) + if (dmax < max_diff) break; } - return params2; + return params; } /** @@ -193,16 +178,16 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, */ template -inline std::vector gaus_interp_vect(const std::vector& xs, - const std::vector& ys, - const std::vector& xs_new) { +inline std::vector gaus_interp(const std::vector& xs, + const std::vector& ys, + const std::vector& xs_new) { int n_interp = xs_new.size(); std::vector ys_new(n_interp); // create interpolation - std::vector params2 = gaus_interp_precomp(xs, ys); + std::vector params = gaus_interp_precomp(xs, ys); for (int i = 0; i < n_interp; i++) { - ys_new[i] = gaus_interp(xs, ys, params2, xs_new[i]); + ys_new[i] = gaus_interp(xs, ys, params, xs_new[i]); } return ys_new; } diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/gaus_interp_test.cpp index 5566e382a28..3571792b246 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/gaus_interp_test.cpp @@ -6,7 +6,6 @@ double ABS_TOL = 1e-12; TEST(mathPrimGausInterp, throwing) { using stan::math::gaus_interp; - using stan::math::gaus_interp_params; using stan::math::gaus_interp_precomp; double nan = std::numeric_limits::quiet_NaN(); double x; @@ -31,7 +30,6 @@ TEST(mathPrimGausInterp, throwing) { // check that error throws when trying to interpolate out of range or nan xs = {0, 1}; ys = {0, 2}; - //gaus_interp_params params = gaus_interp_precomp(xs, ys); std::vector params = gaus_interp_precomp(xs, ys); x = 1.1; EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); @@ -65,9 +63,7 @@ TEST(mathPrimGausInterp, throwing) { TEST(mathPrimGausInterp, interp_line) { using stan::math::gaus_interp; - using stan::math::gaus_interp_params; using stan::math::gaus_interp_precomp; - using stan::math::gaus_interp_vect; // check that interpolation of line returns the same function // generate function tabulation @@ -90,14 +86,13 @@ TEST(mathPrimGausInterp, interp_line) { // create interpolation using precomp std::vector ys_new_gaus(n_interp); - //gaus_interp_params params = gaus_interp_precomp(xs, ys); std::vector params = gaus_interp_precomp(xs, ys); for (int i = 0; i < n_interp; i++) { ys_new_gaus[i] = gaus_interp(xs, ys, params, xs_new[i]); } // create interpolation without precomp - std::vector ys_new_gaus2 = gaus_interp_vect(xs, ys, xs_new); + std::vector ys_new_gaus2 = gaus_interp(xs, ys, xs_new); // test points double tmp, y; @@ -109,9 +104,31 @@ TEST(mathPrimGausInterp, interp_line) { } } +TEST(mathPrimGausInterp, matching_reference_interp_pts) { + using stan::math::gaus_interp; + + // check that interpolation returns the same function + // when interpolation points are the same as reference points + + // generate function tabulation + int n = 3; + std::vector xs = {0, 1, 2}; + std::vector ys = {0, 2, 1}; + + // create interpolation points + std::vector xs_new = xs; + + // create interpolation + std::vector ys_new = gaus_interp(xs, ys, xs_new); + + // test points + for (int i = 0; i < n; i++) { + ASSERT_NEAR(ys_new[i], ys[i], 1e-8); + } +} + TEST(mathPrimGausInterp, gaus_and_lin_interp) { using stan::math::gaus_interp; - using stan::math::gaus_interp_params; using stan::math::gaus_interp_precomp; using stan::math::lin_interp; @@ -150,7 +167,6 @@ TEST(mathPrimGausInterp, gaus_and_lin_interp) { } // gaus interpolation - //gaus_interp_params params = gaus_interp_precomp(xs, ys); std::vector params = gaus_interp_precomp(xs, ys); for (int i = 0; i < n_interp; i++) { ys_new_gaus[i] = gaus_interp(xs, ys, params, xs_new[i]); From c1cc956dd2e15ea65fc230e417266791ec431694 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 2 Sep 2020 08:31:48 -0400 Subject: [PATCH 46/57] various fixes --- stan/math/prim/fun/gaus_interp.hpp | 2 +- stan/math/prim/fun/lin_interp.hpp | 37 ++++++++++++++++++--- stan/math/rev/fun/conv_gaus_line.hpp | 3 +- test/unit/math/mix/fun/gaus_interp_test.cpp | 1 + test/unit/math/prim/fun/lin_interp_test.cpp | 24 +++++++++++++ 5 files changed, 59 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 6ed05b2ead3..cc10df09d12 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -167,7 +167,7 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, /** * This function combines gaus_interp_precomp and gaus_interp. * It takes as input two vectors of reference points (xs and ys) - * in addition to a vector, xs_new, of points at which the + * in addition to a vector, xs_new, of points at which this * function will evaluate the interpolation through those reference * points. * diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/lin_interp.hpp index f992d95db71..bcc002077b7 100644 --- a/stan/math/prim/fun/lin_interp.hpp +++ b/stan/math/prim/fun/lin_interp.hpp @@ -42,22 +42,49 @@ inline double lin_interp(const std::vector& xs, } // find in between which points the input, x, lives - auto ub = upper_bound(xs.begin(), xs.end(), x); - int ind = distance(xs.begin(), ub); + auto ub = std::upper_bound(xs.begin(), xs.end(), x); + auto ind = std::distance(xs.begin(), ub); // check if the interpolation point falls on a reference point - if (x == *(ub - 1)) + if (x == xs[ind - 1]) { return ys[ind - 1]; + } // do linear interpolation - double x1 = *(ub - 1); - double x2 = *ub; + double x1 = xs[ind-1]; + double x2 = xs[ind]; double y1 = ys[ind - 1]; double y2 = ys[ind]; return y1 + (x - x1) * (y2 - y1) / (x2 - x1); } +/** + * This function takes as input two vectors of reference points (xs and ys) + * in addition to a vector, xs_new, of points at which this + * function will evaluate a linear interpolation through those reference + * points. + * + * @param xs vector of independent variable of reference points + * @param ys vector of dependent variable of reference points + * @param xs_new vector of point at which to evaluate interpolation + * @return vector of interpolation values + */ + +template +inline std::vector lin_interp(const std::vector& xs, + const std::vector& ys, + const std::vector& xs_new) { + int n_interp = xs_new.size(); + std::vector ys_new(n_interp); + + // create interpolation + for (int i = 0; i < n_interp; i++) { + ys_new[i] = lin_interp(xs, ys, xs_new[i]); + } + return ys_new; +} + } // namespace math } // namespace stan diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp index 4fd1fb151b7..0551e9e957e 100644 --- a/stan/math/rev/fun/conv_gaus_line.hpp +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -31,9 +31,8 @@ inline double der_conv_gaus_line(double t0, double t1, double a, double b, using std::exp; using std::pow; using std::sqrt; - const double pi = stan::math::pi(); const double sig = sqrt(sig2); - const double alpha = sqrt(2 * pi * sig2); + const double alpha = sqrt(2 * pi() * sig2); double y = (a * x0 + b) / alpha * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index c678fb80764..0c5dee69e70 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -29,6 +29,7 @@ TEST(mathMixGausInterp, derivs) { int n_interp = 100; // add a cushion to each endpoint so that we don't leave interval + // when taking finite differences t0 = xs[0] + 0.01; t1 = xmax - 0.01; for (int i = 0; i < n_interp; i++) { diff --git a/test/unit/math/prim/fun/lin_interp_test.cpp b/test/unit/math/prim/fun/lin_interp_test.cpp index 1342be5aa90..48c5464d474 100644 --- a/test/unit/math/prim/fun/lin_interp_test.cpp +++ b/test/unit/math/prim/fun/lin_interp_test.cpp @@ -85,3 +85,27 @@ TEST(mathPrimLinInterp, interp_line) { double x = 0.5; ASSERT_NEAR(lin_interp(xs, ys, x), lin_interp(xs2, ys2, x), ABS_TOL); } + +TEST(mathPrimLinInterp, matching_reference_interp_pts) { + using stan::math::lin_interp; + + // check that interpolation returns the same function + // when interpolation points are the same as reference points + + // generate function tabulation + int n = 3; + std::vector xs = {0, 1, 2}; + std::vector ys = {0, 2, 1}; + + // create interpolation points, same as reference points + std::vector xs_new = xs; + + // create interpolation + std::vector ys_new = lin_interp(xs, ys, xs_new); + + // test points + for (int i = 0; i < n; i++) { + ASSERT_NEAR(ys_new[i], ys[i], ABS_TOL); + } +} + From 09829bec5e7f14e27dc8ea224c7c71fe5190216d Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 2 Sep 2020 08:53:32 -0400 Subject: [PATCH 47/57] change distance to pointer arithmetic --- stan/math/prim/fun/lin_interp.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/lin_interp.hpp index bcc002077b7..7a612291095 100644 --- a/stan/math/prim/fun/lin_interp.hpp +++ b/stan/math/prim/fun/lin_interp.hpp @@ -43,7 +43,7 @@ inline double lin_interp(const std::vector& xs, // find in between which points the input, x, lives auto ub = std::upper_bound(xs.begin(), xs.end(), x); - auto ind = std::distance(xs.begin(), ub); + int ind = &(*ub) - &(xs[0]); // check if the interpolation point falls on a reference point if (x == xs[ind - 1]) { From 91199f830759edc22f872a46c053faf1195fbb26 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 2 Sep 2020 09:04:29 -0400 Subject: [PATCH 48/57] [Jenkins] auto-formatting by clang-format version 6.0.1-14 (tags/RELEASE_601/final) --- stan/math/prim/fun/conv_gaus_line.hpp | 24 ++++++++++----------- stan/math/prim/fun/gaus_interp.hpp | 12 +++++------ stan/math/prim/fun/lin_interp.hpp | 6 +++--- test/unit/math/prim/fun/lin_interp_test.cpp | 1 - 4 files changed, 21 insertions(+), 22 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 5fac34e016b..84bc020778a 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -41,9 +41,8 @@ inline double conv_gaus_line(double t0, double t1, double a, double b, return y; } - /** - * Evaluate the convolution of a piecewise linear function with a + * Evaluate the convolution of a piecewise linear function with a * Gaussian kernel. * * \f$\int_{x_0}^{x_1} (a_1t + b_1) e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ @@ -52,32 +51,33 @@ inline double conv_gaus_line(double t0, double t1, double a, double b, * * @param x point at which convolution is evaluated * @param xs increasing vector with endpoints of each piecewise linear function - * @param params vector containing slopes and intercepts of peicewise linear function and width of Gaussian kernel + * @param params vector containing slopes and intercepts of peicewise linear + * function and width of Gaussian kernel * @param ind_start * @param ind_end * @return The value of the derivative */ template inline double conv_gaus_line_sum(const Tx& x, std::vector xs, - std::vector params, - int ind_start, - int ind_end) { + std::vector params, int ind_start, + int ind_end) { using stan::math::normal_cdf; using std::exp; using std::pow; using std::sqrt; const double pi = stan::math::pi(); int n = xs.size(); - double sig2 = params[2*n-2]; + double sig2 = params[2 * n - 2]; double sig = std::sqrt(sig2); - double y=0; + double y = 0; for (int i = ind_start; i < ind_end; i++) { - y += (params[i] * value_of(x) + params[n-1 + i]) - * (normal_cdf(xs[i+1], value_of(x), sig) - normal_cdf(xs[i], value_of(x), sig)); + y += (params[i] * value_of(x) + params[n - 1 + i]) + * (normal_cdf(xs[i + 1], value_of(x), sig) + - normal_cdf(xs[i], value_of(x), sig)); y += -params[i] * sig2 / sqrt(2 * pi * sig2) - * (exp(-pow(xs[i+1] - value_of(x), 2) / (2 * sig2)) - - exp(-pow(xs[i] - value_of(x), 2) / (2 * sig2))); + * (exp(-pow(xs[i + 1] - value_of(x), 2) / (2 * sig2)) + - exp(-pow(xs[i] - value_of(x), 2) / (2 * sig2))); } return y; } diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index cc10df09d12..da9d6486fd7 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -77,7 +77,7 @@ inline return_type_t gaus_interp(const std::vector& xs, std::vector xs2 = xs; // extend out first and last lines for convolution - double sig = std::sqrt(params[2*n-2]); + double sig = std::sqrt(params[2 * n - 2]); xs2[0] = xs[0] - NSTDS * sig; xs2[n - 1] = xs[n - 1] + NSTDS * sig; @@ -122,7 +122,7 @@ inline return_type_t gaus_interp(const std::vector& xs, * @return struct containing slopes, intercepts, and width of kernel */ inline std::vector gaus_interp_precomp(const std::vector& xs, - const std::vector& ys) { + const std::vector& ys) { static char const* function = "gaus_interp_precomp"; check_not_nan(function, "xs", xs); check_not_nan(function, "ys", ys); @@ -136,8 +136,8 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, // create the vector to be returned that consists of as, bs, sig2 std::vector params; - params.resize(2*n-1); - params[2*n-2] = square(min_diff(n, xs) * sig2_scale); + params.resize(2 * n - 1); + params[2 * n - 2] = square(min_diff(n, xs) * sig2_scale); // copy ys into a new vector that will be changed std::vector y2s = ys; @@ -179,8 +179,8 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, template inline std::vector gaus_interp(const std::vector& xs, - const std::vector& ys, - const std::vector& xs_new) { + const std::vector& ys, + const std::vector& xs_new) { int n_interp = xs_new.size(); std::vector ys_new(n_interp); diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/lin_interp.hpp index 7a612291095..efce23cb855 100644 --- a/stan/math/prim/fun/lin_interp.hpp +++ b/stan/math/prim/fun/lin_interp.hpp @@ -51,7 +51,7 @@ inline double lin_interp(const std::vector& xs, } // do linear interpolation - double x1 = xs[ind-1]; + double x1 = xs[ind - 1]; double x2 = xs[ind]; double y1 = ys[ind - 1]; double y2 = ys[ind]; @@ -73,8 +73,8 @@ inline double lin_interp(const std::vector& xs, template inline std::vector lin_interp(const std::vector& xs, - const std::vector& ys, - const std::vector& xs_new) { + const std::vector& ys, + const std::vector& xs_new) { int n_interp = xs_new.size(); std::vector ys_new(n_interp); diff --git a/test/unit/math/prim/fun/lin_interp_test.cpp b/test/unit/math/prim/fun/lin_interp_test.cpp index 48c5464d474..a19590ec62f 100644 --- a/test/unit/math/prim/fun/lin_interp_test.cpp +++ b/test/unit/math/prim/fun/lin_interp_test.cpp @@ -108,4 +108,3 @@ TEST(mathPrimLinInterp, matching_reference_interp_pts) { ASSERT_NEAR(ys_new[i], ys[i], ABS_TOL); } } - From 07060e70bf4c5e99b3859dc57816daf3ccfeb9ce Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Thu, 3 Sep 2020 08:23:33 -0400 Subject: [PATCH 49/57] revert conv_gaus_line --- stan/math/prim/fun/conv_gaus_line.hpp | 57 +++------------------ stan/math/prim/fun/gaus_interp.hpp | 6 +-- test/unit/math/mix/fun/gaus_interp_test.cpp | 12 ++--- 3 files changed, 16 insertions(+), 59 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 5fac34e016b..aaea6a39c88 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -23,62 +23,21 @@ namespace math { * @return The value of the derivative */ template -inline double conv_gaus_line(double t0, double t1, double a, double b, - const Tx& x, double sig2) { +inline return_type_t conv_gaus_line(double t0, double t1, double a, double b, + const Tx& x, double sig2) { using stan::math::normal_cdf; using std::exp; using std::pow; using std::sqrt; - const double pi = stan::math::pi(); const double sig = sqrt(sig2); - - double y - = (a * value_of(x) + b) - * (normal_cdf(t1, value_of(x), sig) - normal_cdf(t0, value_of(x), sig)); - y += -a * sig2 / sqrt(2 * pi * sig2) - * (exp(-pow(t1 - value_of(x), 2) / (2 * sig2)) - - exp(-pow(t0 - value_of(x), 2) / (2 * sig2))); - return y; -} - - -/** - * Evaluate the convolution of a piecewise linear function with a - * Gaussian kernel. - * - * \f$\int_{x_0}^{x_1} (a_1t + b_1) e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ - * \f$ + ... + $\int_{x_{n-1}}^{x_{n}} (a_{n-1}t + b_{n-1}) \f$ - * \f$ e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ - * - * @param x point at which convolution is evaluated - * @param xs increasing vector with endpoints of each piecewise linear function - * @param params vector containing slopes and intercepts of peicewise linear function and width of Gaussian kernel - * @param ind_start - * @param ind_end - * @return The value of the derivative - */ -template -inline double conv_gaus_line_sum(const Tx& x, std::vector xs, - std::vector params, - int ind_start, - int ind_end) { - using stan::math::normal_cdf; - using std::exp; - using std::pow; - using std::sqrt; const double pi = stan::math::pi(); - int n = xs.size(); - double sig2 = params[2*n-2]; - double sig = std::sqrt(sig2); - double y=0; - for (int i = ind_start; i < ind_end; i++) { - y += (params[i] * value_of(x) + params[n-1 + i]) - * (normal_cdf(xs[i+1], value_of(x), sig) - normal_cdf(xs[i], value_of(x), sig)); - y += -params[i] * sig2 / sqrt(2 * pi * sig2) - * (exp(-pow(xs[i+1] - value_of(x), 2) / (2 * sig2)) - - exp(-pow(xs[i] - value_of(x), 2) / (2 * sig2))); - } + return_type_t y + = (a * x + b) + * (normal_cdf(t1, x, sig) - normal_cdf(t0, x, sig)); + y += -a * sig2 / sqrt(2 * pi * sig2) + * (exp(-pow(t1 - x, 2) / (2 * sig2)) + - exp(-pow(t0 - x, 2) / (2 * sig2))); return y; } diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index cc10df09d12..bd77953c2e2 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -92,15 +92,11 @@ inline return_type_t gaus_interp(const std::vector& xs, ind_end = std::min(n - 1, ind_end); // sum convolutions over intervals - /* return_type_t y = 0; for (int i = ind_start; i < ind_end; i++) { y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(n-1) + i], x, params[2*n-2]); } - */ - return_type_t y = conv_gaus_line_sum(x, xs2, params, ind_start, ind_end); - return y; } @@ -131,6 +127,8 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, using internal::min_diff; static const double max_diff = 1e-8; + // set Gaussian kernel to sig2_scale times smallest difference between + // successive points static const double sig2_scale = 0.1; int n = xs.size(); diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index 0c5dee69e70..c6a5d26236d 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -1,8 +1,8 @@ #include TEST(mathMixGausInterp, derivs) { - using stan::math::gaus_interp_vect; using stan::math::var; + using stan::math::gaus_interp; std::vector xs, ys, ts; double xmin, xmax, x, y, t, t0, t1, dder, dder2, x0, x1; int n; @@ -40,11 +40,11 @@ TEST(mathMixGausInterp, derivs) { } std::vector ys_new, ys_new_p, ys_new_n, xs_new_p, xs_new_n; - ys_new = gaus_interp_vect(xs, ys, xs_new); + ys_new = gaus_interp(xs, ys, xs_new); // autodiff at each interpolation pt std::vector ys_new_v; - ys_new_v = gaus_interp_vect(xs, ys, xs_new_v); + ys_new_v = gaus_interp(xs, ys, xs_new_v); std::vector ys_new_dder; for (int i = 0; i < n_interp; i++) { @@ -54,14 +54,14 @@ TEST(mathMixGausInterp, derivs) { } // take derivative of interpolation using finite differencing - double h = 1e-8; + double h = 1e-6; xs_new_p = xs_new; xs_new_n = xs_new; for (int i = 0; i < n_interp; i++) { xs_new_p[i] += -h; xs_new_n[i] += h; - ys_new_p = gaus_interp_vect(xs, ys, xs_new_p); - ys_new_n = gaus_interp_vect(xs, ys, xs_new_n); + ys_new_p = gaus_interp(xs, ys, xs_new_p); + ys_new_n = gaus_interp(xs, ys, xs_new_n); dder2 = (ys_new_n[i] - ys_new_p[i]) / (2 * h); ASSERT_NEAR(ys_new_dder[i], dder2, 1e-5); } From 57c7b0da1d7bcdd4fe387d2f5d71a68994750fec Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 3 Sep 2020 12:36:13 +0000 Subject: [PATCH 50/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/conv_gaus_line.hpp | 7 +++---- stan/math/prim/fun/gaus_interp.hpp | 6 +++--- test/unit/math/mix/fun/gaus_interp_test.cpp | 2 +- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 9ae581d583f..25c72e5414d 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -23,8 +23,8 @@ namespace math { * @return The value of the derivative */ template -inline return_type_t conv_gaus_line(double t0, double t1, double a, double b, - const Tx& x, double sig2) { +inline return_type_t conv_gaus_line(double t0, double t1, double a, + double b, const Tx& x, double sig2) { using stan::math::normal_cdf; using std::exp; using std::pow; @@ -33,8 +33,7 @@ inline return_type_t conv_gaus_line(double t0, double t1, double a, double b const double pi = stan::math::pi(); return_type_t y - = (a * x + b) - * (normal_cdf(t1, x, sig) - normal_cdf(t0, x, sig)); + = (a * x + b) * (normal_cdf(t1, x, sig) - normal_cdf(t0, x, sig)); y += -a * sig2 / sqrt(2 * pi * sig2) * (exp(-pow(t1 - x, 2) / (2 * sig2)) - exp(-pow(t0 - x, 2) / (2 * sig2))); diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 0ac9ecd323d..5cc99efdb2a 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -94,8 +94,8 @@ inline return_type_t gaus_interp(const std::vector& xs, // sum convolutions over intervals return_type_t y = 0; for (int i = ind_start; i < ind_end; i++) { - y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(n-1) + i], x, - params[2*n-2]); + y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(n - 1) + i], x, + params[2 * n - 2]); } return y; } @@ -127,7 +127,7 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, using internal::min_diff; static const double max_diff = 1e-8; - // set Gaussian kernel to sig2_scale times smallest difference between + // set Gaussian kernel to sig2_scale times smallest difference between // successive points static const double sig2_scale = 0.1; int n = xs.size(); diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/gaus_interp_test.cpp index c6a5d26236d..a2adf4e2365 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/gaus_interp_test.cpp @@ -1,8 +1,8 @@ #include TEST(mathMixGausInterp, derivs) { - using stan::math::var; using stan::math::gaus_interp; + using stan::math::var; std::vector xs, ys, ts; double xmin, xmax, x, y, t, t0, t1, dder, dder2, x0, x1; int n; From fa5067fb74e6bcccd07cf32879e06e786e39db52 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Tue, 27 Oct 2020 17:37:31 -0400 Subject: [PATCH 51/57] fix issues --- stan/math/prim/fun/conv_gaus_line.hpp | 4 +- stan/math/prim/fun/gaus_interp.hpp | 53 ++++++++++++++------------- stan/math/prim/fun/lin_interp.hpp | 2 +- stan/math/rev/fun/conv_gaus_line.hpp | 10 ++--- 4 files changed, 34 insertions(+), 35 deletions(-) diff --git a/stan/math/prim/fun/conv_gaus_line.hpp b/stan/math/prim/fun/conv_gaus_line.hpp index 9ae581d583f..e2009a337e1 100644 --- a/stan/math/prim/fun/conv_gaus_line.hpp +++ b/stan/math/prim/fun/conv_gaus_line.hpp @@ -14,6 +14,7 @@ namespace math { * * \f$\int_{t_0}^{t_1} (at + b) e^{\frac{-(t-x)^2}{2\sigma^2}} dt \f$ * + * @tparam Tx type of x * @param t0 lower integration bound * @param t1 upper integration bound * @param a coefficient of t in line @@ -30,12 +31,11 @@ inline return_type_t conv_gaus_line(double t0, double t1, double a, double b using std::pow; using std::sqrt; const double sig = sqrt(sig2); - const double pi = stan::math::pi(); return_type_t y = (a * x + b) * (normal_cdf(t1, x, sig) - normal_cdf(t0, x, sig)); - y += -a * sig2 / sqrt(2 * pi * sig2) + y += -a * sig2 / sqrt(2 * stan::math::pi() * sig2) * (exp(-pow(t1 - x, 2) / (2 * sig2)) - exp(-pow(t0 - x, 2) / (2 * sig2))); diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 0ac9ecd323d..9ddff6515bb 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -16,9 +16,10 @@ namespace internal { * find the smallest difference between successive elements in a sorted vector */ template -inline double min_diff(int n, const std::vector& xs) { +inline double min_diff(const std::vector& xs) { double dmin = value_of(xs[1]) - value_of(xs[0]); - for (int i = 1; i < n - 1; i++) { + int N = xs.size(); + for (int i = 1; i < N - 1; i++) { if (value_of(xs[i + 1]) - value_of(xs[i]) < dmin) { dmin = value_of(xs[i + 1]) - value_of(xs[i]); } @@ -31,9 +32,9 @@ inline double min_diff(int n, const std::vector& xs) { /** * Given a set of reference points \f$(xs_i, ys_i)\f$, create a mollifier * that intersects the reference points. This function requires as input - * a struct created by the function gaus_interp_precomp. The algorithm - * used to create the mollifier is an iterative algorithm that works - * as follows. First a linear + * a vector, params, created by the function gaus_interp_precomp. The + * algorithm used to create the mollifier is an iterative algorithm that + * works as follows. First a linear * interpolation is created through the reference points. Then, the * linear interpolation is convolved with a Gaussian whose width is * proportional the smallest distance between successive points @@ -45,6 +46,7 @@ inline double min_diff(int n, const std::vector& xs) { * Note: This interpolation scheme should be used when the function * to be interpolated is well-resolved by the reference points. * + * @tparam Tx type of x * @param xs vector of independent variable of reference points * @param ys vector of dependent variable of reference points * @param params a vector created by gaus_interp_precomp @@ -69,7 +71,7 @@ inline return_type_t gaus_interp(const std::vector& xs, // number of standard deviations to extend endpoints for convolution const double NSTDS = 10; - int n = xs.size(); + int N = xs.size(); // params is vector of the form (as, bs, sig2) @@ -77,25 +79,25 @@ inline return_type_t gaus_interp(const std::vector& xs, std::vector xs2 = xs; // extend out first and last lines for convolution - double sig = std::sqrt(params[2 * n - 2]); + double sig = std::sqrt(params[2 * N - 2]); xs2[0] = xs[0] - NSTDS * sig; - xs2[n - 1] = xs[n - 1] + NSTDS * sig; + xs2[N - 1] = xs[N - 1] + NSTDS * sig; // no need to convolve far from center of gaussian, so // get lower and upper indexes for integration bounds auto lb = std::lower_bound(xs.begin(), xs.end(), x - NSTDS * sig); - int ind_start = &(*lb) - &xs[0] - 1; + int ind_start = std::distance(xs.begin(), lb) - 1; ind_start = std::max(0, ind_start); auto ub = std::upper_bound(xs.begin(), xs.end(), x + NSTDS * sig); - int ind_end = &(*ub) - &xs[0]; - ind_end = std::min(n - 1, ind_end); + int ind_end = std::distance(xs.begin(), ub); + ind_end = std::min(N-1, ind_end); // sum convolutions over intervals return_type_t y = 0; for (int i = ind_start; i < ind_end; i++) { - y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(n-1) + i], x, - params[2*n-2]); + y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(N-1) + i], x, + params[2*N-2]); } return y; } @@ -112,10 +114,9 @@ inline return_type_t gaus_interp(const std::vector& xs, * A tolerance for the maximum distance between the interpolation and * all reference points is also set manually and is not an input. * - * * @param xs vector of independent variable of reference points * @param ys vector of dependent variable of reference points - * @return struct containing slopes, intercepts, and width of kernel + * @return vector containing slopes, intercepts, and width of kernel */ inline std::vector gaus_interp_precomp(const std::vector& xs, const std::vector& ys) { @@ -130,34 +131,34 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, // set Gaussian kernel to sig2_scale times smallest difference between // successive points static const double sig2_scale = 0.1; - int n = xs.size(); + int N = xs.size(); // create the vector to be returned that consists of as, bs, sig2 - std::vector params; - params.resize(2 * n - 1); - params[2 * n - 2] = square(min_diff(n, xs) * sig2_scale); + std::vector params(2 * N - 1, 0.0); + params[2 * N - 2] = square(min_diff(xs) * sig2_scale); // copy ys into a new vector that will be changed std::vector y2s = ys; // interatively find interpolation that coincides with ys at xs - int max_iters = 50; + int max_iters = 100; double dd; for (int j = 0; j < max_iters; j++) { // find slope and intercept of line between each point - for (size_t i = 0; i < n - 1; i++) { + for (size_t i = 0; i < N - 1; i++) { params[i] = (y2s[i + 1] - y2s[i]) / (xs[i + 1] - xs[i]); - params[(n - 1) + i] = -xs[i] * params[i] + y2s[i]; + params[(N - 1) + i] = -xs[i] * params[i] + y2s[i]; } double dmax = 0; - for (size_t i = 0; i < n; i++) { + for (size_t i = 0; i < N; i++) { dd = ys[i] - gaus_interp(xs, y2s, params, xs[i]); y2s[i] += dd; dmax = std::max(std::abs(dd), dmax); } - if (dmax < max_diff) + if (dmax < max_diff) { break; + } } return params; } @@ -166,15 +167,15 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, * This function combines gaus_interp_precomp and gaus_interp. * It takes as input two vectors of reference points (xs and ys) * in addition to a vector, xs_new, of points at which this - * function will evaluate the interpolation through those reference + * function will evaluate the interpolation at those reference * points. * + * @tparam Tx type of xs_new * @param xs vector of independent variable of reference points * @param ys vector of dependent variable of reference points * @param xs_new vector of point at which to evaluate interpolation * @return vector of interpolation values */ - template inline std::vector gaus_interp(const std::vector& xs, const std::vector& ys, diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/lin_interp.hpp index efce23cb855..4315c659d2d 100644 --- a/stan/math/prim/fun/lin_interp.hpp +++ b/stan/math/prim/fun/lin_interp.hpp @@ -43,7 +43,7 @@ inline double lin_interp(const std::vector& xs, // find in between which points the input, x, lives auto ub = std::upper_bound(xs.begin(), xs.end(), x); - int ind = &(*ub) - &(xs[0]); + int ind = std::distance(xs.begin(), ub); // check if the interpolation point falls on a reference point if (x == xs[ind - 1]) { diff --git a/stan/math/rev/fun/conv_gaus_line.hpp b/stan/math/rev/fun/conv_gaus_line.hpp index 0551e9e957e..640b367bcd6 100644 --- a/stan/math/rev/fun/conv_gaus_line.hpp +++ b/stan/math/rev/fun/conv_gaus_line.hpp @@ -34,13 +34,11 @@ inline double der_conv_gaus_line(double t0, double t1, double a, double b, const double sig = sqrt(sig2); const double alpha = sqrt(2 * pi() * sig2); - double y = (a * x0 + b) / alpha - * (-exp(-pow(t1 - x0, 2) / (2 * sig2)) - + exp(-pow(t0 - x0, 2) / (2 * sig2))); + double term1 = exp(-pow(t1 - x0, 2) / (2 * sig2)); + double term2 = exp(-pow(t0 - x0, 2) / (2 * sig2)); + double y = (a * x0 + b) / alpha * (-term1 + term2); y += a * (normal_cdf(t1, x0, sig) - normal_cdf(t0, x0, sig)); - y -= a / alpha - * ((t1 - x0) * exp(-pow(t1 - x0, 2) / (2 * sig2)) - - (t0 - x0) * exp(-pow(t0 - x0, 2) / (2 * sig2))); + y -= a / alpha * ((t1 - x0) * term1 - (t0 - x0) * term2); return y; } From 62ea0833caf5016e904c3ed67a7b17970a0163c1 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 27 Oct 2020 22:14:14 +0000 Subject: [PATCH 52/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/gaus_interp.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index a7e5559b9c3..89637b2bcf0 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -32,8 +32,8 @@ inline double min_diff(const std::vector& xs) { /** * Given a set of reference points \f$(xs_i, ys_i)\f$, create a mollifier * that intersects the reference points. This function requires as input - * a vector, params, created by the function gaus_interp_precomp. The - * algorithm used to create the mollifier is an iterative algorithm that + * a vector, params, created by the function gaus_interp_precomp. The + * algorithm used to create the mollifier is an iterative algorithm that * works as follows. First a linear * interpolation is created through the reference points. Then, the * linear interpolation is convolved with a Gaussian whose width is @@ -91,14 +91,14 @@ inline return_type_t gaus_interp(const std::vector& xs, auto ub = std::upper_bound(xs.begin(), xs.end(), x + NSTDS * sig); int ind_end = std::distance(xs.begin(), ub); - ind_end = std::min(N-1, ind_end); + ind_end = std::min(N - 1, ind_end); // sum convolutions over intervals return_type_t y = 0; for (int i = ind_start; i < ind_end; i++) { <<<<<<< HEAD - y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(N-1) + i], x, - params[2*N-2]); + y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(N - 1) + i], x, + params[2 * N - 2]); ======= y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(n - 1) + i], x, params[2 * n - 2]); From b5f69eb53423b358915dc177757e97103fc6c922 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 28 Oct 2020 07:31:49 -0400 Subject: [PATCH 53/57] merge conflicts --- stan/math/prim/fun/gaus_interp.hpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index a7e5559b9c3..10e877bc284 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -96,13 +96,8 @@ inline return_type_t gaus_interp(const std::vector& xs, // sum convolutions over intervals return_type_t y = 0; for (int i = ind_start; i < ind_end; i++) { -<<<<<<< HEAD y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(N-1) + i], x, params[2*N-2]); -======= - y += conv_gaus_line(xs2[i], xs2[i + 1], params[i], params[(n - 1) + i], x, - params[2 * n - 2]); ->>>>>>> 57c7b0da1d7bcdd4fe387d2f5d71a68994750fec } return y; } From a992cad820bd5e855cab2876391fe2f6648aac02 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Wed, 11 Nov 2020 13:51:52 -0500 Subject: [PATCH 54/57] rerun tests --- stan/math/prim/fun/gaus_interp.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 684b70d0243..3ec373591b1 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -34,8 +34,8 @@ inline double min_diff(const std::vector& xs) { * that intersects the reference points. This function requires as input * a vector, params, created by the function gaus_interp_precomp. The * algorithm used to create the mollifier is an iterative algorithm that - * works as follows. First a linear - * interpolation is created through the reference points. Then, the + * works as follows. First a linear interpolation is created through the + * reference points. Then, the * linear interpolation is convolved with a Gaussian whose width is * proportional the smallest distance between successive points * \f$xs_i\f$ and \f$xs_{i+1}\f$. Since the convolution is unlikely to From 57c2f8bf298f69c6dcd2ff287f9e4faa4b40a173 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 11 Nov 2020 21:05:02 +0000 Subject: [PATCH 55/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/gaus_interp.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/gaus_interp.hpp index 3ec373591b1..f18fd7ac9ac 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/gaus_interp.hpp @@ -34,7 +34,7 @@ inline double min_diff(const std::vector& xs) { * that intersects the reference points. This function requires as input * a vector, params, created by the function gaus_interp_precomp. The * algorithm used to create the mollifier is an iterative algorithm that - * works as follows. First a linear interpolation is created through the + * works as follows. First a linear interpolation is created through the * reference points. Then, the * linear interpolation is convolved with a Gaussian whose width is * proportional the smallest distance between successive points From 2cbbad1f5f1797ec70fb65c682c899f94b0901d6 Mon Sep 17 00:00:00 2001 From: Philip Greengard Date: Tue, 21 Sep 2021 19:58:15 -0400 Subject: [PATCH 56/57] renaming --- stan/math/prim/fun.hpp | 4 +- .../fun/{gaus_interp.hpp => interp_gauss.hpp} | 40 ++++++------- .../fun/{lin_interp.hpp => interp_lin.hpp} | 12 ++-- ..._interp_test.cpp => interp_gauss_test.cpp} | 10 ++-- ..._interp_test.cpp => interp_gauss_test.cpp} | 60 +++++++++---------- ...in_interp_test.cpp => interp_lin_test.cpp} | 32 +++++----- 6 files changed, 79 insertions(+), 79 deletions(-) rename stan/math/prim/fun/{gaus_interp.hpp => interp_gauss.hpp} (83%) rename stan/math/prim/fun/{lin_interp.hpp => interp_lin.hpp} (89%) rename test/unit/math/mix/fun/{gaus_interp_test.cpp => interp_gauss_test.cpp} (87%) rename test/unit/math/prim/fun/{gaus_interp_test.cpp => interp_gauss_test.cpp} (66%) rename test/unit/math/prim/fun/{lin_interp_test.cpp => interp_lin_test.cpp} (72%) diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index e92ad058463..1cec9e8b746 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -103,7 +103,6 @@ #include #include #include -#include #include #include #include @@ -130,6 +129,8 @@ #include #include #include +#include +#include #include #include #include @@ -155,7 +156,6 @@ #include #include #include -#include #include #include #include diff --git a/stan/math/prim/fun/gaus_interp.hpp b/stan/math/prim/fun/interp_gauss.hpp similarity index 83% rename from stan/math/prim/fun/gaus_interp.hpp rename to stan/math/prim/fun/interp_gauss.hpp index 3ec373591b1..3bc4579e436 100644 --- a/stan/math/prim/fun/gaus_interp.hpp +++ b/stan/math/prim/fun/interp_gauss.hpp @@ -1,5 +1,5 @@ -#ifndef STAN_MATH_PRIM_FUN_GAUS_INTERP -#define STAN_MATH_PRIM_FUN_GAUS_INTERP +#ifndef STAN_MATH_PRIM_FUN_INTERP_GAUSS +#define STAN_MATH_PRIM_FUN_INTERP_GAUSS #include #include @@ -32,7 +32,7 @@ inline double min_diff(const std::vector& xs) { /** * Given a set of reference points \f$(xs_i, ys_i)\f$, create a mollifier * that intersects the reference points. This function requires as input - * a vector, params, created by the function gaus_interp_precomp. The + * a vector, params, created by the function interp_gauss_precomp. The * algorithm used to create the mollifier is an iterative algorithm that * works as follows. First a linear interpolation is created through the * reference points. Then, the @@ -49,18 +49,18 @@ inline double min_diff(const std::vector& xs) { * @tparam Tx type of x * @param xs vector of independent variable of reference points * @param ys vector of dependent variable of reference points - * @param params a vector created by gaus_interp_precomp + * @param params a vector created by interp_gauss_precomp * @param x the point at which to evaluate the interpolation * @return value of the interpolation at x */ template -inline return_type_t gaus_interp(const std::vector& xs, - const std::vector& ys, - const std::vector& params, - const Tx& x) { +inline return_type_t interp_gauss(const std::vector& xs, + const std::vector& ys, + const std::vector& params, + const Tx& x) { // enforce that interpolation point is between smallest and largest // reference point - static char const* function = "gaus_interp"; + static char const* function = "interp_gauss"; check_less_or_equal(function, "Interpolation point", x, xs.back()); check_greater_or_equal(function, "Interpolation point", x, xs.front()); check_ordered(function, "xs", xs); @@ -103,7 +103,7 @@ inline return_type_t gaus_interp(const std::vector& xs, } /** - * This function was written to be used with gaus_interp. This function + * This function was written to be used with interp_gauss. This function * computes the shifted y-values of the reference points of an interpolation * in such a way that when that piecewise linear function is convolved * with a Gaussian kernel, the resulting function coincides with the @@ -118,9 +118,9 @@ inline return_type_t gaus_interp(const std::vector& xs, * @param ys vector of dependent variable of reference points * @return vector containing slopes, intercepts, and width of kernel */ -inline std::vector gaus_interp_precomp(const std::vector& xs, - const std::vector& ys) { - static char const* function = "gaus_interp_precomp"; +inline std::vector interp_gauss_precomp(const std::vector& xs, + const std::vector& ys) { + static char const* function = "interp_gauss_precomp"; check_not_nan(function, "xs", xs); check_not_nan(function, "ys", ys); check_ordered(function, "xs", xs); @@ -152,7 +152,7 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, double dmax = 0; for (size_t i = 0; i < N; i++) { - dd = ys[i] - gaus_interp(xs, y2s, params, xs[i]); + dd = ys[i] - interp_gauss(xs, y2s, params, xs[i]); y2s[i] += dd; dmax = std::max(std::abs(dd), dmax); } @@ -164,7 +164,7 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, } /** - * This function combines gaus_interp_precomp and gaus_interp. + * This function combines interp_gauss_precomp and interp_gauss. * It takes as input two vectors of reference points (xs and ys) * in addition to a vector, xs_new, of points at which this * function will evaluate the interpolation at those reference @@ -177,16 +177,16 @@ inline std::vector gaus_interp_precomp(const std::vector& xs, * @return vector of interpolation values */ template -inline std::vector gaus_interp(const std::vector& xs, - const std::vector& ys, - const std::vector& xs_new) { +inline std::vector interp_gauss(const std::vector& xs, + const std::vector& ys, + const std::vector& xs_new) { int n_interp = xs_new.size(); std::vector ys_new(n_interp); // create interpolation - std::vector params = gaus_interp_precomp(xs, ys); + std::vector params = interp_gauss_precomp(xs, ys); for (int i = 0; i < n_interp; i++) { - ys_new[i] = gaus_interp(xs, ys, params, xs_new[i]); + ys_new[i] = interp_gauss(xs, ys, params, xs_new[i]); } return ys_new; } diff --git a/stan/math/prim/fun/lin_interp.hpp b/stan/math/prim/fun/interp_lin.hpp similarity index 89% rename from stan/math/prim/fun/lin_interp.hpp rename to stan/math/prim/fun/interp_lin.hpp index 4315c659d2d..1decb8f36a5 100644 --- a/stan/math/prim/fun/lin_interp.hpp +++ b/stan/math/prim/fun/interp_lin.hpp @@ -1,5 +1,5 @@ -#ifndef STAN_MATH_PRIM_FUN_LIN_INTERP -#define STAN_MATH_PRIM_FUN_LIN_INTERP +#ifndef STAN_MATH_PRIM_FUN_INTERP_LIN +#define STAN_MATH_PRIM_FUN_INTERP_LIN #include #include @@ -23,9 +23,9 @@ namespace math { * @param x the point at which to evaluate the interpolation * @return value of linear interpolation at x */ -inline double lin_interp(const std::vector& xs, +inline double interp_lin(const std::vector& xs, const std::vector& ys, double x) { - static char const* function = "lin_interp"; + static char const* function = "interp_lin"; check_ordered(function, "xs", xs); check_not_nan(function, "xs", xs); check_not_nan(function, "ys", ys); @@ -72,7 +72,7 @@ inline double lin_interp(const std::vector& xs, */ template -inline std::vector lin_interp(const std::vector& xs, +inline std::vector interp_lin(const std::vector& xs, const std::vector& ys, const std::vector& xs_new) { int n_interp = xs_new.size(); @@ -80,7 +80,7 @@ inline std::vector lin_interp(const std::vector& xs, // create interpolation for (int i = 0; i < n_interp; i++) { - ys_new[i] = lin_interp(xs, ys, xs_new[i]); + ys_new[i] = interp_lin(xs, ys, xs_new[i]); } return ys_new; } diff --git a/test/unit/math/mix/fun/gaus_interp_test.cpp b/test/unit/math/mix/fun/interp_gauss_test.cpp similarity index 87% rename from test/unit/math/mix/fun/gaus_interp_test.cpp rename to test/unit/math/mix/fun/interp_gauss_test.cpp index a2adf4e2365..456ff0c8f92 100644 --- a/test/unit/math/mix/fun/gaus_interp_test.cpp +++ b/test/unit/math/mix/fun/interp_gauss_test.cpp @@ -1,7 +1,7 @@ #include TEST(mathMixGausInterp, derivs) { - using stan::math::gaus_interp; + using stan::math::interp_gauss; using stan::math::var; std::vector xs, ys, ts; double xmin, xmax, x, y, t, t0, t1, dder, dder2, x0, x1; @@ -40,11 +40,11 @@ TEST(mathMixGausInterp, derivs) { } std::vector ys_new, ys_new_p, ys_new_n, xs_new_p, xs_new_n; - ys_new = gaus_interp(xs, ys, xs_new); + ys_new = interp_gauss(xs, ys, xs_new); // autodiff at each interpolation pt std::vector ys_new_v; - ys_new_v = gaus_interp(xs, ys, xs_new_v); + ys_new_v = interp_gauss(xs, ys, xs_new_v); std::vector ys_new_dder; for (int i = 0; i < n_interp; i++) { @@ -60,8 +60,8 @@ TEST(mathMixGausInterp, derivs) { for (int i = 0; i < n_interp; i++) { xs_new_p[i] += -h; xs_new_n[i] += h; - ys_new_p = gaus_interp(xs, ys, xs_new_p); - ys_new_n = gaus_interp(xs, ys, xs_new_n); + ys_new_p = interp_gauss(xs, ys, xs_new_p); + ys_new_n = interp_gauss(xs, ys, xs_new_n); dder2 = (ys_new_n[i] - ys_new_p[i]) / (2 * h); ASSERT_NEAR(ys_new_dder[i], dder2, 1e-5); } diff --git a/test/unit/math/prim/fun/gaus_interp_test.cpp b/test/unit/math/prim/fun/interp_gauss_test.cpp similarity index 66% rename from test/unit/math/prim/fun/gaus_interp_test.cpp rename to test/unit/math/prim/fun/interp_gauss_test.cpp index 3571792b246..efdce502dfc 100644 --- a/test/unit/math/prim/fun/gaus_interp_test.cpp +++ b/test/unit/math/prim/fun/interp_gauss_test.cpp @@ -4,9 +4,9 @@ double ABS_TOL = 1e-12; -TEST(mathPrimGausInterp, throwing) { - using stan::math::gaus_interp; - using stan::math::gaus_interp_precomp; +TEST(mathPrimInterpGauss, throwing) { + using stan::math::interp_gauss; + using stan::math::interp_gauss_precomp; double nan = std::numeric_limits::quiet_NaN(); double x; std::vector xs, ys; @@ -15,55 +15,55 @@ TEST(mathPrimGausInterp, throwing) { int n = 2; xs = {1, 1}; ys = {0, 2}; - EXPECT_THROW(gaus_interp_precomp(xs, ys), std::domain_error); + EXPECT_THROW(interp_gauss_precomp(xs, ys), std::domain_error); // check that when xs contain a nan xs = {nan, 1}; ys = {0, 2}; - EXPECT_THROW(gaus_interp_precomp(xs, ys), std::domain_error); + EXPECT_THROW(interp_gauss_precomp(xs, ys), std::domain_error); // xs must contain at least two elements xs = {1}; ys = {0, 2}; - EXPECT_THROW(gaus_interp_precomp(xs, ys), std::domain_error); + EXPECT_THROW(interp_gauss_precomp(xs, ys), std::domain_error); // check that error throws when trying to interpolate out of range or nan xs = {0, 1}; ys = {0, 2}; - std::vector params = gaus_interp_precomp(xs, ys); + std::vector params = interp_gauss_precomp(xs, ys); x = 1.1; - EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + EXPECT_THROW(interp_gauss(xs, ys, params, x), std::domain_error); x = -0.1; - EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + EXPECT_THROW(interp_gauss(xs, ys, params, x), std::domain_error); x = nan; - EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + EXPECT_THROW(interp_gauss(xs, ys, params, x), std::domain_error); // ys can't contain nan xs = {0, 1}; ys = {0, nan}; x = 0.5; - EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + EXPECT_THROW(interp_gauss(xs, ys, params, x), std::domain_error); // xs can't contain nan xs = {0, nan}; ys = {0, 2}; x = 0.5; - EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + EXPECT_THROW(interp_gauss(xs, ys, params, x), std::domain_error); // xs must be increasing xs = {1, 1}; ys = {0, 2}; - EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + EXPECT_THROW(interp_gauss(xs, ys, params, x), std::domain_error); // xs must contain at least two elements xs = {1}; ys = {0, 2}; - EXPECT_THROW(gaus_interp(xs, ys, params, x), std::domain_error); + EXPECT_THROW(interp_gauss(xs, ys, params, x), std::domain_error); } -TEST(mathPrimGausInterp, interp_line) { - using stan::math::gaus_interp; - using stan::math::gaus_interp_precomp; +TEST(mathPrimInterpGauss, interp_line) { + using stan::math::interp_gauss; + using stan::math::interp_gauss_precomp; // check that interpolation of line returns the same function // generate function tabulation @@ -86,13 +86,13 @@ TEST(mathPrimGausInterp, interp_line) { // create interpolation using precomp std::vector ys_new_gaus(n_interp); - std::vector params = gaus_interp_precomp(xs, ys); + std::vector params = interp_gauss_precomp(xs, ys); for (int i = 0; i < n_interp; i++) { - ys_new_gaus[i] = gaus_interp(xs, ys, params, xs_new[i]); + ys_new_gaus[i] = interp_gauss(xs, ys, params, xs_new[i]); } // create interpolation without precomp - std::vector ys_new_gaus2 = gaus_interp(xs, ys, xs_new); + std::vector ys_new_gaus2 = interp_gauss(xs, ys, xs_new); // test points double tmp, y; @@ -104,8 +104,8 @@ TEST(mathPrimGausInterp, interp_line) { } } -TEST(mathPrimGausInterp, matching_reference_interp_pts) { - using stan::math::gaus_interp; +TEST(mathPrimInterpGauss, matching_reference_interp_pts) { + using stan::math::interp_gauss; // check that interpolation returns the same function // when interpolation points are the same as reference points @@ -119,7 +119,7 @@ TEST(mathPrimGausInterp, matching_reference_interp_pts) { std::vector xs_new = xs; // create interpolation - std::vector ys_new = gaus_interp(xs, ys, xs_new); + std::vector ys_new = interp_gauss(xs, ys, xs_new); // test points for (int i = 0; i < n; i++) { @@ -127,10 +127,10 @@ TEST(mathPrimGausInterp, matching_reference_interp_pts) { } } -TEST(mathPrimGausInterp, gaus_and_lin_interp) { - using stan::math::gaus_interp; - using stan::math::gaus_interp_precomp; - using stan::math::lin_interp; +TEST(mathPrimInterpGauss, interp_gauss_and_lin) { + using stan::math::interp_gauss; + using stan::math::interp_gauss_precomp; + using stan::math::interp_lin; // check that interpolation of line returns the same function // generate function tabulation @@ -163,13 +163,13 @@ TEST(mathPrimGausInterp, gaus_and_lin_interp) { // linear interpolation ys_new_lin.resize(n_interp); for (int i = 0; i < n_interp; i++) { - ys_new_lin[i] = lin_interp(xs, ys, xs_new[i]); + ys_new_lin[i] = interp_lin(xs, ys, xs_new[i]); } // gaus interpolation - std::vector params = gaus_interp_precomp(xs, ys); + std::vector params = interp_gauss_precomp(xs, ys); for (int i = 0; i < n_interp; i++) { - ys_new_gaus[i] = gaus_interp(xs, ys, params, xs_new[i]); + ys_new_gaus[i] = interp_gauss(xs, ys, params, xs_new[i]); } // test points diff --git a/test/unit/math/prim/fun/lin_interp_test.cpp b/test/unit/math/prim/fun/interp_lin_test.cpp similarity index 72% rename from test/unit/math/prim/fun/lin_interp_test.cpp rename to test/unit/math/prim/fun/interp_lin_test.cpp index a19590ec62f..98f6a2fd240 100644 --- a/test/unit/math/prim/fun/lin_interp_test.cpp +++ b/test/unit/math/prim/fun/interp_lin_test.cpp @@ -4,8 +4,8 @@ double ABS_TOL = 1e-12; -TEST(mathPrimLinInterp, throwing) { - using stan::math::lin_interp; +TEST(mathPrimInterpLin, throwing) { + using stan::math::interp_lin; double nan = std::numeric_limits::quiet_NaN(); double x = 0.5; std::vector xs, ys; @@ -14,33 +14,33 @@ TEST(mathPrimLinInterp, throwing) { int n = 2; xs = {1, 1}; ys = {0, 2}; - EXPECT_THROW(lin_interp(xs, ys, x), std::domain_error); + EXPECT_THROW(interp_lin(xs, ys, x), std::domain_error); // check when xs contain a nan xs = {nan, 1}; ys = {0, 2}; - EXPECT_THROW(lin_interp(xs, ys, x), std::domain_error); + EXPECT_THROW(interp_lin(xs, ys, x), std::domain_error); // xs must contain at least two elements xs = {1}; ys = {0, 2}; - EXPECT_THROW(lin_interp(xs, ys, x), std::domain_error); + EXPECT_THROW(interp_lin(xs, ys, x), std::domain_error); // ys can't contain nan xs = {0, 1}; ys = {0, nan}; x = 0.5; - EXPECT_THROW(lin_interp(xs, ys, x), std::domain_error); + EXPECT_THROW(interp_lin(xs, ys, x), std::domain_error); // x can't be nan xs = {0, 1}; ys = {0, 2}; x = nan; - EXPECT_THROW(lin_interp(xs, ys, x), std::domain_error); + EXPECT_THROW(interp_lin(xs, ys, x), std::domain_error); } -TEST(mathPrimLinInterp, interp_line) { - using stan::math::lin_interp; +TEST(mathPrimInterpLin, interp_line) { + using stan::math::interp_lin; // check that interpolation of line returns the same function // generate function tabulation @@ -64,7 +64,7 @@ TEST(mathPrimLinInterp, interp_line) { // interpolate std::vector ys_new(n_interp); for (int i = 0; i < n_interp; i++) { - ys_new[i] = lin_interp(xs, ys, xs_new[i]); + ys_new[i] = interp_lin(xs, ys, xs_new[i]); } // test points @@ -76,18 +76,18 @@ TEST(mathPrimLinInterp, interp_line) { } // check values outside of range of reference points - ASSERT_NEAR(lin_interp(xs, ys, -1), ys[0], ABS_TOL); - ASSERT_NEAR(lin_interp(xs, ys, 100), ys[1], ABS_TOL); + ASSERT_NEAR(interp_lin(xs, ys, -1), ys[0], ABS_TOL); + ASSERT_NEAR(interp_lin(xs, ys, 100), ys[1], ABS_TOL); // xs with more than 2 elements std::vector xs2 = {0, 1, 2, 3, 4, 5, 6}; std::vector ys2 = {0, 2, 5, 2, 3, 2, 2}; double x = 0.5; - ASSERT_NEAR(lin_interp(xs, ys, x), lin_interp(xs2, ys2, x), ABS_TOL); + ASSERT_NEAR(interp_lin(xs, ys, x), interp_lin(xs2, ys2, x), ABS_TOL); } -TEST(mathPrimLinInterp, matching_reference_interp_pts) { - using stan::math::lin_interp; +TEST(mathPrimInterpLin, matching_reference_interp_pts) { + using stan::math::interp_lin; // check that interpolation returns the same function // when interpolation points are the same as reference points @@ -101,7 +101,7 @@ TEST(mathPrimLinInterp, matching_reference_interp_pts) { std::vector xs_new = xs; // create interpolation - std::vector ys_new = lin_interp(xs, ys, xs_new); + std::vector ys_new = interp_lin(xs, ys, xs_new); // test points for (int i = 0; i < n; i++) { From c19300c66a443e67ebb19f922e9bdd88ff3bfff4 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 28 Sep 2021 11:21:59 +0000 Subject: [PATCH 57/57] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/interp_gauss.hpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/fun/interp_gauss.hpp b/stan/math/prim/fun/interp_gauss.hpp index 70dfc813016..5f39ff5fd0a 100644 --- a/stan/math/prim/fun/interp_gauss.hpp +++ b/stan/math/prim/fun/interp_gauss.hpp @@ -55,9 +55,9 @@ inline double min_diff(const std::vector& xs) { */ template inline return_type_t interp_gauss(const std::vector& xs, - const std::vector& ys, - const std::vector& params, - const Tx& x) { + const std::vector& ys, + const std::vector& params, + const Tx& x) { // enforce that interpolation point is between smallest and largest // reference point static char const* function = "interp_gauss"; @@ -119,7 +119,7 @@ inline return_type_t interp_gauss(const std::vector& xs, * @return vector containing slopes, intercepts, and width of kernel */ inline std::vector interp_gauss_precomp(const std::vector& xs, - const std::vector& ys) { + const std::vector& ys) { static char const* function = "interp_gauss_precomp"; check_not_nan(function, "xs", xs); check_not_nan(function, "ys", ys); @@ -178,8 +178,8 @@ inline std::vector interp_gauss_precomp(const std::vector& xs, */ template inline std::vector interp_gauss(const std::vector& xs, - const std::vector& ys, - const std::vector& xs_new) { + const std::vector& ys, + const std::vector& xs_new) { int n_interp = xs_new.size(); std::vector ys_new(n_interp);