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             : #pragma once
       6             : #include <optional>
       7             : 
       8             : #include "scipp/common/initialization.h"
       9             : #include "scipp/common/numeric.h"
      10             : #include "scipp/core/dimensions.h"
      11             : #include "scipp/core/eigen.h"
      12             : #include "scipp/core/element_array_view.h"
      13             : #include "scipp/core/except.h"
      14             : #include "scipp/units/unit.h"
      15             : #include "scipp/variable/except.h"
      16             : #include "scipp/variable/transform.h"
      17             : #include "scipp/variable/variable_concept.h"
      18             : 
      19             : namespace scipp::variable {
      20             : 
      21             : template <class T> struct is_span : std::false_type {};
      22             : template <class T> struct is_span<scipp::span<T>> : std::true_type {};
      23             : template <class T> inline constexpr bool is_span_v = is_span<T>::value;
      24             : 
      25             : template <class T1, class T2>
      26      549703 : bool equals_impl(const T1 &view1, const T2 &view2) {
      27             :   // TODO Use optimizations in case of contiguous views (instead of slower
      28             :   // ElementArrayView iteration). Add multi threading?
      29             :   if constexpr (is_span_v<typename T1::value_type>)
      30           0 :     return std::equal(
      31             :         view1.begin(), view1.end(), view2.begin(), view2.end(),
      32           0 :         [](const auto &a, const auto &b) { return equals_impl(a, b); });
      33             :   else
      34      549703 :     return std::equal(view1.begin(), view1.end(), view2.begin(), view2.end());
      35             : }
      36             : 
      37        7734 : template <class T> bool equals_nan(const T &a, const T &b) {
      38             :   if constexpr (std::is_floating_point_v<T>) {
      39             :     using numeric::isnan;
      40        6845 :     if (isnan(a) && isnan(b))
      41          63 :       return true;
      42             :   }
      43        7671 :   return a == b;
      44             : }
      45             : 
      46             : template <class T1, class T2>
      47        2225 : bool equals_nan_impl(const T1 &view1, const T2 &view2) {
      48             :   // TODO Use optimizations in case of contiguous views (instead of slower
      49             :   // ElementArrayView iteration). Add multi threading?
      50             :   if constexpr (is_span_v<typename T1::value_type>)
      51           0 :     return std::equal(
      52             :         view1.begin(), view1.end(), view2.begin(), view2.end(),
      53           0 :         [](const auto &a, const auto &b) { return equals_nan_impl(a, b); });
      54             :   else
      55        2225 :     return std::equal(
      56             :         view1.begin(), view1.end(), view2.begin(), view2.end(),
      57       10058 :         [](const auto &x, const auto &y) { return equals_nan(x, y); });
      58             : }
      59             : 
      60             : /// Implementation of VariableConcept that holds an array with element type T.
      61             : template <class T> class ElementArrayModel : public VariableConcept {
      62             : public:
      63             :   static_assert(!core::is_structured(core::template dtype<T>));
      64             :   using value_type = T;
      65             : 
      66             :   ElementArrayModel(const scipp::index size, const units::Unit &unit,
      67             :                     element_array<T> model,
      68             :                     std::optional<element_array<T>> variances = std::nullopt);
      69             : 
      70          15 :   static DType static_dtype() noexcept { return scipp::dtype<T>; }
      71    53571942 :   DType dtype() const noexcept override { return scipp::dtype<T>; }
      72      390556 :   scipp::index size() const override { return m_values.size(); }
      73             : 
      74             :   VariableConceptHandle
      75             :   makeDefaultFromParent(const scipp::index size) const override;
      76             : 
      77             :   VariableConceptHandle
      78           0 :   makeDefaultFromParent(const Variable &shape) const override {
      79           0 :     return makeDefaultFromParent(shape.dims().volume());
      80             :   }
      81             : 
      82             :   bool equals(const Variable &a, const Variable &b) const override;
      83             :   bool equals_nan(const Variable &a, const Variable &b) const override;
      84             :   void copy(const Variable &src, Variable &dest) const override;
      85             :   void copy(const Variable &src, Variable &&dest) const override;
      86             :   void assign(const VariableConcept &other) override;
      87             : 
      88             :   void setVariances(const Variable &variances) override;
      89             : 
      90             :   VariableConceptHandle clone() const override;
      91             : 
      92     8233331 :   bool has_variances() const noexcept override {
      93     8233331 :     return m_variances.has_value();
      94             :   }
      95             : 
      96     3736026 :   auto values(const core::ElementArrayViewParams &base) const {
      97     3736026 :     return ElementArrayView(base, m_values.data());
      98             :   }
      99     2619361 :   auto values(const core::ElementArrayViewParams &base) {
     100     2619361 :     return ElementArrayView(base, m_values.data());
     101             :   }
     102      318568 :   auto variances(const core::ElementArrayViewParams &base) const {
     103      318568 :     expect_has_variances();
     104      318567 :     return ElementArrayView(base, m_variances->data());
     105             :   }
     106       55219 :   auto variances(const core::ElementArrayViewParams &base) {
     107       55219 :     expect_has_variances();
     108       55217 :     return ElementArrayView(base, m_variances->data());
     109             :   }
     110             : 
     111        1486 :   scipp::index dtype_size() const override { return sizeof(T); }
     112        1517 :   scipp::index object_size() const override { return sizeof(*this); }
     113           0 :   const VariableConceptHandle &bin_indices() const override {
     114           0 :     throw except::TypeError("This data type does not have bin indices.");
     115             :   }
     116             : 
     117      681889 :   scipp::span<const T> values() const {
     118      681889 :     return {m_values.data(), m_values.data() + m_values.size()};
     119             :   }
     120             : 
     121      419253 :   scipp::span<T> values() {
     122      419253 :     return {m_values.data(), m_values.data() + m_values.size()};
     123             :   }
     124             : 
     125             : private:
     126      373787 :   void expect_has_variances() const {
     127      373787 :     if (!has_variances())
     128           3 :       throw except::VariancesError("Variable does not have variances.");
     129      373784 :   }
     130             :   element_array<T> m_values;
     131             :   std::optional<element_array<T>> m_variances;
     132             : };
     133             : 
     134             : namespace {
     135   110452144 : template <class T> auto copy(const T &x) { return x; }
     136   110452193 : constexpr auto do_copy = [](auto &a, const auto &b) { a = copy(b); };
     137             : } // namespace
     138             : 
     139             : /// Helper for implementing Variable(View) copy operations.
     140             : ///
     141             : /// This method is using virtual dispatch as a trick to obtain T, such that
     142             : /// transform can be called with any T.
     143             : template <class T>
     144      332456 : void ElementArrayModel<T>::copy(const Variable &src, Variable &dest) const {
     145      332456 :   transform_in_place<T>(
     146             :       dest, src,
     147             :       overloaded{core::transform_flags::expect_in_variance_if_out_variance,
     148             :                  core::transform_flags::force_variance_broadcast, do_copy},
     149             :       "copy");
     150      332431 : }
     151             : template <class T>
     152        5238 : void ElementArrayModel<T>::copy(const Variable &src, Variable &&dest) const {
     153        5238 :   copy(src, dest);
     154        5234 : }
     155             : 
     156             : } // namespace scipp::variable