[WIP] Tools for testing identities

stan/math/rev/fun/log_sum_exp.hpp

@@ -8,6 +8,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/inv_logit.hpp>
 #include <stan/math/prim/fun/log_sum_exp.hpp>
 #include <cmath>
 #include <vector>
@@ -22,8 +23,10 @@ class log_sum_exp_vv_vari : public op_vv_vari {
   log_sum_exp_vv_vari(vari* avi, vari* bvi)
       : op_vv_vari(log_sum_exp(avi->val_, bvi->val_), avi, bvi) {}
   void chain() {
-    avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_);
-    bvi_->adj_ += adj_ * calculate_chain(bvi_->val_, val_);
+    // avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_);
+    // bvi_->adj_ += adj_ * calculate_chain(bvi_->val_, val_);
+    avi_->adj_ += adj_ * inv_logit(avi_->val_ - bvi_->val_);
+    bvi_->adj_ += adj_ * inv_logit(bvi_->val_ - avi_->val_);
   }
 };
 class log_sum_exp_vd_vari : public op_vd_vari {

stan/math/rev/functor.hpp

@@ -20,5 +20,6 @@
 #include <stan/math/rev/functor/kinsol_solve.hpp>
 #include <stan/math/rev/functor/map_rect_concurrent.hpp>
 #include <stan/math/rev/functor/map_rect_reduce.hpp>
+#include <stan/math/rev/functor/pack_params.hpp>
 
 #endif

stan/math/rev/functor/pack_params.hpp

@@ -0,0 +1,47 @@
+#ifndef STAN_MATH_REV_FUNCTOR_PACK_PARAMS_HPP
+#define STAN_MATH_REV_FUNCTOR_PACK_PARAMS_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stdexcept>
+
+namespace stan {
+namespace math {
+
+template <size_t N, typename F, typename... Ts>
+class invoke_packed {
+ public:
+  static var invoke(const F& f,
+                    const Eigen::Matrix<var, Eigen::Dynamic, 1>& x_arr,
+                    Ts... xs) {
+    return invoke_packed<N - 1, F, const var&, Ts...>::invoke(
+        f, x_arr, x_arr(x_arr.size() - N), xs...);
+  }
+};
+
+template <typename F, typename... Ts>
+class invoke_packed<0, F, Ts...> {
+ public:
+  static var invoke(const F& f, const Eigen::Matrix<var, Eigen::Dynamic, 1>&,
+                    Ts... xs) {
+    return f(xs...);
+  }
+};
+
+template <size_t N, typename F>
+class pack_params {
+ public:
+  explicit pack_params(const F& f) : f_(f) {}
+
+  var operator()(const Eigen::Matrix<var, Eigen::Dynamic, 1>& x) const {
+    return invoke_packed<N, F>::invoke(f_, x);
+  }
+
+ private:
+  const F& f_;
+};
+
+}  // namespace math
+}  // namespace stan
+#endif

test/unit/math/ad_tolerances.hpp

@@ -1,6 +1,8 @@
 #ifndef TEST_UNIT_MATH_AD_TOLERANCES_HPP
 #define TEST_UNIT_MATH_AD_TOLERANCES_HPP
 
+#include <test/unit/math/relative_tolerance.hpp>
+
 namespace stan {
 namespace test {
 
@@ -13,34 +15,35 @@ namespace test {
  * if the function is implemented using only forward mode, and end
  * with the quantity being calculated; for example,
  * `hessian_fvar_grad_` is the gradient calculated by the `hessian`
- * function using forward-mode autodiff.
+ * function using forward-mode autodiff. Those get interpreted as
+ * relative tolerances, see `relative_tolerance` class more details.
  *
  * `gradient_val_`: 1e-8;  `gradient_grad_`: 1e-4
  *
  * `gradient_fvar_val_`: 1e-8;  `gradient_fvar_grad_`: 1e-4
  *
- * `hessian_val_` : 1e-8; `hessian_grad_`: 1e-4; `hessian_hessian_`: 1e-3
+ * `hessian_val_` : 1e-8; `hessian_grad_`: 1e-4; `hessian_hessian_`: (1e-4,1e-3)
  *
  * `hessian_fvar_val_` : 1e-8; `hessian_fvar_grad_`: 1e-4;
- * `hessian_fvar_hessian_`: 1e-3
+ * `hessian_fvar_hessian_`: (1e-4,1e-3)
  *
  * `grad_hessian_val_` : 1e-8; `grad_hessian_hessian_`: 1e-3;
  * `grad_hessian_grad_hessian_`: 1e-2
  */
 struct ad_tolerances {
-  double gradient_val_;
-  double gradient_grad_;
-  double gradient_fvar_val_;
-  double gradient_fvar_grad_;
-  double hessian_val_;
-  double hessian_grad_;
-  double hessian_hessian_;
-  double hessian_fvar_val_;
-  double hessian_fvar_grad_;
-  double hessian_fvar_hessian_;
-  double grad_hessian_val_;
-  double grad_hessian_hessian_;
-  double grad_hessian_grad_hessian_;
+  relative_tolerance gradient_val_;
+  relative_tolerance gradient_grad_;
+  relative_tolerance gradient_fvar_val_;
+  relative_tolerance gradient_fvar_grad_;
+  relative_tolerance hessian_val_;
+  relative_tolerance hessian_grad_;
+  relative_tolerance hessian_hessian_;
+  relative_tolerance hessian_fvar_val_;
+  relative_tolerance hessian_fvar_grad_;
+  relative_tolerance hessian_fvar_hessian_;
+  relative_tolerance grad_hessian_val_;
+  relative_tolerance grad_hessian_hessian_;
+  relative_tolerance grad_hessian_grad_hessian_;
   ad_tolerances()
       : gradient_val_(1e-8),
         gradient_grad_(1e-4),
@@ -50,11 +53,11 @@ struct ad_tolerances {
 
         hessian_val_(1e-8),
         hessian_grad_(1e-4),
-        hessian_hessian_(1e-3),
+        hessian_hessian_(1e-4, 1e-3),
 
         hessian_fvar_val_(1e-8),
         hessian_fvar_grad_(1e-4),
-        hessian_fvar_hessian_(1e-3),
+        hessian_fvar_hessian_(1e-4, 1e-3),
 
         grad_hessian_val_(1e-8),
         grad_hessian_hessian_(1e-3),

test/unit/math/expect_near_rel.hpp

@@ -3,6 +3,7 @@
 
 #include <stan/math.hpp>
 #include <gtest/gtest.h>
+#include <test/unit/math/relative_tolerance.hpp>
 #include <string>
 #include <vector>
 
@@ -15,13 +16,7 @@ namespace internal {
  * Test that the specified values are within the specified tolerance
  * on relative error, and if not, fail the embedded google test.
  *
- * <p>Relative error is defined to be the error `u - v` rescaled by the
- * average absolute value,
- * `rel_err(u, v) = (u - v) / (0.5 * (abs(u) * + abs(v))).`
- *
- * <p>If at least one of `u` or `v` is zero, the absolute error is
- * tested at the specified tolerance, because the relative error
- * reduces to a constant.
+ * Uses relative_tolerance::inexact to compute the tolerance.
  *
  * @tparam T1 type of first argument
  * @tparam T2 type of second argument
@@ -32,23 +27,11 @@ namespace internal {
  */
 template <typename T1, typename T2, require_all_stan_scalar_t<T1, T2>...>
 void expect_near_rel_finite(const std::string& msg, const T1& x1, const T2& x2,
-                            double tol = 1e-8) {
-  using stan::math::fabs;
-  // if either arg near zero, must use absolute tolerance as rel tol -> 2
-  if (fabs(x1) < tol || fabs(x2) < tol) {
-    EXPECT_NEAR(x1, x2, tol) << "expect_near_rel_finite(" << x1 << ", " << x2
-                             << ", absolute tolerance = " << tol << ")"
-                             << "; absolute diff = " << fabs(x1 - x2)
-                             << "    in: " << msg << std::endl;
-    return;
-  }
-  auto avg = 0.5 * (fabs(x1) + fabs(x2));
-  auto relative_diff = (x1 - x2) / avg;
-  EXPECT_NEAR(0, relative_diff, tol)
-      << "expect_near_rel_finite(" << x1 << ", " << x2
-      << ", relative tolerance = " << tol << ")"
-      << "; relative diff = " << relative_diff << std::endl
-      << "    in: " << msg << std::endl;
+                            const relative_tolerance tol
+                            = relative_tolerance()) {
+  double tol_val = tol.inexact(x1, x2);
+  EXPECT_NEAR(x1, x2, tol_val)
+      << "expect_near_rel_finite in: " << msg << std::endl;
 }
 
 template <typename EigMat1, typename EigMat2,
@@ -89,7 +72,7 @@ void expect_near_rel_finite(const std::string& msg, const std::vector<T1>& x1,
  */
 template <typename T1, typename T2, require_all_stan_scalar_t<T1, T2>...>
 void expect_near_rel(const std::string& msg, const T1& x1, const T2& x2,
-                     double tol = 1e-8) {
+                     relative_tolerance tol = relative_tolerance()) {
   if (stan::math::is_nan(x1) || stan::math::is_nan(x2))
     EXPECT_TRUE(stan::math::is_nan(x1) && stan::math::is_nan(x2))
         << "expect_near_rel(" << x1 << ", " << x2 << ")" << std::endl
@@ -120,7 +103,7 @@ void expect_near_rel(const std::string& msg, const T1& x1, const T2& x2,
 template <typename EigMat1, typename EigMat2,
           require_all_eigen_t<EigMat1, EigMat2>...>
 void expect_near_rel(const std::string& msg, EigMat1&& x1, EigMat2&& x2,
-                     double tol = 1e-8) {
+                     relative_tolerance tol = relative_tolerance()) {
   EXPECT_EQ(x1.rows(), x2.rows()) << "expect_near_rel (Eigen::Matrix)"
                                   << " rows must be same size."
                                   << " x1.rows() = " << x1.rows()
@@ -141,7 +124,8 @@ void expect_near_rel(const std::string& msg, EigMat1&& x1, EigMat2&& x2,
 
 template <typename T1, typename T2>
 void expect_near_rel(const std::string& msg, const std::vector<T1>& x1,
-                     const std::vector<T2>& x2, double tol = 1e-8) {
+                     const std::vector<T2>& x2,
+                     relative_tolerance tol = relative_tolerance()) {
   EXPECT_EQ(x1.size(), x2.size()) << "expect_near_rel (std::vector):"
                                   << " vectors must be same size."
                                   << " x1.size() = " << x1.size()

test/unit/math/expect_near_rel_test.cpp

@@ -9,10 +9,10 @@ TEST(testUnitMath, ExpectNearRelScalar) {
   expect_near_rel("test A1", 0, 0, 1e-16);
   expect_near_rel("test A2", 0, 0.0, 1e-16);
   expect_near_rel("test A3", 0.0, 0, 1e-16);
-  expect_near_rel("test B1", 0, 1e-8, 1e-5);
-  expect_near_rel("test B2", 0.0, 1e-8, 1e-5);
-  expect_near_rel("test C1", 1e-8, 0, 1e-5);
-  expect_near_rel("test C2", 1e-8, 0.0, 1e-5);
+  expect_near_rel("test B1", 0, 1e-8, 1e-4);
+  expect_near_rel("test B2", 0.0, 1e-8, 1e-4);
+  expect_near_rel("test C1", 1e-8, 0, 1e-4);
+  expect_near_rel("test C2", 1e-8, 0.0, 1e-4);
 
   // non-zero examples
   expect_near_rel("test D1", 1, 1, 1e-16);
@@ -33,6 +33,7 @@ TEST(testUnitMath, ExpectNearRelScalar) {
 TEST(testUnitMath, ExpectNearRelMatrix) {
   using Eigen::Matrix;
   using stan::test::expect_near_rel;
+  using stan::test::relative_tolerance;
   typedef Matrix<double, -1, 1> v_t;
   typedef Matrix<double, 1, -1> rv_t;
   typedef Matrix<double, -1, -1> m_t;
@@ -57,7 +58,7 @@ TEST(testUnitMath, ExpectNearRelMatrix) {
   m_t d2(2, 3);
   d1 << 1, 2, 3, 0, 0, 0 - 1e-8;
   d2 << 1 + 1e-8, 2 - 1e-8, 3, 0, 0 + 1e-8, 0;
-  expect_near_rel("test D", d1, d2, 1e-6);
+  expect_near_rel("test D", d1, d2, relative_tolerance(1e-6, 1e-8));
 
   // these will fail
   // v_t e1(1);
@@ -82,7 +83,7 @@ TEST(testUnitMath, ExpectNearRelVector) {
 
   expect_near_rel("test C", v_t{1, 1, 1}, v_t{1 + 1e-8, 1 - 1e-9, 1}, 1e-6);
 
-  expect_near_rel("test D", v_t{0, 0, 0}, v_t{0, 0 + 1e-6, 0 - 1e-6}, 1e-4);
+  expect_near_rel("test D", v_t{0, 0, 0}, v_t{0, 0 + 9e-9, 0 - 9e-9}, 1e-4);
 
   // ones after here fail
   // expect_near_rel("test E", v_t{1}, v_t{1, 2}, 1e-6);
@@ -91,6 +92,7 @@ TEST(testUnitMath, ExpectNearRelVector) {
 
 TEST(testUnitMath, ExpectNearRelVectorNesting) {
   using stan::test::expect_near_rel;
+  using stan::test::relative_tolerance;
   using std::vector;
   typedef vector<double> v_t;
   typedef vector<v_t> vv_t;
@@ -102,19 +104,19 @@ TEST(testUnitMath, ExpectNearRelVectorNesting) {
   expect_near_rel("test A", vv_t{}, vv_t{}, 1e-10);
 
   expect_near_rel("test B", vv_t{v_t{1, 2, 3}, v_t{0, 0, 0}},
-                  vv_t{v_t{1, 2, 3}, v_t{0, 0 + 1e-6, 0 - 1e-6}}, 1e-5);
+                  vv_t{v_t{1, 2, 3}, v_t{0, 0 + 1e-6, 0 - 1e-6}}, 1e-3);
 
   expect_near_rel(
       "test C",
       vvv_t{vv_t{v_t{1, 2, 3}, v_t{0, 0, 0}}, vv_t{v_t{1, 2, 3}, v_t{0, 0, 0}}},
       vvv_t{vv_t{v_t{1, 2, 3}, v_t{0, 0 + 1e-6, 0 - 1e-6}},
             vv_t{v_t{1, 2, 3}, v_t{0, 0 + 1e-6, 0 - 1e-6}}},
-      1e-5);
+      relative_tolerance(1e-5, 1e-6));
 
   ev_t d1(3);
   ev_t d2(3);
   d1 << 1, 0, 0;
-  d2 << 1 + 1e-8, 0 + 1e-8, 0;
+  d2 << 1 + 1e-8, 0 + 1e-12, 0;
   vev_t e1{d1, d2};
   vev_t e2{d2, d1};
   expect_near_rel("test E", e1, e2, 1e-6);

test/unit/math/mix/fun/add_diag_test.cpp

@@ -1,6 +1,8 @@
 #include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/ad_tolerances.hpp>
 
 TEST(MathMixMatFun, addDiag) {
+  using stan::test::relative_tolerance;
   auto f
       = [](const auto& x, const auto& y) { return stan::math::add_diag(x, y); };
 
@@ -22,18 +24,21 @@ TEST(MathMixMatFun, addDiag) {
   Eigen::MatrixXd m23(2, 3);
   m23 << 1, 1, 1, 1, 1, 1;
 
+  stan::test::ad_tolerances tol;
+  tol.hessian_hessian_ = relative_tolerance(1e-4, 1e-3);
+  tol.hessian_fvar_hessian_ = relative_tolerance(1e-4, 1e-3);
   // these are OK calls
-  stan::test::expect_ad(f, m00, d);
-  stan::test::expect_ad(f, m00, v0);
-  stan::test::expect_ad(f, m11, d);
-  stan::test::expect_ad(f, m11, v1);
-  stan::test::expect_ad(f, m22, d);
-  stan::test::expect_ad(f, m22, v2);
-  stan::test::expect_ad(f, m23, d);
-  stan::test::expect_ad(f, m23, v2);
+  stan::test::expect_ad(tol, f, m00, d);
+  stan::test::expect_ad(tol, f, m00, v0);
+  stan::test::expect_ad(tol, f, m11, d);
+  stan::test::expect_ad(tol, f, m11, v1);
+  stan::test::expect_ad(tol, f, m22, d);
+  stan::test::expect_ad(tol, f, m22, v2);
+  stan::test::expect_ad(tol, f, m23, d);
+  stan::test::expect_ad(tol, f, m23, v2);
 
   // these throw
-  stan::test::expect_ad(f, m11, v2);
-  stan::test::expect_ad(f, m22, v1);
-  stan::test::expect_ad(f, m23, v1);
+  stan::test::expect_ad(tol, f, m11, v2);
+  stan::test::expect_ad(tol, f, m22, v1);
+  stan::test::expect_ad(tol, f, m23, v1);
 }

test/unit/math/mix/fun/cov_matrix_constrain_test.cpp

@@ -33,9 +33,10 @@ typename stan::scalar_type<T>::type g3(const T& x) {
 
 template <typename T>
 void expect_cov_matrix_transform(const T& x) {
+  using stan::test::relative_tolerance;
   stan::test::ad_tolerances tols;
-  tols.hessian_hessian_ = 1e-2;
-  tols.hessian_fvar_hessian_ = 1e-2;
+  tols.hessian_hessian_ = relative_tolerance(1e-3, 1e-3);
+  tols.hessian_fvar_hessian_ = relative_tolerance(1e-3, 1e-3);
 
   auto f1 = [](const auto& x) { return g1(x); };
   auto f2 = [](const auto& x) { return g2(x); };

test/unit/math/mix/fun/determinant_test.cpp

@@ -3,11 +3,12 @@
 #include <vector>
 
 TEST(MathMixMatFun, determinant) {
+  using stan::test::relative_tolerance;
   auto f = [](const auto& y) { return stan::math::determinant(y); };
 
   stan::test::ad_tolerances tols;
-  tols.hessian_hessian_ = 1e-2;       // default 1e-3
-  tols.hessian_fvar_hessian_ = 1e-2;  // default 1e-3
+  tols.hessian_hessian_ = relative_tolerance(1e-4, 1e-2);
+  tols.hessian_fvar_hessian_ = relative_tolerance(1e-4, 1e-2);
 
   Eigen::MatrixXd z(0, 0);
 

test/unit/math/mix/fun/diag_post_multiply_test.cpp

@@ -17,6 +17,7 @@ void expect_diag_post_multiply(const Eigen::MatrixXd& a,
 }
 
 TEST(MathMixMatFun, diagPostMultiply) {
+  using stan::test::relative_tolerance;
   // 0 x 0
   Eigen::MatrixXd a00(0, 0);
   Eigen::VectorXd u0(0);
@@ -50,8 +51,8 @@ TEST(MathMixMatFun, diagPostMultiply) {
   expect_diag_post_multiply(a33c, u3c);
 
   stan::test::ad_tolerances tols;
-  tols.hessian_hessian_ = 2e-2;       // default 1e-3
-  tols.hessian_fvar_hessian_ = 2e-2;  // default 1e-3
+  tols.hessian_hessian_ = relative_tolerance(1e-4, 2e-2);
+  tols.hessian_fvar_hessian_ = relative_tolerance(1e-4, 2e-2);
 
   Eigen::MatrixXd a33(3, 3);
   a33 << 1, 10, 100, 1000, 2, -4, 8, -16, 32;