[WIP] Adds boost root finders with reverse mode specializations

stan/math/fwd/fun.hpp

@@ -38,6 +38,7 @@
 #include <stan/math/fwd/fun/fmax.hpp>
 #include <stan/math/fwd/fun/fmin.hpp>
 #include <stan/math/fwd/fun/fmod.hpp>
+#include <stan/math/fwd/fun/frexp.hpp>
 #include <stan/math/fwd/fun/gamma_p.hpp>
 #include <stan/math/fwd/fun/gamma_q.hpp>
 #include <stan/math/fwd/fun/grad_inc_beta.hpp>
@@ -103,6 +104,7 @@
 #include <stan/math/fwd/fun/read_fvar.hpp>
 #include <stan/math/fwd/fun/rising_factorial.hpp>
 #include <stan/math/fwd/fun/round.hpp>
+#include <stan/math/fwd/fun/sign.hpp>
 #include <stan/math/fwd/fun/sin.hpp>
 #include <stan/math/fwd/fun/sinh.hpp>
 #include <stan/math/fwd/fun/softmax.hpp>

stan/math/fwd/fun/abs.hpp

@@ -6,6 +6,7 @@
 #include <stan/math/prim/fun/abs.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
 #include <complex>
 
 namespace stan {

stan/math/fwd/fun/frexp.hpp

@@ -0,0 +1,17 @@
+#ifndef STAN_MATH_FWD_FUN_FREXP_HPP
+#define STAN_MATH_FWD_FUN_FREXP_HPP
+
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
+
+namespace stan {
+namespace math {
+
+template <typename T>
+inline auto frexp(const fvar<T>& x, int* exponent) noexcept {
+  return std::frexp(value_of_rec(x), exponent);
+}
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/fwd/fun/sign.hpp

@@ -0,0 +1,19 @@
+#ifndef STAN_MATH_FWD_FUN_SIGN_HPP
+#define STAN_MATH_FWD_FUN_SIGN_HPP
+
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
+
+namespace stan {
+namespace math {
+
+template <typename T>
+inline auto sign(const fvar<T>& x) {
+  double x_val = value_of_rec(x);
+  return (0. < x_val) - (x_val < 0.);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/prim/functor.hpp

@@ -26,5 +26,5 @@
 #include <stan/math/prim/functor/operands_and_partials.hpp>
 #include <stan/math/prim/functor/reduce_sum.hpp>
 #include <stan/math/prim/functor/reduce_sum_static.hpp>
-
+#include <stan/math/prim/functor/root_finder.hpp>
 #endif

stan/math/prim/functor/root_finder.hpp

@@ -0,0 +1,181 @@
+#ifndef STAN_MATH_PRIM_FUNCTOR_ROOT_FINDER_HPP
+#define STAN_MATH_PRIM_FUNCTOR_ROOT_FINDER_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err/check_bounded.hpp>
+#include <stan/math/prim/err/check_positive.hpp>
+#include <stan/math/prim/functor/apply.hpp>
+#include <boost/math/tools/roots.hpp>
+#include <tuple>
+#include <utility>
+
+namespace stan {
+namespace math {
+namespace internal {
+template <typename Tuple, typename... Args>
+inline auto func_with_derivs(Tuple&& f_tuple, Args&&... args) {
+  return stan::math::apply(
+      [&args...](auto&&... funcs) {
+        return [&args..., &funcs...](auto&& g) {
+          return std::make_tuple(funcs(g, args...)...);
+        };
+      },
+      f_tuple);
+}
+}  // namespace internal
+
+/**
+ * Solve for root using Boost's Halley method
+ * @tparam FTuple A tuple holding functors whose signatures all match
+ *  `(GuessScalar g, Types&&... Args)`.
+ * @tparam GuessScalar Scalar type
+ * @tparam MinScalar Scalar type
+ * @tparam MaxScalar Scalar type
+ * @tparam Types Arg types to pass to functors in `f_tuple`
+ * @param f_tuple A tuple of functors to calculate the value and any derivates
+ * needed.
+ * @param guess An initial guess at the root value
+ * @param min The minimum possible value for the result, this is used as an
+ * initial lower bracket
+ * @param max The maximum possible value for the result, this is used as an
+ * initial upper bracket
+ * @param digits The desired number of binary digits precision
+ * @param max_iter An optional maximum number of iterations to perform. On exit,
+ * this is updated to the actual number of iterations performed
+ * @param args Parameter pack of arguments to pass the the functors in `f_tuple`
+ */
+template <typename SolverFun, typename FTuple, typename GuessScalar,
+          typename MinScalar, typename MaxScalar, typename... Types,
+          require_all_not_st_var<GuessScalar, MinScalar, MaxScalar,
+                                 Types...>* = nullptr>
+auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple,
+                     const GuessScalar guess, const MinScalar min,
+                     const MaxScalar max, const int digits,
+                     std::uintmax_t& max_iter, Types&&... args) {
+  check_bounded("root_finder", "initial guess", guess, min, max);
+  check_positive("root_finder", "digits", digits);
+  check_positive("root_finder", "max_iter", max_iter);
+  using ret_t = return_type_t<GuessScalar, MinScalar, MaxScalar, Types...>;
+  ret_t ret = 0;
+  auto f_plus_div = internal::func_with_derivs(f_tuple, args...);
+  try {
+    ret = f_solver(f_plus_div, ret_t(guess), ret_t(min), ret_t(max), digits,
+                   max_iter);
+  } catch (const std::exception& e) {
+    std::cout << "err: \n" << e.what() << "\n";
+    throw e;
+  }
+  return ret;
+}
+
+template <typename FTuple, typename GuessScalar, typename MinScalar,
+          typename MaxScalar, typename... Types>
+auto root_finder_halley_tol(FTuple&& f_tuple, const GuessScalar guess,
+                            const MinScalar min, const MaxScalar max,
+                            const int digits, std::uintmax_t& max_iter,
+                            Types&&... args) {
+  return root_finder_tol(
+      [](auto&&... args) {
+        return boost::math::tools::halley_iterate(args...);
+      },
+      std::forward<FTuple>(f_tuple), guess, min, max, digits, max_iter,
+      std::forward<Types>(args)...);
+}
+
+template <typename FTuple, typename GuessScalar, typename MinScalar,
+          typename MaxScalar, typename... Types>
+auto root_finder_newton_raphson_tol(FTuple&& f_tuple, const GuessScalar guess,
+                                    const MinScalar min, const MaxScalar max,
+                                    const int digits, std::uintmax_t& max_iter,
+                                    Types&&... args) {
+  return root_finder_tol(
+      [](auto&&... args) {
+        return boost::math::tools::newton_raphson_iterate(args...);
+      },
+      std::forward<FTuple>(f_tuple), guess, min, max, digits, max_iter,
+      std::forward<Types>(args)...);
+}
+
+template <typename FTuple, typename GuessScalar, typename MinScalar,
+          typename MaxScalar, typename... Types>
+auto root_finder_schroder_tol(FTuple&& f_tuple, const GuessScalar guess,
+                              const MinScalar min, const MaxScalar max,
+                              const int digits, std::uintmax_t& max_iter,
+                              Types&&... args) {
+  return root_finder_tol(
+      [](auto&&... args) {
+        return boost::math::tools::schroder_iterate(args...);
+      },
+      std::forward<FTuple>(f_tuple), guess, min, max, digits, max_iter,
+      std::forward<Types>(args)...);
+}
+
+/**
+ * Solve for root using Boost Halley method with default values for the
+ * tolerances
+ * @tparam FTuple A tuple holding functors whose signatures all match
+ *  `(GuessScalar g, Types&&... Args)`.
+ * @tparam GuessScalar Scalar type
+ * @tparam MinScalar Scalar type
+ * @tparam MaxScalar Scalar type
+ * @tparam Types Arg types to pass to functors in `f_tuple`
+ * @param f_tuple A tuple of functors to calculate the value and any derivates
+ * needed.
+ * @param guess An initial guess at the root value
+ * @param min The minimum possible value for the result, this is used as an
+ * initial lower bracket
+ * @param max The maximum possible value for the result, this is used as an
+ * initial upper bracket
+ * @param args Parameter pack of arguments to pass the the functors in `f_tuple`
+ */
+template <typename SolverFun, typename FTuple, typename GuessScalar,
+          typename MinScalar, typename MaxScalar, typename... Types>
+auto root_finder(SolverFun&& f_solver, FTuple&& f_tuple,
+                 const GuessScalar guess, const MinScalar min,
+                 const MaxScalar max, Types&&... args) {
+  constexpr int digits = 16;
+  std::uintmax_t max_iter = std::numeric_limits<std::uintmax_t>::max();
+  return root_finder_tol(std::forward<SolverFun>(f_solver),
+                         std::forward<FTuple>(f_tuple), guess, min, max, digits,
+                         max_iter, std::forward<Types>(args)...);
+}
+
+template <typename FTuple, typename GuessScalar, typename MinScalar,
+          typename MaxScalar, typename... Types>
+auto root_finder_hailey(FTuple&& f_tuple, const GuessScalar guess,
+                        const MinScalar min, const MaxScalar max,
+                        Types&&... args) {
+  constexpr int digits = 16;
+  std::uintmax_t max_iter = std::numeric_limits<std::uintmax_t>::max();
+  return root_finder_halley_tol(std::forward<FTuple>(f_tuple), guess, min, max,
+                                digits, max_iter, std::forward<Types>(args)...);
+}
+
+template <typename FTuple, typename GuessScalar, typename MinScalar,
+          typename MaxScalar, typename... Types>
+auto root_finder_newton_raphson(FTuple&& f_tuple, const GuessScalar guess,
+                                const MinScalar min, const MaxScalar max,
+                                Types&&... args) {
+  constexpr int digits = 16;
+  std::uintmax_t max_iter = std::numeric_limits<std::uintmax_t>::max();
+  return root_finder_newton_raphson_tol(std::forward<FTuple>(f_tuple), guess,
+                                        min, max, digits, max_iter,
+                                        std::forward<Types>(args)...);
+}
+
+template <typename FTuple, typename GuessScalar, typename MinScalar,
+          typename MaxScalar, typename... Types>
+auto root_finder_schroder(FTuple&& f_tuple, const GuessScalar guess,
+                          const MinScalar min, const MaxScalar max,
+                          Types&&... args) {
+  constexpr int digits = 16;
+  std::uintmax_t max_iter = std::numeric_limits<std::uintmax_t>::max();
+  return root_finder_schroder_tol(std::forward<FTuple>(f_tuple), guess, min,
+                                  max, digits, max_iter,
+                                  std::forward<Types>(args)...);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/rev/fun.hpp

@@ -67,6 +67,7 @@
 #include <stan/math/rev/fun/fmax.hpp>
 #include <stan/math/rev/fun/fmin.hpp>
 #include <stan/math/rev/fun/fmod.hpp>
+#include <stan/math/rev/fun/frexp.hpp>
 #include <stan/math/rev/fun/from_var_value.hpp>
 #include <stan/math/rev/fun/gamma_p.hpp>
 #include <stan/math/rev/fun/gamma_q.hpp>
@@ -161,6 +162,7 @@
 #include <stan/math/rev/fun/singular_values.hpp>
 #include <stan/math/rev/fun/svd_U.hpp>
 #include <stan/math/rev/fun/svd_V.hpp>
+#include <stan/math/rev/fun/sign.hpp>
 #include <stan/math/rev/fun/sinh.hpp>
 #include <stan/math/rev/fun/softmax.hpp>
 #include <stan/math/rev/fun/sqrt.hpp>

stan/math/rev/fun/frexp.hpp

@@ -0,0 +1,14 @@
+#ifndef STAN_MATH_REV_FUN_FREXP_HPP
+#define STAN_MATH_REV_FUN_FREXP_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+
+namespace stan {
+namespace math {
+inline auto frexp(stan::math::var x, int* exponent) noexcept {
+  return std::frexp(x.val(), exponent);
+}
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/rev/fun/sign.hpp

@@ -0,0 +1,12 @@
+#ifndef STAN_MATH_REV_FUN_SIGN_HPP
+#define STAN_MATH_REV_FUN_SIGN_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+
+namespace stan {
+namespace math {
+inline int sign(stan::math::var z) { return (z == 0) ? 0 : z < 0 ? -1 : 1; }
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/rev/functor.hpp

@@ -27,6 +27,7 @@
 #include <stan/math/rev/functor/map_rect_reduce.hpp>
 #include <stan/math/rev/functor/operands_and_partials.hpp>
 #include <stan/math/rev/functor/reduce_sum.hpp>
+#include <stan/math/rev/functor/root_finder.hpp>
 #include <stan/math/rev/functor/finite_diff_hessian_auto.hpp>
 
 #endif