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      198733 : bool is_default_dimensionless(DType type, std::tuple<Ts...>) {
      28      198733 :   return ((type == dtype<Ts>) || ...);
      29             : }
      30             : } // namespace
      31             : 
      32      198733 : units::Unit default_unit_for(const DType type) {
      33      198733 :   return is_default_dimensionless(type, default_dimensionless_dtypes)
      34      198733 :              ? units::dimensionless
      35      198733 :              : 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       41971 : Variable::Variable(const Variable &parent, const Dimensions &dims)
      42       41971 :     : m_dims(dims), m_strides(dims),
      43       41971 :       m_object(parent.data().makeDefaultFromParent(dims.volume())) {}
      44             : 
      45             : // TODO there is no size check here!
      46     1642461 : Variable::Variable(const Dimensions &dims, VariableConceptHandle data)
      47     1642461 :     : 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     1732977 : void check_nested_in_assign(const Variable &lhs, const Variable &rhs) {
      54     1732977 :   if (!rhs.is_valid() || rhs.dtype() != dtype<Variable>) {
      55     1732948 :     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         879 : Variable &Variable::operator=(const Variable &other) {
      72         879 :   return *this = Variable(other);
      73             : }
      74             : 
      75     1732952 : Variable &Variable::operator=(Variable &&other) {
      76     1732952 :   check_nested_in_assign(*this, other);
      77     1732942 :   m_dims = other.m_dims;
      78     1732942 :   m_strides = other.m_strides;
      79     1732942 :   m_offset = other.m_offset;
      80     1732942 :   m_object = std::move(other.m_object);
      81     1732942 :   m_readonly = other.m_readonly;
      82     1732942 :   m_aligned = other.m_aligned;
      83     1732942 :   return *this;
      84             : }
      85             : 
      86      159175 : void Variable::setDataHandle(VariableConceptHandle object) {
      87      159175 :   if (object->size() != m_object->size())
      88           1 :     throw except::DimensionError("Cannot replace by model of different size.");
      89      159174 :   m_object = object;
      90      159174 : }
      91             : 
      92    22191342 : const Dimensions &Variable::dims() const {
      93    22191342 :   if (!is_valid())
      94           1 :     throw std::runtime_error("invalid variable");
      95    22191341 :   return m_dims;
      96             : }
      97             : 
      98      116807 : scipp::index Variable::ndim() const {
      99      116807 :   if (!is_valid())
     100           0 :     throw std::runtime_error("invalid variable");
     101      116807 :   return m_dims.ndim();
     102             : }
     103             : 
     104    51877302 : DType Variable::dtype() const { return data().dtype(); }
     105             : 
     106     9031936 : bool Variable::has_variances() const { return data().has_variances(); }
     107             : 
     108     1319689 : void Variable::expect_can_set_unit(const units::Unit &unit) const {
     109     1319689 :   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     1319672 : }
     113             : 
     114     6276348 : const units::Unit &Variable::unit() const { return m_object->unit(); }
     115             : 
     116      855892 : void Variable::setUnit(const units::Unit &unit) {
     117      855892 :   expect_writable();
     118      855881 :   expect_can_set_unit(unit);
     119      855871 :   m_object->setUnit(unit);
     120      855867 : }
     121             : 
     122       40596 : Dim Variable::dim() const {
     123       40596 :   core::expect::ndim_is(dims(), 1);
     124       40589 :   return dims().inner();
     125             : }
     126             : 
     127             : namespace {
     128      423550 : bool compare(const Variable &a, const Variable &b, bool equal_nan) {
     129      423550 :   if (equal_nan && a.is_same(b))
     130        1377 :     return true;
     131      422173 :   if (!a.is_valid() || !b.is_valid())
     132           2 :     return a.is_valid() == b.is_valid();
     133             :   // Note: Not comparing strides
     134      422171 :   if (a.unit() != b.unit())
     135         244 :     return false;
     136      421927 :   if (a.dims() != b.dims())
     137         557 :     return false;
     138      421370 :   if (a.dtype() != b.dtype())
     139         162 :     return false;
     140      421208 :   if (a.has_variances() != b.has_variances())
     141         157 :     return false;
     142      421051 :   if (a.dims().volume() == 0 && a.dims() == b.dims())
     143        7011 :     return true;
     144      828080 :   return a.dims() == b.dims() &&
     145      828080 :          (equal_nan ? a.data().equals_nan(a, b) : a.data().equals(a, b));
     146             : }
     147             : } // namespace
     148             : 
     149      419722 : bool Variable::operator==(const Variable &other) const {
     150      419722 :   return compare(*this, other, false);
     151             : }
     152             : 
     153        3828 : bool equals_nan(const Variable &a, const Variable &b) {
     154        3828 :   return compare(a, b, true);
     155             : }
     156             : 
     157        7485 : bool Variable::operator!=(const Variable &other) const {
     158        7485 :   return !(*this == other);
     159             : }
     160             : 
     161    68377021 : const VariableConcept &Variable::data() const & { return *m_object; }
     162             : 
     163     3519736 : VariableConcept &Variable::data() & {
     164     3519736 :   expect_writable();
     165     3519718 :   return *m_object;
     166             : }
     167             : 
     168      207615 : const VariableConceptHandle &Variable::data_handle() const { return m_object; }
     169             : 
     170      425723 : scipp::span<const scipp::index> Variable::strides() const {
     171      425723 :   return scipp::span<const scipp::index>(&*m_strides.begin(),
     172     1277169 :                                          &*m_strides.begin() + dims().ndim());
     173             : }
     174             : 
     175      800245 : scipp::index Variable::stride(const Dim dim) const {
     176      800245 :   return m_strides[dims().index(dim)];
     177             : }
     178             : 
     179           9 : scipp::index Variable::offset() const { return m_offset; }
     180             : 
     181     8012813 : core::ElementArrayViewParams Variable::array_params() const {
     182     8012813 :   return {m_offset, dims(), m_strides, {}};
     183             : }
     184             : 
     185      804433 : Variable Variable::slice(const Slice params) const {
     186      804433 :   core::expect::validSlice(dims(), params);
     187      804333 :   Variable out(*this);
     188      804333 :   if (params == Slice{})
     189        3751 :     return out;
     190      800582 :   const auto dim = params.dim();
     191      800582 :   const auto begin = params.begin();
     192      800582 :   const auto end = params.end();
     193      800582 :   const auto stride = params.stride();
     194      800582 :   const auto index = out.m_dims.index(dim);
     195      800582 :   out.m_offset += begin * m_strides[index];
     196      800582 :   if (end == -1) {
     197      115440 :     out.m_strides.erase(index);
     198      115440 :     out.m_dims.erase(dim);
     199             :   } else {
     200      685142 :     static_cast<Sizes &>(out.m_dims) = out.m_dims.slice(params);
     201      685142 :     out.m_strides[index] *= stride;
     202             :   }
     203      800582 :   return out;
     204           0 : }
     205             : 
     206        3987 : void Variable::validateSlice(const Slice &s, const Variable &data) const {
     207        3987 :   core::expect::validSlice(this->dims(), s);
     208        3987 :   if (variableFactory().has_variances(data) !=
     209        3987 :       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        3986 :   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        3985 :   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        3984 : }
     231             : 
     232        3926 : Variable &Variable::setSlice(const Slice params, const Variable &data) {
     233        3926 :   validateSlice(params, data);
     234        3926 :   copy(data, slice(params));
     235        3912 :   return *this;
     236             : }
     237             : 
     238       80972 : Variable Variable::broadcast(const Dimensions &target) const {
     239       80972 :   expect::includes(target, dims());
     240       80970 :   auto out = target.volume() == dims().volume() ? *this : as_const();
     241       80970 :   out.m_dims = target;
     242       80970 :   out.m_strides.clear();
     243      142809 :   for (const auto &d : target.labels())
     244       61839 :     out.m_strides.push_back(dims().contains(d) ? m_strides[dims().index(d)]
     245             :                                                : 0);
     246       80970 :   return out;
     247           0 : }
     248             : 
     249         247 : Variable Variable::fold(const Dim dim, const Dimensions &target) const {
     250         247 :   auto out(*this);
     251         247 :   out.m_dims = core::fold(dims(), dim, target);
     252         246 :   out.m_strides.clear();
     253         246 :   const Strides substrides(target);
     254         531 :   for (scipp::index i_in = 0; i_in < dims().ndim(); ++i_in) {
     255         285 :     if (dims().label(i_in) == dim)
     256         746 :       for (scipp::index i_target = 0; i_target < target.ndim(); ++i_target)
     257         500 :         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         492 :   return out;
     262         247 : }
     263             : 
     264        5889 : Variable Variable::transpose(const scipp::span<const Dim> order) const {
     265        5889 :   auto transposed(*this);
     266        5897 :   transposed.m_strides = core::transpose(m_strides, dims(), order);
     267        5881 :   transposed.m_dims = core::transpose(dims(), order);
     268        5881 :   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        1957 : Variable Variable::rename_dims(const std::vector<std::pair<Dim, Dim>> &names,
     276             :                                const bool fail_on_unknown) const {
     277        1957 :   auto out(*this);
     278             :   // names is a vector for predictable order
     279        1957 :   out.m_dims = out.m_dims.rename_dims(names, fail_on_unknown);
     280        1956 :   return out;
     281           1 : }
     282             : 
     283    25414926 : bool Variable::is_valid() const noexcept { return m_object.operator bool(); }
     284             : 
     285       71406 : bool Variable::is_slice() const {
     286             :   // TODO Is this condition sufficient?
     287       71406 :   return m_offset != 0 || m_dims.volume() != data().size();
     288             : }
     289             : 
     290      133826 : bool Variable::is_readonly() const noexcept { return m_readonly; }
     291             : 
     292        6754 : bool Variable::is_same(const Variable &other) const noexcept {
     293        6754 :   return std::tie(m_dims, m_strides, m_offset, m_object) ==
     294        6754 :          std::tie(other.m_dims, other.m_strides, other.m_offset,
     295       13508 :                   other.m_object);
     296             : }
     297             : 
     298      355959 : bool Variable::is_aligned() const noexcept { return m_aligned; }
     299             : 
     300      401407 : void Variable::set_aligned(bool aligned) noexcept { m_aligned = aligned; }
     301             : 
     302        3571 : void Variable::setVariances(const Variable &v) {
     303        3571 :   expect_writable();
     304        3560 :   if (is_slice())
     305          11 :     throw except::VariancesError(
     306          22 :         "Cannot add variances via sliced view of Variable.");
     307        3549 :   if (v.is_valid()) {
     308        3524 :     core::expect::equals(unit(), v.unit());
     309        3520 :     core::expect::equals(dims(), v.dims());
     310             :   }
     311        3545 :   data().setVariances(v);
     312        3540 : }
     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      938260 : Variable Variable::bin_indices() const {
     323      938260 :   auto out{*this};
     324      938260 :   out.m_object = data().bin_indices();
     325      938260 :   return out;
     326           0 : }
     327             : 
     328       69523 : Variable Variable::as_const() const {
     329       69523 :   Variable out(*this);
     330       69523 :   out.m_readonly = true;
     331       69523 :   return out;
     332             : }
     333             : 
     334     4379199 : void Variable::expect_writable() const {
     335     4379199 :   if (m_readonly)
     336          40 :     throw except::VariableError("Read-only flag is set, cannot mutate data.");
     337     4379159 : }
     338             : } // namespace scipp::variable