added a mapper that maps the indexes of a view to the indices of a container

include/xtensor/views/index_mapper.hpp

@@ -0,0 +1,236 @@
+/***************************************************************************
+ * Copyright (c) Johan Mabille, Sylvain Corlay and Wolf Vollprecht          *
+ * Copyright (c) QuantStack                                                 *
+ *                                                                          *
+ * Distributed under the terms of the BSD 3-Clause License.                 *
+ *                                                                          *
+ * The full license is in the file LICENSE, distributed with this software. *
+ ****************************************************************************/
+
+#ifndef XTENSOR_INDEX_MAPPER_HPP
+#define XTENSOR_INDEX_MAPPER_HPP
+
+#include "xview.hpp"
+
+namespace xt
+{
+
+template<class UndefinedView> struct index_mapper; 
+
+/**
+ * @class index_mapper
+ * @brief A helper class for mapping indices between views and their underlying containers.
+ * 
+ * The `index_mapper` class provides functionality to convert indices from a view's coordinate system
+ * to the corresponding indices in the underlying container. This is particularly useful for views
+ * that contain integral slices (fixed indices), as these slices reduce the dimensionality of the view.
+ * 
+ * @tparam UndefinedView The primary template parameter, specialized for `xt::xview` types.
+ * 
+ * @note This class is specialized for `xt::xview<UnderlyingContainer, Slices...>` types only.
+ *       Other view types will trigger a compilation error.
+ * 
+ * @example
+ * @code
+ * xt::xarray<double> a = xt::arange(24).reshape({2, 3, 4});
+ * auto view1 = xt::view(a, 1, xt::all(), xt::all());  // Fixed first dimension
+ * index_mapper<decltype(view1)> mapper;
+ * 
+ * // Map view indices (i,j) to container indices (1,i,j)
+ * double val = mapper.map(a, view1, 0, 0);  // Returns a(1, 0, 0)
+ * double val2 = mapper.map(a, view1, 1, 2); // Returns a(1, 1, 2)
+ * @endcode
+ */
+template<class UnderlyingContainer, class... Slices>
+class index_mapper< xt::xview<UnderlyingContainer, Slices...> >
+{
+    static constexpr size_t n_slices = sizeof...(Slices);                      ///<  @brief Total number of slices in the view
+    static constexpr size_t n_free   = ((!std::is_integral_v<Slices>) + ... ); ///<  @brief Number of free (non-integral) slices
+public:
+    using view_type = xt::xview<UnderlyingContainer, Slices...>;  ///<   @brief The view type this mapper works with
+
+    using value_type = typename xt::xview<UnderlyingContainer, Slices...>::value_type; ///<  @brief Value type of the underlying container
+private:
+	/// @brief Helper type alias for the I-th slice type
+    template<size_t I> using ith_slice_type = std::tuple_element_t<I, std::tuple<Slices...> >;
+
+	/**
+     * @brief Helper metafunction to generate an index sequence excluding newaxis slices.
+     * 
+     * This recursive template builds an `std::index_sequence` containing indices of slices
+     * that are not `xt::newaxis`. Newaxis slices increase dimensionality but don't correspond
+     * to actual dimensions in the underlying container.
+     * 
+     * @tparam first Current slice index being processed.
+     * @tparam indices... Accumulated indices of non-newaxis slices.
+     */
+    template<size_t first, size_t... indices>
+    struct indices_sequence_helper
+    {
+        using not_new_axis_type = typename indices_sequence_helper<first + 1, indices..., first>::Type; // we add the current axis
+        using new_axis_type     = typename indices_sequence_helper<first + 1, indices...>::Type;        // we skip the current axis
+
+        using Type = std::conditional_t< xt::detail::is_newaxis< ith_slice_type<first> >::value , new_axis_type, not_new_axis_type>;
+    };
+    /// @brief Base case: recursion termination
+    template<size_t... indices>
+    struct indices_sequence_helper<n_slices, indices...>
+    {
+        using Type = std::index_sequence<indices...>;
+    };
+
+    using indices_sequence  = indices_sequence_helper<0>::Type; ///<  @brief Index sequence of non-newaxis slices
+
+	/**
+     * @brief Maps an index for a specific slice to the corresponding index in the underlying container.
+     * 
+     * For integral slices (fixed indices), returns the fixed index value.
+     * For non-integral slices (like `xt::all()`), applies the slice transformation to the index.
+     * 
+     * @tparam I The slice index to map.
+     * @tparam Index Type of the index (must be integral).
+     * @param view The view object containing slice information.
+     * @param i The index within the slice to map.
+     * @return size_t The mapped index in the underlying container.
+     * 
+     * @throws Assertion failure if `i != 0` for integral slices.
+     * @throws Assertion failure if `i >= slice.size()` for non-integral slices.
+     */
+    template<size_t I, std::integral Index>
+    size_t map_ith_index(const view_type& view, const Index i) const;
+
+	/**
+     * @brief Maps all indices and accesses the container.
+     * 
+     * @tparam Is Index sequence for parameter pack expansion.
+     * @param container The underlying container to access.
+     * @param view The view providing slice information.
+     * @param indices Array of indices for all slices.
+     * @return value_type The value at the mapped location in the container.
+     */
+    template<size_t... Is>
+    value_type map_all_indices(const UnderlyingContainer& container, const view_type& view, std::index_sequence<Is...>, const std::array<size_t, n_slices>& indices) const;
+
+	/**
+     * @brief Maps all indices and accesses the container with bounds checking.
+     * 
+     * Same as `map_all_indices` but uses `container.at()` which performs bounds checking.
+     * 
+     * @tparam Is Index sequence for parameter pack expansion.
+     * @param container The underlying container to access.
+     * @param view The view providing slice information.
+     * @param indices Array of indices for all slices.
+     * @return value_type The value at the mapped location in the container.
+     * 
+     * @throws std::out_of_range if any index is out of bounds.
+     */
+    template<size_t... Is>
+    value_type map_at_all_indices(const UnderlyingContainer& container, const view_type& view, std::index_sequence<Is...>, const std::array<size_t, n_slices>& indices) const;
+public:
+	/**
+     * @brief Maps view indices to container indices and returns the value.
+     * 
+     * Converts the provided indices (for the free dimensions of the view) to
+     * the corresponding indices in the underlying container and returns the value.
+     * 
+     * @tparam Indices Types of the indices (must be integral).
+     * @param container The underlying container to access.
+     * @param view The view providing slice information.
+     * @param indices The indices for the free dimensions of the view.
+     * @return value_type The value at the mapped location in the container.
+     * 
+     * @note The number of provided indices must equal `n_free` (number of non-integral slices).
+     * 
+     * @example
+     * @code
+     * // For view(a, 1, all(), all()):
+     * // n_free = 2 (two all() slices)
+     * mapper.map(a, view, i, j);  // Maps to a(1, i, j)
+     * @endcode
+     */
+    template<std::integral... Indices> 
+    value_type map(const UnderlyingContainer& container, const view_type& view, const Indices... indices) const 
+        requires(sizeof...(Indices) == n_free);
+
+	/**
+     * @brief Maps view indices to container indices with bounds checking.
+     * 
+     * Same as `map()` but uses bounds-checked access via `container.at()`.
+     * 
+     * @tparam Indices Types of the indices (must be integral).
+     * @param container The underlying container to access.
+     * @param view The view providing slice information.
+     * @param indices The indices for the free dimensions of the view.
+     * @return value_type The value at the mapped location in the container.
+     * 
+     * @throws std::out_of_range if any mapped index is out of bounds.
+     */
+    template<std::integral... Indices> 
+    value_type map_at(const UnderlyingContainer& container, const view_type& view, const Indices... indices) const 
+        requires(sizeof...(Indices) == n_free);
+
+    constexpr size_t dimension() const { return n_free; }
+};
+
+/*******************************
+ * index_mapper implementation *
+ *******************************/
+
+template<class UnderlyingContainer, class... Slices> template<std::integral... Indices> 
+auto index_mapper< xt::xview<UnderlyingContainer, Slices...> >::map(
+    const UnderlyingContainer& container, 
+    const view_type&           view, 
+    const Indices...           indices) const -> value_type 
+    requires(sizeof...(Indices) == n_free)
+{
+    std::array<size_t, sizeof...(indices)> args{ size_t(indices)...};
+
+	auto it = std::cbegin(args);
+	std::array<size_t, n_slices> args_full{ (std::is_integral_v<Slices> ? size_t(0) : *it++)... };
+
+    return map_all_indices(container, view, indices_sequence{}, args_full);
+}
+
+template<class UnderlyingContainer, class... Slices> template<std::integral... Indices> 
+auto index_mapper< xt::xview<UnderlyingContainer, Slices...> >::map_at(
+	const UnderlyingContainer& container, 
+    const view_type&           view, 
+    const Indices...           indices) const -> value_type 
+    requires(sizeof...(Indices) == n_free)
+{
+    std::array<size_t, sizeof...(indices)> args{ size_t(indices)...};
+
+    auto it = std::cbegin(args);
+    std::array<size_t, n_slices> args_full{ (std::is_integral_v<Slices> ? size_t(0) : *it++)... };
+
+    return map_at_all_indices(container, view, indices_sequence{}, args_full);
+}
+
+template<class UnderlyingContainer, class... Slices> template<size_t... Is>
+auto index_mapper< xt::xview<UnderlyingContainer, Slices...> >::map_all_indices(const UnderlyingContainer& container, const view_type& view, std::index_sequence<Is...>, const std::array<size_t, n_slices>& indices) const -> value_type
+{		
+	return container(map_ith_index<Is>(view, indices[Is])...);
+}
+
+template<class UnderlyingContainer, class... Slices> template<size_t... Is>
+auto index_mapper< xt::xview<UnderlyingContainer, Slices...> >::map_at_all_indices(const UnderlyingContainer& container, const view_type& view, std::index_sequence<Is...>, const std::array<size_t, n_slices>& indices) const -> value_type
+{
+	return container.at(map_ith_index<Is>(view, indices[Is])...);
+}
+
+template<class UnderlyingContainer, class... Slices> template<size_t I, std::integral Index>
+auto index_mapper< xt::xview<UnderlyingContainer, Slices...> >::map_ith_index(const view_type& view, const Index i) const -> size_t
+{
+    using current_slice = std::tuple_element_t<I, std::tuple<Slices...>>;
+
+    static_assert(not xt::detail::is_newaxis<current_slice>::value);
+
+    const auto& slice = std::get<I>(view.slices());
+
+    if constexpr (std::is_integral_v<current_slice>) { assert(i == 0);           return size_t(slice);    }
+    else                                             { assert(i < slice.size()); return size_t(slice(i)); }
+}
+
+} // namespace xt
+
+#endif // XTENSOR_INDEX_MAPPER_HPP

test/test_xview.cpp

@@ -30,6 +30,7 @@
 #include "xtensor/misc/xmanipulation.hpp"
 #include "xtensor/views/xstrided_view.hpp"
 #include "xtensor/views/xview.hpp"
+#include "xtensor/views/index_mapper.hpp"
 
 namespace xt
 {
@@ -143,6 +144,57 @@ namespace xt
         }
     }
 
+    TEST(xview_mapping, simple)
+    {
+        view_shape_type shape = {3, 4};
+        xarray<double> a(shape);
+        std::vector<double> data = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
+        std::copy(data.cbegin(), data.cend(), a.template begin<layout_type::row_major>());
+
+        auto view1 = view(a, 1, range(1, 4));
+
+        index_mapper<decltype(view1)> mapper1;
+
+        EXPECT_EQ(a(1, 1), mapper1.map(a, view1, 0));
+        EXPECT_EQ(a(1, 2), mapper1.map(a, view1, 1));
+        EXPECT_EQ(size_t(1), mapper1.dimension());
+        //~ XT_EXPECT_ANY_THROW(mapper1.map_at(a, view1, 10));
+
+        auto view0 = view(a, 0, range(0, 3));
+        index_mapper<decltype(view0)> mapper0;
+
+        EXPECT_EQ(a(0, 0), mapper0.map(a, view0, 0));
+        EXPECT_EQ(a(0, 1), mapper0.map(a, view0, 1));
+        EXPECT_EQ(size_t(1), mapper0.dimension());
+
+        auto view2 = view(a, range(0, 2), 2);
+        index_mapper<decltype(view2)> mapper2;
+        EXPECT_EQ(a(0, 2), mapper2.map(a, view2, 0));
+        EXPECT_EQ(a(1, 2), mapper2.map(a, view2, 1));
+        EXPECT_EQ(size_t(1), mapper2.dimension());
+
+        //~ auto view4 = view(a, 1);
+        //~ index_mapper<decltype(view4)> mapper4;
+        //~ EXPECT_EQ(size_t(1), mapper4.dimension());
+        //~ 
+        //~ auto view5 = view(view4, 1);
+        //~ index_mapper<decltype(view5)> mapper5;
+        //~ EXPECT_EQ(size_t(0), mapper5.dimension());
+
+        auto view6 = view(a, 1, all());
+        index_mapper<decltype(view6)> mapper6;
+        EXPECT_EQ(a(1, 0), mapper6.map(a, view6, 0));
+        EXPECT_EQ(a(1, 1), mapper6.map(a, view6, 1));
+        EXPECT_EQ(a(1, 2), mapper6.map(a, view6, 2));
+        EXPECT_EQ(a(1, 3), mapper6.map(a, view6, 3));
+
+        auto view7 = view(a, all(), 2);
+        index_mapper<decltype(view7)> mapper7;
+        EXPECT_EQ(a(0, 2), mapper7.map(a, view7, 0));
+        EXPECT_EQ(a(1, 2), mapper7.map(a, view7, 1));
+        EXPECT_EQ(a(2, 2), mapper7.map(a, view7, 2));
+    }
+
     TEST(xview, negative_index)
     {
         view_shape_type shape = {3, 4};
@@ -269,6 +321,29 @@ namespace xt
         EXPECT_EQ(a(1, 1, 1), view1.element(idx.cbegin(), idx.cend()));
     }
 
+    TEST(xview_mapping, three_dimensional)
+    {
+        view_shape_type shape = {3, 4, 2};
+        std::vector<double> data = {1,  2,  3,  4,  5,  6,   7,   8,
+
+                                    9,  10, 11, 12, 21, 22,  23,  24,
+
+                                    25, 26, 27, 28, 29, 210, 211, 212};
+        xarray<double> a(shape);
+        std::copy(data.cbegin(), data.cend(), a.template begin<layout_type::row_major>());
+
+        auto view1 = view(a, 1, all(), all());
+        index_mapper<decltype(view1)> mapper1;
+
+        EXPECT_EQ(size_t(2), mapper1.dimension());
+        std::cout << "===================================" << std::endl;
+        EXPECT_EQ(a(1, 0, 0), mapper1.map(a, view1, 0, 0));
+        EXPECT_EQ(a(1, 0, 1), mapper1.map(a, view1, 0, 1));
+        //~ EXPECT_EQ(a(1, 1, 0), mapper1.map(a, view1, 1, 0));
+        //~ EXPECT_EQ(a(1, 1, 1), mapper1.map(a, view1, 1, 1));
+        //~ XT_EXPECT_ANY_THROW(mapper1.map_at(a, view1, 10, 10));
+    }
+
     TEST(xview, integral_count)
     {
         size_t squeeze1 = integral_count<size_t, size_t, size_t, xrange<size_t>>();