From f55afe67215a386029490d14b87fd70d0d904eb9 Mon Sep 17 00:00:00 2001
From: Stefan Wunsch <stefan.wunsch@cern.ch>
Date: Wed, 19 Dec 2018 18:53:44 +0100
Subject: [PATCH] [TMVA::Experimental] Introduce RTensor class

---
 tmva/tmva/CMakeLists.txt         |   3 +
 tmva/tmva/inc/TMVA/RTensor.hxx   | 429 +++++++++++++++++++++++++++++++
 tmva/tmva/test/CMakeLists.txt    |   1 +
 tmva/tmva/test/rtensor.cxx       | 353 +++++++++++++++++++++++++
 tutorials/tmva/tmva001_RTensor.C |  54 ++++
 5 files changed, 840 insertions(+)
 create mode 100644 tmva/tmva/inc/TMVA/RTensor.hxx
 create mode 100644 tmva/tmva/test/rtensor.cxx
 create mode 100644 tutorials/tmva/tmva001_RTensor.C

diff --git a/tmva/tmva/CMakeLists.txt b/tmva/tmva/CMakeLists.txt
index dcd496b581f..e07319d6abc 100644
--- a/tmva/tmva/CMakeLists.txt
+++ b/tmva/tmva/CMakeLists.txt
@@ -178,6 +178,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(TMVA
     TMVA/VarTransformHandler.h
     TMVA/Version.h
     TMVA/Volume.h
+    TMVA/RTensor.hxx
   SOURCES
     src/BDTEventWrapper.cxx
     src/BinarySearchTree.cxx
@@ -355,6 +356,8 @@ ROOT_STANDARD_LIBRARY_PACKAGE(TMVA
   DEPENDENCIES
     TreePlayer
     Tree
+    ROOTDataFrame
+    ROOTVecOps
     Hist
     Matrix
     Minuit
diff --git a/tmva/tmva/inc/TMVA/RTensor.hxx b/tmva/tmva/inc/TMVA/RTensor.hxx
new file mode 100644
index 00000000000..e0c6cc86682
--- /dev/null
+++ b/tmva/tmva/inc/TMVA/RTensor.hxx
@@ -0,0 +1,429 @@
+#ifndef TMVA_RTENSOR
+#define TMVA_RTENSOR
+
+#include <vector>
+#include <cstddef>     // std::size_t
+#include <stdexcept>   // std::runtime_error
+#include <sstream>     // std::stringstream
+#include <memory>      // std::shared_ptr
+#include <type_traits> // std::is_convertible
+#include <algorithm>   // std::reverse
+
+namespace TMVA {
+namespace Experimental {
+
+/// Memory layout type
+enum class MemoryLayout : uint8_t {
+   RowMajor = 0x01,
+   ColumnMajor = 0x02
+};
+
+namespace Internal {
+
+/// \brief Get size of tensor from shape vector
+/// \param[in] shape Shape vector
+/// \return Size of contiguous memory
+template <typename T>
+std::size_t GetSizeFromShape(const T &shape)
+{
+   if (shape.size() == 0)
+      return 0;
+   std::size_t size = 1;
+   for (auto &s : shape)
+      size *= s;
+   return size;
+}
+
+/// \brief Compute strides from shape vector.
+/// \param[in] shape Shape vector
+/// \param[in] layout Memory layout
+/// \return Size of contiguous memory
+///
+/// This information is needed for the multi-dimensional indexing. See here:
+/// https://en.wikipedia.org/wiki/Row-_and_column-major_order
+/// https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.strides.html
+template <typename T>
+std::vector<std::size_t> ComputeStridesFromShape(const T &shape, MemoryLayout layout)
+{
+   const auto size = shape.size();
+   T strides(size);
+   if (layout == MemoryLayout::RowMajor) {
+      for (std::size_t i = 0; i < size; i++) {
+         if (i == 0) {
+            strides[size - 1 - i] = 1;
+         } else {
+            strides[size - 1 - i] = strides[size - 1 - i + 1] * shape[size - 1 - i + 1];
+         }
+      }
+   } else if (layout == MemoryLayout::ColumnMajor) {
+      for (std::size_t i = 0; i < size; i++) {
+         if (i == 0) {
+            strides[i] = 1;
+         } else {
+            strides[i] = strides[i - 1] * shape[i - 1];
+         }
+      }
+   } else {
+      std::stringstream ss;
+      ss << "Memory layout type is not valid for calculating strides.";
+      throw std::runtime_error(ss.str());
+   }
+   return strides;
+}
+
+/// \brief Type checking for all types of a parameter pack, e.g., used in combination with std::is_convertible
+template <class... Ts>
+struct and_types : std::true_type {
+};
+
+template <class T0, class... Ts>
+struct and_types<T0, Ts...> : std::integral_constant<bool, T0{} && and_types<Ts...>{}> {
+};
+
+} // namespace TMVA::Experimental::Internal
+
+/// \class TMVA::Experimental::RTensor
+/// \brief RTensor is a container with contiguous memory and shape information.
+/// \tparam T Data-type of the tensor
+///
+/// An RTensor is a vector-like container, which has additional shape information.
+/// The elements of the multi-dimensional container can be accessed by their
+/// indices in a coherent way without taking care about the one-dimensional memory
+/// layout of the contiguous storage. This also allows to manipulate the shape
+/// of the container without moving the actual elements in memory. Another feature
+/// is that an RTensor can own the underlying contiguous memory but can also represent
+/// only a view on existing data without owning it.
+template <typename V, typename C = std::vector<V>>
+class RTensor {
+public:
+   // Typedefs
+   using Value_t = V;
+   using Shape_t = std::vector<std::size_t>;
+   using Index_t = Shape_t;
+   using Slice_t = std::vector<Shape_t>;
+   using Container_t = C;
+
+private:
+   Shape_t fShape;
+   Shape_t fStrides;
+   std::size_t fSize;
+   MemoryLayout fLayout;
+   Value_t *fData;
+   std::shared_ptr<Container_t> fContainer;
+
+public:
+   // Constructors
+
+   /// \brief Construct a tensor as view on data
+   /// \param[in] data Pointer to data contiguous in memory
+   /// \param[in] shape Shape vector
+   /// \param[in] layout Memory layout
+   RTensor(Value_t *data, Shape_t shape, MemoryLayout layout = MemoryLayout::RowMajor)
+      : fShape(shape), fLayout(layout), fData(data), fContainer(NULL)
+   {
+      fSize = Internal::GetSizeFromShape(shape);
+      fStrides = Internal::ComputeStridesFromShape(shape, layout);
+   }
+
+   /// \brief Construct a tensor owning externally provided data
+   /// \param[in] container Shared pointer to data container
+   /// \param[in] shape Shape vector
+   /// \param[in] layout Memory layout
+   RTensor(std::shared_ptr<Container_t> container, Shape_t shape,
+           MemoryLayout layout = MemoryLayout::RowMajor)
+      : fShape(shape), fLayout(layout), fContainer(container)
+   {
+      fSize = Internal::GetSizeFromShape(shape);
+      fStrides = Internal::ComputeStridesFromShape(shape, layout);
+      fData = &(*container->begin());
+   }
+
+   /// \brief Construct a tensor owning data initialized with new container
+   /// \param[in] shape Shape vector
+   /// \param[in] layout Memory layout
+   RTensor(Shape_t shape, MemoryLayout layout = MemoryLayout::RowMajor)
+      : fShape(shape), fLayout(layout)
+   {
+      // TODO: Document how data pointer is determined using STL iterator interface.
+      // TODO: Sanitize given container type with type traits
+      fSize = Internal::GetSizeFromShape(shape);
+      fStrides = Internal::ComputeStridesFromShape(shape, layout);
+      fContainer = std::make_shared<Container_t>(fSize);
+      fData = &(*fContainer->begin());
+   }
+
+   // Access elements
+   Value_t &operator()(const Index_t &idx);
+   template <typename... Idx> Value_t &operator()(Idx... idx);
+
+   // Access properties
+   std::size_t GetSize() { return fSize; }
+   Shape_t GetShape() { return fShape; }
+   Shape_t GetStrides() { return fStrides; }
+   Value_t *GetData() { return fData; }
+   std::shared_ptr<Container_t> GetContainer() { return fContainer; }
+   MemoryLayout GetMemoryLayout() { return fLayout; }
+   bool IsView() { return fContainer == NULL; }
+   bool IsOwner() { return !IsView(); }
+
+   // Transformations
+   RTensor<Value_t, Container_t> Transpose();
+   RTensor<Value_t, Container_t> Squeeze();
+   RTensor<Value_t, Container_t> ExpandDims(int idx);
+   RTensor<Value_t, Container_t> Reshape(const Shape_t &shape);
+   RTensor<Value_t, Container_t> Slice(const Slice_t &slice);
+};
+
+/// \brief Access elements
+/// \param[in] idx Index vector
+/// \return Reference to element
+template <typename Value_t, typename Container_t>
+inline Value_t &RTensor<Value_t, Container_t>::operator()(const Index_t &idx)
+{
+   std::size_t globalIndex = 0;
+   const auto size = idx.size();
+   for (std::size_t i = 0; i < size; i++) {
+      globalIndex += fStrides[size - 1 - i] * idx[size - 1 - i];
+   }
+   return *(fData + globalIndex);
+}
+
+/// \brief Access elements
+/// \param[in] idx Indices
+/// \return Reference to element
+template <typename Value_t, typename Container_t>
+template <typename... Idx>
+Value_t &RTensor<Value_t, Container_t>::operator()(Idx... idx)
+{
+   static_assert(Internal::and_types<std::is_convertible<Idx, std::size_t>...>{},
+                 "Indices are not convertible to std::size_t.");
+   return this->operator()({static_cast<std::size_t>(idx)...});
+}
+
+/// \brief Transpose
+/// \returns New RTensor
+/// The tensor is transposed by inverting the associated memory layout from row-
+/// major to column-major and vice versa. Therefore, the underlying data is not
+/// touched.
+template <typename Value_t, typename Container_t>
+inline RTensor<Value_t, Container_t> RTensor<Value_t, Container_t>::Transpose()
+{
+   // Transpose by inverting memory layout
+   if (fLayout == MemoryLayout::RowMajor) {
+      fLayout = MemoryLayout::ColumnMajor;
+   } else if (fLayout == MemoryLayout::ColumnMajor) {
+      fLayout = MemoryLayout::RowMajor;
+   } else {
+      std::runtime_error("Memory layout is not known.");
+   }
+
+   // Create copy of container
+   RTensor<Value_t, Container_t> x(*this);
+
+   // Reverse shape
+   std::reverse(x.fShape.begin(), x.fShape.end());
+
+   // Reverse strides
+   std::reverse(x.fStrides.begin(), x.fStrides.end());
+
+   return x;
+}
+
+/// \brief Squeeze dimensions
+/// \returns New RTensor
+/// Squeeze removes the dimensions of size one from the shape.
+template <typename Value_t, typename Container_t>
+inline RTensor<Value_t, Container_t> RTensor<Value_t, Container_t>::Squeeze()
+{
+   // Remove dimensions of one and associated strides
+   Shape_t shape;
+   Shape_t strides;
+   for (std::size_t i = 0; i < fShape.size(); i++) {
+      if (fShape[i] != 1) {
+         shape.emplace_back(fShape[i]);
+         strides.emplace_back(fStrides[i]);
+      }
+   }
+
+   // If all dimensions are 1, we need to keep one.
+   // This does not apply if the inital shape is already empty. Then, return
+   // the empty shape.
+   if (shape.size() == 0 && fShape.size() != 0) {
+      shape.emplace_back(1);
+      strides.emplace_back(1);
+   }
+
+   // Create copy, attach new shape and strides and return
+   RTensor<Value_t, Container_t> x(*this);
+   x.fShape = shape;
+   x.fStrides = strides;
+   return x;
+}
+
+/// \brief Expand dimensions
+/// \param[in] idx Index in shape vector where dimension is added
+/// \returns New RTensor
+/// Inserts a dimension of one into the shape.
+template <typename Value_t, typename Container_t>
+inline RTensor<Value_t, Container_t> RTensor<Value_t, Container_t>::ExpandDims(int idx)
+{
+   // Compose shape vector with additional dimensions and adjust strides
+   const int len = fShape.size();
+   auto shape = fShape;
+   auto strides = fStrides;
+   if (idx < 0) {
+      if (len + idx + 1 < 0) {
+         throw std::runtime_error("Given negative index is invalid.");
+      }
+      shape.insert(shape.end() + 1 + idx, 1);
+      strides.insert(strides.begin() + 1 + idx, 1);
+   } else {
+      if (idx > len) {
+         throw std::runtime_error("Given index is invalid.");
+      }
+      shape.insert(shape.begin() + idx, 1);
+      strides.insert(strides.begin() + idx, 1);
+   }
+
+   // Create copy, attach new shape and strides and return
+   RTensor<Value_t, Container_t> x(*this);
+   x.fShape = shape;
+   x.fStrides = strides;
+   return x;
+}
+
+/// \brief Reshape tensor
+/// \param[in] shape Shape vector
+/// \returns New RTensor
+/// Reshape tensor without changing the overall size
+template <typename Value_t, typename Container_t>
+inline RTensor<Value_t, Container_t> RTensor<Value_t, Container_t>::Reshape(const Shape_t &shape)
+{
+   const auto size = Internal::GetSizeFromShape(shape);
+   if (size != fSize) {
+      std::stringstream ss;
+      ss << "Cannot reshape tensor with size " << fSize << " into shape { ";
+      for (std::size_t i = 0; i < shape.size(); i++) {
+         if (i != shape.size() - 1) {
+            ss << shape[i] << ", ";
+         } else {
+            ss << shape[i] << " }.";
+         }
+      }
+      throw std::runtime_error(ss.str());
+   }
+
+   // Compute new strides from shape
+   auto strides = Internal::ComputeStridesFromShape(shape, fLayout);
+
+   // Create copy, attach new shape and strides and return
+   RTensor<Value_t, Container_t> x(*this);
+   x.fShape = shape;
+   x.fStrides = strides;
+   return x;
+}
+
+/// \brief Create a slice of the tensor
+/// \param[in] slice Slice vector
+/// \returns New RTensor
+/// A slice is a subset of the tensor defined by a vector of pairs of indices.
+template <typename Value_t, typename Container_t>
+inline RTensor<Value_t, Container_t> RTensor<Value_t, Container_t>::Slice(const Slice_t &slice)
+{
+   // Sanitize size of slice
+   const auto sliceSize = slice.size();
+   const auto shapeSize = fShape.size();
+   if (sliceSize != shapeSize) {
+      std::stringstream ss;
+      ss << "Size of slice (" << sliceSize << ") is unequal number of dimensions (" << shapeSize << ").";
+      throw std::runtime_error(ss.str());
+   }
+
+   // Sanitize slice indices
+   // TODO: Sanitize slice indices
+   /*
+   for (std::size_t i = 0; i < sliceSize; i++) {
+   }
+   */
+
+   // Convert -1 in slice to proper pair of indices
+   // TODO
+
+   // Recompute shape and size
+   Shape_t shape(sliceSize);
+   for (std::size_t i = 0; i < sliceSize; i++) {
+      shape[i] = slice[i][1] - slice[i][0];
+   }
+   auto size = Internal::GetSizeFromShape(shape);
+
+   // Determine first element contributing to the slice and get the data pointer
+   Value_t *data;
+   Shape_t idx(sliceSize);
+   for (std::size_t i = 0; i < sliceSize; i++) {
+      idx[i] = slice[i][0];
+   }
+   data = &this->operator()(idx);
+
+   // Create copy and modify properties
+   RTensor<Value_t, Container_t> x(*this);
+   x.fData = data;
+   x.fShape = shape;
+   x.fSize = size;
+
+   // Squeeze tensor and return
+   return x.Squeeze();
+}
+
+/// \brief Pretty printing
+/// \param[in] os Output stream
+/// \param[in] x RTensor
+/// \return Modified output stream
+template <typename T>
+std::ostream &operator<<(std::ostream &os, RTensor<T> &x)
+{
+   const auto shapeSize = x.GetShape().size();
+   if (shapeSize == 1) {
+      os << "{ ";
+      const auto size = x.GetSize();
+      for (std::size_t i = 0; i < size; i++) {
+         os << x({i});
+         if (i != size - 1)
+            os << ", ";
+      }
+      os << " }";
+   } else if (shapeSize == 2) {
+      os << "{";
+      const auto shape = x.GetShape();
+      for (std::size_t i = 0; i < shape[0]; i++) {
+         os << " { ";
+         for (std::size_t j = 0; j < shape[1]; j++) {
+            os << x({i, j});
+            if (j < shape[1] - 1) {
+               os << ", ";
+            } else {
+               os << " ";
+            }
+         }
+         os << "}";
+      }
+      os << " }";
+   } else {
+      os << "{ printing not yet implemented for this rank }";
+   }
+   return os;
+}
+
+} // namespace TMVA::Experimental
+} // namespace TMVA
+
+namespace cling {
+template <typename T>
+std::string printValue(TMVA::Experimental::RTensor<T> *x)
+{
+   std::stringstream ss;
+   ss << *x;
+   return ss.str();
+}
+} // namespace cling
+
+#endif // TMVA_RTENSOR
diff --git a/tmva/tmva/test/CMakeLists.txt b/tmva/tmva/test/CMakeLists.txt
index 8b429cbbf5c..058a5b98909 100644
--- a/tmva/tmva/test/CMakeLists.txt
+++ b/tmva/tmva/test/CMakeLists.txt
@@ -10,6 +10,7 @@ include_directories(${ROOT_INCLUDE_DIRS})
 ROOT_ADD_GTEST(TestRandomGenerator
                TestRandomGenerator.cxx
                LIBRARIES ${Libraries})
+ROOT_ADD_GTEST(rtensor rtensor.cxx LIBRARIES ROOTVecOps TMVA)
 
 project(tmva-tests)
 find_package(ROOT REQUIRED)
diff --git a/tmva/tmva/test/rtensor.cxx b/tmva/tmva/test/rtensor.cxx
new file mode 100644
index 00000000000..432a029847b
--- /dev/null
+++ b/tmva/tmva/test/rtensor.cxx
@@ -0,0 +1,353 @@
+#include <gtest/gtest.h>
+#include <TMVA/RTensor.hxx>
+
+using namespace TMVA::Experimental;
+
+TEST(RTensor, GetElement)
+{
+   float data[6] = {0, 1, 2, 3, 4, 5};
+   RTensor<float> x(data, {2, 3});
+   auto shape = x.GetShape();
+   float count = 0.0;
+   for (size_t i = 0; i < shape[0]; i++) {
+      for (size_t j = 0; j < shape[1]; j++) {
+         EXPECT_EQ(count, x({i, j}));
+         EXPECT_EQ(count, x(i, j));
+         count++;
+      }
+   }
+}
+
+TEST(RTensor, SetElement)
+{
+   RTensor<float> x({2, 3});
+   auto shape = x.GetShape();
+   float count = 0.0;
+   for (size_t i = 0; i < shape[0]; i++) {
+      for (size_t j = 0; j < shape[1]; j++) {
+         x({i, j}) = count;
+         x(i, j) = count;
+         count++;
+      }
+   }
+   auto data = x.GetData();
+   for (size_t i = 0; i < shape.size(); i++)
+      EXPECT_EQ((float)i, *(data + i));
+}
+
+TEST(RTensor, AdoptMemory)
+{
+   float data[4] = {0, 0, 0, 0};
+   RTensor<float> x(data, {4});
+   for (size_t i = 0; i < 4; i++)
+      x(i) = (float)i;
+   for (size_t i = 0; i < 4; i++)
+      EXPECT_EQ((float)i, data[i]);
+}
+
+TEST(RTensor, GetShape)
+{
+   RTensor<int> x({2, 3});
+   const auto s = x.GetShape();
+   EXPECT_EQ(s.size(), 2u);
+   EXPECT_EQ(s[0], 2u);
+   EXPECT_EQ(s[1], 3u);
+}
+
+TEST(RTensor, Reshape)
+{
+   RTensor<int> x({2, 3});
+   const auto s = x.GetShape();
+   EXPECT_EQ(s.size(), 2u);
+   EXPECT_EQ(s[0], 2u);
+   EXPECT_EQ(s[1], 3u);
+
+   auto x2 = x.Reshape({1, 6});
+   const auto s2 = x2.GetShape();
+   EXPECT_EQ(s2.size(), 2u);
+   EXPECT_EQ(s2[0], 1u);
+   EXPECT_EQ(s2[1], 6u);
+
+   auto x3 = x.Reshape({6});
+   const auto s3 = x3.GetShape();
+   EXPECT_EQ(s3.size(), 1u);
+   EXPECT_EQ(s3[0], 6u);
+}
+
+TEST(RTensor, ExpandDims)
+{
+   RTensor<int> x({2, 3});
+   auto xa = x.ExpandDims(0);
+   const auto s = xa.GetShape();
+   EXPECT_EQ(s.size(), 3u);
+   EXPECT_EQ(s[0], 1u);
+   EXPECT_EQ(s[1], 2u);
+   EXPECT_EQ(s[2], 3u);
+
+   RTensor<int> x1({2, 3});
+   auto xb = x1.ExpandDims(1);
+   const auto s1 = xb.GetShape();
+   EXPECT_EQ(s1.size(), 3u);
+   EXPECT_EQ(s1[0], 2u);
+   EXPECT_EQ(s1[1], 1u);
+   EXPECT_EQ(s1[2], 3u);
+
+   RTensor<int> x2({2, 3});
+   auto xc = x2.ExpandDims(-1);
+   const auto s2 = xc.GetShape();
+   EXPECT_EQ(s2.size(), 3u);
+   EXPECT_EQ(s2[0], 2u);
+   EXPECT_EQ(s2[1], 3u);
+   EXPECT_EQ(s2[2], 1u);
+}
+
+TEST(RTensor, Squeeze)
+{
+   RTensor<int> x({1, 2, 3});
+   auto xa = x.Squeeze();
+   const auto s = xa.GetShape();
+   EXPECT_EQ(s.size(), 2u);
+   EXPECT_EQ(s[0], 2u);
+   EXPECT_EQ(s[1], 3u);
+
+   RTensor<int> x1({1, 2, 1, 3, 1});
+   auto xb = x1.Squeeze();
+   const auto s1 = xb.GetShape();
+   EXPECT_EQ(s1.size(), 2u);
+   EXPECT_EQ(s1[0], 2u);
+   EXPECT_EQ(s1[1], 3u);
+
+   RTensor<int> x2({1, 1, 1});
+   auto xc = x2.Squeeze();
+   const auto s2 = xc.GetShape();
+   EXPECT_EQ(s2.size(), 1u);
+   EXPECT_EQ(s2[0], 1u);
+
+   RTensor<int> x3({});
+   auto xd = x3.Squeeze();
+   const auto s3 = xd.GetShape();
+   EXPECT_EQ(s3.size(), 0u);
+}
+
+TEST(RTensor, Transpose)
+{
+   // Layout (physical):
+   // 0, 1, 2, 3, 4, 5
+   float data[6] = {0, 1, 2, 3, 4, 5};
+   RTensor<float> x(data, {2, 3});
+
+   // Layout (logical):
+   // 0, 1, 2
+   // 3, 4, 5
+   float count = 0;
+   for (size_t i = 0; i < 2; i++) {
+      for (size_t j = 0; j < 3; j++) {
+         x(i, j) = count;
+         count++;
+      }
+   }
+
+   // Layout (logical):
+   // 0, 3
+   // 1, 4
+   // 2, 5
+   auto x2 = x.Transpose();
+   EXPECT_EQ(x2.GetShape().size(), 2u);
+   EXPECT_EQ(x2.GetShape()[0], 3u);
+   EXPECT_EQ(x2.GetShape()[1], 2u);
+   count = 0;
+   for (size_t i = 0; i < 2; i++) {
+      for (size_t j = 0; j < 3; j++) {
+         x2(j, i) = count;
+         count++;
+      }
+   }
+}
+
+TEST(RTensor, InitWithZeros)
+{
+   RTensor<float> x({2, 3});
+   for (size_t i = 0; i < 2; i++) {
+      for (size_t j = 0; j < 3; j++) {
+         EXPECT_EQ(0.0, x(i, j));
+      }
+   }
+}
+
+TEST(RTensor, SetElementRowMajor)
+{
+   // Layout (logical):
+   // 0, 1, 2
+   // 3, 4, 5
+   RTensor<float> x({2, 3}, MemoryLayout::RowMajor);
+   float count = 0;
+   for (size_t i = 0; i < 2; i++) {
+      for (size_t j = 0; j < 3; j++) {
+         x(i, j) = count;
+         count++;
+      }
+   }
+
+   // Layout (physical):
+   // 0, 1, 2, 3, 4, 5
+   float ref[6] = {0, 1, 2, 3, 4, 5};
+   float *data = x.GetData();
+   for (size_t i = 0; i < 6; i++) {
+      EXPECT_EQ(data[i], ref[i]);
+   }
+}
+
+TEST(RTensor, SetElementColumnMajor)
+{
+   // Layout (logical):
+   // 0, 1, 2
+   // 3, 4, 5
+   RTensor<float> x({2, 3}, MemoryLayout::ColumnMajor);
+   float count = 0;
+   for (size_t i = 0; i < 2; i++) {
+      for (size_t j = 0; j < 3; j++) {
+         x(i, j) = count;
+         count++;
+      }
+   }
+
+   // Layout (physical):
+   // 0, 3, 1, 4, 2, 5
+   float ref[6] = {0, 3, 1, 4, 2, 5};
+   float *data = x.GetData();
+   for (size_t i = 0; i < 6; i++) {
+      EXPECT_EQ(data[i], ref[i]);
+   }
+}
+
+TEST(RTensor, GetElementRowMajor)
+{
+   // Layout (physical):
+   // 0, 1, 2, 3, 4, 5
+   float data[6] = {0, 1, 2, 3, 4, 5};
+   RTensor<float> x(data, {2, 3}, MemoryLayout::RowMajor);
+
+   // Layout (logical):
+   // 0, 1, 2
+   // 3, 4, 5
+   float count = 0;
+   for (size_t i = 0; i < 2; i++) {
+      for (size_t j = 0; j < 3; j++) {
+         EXPECT_EQ(x(i, j), count);
+         count++;
+      }
+   }
+}
+
+TEST(RTensor, GetElementColumnMajor)
+{
+   // Layout (physical):
+   // 0, 3, 1, 4, 2, 5
+   float data[6] = {0, 3, 1, 4, 2, 5};
+   RTensor<float> x(data, {2, 3}, MemoryLayout::ColumnMajor);
+
+   // Layout (logical):
+   // 0, 1, 2
+   // 3, 4, 5
+   float count = 0;
+   for (size_t i = 0; i < 2; i++) {
+      for (size_t j = 0; j < 3; j++) {
+         EXPECT_EQ(x(i, j), count);
+         count++;
+      }
+   }
+}
+
+TEST(RTensor, SliceRowMajor)
+{
+   // Data layout:
+   // [ 0, 1, 2, 3, 4, 5 ] ]
+   RTensor<float> x({2, 3}, MemoryLayout::RowMajor);
+   float c = 0.f;
+   for (auto i = 0; i < 2; i++) {
+      for (auto j = 0; j < 3; j++) {
+         x(i, j) = c;
+         c++;
+      }
+   }
+
+   // Slice:
+   // [ [ 0, 1 ], [ 3, 4 ] ]
+   auto s1 = x.Slice({{0, 2}, {0, 2}});
+   EXPECT_EQ(s1.GetSize(), 4u);
+   EXPECT_EQ(s1.GetShape().size(), 2u);
+   EXPECT_EQ(s1.GetShape()[0], 2u);
+   EXPECT_EQ(s1.GetShape()[1], 2u);
+   EXPECT_EQ(s1(0, 0), 0.f);
+   EXPECT_EQ(s1(0, 1), 1.f);
+   EXPECT_EQ(s1(1, 0), 3.f);
+   EXPECT_EQ(s1(1, 1), 4.f);
+
+   // Slice:
+   // [ 5 ]
+   auto s2 = x.Slice({{1, 2}, {2, 3}});
+   EXPECT_EQ(s2.GetSize(), 1u);
+   EXPECT_EQ(s2.GetShape().size(), 1u);
+   EXPECT_EQ(s2.GetShape()[0], 1u);
+   EXPECT_EQ(s2(0), 5.f);
+
+   // Slice:
+   // [ [ 0, 1, 2 ], [ 3, 4, 5 ] ]
+   auto s3 = x.Slice({{0, 2}, {0, 3}});
+   EXPECT_EQ(s3.GetSize(), 6u);
+   EXPECT_EQ(s3.GetShape().size(), 2u);
+   EXPECT_EQ(s3.GetShape()[0], 2u);
+   EXPECT_EQ(s3.GetShape()[1], 3u);
+   for(auto i = 0; i < 3; i ++) {
+      for(auto j = 0; j < 2; j ++) {
+         EXPECT_EQ(x(i, j), s3(i, j));
+      }
+   }
+}
+
+TEST(RTensor, SliceColumnMajor)
+{
+   // Data layout:
+   // [ 0, 3, 1, 4, 2, 5  ]
+   RTensor<float> x({2, 3}, MemoryLayout::ColumnMajor);
+   float c = 0.f;
+   for (auto i = 0; i < 2; i++) {
+      for (auto j = 0; j < 3; j++) {
+         x(i, j) = c;
+         c++;
+      }
+   }
+
+   // Slice:
+   // [ [ 0, 1 ], [ 3, 4 ] ]
+   auto s1 = x.Slice({{0, 2}, {0, 2}});
+   EXPECT_EQ(s1.GetSize(), 4u);
+   EXPECT_EQ(s1.GetShape().size(), 2u);
+   EXPECT_EQ(s1.GetShape()[0], 2u);
+   EXPECT_EQ(s1.GetShape()[1], 2u);
+   EXPECT_EQ(s1(0, 0), 0.f);
+   EXPECT_EQ(s1(0, 1), 1.f);
+   EXPECT_EQ(s1(1, 0), 3.f);
+   EXPECT_EQ(s1(1, 1), 4.f);
+
+   // Slice:
+   // [ 5 ]
+   auto s2 = x.Slice({{1, 2}, {2, 3}});
+   EXPECT_EQ(s2.GetSize(), 1u);
+   EXPECT_EQ(s2.GetShape().size(), 1u);
+   EXPECT_EQ(s2.GetShape()[0], 1u);
+   EXPECT_EQ(s2(0), 5.f);
+
+   // Slice:
+   // [ [ 0, 1, 2 ], [ 3, 4, 5 ] ]
+   auto s3 = x.Slice({{0, 2}, {0, 3}});
+   EXPECT_EQ(s3.GetSize(), 6u);
+   EXPECT_EQ(s3.GetShape().size(), 2u);
+   EXPECT_EQ(s3.GetShape()[0], 2u);
+   EXPECT_EQ(s3.GetShape()[1], 3u);
+   for(auto i = 0; i < 3; i ++) {
+      for(auto j = 0; j < 2; j ++) {
+         EXPECT_EQ(x(i, j), s3(i, j));
+      }
+   }
+}
diff --git a/tutorials/tmva/tmva001_RTensor.C b/tutorials/tmva/tmva001_RTensor.C
new file mode 100644
index 00000000000..c8c764113ef
--- /dev/null
+++ b/tutorials/tmva/tmva001_RTensor.C
@@ -0,0 +1,54 @@
+/// \file
+/// \ingroup tutorial_dataframe
+/// \notebook -nodraw
+/// This tutorial illustrates the basic features of the RTensor class,
+/// RTensor is a std::vector-like container with additional shape information.
+/// The class serves as an interface in C++ between multi-dimensional data and
+/// the algorithm such as in machine learning workflows. The interface is similar
+/// to Numpy arrays and provides a subset of the functionality.
+///
+/// \macro_code
+/// \macro_output
+///
+/// \date December 2018
+/// \author Stefan Wunsch
+
+using namespace TMVA::Experimental;
+using namespace ROOT::VecOps;
+
+void tmva001_RTensor()
+{
+   // Create RTensor from scratch
+   RTensor<float> x({2, 2});
+   cout << x << endl;
+
+   // Assign some data
+   x({0, 0}) = 1;
+   x({0, 1}) = 2;
+   x({1, 0}) = 3;
+   x({1, 1}) = 4;
+
+   // Apply transformations
+   auto x2 = x.Reshape({1, 4}).Squeeze();
+   cout << x2 << endl;
+
+   // Slice
+   auto x3 = x.Reshape({2, 2}).Slice({{0, 2}, {0, 1}});
+   cout << x3 << endl;
+
+   // Create tensor as view on data without ownership
+   float data[] = {5, 6, 7, 8};
+   RTensor<float> y(data, {2, 2});
+   cout << y << endl;
+
+   // Create tensor as view on data with ownership
+   auto data2 = std::make_shared<std::vector<float>>(4);
+   float c = 9;
+   for (auto &v : *data2) {
+      v = c;
+      c++;
+   }
+
+   RTensor<float> z(data2, {2, 2});
+   cout << z << endl;
+}
-- 
GitLab