diff --git a/stan/math/prim/fun/max.hpp b/stan/math/prim/fun/max.hpp index 2f0af8187fc..28e4dad67c6 100644 --- a/stan/math/prim/fun/max.hpp +++ b/stan/math/prim/fun/max.hpp @@ -21,7 +21,8 @@ namespace math { * @param y second argument * @return maximum value of the two arguments */ -template * = nullptr> +template * = nullptr, + require_all_not_var_t* = nullptr> inline auto max(T1 x, T2 y) { return std::max(x, y); } @@ -38,7 +39,8 @@ inline auto max(T1 x, T2 y) { * @throws std::invalid_argument if the vector is size zero and the * scalar type in the container is integer */ -template * = nullptr> +template * = nullptr, + require_not_st_var* = nullptr> inline value_type_t max(T&& m) { if constexpr (std::is_integral>::value) { check_nonzero_size("max", "int vector", m); diff --git a/stan/math/prim/fun/min.hpp b/stan/math/prim/fun/min.hpp index a6e006f4c74..fd724a03249 100644 --- a/stan/math/prim/fun/min.hpp +++ b/stan/math/prim/fun/min.hpp @@ -37,7 +37,8 @@ inline auto min(T1 x, T2 y) { * @throws std::invalid_argument if the vector is size zero and the * scalar type in the container is integer */ -template * = nullptr> +template * = nullptr, + require_not_st_var* = nullptr> inline value_type_t min(T&& m) { if constexpr (std::is_integral>::value) { check_nonzero_size("min", "int vector", m); diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index b9a8bfa3e2c..0901797d37a 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -118,10 +118,12 @@ #include #include #include +#include #include #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/max.hpp b/stan/math/rev/fun/max.hpp new file mode 100644 index 00000000000..5ee54d382d7 --- /dev/null +++ b/stan/math/rev/fun/max.hpp @@ -0,0 +1,129 @@ +#ifndef STAN_MATH_REV_FUN_MAX_HPP +#define STAN_MATH_REV_FUN_MAX_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Returns the maximum value of the two specified scalar arguments, + * where at least one argument is a 'var'. + * + * @tparam Scal1 type of first argument (must be 'var' if Scal2 is not) + * @tparam Scal2 type of second argument (must be 'var' if Scal1 is not) + * @param[in] x first argument + * @param[in] y second argument + * @return maximum value of the two arguments + */ +template * = nullptr> +inline var max(Scal1 x, Scal2 y) { + double x_val = stan::math::value_of(x); + double y_val = stan::math::value_of(y); + if (unlikely(is_nan(x_val))) { + if (unlikely(is_nan(y_val))) { + return make_callback_var(NOT_A_NUMBER, [x, y](auto& vi) mutable { + if constexpr (is_var::value) + x.adj() = NOT_A_NUMBER; + if constexpr (is_var::value) + y.adj() = NOT_A_NUMBER; + }); + } + return y; + } + if (unlikely(is_nan(y_val))) { + return x; + } + double max_val = std::fmax(x_val, y_val); + return make_callback_var(max_val, [x, y, x_val, y_val](auto& vi) mutable { + if (x_val > y_val) { + if constexpr (is_var::value) + x.adj() += vi.adj(); + } else if (y_val > x_val) { + if constexpr (is_var::value) + y.adj() += vi.adj(); + } else { + if constexpr (is_var::value) + x.adj() += vi.adj() * 0.5; + if constexpr (is_var::value) + y.adj() += vi.adj() * 0.5; + } + }); +} +/** + * Return the maximum value in a container of 'var'. + * @tparam T Type of container with 'var' scalar type + * @param[in] x container (Eigen matrix, vector, row vector or std vector) + * @return maximum value in the container, or -infinity if size is zero + */ +template * = nullptr, + require_container_t* = nullptr, + require_not_var_matrix_t* = nullptr> +inline var max(T&& x) { + if (unlikely(x.size() == 0)) { + return NEGATIVE_INFTY; + } + return apply_vector_unary::apply(std::forward(x), [](auto&& v) { + const auto& x_ref = to_ref(v); + arena_t> x_arena(x_ref); + double max_val = x_arena.val().maxCoeff(); + return make_callback_var(max_val, [x_arena, max_val](auto& vi) mutable { + if (unlikely(is_nan(max_val))) { + for (Eigen::Index i = 0; i < x_arena.size(); ++i) { + if (is_nan(x_arena.val().coeff(i))) { + x_arena.adj().coeffRef(i) = NOT_A_NUMBER; + } + } + return; + } + auto is_max = (x_arena.val().array() == max_val); + double count = is_max.template cast().sum(); + double adj_to_propagate = vi.adj() / count; + x_arena.adj().array() + += is_max.template cast() * adj_to_propagate; + }); + }); +} + +/** + * Return the maximum value in a `var_value` containing an Eigen matrix, + * vector, or row vector. + * + * @tparam T Type of `var_value` with an Eigen inner type + * @param[in] x `var_value` matrix + * @return maximum value in the matrix, or -infinity if size is zero + */ +template * = nullptr> +inline var max(const T& x) { + if (unlikely(x.size() == 0)) { + return NEGATIVE_INFTY; + } + double max_val = x.val().maxCoeff(); + return make_callback_var(max_val, [x, max_val](auto& vi) mutable { + if (unlikely(is_nan(max_val))) { + for (Eigen::Index i = 0; i < x.size(); ++i) { + if (is_nan(x.val().coeff(i))) { + x.adj().coeffRef(i) = NOT_A_NUMBER; + } + } + return; + } + auto is_max = (x.val().array() == max_val); + double count = is_max.template cast().sum(); + double adj_to_propagate = vi.adj() / count; + x.adj().array() += is_max.template cast() * adj_to_propagate; + }); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun/min.hpp b/stan/math/rev/fun/min.hpp new file mode 100644 index 00000000000..866c338e69f --- /dev/null +++ b/stan/math/rev/fun/min.hpp @@ -0,0 +1,131 @@ +#ifndef STAN_MATH_REV_FUN_MIN_HPP +#define STAN_MATH_REV_FUN_MIN_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Returns the minimum value of the two specified scalar arguments, + * where at least one argument is a 'var'. + * + * @tparam Scal1 type of first argument (must be 'var' if Scal2 is not) + * @tparam Scal2 type of second argument (must be 'var' if Scal1 is not) + * @param[in] x first argument + * @param[in] y second argument + * @return minimum value of the two arguments + */ +template * = nullptr> +inline var min(Scal1 x, Scal2 y) { + double x_val = stan::math::value_of(x); + double y_val = stan::math::value_of(y); + if (unlikely(is_nan(x_val))) { + if (unlikely(is_nan(y_val))) { + return make_callback_var(NOT_A_NUMBER, [x, y](auto& vi) mutable { + if constexpr (is_var::value) + x.adj() = NOT_A_NUMBER; + if constexpr (is_var::value) + y.adj() = NOT_A_NUMBER; + }); + } + return y; + } + if (unlikely(is_nan(y_val))) { + return x; + } + double min_val = std::fmin(x_val, y_val); + return make_callback_var(min_val, [x, y, x_val, y_val](auto& vi) mutable { + if (x_val < y_val) { + if constexpr (is_var::value) + x.adj() += vi.adj(); + } else if (y_val < x_val) { + if constexpr (is_var::value) + y.adj() += vi.adj(); + } else { + if constexpr (is_var::value) + x.adj() += vi.adj() * 0.5; + if constexpr (is_var::value) + y.adj() += vi.adj() * 0.5; + } + }); +} + +/** + * Return the minimum value in a container of 'var'. + * + * @tparam T Type of container with 'var' scalar type + * @param[in] x container (Eigen matrix, vector, row vector or std vector) + * @return minimum value in the container, or positive infinity if size is zero + */ +template * = nullptr, + require_container_t* = nullptr, + require_not_var_matrix_t* = nullptr> +inline var min(T&& x) { + if (unlikely(x.size() == 0)) { + return INFTY; + } + return apply_vector_unary::apply(std::forward(x), [](auto&& v) { + const auto& x_ref = to_ref(v); + arena_t> x_arena(x_ref); + double min_val = x_arena.val().minCoeff(); + return make_callback_var(min_val, [x_arena, min_val](auto& vi) mutable { + if (unlikely(is_nan(min_val))) { + for (Eigen::Index i = 0; i < x_arena.size(); ++i) { + if (is_nan(x_arena.val().coeff(i))) { + x_arena.adj().coeffRef(i) = NOT_A_NUMBER; + } + } + return; + } + auto is_min = (x_arena.val().array() == min_val); + double count = is_min.template cast().sum(); + double adj_to_propagate = vi.adj() / count; + x_arena.adj().array() + += is_min.template cast() * adj_to_propagate; + }); + }); +} + +/** + * Return the minimum value in a `var_value` containing an Eigen matrix, + * vector, or row vector. + * + * @tparam T Type of `var_value` with an Eigen inner type + * @param[in] x `var_value` matrix + * @return minimum value in the matrix, or infinity if size is zero + */ +template * = nullptr> +inline var min(const T& x) { + if (unlikely(x.size() == 0)) { + return INFTY; + } + double min_val = x.val().minCoeff(); + return make_callback_var(min_val, [x, min_val](auto& vi) mutable { + if (unlikely(is_nan(min_val))) { + for (Eigen::Index i = 0; i < x.size(); ++i) { + if (is_nan(x.val().coeff(i))) { + x.adj().coeffRef(i) = NOT_A_NUMBER; + } + } + return; + } + auto is_min = (x.val().array() == min_val); + double count = is_min.template cast().sum(); + double adj_to_propagate = vi.adj() / count; + x.adj().array() += is_min.template cast() * adj_to_propagate; + }); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/mix/fun/max_test.cpp b/test/unit/math/mix/fun/max_test.cpp index fe60483fa67..c0ecb32532a 100644 --- a/test/unit/math/mix/fun/max_test.cpp +++ b/test/unit/math/mix/fun/max_test.cpp @@ -13,25 +13,48 @@ inline void expect_max(const T& m) { stan::test::expect_ad(f, v); stan::test::expect_ad(f, rv); stan::test::expect_ad(f, m); + stan::test::expect_ad_matvar(f, v); + stan::test::expect_ad_matvar(f, rv); + stan::test::expect_ad_matvar(f, m); } TEST(MathMixMatFun, max) { Eigen::MatrixXd a(0, 0); expect_max(a); - Eigen::MatrixXd b(1, 1); b << 12; expect_max(b); - Eigen::MatrixXd c(3, 1); c << 100, 0, -3; expect_max(c); - Eigen::MatrixXd d(1, 3); d << 100, 0, -3; expect_max(d); - Eigen::MatrixXd e(3, 2); e << -100, 0, 1, 20, -40, 2; expect_max(e); + Eigen::MatrixXd ties(2, 2); + ties << 10.5, 1.0, 10.5, 10.5; + expect_max(ties); + double inf = std::numeric_limits::infinity(); + Eigen::VectorXd o(3); + o << -inf, 5.0, inf; + expect_max(o); + double nan = std::numeric_limits::quiet_NaN(); + Eigen::VectorXd n1(3); + n1 << 1.0, nan, 2.0; + expect_max(n1); + Eigen::MatrixXd n2(2, 2); + n2 << nan, nan, nan, nan; + expect_max(n2); +} +TEST(MathMixMatFun, max_binary) { + auto f = [](const auto& x, const auto& y) { return stan::math::max(x, y); }; + stan::test::expect_ad(f, 1.0, 2.0); + stan::test::expect_ad(f, 2.0, 1.0); + stan::test::expect_ad(f, 3.0, 3.0); + double nan = std::numeric_limits::quiet_NaN(); + stan::test::expect_ad(f, nan, 1.0); + stan::test::expect_ad(f, 1.0, nan); + stan::test::expect_ad(f, nan, nan); } diff --git a/test/unit/math/rev/fun/max_test.cpp b/test/unit/math/rev/fun/max_test.cpp new file mode 100644 index 00000000000..fda5d6c3782 --- /dev/null +++ b/test/unit/math/rev/fun/max_test.cpp @@ -0,0 +1,61 @@ +#include +#include +#include + +TEST(MathRev, max_scalar_permutations) { + auto f = [](const auto& x, const auto& y) { return stan::math::max(x, y); }; + stan::test::expect_ad(f, 1.0, 2.0); + stan::test::expect_ad(f, 2.0, 1.0); + stan::test::expect_ad(f, -5.0, 0.0); +} +TEST(MathRev, max_scalar_tie) { + auto f = [](const auto& x, const auto& y) { return stan::math::max(x, y); }; + stan::test::expect_ad(f, 3.0, 3.0); +} +TEST(MathRev, max_scalar_boundary) { + auto f = [](const auto& x, const auto& y) { return stan::math::max(x, y); }; + double nan = std::numeric_limits::quiet_NaN(); + double inf = std::numeric_limits::infinity(); + stan::test::expect_ad(f, nan, 1.0); + stan::test::expect_ad(f, 1.0, nan); + stan::test::expect_ad(f, inf, -inf); +} +TEST(MathRev, max_container_types) { + auto f = [](const auto& x) { return stan::math::max(x); }; + Eigen::VectorXd v(3); + v << 1, 2, 3; + Eigen::RowVectorXd rv(3); + rv << 1, 2, 3; + std::vector sv{1, 2, 3}; + Eigen::MatrixXd m(2, 2); + m << 1, 5, 2, 3; + stan::test::expect_ad(f, v); + stan::test::expect_ad(f, rv); + stan::test::expect_ad(f, sv); + stan::test::expect_ad(f, m); +} +TEST(MathRev, max_container_ties) { + auto f = [](const auto& x) { return stan::math::max(x); }; + Eigen::VectorXd v(4); + v << 1.5, 1.5, 1.5, 1.5; + stan::test::ad_tolerances tols; + tols.gradient_grad_ = 1e-3; + stan::test::expect_ad(tols, f, v); +} +TEST(MathRev, max_container_edges) { + Eigen::VectorXd v_single(1); + v_single << 42.0; + EXPECT_FLOAT_EQ(42.0, stan::math::value_of(stan::math::max(v_single))); + Eigen::VectorXd v_empty(0); + EXPECT_FLOAT_EQ(stan::math::NEGATIVE_INFTY, stan::math::max(v_empty)); +} +TEST(MathRev, max_container_nan) { + auto f = [](const auto& x) { return stan::math::max(x); }; + double nan = std::numeric_limits::quiet_NaN(); + Eigen::VectorXd v1(3); + v1 << 1.0, nan, 2.0; + stan::test::expect_ad(f, v1); + Eigen::VectorXd v2(2); + v2 << nan, nan; + stan::test::expect_ad(f, v2); +} diff --git a/test/unit/math/rev/fun/min_test.cpp b/test/unit/math/rev/fun/min_test.cpp new file mode 100644 index 00000000000..b79f1b9f960 --- /dev/null +++ b/test/unit/math/rev/fun/min_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +TEST(MathRev, min_scalar_permutations) { + auto f = [](const auto& x, const auto& y) { return stan::math::min(x, y); }; + stan::test::expect_ad(f, 1.0, 2.0); + stan::test::expect_ad(f, 2.0, 1.0); +} +TEST(MathRev, min_scalar_tie) { + auto f = [](const auto& x, const auto& y) { return stan::math::min(x, y); }; + stan::test::expect_ad(f, 0.0, 0.0); +} +TEST(MathRev, min_scalar_boundary) { + auto f = [](const auto& x, const auto& y) { return stan::math::min(x, y); }; + double nan = std::numeric_limits::quiet_NaN(); + stan::test::expect_ad(f, nan, 1.0); + stan::test::expect_ad(f, 1.0, nan); +} +TEST(MathRev, min_container_types) { + auto f = [](const auto& x) { return stan::math::min(x); }; + Eigen::VectorXd v(3); + v << 3, 2, 1; + Eigen::RowVectorXd rv(3); + rv << 3, 2, 1; + std::vector sv{3, 2, 1}; + Eigen::MatrixXd m(2, 2); + m << 1, 0.5, 2, 3; + stan::test::expect_ad(f, v); + stan::test::expect_ad(f, rv); + stan::test::expect_ad(f, sv); + stan::test::expect_ad(f, m); +} +TEST(MathRev, min_container_ties) { + auto f = [](const auto& x) { return stan::math::min(x); }; + Eigen::VectorXd v(4); + v << -1.1, -1.1, -1.1, -1.1; + stan::test::ad_tolerances tols; + tols.gradient_grad_ = 1e-3; + stan::test::expect_ad(tols, f, v); +} +TEST(MathRev, min_container_edges) { + Eigen::VectorXd v_single(1); + v_single << -42.0; + EXPECT_FLOAT_EQ(-42.0, stan::math::value_of(stan::math::min(v_single))); + Eigen::VectorXd v_empty(0); + EXPECT_FLOAT_EQ(stan::math::INFTY, stan::math::min(v_empty)); +} +TEST(MathRev, min_container_nan) { + auto f = [](const auto& x) { return stan::math::min(x); }; + double nan = std::numeric_limits::quiet_NaN(); + Eigen::VectorXd v1(3); + v1 << 10.0, nan, -5.0; + stan::test::expect_ad(f, v1); + Eigen::VectorXd v2(2); + v2 << nan, nan; + stan::test::expect_ad(f, v2); +}