MDLeanEvent.h

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#pragma once
 
#include "MantidDataObjects/DllConfig.h"
#include "MantidDataObjects/MortonIndex/BitInterleaving.h"
#include "MantidDataObjects/MortonIndex/CoordinateConversion.h"
#include "MantidGeometry/MDGeometry/MDTypes.h"
#include <algorithm>
#include <cmath>
#include <cstring>
#include <numeric>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace Mantid {
namespace DataObjects {
 
template <size_t nd> class MDLeanEvent;
 
template <size_t nd> void swap(MDLeanEvent<nd> &first, MDLeanEvent<nd> &second);
 
struct EventAccessor {};
 
template <size_t nd> class MANTID_DATAOBJECTS_DLL MDLeanEvent {
public:
  using IntT = typename morton_index::IndexTypes<nd, coord_t>::IntType;
  using MortonT = typename morton_index::IndexTypes<nd, coord_t>::MortonType;
 
  template <class Accessor>
  struct AccessFor {
    static std::enable_if_t<std::is_same<EventAccessor, typename Accessor::EventAccessType>::value>
    convertToCoordinates(MDLeanEvent<nd> &event, const morton_index::MDSpaceBounds<nd> &space) {
      event.convertToCoordinates(space);
    }
    static std::enable_if_t<std::is_same<EventAccessor, typename Accessor::EventAccessType>::value>
    convertToIndex(MDLeanEvent<nd> &event, const morton_index::MDSpaceBounds<nd> &space) {
      event.convertToIndex(space);
    }
    static
        typename std::enable_if<std::is_same<EventAccessor, typename Accessor::EventAccessType>::value, MortonT>::type
        getIndex(const MDLeanEvent<nd> &event) {
      return event.index;
    }
  };
  template <class Accessor> friend struct AccessFor;
 
protected:
  float signal;
 
  float errorSquared;
 
#pragma pack(push, 1)
  union {
    coord_t center[nd];
    MortonT index = 0;
  };
#pragma pack(pop)
 
private:
  void convertToIndex(const morton_index::MDSpaceBounds<nd> &space) {
    index = morton_index::coordinatesToIndex<nd, IntT, MortonT>(center, space);
  }
 
  void convertToCoordinates(const morton_index::MDSpaceBounds<nd> &space) {
    auto coords = morton_index::indexToCoordinates<nd, IntT, MortonT>(index, space);
    for (unsigned i = 0; i < nd; ++i)
      center[i] = coords[i];
  }
 
public:
  /* Will be keeping functions inline for (possible?) performance improvements
   */
 
  // Enum to flag this templated type as NOT a full md event type.
  enum { is_full_mdevent = false };
 
  //---------------------------------------------------------------------------------------------
  MDLeanEvent() : signal(1.0), errorSquared(1.0) {}
 
  //---------------------------------------------------------------------------------------------
  MDLeanEvent(const double signal, const double errorSquared)
      : signal(float(signal)), errorSquared(float(errorSquared)) {}
 
  //---------------------------------------------------------------------------------------------
  MDLeanEvent(const float signal, const float errorSquared) : signal(signal), errorSquared(errorSquared) {}
 
  //---------------------------------------------------------------------------------------------
  MDLeanEvent(const float signal, const float errorSquared, const coord_t *centers)
      : signal(signal), errorSquared(errorSquared) {
    for (size_t i = 0; i < nd; i++)
      center[i] = centers[i];
  }
  //---------------------------------------------------------------------------------------------
  MDLeanEvent(const double signal, const double errorSquared, const coord_t *centers)
      : signal(float(signal)), errorSquared(float(errorSquared)) {
    for (size_t i = 0; i < nd; i++)
      center[i] = centers[i];
  }
 
#ifdef COORDT_IS_FLOAT
  //---------------------------------------------------------------------------------------------
  MDLeanEvent(const float signal, const float errorSquared, const double *centers)
      : signal(signal), errorSquared(errorSquared) {
    for (size_t i = 0; i < nd; i++)
      center[i] = static_cast<coord_t>(centers[i]);
  }
#endif
  //---------------------------------------------------------------------------------------------
  MDLeanEvent(const MDLeanEvent &rhs) : signal(rhs.signal), errorSquared(rhs.errorSquared) {
    for (size_t i = 0; i < nd; i++)
      center[i] = rhs.center[i];
  }
 
  //---------------------------------------------------------------------------------------------
  friend void swap<nd>(MDLeanEvent &first, MDLeanEvent &second);
 
  MDLeanEvent &operator=(MDLeanEvent other) {
    swap(*this, other);
    return *this;
  }
 
  coord_t getCenter(const size_t n) const { return center[n]; }
 
  //---------------------------------------------------------------------------------------------
  const coord_t *getCenter() const { return center; }
 
  //---------------------------------------------------------------------------------------------
  coord_t *getCenterNonConst() { return center; }
 
  //---------------------------------------------------------------------------------------------
  void setCenter(const size_t n, const coord_t value) { center[n] = value; }
 
#ifdef COORDT_IS_FLOAT
  //---------------------------------------------------------------------------------------------
  void setCenter(const size_t n, const double value) { center[n] = static_cast<coord_t>(value); }
#endif
 
  //---------------------------------------------------------------------------------------------
  void setCoords(const coord_t *centers) {
    for (size_t i = 0; i < nd; i++)
      center[i] = centers[i];
  }
 
  //---------------------------------------------------------------------------------------------
  size_t getNumDims() const { return nd; }
 
  //---------------------------------------------------------------------------------------------
  float getSignal() const { return signal; }
 
  //---------------------------------------------------------------------------------------------
  float getErrorSquared() const { return errorSquared; }
 
  //---------------------------------------------------------------------------------------------
  float getError() const { return float(sqrt(errorSquared)); }
 
  //---------------------------------------------------------------------------------------------
  void setSignal(const float newSignal) { signal = newSignal; }
 
  //---------------------------------------------------------------------------------------------
  void setErrorSquared(const float newerrorSquared) { errorSquared = newerrorSquared; }
 
  //---------------------------------------------------------------------------------------------
  static std::string getTypeName() { return "MDLeanEvent"; }
 
  //---------------------------------------------------------------------------------------------
  uint16_t getExpInfoIndex() const { return 0; }
 
  //---------------------------------------------------------------------------------------------
  uint16_t getGoniometerIndex() const { return 0; }
 
  //---------------------------------------------------------------------------------------------
  int32_t getDetectorID() const { return 0; }
 
  /* static method used to convert vector of lean events into vector of their
   coordinates & signal and error
   @param events    -- vector of events
   @return data     -- vector of events coordinates, their signal and error
   casted to coord_t type
   @return ncols    -- the number of colunts  in the data (it is nd+2 here but
   may be different for other data types)
   @return totalSignal -- total signal in the vector of events
   @return totalErr   -- total error corresponting to the vector of events
  */
  static inline void eventsToData(const std::vector<MDLeanEvent<nd>> &mdLeanEvents, std::vector<coord_t> &data,
                                  size_t &ncols, double &totalSignal, double &totalErrSq) {
    ncols = nd + 2;
    size_t nEvents = mdLeanEvents.size();
    data.resize(nEvents * ncols);
 
    totalSignal = 0;
    totalErrSq = 0;
 
    size_t dataIndex(0);
    for (const MDLeanEvent<nd> &mdLeanEvent : mdLeanEvents) {
      data[dataIndex++] = static_cast<coord_t>(mdLeanEvent.signal);
      data[dataIndex++] = static_cast<coord_t>(mdLeanEvent.errorSquared);
      for (size_t d = 0; d < nd; d++)
        data[dataIndex++] = mdLeanEvent.center[d];
      // Track the total signal
      totalSignal += signal_t(mdLeanEvent.signal);
      totalErrSq += signal_t(mdLeanEvent.errorSquared);
    }
  }
 
  /* static method used to convert vector of data into vector of lean events
   @return coord    -- vector of events coordinates, their signal and error
   casted to coord_t type
   @param events    -- vector of events
   @param reserveMemory -- reserve memory for events copying. Set to false if
   one wants to add new events to the existing one.
  */
  static inline void dataToEvents(const std::vector<coord_t> &coord, std::vector<MDLeanEvent<nd>> &events,
                                  bool reserveMemory = true) {
    // Number of columns = number of dimensions + 2 (signal/error)
    size_t numColumns = (nd + 2);
    size_t numEvents = coord.size() / numColumns;
    if (numEvents * numColumns != coord.size())
      throw(std::invalid_argument("wrong input array of data to convert to "
                                  "lean events, suspected column data for "
                                  "different dimensions/(type of) events "));
 
    if (reserveMemory) // Reserve the amount of space needed. Significant speed
                       // up (~30% thanks to this)
    {
      events.clear();
      events.reserve(numEvents);
    }
    for (size_t i = 0; i < numEvents; i++) {
      // Index into the data array
      size_t ii = i * numColumns;
 
      // Point directly into the data block for the centers.
      const coord_t *centers = &(coord[ii + 2]);
 
      // Create the event with signal, error squared, and the centers
      events.emplace_back(static_cast<signal_t>(coord[ii]), static_cast<signal_t>(coord[ii + 1]), centers);
    }
  }
};
 
template <size_t ND> void swap(MDLeanEvent<ND> &first, MDLeanEvent<ND> &second) {
  std::swap(first.signal, second.signal);
  std::swap(first.errorSquared, second.errorSquared);
  // Index always of the same size as center
  std::swap(first.index, second.index);
}
 
} // namespace DataObjects
} // namespace Mantid