Line data    Source code

       1             : // SPDX-License-Identifier: BSD-3-Clause
       2             : // Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
       3             : /// @file
       4             : /// @author Simon Heybrock
       5             : #include <algorithm>
       6             : #include <numeric>
       7             : 
       8             : #include "rename.h"
       9             : #include "scipp/core/dimensions.h"
      10             : #include "scipp/core/except.h"
      11             : #include "scipp/core/sizes.h"
      12             : 
      13             : namespace scipp::core {
      14       74037 : Dimensions::Dimensions(const scipp::span<const Dim> labels,
      15       74037 :                        const scipp::span<const scipp::index> shape) {
      16       74037 :   if (labels.size() != shape.size())
      17           0 :     throw except::DimensionError(
      18             :         "Constructing Dimensions: Number of dimensions "
      19           0 :         "labels (" +
      20           0 :         std::to_string(labels.size()) + ") does not match shape size (" +
      21           0 :         std::to_string(shape.size()) + ").");
      22      214548 :   for (scipp::index i = 0; i < scipp::size(shape); ++i)
      23      140511 :     addInner(labels[i], shape[i]);
      24       74037 : }
      25             : 
      26     2035629 : Dim Dimensions::label(const scipp::index i) const { return labels()[i]; }
      27             : 
      28     7490975 : scipp::index Dimensions::size(const scipp::index i) const { return shape()[i]; }
      29             : 
      30             : /// Return the offset of elements along this dimension in a multi-dimensional
      31             : /// array defined by this.
      32         118 : scipp::index Dimensions::offset(const Dim label) const {
      33         118 :   scipp::index offset{1};
      34         158 :   for (int32_t i = index(label) + 1; i < ndim(); ++i)
      35          40 :     offset *= shape()[i];
      36         118 :   return offset;
      37             : }
      38             : 
      39             : /// Add a new dimension, which will be the outermost dimension.
      40         240 : void Dimensions::add(const Dim label, const scipp::index size) {
      41         240 :   expect::validDim(label);
      42         240 :   expect::validExtent(size);
      43         240 :   insert_left(label, size);
      44         237 : }
      45             : 
      46             : /// Add a new dimension, which will be the innermost dimension.
      47     1455120 : void Dimensions::addInner(const Dim label, const scipp::index size) {
      48     1455120 :   expect::validDim(label);
      49     1455120 :   expect::validExtent(size);
      50     1455120 :   insert_right(label, size);
      51     1455112 : }
      52             : 
      53             : /// Return the innermost dimension. Returns Dim::Invalid if *this is empty
      54      328215 : Dim Dimensions::inner() const noexcept {
      55      328215 :   if (empty())
      56          72 :     return Dim::Invalid;
      57      328143 :   return labels().back();
      58             : }
      59             : 
      60             : Dimensions
      61        1957 : Dimensions::rename_dims(const std::vector<std::pair<Dim, Dim>> &names,
      62             :                         const bool fail_on_unknown) const {
      63        1957 :   return detail::rename_dims(*this, names, fail_on_unknown);
      64             : }
      65             : 
      66      527887 : Dimensions merge(const Dimensions &a) { return a; }
      67             : 
      68             : /// Return the direct sum, i.e., the combination of dimensions in a and b.
      69             : ///
      70             : /// Throws if there is a mismatching dimension extent.
      71             : /// The implementation "favors" the order of the first argument if both
      72             : /// inputs have the same number of dimension. Transposing is avoided where
      73             : /// possible, which is crucial for accumulate performance.
      74     1906590 : Dimensions merge(const Dimensions &a, const Dimensions &b) {
      75     1906590 :   Dimensions out;
      76     1906590 :   auto it = b.labels().begin();
      77     1906590 :   auto end = b.labels().end();
      78     2881594 :   for (const auto &dim : a.labels()) {
      79             :     // add any labels appearing *before* dim
      80      975026 :     if (b.contains(dim)) {
      81      868654 :       if (a[dim] != b[dim])
      82          22 :         throw except::DimensionError(
      83          22 :             "Cannot merge dimensions with mismatching extent in '" +
      84          66 :             to_string(dim) + "': " + to_string(a) + " and " + to_string(b));
      85     1183721 :       while (it != end && *it != dim) {
      86      315089 :         if (!a.contains(*it))
      87        3032 :           out.addInner(*it, b[*it]);
      88      315089 :         ++it;
      89             :       }
      90             :     }
      91      975004 :     out.addInner(dim, a[dim]);
      92             :   }
      93             :   // add remaining labels appearing after last of a's labels
      94     2653896 :   while (it != end) {
      95      747328 :     if (!a.contains(*it))
      96      190753 :       out.addInner(*it, b[*it]);
      97      747328 :     ++it;
      98             :   }
      99     1906568 :   return out;
     100          22 : }
     101             : 
     102             : /// Return the dimensions contained in both a and b (dimension order is not
     103             : /// checked).
     104             : /// The convention is the same as for merge: we favor the dimension order in a
     105             : /// for dimensions found both in a and b.
     106        1571 : Dimensions intersection(const Dimensions &a, const Dimensions &b) {
     107        1571 :   Dimensions out;
     108        1571 :   Dimensions m = merge(a, b);
     109        4735 :   for (const auto &dim : m.labels()) {
     110        3164 :     if (a.contains(dim) && b.contains(dim))
     111        3133 :       out.addInner(dim, m[dim]);
     112             :   }
     113        3142 :   return out;
     114        1571 : }
     115             : 
     116             : namespace {
     117       11793 : Dimensions transpose_impl(const Dimensions &dims,
     118             :                           const scipp::span<const Dim> labels) {
     119       11793 :   if (scipp::size(labels) != dims.ndim())
     120          11 :     throw except::DimensionError("Cannot transpose: Requested new dimension "
     121          22 :                                  "order contains different number of labels.");
     122       11782 :   std::vector<scipp::index> shape(labels.size());
     123       11782 :   std::transform(labels.begin(), labels.end(), shape.begin(),
     124       27217 :                  [&dims](const auto &dim) { return dims[dim]; });
     125       23556 :   return {labels, shape};
     126       11782 : }
     127             : } // namespace
     128             : 
     129       11793 : Dimensions transpose(const Dimensions &dims,
     130             :                      const scipp::span<const Dim> labels) {
     131       11793 :   if (labels.empty()) {
     132        2099 :     std::vector<Dim> default_labels{dims.labels().rbegin(),
     133        2099 :                                     dims.labels().rend()};
     134        2099 :     return transpose_impl(dims, default_labels);
     135        2099 :   }
     136        9694 :   return transpose_impl(dims, labels);
     137             : }
     138             : 
     139             : /// Fold one dim into multiple dims
     140             : ///
     141             : /// Go through the old dims and:
     142             : /// - if the dim does not equal the dim that is being stacked, copy dim/shape
     143             : /// - if the dim equals the dim to be stacked, replace by stack of new dims
     144             : ///
     145             : /// Note that addInner will protect against inserting new dims that already
     146             : /// exist in the old dims.
     147         253 : Dimensions fold(const Dimensions &old_dims, const Dim from_dim,
     148             :                 const Dimensions &to_dims) {
     149         253 :   scipp::expect::contains(old_dims, from_dim);
     150         251 :   Dimensions new_dims;
     151         542 :   for (const auto &dim : old_dims.labels())
     152         291 :     if (dim != from_dim)
     153          40 :       new_dims.addInner(dim, old_dims[dim]);
     154             :     else
     155         762 :       for (const auto &lab : to_dims.labels())
     156         511 :         new_dims.addInner(lab, to_dims[lab]);
     157         251 :   if (old_dims.volume() != new_dims.volume())
     158           2 :     throw except::DimensionError(
     159           4 :         "Sizes " + to_string(to_dims) +
     160           4 :         " provided to `fold` not compatible with length '" +
     161           8 :         std::to_string(old_dims[from_dim]) + "' of dimension '" +
     162           8 :         from_dim.name() + "' being folded.");
     163         249 :   return new_dims;
     164           2 : }
     165             : 
     166             : } // namespace scipp::core