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       77254 : Dimensions::Dimensions(const scipp::span<const Dim> labels,
      15       77254 :                        const scipp::span<const scipp::index> shape) {
      16       77254 :   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      226451 :   for (scipp::index i = 0; i < scipp::size(shape); ++i)
      23      149197 :     addInner(labels[i], shape[i]);
      24       77254 : }
      25             : 
      26     2133629 : Dim Dimensions::label(const scipp::index i) const { return labels()[i]; }
      27             : 
      28     7599178 : 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         245 : void Dimensions::add(const Dim label, const scipp::index size) {
      41         245 :   expect::validDim(label);
      42         245 :   expect::validExtent(size);
      43         245 :   insert_left(label, size);
      44         242 : }
      45             : 
      46             : /// Add a new dimension, which will be the innermost dimension.
      47     1526448 : void Dimensions::addInner(const Dim label, const scipp::index size) {
      48     1526448 :   expect::validDim(label);
      49     1526448 :   expect::validExtent(size);
      50     1526448 :   insert_right(label, size);
      51     1526440 : }
      52             : 
      53             : /// Return the innermost dimension. Returns Dim::Invalid if *this is empty
      54      330431 : Dim Dimensions::inner() const noexcept {
      55      330431 :   if (empty())
      56          72 :     return Dim::Invalid;
      57      330359 :   return labels().back();
      58             : }
      59             : 
      60             : Dimensions
      61        1990 : Dimensions::rename_dims(const std::vector<std::pair<Dim, Dim>> &names,
      62             :                         const bool fail_on_unknown) const {
      63        1990 :   return detail::rename_dims(*this, names, fail_on_unknown);
      64             : }
      65             : 
      66      566771 : 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     1948240 : Dimensions merge(const Dimensions &a, const Dimensions &b) {
      75     1948240 :   Dimensions out;
      76     1948240 :   auto it = b.labels().begin();
      77     1948240 :   auto end = b.labels().end();
      78     2967573 :   for (const auto &dim : a.labels()) {
      79             :     // add any labels appearing *before* dim
      80     1019355 :     if (b.contains(dim)) {
      81      906711 :       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     1243059 :       while (it != end && *it != dim) {
      86      336370 :         if (!a.contains(*it))
      87        3978 :           out.addInner(*it, b[*it]);
      88      336370 :         ++it;
      89             :       }
      90             :     }
      91     1019333 :     out.addInner(dim, a[dim]);
      92             :   }
      93             :   // add remaining labels appearing after last of a's labels
      94     2725141 :   while (it != end) {
      95      776923 :     if (!a.contains(*it))
      96      202626 :       out.addInner(*it, b[*it]);
      97      776923 :     ++it;
      98             :   }
      99     1948218 :   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        2016 : Dimensions intersection(const Dimensions &a, const Dimensions &b) {
     107        2016 :   Dimensions out;
     108        2016 :   Dimensions m = merge(a, b);
     109        6394 :   for (const auto &dim : m.labels()) {
     110        4378 :     if (a.contains(dim) && b.contains(dim))
     111        4347 :       out.addInner(dim, m[dim]);
     112             :   }
     113        4032 :   return out;
     114        2016 : }
     115             : 
     116             : namespace {
     117       14921 : Dimensions transpose_impl(const Dimensions &dims,
     118             :                           const scipp::span<const Dim> labels) {
     119       14921 :   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       14910 :   std::vector<scipp::index> shape(labels.size());
     123       14910 :   std::transform(labels.begin(), labels.end(), shape.begin(),
     124       35783 :                  [&dims](const auto &dim) { return dims[dim]; });
     125       29812 :   return {labels, shape};
     126       14910 : }
     127             : } // namespace
     128             : 
     129       14921 : Dimensions transpose(const Dimensions &dims,
     130             :                      const scipp::span<const Dim> labels) {
     131       14921 :   if (labels.empty()) {
     132        2659 :     std::vector<Dim> default_labels{dims.labels().rbegin(),
     133        2659 :                                     dims.labels().rend()};
     134        2659 :     return transpose_impl(dims, default_labels);
     135        2659 :   }
     136       12262 :   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         265 : Dimensions fold(const Dimensions &old_dims, const Dim from_dim,
     148             :                 const Dimensions &to_dims) {
     149         265 :   scipp::expect::contains(old_dims, from_dim);
     150         263 :   Dimensions new_dims;
     151         566 :   for (const auto &dim : old_dims.labels())
     152         303 :     if (dim != from_dim)
     153          40 :       new_dims.addInner(dim, old_dims[dim]);
     154             :     else
     155         823 :       for (const auto &lab : to_dims.labels())
     156         560 :         new_dims.addInner(lab, to_dims[lab]);
     157         263 :   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         261 :   return new_dims;
     164           2 : }
     165             : 
     166             : } // namespace scipp::core