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 <units/units.hpp>
       6             : 
       7             : #include "scipp/variable/variable.h"
       8             : 
       9             : #include "scipp/core/dtype.h"
      10             : #include "scipp/core/eigen.h"
      11             : #include "scipp/variable/arithmetic.h"
      12             : #include "scipp/variable/except.h"
      13             : #include "scipp/variable/variable_concept.h"
      14             : #include "scipp/variable/variable_factory.h"
      15             : 
      16             : namespace scipp::variable {
      17             : 
      18             : namespace {
      19             : /// Types that default to unit=units::dimensionless. Everything else is
      20             : /// units::none.
      21             : const std::tuple<double, float, int64_t, int32_t, core::time_point,
      22             :                  Eigen::Vector3d, Eigen::Matrix3d, Eigen::Affine3d,
      23             :                  core::Quaternion, core::Translation>
      24             :     default_dimensionless_dtypes;
      25             : 
      26             : template <class... Ts>
      27      202585 : bool is_default_dimensionless(DType type, std::tuple<Ts...>) {
      28      202585 :   return ((type == dtype<Ts>) || ...);
      29             : }
      30             : } // namespace
      31             : 
      32      202585 : units::Unit default_unit_for(const DType type) {
      33      202585 :   return is_default_dimensionless(type, default_dimensionless_dtypes)
      34      202585 :              ? units::dimensionless
      35      202585 :              : units::none;
      36             : }
      37             : 
      38             : /// Construct from parent with same dtype, unit, and has_variances but new dims.
      39             : ///
      40             : /// In the case of bucket variables the buffer size is set to zero.
      41       43329 : Variable::Variable(const Variable &parent, const Dimensions &dims)
      42       43329 :     : m_dims(dims), m_strides(dims),
      43       43329 :       m_object(parent.data().makeDefaultFromParent(dims.volume())) {}
      44             : 
      45             : // TODO there is no size check here!
      46     1714731 : Variable::Variable(const Dimensions &dims, VariableConceptHandle data)
      47     1714731 :     : m_dims(dims), m_strides(dims), m_object(std::move(data)) {}
      48             : 
      49           2 : Variable::Variable(const llnl::units::precise_measurement &m)
      50           2 :     : Variable(m.value() * units::Unit(m.units())) {}
      51             : 
      52             : namespace {
      53     1804412 : void check_nested_in_assign(const Variable &lhs, const Variable &rhs) {
      54     1804412 :   if (!rhs.is_valid() || rhs.dtype() != dtype<Variable>) {
      55     1804383 :     return;
      56             :   }
      57             :   // In principle we should also check when the RHS contains DataArrays or
      58             :   // Datasets. But those are copied when stored in Variables,
      59             :   // so no check needed here.
      60          49 :   for (const auto &nested : rhs.values<Variable>()) {
      61          35 :     if (&lhs == &nested) {
      62          10 :       throw std::invalid_argument("Cannot assign Variable, the right hand side "
      63             :                                   "contains a reference to the left hand side. "
      64          20 :                                   "Reference cycles are not allowed.");
      65             :     }
      66          25 :     check_nested_in_assign(lhs, nested);
      67          29 :   }
      68             : }
      69             : } // namespace
      70             : 
      71         884 : Variable &Variable::operator=(const Variable &other) {
      72         884 :   return *this = Variable(other);
      73             : }
      74             : 
      75     1804387 : Variable &Variable::operator=(Variable &&other) {
      76     1804387 :   check_nested_in_assign(*this, other);
      77     1804377 :   m_dims = other.m_dims;
      78     1804377 :   m_strides = other.m_strides;
      79     1804377 :   m_offset = other.m_offset;
      80     1804377 :   m_object = std::move(other.m_object);
      81     1804377 :   m_readonly = other.m_readonly;
      82     1804377 :   m_aligned = other.m_aligned;
      83     1804377 :   return *this;
      84             : }
      85             : 
      86      164276 : void Variable::setDataHandle(VariableConceptHandle object) {
      87      164276 :   if (object->size() != m_object->size())
      88           1 :     throw except::DimensionError("Cannot replace by model of different size.");
      89      164275 :   m_object = object;
      90      164275 : }
      91             : 
      92    22886633 : const Dimensions &Variable::dims() const {
      93    22886633 :   if (!is_valid())
      94           1 :     throw std::runtime_error("invalid variable");
      95    22886632 :   return m_dims;
      96             : }
      97             : 
      98      170939 : scipp::index Variable::ndim() const {
      99      170939 :   if (!is_valid())
     100           0 :     throw std::runtime_error("invalid variable");
     101      170939 :   return m_dims.ndim();
     102             : }
     103             : 
     104    53704129 : DType Variable::dtype() const { return data().dtype(); }
     105             : 
     106     9336201 : bool Variable::has_variances() const { return data().has_variances(); }
     107             : 
     108     1337596 : void Variable::expect_can_set_unit(const units::Unit &unit) const {
     109     1337596 :   if (this->unit() != unit && is_slice())
     110          17 :     throw except::UnitError("Partial view on data of variable cannot be used "
     111          34 :                             "to change the unit.");
     112     1337579 : }
     113             : 
     114     6420701 : const units::Unit &Variable::unit() const { return m_object->unit(); }
     115             : 
     116      867315 : void Variable::setUnit(const units::Unit &unit) {
     117      867315 :   expect_writable();
     118      867304 :   expect_can_set_unit(unit);
     119      867294 :   m_object->setUnit(unit);
     120      867290 : }
     121             : 
     122       40504 : Dim Variable::dim() const {
     123       40504 :   core::expect::ndim_is(dims(), 1);
     124       40497 :   return dims().inner();
     125             : }
     126             : 
     127             : namespace {
     128      424263 : bool compare(const Variable &a, const Variable &b, bool equal_nan) {
     129      424263 :   if (equal_nan && a.is_same(b))
     130        1402 :     return true;
     131      422861 :   if (!a.is_valid() || !b.is_valid())
     132           2 :     return a.is_valid() == b.is_valid();
     133             :   // Note: Not comparing strides
     134      422859 :   if (a.unit() != b.unit())
     135         244 :     return false;
     136      422615 :   if (a.dims() != b.dims())
     137         557 :     return false;
     138      422058 :   if (a.dtype() != b.dtype())
     139         162 :     return false;
     140      421896 :   if (a.has_variances() != b.has_variances())
     141         157 :     return false;
     142      421739 :   if (a.dims().volume() == 0 && a.dims() == b.dims())
     143        7059 :     return true;
     144      829360 :   return a.dims() == b.dims() &&
     145      829360 :          (equal_nan ? a.data().equals_nan(a, b) : a.data().equals(a, b));
     146             : }
     147             : } // namespace
     148             : 
     149      420404 : bool Variable::operator==(const Variable &other) const {
     150      420404 :   return compare(*this, other, false);
     151             : }
     152             : 
     153        3859 : bool equals_nan(const Variable &a, const Variable &b) {
     154        3859 :   return compare(a, b, true);
     155             : }
     156             : 
     157        7566 : bool Variable::operator!=(const Variable &other) const {
     158        7566 :   return !(*this == other);
     159             : }
     160             : 
     161    70699536 : const VariableConcept &Variable::data() const & { return *m_object; }
     162             : 
     163     3695273 : VariableConcept &Variable::data() & {
     164     3695273 :   expect_writable();
     165     3695255 :   return *m_object;
     166             : }
     167             : 
     168      211189 : const VariableConceptHandle &Variable::data_handle() const { return m_object; }
     169             : 
     170      442708 : scipp::span<const scipp::index> Variable::strides() const {
     171      442708 :   return scipp::span<const scipp::index>(&*m_strides.begin(),
     172      885416 :                                          &*m_strides.begin() + dims().ndim());
     173             : }
     174             : 
     175      803813 : scipp::index Variable::stride(const Dim dim) const {
     176      803813 :   return m_strides[dims().index(dim)];
     177             : }
     178             : 
     179           9 : scipp::index Variable::offset() const { return m_offset; }
     180             : 
     181     8304311 : core::ElementArrayViewParams Variable::array_params() const {
     182     8304311 :   return {m_offset, dims(), m_strides, {}};
     183             : }
     184             : 
     185      809378 : Variable Variable::slice(const Slice params) const {
     186      809378 :   core::expect::validSlice(dims(), params);
     187      809278 :   Variable out(*this);
     188      809278 :   if (params == Slice{})
     189        3751 :     return out;
     190      805527 :   const auto dim = params.dim();
     191      805527 :   const auto begin = params.begin();
     192      805527 :   const auto end = params.end();
     193      805527 :   const auto stride = params.stride();
     194      805527 :   const auto index = out.m_dims.index(dim);
     195      805527 :   out.m_offset += begin * m_strides[index];
     196      805527 :   if (end == -1) {
     197      118870 :     out.m_strides.erase(index);
     198      118870 :     out.m_dims.erase(dim);
     199             :   } else {
     200      686657 :     static_cast<Sizes &>(out.m_dims) = out.m_dims.slice(params);
     201      686657 :     out.m_strides[index] *= stride;
     202             :   }
     203      805527 :   return out;
     204           0 : }
     205             : 
     206        3990 : void Variable::validateSlice(const Slice &s, const Variable &data) const {
     207        3990 :   core::expect::validSlice(this->dims(), s);
     208        3990 :   if (variableFactory().has_variances(data) !=
     209        3990 :       variableFactory().has_variances(*this)) {
     210           2 :     auto variances_message = [](const auto &variable) {
     211             :       return "does" +
     212           4 :              std::string(variableFactory().has_variances(variable) ? ""
     213             :                                                                    : " NOT") +
     214           4 :              " have variances.";
     215             :     };
     216           2 :     throw except::VariancesError("Invalid slice operation. Slice " +
     217           3 :                                  variances_message(data) + "Variable " +
     218           2 :                                  variances_message(*this));
     219             :   }
     220        3989 :   if (variableFactory().elem_unit(data) != variableFactory().elem_unit(*this))
     221           1 :     throw except::UnitError(
     222           1 :         "Invalid slice operation. Slice has unit: " +
     223           3 :         to_string(variableFactory().elem_unit(data)) +
     224           4 :         " Variable has unit: " + to_string(variableFactory().elem_unit(*this)));
     225        3988 :   if (variableFactory().elem_dtype(data) != variableFactory().elem_dtype(*this))
     226           2 :     throw except::TypeError("Invalid slice operation. Slice has dtype " +
     227           3 :                             to_string(variableFactory().elem_dtype(data)) +
     228           2 :                             ". Variable has dtype " +
     229           2 :                             to_string(variableFactory().elem_dtype(*this)));
     230        3987 : }
     231             : 
     232        3929 : Variable &Variable::setSlice(const Slice params, const Variable &data) {
     233        3929 :   validateSlice(params, data);
     234        3929 :   copy(data, slice(params));
     235        3915 :   return *this;
     236             : }
     237             : 
     238       84716 : Variable Variable::broadcast(const Dimensions &target) const {
     239       84716 :   expect::includes(target, dims());
     240       84714 :   auto out = target.volume() == dims().volume() ? *this : as_const();
     241       84714 :   out.m_dims = target;
     242       84714 :   out.m_strides.clear();
     243      150618 :   for (const auto &d : target.labels())
     244       65904 :     out.m_strides.push_back(dims().contains(d) ? m_strides[dims().index(d)]
     245             :                                                : 0);
     246       84714 :   return out;
     247           0 : }
     248             : 
     249         259 : Variable Variable::fold(const Dim dim, const Dimensions &target) const {
     250         259 :   auto out(*this);
     251         259 :   out.m_dims = core::fold(dims(), dim, target);
     252         258 :   out.m_strides.clear();
     253         258 :   const Strides substrides(target);
     254         555 :   for (scipp::index i_in = 0; i_in < dims().ndim(); ++i_in) {
     255         297 :     if (dims().label(i_in) == dim)
     256         807 :       for (scipp::index i_target = 0; i_target < target.ndim(); ++i_target)
     257         549 :         out.m_strides.push_back(m_strides[i_in] * substrides[i_target]);
     258             :     else
     259          39 :       out.m_strides.push_back(m_strides[i_in]);
     260             :   }
     261         516 :   return out;
     262         259 : }
     263             : 
     264        7453 : Variable Variable::transpose(const scipp::span<const Dim> order) const {
     265        7453 :   auto transposed(*this);
     266        7461 :   transposed.m_strides = core::transpose(m_strides, dims(), order);
     267        7445 :   transposed.m_dims = core::transpose(dims(), order);
     268        7445 :   return transposed;
     269           8 : }
     270             : 
     271             : /// Return new variable with renamed dimensions.
     272             : ///
     273             : /// The `fail_on_unknown` option is used internally for implementing rename_dims
     274             : /// in DataArray and related classes.
     275        1990 : Variable Variable::rename_dims(const std::vector<std::pair<Dim, Dim>> &names,
     276             :                                const bool fail_on_unknown) const {
     277        1990 :   auto out(*this);
     278             :   // names is a vector for predictable order
     279        1990 :   out.m_dims = out.m_dims.rename_dims(names, fail_on_unknown);
     280        1989 :   return out;
     281           1 : }
     282             : 
     283    26246831 : bool Variable::is_valid() const noexcept { return m_object.operator bool(); }
     284             : 
     285       73207 : bool Variable::is_slice() const {
     286             :   // TODO Is this condition sufficient?
     287       73207 :   return m_offset != 0 || m_dims.volume() != data().size();
     288             : }
     289             : 
     290      137723 : bool Variable::is_readonly() const noexcept { return m_readonly; }
     291             : 
     292        7033 : bool Variable::is_same(const Variable &other) const noexcept {
     293        7033 :   return std::tie(m_dims, m_strides, m_offset, m_object) ==
     294        7033 :          std::tie(other.m_dims, other.m_strides, other.m_offset,
     295       14066 :                   other.m_object);
     296             : }
     297             : 
     298      364724 : bool Variable::is_aligned() const noexcept { return m_aligned; }
     299             : 
     300      414351 : void Variable::set_aligned(bool aligned) noexcept { m_aligned = aligned; }
     301             : 
     302        3585 : void Variable::setVariances(const Variable &v) {
     303        3585 :   expect_writable();
     304        3574 :   if (is_slice())
     305          11 :     throw except::VariancesError(
     306          22 :         "Cannot add variances via sliced view of Variable.");
     307        3563 :   if (v.is_valid()) {
     308        3530 :     core::expect::equals(unit(), v.unit());
     309        3526 :     core::expect::equals(dims(), v.dims());
     310             :   }
     311        3559 :   data().setVariances(v);
     312        3554 : }
     313             : 
     314             : namespace detail {
     315           0 : void throw_keyword_arg_constructor_bad_dtype(const DType dtype) {
     316           0 :   throw except::TypeError("Cannot create the Variable with type " +
     317           0 :                           to_string(dtype) +
     318           0 :                           " with such values and/or variances.");
     319             : }
     320             : } // namespace detail
     321             : 
     322      968073 : Variable Variable::bin_indices() const {
     323      968073 :   auto out{*this};
     324      968073 :   out.m_object = data().bin_indices();
     325      968073 :   return out;
     326           0 : }
     327             : 
     328       71992 : Variable Variable::as_const() const {
     329       71992 :   Variable out(*this);
     330       71992 :   out.m_readonly = true;
     331       71992 :   return out;
     332             : }
     333             : 
     334     4566173 : void Variable::expect_writable() const {
     335     4566173 :   if (m_readonly)
     336          40 :     throw except::VariableError("Read-only flag is set, cannot mutate data.");
     337     4566133 : }
     338             : } // namespace scipp::variable