diff options
-rw-r--r-- | unsupported/Eigen/CXX11/Tensor | 1 | ||||
-rw-r--r-- | unsupported/Eigen/CXX11/src/Tensor/README.md | 38 | ||||
-rw-r--r-- | unsupported/Eigen/CXX11/src/Tensor/TensorBase.h | 15 | ||||
-rw-r--r-- | unsupported/Eigen/CXX11/src/Tensor/TensorForwardDeclarations.h | 1 | ||||
-rw-r--r-- | unsupported/Eigen/CXX11/src/Tensor/TensorScan.h | 197 | ||||
-rw-r--r-- | unsupported/test/CMakeLists.txt | 1 | ||||
-rw-r--r-- | unsupported/test/cxx11_tensor_scan.cpp | 98 |
7 files changed, 351 insertions, 0 deletions
diff --git a/unsupported/Eigen/CXX11/Tensor b/unsupported/Eigen/CXX11/Tensor index 77431cfc9..859147404 100644 --- a/unsupported/Eigen/CXX11/Tensor +++ b/unsupported/Eigen/CXX11/Tensor @@ -114,6 +114,7 @@ typedef unsigned __int64 uint64_t; #include "src/Tensor/TensorForcedEval.h" #include "src/Tensor/TensorGenerator.h" #include "src/Tensor/TensorAssign.h" +#include "src/Tensor/TensorScan.h" #include "src/Tensor/TensorExecutor.h" #include "src/Tensor/TensorDevice.h" diff --git a/unsupported/Eigen/CXX11/src/Tensor/README.md b/unsupported/Eigen/CXX11/src/Tensor/README.md index eeca2f69e..fda33edda 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/README.md +++ b/unsupported/Eigen/CXX11/src/Tensor/README.md @@ -1168,6 +1168,44 @@ Reduce a tensor using a user-defined reduction operator. See ```SumReducer``` in TensorFunctors.h for information on how to implement a reduction operator. +## Scan Operations + +A *Scan* operation returns a tensor with the same dimensions as the original +tensor. The operation performs an inclusive scan along the specified +axis, which means it computes a running total along the axis for a given +reduction operation. +If the reduction operation corresponds to summation, then this computes the +prefix sum of the tensor along the given axis. + +Example: +dd a comment to this line + + // Create a tensor of 2 dimensions + Eigen::Tensor<int, 2> a(2, 3); + a.setValues({{1, 2, 3}, {4, 5, 6}}); + // Scan it along the second dimension (1) using summation + Eigen::Tensor<int, 2> b = a.cumsum(1); + // The result is a tensor with the same size as the input + cout << "a" << endl << a << endl << endl; + cout << "b" << endl << b << endl << endl; + => + a + 1 2 3 + 6 5 4 + + b + 1 3 6 + 4 9 15 + +### <Operation> cumsum(const Index& axis) + +Perform a scan by summing consecutive entries. + +### <Operation> cumprod(const Index& axis) + +Perform a scan by multiplying consecutive entries. + + ## Convolutions ### <Operation> convolve(const Kernel& kernel, const Dimensions& dims) diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h b/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h index 07dcfa556..1eaa8d4fc 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h @@ -453,6 +453,21 @@ class TensorBase<Derived, ReadOnlyAccessors> return TensorFFTOp<const FFT, const Derived, FFTDataType, FFTDirection>(derived(), fft); } + // Scan. + typedef TensorScanOp<internal::SumReducer<CoeffReturnType>, const Derived> TensorScanSumOp; + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const TensorScanSumOp + cumsum(const Index& axis) const { + return TensorScanSumOp(derived(), axis); + } + + typedef TensorScanOp<internal::ProdReducer<CoeffReturnType>, const Derived> TensorScanProdOp; + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const TensorScanProdOp + cumprod(const Index& axis) const { + return TensorScanProdOp(derived(), axis); + } + // Reductions. template <typename Dims> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const TensorReductionOp<internal::SumReducer<CoeffReturnType>, const Dims, const Derived> diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorForwardDeclarations.h b/unsupported/Eigen/CXX11/src/Tensor/TensorForwardDeclarations.h index 5e59c7dee..a1a18d938 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorForwardDeclarations.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorForwardDeclarations.h @@ -46,6 +46,7 @@ template<typename StartIndices, typename StopIndices, typename Strides, typename template<typename Strides, typename XprType> class TensorInflationOp; template<typename Generator, typename XprType> class TensorGeneratorOp; template<typename LeftXprType, typename RightXprType> class TensorAssignOp; +template<typename Op, typename XprType> class TensorScanOp; template<typename CustomUnaryFunc, typename XprType> class TensorCustomUnaryOp; template<typename CustomBinaryFunc, typename LhsXprType, typename RhsXprType> class TensorCustomBinaryOp; diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorScan.h b/unsupported/Eigen/CXX11/src/Tensor/TensorScan.h new file mode 100644 index 000000000..031dbf6f2 --- /dev/null +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorScan.h @@ -0,0 +1,197 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Igor Babuschkin <igor@babuschk.in> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_CXX11_TENSOR_TENSOR_SCAN_H +#define EIGEN_CXX11_TENSOR_TENSOR_SCAN_H +namespace Eigen { + +namespace internal { +template <typename Op, typename XprType> +struct traits<TensorScanOp<Op, XprType> > + : public traits<XprType> { + typedef typename XprType::Scalar Scalar; + typedef traits<XprType> XprTraits; + typedef typename XprTraits::StorageKind StorageKind; + typedef typename XprType::Nested Nested; + typedef typename remove_reference<Nested>::type _Nested; + static const int NumDimensions = XprTraits::NumDimensions; + static const int Layout = XprTraits::Layout; +}; + +template<typename Op, typename XprType> +struct eval<TensorScanOp<Op, XprType>, Eigen::Dense> +{ + typedef const TensorScanOp<Op, XprType>& type; +}; + +template<typename Op, typename XprType> +struct nested<TensorScanOp<Op, XprType>, 1, + typename eval<TensorScanOp<Op, XprType> >::type> +{ + typedef TensorScanOp<Op, XprType> type; +}; +} // end namespace internal + +/** \class TensorScan + * \ingroup CXX11_Tensor_Module + * + * \brief Tensor scan class. + * + */ + +template <typename Op, typename XprType> +class TensorScanOp + : public TensorBase<TensorScanOp<Op, XprType>, ReadOnlyAccessors> { +public: + typedef typename Eigen::internal::traits<TensorScanOp>::Scalar Scalar; + typedef typename Eigen::NumTraits<Scalar>::Real RealScalar; + typedef typename XprType::CoeffReturnType CoeffReturnType; + typedef typename Eigen::internal::nested<TensorScanOp>::type Nested; + typedef typename Eigen::internal::traits<TensorScanOp>::StorageKind StorageKind; + typedef typename Eigen::internal::traits<TensorScanOp>::Index Index; + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorScanOp( + const XprType& expr, const Index& axis, const Op& op = Op()) + : m_expr(expr), m_axis(axis), m_accumulator(op) {} + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const Index axis() const { return m_axis; } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const XprType& expression() const { return m_expr; } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const Op accumulator() const { return m_accumulator; } + +protected: + typename XprType::Nested m_expr; + const Index m_axis; + const Op m_accumulator; +}; + +// Eval as rvalue +template <typename Op, typename ArgType, typename Device> +struct TensorEvaluator<const TensorScanOp<Op, ArgType>, Device> { + + typedef TensorScanOp<Op, ArgType> XprType; + typedef typename XprType::Index Index; + static const int NumDims = internal::array_size<typename TensorEvaluator<ArgType, Device>::Dimensions>::value; + typedef DSizes<Index, NumDims> Dimensions; + typedef typename XprType::Scalar Scalar; + typedef typename XprType::CoeffReturnType CoeffReturnType; + typedef typename PacketType<CoeffReturnType, Device>::type PacketReturnType; + + enum { + IsAligned = false, + PacketAccess = (internal::packet_traits<Scalar>::size > 1), + BlockAccess = false, + Layout = TensorEvaluator<ArgType, Device>::Layout, + CoordAccess = false, + RawAccess = true + }; + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorEvaluator(const XprType& op, + const Device& device) + : m_impl(op.expression(), device), + m_device(device), + m_axis(op.axis()), + m_accumulator(op.accumulator()), + m_dimensions(m_impl.dimensions()), + m_size(m_dimensions[m_axis]), + m_stride(1), + m_output(NULL) { + + // Accumulating a scalar isn't supported. + EIGEN_STATIC_ASSERT(NumDims > 0, YOU_MADE_A_PROGRAMMING_MISTAKE); + eigen_assert(m_axis >= 0 && m_axis < NumDims); + + // Compute stride of scan axis + if (static_cast<int>(Layout) == static_cast<int>(ColMajor)) { + for (int i = 0; i < m_axis; ++i) { + m_stride = m_stride * m_dimensions[i]; + } + } else { + for (int i = NumDims - 1; i > m_axis; --i) { + m_stride = m_stride * m_dimensions[i]; + } + } + } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Dimensions& dimensions() const { + return m_dimensions; + } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool evalSubExprsIfNeeded(Scalar* data) { + m_impl.evalSubExprsIfNeeded(NULL); + if (data) { + accumulateTo(data); + return false; + } else { + m_output = static_cast<CoeffReturnType*>(m_device.allocate(dimensions().TotalSize() * sizeof(Scalar))); + accumulateTo(m_output); + return true; + } + } + + template<int LoadMode> + EIGEN_DEVICE_FUNC PacketReturnType packet(Index index) const { + return internal::ploadt<PacketReturnType, LoadMode>(m_output + index); + } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType* data() const + { + return m_output; + } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const + { + return m_output[index]; + } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void cleanup() { + if (m_output != NULL) { + m_device.deallocate(m_output); + m_output = NULL; + } + m_impl.cleanup(); + } + +protected: + TensorEvaluator<ArgType, Device> m_impl; + const Device& m_device; + const Index m_axis; + Op m_accumulator; + const Dimensions& m_dimensions; + const Index& m_size; + Index m_stride; + CoeffReturnType* m_output; + + // TODO(ibab) Parallelize this single-threaded implementation if desired + EIGEN_DEVICE_FUNC void accumulateTo(Scalar* data) { + // We fix the index along the scan axis to 0 and perform an + // scan per remaining entry. The iteration is split into two nested + // loops to avoid an integer division by keeping track of each idx1 and idx2. + for (Index idx1 = 0; idx1 < dimensions().TotalSize() / m_size; idx1 += m_stride) { + for (Index idx2 = 0; idx2 < m_stride; idx2++) { + // Calculate the starting offset for the scan + Index offset = idx1 * m_size + idx2; + + // Compute the prefix sum along the axis, starting at the calculated offset + CoeffReturnType accum = m_accumulator.initialize(); + for (Index idx3 = 0; idx3 < m_size; idx3++) { + Index curr = offset + idx3 * m_stride; + m_accumulator.reduce(m_impl.coeff(curr), &accum); + data[curr] = m_accumulator.finalize(accum); + } + } + } + } +}; + +} // end namespace Eigen + +#endif // EIGEN_CXX11_TENSOR_TENSOR_SCAN_H diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index 70a3e4565..2d65eb0cd 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -176,6 +176,7 @@ if(EIGEN_TEST_CXX11) ei_add_test(cxx11_tensor_custom_index) ei_add_test(cxx11_tensor_fft) ei_add_test(cxx11_tensor_ifft) + ei_add_test(cxx11_tensor_scan) endif() diff --git a/unsupported/test/cxx11_tensor_scan.cpp b/unsupported/test/cxx11_tensor_scan.cpp new file mode 100644 index 000000000..dbd3023d7 --- /dev/null +++ b/unsupported/test/cxx11_tensor_scan.cpp @@ -0,0 +1,98 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Igor Babuschkin <igor@babuschk.in> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include <limits> +#include <numeric> +#include <Eigen/CXX11/Tensor> + +using Eigen::Tensor; + +template <int DataLayout, typename Type=float> +static void test_1d_scan() +{ + int size = 50; + Tensor<Type, 1, DataLayout> tensor(size); + tensor.setRandom(); + Tensor<Type, 1, DataLayout> result = tensor.cumsum(0); + + VERIFY_IS_EQUAL(tensor.dimension(0), result.dimension(0)); + + float accum = 0; + for (int i = 0; i < size; i++) { + accum += tensor(i); + VERIFY_IS_EQUAL(result(i), accum); + } + + accum = 1; + result = tensor.cumprod(0); + for (int i = 0; i < size; i++) { + accum *= tensor(i); + VERIFY_IS_EQUAL(result(i), accum); + } +} + +template <int DataLayout, typename Type=float> +static void test_4d_scan() +{ + int size = 5; + Tensor<Type, 4, DataLayout> tensor(size, size, size, size); + tensor.setRandom(); + + Tensor<Type, 4, DataLayout> result(size, size, size, size); + + result = tensor.cumsum(0); + float accum = 0; + for (int i = 0; i < size; i++) { + accum += tensor(i, 0, 0, 0); + VERIFY_IS_EQUAL(result(i, 0, 0, 0), accum); + } + result = tensor.cumsum(1); + accum = 0; + for (int i = 0; i < size; i++) { + accum += tensor(0, i, 0, 0); + VERIFY_IS_EQUAL(result(0, i, 0, 0), accum); + } + result = tensor.cumsum(2); + accum = 0; + for (int i = 0; i < size; i++) { + accum += tensor(0, 0, i, 0); + VERIFY_IS_EQUAL(result(0, 0, i, 0), accum); + } + result = tensor.cumsum(3); + accum = 0; + for (int i = 0; i < size; i++) { + accum += tensor(0, 0, 0, i); + VERIFY_IS_EQUAL(result(0, 0, 0, i), accum); + } +} + +template <int DataLayout> +static void test_tensor_maps() { + int inputs[20]; + TensorMap<Tensor<int, 1, DataLayout> > tensor_map(inputs, 20); + tensor_map.setRandom(); + + Tensor<int, 1, DataLayout> result = tensor_map.cumsum(0); + + int accum = 0; + for (int i = 0; i < 20; ++i) { + accum += tensor_map(i); + VERIFY_IS_EQUAL(result(i), accum); + } +} + +void test_cxx11_tensor_scan() { + CALL_SUBTEST(test_1d_scan<ColMajor>()); + CALL_SUBTEST(test_1d_scan<RowMajor>()); + CALL_SUBTEST(test_4d_scan<ColMajor>()); + CALL_SUBTEST(test_4d_scan<RowMajor>()); + CALL_SUBTEST(test_tensor_maps<ColMajor>()); + CALL_SUBTEST(test_tensor_maps<RowMajor>()); +} |