Mantid
Loading...
Searching...
No Matches
IMDHistoWorkspace.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
13
14#include <boost/python/class.hpp>
15#include <boost/python/copy_non_const_reference.hpp>
16#define PY_ARRAY_UNIQUE_SYMBOL API_ARRAY_API
17#define NO_IMPORT_ARRAY
18#include <numpy/arrayobject.h>
19
20using namespace Mantid::API;
24using namespace boost::python;
25
27
29extern template int NDArrayTypeIndex<float>::typenum;
30extern template int NDArrayTypeIndex<double>::typenum;
31} // namespace Mantid::PythonInterface::Converters
32
33namespace {
41PyObject *WrapReadOnlyNumpyFArray(const Mantid::signal_t *arr, std::vector<Py_intptr_t> dims) {
43#if NPY_API_VERSION >= 0x00000007 //(1.7)
44 auto *nparray = reinterpret_cast<PyArrayObject *>(
45 PyArray_New(&PyArray_Type, static_cast<int>(dims.size()), &dims[0], datatype, nullptr,
46 static_cast<void *>(const_cast<double *>(arr)), 0, NPY_ARRAY_FARRAY, nullptr));
47 PyArray_CLEARFLAGS(nparray, NPY_ARRAY_WRITEABLE);
48#else
49 PyArrayObject *nparray =
50 (PyArrayObject *)PyArray_New(&PyArray_Type, static_cast<int>(dims.size()), &dims[0], datatype, nullptr,
51 static_cast<void *>(const_cast<double *>(arr)), 0, NPY_FARRAY, nullptr);
52 nparray->flags &= ~NPY_WRITEABLE;
53#endif
54 return reinterpret_cast<PyObject *>(nparray);
55}
56
62std::vector<Py_intptr_t> countDimensions(const IMDHistoWorkspace &self) {
63 size_t ndims = self.getNumDims();
64 std::vector<size_t> nd;
65 nd.reserve(ndims);
66
67 // invert dimensions in C way, e.g. slowest changing ndim goes first
68 for (size_t i = 0; i < ndims; ++i) {
69 nd.emplace_back(self.getDimension(i)->getNBins());
70 }
71
72 ndims = nd.size();
73 std::vector<Py_intptr_t> dims(ndims);
74 for (size_t i = 0; i < ndims; ++i)
75 dims[i] = static_cast<Py_intptr_t>(nd[i]);
76
77 if (dims.empty()) {
78 throw std::runtime_error("Workspace has zero dimensions!");
79 } else {
80 return dims;
81 }
82}
83
88PyObject *getSignalArrayAsNumpyArray(const IMDHistoWorkspace &self) {
89 auto dims = countDimensions(self);
90 return WrapReadOnlyNumpyFArray(self.getSignalArray(), dims);
91}
92
97PyObject *getErrorSquaredArrayAsNumpyArray(const IMDHistoWorkspace &self) {
98 auto dims = countDimensions(self);
99 return WrapReadOnlyNumpyFArray(self.getErrorSquaredArray(), dims);
100}
101
106PyObject *getNumEventsArrayAsNumpyArray(const IMDHistoWorkspace &self) {
107 auto dims = countDimensions(self);
108 return WrapReadOnlyNumpyFArray(self.getNumEventsArray(), dims);
109}
110
118void throwIfSizeIncorrect(const IMDHistoWorkspace &self, const NDArray &signal, const std::string &fnLabel) {
119 auto wsShape = countDimensions(self);
120 const size_t ndims = wsShape.size();
121 auto arrShape = signal.attr("shape");
122 if (ndims != static_cast<size_t>(len(arrShape))) {
123 std::ostringstream os;
124 os << fnLabel
125 << ": The number of dimensions doe not match the current "
126 "workspace size. Workspace="
127 << ndims << " array=" << len(arrShape);
128 throw std::invalid_argument(os.str());
129 }
130
131 for (size_t i = 0; i < ndims; ++i) {
132 int arrDim = extract<int>(arrShape[i])();
133 if (wsShape[i] != arrDim) {
134 std::ostringstream os;
135 os << fnLabel << ": The dimension size for the " << std::to_string(i) << "th dimension do not match. "
136 << "Workspace dimension size=" << wsShape[i] << ", array size=" << arrDim;
137 throw std::invalid_argument(os.str());
138 }
139 }
140}
141
149void setSignalArray(IMDHistoWorkspace &self, const NDArray &signalValues) {
150 throwIfSizeIncorrect(self, signalValues, "setSignalArray");
151 object rav = signalValues.attr("ravel")("F");
152 object flattened = rav.attr("flat");
153 auto length = len(flattened);
154 for (auto i = 0; i < length; ++i) {
155 self.setSignalAt(i, extract<double>(flattened[i])());
156 }
157}
158
166void setErrorSquaredArray(IMDHistoWorkspace &self, const NDArray &errorSquared) {
167 throwIfSizeIncorrect(self, errorSquared, "setErrorSquaredArray");
168 object rav = errorSquared.attr("ravel")("F");
169 object flattened = rav.attr("flat");
170 auto length = len(flattened);
171 for (auto i = 0; i < length; ++i) {
172 self.setErrorSquaredAt(i, extract<double>(flattened[i])());
173 }
174}
175
179void setSignalAt(IMDHistoWorkspace &self, const size_t index, const double value) {
180 if (index >= self.getNPoints())
181 throw std::invalid_argument("setSignalAt: The index is greater than the "
182 "number of bins in the workspace");
183
184 self.setSignalAt(index, value);
185}
186} // namespace
187
189 // IMDHistoWorkspace class
190 class_<IMDHistoWorkspace, bases<IMDWorkspace, MultipleExperimentInfos>, boost::noncopyable>("IMDHistoWorkspace",
191 no_init)
192 .def("getSignalArray", &getSignalArrayAsNumpyArray, arg("self"),
193 "Returns a read-only numpy array containing the signal values")
194
195 .def("getErrorSquaredArray", &getErrorSquaredArrayAsNumpyArray, arg("self"),
196 "Returns a read-only numpy array containing the square of the error "
197 "values")
198
199 .def("getNumEventsArray", &getNumEventsArrayAsNumpyArray, arg("self"),
200 "Returns a read-only numpy array containing the number of MD events "
201 "in each bin")
202
203 .def("signalAt", &IMDHistoWorkspace::signalAt, (arg("self"), arg("index")),
204 return_value_policy<copy_non_const_reference>(), "Return a reference to the signal at the linear index")
205
206 .def("errorSquaredAt", &IMDHistoWorkspace::errorSquaredAt, (arg("self"), arg("index")),
207 return_value_policy<copy_non_const_reference>(), "Return the squared-errors at the linear index")
208
209 .def("setSignalAt", &setSignalAt, (arg("self"), arg("index"), arg("value")),
210 "Sets the signal at the specified index.")
211
212 .def("setErrorSquaredAt", &IMDHistoWorkspace::setErrorSquaredAt, (arg("self"), arg("index"), arg("value")),
213 "Sets the squared-error at the specified index.")
214
215 .def("setSignalArray", &setSignalArray, (arg("self"), arg("signalValues")),
216 "Sets the signal from a numpy array. The sizes must match the "
217 "current workspace sizes. A ValueError is thrown if not")
218
219 .def("setErrorSquaredArray", &setErrorSquaredArray, (arg("self"), arg("errorSquared")),
220 "Sets the square of the errors from a numpy array. The sizes must "
221 "match the current workspace sizes. A ValueError is thrown if not")
222
223 .def("setTo", &IMDHistoWorkspace::setTo, (arg("self"), arg("signal"), arg("error_squared"), arg("num_events")),
224 "Sets all signals/errors in the workspace to the given values")
225
226 .def("getInverseVolume", &IMDHistoWorkspace::getInverseVolume, arg("self"),
227 return_value_policy<return_by_value>(), "Return the inverse of volume of EACH cell in the workspace.")
228
229 .def("getLinearIndex", (size_t(IMDHistoWorkspace::*)(size_t, size_t) const) & IMDHistoWorkspace::getLinearIndex,
230 (arg("self"), arg("index1"), arg("index2")), return_value_policy<return_by_value>(),
231 "Get the 1D linear index from the 2D array")
232
233 .def("getLinearIndex",
234 (size_t(IMDHistoWorkspace::*)(size_t, size_t, size_t) const) & IMDHistoWorkspace::getLinearIndex,
235 (arg("self"), arg("index1"), arg("index2"), arg("index3")), return_value_policy<return_by_value>(),
236 "Get the 1D linear index from the 3D array")
237
238 .def("getLinearIndex",
239 (size_t(IMDHistoWorkspace::*)(size_t, size_t, size_t, size_t) const) & IMDHistoWorkspace::getLinearIndex,
240 (arg("self"), arg("index1"), arg("index2"), arg("index3"), arg("index4")),
241 return_value_policy<return_by_value>(), "Get the 1D linear index from the 4D array")
242
243 .def("getCenter", &IMDHistoWorkspace::getCenter, (arg("self"), arg("linear_index")),
244 return_value_policy<return_by_value>(), "Return the position of the center of a bin at a given position")
245
246 .def("setDisplayNormalization", &IMDHistoWorkspace::setDisplayNormalization, (arg("self"), arg("normalization")),
247 "Sets the visual normalization of"
248 " the workspace.");
249
250 //-------------------------------------------------------------------------------------------------
251
253}
double value
The value of the point.
Definition: FitMW.cpp:51
#define GET_POINTER_SPECIALIZATION(TYPE)
Definition: GetPointer.h:17
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
tagPyArrayObject PyArrayObject
void export_IMDHistoWorkspace()
Abstract interface to MDHistoWorkspace, for use in exposing to Python.
virtual signal_t & errorSquaredAt(size_t index)=0
virtual Mantid::Kernel::VMD getCenter(size_t linearIndex) const =0
virtual const signal_t * getErrorSquaredArray() const =0
virtual const signal_t * getNumEventsArray() const =0
virtual signal_t & signalAt(size_t index)=0
virtual void setErrorSquaredAt(size_t index, signal_t value)=0
virtual const signal_t * getSignalArray() const =0
virtual void setDisplayNormalization(const Mantid::API::MDNormalization &preferredNormalization)=0
virtual void setSignalAt(size_t index, signal_t value)=0
virtual coord_t getInverseVolume() const =0
See the MDHistoWorkspace definition for descriptions of these.
virtual size_t getLinearIndex(size_t index1, size_t index2) const =0
virtual void setTo(signal_t signal, signal_t errorSquared, signal_t numEvents)=0
virtual uint64_t getNPoints() const =0
Get the number of points associated with the workspace.
virtual std::shared_ptr< const Mantid::Geometry::IMDDimension > getDimension(size_t index) const
Get a dimension.
Definition: MDGeometry.cpp:161
virtual size_t getNumDims() const
Definition: MDGeometry.cpp:148
Thin object wrapper around a numpy array.
Definition: NDArray.h:31
double signal_t
Typedef for the signal recorded in a MDBox, etc.
Definition: MDTypes.h:36
std::string to_string(const wide_integer< Bits, Signed > &n)
Defines a mapping between C++ type given by the template parameter and numpy type enum NPY_TYPES.
Encapsulates the registration required for an interface type T that sits on top of a Kernel::DataItem...