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 "scipp/variable/operations.h"
       6             : 
       7             : #include "scipp/dataset/counts.h"
       8             : #include "scipp/dataset/dataset.h"
       9             : #include "scipp/dataset/except.h"
      10             : 
      11             : namespace scipp::dataset {
      12             : 
      13             : namespace counts {
      14             : 
      15           6 : std::vector<Variable> getBinWidths(const Coords &c,
      16             :                                    const std::vector<Dim> &dims) {
      17           6 :   std::vector<Variable> binWidths;
      18          12 :   for (const auto &dim : dims) {
      19           6 :     const auto &coord = c[dim];
      20           6 :     if (coord.unit() == units::dimensionless)
      21           0 :       throw std::runtime_error("Dimensionless axis cannot be used for "
      22           0 :                                "conversion from or to density");
      23           6 :     binWidths.emplace_back(coord.slice({dim, 1, coord.dims()[dim]}) -
      24          12 :                            coord.slice({dim, 0, coord.dims()[dim] - 1}));
      25             :   }
      26           6 :   return binWidths;
      27           0 : }
      28             : 
      29           3 : void toDensity(DataArray &data, const std::vector<Variable> &binWidths) {
      30           3 :   if (data.unit().isCounts()) {
      31           6 :     for (const auto &binWidth : binWidths)
      32           3 :       data /= binWidth;
      33           0 :   } else if (data.unit().isCountDensity()) {
      34             :     // This error implies that conversion to multi-dimensional densities
      35             :     // must be done in one step, e.g., counts -> counts/(m*m*s). We cannot
      36             :     // do counts -> counts/m -> counts/(m*m) -> counts/(m*m*s) since the
      37             :     // unit-based distinction between counts and counts-density cannot tell
      38             :     // apart different length dimensions such as X and Y, so we would not be
      39             :     // able to prevent converting to density using Dim::X twice.
      40           0 :     throw std::runtime_error("Cannot convert counts-variable to density, "
      41             :                              "it looks like it has already been "
      42           0 :                              "converted.");
      43             :   }
      44             :   // No `else`, variables that do not contain a `counts` factor are left
      45             :   // unchanged.
      46           3 : }
      47             : 
      48           2 : Dataset toDensity(const Dataset &d, const Dim dim) {
      49           4 :   return toDensity(d, std::vector<Dim>{dim});
      50             : }
      51             : 
      52           2 : Dataset toDensity(const Dataset &d, const std::vector<Dim> &dims) {
      53           2 :   const auto binWidths = getBinWidths(d.coords(), dims);
      54           4 :   for (auto &&data : d)
      55           2 :     toDensity(data, binWidths);
      56           4 :   return d;
      57           2 : }
      58             : 
      59           1 : DataArray toDensity(const DataArray &a, const Dim dim) {
      60           2 :   return toDensity(a, std::vector<Dim>{dim});
      61             : }
      62             : 
      63           1 : DataArray toDensity(const DataArray &a, const std::vector<Dim> &dims) {
      64           1 :   const auto binWidths = getBinWidths(a.coords(), dims);
      65           1 :   auto out = copy(a);
      66           1 :   toDensity(out, binWidths);
      67           2 :   return out;
      68           1 : }
      69             : 
      70           3 : void fromDensity(DataArray &data, const std::vector<Variable> &binWidths) {
      71           3 :   if (data.unit().isCounts()) {
      72             :     // Do nothing, but do not fail either.
      73           3 :   } else if (data.unit().isCountDensity()) {
      74           6 :     for (const auto &binWidth : binWidths)
      75           3 :       data *= binWidth;
      76             :   }
      77           3 : }
      78             : 
      79           2 : Dataset fromDensity(const Dataset &d, const Dim dim) {
      80           4 :   return fromDensity(d, std::vector<Dim>{dim});
      81             : }
      82             : 
      83           2 : Dataset fromDensity(const Dataset &d, const std::vector<Dim> &dims) {
      84           2 :   const auto binWidths = getBinWidths(d.coords(), dims);
      85           4 :   for (auto &&data : d)
      86           2 :     fromDensity(data, binWidths);
      87           4 :   return d;
      88           2 : }
      89             : 
      90           1 : DataArray fromDensity(const DataArray &a, const Dim dim) {
      91           2 :   return fromDensity(a, std::vector<Dim>{dim});
      92             : }
      93             : 
      94           1 : DataArray fromDensity(const DataArray &a, const std::vector<Dim> &dims) {
      95           1 :   const auto binWidths = getBinWidths(a.coords(), dims);
      96           1 :   auto out = copy(a);
      97           1 :   fromDensity(out, binWidths);
      98           2 :   return out;
      99           1 : }
     100             : 
     101             : } // namespace counts
     102             : } // namespace scipp::dataset