Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closures #2384

Closed
wants to merge 35 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
35 commits
Select commit Hold shift + click to select a range
01f0a03
add closure support and a few tests
nhuurre Feb 19, 2021
7382327
some docs
nhuurre Feb 21, 2021
ee21600
fix headers
nhuurre Feb 21, 2021
a595b43
add from_lambda() for suffix functions, and fix integrate_ode adapter
nhuurre Feb 23, 2021
29c165f
add some minimal docs
nhuurre Mar 6, 2021
5dce92a
Merge branch 'develop' into feature/closures-v2
nhuurre Mar 20, 2021
bbabc92
fix reduce_sum off-by-one index
nhuurre Mar 20, 2021
a4032be
Merge branch 'feature/variadic-integrate-1d' into feature/closures-v2
nhuurre Mar 20, 2021
76a8991
closure support for integrate_1d
nhuurre Mar 20, 2021
a55ddb3
Merge branch 'develop' into feature/closures-v2
nhuurre Apr 20, 2021
990c070
move ode_closure_adapter out of internal namespace
nhuurre Jun 15, 2021
cbf48fa
refactor tests
nhuurre Jun 15, 2021
977f54d
Merge branch 'develop' into feature/closures-v2
nhuurre Jun 15, 2021
c043511
[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…
stan-buildbot Jun 15, 2021
e72e6f4
fix cpplint
nhuurre Jun 15, 2021
e120064
fix test
nhuurre Jun 15, 2021
1245fa6
fix promote_scalar_type
nhuurre Jun 15, 2021
9c25817
fix
nhuurre Jun 15, 2021
a749f61
Merge branch 'develop' into feature/closures-v2
nhuurre Jul 14, 2021
6990768
generalize base_closure
nhuurre Jul 19, 2021
2e3180a
non-empty suffix closures
nhuurre Jul 27, 2021
610af2d
cleanup
nhuurre Jul 27, 2021
880e270
remove empty statement
nhuurre Jul 27, 2021
4f1f6eb
fix includes
nhuurre Jul 28, 2021
f883e42
Merge branch 'develop' into feature/closures-v2
nhuurre Jul 28, 2021
75f6d30
missing include
nhuurre Jul 28, 2021
0a609db
fixup names and references logic for closure classes
SteveBronder Aug 10, 2021
e0f6145
remove fn_return_type
SteveBronder Aug 10, 2021
9eca190
Merge commit '02dc560b022f328f917af80a1b9b7f1feb249ee4' into HEAD
yashikno Aug 12, 2021
9436a18
[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…
stan-buildbot Aug 12, 2021
2b2bee2
fix headers
SteveBronder Aug 13, 2021
ecec96a
Merge branch 'develop' into feature/closures-v2
nhuurre Sep 7, 2021
03ab504
some docs, a test
nhuurre Sep 7, 2021
9fc569f
Merge branch 'develop' into feature/closures-v2
nhuurre Sep 23, 2021
9fb3740
Merge branch 'develop' into feature/closures-v2
nhuurre Oct 17, 2021
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions stan/math/prim/err/elementwise_check.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,13 @@ inline void elementwise_check(const F& is_good, const char* function,
}();
}
}
template <typename F, typename T, typename... Indexings,
require_stan_closure_t<T>* = nullptr>
inline void elementwise_check(const F& is_good, const char* function,
const char* name, const T& x, const char* must_be,
const Indexings&... indexings) {
// XXX skip closures
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch. We'll need to implement these.

Copy link
Member

@bbbales2 bbbales2 Mar 4, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, thinking more, I don't think there's a big advantage to implementing this.

We could expose the variables captured by a closure to checks, but the Math checks wouldn't know in what order its getting them, and then depending on which function was accepting closures it would need to decide which checks to do on which inputs.

I think instead in the ODE solves we check only the arguments passed in explicitly (which this is effectively doing) or we get rid of the infinity checks on the inputs to the ODE solves. I'll make an issue and see if getting rid of the checks altogether is an option. (Edit: Issue #2406)

}
/**
* Check that the predicate holds for all elements of the value of `x`. This
* overload works on Eigen types that support linear indexing.
Expand Down
19 changes: 19 additions & 0 deletions stan/math/prim/fun/value_of.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@

#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
#include <stan/math/prim/fun/eval.hpp>
#include <stan/math/prim/functor/apply.hpp>
#include <cstddef>
#include <vector>

Expand Down Expand Up @@ -77,6 +79,23 @@ inline auto value_of(EigMat&& M) {
std::forward<EigMat>(M));
}

/**
* Closures that capture non-arithmetic types have value_of__() method.
*
* @tparam F Input element type
* @param[in] f Input closure
* @return closure
**/
template <typename F, require_stan_closure_t<F>* = nullptr,
require_not_st_arithmetic<F>* = nullptr>
inline auto value_of(const F& f) {
return apply(
[&f](const auto&... s) {
return typename F::partials_closure_t_(f.f_, eval(value_of(s))...);
},
f.captures_);
}

} // namespace math
} // namespace stan

Expand Down
1 change: 1 addition & 0 deletions stan/math/prim/functor.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#include <stan/math/prim/functor/apply_scalar_binary.hpp>
#include <stan/math/prim/functor/apply_vector_unary.hpp>
#include <stan/math/prim/functor/coupled_ode_system.hpp>
#include <stan/math/prim/functor/closure_adapter.hpp>
#include <stan/math/prim/functor/finite_diff_gradient.hpp>
#include <stan/math/prim/functor/finite_diff_gradient_auto.hpp>
#include <stan/math/prim/functor/for_each.hpp>
Expand Down
220 changes: 220 additions & 0 deletions stan/math/prim/functor/closure_adapter.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,220 @@
#ifndef STAN_MATH_PRIM_FUNCTOR_CLOSURE_ADAPTER_HPP
#define STAN_MATH_PRIM_FUNCTOR_CLOSURE_ADAPTER_HPP

#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/functor/apply.hpp>
#include <ostream>

namespace stan {
namespace math {
namespace internal {

/**
* A closure that wraps a C++ lambda and captures values.
*
* @tparam Ref if true values are captured by reference
* @tparam F the lambda functor type
* @tparam Ts types of the captured values
*/
template <bool Ref, typename F, typename... Ts>
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This all needs docs for template parameters etc.

struct base_closure {
using return_scalar_t_ = return_type_t<Ts...>;
/*The base closure with `Ts` as the non-expression partials of `Ts`*/
using partials_closure_t_
= base_closure<false, F, decltype(eval(value_of(std::declval<Ts>())))...>;
using Base_ = base_closure<false, F, Ts...>;
std::decay_t<F> f_;
std::tuple<closure_return_type_t<Ts, Ref>...> captures_;
template <typename FF, require_same_t<FF, F>* = nullptr, typename... Args>
explicit base_closure(FF&& f, Args&&... args)
: f_(std::forward<FF>(f)), captures_(std::forward<Args>(args)...) {}

template <typename... Args>
auto operator()(std::ostream* msgs, const Args&... args) const {
return apply(
[this, msgs, &args...](const auto&... s) {
return this->f_(s..., args..., msgs);
},
captures_);
}
};

/**
* A closure that takes rng argument.
*
* @tparam Ref if true values are captured by reference
* @tparam F the lambda functor type
* @tparam Ts types of the captured values
*/
template <bool Ref, typename F, typename... Ts>
struct closure_rng {
using return_scalar_t_ = double;
using partials_closure_t_ = closure_rng<false, F, Ts...>;
using Base_ = closure_rng<false, F, Ts...>;
std::decay_t<F> f_;
std::tuple<closure_return_type_t<Ts, Ref>...> captures_;

template <typename FF, require_same_t<FF, F>* = nullptr, typename... Args>
explicit closure_rng(FF&& f, Args&&... args)
: f_(std::forward<FF>(f)), captures_(std::forward<Args>(args)...) {}

template <typename Rng, typename... Args>
auto operator()(Rng& rng, std::ostream* msgs, const Args&... args) const {
return apply(
[this, &rng, msgs, &args...](const auto&... s) {
return this->f_(s..., args..., rng, msgs);
},
captures_);
}
};

/**
* A closure that may compute an unnormalized propability density.
*
* @tparam Propto if true the function is unnormalized
* @tparam Ref if true values are captured by reference
* @tparam F the lambda functor type
* @tparam Ts types of the captured values
*/
template <bool Propto, bool Ref, typename F, typename... Ts>
struct closure_lpdf {
using return_scalar_t_ = return_type_t<Ts...>;
using partials_closure_t_ = closure_lpdf<Propto, false, F, Ts...>;
using Base_ = closure_lpdf<Propto, false, F, Ts...>;
std::decay_t<F> f_;
std::tuple<closure_return_type_t<Ts, Ref>...> captures_;

template <typename FF, require_same_t<FF, F>* = nullptr, typename... Args>
explicit closure_lpdf(FF&& f, Args&&... args)
: f_(std::forward<FF>(f)), captures_(std::forward<Args>(args)...) {}

template <bool propto>
auto with_propto() {
Comment on lines +91 to +92
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We use camelcase for template parameters. How is this propto different from Propto in the template parameters of the class?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this for like lupdf or something?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it's an lupdf hack. This normalization nonsense is rather confusing so here's an example:

parameters {
  real y[100];
}
model {
  function
  real higher_lpdf(real[] x, real(real[], int, int) f_lpdf) {
    real lp = 0;
    lp += reduce_sum(f_lpdf, x, 1);  // <-- A
    lp += reduce_sum(f_lupdf, x, 1); // <-- B
    return lp;
  }
  function
  real partial_lpdf(real[] x, int s, int e) {
    return std_normal_lupdf(x|);
  }
  target += higher_lpdf( y| partial_lpdf);  // <-- 1
  target += higher_lupdf(y| partial_lupdf); // <-- 2
}

Using lpdf instead of lupdf anywhere makes the return value normalized.
In the example reduce_sum is called four times, at (1A), (1B), (2A), (2B), and only (2B) is unnormalized.

The above compiles to C++ that looks something like

auto higher_lpdf = from_lambda([&](auto f_lpdf) {
  var lp = 0;
  lp += reduce_sum(f_lpdf.with_propto<false>(), x, 1); // <-- A
  lp += reduce_sum(f_lpdf.with_propto<true>(), x, 1);  // <-- B
  return lp;
});
auto partial_lpdf = from_lambda([]<bool propto>(auto x, int s, int e) {
  return std_normal_lpdf<propto>(x);
});
lp_accum__.add(higher_lpdf(y, partial_lpdf.with_propto<false>()); // <-- 1
lp_accum__.add(higher_lpdf(y, partial_lpdf.with_propto<true>());  // <-- 2

Every time the closure object is passed to a higher-order function Propto records if it's in its lpdf or lupdf form.

return apply(
[this](const auto&... args) {
return closure_lpdf < Propto && propto, true, F,
Ts... > (this->f_, args...);
},
captures_);
}

template <bool propto = false, typename... Args>
auto operator()(std::ostream* msgs, const Args&... args) const {
return apply(
[this, msgs, &args...](const auto&... s) {
return this->f_.template operator()<propto>(s..., args..., msgs);
},
captures_);
}
};

/**
* A closure that accesses logprob accumulator.
*
* @tparam Propto if true the logprob is unnormalized
* @tparam Ref if true values are captured by reference
* @tparam F the lambda functor type
* @tparam Ts types of the captured values
*/
template <bool Propto, bool Ref, typename F, typename... Ts>
struct closure_lp {
using return_scalar_t_ = return_type_t<Ts...>;
using partials_closure_t_ = closure_lp<Propto, true, F, Ts...>;
using Base_ = closure_lp<Propto, true, F, Ts...>;
std::decay_t<F> f_;
std::tuple<closure_return_type_t<Ts, Ref>...> captures_;

template <typename FF, require_same_t<FF, F>* = nullptr, typename... Args>
explicit closure_lp(FF&& f, Args&&... args)
: f_(std::forward<FF>(f)), captures_(std::forward<Args>(args)...) {}

template <bool propto = false, typename T_lp, typename T_lp_accum,
typename... Args>
auto operator()(T_lp& lp, T_lp_accum& lp_accum, std::ostream* msgs,
const Args&... args) const {
return apply(
[this, &lp, &lp_accum, msgs, &args...](const auto&... s) {
return this->f_.template operator()<propto>(s..., args..., lp,
lp_accum, msgs);
},
captures_);
}
};

} // namespace internal

/**
* Higher-order functor suitable for calling a closure inside variadic ODE
* solvers.
*/
struct ode_closure_adapter {
template <typename F, typename T0, typename T1, typename... Args>
auto operator()(const T0& t, const T1& y, std::ostream* msgs, F&& f,
Args&&... args) const {
return std::forward<F>(f)(msgs, t, y, std::forward<Args>(args)...);
}
};

struct integrate_ode_closure_adapter {
template <typename F, typename T0, typename T1, typename... Args>
auto operator()(const T0& t, const T1& y, std::ostream* msgs, F&& f,
Args&&... args) const {
return to_vector(std::forward<F>(f)(msgs, t, to_array_1d(y),
std::forward<Args>(args)...));
}
};

/**
* Create a closure from a C++ lambda and captures.
*/
template <typename F, typename... Args>
auto from_lambda(F&& f, Args&&... args) {
return internal::base_closure<true, F, Args...>(std::forward<F>(f),
std::forward<Args>(args)...);
}

/**
* Create a closure from an rng functor.
*/
template <typename F, typename... Args>
auto rng_from_lambda(F&& f, Args&&... args) {
return internal::closure_rng<true, F, Args...>(std::forward<F>(f),
std::forward<Args>(args)...);
}

/**
* Create a closure from an lpdf functor.
*/
template <bool propto, typename F, typename... Args>
auto lpdf_from_lambda(F&& f, Args&&... args) {
return internal::closure_lpdf<propto, true, F, Args...>(
std::forward<F>(f), std::forward<Args>(args)...);
}

/**
* Create a closure from a functor that needs access to logprob accumulator.
*/
template <bool Propto, typename F, typename... Args>
auto lp_from_lambda(F&& f, Args&&... args) {
return internal::closure_lp<Propto, true, F, Args...>(
std::forward<F>(f), std::forward<Args>(args)...);
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So for each kind of function we might have in Stan, we can also have a closure version of that.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Each function kind follows a different calling convention so each kind needs its own adapter closure. These aren't used in math library but stanc3 allows userdefined higher order functions that might need them.


/**
* Higher-order functor that invokes a closure inside a reduce_sum call.
*/
struct reduce_sum_closure_adapter {
template <typename F, typename T, typename... Args>
auto operator()(const std::vector<T>& sub_slice, std::size_t start,
std::size_t end, std::ostream* msgs, F&& f,
Args&&... args) const {
return std::forward<F>(f)(msgs, sub_slice, start + error_index::value,
end + error_index::value,
std::forward<Args>(args)...);
}
};

} // namespace math
} // namespace stan

#endif
14 changes: 13 additions & 1 deletion stan/math/prim/functor/integrate_1d.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -236,7 +236,7 @@ inline double integrate_1d_impl(const F& f, double a, double b,
* @param relative_tolerance tolerance passed to Boost quadrature
* @return numeric integral of function f
*/
template <typename F>
template <typename F, require_not_stan_closure_t<F>* = nullptr>
inline double integrate_1d(const F& f, double a, double b,
const std::vector<double>& theta,
const std::vector<double>& x_r,
Expand All @@ -247,6 +247,18 @@ inline double integrate_1d(const F& f, double a, double b,
msgs, theta, x_r, x_i);
}

template <typename F, require_stan_closure_t<F>* = nullptr,
require_arithmetic_t<return_type_t<F>>* = nullptr>
inline double integrate_1d(const F& f, double a, double b,
const std::vector<double>& theta,
const std::vector<double>& x_r,
const std::vector<int>& x_i, std::ostream* msgs,
const double relative_tolerance
= std::sqrt(EPSILON)) {
return integrate_1d_impl(integrate_1d_closure_adapter(), a, b,
relative_tolerance, msgs, f, theta, x_r, x_i);
}
Comment on lines +250 to +260
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[Q] Think I'm just missing some context, why is f being passed as an argument here? instead of being pushed to integrate_1d_adapter?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

integrate_1d_impl already knows how to handle the autodiff variables passed as arguments but would need new logic for extracting them from the functor.


} // namespace math
} // namespace stan

Expand Down
15 changes: 15 additions & 0 deletions stan/math/prim/functor/integrate_1d_adapter.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,4 +25,19 @@ struct integrate_1d_adapter {
}
};

/**
* Call a closure object from integrate_1d
*/
struct integrate_1d_closure_adapter {
integrate_1d_closure_adapter() {}

template <typename F, typename T_a, typename T_b, typename T_theta>
auto operator()(const T_a& x, const T_b& xc, std::ostream* msgs, const F& f,
const std::vector<T_theta>& theta,
const std::vector<double>& x_r,
const std::vector<int>& x_i) const {
return f(msgs, x, xc, theta, x_r, x_i);
}
};

#endif
42 changes: 37 additions & 5 deletions stan/math/prim/functor/integrate_ode_rk45.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#define STAN_MATH_PRIM_FUNCTOR_INTEGRATE_ODE_RK45_HPP

#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/functor/closure_adapter.hpp>
#include <stan/math/prim/functor/integrate_ode_std_vector_interface_adapter.hpp>
#include <stan/math/prim/functor/ode_rk45.hpp>
#include <ostream>
Expand All @@ -10,6 +11,38 @@
namespace stan {
namespace math {

namespace internal {

template <typename F, typename T_y0, typename T_param, typename T_t0,
typename T_ts, require_not_stan_closure_t<F>* = nullptr>
inline auto integrate_ode_rk45_impl(
const F& f, const std::vector<T_y0>& y0, const T_t0& t0,
const std::vector<T_ts>& ts, const std::vector<T_param>& theta,
const std::vector<double>& x, const std::vector<int>& x_int,
std::ostream* msgs, double relative_tolerance, double absolute_tolerance,
int max_num_steps) {
internal::integrate_ode_std_vector_interface_adapter<F> f_adapted(f);
return ode_rk45_tol_impl("integrate_ode_rk45", f_adapted, to_vector(y0), t0,
ts, relative_tolerance, absolute_tolerance,
max_num_steps, msgs, theta, x, x_int);
}

template <typename F, typename T_y0, typename T_param, typename T_t0,
typename T_ts, require_stan_closure_t<F>* = nullptr>
inline auto integrate_ode_rk45_impl(
const F& f, const std::vector<T_y0>& y0, const T_t0& t0,
const std::vector<T_ts>& ts, const std::vector<T_param>& theta,
const std::vector<double>& x, const std::vector<int>& x_int,
std::ostream* msgs, double relative_tolerance, double absolute_tolerance,
int max_num_steps) {
return ode_rk45_tol_impl("integrate_ode_rk45",
integrate_ode_closure_adapter(), to_vector(y0), t0,
ts, relative_tolerance, absolute_tolerance,
max_num_steps, msgs, f, theta, x, x_int);
}

} // namespace internal

/**
* @deprecated use <code>ode_rk45</code>
*/
Expand All @@ -21,12 +54,11 @@ inline auto integrate_ode_rk45(
const std::vector<double>& x, const std::vector<int>& x_int,
std::ostream* msgs = nullptr, double relative_tolerance = 1e-6,
double absolute_tolerance = 1e-6, int max_num_steps = 1e6) {
internal::integrate_ode_std_vector_interface_adapter<F> f_adapted(f);
auto y = ode_rk45_tol_impl("integrate_ode_rk45", f_adapted, to_vector(y0), t0,
ts, relative_tolerance, absolute_tolerance,
max_num_steps, msgs, theta, x, x_int);
auto y = internal::integrate_ode_rk45_impl(f, y0, t0, ts, theta, x, x_int,
msgs, relative_tolerance,
absolute_tolerance, max_num_steps);

std::vector<std::vector<return_type_t<T_y0, T_param, T_t0, T_ts>>>
std::vector<std::vector<return_type_t<F, T_y0, T_param, T_t0, T_ts>>>
y_converted;
y_converted.reserve(y.size());
for (size_t i = 0; i < y.size(); ++i)
Expand Down
Loading