diff --git a/.gitignore b/.gitignore index 936f0405..2df4161c 100644 --- a/.gitignore +++ b/.gitignore @@ -20,6 +20,7 @@ # These are directories used by IDEs for storing settings .idea/ .vscode/ +.cache/ # These are common Python virtual enviornment directory names venv/ diff --git a/docs/source/developer/design/assets/shape_design_v2.png b/docs/source/developer/design/assets/shape_design_v2.png new file mode 100644 index 00000000..d8edb4ce Binary files /dev/null and b/docs/source/developer/design/assets/shape_design_v2.png differ diff --git a/docs/source/developer/design/index.rst b/docs/source/developer/design/index.rst index 33cda52a..81b3540c 100644 --- a/docs/source/developer/design/index.rst +++ b/docs/source/developer/design/index.rst @@ -33,3 +33,9 @@ Design of TensorWrapper expression op_graph sparse_maps + +.. toctree:: + :maxdepth: 2 + :caption: Old Designs + + shapev1 diff --git a/docs/source/developer/design/shape.rst b/docs/source/developer/design/shape.rst index ef319886..3804121a 100644 --- a/docs/source/developer/design/shape.rst +++ b/docs/source/developer/design/shape.rst @@ -153,603 +153,33 @@ Memory allocation. Shape Design ************ -In designing the class hierarchy we note the following: +TODO: Express better. -- An algorithm which works for a jagged shape should work for a smooth shape - as well. The reverse, smooth algorithms with jagged shapes, will in general - not work. -- Nestings are logically imposed over an existing shape. The resulting nested - shape is still an instance of the underlying shape. -- Tiled shapes are a subcategory of jagged shapes. +.. note:: -.. _fig_shape_design: + This is a proposal. It has not been fully implemented yet. -.. figure:: assets/shape.png - :align: center - - The architecture of TensorWrapper's Shape component. - -:numref:`fig_shape_design` shows the classes primarily responsible for -implementing the shape component. Most end users will deal with the ``Shape`` -class. - -ShapeBase -========= - -The unifying features of all shapes were summarized in the -:ref:`shape_rank_and_extents` consideration. ``ShapeBase`` provides the API -that all shapes must minimally satisfy because they are shapes. The actual -class serves primarily as code factorization. - -Shape -===== - -The ``Shape`` class describes a (smooth) hyper-rectangular array of data and -can be used for "traditional" tensors (those which are not nested or jagged). -Most end users will simply create ``Shape`` objects and pass them on to -``TensorWrapper``. We expect that manipulations of ``Shape`` objects will be -aimed at TensorWrapper developers. - -JaggedShape -=========== - -To satisfy the :ref:`shape_jagged` consideration we introduce ``JaggedShape``. -``JaggedShape`` is similar to ``Shape`` except that users must explicitly -provide the shape of the slices. Generally speaking ``JaggedShape`` objects of -rank |r| will contain a series of rank :math:`s` ``Shape`` objects. The -actual ``JaggedShape`` object serves as a map from an index with :math:`(r-s)` -indices to the ``Shape`` of that slice. Like ``Shape`` we expect users to -primarily be concerned with construction. Again, manipulations of the -``JaggedShape`` will be primarily of interest to TensorWrapper developers. - -TiledShape -========== - -Introduced primarily as a convenience for constructing ``JaggedShape`` objects -by tiling. - -Nested -====== - -To address the :ref:`shape_nested` consideration, we have added a ``Nested`` -class which is templated on the type of the tensor's overall shape. - -With objects like ``Shape`` and ``JaggedShape``, TensorWrapper can't tell how -the user is thinking of the tensor. For example, they could be thinking of a -matrix as a matrix or as a vector of vectors. The point of the ``Nested`` -object is to partition the ranks of the tensor into layers, so we know how many -layers the user is viewing the tensor as, and how many ranks each layer has. -Mathematically the various ways of a viewing a tensor do not change the -properties of the tensor; however, when we are physically laying the tensor out -on the computer, how we view the tensor can affect physical layout. - -IndexedShape -============ - -Consideration :ref:`shape_combining_shapes` requires us to be able to compose -the various shape objects. To do this, we rely on the same mechanism used for -``TensorWrapper``, *i.e.*, an expression layer. More specifically, -``IndexedShape`` objects result from indexing a shape like ``s("i,j,k")``. While -``IndexedShape`` is technically exposed to the user, user can be somewhat -oblivious to its existence. - -******************* -Proposed Shape APIs -******************* - -Constructing a ``Shape`` -======================== - -Creating a non-nested shape just requires knowing the extent of each mode: - -.. code-block:: c++ - - Shape null_shape; // No rank and no elements - Shape rank0_shape{}; // A scalar - Shape rank1_shape{10}; // 10 element vector - Shape rank2_shape{10, 20}; // 10 by 20 matrix - Shape rank3_shape{10, 20, 30}; // 10 by 20 by 30 tensor - -Note that following usual C++ rules the first two lines actually call -different constructors (default ctor vs. initializer list). Using an initializer -list requires us to know the rank at compile time. If we want to determine the -rank at runtime we can use iterators: - -.. code-block:: c++ - - // Somehow create a vector of extents - using size_type = Shape::size_type; - std::vector extents = get_extents(); - - // Construct Shape from iterator pair - Shape runtime_rank_shape(extents.begin(), extents.end()); - -Jagged Shape Construction -========================= - -For a ``Shape`` we need to specify the extents of each mode. ``JaggedShape`` -declaration is done in terms of ``Shape`` objects and looks like: - -.. code-block:: c++ - - // For brevity define variables - Shape s10{10}, s20{20}, s30{30}; - Shape s10_20{10, 20}, s30_40{30, 40}, s50_60{50, 60}; - Shape s10_20_30{10, 20, 30}, s40_50_60{40, 50, 60}; - - // No elements, no rank - JaggedShape null_shape; - - // Smooth scalar viewed as a JaggedShape (note () not {}) - JaggedShape smooth0_shape(Shape{}); - - // Smooth vector viewed as a JaggedShape (note () not {}) - JaggedShape smooth1_shape(Shape{10}); - - // Smooth matrix viewed as a JaggedShape (note () not {}) - JaggedShape smooth2_shape(Shape{10, 20}); - - // A "jagged" vector with no elements - JaggedShape rank0_shape{}; - - // A jagged matrix with 1 row, note the {} - JaggedShape rank2_shape{s10}; - - // A jagged matrix with 3 rows; row 0 has 10 elements, row 1 has 20, row 2 30 - JaggedShape rank2_shape2{s10, s20, s30}; - - // A jagged rank 3 tensor with smooth matrices. Matrix 0 is 10 by 20, - // matrix 1 is 30 by 40, and matrix 2 is 50 by 60 - JaggedShape rank3_shape{s10_20, s30_40, s50_60}; - - // A jagged rank 3 tensor where elements are jagged matrices. Matrix 0 is - // 1 by 10, matrix 2 has 20 columns in row 0 and 30 columns in row 2, and - // matrix 3 has 30 columns in row 0, 10 columns in row 1, and 20 columns in - // row 2 - JaggedShape rank3_shape2{{s10}, - {s20, s30}, - {s30, s10, s20}}; - - // A jagged rank 4 tensor where the 0-th element of the 0-th mode is a - // 10 by 20 by 30 smooth tensor and the 1-st element is a 40 by 50 by 60 - // smooth tensor - JaggedShape rank4_shape{s10_20_30, s40_50_60}; - - // A jagged rank 4 tensor where the elements are jagged rank 3 tensors. - // Taking slices along the 0 and 1-st modes, the (0,0)-th slice is a 10 by 20 - // matrix, the (0,1)-th slice is a 30 by 40 matrix, the (1,0)-th slice is - // a 30 by 40 matrix, the (1,1)-th slice is a 10 by 20 matrix, and the - // (1,2)-th slice is a 50 by 60 matrix - JaggedShape rank4_shape2{{s10_20, s30_40}, - {s30_40, s10_20, s50_60}}; - - // A jagged rank 4 tensors where the elements are jagged rank 3 tensors, - // which have jagged matrices for elements. Taking slices along the 0, 1, and - // 2 modes we have: - // - (0,0,0) is a 10 element vector, - // - (0,1,0) is a 20 element vector, - // - (0,1,1) is a 30 element vector, - // - (1,0,0) is a 10 element vector, - // - (1,0,1) is a 30 element vector, - // - (1,1,0) is a 20 element vector, - // - (1,2,0) is a 10 element vector, - // - (1,2,1) is a 20 element vector, - // - (1,2,2) is a 30 element vector - JaggedShape rank4_shape{{{s10}, {s20, s30}}, - {{s10, s30}, {s20}, {s10, s20, s30}}}; - - -Consider the shape of the (0,1) slice of ``rank4_shape``. This slice is a -vector of vectors where the outer vector has two elements, element 0 of the -outer vector is a 10-element vector and element 1 is a 30-element vector. In -other words the shape of the (0,1) slice of ``rank4_shape`` describes a jagged -matrix, which could have been initialized by ``JaggedShape{s20, s30}``. In turn -the above construction of ``rank4_shape`` is actually equivalent to: - -.. code-block:: c++ - - JaggedShape e00{s10}; - JaggedShape e01{s20, s30}; - JaggedShape e10{s10, s30}; - JaggedShape e11{s20}; - JaggedShape e12{s10, s20, s30}; - JaggedShape e0{e00, e01}; - JaggedShape e1{e10, e11, e12}; - JaggedShape rank4_shape{e0, e1}; - -And we see that ``JaggedShape`` is a recursive structure and thus the runtime -mechanism for initializing a ``JaggedShape`` is with iterators running over -``JaggedShape`` objects: - -.. code-block:: c++ - - std::vector slice_shapes = get_slices(); - JaggedShape shape(slice_shapes.begin(), slice_shapes.end()); - -So far we have focused on the most general way to create a ``JaggedShape`` one -of the most common ways to form a ``JaggedShape`` is by tiling. Consider a -30 by 30 matrix where we tile each mode into 5, 15, and 10 element chunks. -Using ``JaggedShape`` this can be done by: - -.. code-block:: c++ - - JaggedShape js{{Shape{5,5}, Shape{5, 15}, Shape{5,10}}, - {Shape{15,5}, Shape{15,15}, Shape{15,10}}, - {Shape{10,5}, Shape{10,15}, Shape{10,10}}}; - -This is an admittedly verbose declaration. Thus for the special case of crating -``JaggedShape`` objects which result from tiling smooth ``Shape`` objects we -introduce the ``TiledShape`` class. Using ``TiledShape`` the same shape could -be declared via: - -.. code-block:: c++ - - TiledShape s{{5, 10, 15}, {5, 10, 15}}; - - -Constructing Nested Shapes -========================== - -Creating a ``Nested`` object requires knowing the shape of the tensor and how -the indices are partitioned into layers. - -.. code-block:: c++ - - // Zero layer scalar - Nested s({}, Shape{}); - - // One layer scalar - Nested s0({0}, Shape{}); - - // Two layer scalar - Nested s0_0({0, 0}, Shape{}); - - // One layer vector - Nested s1({1}, Shape{10}); - - // Two layer vector (mode in layer 0) - Nested s1_0({1, 0}, Shape{10}); - - // Two layer vector (mode in layer 1) - Nested s0_1({0, 1}, Shape{10}); - - // One layer matrix - Nested s1({2}, Shape{10, 20}); - - // Two layer matrix (both modes in layer 0) - Nested s2_0({2, 0}, Shape{10, 20}); - - // Two layer matrix (one mode per layer) - Nested s1_1({1, 1}, Shape{10, 20}); - - // Two layer matrix (both modes in layer 1) - Nested s0_2({0, 2}, Shape{10, 20}); - - // One layer rank 3 - Nested s3({3}, Shape{10, 20, 30}); - - // Two layer rank 3 one mode in layer 0 two in layer 1 - Nested s1_2({1, 2}, Shape{10, 20, 30}); - - // Three layer rank 3, one mode per layer - Nested s1_1_1({1, 1, 1}, Shape{10, 20, 30}); - - // A two-layer shape where modes 0 and 1 are in layer 0 and modes 2 and 3 - // are layer 1 - Nested s({2, 2}, Shape{5, 10, 15, 20}); - -The general syntax for an |n| layer tensor is an |n| element -container where the :math:`i`-th element is the number of ranks in that -layer (ranks from the shape object are assigned to layers left to right; so -permutations may be needed to line up with layering). - -Basic Operations -================ - -All shapes know their total rank and the total number of scalar elements: - -.. code-block:: c++ - - Shape s{10, 20, 30}; - JaggedShape js{s, Shape{10, 20}}; - - // Total rank of the tensor - assert(s.rank() == 3); - assert(js.rank() == 3); - - // Total number of elements in the tensor - assert(s.size() == 6000); // 10 * 20 * 30 = 6000 - assert(js.size() == 6200); // 6000 + (10*20) = 6200; - - -``Nested`` additionally allows you to get this information per layer: - -.. code-block:: c++ - - Nested s1_2({1, 2}, s); - Nested js1_2({1, 2}, js); - - assert(s1_2.n_layers() == 2); - assert(js1_2.n_layers() == 2); - - assert(s1_2.rank_layer(0) == 1); - assert(s1_2.rank_layer(1) == 2); - assert(js1_2.rank_layer(0) == 1); - assert(js1_2.rank_layer(1) == 2); - - assert(s1_2.elements_in_layer(0) == 10); - assert(s1_2.elements_in_layer(1) == 6000); - assert(js1_2.elements_in_layer(0) == 2); - assert(js1_2.elements_in_layer(1) == 6200); +The last design, :ref:`shape_v1_design` had problems: - // Get the shape of the 0,0-th element (returns a std::variant) - assert(s3_3({0, 0}) == s); +- Passing things around as ``ShapeBase`` involved holding pointers and needing + to downcast them. Pointers are "unnatural". +- There was a lot of code duplication between the view and value objects. +- Developers had to manually synchronize the APIs of views/values, i.e., there + was no check to ensure they interfaces remained equivalent. +- Trying to introduce code factorization became difficult because the value + inherited from ``ShapeBase``, but the views didn't. +- Exposing the actual polymorphic objects meant users had to be careful to not + slice the objects. -Shape Composition -================= +.. _fig_shape_designv2: -``Shape`` and ``JaggedShape`` objects are composed similarly (with -``JaggedShape`` objects having many more checks to ensure slices are of -compatible sizes). - -.. code-block:: c++ - - Shape s0{10, 20, 30}, s1; - JaggedShape js0{Shape{10}, Shape{20}}, js1; - - // Since addition, subtraction, and element-wise multiplication only work out - // the shape of the result, they often amount to copying the state on the - // right side of the assignment operator (possibly with a permutation) - s1("i,j,k") = s0("i,j,k") + s0("i,j,k"); - assert(s1 == s0); - - js1("i,j") = js0("i,j") + js0("i,j"); - assert(js1 == js0); - - // Permuting modes - s1("j,i,k") = s0("i,j,k") + s0("i,j,k"); - assert(s1 == Shape{20,10,30}); - - js1("j,i") = js0("i,j") + js0("i,j"); - assert(js1 == JaggedShape{Shape{20}, Shape{10}}); - - // Contraction - s1("i,k") = s0("i,j,k") * s0("i,j,k"); - assert(s1 == Shape{10, 30}); - - // js0 is a jagged matrix with 2 rows, contracting over the variable number - // of columns gives a 2 by 2 matrix (represented as jagged matrix even though - // it's smooth) - js1("i,k") = js0("i,j") * js0("k,j"); - assert(js1 == JaggedShape{Shape{2}, Shape{2}}); - - // These would throw since contracted modes aren't the same length - // s1("i,k") = s0("j,i,k") * s0("i,j,k"); - - // js1("i,k") = js0("i,j") * js0("j,k"); - - // Direct product - s1("i,j,k,l") = s0("i,j,k") * s0("i,j,l"); - assert(s1 == Shape{10, 20, 30, 30}); - - js1("i,j,k") = js0("i,j") * js0("i,k"); - assert(js1 == JaggedShape{Shape{10,10}, Shape{20,20}}); - -Combining ``Nested`` objects is conceptually done layer-by-layer. In practice -we just combine the underlying ``T`` objects while preserving the layer -assignments and ensuring layer shapes are compatible: - -.. code-block:: c++ - - Shape s{10, 20, 30}; - Nested s1_2({1, 2}, s), s2_1({2,1}, s), result; - - result("i,j,k") = s1_2("i,j,k") + s1_2("i,j,k"); - assert(result == s1_2); - - // Not allowed because we can't add rank 1 tensors to rank 2 tensors - // result("i,j,k") = s1_2("i,j,k") + s2_1("i,j,k"); - - result("i,j") = s1_2("i,j,k") * s1_2("i,j,k"); - assert(result == Nested({1, 1}, Shape{10, 20})); - - result("j,k") = s1_2("i,j,k") * s1_2("i,j,k"); - assert(result == Nested({0,2}, Shape{20, 30})); - - // Layers only need compatible, not identical, shapes - result("j,k") = s1_2("i,j,k") * s2_1("i,j,k"); - assert(result == Nested({1, 1}, Shape{20, 30})); - - -We note that it's quite likely that scenarios will arise where the user will -want the result to be layered different than the default behavior provides. In -practice re-layering a shape is a trivial operation (swapping two small -vectors of integers). - -Slicing and Chipping -===================== - -Slices of a shape have the same rank, chips have different ranks: - -.. code-block:: c++ - - Shape s{10, 20}; - - // Get the shape of row 0 as a matrix - auto s0 = s.slice(0); - assert(s0 == Shape{1, 20}); - - // Get the shape of column 0 as a matrix - auto sx0 = s.slice({0, 0}, {10, 1}); - assert(sx0 == Shape{10, 1}); - - // Get the shape of the first five columns of the first five rows... - auto s05_05 = s.slice({0,0}, {5,5}); - assert(s05_05 == Shape{5, 5}); - - // Note that this shape still refers to a rank 2 tensor even though the - // first mode has a single element - auto s01_05 = s.slice({0, 0}, {1, 5}); - assert(s01_05 == Shape{1, 5}); - - // Get the shape of row 2 as a vector - auto s2 = s.chip(2); - assert(s2 == Shape{20}); - - // Get the shape of column 2 as a vector - auto sx2 = s.chip({0, 2}, {10, 3}); - assert(sx2 == Shape{10}); - -For a rank |r| tensor, the general overload of ``slice`` and ``chip`` takes -two |r|-element vectors. The first vector is the first element in the -slice/chip and the second vector is the first element not in the slice/chip. -For convenience we also provide an overload where the user may provide up to -|r| integers. This overload pins the :math:`i`-th mode to the :math:`i`-th -integer all other modes run their entire span. - -Slicing and chipping ``JaggedShape`` objects is largely the same: - -.. code-block:: c++ - - JaggedShape js0{Shape{10}, Shape{20}}; - - auto j0 = js0.chip(0); - assert(j0 == JaggedShape{Shape10}); - - auto j1 = js0.slice(0); - assert(j1 == JaggedShape({Shape{10}}); - -Because chipping selects a single element per mode per layer, chipping a -``Nested`` object is fairly straightforward: - -.. code-block:: c++ - - Nested s2_2({2, 2}, Shape{2, 2, 10, 10}); - assert(s2_2.chip(0) == Nested({1, 2}, Shape{2, 10, 10})); - assert(s2_2.chip(1,2) == Shape{10, 10}); - assert(s2_2.chip(1,2,3) == Shape{10});e - assert(s2_2.chip(1,2,3,4) == Shape{}); - -Taking arbitrary slices of a ``Nested`` object is significantly more -complicated on account of the fact that slice requests for any of the inner -layers will in general be slicing multiple tensors simultaneously. For example -consider ``s2_2`` from the previous code snippet. Slicing layer 0 is -straightforward, asking for say ``{1,0}, {2,2}`` selects row 1 of the outer -matrix. Generalizing, something like ``{1,0,5,5}, {2,2,10,10}`` would grab -row 1 of layer 0, and rows 5 through 9 (inclusive) for each inner matrix. What -if we wanted the first 5 rows of outer element ``{1,0}`` and the last 3 rows -of outer element ``{1,1}``? This request requires more than just a block range, -it requires a ``JaggedShape``. The above request can be requested by: - -.. code-block:: c++ - - // The shape resulting from taking the first 5 rows of a 10 by 10 matrix - Shape e10({5,10}, {0,0}); - // The shape resulting from taking the last 5 rows from a 10 by 10 matrix - Shape e11({3,10}, {7,0}); - // A 1 row matrix with 2 columns whose elements are a 5 by 10 and a - // 3 by 10 matrix, the origin of the outer tensor is {1,0} - JaggedShape slice({{e10, e11}}, {1,0}); - - auto requested_slice = s2_2.slice(slice); - - assert(requested_slice == slice); - -As this also shows, requesting such slices also completely negates the point of -the ``slice`` member because the input is the result. As a result, we have not -designed such an API. Instead the slicing APIs for a ``Nested`` object are: - -.. code-block:: c++ - - Nested s2_2({2, 2}, Shape{2, 2, 10, 10}); - Shape s01({1, 1, 10, 10}, {0, 1, 0, 0}); - - // Grabs the 0,1 element of the outer matrix preserving the rank - assert(s2_2.slice(0, 1) == Nested({2, 2}, s01)); - - // Grabs the bottom row of the outer matrix, and the bottom 5 rows of - // the inner matrices - Nested s1050({1, 2, 5, 5}, {1, 0, 5, 0}); - s2_2.slice({1, 0, 5, 0}, {2, 2, 10, 10}); - - -Ultimately, bear in mind, chipping and slicing are little more than convenience -functions for working out the shapes resulting from slicing/chipping the tensor; -for complicated selections it should always be possible to build the resulting -shape manually. - -Iterating -========= - -By default the origin of a freshly constructed shape is the zero vector. For -slices and chips, the origin is the first element in the slice or chip (note -that in the previous section we conveniently chose our slices/chips so the -origin was the zero vector). By default, when iterating over a shape indices are -returned as offsets from the origin, in lexicographical order. For example: - -.. code-block:: c++ - - auto print_shape = [](auto&& s){ - for(const auto& index : s) - std::cout << "{" << index[0] << "," << index[1] << "} "; - }; - - Shape s{2, 3}; - print_shape(s); // prints {0,0} {0,1} {0,2} {1,0} {1,1} {1,2} - - auto s01_13 = s.slice({0, 1}, {1, 3}); - print_shape(s01_13); // prints {0,1}, {0,2} NOT {0,0} {0,1} - - // If we wanted {0,0} {0,1} - print_shape(s01_13.offsets()); - - // We can move the origin - s.set_origin({10, 10}); - print_shape(s); // prints {10,10} {10,11} {10,12} {11,10} {11,11} {11,12} - -For completeness we define overloads of ``Shape`` and ``JaggedShape`` which -also take an origin. For ``JaggedShape`` the origin needs to be specified for -the internal shapes and the explicitly unrolled ranks. - -.. code-block:: c++ - - // Makes a shape for a 2 by 3 matrix whose first element is {10, 10} - Shape s({2, 3}, {10, 10}); - - // Outer vector starts at 10, element 11 of the outer vector starts at 5, - // element 12 of the outer vector starts at 6 - JaggedShape js({{Shape({10}, {5}), Shape({20}, {6})}, {10}); - - // Outer vector starts at 10, inner vectors start at 10. - Nested s1_1({1, 1}, s); - -******* -Summary -******* - -The design of the shape component satisfies the considerations raised above -by: - -:ref:`shape_rank_and_extents` - The ``ShapeBase`` class will provide a common API for getting/setting basic - information and performing common operations. - -:ref:`shape_nested` - The ``Nested`` class tracks how the modes of a tensor are layered. - -:ref:`shape_jagged` - The ``JaggedShape`` class is used to represent jagged shapes. - -:ref:`shape_combining_shapes` - The ``IndexedShape`` class allows us to easily compose shapes. +.. figure:: assets/shape_design_v2.png + :align: center -:ref:`shape_iterable` - The various classes define iterators which allow users to iterate over the - indices contained in the shape. + The architecture of TensorWrapper's Shape component. -**************** -Additional Notes -**************** +What this design changes: -- Slicing and chipping assume contiguous sub-tensors. For grabbing noncontiguous - sub-blocks and using them as if they were contiguous, one needs a mask. +- Moves to a "type-erased" architecture. +- Uses CRTP to factor out common APIs. +- Better separation of user-API and performance details. diff --git a/docs/source/developer/design/shapev1.rst b/docs/source/developer/design/shapev1.rst new file mode 100644 index 00000000..922975b5 --- /dev/null +++ b/docs/source/developer/design/shapev1.rst @@ -0,0 +1,624 @@ +.. Copyright 2026 NWChemEx-Project +.. +.. Licensed under the Apache License, Version 2.0 (the "License"); +.. you may not use this file except in compliance with the License. +.. You may obtain a copy of the License at +.. +.. http://www.apache.org/licenses/LICENSE-2.0 +.. +.. Unless required by applicable law or agreed to in writing, software +.. distributed under the License is distributed on an "AS IS" BASIS, +.. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +.. See the License for the specific language governing permissions and +.. limitations under the License. + +.. _shape_v1_design: + +############### +Shape V1 Design +############### + +.. |n| replace:: :math:`n` +.. |r| replace:: :math:`r` + +In designing the class hierarchy we note the following: + +- An algorithm which works for a jagged shape should work for a smooth shape + as well. The reverse, smooth algorithms with jagged shapes, will in general + not work. +- Nestings are logically imposed over an existing shape. The resulting nested + shape is still an instance of the underlying shape. +- Tiled shapes are a subcategory of jagged shapes. + +.. _fig_shape_design: + +.. figure:: assets/shape.png + :align: center + + The architecture of TensorWrapper's Shape component. + +Hi there +:numref:`fig_shape_design` shows the classes primarily responsible for +implementing the shape component. Most end users will deal with the ``Shape`` +class. + +ShapeBase +========= + +The unifying features of all shapes were summarized in the +:ref:`shape_rank_and_extents` consideration. ``ShapeBase`` provides the API +that all shapes must minimally satisfy because they are shapes. The actual +class serves primarily as code factorization. + +Shape +===== + +The ``Shape`` class describes a (smooth) hyper-rectangular array of data and +can be used for "traditional" tensors (those which are not nested or jagged). +Most end users will simply create ``Shape`` objects and pass them on to +``TensorWrapper``. We expect that manipulations of ``Shape`` objects will be +aimed at TensorWrapper developers. + +JaggedShape +=========== + +To satisfy the :ref:`shape_jagged` consideration we introduce ``JaggedShape``. +``JaggedShape`` is similar to ``Shape`` except that users must explicitly +provide the shape of the slices. Generally speaking ``JaggedShape`` objects of +rank |r| will contain a series of rank :math:`s` ``Shape`` objects. The +actual ``JaggedShape`` object serves as a map from an index with :math:`(r-s)` +indices to the ``Shape`` of that slice. Like ``Shape`` we expect users to +primarily be concerned with construction. Again, manipulations of the +``JaggedShape`` will be primarily of interest to TensorWrapper developers. + +TiledShape +========== + +Introduced primarily as a convenience for constructing ``JaggedShape`` objects +by tiling. + +Nested +====== + +To address the :ref:`shape_nested` consideration, we have added a ``Nested`` +class which is templated on the type of the tensor's overall shape. + +With objects like ``Shape`` and ``JaggedShape``, TensorWrapper can't tell how +the user is thinking of the tensor. For example, they could be thinking of a +matrix as a matrix or as a vector of vectors. The point of the ``Nested`` +object is to partition the ranks of the tensor into layers, so we know how many +layers the user is viewing the tensor as, and how many ranks each layer has. +Mathematically the various ways of a viewing a tensor do not change the +properties of the tensor; however, when we are physically laying the tensor out +on the computer, how we view the tensor can affect physical layout. + +IndexedShape +============ + +Consideration :ref:`shape_combining_shapes` requires us to be able to compose +the various shape objects. To do this, we rely on the same mechanism used for +``TensorWrapper``, *i.e.*, an expression layer. More specifically, +``IndexedShape`` objects result from indexing a shape like ``s("i,j,k")``. While +``IndexedShape`` is technically exposed to the user, user can be somewhat +oblivious to its existence. + +******************* +Proposed Shape APIs +******************* + +Constructing a ``Shape`` +======================== + +Creating a non-nested shape just requires knowing the extent of each mode: + +.. code-block:: c++ + + Shape null_shape; // No rank and no elements + Shape rank0_shape{}; // A scalar + Shape rank1_shape{10}; // 10 element vector + Shape rank2_shape{10, 20}; // 10 by 20 matrix + Shape rank3_shape{10, 20, 30}; // 10 by 20 by 30 tensor + +Note that following usual C++ rules the first two lines actually call +different constructors (default ctor vs. initializer list). Using an initializer +list requires us to know the rank at compile time. If we want to determine the +rank at runtime we can use iterators: + +.. code-block:: c++ + + // Somehow create a vector of extents + using size_type = Shape::size_type; + std::vector extents = get_extents(); + + // Construct Shape from iterator pair + Shape runtime_rank_shape(extents.begin(), extents.end()); + +Jagged Shape Construction +========================= + +For a ``Shape`` we need to specify the extents of each mode. ``JaggedShape`` +declaration is done in terms of ``Shape`` objects and looks like: + +.. code-block:: c++ + + // For brevity define variables + Shape s10{10}, s20{20}, s30{30}; + Shape s10_20{10, 20}, s30_40{30, 40}, s50_60{50, 60}; + Shape s10_20_30{10, 20, 30}, s40_50_60{40, 50, 60}; + + // No elements, no rank + JaggedShape null_shape; + + // Smooth scalar viewed as a JaggedShape (note () not {}) + JaggedShape smooth0_shape(Shape{}); + + // Smooth vector viewed as a JaggedShape (note () not {}) + JaggedShape smooth1_shape(Shape{10}); + + // Smooth matrix viewed as a JaggedShape (note () not {}) + JaggedShape smooth2_shape(Shape{10, 20}); + + // A "jagged" vector with no elements + JaggedShape rank0_shape{}; + + // A jagged matrix with 1 row, note the {} + JaggedShape rank2_shape{s10}; + + // A jagged matrix with 3 rows; row 0 has 10 elements, row 1 has 20, row 2 30 + JaggedShape rank2_shape2{s10, s20, s30}; + + // A jagged rank 3 tensor with smooth matrices. Matrix 0 is 10 by 20, + // matrix 1 is 30 by 40, and matrix 2 is 50 by 60 + JaggedShape rank3_shape{s10_20, s30_40, s50_60}; + + // A jagged rank 3 tensor where elements are jagged matrices. Matrix 0 is + // 1 by 10, matrix 2 has 20 columns in row 0 and 30 columns in row 2, and + // matrix 3 has 30 columns in row 0, 10 columns in row 1, and 20 columns in + // row 2 + JaggedShape rank3_shape2{{s10}, + {s20, s30}, + {s30, s10, s20}}; + + // A jagged rank 4 tensor where the 0-th element of the 0-th mode is a + // 10 by 20 by 30 smooth tensor and the 1-st element is a 40 by 50 by 60 + // smooth tensor + JaggedShape rank4_shape{s10_20_30, s40_50_60}; + + // A jagged rank 4 tensor where the elements are jagged rank 3 tensors. + // Taking slices along the 0 and 1-st modes, the (0,0)-th slice is a 10 by 20 + // matrix, the (0,1)-th slice is a 30 by 40 matrix, the (1,0)-th slice is + // a 30 by 40 matrix, the (1,1)-th slice is a 10 by 20 matrix, and the + // (1,2)-th slice is a 50 by 60 matrix + JaggedShape rank4_shape2{{s10_20, s30_40}, + {s30_40, s10_20, s50_60}}; + + // A jagged rank 4 tensors where the elements are jagged rank 3 tensors, + // which have jagged matrices for elements. Taking slices along the 0, 1, and + // 2 modes we have: + // - (0,0,0) is a 10 element vector, + // - (0,1,0) is a 20 element vector, + // - (0,1,1) is a 30 element vector, + // - (1,0,0) is a 10 element vector, + // - (1,0,1) is a 30 element vector, + // - (1,1,0) is a 20 element vector, + // - (1,2,0) is a 10 element vector, + // - (1,2,1) is a 20 element vector, + // - (1,2,2) is a 30 element vector + JaggedShape rank4_shape{{{s10}, {s20, s30}}, + {{s10, s30}, {s20}, {s10, s20, s30}}}; + + +Consider the shape of the (0,1) slice of ``rank4_shape``. This slice is a +vector of vectors where the outer vector has two elements, element 0 of the +outer vector is a 10-element vector and element 1 is a 30-element vector. In +other words the shape of the (0,1) slice of ``rank4_shape`` describes a jagged +matrix, which could have been initialized by ``JaggedShape{s20, s30}``. In turn +the above construction of ``rank4_shape`` is actually equivalent to: + +.. code-block:: c++ + + JaggedShape e00{s10}; + JaggedShape e01{s20, s30}; + JaggedShape e10{s10, s30}; + JaggedShape e11{s20}; + JaggedShape e12{s10, s20, s30}; + JaggedShape e0{e00, e01}; + JaggedShape e1{e10, e11, e12}; + JaggedShape rank4_shape{e0, e1}; + +And we see that ``JaggedShape`` is a recursive structure and thus the runtime +mechanism for initializing a ``JaggedShape`` is with iterators running over +``JaggedShape`` objects: + +.. code-block:: c++ + + std::vector slice_shapes = get_slices(); + JaggedShape shape(slice_shapes.begin(), slice_shapes.end()); + +So far we have focused on the most general way to create a ``JaggedShape`` one +of the most common ways to form a ``JaggedShape`` is by tiling. Consider a +30 by 30 matrix where we tile each mode into 5, 15, and 10 element chunks. +Using ``JaggedShape`` this can be done by: + +.. code-block:: c++ + + JaggedShape js{{Shape{5,5}, Shape{5, 15}, Shape{5,10}}, + {Shape{15,5}, Shape{15,15}, Shape{15,10}}, + {Shape{10,5}, Shape{10,15}, Shape{10,10}}}; + +This is an admittedly verbose declaration. Thus for the special case of crating +``JaggedShape`` objects which result from tiling smooth ``Shape`` objects we +introduce the ``TiledShape`` class. Using ``TiledShape`` the same shape could +be declared via: + +.. code-block:: c++ + + TiledShape s{{5, 10, 15}, {5, 10, 15}}; + + +Constructing Nested Shapes +========================== + +Creating a ``Nested`` object requires knowing the shape of the tensor and how +the indices are partitioned into layers. + +.. code-block:: c++ + + // Zero layer scalar + Nested s({}, Shape{}); + + // One layer scalar + Nested s0({0}, Shape{}); + + // Two layer scalar + Nested s0_0({0, 0}, Shape{}); + + // One layer vector + Nested s1({1}, Shape{10}); + + // Two layer vector (mode in layer 0) + Nested s1_0({1, 0}, Shape{10}); + + // Two layer vector (mode in layer 1) + Nested s0_1({0, 1}, Shape{10}); + + // One layer matrix + Nested s1({2}, Shape{10, 20}); + + // Two layer matrix (both modes in layer 0) + Nested s2_0({2, 0}, Shape{10, 20}); + + // Two layer matrix (one mode per layer) + Nested s1_1({1, 1}, Shape{10, 20}); + + // Two layer matrix (both modes in layer 1) + Nested s0_2({0, 2}, Shape{10, 20}); + + // One layer rank 3 + Nested s3({3}, Shape{10, 20, 30}); + + // Two layer rank 3 one mode in layer 0 two in layer 1 + Nested s1_2({1, 2}, Shape{10, 20, 30}); + + // Three layer rank 3, one mode per layer + Nested s1_1_1({1, 1, 1}, Shape{10, 20, 30}); + + // A two-layer shape where modes 0 and 1 are in layer 0 and modes 2 and 3 + // are layer 1 + Nested s({2, 2}, Shape{5, 10, 15, 20}); + +The general syntax for an |n| layer tensor is an |n| element +container where the :math:`i`-th element is the number of ranks in that +layer (ranks from the shape object are assigned to layers left to right; so +permutations may be needed to line up with layering). + +Basic Operations +================ + +All shapes know their total rank and the total number of scalar elements: + +.. code-block:: c++ + + Shape s{10, 20, 30}; + JaggedShape js{s, Shape{10, 20}}; + + // Total rank of the tensor + assert(s.rank() == 3); + assert(js.rank() == 3); + + // Total number of elements in the tensor + assert(s.size() == 6000); // 10 * 20 * 30 = 6000 + assert(js.size() == 6200); // 6000 + (10*20) = 6200; + + +``Nested`` additionally allows you to get this information per layer: + +.. code-block:: c++ + + Nested s1_2({1, 2}, s); + Nested js1_2({1, 2}, js); + + assert(s1_2.n_layers() == 2); + assert(js1_2.n_layers() == 2); + + assert(s1_2.rank_layer(0) == 1); + assert(s1_2.rank_layer(1) == 2); + assert(js1_2.rank_layer(0) == 1); + assert(js1_2.rank_layer(1) == 2); + + assert(s1_2.elements_in_layer(0) == 10); + assert(s1_2.elements_in_layer(1) == 6000); + assert(js1_2.elements_in_layer(0) == 2); + assert(js1_2.elements_in_layer(1) == 6200); + + // Get the shape of the 0,0-th element (returns a std::variant) + assert(s3_3({0, 0}) == s); + +Shape Composition +================= + +``Shape`` and ``JaggedShape`` objects are composed similarly (with +``JaggedShape`` objects having many more checks to ensure slices are of +compatible sizes). + +.. code-block:: c++ + + Shape s0{10, 20, 30}, s1; + JaggedShape js0{Shape{10}, Shape{20}}, js1; + + // Since addition, subtraction, and element-wise multiplication only work out + // the shape of the result, they often amount to copying the state on the + // right side of the assignment operator (possibly with a permutation) + s1("i,j,k") = s0("i,j,k") + s0("i,j,k"); + assert(s1 == s0); + + js1("i,j") = js0("i,j") + js0("i,j"); + assert(js1 == js0); + + // Permuting modes + s1("j,i,k") = s0("i,j,k") + s0("i,j,k"); + assert(s1 == Shape{20,10,30}); + + js1("j,i") = js0("i,j") + js0("i,j"); + assert(js1 == JaggedShape{Shape{20}, Shape{10}}); + + // Contraction + s1("i,k") = s0("i,j,k") * s0("i,j,k"); + assert(s1 == Shape{10, 30}); + + // js0 is a jagged matrix with 2 rows, contracting over the variable number + // of columns gives a 2 by 2 matrix (represented as jagged matrix even though + // it's smooth) + js1("i,k") = js0("i,j") * js0("k,j"); + assert(js1 == JaggedShape{Shape{2}, Shape{2}}); + + // These would throw since contracted modes aren't the same length + // s1("i,k") = s0("j,i,k") * s0("i,j,k"); + + // js1("i,k") = js0("i,j") * js0("j,k"); + + // Direct product + s1("i,j,k,l") = s0("i,j,k") * s0("i,j,l"); + assert(s1 == Shape{10, 20, 30, 30}); + + js1("i,j,k") = js0("i,j") * js0("i,k"); + assert(js1 == JaggedShape{Shape{10,10}, Shape{20,20}}); + +Combining ``Nested`` objects is conceptually done layer-by-layer. In practice +we just combine the underlying ``T`` objects while preserving the layer +assignments and ensuring layer shapes are compatible: + +.. code-block:: c++ + + Shape s{10, 20, 30}; + Nested s1_2({1, 2}, s), s2_1({2,1}, s), result; + + result("i,j,k") = s1_2("i,j,k") + s1_2("i,j,k"); + assert(result == s1_2); + + // Not allowed because we can't add rank 1 tensors to rank 2 tensors + // result("i,j,k") = s1_2("i,j,k") + s2_1("i,j,k"); + + result("i,j") = s1_2("i,j,k") * s1_2("i,j,k"); + assert(result == Nested({1, 1}, Shape{10, 20})); + + result("j,k") = s1_2("i,j,k") * s1_2("i,j,k"); + assert(result == Nested({0,2}, Shape{20, 30})); + + // Layers only need compatible, not identical, shapes + result("j,k") = s1_2("i,j,k") * s2_1("i,j,k"); + assert(result == Nested({1, 1}, Shape{20, 30})); + + +We note that it's quite likely that scenarios will arise where the user will +want the result to be layered different than the default behavior provides. In +practice re-layering a shape is a trivial operation (swapping two small +vectors of integers). + +Slicing and Chipping +===================== + +Slices of a shape have the same rank, chips have different ranks: + +.. code-block:: c++ + + Shape s{10, 20}; + + // Get the shape of row 0 as a matrix + auto s0 = s.slice(0); + assert(s0 == Shape{1, 20}); + + // Get the shape of column 0 as a matrix + auto sx0 = s.slice({0, 0}, {10, 1}); + assert(sx0 == Shape{10, 1}); + + // Get the shape of the first five columns of the first five rows... + auto s05_05 = s.slice({0,0}, {5,5}); + assert(s05_05 == Shape{5, 5}); + + // Note that this shape still refers to a rank 2 tensor even though the + // first mode has a single element + auto s01_05 = s.slice({0, 0}, {1, 5}); + assert(s01_05 == Shape{1, 5}); + + // Get the shape of row 2 as a vector + auto s2 = s.chip(2); + assert(s2 == Shape{20}); + + // Get the shape of column 2 as a vector + auto sx2 = s.chip({0, 2}, {10, 3}); + assert(sx2 == Shape{10}); + +For a rank |r| tensor, the general overload of ``slice`` and ``chip`` takes +two |r|-element vectors. The first vector is the first element in the +slice/chip and the second vector is the first element not in the slice/chip. +For convenience we also provide an overload where the user may provide up to +|r| integers. This overload pins the :math:`i`-th mode to the :math:`i`-th +integer all other modes run their entire span. + +Slicing and chipping ``JaggedShape`` objects is largely the same: + +.. code-block:: c++ + + JaggedShape js0{Shape{10}, Shape{20}}; + + auto j0 = js0.chip(0); + assert(j0 == JaggedShape{Shape10}); + + auto j1 = js0.slice(0); + assert(j1 == JaggedShape({Shape{10}}); + +Because chipping selects a single element per mode per layer, chipping a +``Nested`` object is fairly straightforward: + +.. code-block:: c++ + + Nested s2_2({2, 2}, Shape{2, 2, 10, 10}); + assert(s2_2.chip(0) == Nested({1, 2}, Shape{2, 10, 10})); + assert(s2_2.chip(1,2) == Shape{10, 10}); + assert(s2_2.chip(1,2,3) == Shape{10});e + assert(s2_2.chip(1,2,3,4) == Shape{}); + +Taking arbitrary slices of a ``Nested`` object is significantly more +complicated on account of the fact that slice requests for any of the inner +layers will in general be slicing multiple tensors simultaneously. For example +consider ``s2_2`` from the previous code snippet. Slicing layer 0 is +straightforward, asking for say ``{1,0}, {2,2}`` selects row 1 of the outer +matrix. Generalizing, something like ``{1,0,5,5}, {2,2,10,10}`` would grab +row 1 of layer 0, and rows 5 through 9 (inclusive) for each inner matrix. What +if we wanted the first 5 rows of outer element ``{1,0}`` and the last 3 rows +of outer element ``{1,1}``? This request requires more than just a block range, +it requires a ``JaggedShape``. The above request can be requested by: + +.. code-block:: c++ + + // The shape resulting from taking the first 5 rows of a 10 by 10 matrix + Shape e10({5,10}, {0,0}); + // The shape resulting from taking the last 5 rows from a 10 by 10 matrix + Shape e11({3,10}, {7,0}); + // A 1 row matrix with 2 columns whose elements are a 5 by 10 and a + // 3 by 10 matrix, the origin of the outer tensor is {1,0} + JaggedShape slice({{e10, e11}}, {1,0}); + + auto requested_slice = s2_2.slice(slice); + + assert(requested_slice == slice); + +As this also shows, requesting such slices also completely negates the point of +the ``slice`` member because the input is the result. As a result, we have not +designed such an API. Instead the slicing APIs for a ``Nested`` object are: + +.. code-block:: c++ + + Nested s2_2({2, 2}, Shape{2, 2, 10, 10}); + Shape s01({1, 1, 10, 10}, {0, 1, 0, 0}); + + // Grabs the 0,1 element of the outer matrix preserving the rank + assert(s2_2.slice(0, 1) == Nested({2, 2}, s01)); + + // Grabs the bottom row of the outer matrix, and the bottom 5 rows of + // the inner matrices + Nested s1050({1, 2, 5, 5}, {1, 0, 5, 0}); + s2_2.slice({1, 0, 5, 0}, {2, 2, 10, 10}); + + +Ultimately, bear in mind, chipping and slicing are little more than convenience +functions for working out the shapes resulting from slicing/chipping the tensor; +for complicated selections it should always be possible to build the resulting +shape manually. + +Iterating +========= + +By default the origin of a freshly constructed shape is the zero vector. For +slices and chips, the origin is the first element in the slice or chip (note +that in the previous section we conveniently chose our slices/chips so the +origin was the zero vector). By default, when iterating over a shape indices are +returned as offsets from the origin, in lexicographical order. For example: + +.. code-block:: c++ + + auto print_shape = [](auto&& s){ + for(const auto& index : s) + std::cout << "{" << index[0] << "," << index[1] << "} "; + }; + + Shape s{2, 3}; + print_shape(s); // prints {0,0} {0,1} {0,2} {1,0} {1,1} {1,2} + + auto s01_13 = s.slice({0, 1}, {1, 3}); + print_shape(s01_13); // prints {0,1}, {0,2} NOT {0,0} {0,1} + + // If we wanted {0,0} {0,1} + print_shape(s01_13.offsets()); + + // We can move the origin + s.set_origin({10, 10}); + print_shape(s); // prints {10,10} {10,11} {10,12} {11,10} {11,11} {11,12} + +For completeness we define overloads of ``Shape`` and ``JaggedShape`` which +also take an origin. For ``JaggedShape`` the origin needs to be specified for +the internal shapes and the explicitly unrolled ranks. + +.. code-block:: c++ + + // Makes a shape for a 2 by 3 matrix whose first element is {10, 10} + Shape s({2, 3}, {10, 10}); + + // Outer vector starts at 10, element 11 of the outer vector starts at 5, + // element 12 of the outer vector starts at 6 + JaggedShape js({{Shape({10}, {5}), Shape({20}, {6})}, {10}); + + // Outer vector starts at 10, inner vectors start at 10. + Nested s1_1({1, 1}, s); + +******* +Summary +******* + +The design of the shape component satisfies the considerations raised above +by: + +:ref:`shape_rank_and_extents` + The ``ShapeBase`` class will provide a common API for getting/setting basic + information and performing common operations. + +:ref:`shape_nested` + The ``Nested`` class tracks how the modes of a tensor are layered. + +:ref:`shape_jagged` + The ``JaggedShape`` class is used to represent jagged shapes. + +:ref:`shape_combining_shapes` + The ``IndexedShape`` class allows us to easily compose shapes. + +:ref:`shape_iterable` + The various classes define iterators which allow users to iterate over the + indices contained in the shape. + +**************** +Additional Notes +**************** + +- Slicing and chipping assume contiguous sub-tensors. For grabbing noncontiguous + sub-blocks and using them as if they were contiguous, one needs a mask. diff --git a/include/tensorwrapper/shape/shape_fwd.hpp b/include/tensorwrapper/shape/shape_fwd.hpp index cbb852f5..9486d42e 100644 --- a/include/tensorwrapper/shape/shape_fwd.hpp +++ b/include/tensorwrapper/shape/shape_fwd.hpp @@ -26,6 +26,9 @@ class SmoothViewPIMPL; class ShapeBase; +template +class SmoothCommon; + class Smooth; template diff --git a/include/tensorwrapper/shape/shape_traits.hpp b/include/tensorwrapper/shape/shape_traits.hpp index 000f753a..ff5879bd 100644 --- a/include/tensorwrapper/shape/shape_traits.hpp +++ b/include/tensorwrapper/shape/shape_traits.hpp @@ -39,28 +39,38 @@ struct ShapeTraits { using size_type = std::size_t; }; -template<> -struct ShapeTraits : public ShapeTraits { - using value_type = Smooth; +template +struct ShapeTraits> : public ShapeTraits { + using value_type = Derived; using const_value_type = const value_type; using reference = value_type&; using const_reference = const value_type&; using pointer = value_type*; using const_pointer = const value_type*; + using slice_type = Derived; }; -template<> -struct ShapeTraits : public ShapeTraits { - using value_type = Smooth; +template +struct ShapeTraits> + : public ShapeTraits { + using value_type = Derived; using const_value_type = const value_type; using reference = const value_type&; using const_reference = const value_type&; using pointer = const value_type*; using const_pointer = const value_type*; + using slice_type = Derived; }; +template<> +struct ShapeTraits : public ShapeTraits> {}; + +template<> +struct ShapeTraits + : public ShapeTraits> {}; + template -struct ShapeTraits> { +struct ShapeTraits> : public ShapeTraits> { using smooth_traits = ShapeTraits; using pimpl_type = detail_::SmoothViewPIMPL; using const_pimpl_type = diff --git a/include/tensorwrapper/shape/smooth.hpp b/include/tensorwrapper/shape/smooth.hpp index fd6cc86e..d915e278 100644 --- a/include/tensorwrapper/shape/smooth.hpp +++ b/include/tensorwrapper/shape/smooth.hpp @@ -17,6 +17,8 @@ #pragma once #include #include +#include +#include #include #include @@ -29,7 +31,11 @@ namespace tensorwrapper::shape { * geometric dimension of the (hyper-)rectangle and the number of elements in * the array. */ -class Smooth : public ShapeBase { +class Smooth : public ShapeBase, public SmoothCommon { +private: + using my_type = Smooth; + using common_base = SmoothCommon; + public: // Pull in base class's types using ShapeBase::rank_type; @@ -39,6 +45,15 @@ class Smooth : public ShapeBase { // -- Ctors, assignment, and dtor // ------------------------------------------------------------------------- + /** @brief Creates a shape for a null tensor. + * + * A null tensor is a tensor that contains zero elements. Thus a null + * tensor is not the same as scalar tensor (which contains 1 element). To + * create a scalar tensor use ``Smooth scalar{};`` (i.e., the initializer + * list ctor). + * + * @throw None No throw guarantee. + */ Smooth() noexcept = default; /** @brief Constructs *this with a statically specified number of extents. @@ -87,18 +102,6 @@ class Smooth : public ShapeBase { // -- Accessor methods // ------------------------------------------------------------------------- - /** @brief Returns the extent of the @p i -th mode. - * - * @param[in] i The mode the user wants the extent of. @p i must be in the - * range [0, rank()). - * - * @return The extent of the requested mode. - * - * @throw std::out_of_range if @p i is not in the range [0, range()). - * Strong throw guarantee. - */ - rank_type extent(size_type i) const { return m_extents_.at(i); } - // ------------------------------------------------------------------------- // -- Utility methods // ------------------------------------------------------------------------- @@ -152,6 +155,10 @@ class Smooth : public ShapeBase { } protected: + friend common_base; + + size_type extent_impl(rank_type i) const { return m_extents_.at(i); } + /// Implement clone() by calling copy ctor base_pointer clone_() const override { return std::make_unique(*this); diff --git a/include/tensorwrapper/shape/smooth_common.hpp b/include/tensorwrapper/shape/smooth_common.hpp new file mode 100644 index 00000000..a8d38fbb --- /dev/null +++ b/include/tensorwrapper/shape/smooth_common.hpp @@ -0,0 +1,216 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include + +namespace tensorwrapper::shape { + +/** @brief Code factorization for Smooth/SmoothView. + * + * @tparam Derived The class *this is implementing. Expected to be unqualified + * Smooth or SmoothView. + * + * To use this class the derived class must define: + * - `size_type extent(rank_type i) const` so that it returns the extent of + * mode i. + */ +template +class SmoothCommon { +private: + using traits_type = ShapeTraits; + +public: + using rank_type = typename traits_type::rank_type; + using size_type = typename traits_type::size_type; + using slice_type = typename traits_type::slice_type; + using slice_il_type = std::initializer_list; + + /** @brief Returns the extent of the @p i -th mode. + * + * @param[in] i The mode the user wants the extent of. @p i must be in the + * range [0, rank()). + * + * @return The extent of the requested mode. + * + * @throw std::out_of_range if @p i is not in the range [0, range()). + * Strong throw guarantee. + */ + decltype(auto) extent(rank_type i) const { + return derived().extent_impl(i); + } + + /** @brief Slices a shape given two initializer lists. + * + * C++ doesn't allow templates to work with initializer lists, therefore + * we must provide a special overload for when the input containers are + * initializer lists. This method simply dispatches to the range-based + * method by calling begin()/end() on each initializer list. See the + * description of that method for more details. + * + * @param[in] first_elem An initializer list containing the offsets of + * the first element IN the slice such that + * `first_elem[i]` is the offset along mode i. + * @param[in] last_elem An initializer list containing the offsets of + * the first element NOT IN the slice such that + * `last_elem[i]` is the offset along mode i. + * + * @return The requested slice. + * + * @throws ??? If the range-based method throws. Same throw guarantee. + */ + slice_type slice(slice_il_type first_elem, slice_il_type last_elem) const { + return slice(first_elem.begin(), first_elem.end(), last_elem.begin(), + last_elem.end()); + } + + /** @brief Slices a shape given two containers. + * + * @tparam ContainerType0 The type of first_elem. Assumed to have + * begin()/end() methods. + * @tparam ContainerType1 The type of last_elem. Assumed to have + * begin()/end() methods. + * + * Element indices are usually stored in containers. This overload is a + * convenience method for calling begin()/end() on the containers before + * dispatching to the range-based overload. See the documentation for the + * range-based overload for more details. + * + * @param[in] first_elem A container containing the offsets of + * the first element IN the slice such that + * `first_elem[i]` is the offset along mode i. + * @param[in] last_elem A container containing the offsets of + * the first element NOT IN the slice such that + * `last_elem[i]` is the offset along mode i. + * + * @return The requested slice. + * + * @throws ??? If the range-based method throws. Same throw guarantee. + */ + template + slice_type slice(ContainerType0&& first_elem, + ContainerType1&& last_elem) const { + return slice(first_elem.begin(), first_elem.end(), last_elem.begin(), + last_elem.end()); + } + + /** @brief Implements slicing given two ranges. + * + * @tparam BeginItr The type of the iterators pointing to offsets in the + * container holding the first element of the slice. + * @tparam EndItr The type of the iterators pointing to the offsets in + * the container holding the first element NOT in the + * slice. + * + * All other slice functions dispatch to this method. + * + * Slices are assumed to be contiguous, meaning we can uniquely specify + * the slice by providing the first element IN the slice and the first + * element NOT IN the slice. + * + * Specifying an element of a rank @f$r@f$ tensor requires providing + * @f$r@f$ offsets (one for each mode). Generally speaking, this requires + * the offsets to be in a container. This method takes iterators to those + * containers such that the @f$r@f$ elements in the range + * [first_elem_begin, first_elem_end) are the offsets of first element IN + * the slice and [last_elem_begin, last_elem_end) are the offsets of the + * first element NOT IN the slice. + * + * @note Both [first_elem_begin, first_elem_end) and + * [last_elem_begin, last_elem_end) being empty is allowed as long + * as *this is null or for a scalar. In these cases you will get back + * the only slice possible, which is the entire shape, i.e. a copy of + * *this. + * + * @param[in] first_elem_begin An iterator to the offset along mode 0 for + * the first element in the slice. + * @param[in] first_elem_end An iterator pointing to just past the offset + * along mode "r-1" (r being the rank of *this) for the first + * element in the slice. + * @param[in] last_elem_begin An iterator to the offset along mode 0 for + * the first element NOT in the slice. + * @param[in] last_elem_end An iterator pointing to just past the offset + * along mode "r-1" (r being the rank of *this) for the first + * element NOT in the slice. + * + * @return The requested slice. + * + * @throw std::runtime_error if the range + * [first_elem_begin, first_elem_end) does not contain the same + * number of elements as [last_elem_begin, last_elem_end). + * Strong throw guarantee. + * @throw std::runtime_error if the offsets in the range + * [first_elem_begin, first_elem_end) do not come before the + * offsets in [last_elem_begin, last_elem_end). Strong throw + * guarantee. + * + */ + template + slice_type slice(BeginItr first_elem_begin, BeginItr first_elem_end, + EndItr last_elem_begin, EndItr last_elem_end) const; + +private: + // Downcasts *this to mutable Derived object + decltype(auto) derived() { return static_cast(*this); } + + // Downcasts *this to read-only Derived object + decltype(auto) derived() const { + return static_cast(*this); + } +}; + +// ----------------------------------------------------------------------------- +// -- Out of line implementations +// ----------------------------------------------------------------------------- + +template +template +auto SmoothCommon::slice(BeginItr first_elem_begin, + BeginItr first_elem_end, + EndItr last_elem_begin, + EndItr last_elem_end) const -> slice_type { + std::vector new_extents; + + auto first_done = [&]() { return first_elem_begin == first_elem_end; }; + auto last_done = [&]() { return last_elem_begin == last_elem_end; }; + + if(first_done() && last_done()) { + // TODO: Assert rank 0 + return slice_type{}; + } else if(first_done() || last_done()) { + throw std::runtime_error("Ranges were NOT equal"); + } + + for(; !first_done(); ++first_elem_begin, ++last_elem_begin) { + if(last_done()) // Last ended before first + throw std::runtime_error("Ranges were NOT equal"); + + const auto ei = *last_elem_begin; + const auto bi = *first_elem_begin; + if(bi >= ei) + throw std::runtime_error("begin index must be < end index"); + + new_extents.push_back(ei - bi); + } + if(!last_done()) { throw std::runtime_error("Ranges were NOT equal"); } + + // TODO: assert rank == new_extents.size() + return slice_type(new_extents.begin(), new_extents.end()); +} + +} // namespace tensorwrapper::shape diff --git a/include/tensorwrapper/shape/smooth_view.hpp b/include/tensorwrapper/shape/smooth_view.hpp index cf2fc522..de40c0ac 100644 --- a/include/tensorwrapper/shape/smooth_view.hpp +++ b/include/tensorwrapper/shape/smooth_view.hpp @@ -15,6 +15,7 @@ */ #pragma once +#include #include #include @@ -30,11 +31,14 @@ namespace tensorwrapper::shape { * API to the existing state. */ template -class SmoothView { +class SmoothView : public SmoothCommon> { private: /// Type of *this using my_type = SmoothView; + /// Type of the common base + using my_base = SmoothCommon; + /// Type defining the traits for *this using traits_type = ShapeTraits; @@ -152,20 +156,6 @@ class SmoothView { /// Nothrow defaulted dtor ~SmoothView() noexcept; - /** @brief What is the extent of the i-th mode of the tensor with the - * aliased shape? - * - * @param[in] i The offset of the requested mode. @p i must be in the - * range [0, size()). - * - * @return The length of the @p i-th mode in a tensor with the aliased - * shape. - * - * @throw std::out_of_range if @p i is not in the range [0, size()). Strong - * throw guarantee. - */ - rank_type extent(size_type i) const; - /** @brief What is the rank of the tensor the aliased shape describes? * * @return The rank of the tensor with the aliased shape. @@ -227,6 +217,9 @@ class SmoothView { template friend class SmoothView; + friend my_base; + size_type extent_impl(rank_type i) const; + private: /// Type of a pointer to the PIMPL using pimpl_pointer = typename traits_type::pimpl_pointer; diff --git a/src/tensorwrapper/backends/eigen/eigen_tensor_impl.hpp b/src/tensorwrapper/backends/eigen/eigen_tensor_impl.hpp index f3d0b60c..d9015944 100644 --- a/src/tensorwrapper/backends/eigen/eigen_tensor_impl.hpp +++ b/src/tensorwrapper/backends/eigen/eigen_tensor_impl.hpp @@ -131,7 +131,8 @@ class EigenTensorImpl : public EigenTensor { auto make_from_shape_(std::span data, const_shape_reference shape, std::index_sequence) { - return eigen_data_type(data.data(), shape.extent(I)...); + using tensorwrapper::detail_::to_long; + return eigen_data_type(data.data(), to_long(shape.extent(I))...); } // Gets an element from the Eigen Tensor by unwrapping a std::vector diff --git a/src/tensorwrapper/shape/smooth_view.cpp b/src/tensorwrapper/shape/smooth_view.cpp index d5398c27..1f7e4c78 100644 --- a/src/tensorwrapper/shape/smooth_view.cpp +++ b/src/tensorwrapper/shape/smooth_view.cpp @@ -49,11 +49,6 @@ SMOOTH_VIEW& SMOOTH_VIEW::operator=(SmoothView&& rhs) noexcept = default; TPARAMS SMOOTH_VIEW::~SmoothView() noexcept = default; -TPARAMS -typename SMOOTH_VIEW::rank_type SMOOTH_VIEW::extent(size_type i) const { - return m_pimpl_->extent(i); -} - TPARAMS typename SMOOTH_VIEW::rank_type SMOOTH_VIEW::rank() const noexcept { return m_pimpl_->rank(); @@ -79,10 +74,15 @@ bool SMOOTH_VIEW::operator==(const SmoothView& rhs) const noexcept { TPARAMS auto SMOOTH_VIEW::make_smooth() const -> smooth_type { std::vector extents(rank()); - for(rank_type i = 0; i < rank(); ++i) { extents[i] = extent(i); } + for(rank_type i = 0; i < rank(); ++i) { extents[i] = this->extent(i); } return smooth_type(extents.begin(), extents.end()); } +TPARAMS +auto SMOOTH_VIEW::extent_impl(rank_type i) const -> size_type { + return m_pimpl_->extent(i); +} + TPARAMS bool SMOOTH_VIEW::has_pimpl_() const noexcept { return static_cast(m_pimpl_); diff --git a/tests/cxx/unit_tests/tensorwrapper/shape/smooth.cpp b/tests/cxx/unit_tests/tensorwrapper/shape/smooth.cpp index 208f72ee..60331521 100644 --- a/tests/cxx/unit_tests/tensorwrapper/shape/smooth.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/shape/smooth.cpp @@ -56,22 +56,6 @@ TEST_CASE("Smooth") { test_copy_move_ctor_and_assignment(scalar, vector, matrix, tensor); } - SECTION("extent") { - REQUIRE_THROWS_AS(scalar.extent(0), std::out_of_range); - - REQUIRE(vector.extent(0) == 1); - REQUIRE_THROWS_AS(vector.extent(1), std::out_of_range); - - REQUIRE(matrix.extent(0) == matrix_extents[0]); - REQUIRE(matrix.extent(1) == matrix_extents[1]); - REQUIRE_THROWS_AS(matrix.extent(2), std::out_of_range); - - REQUIRE(tensor.extent(0) == 3); - REQUIRE(tensor.extent(1) == 4); - REQUIRE(tensor.extent(2) == 5); - REQUIRE_THROWS_AS(tensor.extent(3), std::out_of_range); - } - SECTION("Virtual implementations") { SECTION("clone") { REQUIRE(scalar.clone()->are_equal(scalar)); diff --git a/tests/cxx/unit_tests/tensorwrapper/shape/smooth_common.cpp b/tests/cxx/unit_tests/tensorwrapper/shape/smooth_common.cpp new file mode 100644 index 00000000..f542d6ab --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/shape/smooth_common.cpp @@ -0,0 +1,230 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "../testing/testing.hpp" +#include +#include +#include + +namespace { + +/* All of these functions expect a vector of shapes or views of shapes such + * that: + * + * - Element 0 is the empty shape. + * - Element 1 is a scalar shape. + * - Element 2 is a vector shape. + * - Element 3 is a matrix shape. + * - Element 4 is a rank 3 tensor shape. + */ + +template +void test_extent(std::vector shapes) { + SECTION("Empty Shape") { + REQUIRE_THROWS_AS(shapes[0].extent(0), std::out_of_range); + } + + SECTION("Scalar Shape") { + REQUIRE_THROWS_AS(shapes[1].extent(0), std::out_of_range); + } + + SECTION("Vector Shape") { + REQUIRE(shapes[2].extent(0) == 3); + REQUIRE_THROWS_AS(shapes[2].extent(1), std::out_of_range); + } + + SECTION("Matrix Shape") { + REQUIRE(shapes[3].extent(0) == 2); + REQUIRE(shapes[3].extent(1) == 4); + REQUIRE_THROWS_AS(shapes[3].extent(2), std::out_of_range); + } + + SECTION("Tensor Shape") { + REQUIRE(shapes[4].extent(0) == 4); + REQUIRE(shapes[4].extent(1) == 5); + REQUIRE(shapes[4].extent(2) == 6); + REQUIRE_THROWS_AS(shapes[4].extent(3), std::out_of_range); + } +} + +template +void test_slice_il(std::vector shapes) { + // N.b. we just spot check here, full checking happens in the range test + + using smooth_type = tensorwrapper::shape::Smooth; + + smooth_type defaulted_corr; + smooth_type scalar_corr{}; + REQUIRE(shapes[0].slice({}, {}) == defaulted_corr); + + REQUIRE(shapes[1].slice({}, {}) == scalar_corr); + auto vslice = shapes[2].slice({0}, {2}); + REQUIRE(vslice == smooth_type{2}); + + auto mslice = shapes[3].slice({1, 1}, {3, 2}); + REQUIRE(mslice == smooth_type{2, 1}); + + auto tslice = shapes[4].slice({0, 0, 0}, {4, 5, 6}); + REQUIRE(tslice == shapes[4]); +} + +template +void test_slice_container(std::vector shapes) { + // N.b. we just spot check here, full checking happens in the range test + + using smooth_type = tensorwrapper::shape::Smooth; + using size_type = typename smooth_type::size_type; + using size_vector = std::vector; + + size_vector empty; + smooth_type defaulted_corr; + smooth_type scalar_corr{}; + REQUIRE(shapes[0].slice(empty, empty) == defaulted_corr); + + REQUIRE(shapes[1].slice(empty, empty) == scalar_corr); + + size_vector i0{0}, i2{2}; + auto vslice = shapes[2].slice(i0, i2); + REQUIRE(vslice == smooth_type{2}); + + size_vector i11{1, 1}, i32{3, 2}; + auto mslice = shapes[3].slice(i11, i32); + REQUIRE(mslice == smooth_type{2, 1}); + + size_vector i000{0, 0, 0}, i456{4, 5, 6}; + auto tslice = shapes[4].slice(i000, i456); + REQUIRE(tslice == shapes[4]); +} + +template +void test_slice_ranges(std::vector shapes) { + using smooth_type = tensorwrapper::shape::Smooth; + using size_type = typename smooth_type::size_type; + using size_vector = std::vector; + + size_vector empty; + auto eb = empty.begin(); + auto ee = empty.end(); + + size_vector i0{0}, i2{2}; + auto i0b = i0.begin(); + auto i0e = i0.end(); + auto i2b = i2.begin(); + auto i2e = i2.end(); + + size_vector i11{1, 1}, i32{3, 2}; + auto i11b = i11.begin(); + auto i11e = i11.end(); + auto i32b = i32.begin(); + auto i32e = i32.end(); + + size_vector i000{0, 0, 0}, i456{4, 5, 6}; + auto i000b = i000.begin(); + auto i000e = i000.end(); + auto i456b = i456.begin(); + auto i456e = i456.end(); + + using except_t = std::runtime_error; + + smooth_type defaulted_corr; + smooth_type scalar_corr{}; + smooth_type tensor_corr{4, 5, 6}; + + SECTION("defaulted") { + REQUIRE(shapes[0].slice(eb, ee, eb, ee) == defaulted_corr); + } + + SECTION("Scalar") { + REQUIRE(shapes[1].slice(eb, ee, eb, ee) == scalar_corr); + } + + SECTION("Vector") { + REQUIRE(shapes[2].slice(i0b, i0e, i2b, i2e) == smooth_type{2}); + } + + SECTION("matrix") { + REQUIRE(shapes[3].slice(i11b, i11e, i32b, i32e) == smooth_type{2, 1}); + } + + SECTION("tensor") { + REQUIRE(shapes[4].slice(i000b, i000e, i456b, i456e) == tensor_corr); + } + + SECTION("Different size ranges") { + REQUIRE_THROWS_AS(shapes[0].slice(i0b, i0e, eb, ee), except_t); + REQUIRE_THROWS_AS(shapes[1].slice(i0b, i0e, eb, ee), except_t); + + // Catch it in first preliminary check + REQUIRE_THROWS_AS(shapes[2].slice(i0b, i0e, eb, ee), except_t); + + // // Catch it in the loop + REQUIRE_THROWS_AS(shapes[3].slice(i11b, i11e, i2b, i2e), except_t); + + // Catch it after the loop + REQUIRE_THROWS_AS(shapes[4].slice(i0b, i0e, i11b, i11e), except_t); + } + + SECTION("Last element < first element") { + REQUIRE_THROWS_AS(shapes[3].slice(i2b, i2e, i0b, i0e), except_t); + } +} + +} // namespace + +TEST_CASE("SmoothCommon", "shape") { + using smooth_type = tensorwrapper::shape::Smooth; + using smooth_view = tensorwrapper::shape::SmoothView; + using const_view = tensorwrapper::shape::SmoothView; + + std::vector shapes; + smooth_type defaulted; + shapes.push_back(defaulted); + shapes.push_back(smooth_type{}); + shapes.push_back(smooth_type{3}); + shapes.push_back(smooth_type{2, 4}); + shapes.push_back(smooth_type{4, 5, 6}); + + std::vector shape_views; + std::vector const_shape_views; + for(std::size_t i = 0; i < shapes.size(); ++i) { + shape_views.push_back(smooth_view(shapes[i])); + const_shape_views.push_back(const_view(shapes[i])); + } + + SECTION("extent") { + test_extent(shapes); + test_extent(shape_views); + test_extent(const_shape_views); + } + + SECTION("slice(initializer_lists)") { + test_slice_il(shapes); + test_slice_il(shape_views); + test_extent(const_shape_views); + } + + SECTION("slice(containers)") { + test_slice_container(shapes); + test_slice_container(shape_views); + test_slice_container(const_shape_views); + } + + SECTION("slice(ranges)") { + test_slice_ranges(shapes); + test_slice_ranges(shape_views); + test_slice_ranges(const_shape_views); + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/shape/smooth_view.cpp b/tests/cxx/unit_tests/tensorwrapper/shape/smooth_view.cpp index 654557c4..147e1769 100644 --- a/tests/cxx/unit_tests/tensorwrapper/shape/smooth_view.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/shape/smooth_view.cpp @@ -73,13 +73,6 @@ TEMPLATE_LIST_TEST_CASE("SmoothView", "", types2test) { } } - SECTION("extent") { - REQUIRE_THROWS_AS(alias_scalar.extent(0), std::out_of_range); - - REQUIRE(alias_vector.extent(0) == 3); - REQUIRE_THROWS_AS(alias_vector.extent(1), std::out_of_range); - } - SECTION("rank") { REQUIRE(alias_scalar.rank() == rank_type(0)); REQUIRE(alias_vector.rank() == rank_type(1));