IMDHistoWorkspace.cpp

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 +
#include "MantidAPI/IMDHistoWorkspace.h"
#include "MantidGeometry/MDGeometry/IMDDimension.h"
#include "MantidPythonInterface/api/RegisterWorkspacePtrToPython.h"
#include "MantidPythonInterface/core/Converters/NDArrayTypeIndex.h"
#include "MantidPythonInterface/core/GetPointer.h"
#include "MantidPythonInterface/core/NDArray.h"
 
#include <boost/python/class.hpp>
#include <boost/python/copy_non_const_reference.hpp>
#define PY_ARRAY_UNIQUE_SYMBOL API_ARRAY_API
#define NO_IMPORT_ARRAY
#include <numpy/arrayobject.h>
 
#include <algorithm>
#include <stdexcept>
#include <vector>
 
using namespace Mantid::API;
using Mantid::PythonInterface::NDArray;
using Mantid::PythonInterface::Registry::RegisterWorkspacePtrToPython;
namespace Converters = Mantid::PythonInterface::Converters;
using namespace boost::python;
 
GET_POINTER_SPECIALIZATION(IMDHistoWorkspace)
 
namespace Mantid::PythonInterface::Converters {
extern template int NDArrayTypeIndex<float>::typenum;
extern template int NDArrayTypeIndex<double>::typenum;
} // namespace Mantid::PythonInterface::Converters
 
namespace {
PyObject *WrapReadOnlyNumpyFArray(const Mantid::signal_t *arr, std::vector<Py_intptr_t> dims) {
  int datatype = Converters::NDArrayTypeIndex<Mantid::signal_t>::typenum;
#if NPY_API_VERSION >= 0x00000007 //(1.7)
  auto *nparray = reinterpret_cast<PyArrayObject *>(
      PyArray_New(&PyArray_Type, static_cast<int>(dims.size()), &dims[0], datatype, nullptr,
                  static_cast<void *>(const_cast<double *>(arr)), 0, NPY_ARRAY_FARRAY, nullptr));
  PyArray_CLEARFLAGS(nparray, NPY_ARRAY_WRITEABLE);
#else
  PyArrayObject *nparray = reinterpret_cast<PyArrayObject *>(
      PyArray_New(&PyArray_Type, static_cast<int>(dims.size()), &dims[0], datatype, nullptr,
                  static_cast<void *>(const_cast<double *>(arr)), 0, NPY_FARRAY, nullptr));
  nparray->flags &= ~NPY_WRITEABLE;
#endif
  return reinterpret_cast<PyObject *>(nparray);
}
 
std::vector<Py_intptr_t> countDimensions(const IMDHistoWorkspace &self) {
  size_t ndims = self.getNumDims();
  std::vector<size_t> nd;
  nd.reserve(ndims);
 
  // invert dimensions in C way, e.g. slowest changing ndim goes first
  for (size_t i = 0; i < ndims; ++i) {
    nd.emplace_back(self.getDimension(i)->getNBins());
  }
 
  ndims = nd.size();
  std::vector<Py_intptr_t> dims(ndims);
  for (size_t i = 0; i < ndims; ++i)
    dims[i] = static_cast<Py_intptr_t>(nd[i]);
 
  if (dims.empty()) {
    throw std::runtime_error("Workspace has zero dimensions!");
  } else {
    return dims;
  }
}
 
PyObject *getSignalArrayAsNumpyArray(const IMDHistoWorkspace &self) {
  auto dims = countDimensions(self);
  return WrapReadOnlyNumpyFArray(self.getSignalArray(), dims);
}
 
PyObject *getErrorSquaredArrayAsNumpyArray(const IMDHistoWorkspace &self) {
  auto dims = countDimensions(self);
  return WrapReadOnlyNumpyFArray(self.getErrorSquaredArray(), dims);
}
 
PyObject *getNumEventsArrayAsNumpyArray(const IMDHistoWorkspace &self) {
  auto dims = countDimensions(self);
  return WrapReadOnlyNumpyFArray(self.getNumEventsArray(), dims);
}
 
void throwIfSizeIncorrect(const IMDHistoWorkspace &self, const NDArray &signal, const std::string &fnLabel) {
  auto wsShape = countDimensions(self);
  const size_t ndims = wsShape.size();
  auto arrShape = signal.attr("shape");
  if (ndims != static_cast<size_t>(len(arrShape))) {
    std::ostringstream os;
    os << fnLabel
       << ": The number of  dimensions doe not match the current "
          "workspace size. Workspace="
       << ndims << " array=" << len(arrShape);
    throw std::invalid_argument(os.str());
  }
 
  for (size_t i = 0; i < ndims; ++i) {
    int arrDim = extract<int>(arrShape[i])();
    if (wsShape[i] != arrDim) {
      std::ostringstream os;
      os << fnLabel << ": The dimension size for the " << std::to_string(i) << "th dimension do not match. "
         << "Workspace dimension size=" << wsShape[i] << ", array size=" << arrDim;
      throw std::invalid_argument(os.str());
    }
  }
}
 
void copyToBuffer(const NDArray &values, Mantid::signal_t *dest) {
  int flags = NPY_ARRAY_F_CONTIGUOUS | NPY_ARRAY_ALIGNED;
  const int type = PyArray_TYPE(reinterpret_cast<PyArrayObject *>(values.ptr()));
  if (type == NPY_OBJECT || PyTypeNum_ISFLOAT(type))
    flags |= NPY_ARRAY_FORCECAST;
  PyObject *converted = PyArray_FROMANY(values.ptr(), NPY_DOUBLE, 0, 0, flags);
  // Throw directly rather than through throw_error_already_set(), so static analysis sees that execution stops here
  if (!converted)
    throw error_already_set();
  handle<> owner(converted); // releases the converted array, also if std::copy throws
  auto *array = reinterpret_cast<PyArrayObject *>(converted);
  // Raw pointer access below requires an aligned, native float64 array in Fortran order
  if (PyArray_TYPE(array) != NPY_DOUBLE || !PyArray_ISFARRAY_RO(array) || !PyArray_ISNOTSWAPPED(array))
    throw std::runtime_error("numpy did not convert the array to an aligned, native float64 array in Fortran order");
  const auto *source = static_cast<const double *>(PyArray_DATA(array));
  std::copy(source, source + PyArray_SIZE(array), dest);
}
 
void setSignalArray(IMDHistoWorkspace &self, const NDArray &signalValues) {
  throwIfSizeIncorrect(self, signalValues, "setSignalArray");
  copyToBuffer(signalValues, self.mutableSignalArray());
  self.clearOriginalWorkspaces();
}
 
void setErrorSquaredArray(IMDHistoWorkspace &self, const NDArray &errorSquared) {
  throwIfSizeIncorrect(self, errorSquared, "setErrorSquaredArray");
  copyToBuffer(errorSquared, self.mutableErrorSquaredArray());
  self.clearOriginalWorkspaces();
}
 
void setNumEventsArray(IMDHistoWorkspace &self, const NDArray &numEvents) {
  throwIfSizeIncorrect(self, numEvents, "setNumEventsArray");
  copyToBuffer(numEvents, self.mutableNumEventsArray());
  self.clearOriginalWorkspaces();
}
 
void setSignalAt(IMDHistoWorkspace &self, const size_t index, const double value) {
  if (index >= self.getNPoints())
    throw std::invalid_argument("setSignalAt: The index is greater than the "
                                "number of bins in the workspace");
 
  self.setSignalAt(index, value);
}
} // namespace
 
void export_IMDHistoWorkspace() {
  // IMDHistoWorkspace class
  class_<IMDHistoWorkspace, bases<IMDWorkspace, MultipleExperimentInfos>, boost::noncopyable>("IMDHistoWorkspace",
                                                                                              no_init)
      .def("getSignalArray", &getSignalArrayAsNumpyArray, arg("self"),
           "Returns a read-only numpy array containing the signal values")
 
      .def("getErrorSquaredArray", &getErrorSquaredArrayAsNumpyArray, arg("self"),
           "Returns a read-only numpy array containing the square of the error "
           "values")
 
      .def("getNumEventsArray", &getNumEventsArrayAsNumpyArray, arg("self"),
           "Returns a read-only numpy array containing the number of MD events "
           "in each bin")
 
      .def("signalAt", &IMDHistoWorkspace::signalAt, (arg("self"), arg("index")),
           return_value_policy<copy_non_const_reference>(), "Return a reference to the signal at the linear index")
 
      .def("errorSquaredAt", &IMDHistoWorkspace::errorSquaredAt, (arg("self"), arg("index")),
           return_value_policy<copy_non_const_reference>(), "Return the squared-errors at the linear index")
 
      .def("setSignalAt", &setSignalAt, (arg("self"), arg("index"), arg("value")),
           "Sets the signal at the specified index.")
 
      .def("setErrorSquaredAt", &IMDHistoWorkspace::setErrorSquaredAt, (arg("self"), arg("index"), arg("value")),
           "Sets the squared-error at the specified index.")
 
      .def("setSignalArray", &setSignalArray, (arg("self"), arg("signalValues")),
           "Sets the signal from a numpy array. The sizes must match the "
           "current workspace sizes. A ValueError is thrown if not. Any link "
           "to the original MDEventWorkspace(s) is cleared.")
 
      .def("setErrorSquaredArray", &setErrorSquaredArray, (arg("self"), arg("errorSquared")),
           "Sets the square of the errors from a numpy array. The sizes must "
           "match the current workspace sizes. A ValueError is thrown if not. Any "
           "link to the original MDEventWorkspace(s) is cleared.")
 
      .def("setNumEventsArray", &setNumEventsArray, (arg("self"), arg("numEvents")),
           "Sets the number of events from a numpy array. The sizes must match "
           "the current workspace sizes, otherwise a ValueError is thrown. Any link "
           "to the original MDEventWorkspace(s) is cleared.")
 
      .def("setTo", &IMDHistoWorkspace::setTo, (arg("self"), arg("signal"), arg("error_squared"), arg("num_events")),
           "Sets all signals/errors in the workspace to the given values")
 
      .def("getInverseVolume", &IMDHistoWorkspace::getInverseVolume, arg("self"),
           return_value_policy<return_by_value>(), "Return the inverse of volume of EACH cell in the workspace.")
 
      .def("getLinearIndex", (size_t (IMDHistoWorkspace::*)(size_t, size_t) const) & IMDHistoWorkspace::getLinearIndex,
           (arg("self"), arg("index1"), arg("index2")), return_value_policy<return_by_value>(),
           "Get the 1D linear index from the 2D array")
 
      .def("getLinearIndex",
           (size_t (IMDHistoWorkspace::*)(size_t, size_t, size_t) const) & IMDHistoWorkspace::getLinearIndex,
           (arg("self"), arg("index1"), arg("index2"), arg("index3")), return_value_policy<return_by_value>(),
           "Get the 1D linear index from the 3D array")
 
      .def("getLinearIndex",
           (size_t (IMDHistoWorkspace::*)(size_t, size_t, size_t, size_t) const) & IMDHistoWorkspace::getLinearIndex,
           (arg("self"), arg("index1"), arg("index2"), arg("index3"), arg("index4")),
           return_value_policy<return_by_value>(), "Get the 1D linear index from the 4D array")
 
      .def("getCenter", &IMDHistoWorkspace::getCenter, (arg("self"), arg("linear_index")),
           return_value_policy<return_by_value>(), "Return the position of the center of a bin at a given position")
 
      .def("setDisplayNormalization", &IMDHistoWorkspace::setDisplayNormalization, (arg("self"), arg("normalization")),
           "Sets the visual normalization of"
           " the workspace.");
 
  //-------------------------------------------------------------------------------------------------
 
  RegisterWorkspacePtrToPython<IMDHistoWorkspace>();
}