MDHistoWorkspace.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 "MantidAPI/IMDHistoWorkspace.h"
#include "MantidAPI/IMDIterator.h"
#include "MantidAPI/IMDWorkspace.h"
#include "MantidAPI/MDGeometry.h"
#include "MantidDataObjects/DllConfig.h"
#include "MantidDataObjects/WorkspaceSingleValue.h"
#include "MantidGeometry/MDGeometry/IMDDimension.h"
#include "MantidGeometry/MDGeometry/MDHistoDimension.h"
#include "MantidGeometry/MDGeometry/MDImplicitFunction.h"
#include "MantidKernel/Exception.h"
 
namespace Mantid {
namespace DataObjects {
 
class MANTID_DATAOBJECTS_DLL MDHistoWorkspace : public API::IMDHistoWorkspace {
public:
  MDHistoWorkspace(Mantid::Geometry::MDHistoDimension_sptr dimX,
                   Mantid::Geometry::MDHistoDimension_sptr dimY = Mantid::Geometry::MDHistoDimension_sptr(),
                   Mantid::Geometry::MDHistoDimension_sptr dimZ = Mantid::Geometry::MDHistoDimension_sptr(),
                   Mantid::Geometry::MDHistoDimension_sptr dimT = Mantid::Geometry::MDHistoDimension_sptr(),
                   Mantid::API::MDNormalization displayNormalization = Mantid::API::NoNormalization);
 
  MDHistoWorkspace(std::vector<Mantid::Geometry::MDHistoDimension_sptr> &dimensions,
                   Mantid::API::MDNormalization displayNormalization = Mantid::API::NoNormalization);
  MDHistoWorkspace(std::vector<Mantid::Geometry::IMDDimension_sptr> &dimensions,
                   Mantid::API::MDNormalization displayNormalization = Mantid::API::NoNormalization);
  MDHistoWorkspace &operator=(const MDHistoWorkspace &other) = delete;
 
  std::unique_ptr<MDHistoWorkspace> clone() const { return std::unique_ptr<MDHistoWorkspace>(doClone()); }
 
  std::unique_ptr<MDHistoWorkspace> cloneEmpty() const { return std::unique_ptr<MDHistoWorkspace>(doCloneEmpty()); }
 
  void init(std::vector<Mantid::Geometry::MDHistoDimension_sptr> &dimensions);
  void init(std::vector<Mantid::Geometry::IMDDimension_sptr> &dimensions);
 
  void cacheValues();
 
  const std::string id() const override { return "MDHistoWorkspace"; }
 
  size_t getMemorySize() const override;
 
  uint64_t getNPoints() const override { return m_length; }
  uint64_t getNEvents() const override;
  std::vector<std::unique_ptr<Mantid::API::IMDIterator>>
  createIterators(size_t suggestedNumCores = 1,
                  Mantid::Geometry::MDImplicitFunction *function = nullptr) const override;
 
  LinePlot getLinePlot(const Mantid::Kernel::VMD &start, const Mantid::Kernel::VMD &end,
                       Mantid::API::MDNormalization normalize) const override;
 
  LinePlot getLineData(const Mantid::Kernel::VMD &start, const Mantid::Kernel::VMD &end,
                       Mantid::API::MDNormalization normalize) const override;
 
  void checkWorkspaceSize(const MDHistoWorkspace &other, const std::string &operation);
 
  // --------------------------------------------------------------------------------------------
  MDHistoWorkspace &operator+=(const MDHistoWorkspace &b);
  void add(const MDHistoWorkspace &b);
  void add(const signal_t signal, const signal_t error);
 
  MDHistoWorkspace &operator-=(const MDHistoWorkspace &b);
  void subtract(const MDHistoWorkspace &b);
  void subtract(const signal_t signal, const signal_t error);
 
  MDHistoWorkspace &operator*=(const MDHistoWorkspace &b_ws);
  void multiply(const MDHistoWorkspace &b_ws);
  void multiply(const signal_t signal, const signal_t error);
 
  MDHistoWorkspace &operator/=(const MDHistoWorkspace &b_ws);
  void divide(const MDHistoWorkspace &b_ws);
  void divide(const signal_t signal, const signal_t error);
 
  void log(double filler = 0.0);
  void log10(double filler = 0.0);
  void exp();
  void power(double exponent);
 
  // --------------------------------------------------------------------------------------------
 
  MDHistoWorkspace &operator&=(const MDHistoWorkspace &b);
  MDHistoWorkspace &operator|=(const MDHistoWorkspace &b);
  MDHistoWorkspace &operator^=(const MDHistoWorkspace &b);
  void operatorNot();
 
  void lessThan(const MDHistoWorkspace &b);
  void lessThan(const signal_t signal);
  void greaterThan(const MDHistoWorkspace &b);
  void greaterThan(const signal_t signal);
  void equalTo(const MDHistoWorkspace &b, const signal_t tolerance = 1e-5);
  void equalTo(const signal_t signal, const signal_t tolerance = 1e-5);
 
  void setUsingMask(const MDHistoWorkspace &mask, const MDHistoWorkspace &values);
  void setUsingMask(const MDHistoWorkspace &mask, const signal_t signal, const signal_t error);
 
  // --------------------------------------------------------------------------------------------
  const size_t *getIndexMultiplier() const { return indexMultiplier.data(); }
 
  const signal_t *getSignalArray() const override { return m_signals.data(); }
 
  coord_t getInverseVolume() const override { return m_inverseVolume; }
 
  const signal_t *getErrorSquaredArray() const override { return m_errorsSquared.data(); }
 
  const signal_t *getNumEventsArray() const override { return m_numEvents.data(); }
 
  const bool *getMaskArray() const { return m_masks.get(); }
 
  const coord_t *getBinWidths() const { return m_boxLength.data(); }
 
  signal_t *mutableSignalArray() override { return m_signals.data(); }
 
  signal_t *mutableErrorSquaredArray() override { return m_errorsSquared.data(); }
 
  signal_t *mutableNumEventsArray() override { return m_numEvents.data(); }
 
  bool *mutableMaskArray() { return m_masks.get(); }
 
  Kernel::SpecialCoordinateSystem getSpecialCoordinateSystem() const override;
 
  void setCoordinateSystem(const Kernel::SpecialCoordinateSystem coordinateSystem) override;
 
  void setTo(signal_t signal, signal_t errorSquared, signal_t numEvents) override;
 
  void applyImplicitFunction(Mantid::Geometry::MDImplicitFunction *function, signal_t signal, signal_t errorSquared);
 
  std::unique_ptr<coord_t[]> getVertexesArray(size_t linearIndex, size_t &numVertices) const;
 
  Kernel::VMD getCenter(size_t linearIndex) const override;
 
  signal_t getSignalAtCoord(const coord_t *coords, const Mantid::API::MDNormalization &normalization) const override;
 
  // or 0 if masked
  signal_t getSignalWithMaskAtCoord(const coord_t *coords,
                                    const Mantid::API::MDNormalization &normalization) const override;
 
  void setSignalAt(size_t index, signal_t value) override { m_signals[index] = value; }
 
  void setErrorSquaredAt(size_t index, signal_t value) override { m_errorsSquared[index] = value; }
 
  void setNumEventsAt(size_t index, signal_t value) { m_numEvents[index] = value; }
 
  signal_t getNumEventsAt(size_t index) const { return m_numEvents[index]; }
 
  signal_t getErrorAt(size_t index) const override { return std::sqrt(m_errorsSquared[index]); }
 
  signal_t getErrorAt(size_t index1, size_t index2) const override {
    return std::sqrt(m_errorsSquared[index1 + indexMultiplier[0] * index2]);
  }
 
  signal_t getErrorAt(size_t index1, size_t index2, size_t index3) const override {
    return std::sqrt(m_errorsSquared[index1 + indexMultiplier[0] * index2 + indexMultiplier[1] * index3]);
  }
 
  signal_t getErrorAt(size_t index1, size_t index2, size_t index3, size_t index4) const override {
    return std::sqrt(m_errorsSquared[index1 + indexMultiplier[0] * index2 + indexMultiplier[1] * index3 +
                                     indexMultiplier[2] * index4]);
  }
 
  bool getIsMaskedAt(size_t index) const { return m_masks[index]; }
 
  signal_t getSignalAt(size_t index) const override { return m_signals[index]; }
 
  signal_t getSignalAt(size_t index1, size_t index2) const override {
    return m_signals[index1 + indexMultiplier[0] * index2];
  }
 
  signal_t getSignalAt(size_t index1, size_t index2, size_t index3) const override {
    return m_signals[index1 + indexMultiplier[0] * index2 + indexMultiplier[1] * index3];
  }
 
  signal_t getSignalAt(size_t index1, size_t index2, size_t index3, size_t index4) const override {
    return m_signals[index1 + indexMultiplier[0] * index2 + indexMultiplier[1] * index3 + indexMultiplier[2] * index4];
  }
 
  signal_t getSignalNormalizedAt(size_t index) const override { return m_signals[index] * m_inverseVolume; }
 
  signal_t getSignalNormalizedAt(size_t index1, size_t index2) const override {
    return m_signals[index1 + indexMultiplier[0] * index2] * m_inverseVolume;
  }
 
  signal_t getSignalNormalizedAt(size_t index1, size_t index2, size_t index3) const override {
    return m_signals[index1 + indexMultiplier[0] * index2 + indexMultiplier[1] * index3] * m_inverseVolume;
  }
 
  signal_t getSignalNormalizedAt(size_t index1, size_t index2, size_t index3, size_t index4) const override {
    return m_signals[index1 + indexMultiplier[0] * index2 + indexMultiplier[1] * index3 + indexMultiplier[2] * index4] *
           m_inverseVolume;
  }
 
  signal_t getErrorNormalizedAt(size_t index) const override {
    return std::sqrt(m_errorsSquared[index]) * m_inverseVolume;
  }
 
  signal_t getErrorNormalizedAt(size_t index1, size_t index2) const override {
    return getErrorAt(index1, index2) * m_inverseVolume;
  }
 
  signal_t getErrorNormalizedAt(size_t index1, size_t index2, size_t index3) const override {
    return getErrorAt(index1, index2, index3) * m_inverseVolume;
  }
 
  signal_t getErrorNormalizedAt(size_t index1, size_t index2, size_t index3, size_t index4) const override {
    return getErrorAt(index1, index2, index3, index4) * m_inverseVolume;
  }
 
  //---------------------------------------------------------------------------------------------
  signal_t &errorSquaredAt(size_t index) override {
    if (index < m_length)
      return m_errorsSquared[index];
    else
      throw std::invalid_argument("MDHistoWorkspace::array index out of range");
  }
 
  signal_t &signalAt(size_t index) override {
    if (index < m_length)
      return m_signals[index];
    else
      throw std::invalid_argument("MDHistoWorkspace::array index out of range");
  }
 
  //---------------------------------------------------------------------------------------------
  size_t getLinearIndex(size_t index1, size_t index2) const override {
    if (this->getNumDims() != 2)
      throw std::runtime_error("Workspace does not have 2 dimensions!");
    return index1 + indexMultiplier[0] * index2;
  }
 
  size_t getLinearIndex(size_t index1, size_t index2, size_t index3) const override {
    if (this->getNumDims() != 3)
      throw std::runtime_error("Workspace does not have 3 dimensions!");
    return index1 + indexMultiplier[0] * index2 + indexMultiplier[1] * index3;
  }
 
  size_t getLinearIndex(size_t index1, size_t index2, size_t index3, size_t index4) const override {
    if (this->getNumDims() != 4)
      throw std::runtime_error("Workspace does not have 4 dimensions!");
    return index1 + indexMultiplier[0] * index2 + indexMultiplier[1] * index3 + indexMultiplier[2] * index4;
  }
 
  size_t getLinearIndex(size_t *index) const {
    size_t nd = this->getNumDims();
    size_t out = index[0];
    for (size_t d = 1; d < nd; d++)
      out += indexMultiplier[d - 1] * index[d];
    return out;
  }
 
  size_t getLinearIndexAtCoord(const coord_t *coords) const;
 
  signal_t &operator[](const size_t &index) override {
    if (index < m_length)
      return m_signals[index];
    else
      throw std::invalid_argument("MDHistoWorkspace::array index out of range");
  }
 
  virtual std::vector<signal_t> getSignalDataVector() const;
  virtual std::vector<signal_t> getErrorDataVector() const;
 
  void setMDMasking(std::unique_ptr<Mantid::Geometry::MDImplicitFunction> maskingRegion) override;
  void setMDMaskAt(const size_t &index, bool mask);
 
  void clearMDMasking() override;
  uint64_t sumNContribEvents() const;
  void updateSum() { m_nEventsContributed = sumNContribEvents(); }
 
  static size_t sizeOfElement();
 
  Mantid::API::MDNormalization displayNormalization() const override;
 
  Mantid::API::MDNormalization displayNormalizationHisto() const override;
 
  void setDisplayNormalization(const Mantid::API::MDNormalization &preferredNormalization) override;
 
  bool isMDHistoWorkspace() const override { return true; }
 
private:
  MDHistoWorkspace *doClone() const override { return new MDHistoWorkspace(*this); }
 
  MDHistoWorkspace *doCloneEmpty() const override { return new MDHistoWorkspace(nullptr); }
 
  void makeSingleBinWithNaN(std::vector<coord_t> &x, std::vector<signal_t> &y, std::vector<signal_t> &e) const;
 
  void initVertexesArray();
 
  size_t numDimensions;
 
  std::vector<signal_t> m_signals;
 
  std::vector<signal_t> m_errorsSquared;
 
  std::vector<signal_t> m_numEvents;
 
  size_t m_length;
 
  std::vector<size_t> indexMultiplier;
  std::vector<size_t> m_indexMaker;
  std::vector<size_t> m_indexMax;
 
  coord_t m_inverseVolume;
 
  std::vector<coord_t> m_vertexesArray;
 
  std::vector<coord_t> m_boxLength;
 
  std::vector<coord_t> m_origin;
  mutable uint64_t m_nEventsContributed;
 
  Kernel::SpecialCoordinateSystem m_coordSystem;
 
  Mantid::API::MDNormalization m_displayNormalization;
 
  // Get ordered list of boundaries in position-along-the-line coordinates
  std::set<coord_t> getBinBoundariesOnLine(const Kernel::VMD &start, const Kernel::VMD &end, size_t nd,
                                           const Kernel::VMD &dir, coord_t length) const;
 
  signal_t getNormalizationFactor(const API::MDNormalization &normalize, size_t linearIndex) const;
 
protected:
  LinePlot getLinePoints(const Mantid::Kernel::VMD &start, const Mantid::Kernel::VMD &end,
                         Mantid::API::MDNormalization normalize, const bool bin_centres) const;
 
  MDHistoWorkspace(const MDHistoWorkspace &other);
 
  std::unique_ptr<bool[]> m_masks;
};
 
using MDHistoWorkspace_sptr = std::shared_ptr<MDHistoWorkspace>;
 
using MDHistoWorkspace_const_sptr = std::shared_ptr<const MDHistoWorkspace>;
 
} // namespace DataObjects
} // namespace Mantid