xtensor
 
Loading...
Searching...
No Matches
xblockwise_reducer.hpp
1#ifndef XTENSOR_XBLOCKWISE_REDUCER_HPP
2#define XTENSOR_XBLOCKWISE_REDUCER_HPP
3
4#include "../core/xmultiindex_iterator.hpp"
5#include "../core/xshape.hpp"
6#include "../reducers/xblockwise_reducer_functors.hpp"
7#include "../reducers/xreducer.hpp"
8#include "xtl/xclosure.hpp"
9#include "xtl/xsequence.hpp"
10
11namespace xt
12{
13
14 template <class CT, class F, class X, class O>
15 class xblockwise_reducer
16 {
17 public:
18
19 using self_type = xblockwise_reducer<CT, F, X, O>;
20 using raw_options_type = std::decay_t<O>;
21 using keep_dims = xtl::mpl::contains<raw_options_type, xt::keep_dims_type>;
22 using xexpression_type = std::decay_t<CT>;
23 using shape_type = typename xreducer_shape_type<typename xexpression_type::shape_type, std::decay_t<X>, keep_dims>::type;
24 using functor_type = F;
25 using value_type = typename functor_type::value_type;
26 using input_shape_type = typename xexpression_type::shape_type;
27 using input_chunk_index_type = filter_fixed_shape_t<input_shape_type>;
28 using input_grid_strides = filter_fixed_shape_t<input_shape_type>;
29 using axes_type = X;
30 using chunk_shape_type = filter_fixed_shape_t<shape_type>;
31
32
33 template <class E, class BS, class XX, class OO, class FF>
34 xblockwise_reducer(E&& e, BS&& block_shape, XX&& axes, OO&& options, FF&& functor);
35
36 const input_shape_type& input_shape() const;
37 const axes_type& axes() const;
38
39 std::size_t dimension() const;
40
41 const shape_type& shape() const;
42
43 const chunk_shape_type& chunk_shape() const;
44
45 template <class R>
46 void assign_to(R& result) const;
47
48 private:
49
50 using mapping_type = filter_fixed_shape_t<shape_type>;
51 using input_chunked_view_type = xchunked_view<const std::decay_t<CT>&>;
52 using input_const_chunked_iterator_type = typename input_chunked_view_type::const_chunk_iterator;
53 using input_chunk_range_type = std::array<xmultiindex_iterator<input_chunk_index_type>, 2>;
54
55 template <class CI>
56 void assign_to_chunk(CI& result_chunk_iter) const;
57
58 template <class CI>
59 input_chunk_range_type compute_input_chunk_range(CI& result_chunk_iter) const;
60
61 input_const_chunked_iterator_type get_input_chunk_iter(input_chunk_index_type input_chunk_index) const;
62 void init_shapes();
63
64 CT m_e;
65 xchunked_view<const std::decay_t<CT>&> m_e_chunked_view;
66 axes_type m_axes;
67 raw_options_type m_options;
68 functor_type m_functor;
69 shape_type m_result_shape;
70 chunk_shape_type m_result_chunk_shape;
71 mapping_type m_mapping;
72 input_grid_strides m_input_grid_strides;
73 };
74
75 template <class CT, class F, class X, class O>
76 template <class E, class BS, class XX, class OO, class FF>
77 xblockwise_reducer<CT, F, X, O>::xblockwise_reducer(E&& e, BS&& block_shape, XX&& axes, OO&& options, FF&& functor)
78 : m_e(std::forward<E>(e))
79 , m_e_chunked_view(m_e, std::forward<BS>(block_shape))
80 , m_axes(std::forward<XX>(axes))
81 , m_options(std::forward<OO>(options))
82 , m_functor(std::forward<FF>(functor))
83 , m_result_shape()
84 , m_result_chunk_shape()
85 , m_mapping()
86 , m_input_grid_strides()
87 {
88 init_shapes();
89 resize_container(m_input_grid_strides, m_e.dimension());
90 std::size_t stride = 1;
91
92 for (std::size_t i = m_input_grid_strides.size(); i != 0; --i)
93 {
94 m_input_grid_strides[i - 1] = stride;
95 stride *= m_e_chunked_view.grid_shape()[i - 1];
96 }
97 }
98
99 template <class CT, class F, class X, class O>
100 inline auto xblockwise_reducer<CT, F, X, O>::input_shape() const -> const input_shape_type&
101 {
102 return m_e.shape();
103 }
104
105 template <class CT, class F, class X, class O>
106 inline auto xblockwise_reducer<CT, F, X, O>::axes() const -> const axes_type&
107 {
108 return m_axes;
109 }
110
111 template <class CT, class F, class X, class O>
112 inline std::size_t xblockwise_reducer<CT, F, X, O>::dimension() const
113 {
114 return m_result_shape.size();
115 }
116
117 template <class CT, class F, class X, class O>
118 inline auto xblockwise_reducer<CT, F, X, O>::shape() const -> const shape_type&
119 {
120 return m_result_shape;
121 }
122
123 template <class CT, class F, class X, class O>
124 inline auto xblockwise_reducer<CT, F, X, O>::chunk_shape() const -> const chunk_shape_type&
125 {
126 return m_result_chunk_shape;
127 }
128
129 template <class CT, class F, class X, class O>
130 template <class R>
131 inline void xblockwise_reducer<CT, F, X, O>::assign_to(R& result) const
132 {
133 auto result_chunked_view = as_chunked(result, m_result_chunk_shape);
134 for (auto chunk_iter = result_chunked_view.chunk_begin(); chunk_iter != result_chunked_view.chunk_end();
135 ++chunk_iter)
136 {
137 assign_to_chunk(chunk_iter);
138 }
139 }
140
141 template <class CT, class F, class X, class O>
142 auto xblockwise_reducer<CT, F, X, O>::get_input_chunk_iter(input_chunk_index_type input_chunk_index) const
143 -> input_const_chunked_iterator_type
144 {
145 std::size_t chunk_linear_index = 0;
146 for (std::size_t i = 0; i < m_e_chunked_view.dimension(); ++i)
147 {
148 chunk_linear_index += input_chunk_index[i] * m_input_grid_strides[i];
149 }
150 return input_const_chunked_iterator_type(m_e_chunked_view, std::move(input_chunk_index), chunk_linear_index);
151 }
152
153 template <class CT, class F, class X, class O>
154 template <class CI>
155 void xblockwise_reducer<CT, F, X, O>::assign_to_chunk(CI& result_chunk_iter) const
156 {
157 auto result_chunk_view = *result_chunk_iter;
158 auto reduction_variable = m_functor.reduction_variable(result_chunk_view);
159
160 // get the range of input chunks we need to compute the desired ouput chunk
161 auto range = compute_input_chunk_range(result_chunk_iter);
162
163 // iterate over input chunk (indics)
164 auto first = true;
165 // std::for_each(std::get<0>(range), std::get<1>(range), [&](auto && input_chunk_index)
166 auto iter = std::get<0>(range);
167 while (iter != std::get<1>(range))
168 {
169 const auto& input_chunk_index = *iter;
170 // get input chunk iterator from chunk index
171 auto chunked_input_iter = this->get_input_chunk_iter(input_chunk_index);
172 auto input_chunk_view = *chunked_input_iter;
173
174 // compute the per block result
175 auto block_res = m_functor.compute(input_chunk_view, m_axes, m_options);
176
177 // merge
178 m_functor.merge(block_res, first, result_chunk_view, reduction_variable);
179 first = false;
180 ++iter;
181 }
182
183 // finalize (ie smth like normalization)
184 m_functor.finalize(reduction_variable, result_chunk_view, *this);
185 }
186
187 template <class CT, class F, class X, class O>
188 template <class CI>
189 auto xblockwise_reducer<CT, F, X, O>::compute_input_chunk_range(CI& result_chunk_iter) const
190 -> input_chunk_range_type
191 {
192 auto input_chunks_begin = xtl::make_sequence<input_chunk_index_type>(m_e_chunked_view.dimension(), 0);
193 auto input_chunks_end = xtl::make_sequence<input_chunk_index_type>(m_e_chunked_view.dimension());
194
195 XTENSOR_ASSERT(input_chunks_begin.size() == m_e_chunked_view.dimension());
196 XTENSOR_ASSERT(input_chunks_end.size() == m_e_chunked_view.dimension());
197
198 std::copy(
199 m_e_chunked_view.grid_shape().begin(),
200 m_e_chunked_view.grid_shape().end(),
201 input_chunks_end.begin()
202 );
203
204 const auto& chunk_index = result_chunk_iter.chunk_index();
205 for (std::size_t result_ax_index = 0; result_ax_index < m_result_shape.size(); ++result_ax_index)
206 {
207 if (m_result_shape[result_ax_index] != 1)
208 {
209 const auto input_ax_index = m_mapping[result_ax_index];
210 input_chunks_begin[input_ax_index] = chunk_index[result_ax_index];
211 input_chunks_end[input_ax_index] = chunk_index[result_ax_index] + 1;
212 }
213 }
214 return input_chunk_range_type{
215 multiindex_iterator_begin<input_chunk_index_type>(input_chunks_begin, input_chunks_end),
216 multiindex_iterator_end<input_chunk_index_type>(input_chunks_begin, input_chunks_end)
217 };
218 }
219
220 template <class CT, class F, class X, class O>
221 void xblockwise_reducer<CT, F, X, O>::init_shapes()
222 {
223 const auto& shape = m_e.shape();
224 const auto dimension = m_e.dimension();
225 const auto& block_shape = m_e_chunked_view.chunk_shape();
226 if (xtl::mpl::contains<raw_options_type, xt::keep_dims_type>::value)
227 {
228 resize_container(m_result_shape, dimension);
229 resize_container(m_result_chunk_shape, dimension);
230 resize_container(m_mapping, dimension);
231 for (std::size_t i = 0; i < dimension; ++i)
232 {
233 m_mapping[i] = i;
234 if (std::find(m_axes.begin(), m_axes.end(), i) == m_axes.end())
235 {
236 // i not in m_axes!
237 m_result_shape[i] = shape[i];
238 m_result_chunk_shape[i] = block_shape[i];
239 }
240 else
241 {
242 m_result_shape[i] = 1;
243 m_result_chunk_shape[i] = 1;
244 }
245 }
246 }
247 else
248 {
249 const auto result_dim = dimension - m_axes.size();
250 resize_container(m_result_shape, result_dim);
251 resize_container(m_result_chunk_shape, result_dim);
252 resize_container(m_mapping, result_dim);
253
254 for (std::size_t i = 0, idx = 0; i < dimension; ++i)
255 {
256 if (std::find(m_axes.begin(), m_axes.end(), i) == m_axes.end())
257 {
258 // i not in axes!
259 m_result_shape[idx] = shape[i];
260 m_result_chunk_shape[idx] = block_shape[i];
261 m_mapping[idx] = i;
262 ++idx;
263 }
264 }
265 }
266 }
267
268 template <class E, class CS, class A, class O, class FF>
269 inline auto blockwise_reducer(E&& e, CS&& chunk_shape, A&& axes, O&& raw_options, FF&& functor)
270 {
271 using functor_type = std::decay_t<FF>;
272 using closure_type = xtl::const_closure_type_t<E>;
273 using axes_type = std::decay_t<A>;
274
275 return xblockwise_reducer<closure_type, functor_type, axes_type, O>(
276 std::forward<E>(e),
277 std::forward<CS>(chunk_shape),
278 std::forward<A>(axes),
279 std::forward<O>(raw_options),
280 std::forward<FF>(functor)
281 );
282 }
283
284 namespace blockwise
285 {
286
287#define XTENSOR_BLOCKWISE_REDUCER_FUNC(FNAME, FUNCTOR) \
288 template < \
289 class T = void, \
290 class E, \
291 class BS, \
292 class X, \
293 class O = DEFAULT_STRATEGY_REDUCERS, \
294 XTL_REQUIRES(std::negation<is_reducer_options<X>>, std::negation<xtl::is_integral<std::decay_t<X>>>)> \
295 auto FNAME(E&& e, BS&& block_shape, X&& axes, O options = O()) \
296 { \
297 using input_expression_type = std::decay_t<E>; \
298 using functor_type = FUNCTOR<typename input_expression_type::value_type, T>; \
299 return blockwise_reducer( \
300 std::forward<E>(e), \
301 std::forward<BS>(block_shape), \
302 std::forward<X>(axes), \
303 std::forward<O>(options), \
304 functor_type() \
305 ); \
306 } \
307 template < \
308 class T = void, \
309 class E, \
310 class BS, \
311 class X, \
312 class O = DEFAULT_STRATEGY_REDUCERS, \
313 XTL_REQUIRES(xtl::is_integral<std::decay_t<X>>)> \
314 auto FNAME(E&& e, BS&& block_shape, X axis, O options = O()) \
315 { \
316 std::array<X, 1> axes{axis}; \
317 using input_expression_type = std::decay_t<E>; \
318 using functor_type = FUNCTOR<typename input_expression_type::value_type, T>; \
319 return blockwise_reducer( \
320 std::forward<E>(e), \
321 std::forward<BS>(block_shape), \
322 axes, \
323 std::forward<O>(options), \
324 functor_type() \
325 ); \
326 } \
327 template < \
328 class T = void, \
329 class E, \
330 class BS, \
331 class O = DEFAULT_STRATEGY_REDUCERS, \
332 XTL_REQUIRES(is_reducer_options<O>, std::negation<xtl::is_integral<std::decay_t<O>>>)> \
333 auto FNAME(E&& e, BS&& block_shape, O options = O()) \
334 { \
335 using input_expression_type = std::decay_t<E>; \
336 using axes_type = filter_fixed_shape_t<typename input_expression_type::shape_type>; \
337 axes_type axes = xtl::make_sequence<axes_type>(e.dimension()); \
338 XTENSOR_ASSERT(axes.size() == e.dimension()); \
339 std::iota(axes.begin(), axes.end(), 0); \
340 using functor_type = FUNCTOR<typename input_expression_type::value_type, T>; \
341 return blockwise_reducer( \
342 std::forward<E>(e), \
343 std::forward<BS>(block_shape), \
344 std::move(axes), \
345 std::forward<O>(options), \
346 functor_type() \
347 ); \
348 } \
349 template <class T = void, class E, class BS, class I, std::size_t N, class O = DEFAULT_STRATEGY_REDUCERS> \
350 auto FNAME(E&& e, BS&& block_shape, const I(&axes)[N], O options = O()) \
351 { \
352 using input_expression_type = std::decay_t<E>; \
353 using functor_type = FUNCTOR<typename input_expression_type::value_type, T>; \
354 using axes_type = std::array<std::size_t, N>; \
355 auto ax = xt::forward_normalize<axes_type>(e, axes); \
356 return blockwise_reducer( \
357 std::forward<E>(e), \
358 std::forward<BS>(block_shape), \
359 std::move(ax), \
360 std::forward<O>(options), \
361 functor_type() \
362 ); \
363 }
364 XTENSOR_BLOCKWISE_REDUCER_FUNC(sum, xt::detail::blockwise::sum_functor)
365 XTENSOR_BLOCKWISE_REDUCER_FUNC(prod, xt::detail::blockwise::prod_functor)
366 XTENSOR_BLOCKWISE_REDUCER_FUNC(amin, xt::detail::blockwise::amin_functor)
367 XTENSOR_BLOCKWISE_REDUCER_FUNC(amax, xt::detail::blockwise::amax_functor)
368 XTENSOR_BLOCKWISE_REDUCER_FUNC(mean, xt::detail::blockwise::mean_functor)
369 XTENSOR_BLOCKWISE_REDUCER_FUNC(variance, xt::detail::blockwise::variance_functor)
370 XTENSOR_BLOCKWISE_REDUCER_FUNC(stddev, xt::detail::blockwise::stddev_functor)
371
372#undef XTENSOR_BLOCKWISE_REDUCER_FUNC
373
374
375// norm reducers do *not* allow to to pass a template
376// parameter to specifiy the internal computation type
377#define XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC(FNAME, FUNCTOR) \
378 template < \
379 class E, \
380 class BS, \
381 class X, \
382 class O = DEFAULT_STRATEGY_REDUCERS, \
383 XTL_REQUIRES(std::negation<is_reducer_options<X>>, std::negation<xtl::is_integral<std::decay_t<X>>>)> \
384 auto FNAME(E&& e, BS&& block_shape, X&& axes, O options = O()) \
385 { \
386 using input_expression_type = std::decay_t<E>; \
387 using functor_type = FUNCTOR<typename input_expression_type::value_type>; \
388 return blockwise_reducer( \
389 std::forward<E>(e), \
390 std::forward<BS>(block_shape), \
391 std::forward<X>(axes), \
392 std::forward<O>(options), \
393 functor_type() \
394 ); \
395 } \
396 template <class E, class BS, class X, class O = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(xtl::is_integral<std::decay_t<X>>)> \
397 auto FNAME(E&& e, BS&& block_shape, X axis, O options = O()) \
398 { \
399 std::array<X, 1> axes{axis}; \
400 using input_expression_type = std::decay_t<E>; \
401 using functor_type = FUNCTOR<typename input_expression_type::value_type>; \
402 return blockwise_reducer( \
403 std::forward<E>(e), \
404 std::forward<BS>(block_shape), \
405 axes, \
406 std::forward<O>(options), \
407 functor_type() \
408 ); \
409 } \
410 template < \
411 class E, \
412 class BS, \
413 class O = DEFAULT_STRATEGY_REDUCERS, \
414 XTL_REQUIRES(is_reducer_options<O>, std::negation<xtl::is_integral<std::decay_t<O>>>)> \
415 auto FNAME(E&& e, BS&& block_shape, O options = O()) \
416 { \
417 using input_expression_type = std::decay_t<E>; \
418 using axes_type = filter_fixed_shape_t<typename input_expression_type::shape_type>; \
419 axes_type axes = xtl::make_sequence<axes_type>(e.dimension()); \
420 XTENSOR_ASSERT(axes.size() == e.dimension()); \
421 std::iota(axes.begin(), axes.end(), 0); \
422 using functor_type = FUNCTOR<typename input_expression_type::value_type>; \
423 return blockwise_reducer( \
424 std::forward<E>(e), \
425 std::forward<BS>(block_shape), \
426 std::move(axes), \
427 std::forward<O>(options), \
428 functor_type() \
429 ); \
430 } \
431 template <class E, class BS, class I, std::size_t N, class O = DEFAULT_STRATEGY_REDUCERS> \
432 auto FNAME(E&& e, BS&& block_shape, const I(&axes)[N], O options = O()) \
433 { \
434 using input_expression_type = std::decay_t<E>; \
435 using functor_type = FUNCTOR<typename input_expression_type::value_type>; \
436 using axes_type = std::array<std::size_t, N>; \
437 auto ax = xt::forward_normalize<axes_type>(e, axes); \
438 return blockwise_reducer( \
439 std::forward<E>(e), \
440 std::forward<BS>(block_shape), \
441 std::move(ax), \
442 std::forward<O>(options), \
443 functor_type() \
444 ); \
445 }
446 XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC(norm_l0, xt::detail::blockwise::norm_l0_functor)
447 XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC(norm_l1, xt::detail::blockwise::norm_l1_functor)
448 XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC(norm_l2, xt::detail::blockwise::norm_l2_functor)
449 XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC(norm_sq, xt::detail::blockwise::norm_sq_functor)
450 XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC(norm_linf, xt::detail::blockwise::norm_linf_functor)
451
452#undef XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC
453
454
455#define XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC(FNAME, FUNCTOR) \
456 template < \
457 class E, \
458 class BS, \
459 class X, \
460 class O = DEFAULT_STRATEGY_REDUCERS, \
461 XTL_REQUIRES(std::negation<is_reducer_options<X>>, std::negation<xtl::is_integral<std::decay_t<X>>>)> \
462 auto FNAME(E&& e, BS&& block_shape, double p, X&& axes, O options = O()) \
463 { \
464 using input_expression_type = std::decay_t<E>; \
465 using functor_type = FUNCTOR<typename input_expression_type::value_type>; \
466 return blockwise_reducer( \
467 std::forward<E>(e), \
468 std::forward<BS>(block_shape), \
469 std::forward<X>(axes), \
470 std::forward<O>(options), \
471 functor_type(p) \
472 ); \
473 } \
474 template <class E, class BS, class X, class O = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(xtl::is_integral<std::decay_t<X>>)> \
475 auto FNAME(E&& e, BS&& block_shape, double p, X axis, O options = O()) \
476 { \
477 std::array<X, 1> axes{axis}; \
478 using input_expression_type = std::decay_t<E>; \
479 using functor_type = FUNCTOR<typename input_expression_type::value_type>; \
480 return blockwise_reducer( \
481 std::forward<E>(e), \
482 std::forward<BS>(block_shape), \
483 axes, \
484 std::forward<O>(options), \
485 functor_type(p) \
486 ); \
487 } \
488 template < \
489 class E, \
490 class BS, \
491 class O = DEFAULT_STRATEGY_REDUCERS, \
492 XTL_REQUIRES(is_reducer_options<O>, std::negation<xtl::is_integral<std::decay_t<O>>>)> \
493 auto FNAME(E&& e, BS&& block_shape, double p, O options = O()) \
494 { \
495 using input_expression_type = std::decay_t<E>; \
496 using axes_type = filter_fixed_shape_t<typename input_expression_type::shape_type>; \
497 axes_type axes = xtl::make_sequence<axes_type>(e.dimension()); \
498 XTENSOR_ASSERT(axes.size() == e.dimension()); \
499 std::iota(axes.begin(), axes.end(), 0); \
500 using functor_type = FUNCTOR<typename input_expression_type::value_type>; \
501 return blockwise_reducer( \
502 std::forward<E>(e), \
503 std::forward<BS>(block_shape), \
504 std::move(axes), \
505 std::forward<O>(options), \
506 functor_type(p) \
507 ); \
508 } \
509 template <class E, class BS, class I, std::size_t N, class O = DEFAULT_STRATEGY_REDUCERS> \
510 auto FNAME(E&& e, BS&& block_shape, double p, const I(&axes)[N], O options = O()) \
511 { \
512 using input_expression_type = std::decay_t<E>; \
513 using functor_type = FUNCTOR<typename input_expression_type::value_type>; \
514 using axes_type = std::array<std::size_t, N>; \
515 auto ax = xt::forward_normalize<axes_type>(e, axes); \
516 return blockwise_reducer( \
517 std::forward<E>(e), \
518 std::forward<BS>(block_shape), \
519 std::move(ax), \
520 std::forward<O>(options), \
521 functor_type(p) \
522 ); \
523 }
524
525 XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC(norm_lp_to_p, xt::detail::blockwise::norm_lp_to_p_functor);
526 XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC(norm_lp, xt::detail::blockwise::norm_lp_functor);
527
528#undef XTENSOR_BLOCKWISE_NORM_REDUCER_FUNC
529 }
530
531}
532
533#endif
xt::xblockwise_reducer
Definition xblockwise_reducer.hpp:16
xt::xchunked_view
Definition xchunked_view.hpp:37
xt
standard mathematical functions for xexpressions
Definition xchunked_array.hpp:20
xt::range
auto range(A start_val, B stop_val)
Select a range from start_val to stop_val (excluded).
Definition xslice.hpp:744
xt::xreducer_shape_type
Definition xreducer.hpp:1239