17 template <
typename T,
typename Enabled =
void>
18 struct TrapzHelper {};
24 struct TrapzHelper<T, typename std::enable_if<std::is_floating_point<T>::value>::type> {
25 static NdArray::NdArray<double>
trapz(
const NdArray::NdArray<double>& grid,
const std::vector<T>& knots,
size_t axis) {
26 if (knots.
size() == 1) {
37 struct TrapzHelper<T, typename std::enable_if<!std::is_floating_point<T>::value>::type> {
38 static NdArray::NdArray<double>
trapz(
const NdArray::NdArray<double>& grid,
const std::vector<T>&,
size_t axis) {
49 template <
typename TKnot>
50 class NdSampler<TKnot> {
55 : m_inv_cumulative(std::
move(knots), std::vector<double>(grid.
begin(), grid.
end())) {
56 if (grid.shape().size() != 1) {
57 throw Elements::Exception() <<
"Grid with " << grid.shape().size() <<
" axes passed to a 1D sampler";
64 template <
typename Generator,
typename... OKnots>
66 constexpr
std::size_t this_n =
sizeof...(OKnots) - 1;
67 static_assert(
sizeof...(OKnots) >= 1,
"The output tuple must have at least one element");
69 if (output.size() != ndraws) {
70 throw Elements::Exception() <<
"Output area does not match the required shape: expected at least " << ndraws <<
", got "
77 for (
auto& row : output) {
78 auto p = uniform(rng);
79 std::get<this_n>(row) = m_inv_cumulative(p);
83 template <
typename Generator>
86 draw(ndraws, rng, samples);
91 const InverseCumulative<TKnot> m_inv_cumulative;
97 template <std::
size_t Start,
typename Seq>
98 struct _CallHelper {};
116 template <
typename Func,
typename TKnot0,
typename... TKnotN>
118 return func(x0, std::get<Start + 1 + Is>(xs)...);
125 template <
typename TKnot0,
typename... TKnotN>
126 class NdSampler<TKnot0, TKnotN...> {
136 if (values.shape().size() != N) {
137 throw Elements::Exception() <<
"Grid with " << values.shape().size() <<
" axes passed to a " << N <<
"D sampler";
144 m_knots0 =
std::move(std::get<0>(knots));
162 template <
typename Generator>
165 draw(ndraws, rng, output);
176 template <
typename Generator,
typename... OutputKnots>
178 constexpr
std::size_t this_n =
sizeof...(OutputKnots) -
sizeof...(TKnotN) - 1;
180 if (output.size() != ndraws) {
181 throw Elements::Exception() <<
"Output area does not match the required shape: expected at least " << ndraws <<
", got "
189 m_subsampler->draw(ndraws, rng, output);
192 for (
std::size_t draw_i = 0; draw_i < ndraws; ++draw_i) {
194 auto& subsample = output[draw_i];
196 for (
std::size_t i = 0; i < m_knots0.size(); ++i) {
197 pdf[i] = _CallHelper<this_n,
_make_index_sequence<
sizeof...(TKnotN)>>::call(*m_interpolation, m_knots0[i], subsample);
200 auto p = uniform(rng);
201 InverseCumulative<TKnot0> inv_cumulative(m_knots0,
std::move(pdf));
202 std::get<this_n>(subsample) = inv_cumulative(p);
_integer_sequence< std::size_t, Idx...> _index_sequence
static std::tuple< Tn...> Tail(std::tuple< T0, Tn...> &&tuple)
std::vector< std::tuple< TKnot...> > draw(std::size_t ndraws, Generator &rng) const
NdArray< T > trapz(const NdArray< T > &array, const Iterator kbegin, const Iterator kend, int axis)
NdSampler(std::tuple< std::vector< TKnot >...> knots, const NdArray::NdArray< double > &grid)
typename _index_sequence_helper< N >::type _make_index_sequence
std::unique_ptr< T > make_unique(Args &&...args)
Constructs an object of type T and wraps it in a std::unique_ptr using args as the parameter list for...
T sum(const NdArray< T > &array)