fix finite_diff_stepsize and lower tolerance for AD tests on inv_Phi, mdivide_left_test, and mdivide_right_test

Author: SteveBronder

stan/math/mix/functor/wolfe_line_search.hpp

@@ -287,7 +287,7 @@ template <typename Eval, typename WolfeT, typename Option>
 inline bool check_armijo(const Eval& eval, const WolfeT& prev,
                          const Option& opt) {
   return check_armijo(eval.obj(), prev.obj(), eval.alpha(), prev.dir(), opt);
-};
+}
 
 template <typename Option>
 inline auto check_wolfe_curve(double dir_deriv_next, double dir_deriv_init,
@@ -404,7 +404,7 @@ struct Eval {
   inline const auto& dir() const { return dir_; }
   constexpr Eval(double alpha, double obj, double dir)
       : alpha_(alpha), obj_(obj), dir_(dir) {}
-  constexpr explicit Eval() = default;
+  constexpr Eval() = default;
 };
 
 /**

stan/math/prim/fun/finite_diff_stepsize.hpp

@@ -3,15 +3,12 @@
 
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/fabs.hpp>
+#include <limits>
 #include <cmath>
 
 namespace stan {
 namespace math {
 
-#include <cmath>
-#include <limits>
-#include <algorithm>
-
 namespace internal {
 template <std::size_t StencilOrder = 2, typename T>
 inline constexpr auto eps_root_calc() {
@@ -74,13 +71,15 @@ inline constexpr auto eps_root_calc() {
  *       differentiable) behavior along the perturbed coordinate.
  */
 template <std::size_t StencilOrder = 2, typename T>
-inline T finite_diff_stepsize(T u, T c = T(1)) {
-  const T scale = std::max(T(1), std::abs(u));
-  const T eps_root = internal::eps_root_calc<StencilOrder, T>();
-  const T h = c * eps_root * scale;
+inline auto finite_diff_stepsize(T u_, T c = T(1)) {
+  using fp_t = std::conditional_t<std::is_integral_v<T>, double, T>;
+  const fp_t u = static_cast<fp_t>(u_);
+  const fp_t scale = std::max(fp_t{1}, std::abs(u));
+  const fp_t eps_root = internal::eps_root_calc<StencilOrder, fp_t>();
+  const fp_t h = static_cast<fp_t>(c) * eps_root * scale;
   // Ensure perturbation isn’t rounded away at u.
   if (u + h == u) {
-    const T next = std::nextafter(u, std::numeric_limits<T>::infinity());
+    const fp_t next = std::nextafter(u, std::numeric_limits<fp_t>::infinity());
     return std::max(h, next - u);
   }
   return h;

test/unit/math/prim/fun/inv_Phi_test.cpp

@@ -1,5 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
+#include <array>
 #include <limits>
 
 TEST(MathFunctions, inv_Phi) {
@@ -23,14 +24,24 @@ TEST(MathFunctions, inv_Phi) {
 TEST(MathFunctions, Equal) {
   using stan::math::inv_Phi;
   // test output generated with WolframAlpha
-
-  double p[] = {0.000000001, 0.0000001, 0.00001,   0.001,      0.05, 0.15, 0.25,
+  constexpr std::array p_tails{
+    0.000000001,
+    0.0000001,
+    0.99999,
+    0.9999999,
+    0.999999999
+  };
+  constexpr std::array p{0.00001,   0.001,      0.05, 0.15, 0.25,
                 0.35,        0.45,      0.55,      0.65,       0.75, 0.85, 0.95,
-                0.999,       0.99999,   0.9999999, 0.999999999};
-
-  double exact[]
-      = {-5.9978070150076868715623102049115374195951202210145432633059905935,
-         -5.199337582192816931587347266962336866509737160238716454318950531,
+                0.999};
+  constexpr std::array exact_tails{
+    -5.9978070150076868715623102049115374195951202210145432633059905935,
+    -5.199337582192816931587347266962336866509737160238716454318950531,
+    4.2648907939228246284985246989063446293560532226954907262010508062,
+    5.1993375821928169315873472669623368665097371602387164543189505310,
+    5.9978070150076868715623102049115374195951202210145432633059905935
+  };
+  constexpr std::array exact{
          -4.264890793922824628498524698906344629356053222695490726201050806,
          -3.090232306167813541540399830107379205491008491865808855697171108,
          -1.644853626951472714863848907991632136083195744275322071769672094,
@@ -43,15 +54,16 @@ TEST(MathFunctions, Equal) {
          0.6744897501960817432022270145413071853869044150498618956620937885,
          1.0364333894937895797132440746735033661347405959859176279044548663,
          1.6448536269514727148638489079916321360831957442753220717696720944,
-         3.0902323061678135415403998301073792054910084918658088556971711085,
-         4.2648907939228246284985246989063446293560532226954907262010508062,
-         5.1993375821928169315873472669623368665097371602387164543189505310,
-         5.9978070150076868715623102049115374195951202210145432633059905935};
+         3.0902323061678135415403998301073792054910084918658088556971711085};
 
-  int numValues = sizeof(p) / sizeof(double);
 
-  for (int i = 0; i < numValues; ++i) {
-    EXPECT_NEAR(exact[i], inv_Phi(p[i]), 1.5e-15);
+  for (int i = 0; i < p_tails.size(); ++i) {
+    EXPECT_NEAR(exact_tails[i], inv_Phi(p_tails[i]), 4.6e-9) <<
+     "at p_tails(" << i << ") = " << p_tails[i];
+  }
+  for (int i = 0; i < p.size(); ++i) {
+    EXPECT_NEAR(exact[i], inv_Phi(p[i]), 1.5e-15) <<
+     "at p(" << i << ") = " << p[i];
   }
 }
 

test/unit/math/prim/fun/mdivide_left_test.cpp

@@ -7,7 +7,7 @@ TEST(MathMatrixPrim, mdivide_left_val) {
   Ad << 2.0, 3.0, 5.0, 7.0;
 
   stan::math::matrix_d I = Eigen::MatrixXd::Identity(2, 2);
-  EXPECT_MATRIX_FLOAT_EQ(I, stan::math::mdivide_left(Ad, Ad));
+  EXPECT_MATRIX_NEAR(I, stan::math::mdivide_left(Ad, Ad), 1e-15);
 }
 
 TEST(MathMatrixPrim, mdivide_left_val2) {

test/unit/math/prim/fun/mdivide_right_test.cpp

@@ -7,7 +7,7 @@ TEST(MathMatrixPrim, mdivide_right_val) {
   Ad << 2.0, 3.0, 5.0, 7.0;
 
   stan::math::matrix_d I = Eigen::MatrixXd::Identity(2, 2);
-  EXPECT_MATRIX_FLOAT_EQ(I, stan::math::mdivide_left(Ad, Ad));
+  EXPECT_MATRIX_NEAR(I, stan::math::mdivide_left(Ad, Ad), 1e-15);
 }
 
 TEST(MathMatrixPrim, mdivide_right_val2) {