Alexandria  2.27.0
SDC-CH common library for the Euclid project
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
NdArray.icpp
Go to the documentation of this file.
1 
19 #ifdef NDARRAY_IMPL
20 
21 #include <ElementsKernel/Unused.h>
22 #include <algorithm>
23 #include <utility>
24 
25 namespace Euclid {
26 namespace NdArray {
27 
28 template <typename T>
29 template <bool Const>
30 NdArray<T>::Iterator<Const>::Iterator(ContainerInterface* container_ptr, size_t offset,
31  const std::vector<size_t>& shape, const std::vector<size_t>& strides,
32  size_t start)
33  : m_container_ptr(container_ptr)
34  , m_offset(offset)
35  , m_row_size(std::accumulate(shape.begin(), shape.end(), 1, std::multiplies<size_t>()))
36  , m_stride{strides.back()}
37  , m_i{start} {}
38 
39 template <typename T>
40 template <bool Const>
41 NdArray<T>::Iterator<Const>::Iterator(ContainerInterface* container_ptr, size_t offset, size_t row_size, size_t stride,
42  size_t start)
43  : m_container_ptr(container_ptr), m_offset(offset), m_row_size(row_size), m_stride(stride), m_i(start) {}
44 
45 template <typename T>
46 template <bool Const>
47 NdArray<T>::Iterator<Const>::Iterator(const Iterator<false>& other)
48  : m_container_ptr{other.m_container_ptr}
49  , m_offset{other.m_offset}
50  , m_row_size{other.m_row_size}
51  , m_stride{other.m_stride}
52  , m_i{other.m_i} {}
53 
54 template <typename T>
55 template <bool Const>
56 auto NdArray<T>::Iterator<Const>::operator++() -> Iterator& {
57  ++m_i;
58  return *this;
59 }
60 
61 template <typename T>
62 template <bool Const>
63 auto NdArray<T>::Iterator<Const>::operator++(int) -> const Iterator {
64  return Iterator{m_container_ptr, m_offset, m_row_size, m_stride, m_i + 1};
65 }
66 
67 template <typename T>
68 template <bool Const>
69 bool NdArray<T>::Iterator<Const>::operator==(const Iterator& other) const {
70  return m_container_ptr == other.m_container_ptr && m_offset == other.m_offset && m_i == other.m_i;
71 }
72 
73 template <typename T>
74 template <bool Const>
75 bool NdArray<T>::Iterator<Const>::operator!=(const Iterator& other) const {
76  return !(*this == other);
77 }
78 
79 template <typename T>
80 template <bool Const>
81 auto NdArray<T>::Iterator<Const>::operator*() -> value_t& {
82  return operator[](0);
83 }
84 
85 template <typename T>
86 template <bool Const>
87 auto NdArray<T>::Iterator<Const>::operator*() const -> value_t {
88  return operator[](0);
89 }
90 
91 template <typename T>
92 template <bool Const>
93 auto NdArray<T>::Iterator<Const>::operator+=(size_t n) -> Iterator& {
94  m_i += n;
95  return *this;
96 }
97 
98 template <typename T>
99 template <bool Const>
100 auto NdArray<T>::Iterator<Const>::operator+(size_t n) const -> Iterator {
101  return Iterator{m_container_ptr, m_offset, m_row_size, m_stride, m_i + n};
102 }
103 
104 template <typename T>
105 template <bool Const>
106 auto NdArray<T>::Iterator<Const>::operator-=(size_t n) -> Iterator& {
107  assert(n <= m_i);
108  m_i -= n;
109  return *this;
110 }
111 
112 template <typename T>
113 template <bool Const>
114 auto NdArray<T>::Iterator<Const>::operator-(size_t n) const -> Iterator {
115  assert(n <= m_i);
116  return Iterator{m_container_ptr, m_offset, m_row_size, m_stride, m_i - n};
117 }
118 
119 template <typename T>
120 template <bool Const>
121 auto NdArray<T>::Iterator<Const>::operator-(const Iterator& other) -> difference_type {
122  assert(m_container_ptr == other.m_container_ptr);
123  return m_i - other.m_i;
124 }
125 
126 template <typename T>
127 template <bool Const>
128 auto NdArray<T>::Iterator<Const>::operator[](size_t i) -> value_t& {
129  return m_container_ptr->get(m_offset + (m_i + i) * m_stride);
130 }
131 
132 template <typename T>
133 template <bool Const>
134 auto NdArray<T>::Iterator<Const>::operator[](size_t i) const -> value_t {
135  return m_container_ptr->get(m_offset + (m_i + i) * m_stride);
136 }
137 
138 template <typename T>
139 template <bool Const>
140 bool NdArray<T>::Iterator<Const>::operator<(const Iterator& other) {
141  assert(m_container_ptr == other.m_container_ptr);
142  return m_i < other.m_i;
143 }
144 
145 template <typename T>
146 template <bool Const>
147 bool NdArray<T>::Iterator<Const>::operator>(const Iterator& other) {
148  assert(m_container_ptr == other.m_container_ptr);
149  return m_i > other.m_i;
150 }
151 
152 template <typename T>
153 NdArray<T>::NdArray(std::vector<size_t> shape_)
154  : m_details_ptr(new Details{0,
155  std::accumulate(shape_.begin(), shape_.end(), 1u, std::multiplies<size_t>()),
156  0,
157  std::move(shape_),
158  {},
159  {},
160  nullptr}) {
161  m_details_ptr->m_container = std::make_shared<ContainerWrapper<std::vector>>(m_details_ptr->m_size);
162  update_strides();
163 }
164 
165 template <typename T>
166 template <template <class...> class Container>
167 NdArray<T>::NdArray(std::vector<size_t> shape_, const Container<T>& data)
168  : m_details_ptr(new Details{0,
169  std::accumulate(shape_.begin(), shape_.end(), 1u, std::multiplies<size_t>()),
170  0,
171  std::move(shape_),
172  {},
173  {},
174  std::make_shared<ContainerWrapper<Container>>(data)}) {
175  if (m_details_ptr->m_size != m_details_ptr->m_container->size()) {
176  throw std::invalid_argument("Data size does not match the shape");
177  }
178  update_strides();
179 }
180 
181 template <typename T>
182 template <template <class...> class Container>
183 NdArray<T>::NdArray(std::vector<size_t> shape_, Container<T>&& data)
184  : m_details_ptr(new Details{0,
185  std::accumulate(shape_.begin(), shape_.end(), 1u, std::multiplies<size_t>()),
186  0,
187  std::move(shape_),
188  {},
189  {},
190  std::make_shared<ContainerWrapper<Container>>(std::move(data))}) {
191  if (m_details_ptr->m_size != m_details_ptr->m_container->size()) {
192  throw std::invalid_argument("Data size does not match the shape");
193  }
194  update_strides();
195 }
196 
197 template <typename T>
198 template <template <class...> class Container>
199 NdArray<T>::NdArray(std::vector<size_t> shape_, std::vector<size_t> strides_, Container<T>&& data)
200  : m_details_ptr(new Details{0,
201  std::accumulate(shape_.begin(), shape_.end(), 1u, std::multiplies<size_t>()),
202  0,
203  std::move(shape_),
204  std::move(strides_),
205  {},
206  std::make_shared<ContainerWrapper<Container>>(std::move(data))}) {
207  if (m_details_ptr->m_shape.size() != m_details_ptr->m_stride_size.size()) {
208  throw std::invalid_argument("The size of the shape and strides parameters do not match");
209  }
210  if (!std::is_sorted(m_details_ptr->m_stride_size.rbegin(), m_details_ptr->m_stride_size.rend())) {
211  throw std::runtime_error("Only C style arrays are supported");
212  }
213 }
214 
215 template <typename T>
216 template <typename II>
217 NdArray<T>::NdArray(std::vector<size_t> shape_, II ibegin, II iend)
218  : m_details_ptr(new Details{0,
219  std::accumulate(shape_.begin(), shape_.end(), 1u, std::multiplies<size_t>()),
220  0,
221  std::move(shape_),
222  {},
223  {},
224  std::make_shared<ContainerWrapper<std::vector>>(ibegin, iend)}) {
225  if (m_details_ptr->m_size != m_details_ptr->m_container->size()) {
226  throw std::invalid_argument("Data size does not match the shape");
227  }
228  update_strides();
229 }
230 
231 template <typename T>
232 NdArray<T>::NdArray(const self_type* other) : m_details_ptr(make_unique<Details>(*other->m_details_ptr)) {
233  m_details_ptr->m_container = other->m_details_ptr->m_container->copy();
234 }
235 
236 inline std::vector<size_t> appendAttrShape(std::vector<size_t> shape, size_t append) {
237  if (append)
238  shape.push_back(append);
239  return shape;
240 }
241 
242 template <typename T>
243 template <typename... Args>
244 NdArray<T>::NdArray(const std::vector<size_t>& shape_, const std::vector<std::string>& attr_names, Args&&... args)
245  : NdArray(appendAttrShape(shape_, attr_names.size()), std::forward<Args>(args)...) {
246  m_details_ptr->m_attr_names = attr_names;
247 }
248 
249 template <typename T>
250 auto NdArray<T>::reshape(const std::vector<size_t>& new_shape) -> self_type& {
251  if (!m_details_ptr->m_attr_names.empty())
252  throw std::invalid_argument("Can not reshape arrays with attribute names");
253 
254  size_t new_size = std::accumulate(new_shape.begin(), new_shape.end(), 1, std::multiplies<size_t>());
255  if (new_size != m_details_ptr->m_size) {
256  throw std::range_error("New shape does not match the number of contained elements");
257  }
258  m_details_ptr->m_shape = new_shape;
259  update_strides();
260  return *this;
261 }
262 
263 template <typename T>
264 template <typename... D>
265 auto NdArray<T>::reshape(size_t i, D... rest) -> self_type& {
266  std::vector<size_t> acc{i};
267  return reshape_helper(acc, rest...);
268 }
269 
270 template <typename T>
271 const T& NdArray<T>::front() const {
272  return m_details_ptr->m_container->get(m_details_ptr->m_offset);
273 }
274 
275 template <typename T>
276 T& NdArray<T>::at(const std::vector<size_t>& coords) {
277  auto offset = get_offset(coords);
278  // cppcheck-suppress "returnTempReference"
279  return m_details_ptr->m_container->get(offset);
280 }
281 
282 template <typename T>
283 const T& NdArray<T>::at(const std::vector<size_t>& coords) const {
284  auto offset = get_offset(coords);
285  // cppcheck-suppress returnTempReference
286  return m_details_ptr->m_container->get(offset);
287 }
288 
289 template <typename T>
290 T& NdArray<T>::at(const std::vector<size_t>& coords, const std::string& attr) {
291  auto offset = get_offset(coords, attr);
292  // cppcheck-suppress returnTempReference
293  return m_details_ptr->m_container->get(offset);
294 }
295 
296 template <typename T>
297 const T& NdArray<T>::at(const std::vector<size_t>& coords, const std::string& attr) const {
298  auto offset = get_offset(coords, attr);
299  // cppcheck-suppress returnTempReference
300  return m_details_ptr->m_container->get(offset);
301 }
302 
303 template <typename T>
304 template <typename... D>
305 T& NdArray<T>::at(size_t i, D... rest) {
306  return at_helper(m_details_ptr->m_offset, 0, i, rest...);
307 }
308 
309 template <typename T>
310 template <typename... D>
311 const T& NdArray<T>::at(size_t i, D... rest) const {
312  return at_helper(m_details_ptr->m_offset, 0, i, rest...);
313 }
314 
315 template <typename T>
316 auto NdArray<T>::begin() -> iterator {
317  return iterator{m_details_ptr->m_container.get(), m_details_ptr->m_offset, m_details_ptr->m_shape,
318  m_details_ptr->m_stride_size, 0};
319 }
320 
321 template <typename T>
322 auto NdArray<T>::end() -> iterator {
323  return iterator{m_details_ptr->m_container.get(), m_details_ptr->m_offset, m_details_ptr->m_shape,
324  m_details_ptr->m_stride_size, size()};
325 }
326 
327 template <typename T>
328 auto NdArray<T>::begin() const -> const_iterator {
329  return const_iterator{m_details_ptr->m_container.get(), m_details_ptr->m_offset, m_details_ptr->m_shape,
330  m_details_ptr->m_stride_size, 0};
331 }
332 
333 template <typename T>
334 auto NdArray<T>::end() const -> const_iterator {
335  return const_iterator{m_details_ptr->m_container.get(), m_details_ptr->m_offset, m_details_ptr->m_shape,
336  m_details_ptr->m_stride_size, size()};
337 }
338 
339 template <typename T>
340 size_t NdArray<T>::size() const {
341  return m_details_ptr->m_size;
342 }
343 
344 template <typename T>
345 bool NdArray<T>::operator==(const self_type& b) const {
346  if (shape() != b.shape())
347  return false;
348  for (auto ai = begin(), bi = b.begin(); ai != end() && bi != b.end(); ++ai, ++bi) {
349  if (*ai != *bi)
350  return false;
351  }
352  return true;
353 }
354 
355 template <typename T>
356 bool NdArray<T>::operator!=(const self_type& b) const {
357  return !(*this == b);
358 }
359 
360 template <typename T>
361 const std::vector<std::string>& NdArray<T>::attributes() const {
362  return m_details_ptr->m_attr_names;
363 }
364 
365 template <typename T>
366 auto NdArray<T>::concatenate(const self_type& other) -> self_type& {
367  // Verify dimensionality
368  if (m_details_ptr->m_shape.size() != other.m_details_ptr->m_shape.size()) {
369  throw std::length_error("Can not concatenate arrays with different dimensionality");
370  }
371  for (size_t i = 1; i < m_details_ptr->m_shape.size(); ++i) {
372  if (m_details_ptr->m_shape[i] != other.m_details_ptr->m_shape[i])
373  throw std::length_error("The size of all axis except for the first one must match");
374  }
375 
376  // New shape
377  auto old_size = size();
378  auto new_shape = m_details_ptr->m_shape;
379  new_shape[0] += other.m_details_ptr->m_shape[0];
380 
381  // Resize container
382  m_details_ptr->m_container->resize(new_shape);
383 
384  // Copy to the end
385  std::copy(std::begin(other), std::end(other), &m_details_ptr->m_container->get(0) + old_size);
386  // Done!
387  m_details_ptr->m_shape = new_shape;
388  m_details_ptr->m_size += other.m_details_ptr->m_size;
389  return *this;
390 }
391 
392 template <typename T>
393 NdArray<T>::NdArray(std::shared_ptr<ContainerInterface> container, size_t offset, std::vector<size_t> shape_,
395  : m_details_ptr(new Details{offset, std::accumulate(shape_.begin(), shape_.end(), 1ul, std::multiplies<size_t>()),
396  stride.front() * shape_.front(), std::move(shape_), std::move(stride),
397  std::move(attr_names), std::move(container)}) {}
398 
399 template <typename T>
400 NdArray<T>::NdArray(const self_type& other) : m_details_ptr(make_unique<Details>(*other.m_details_ptr)) {}
401 
402 template <typename T>
403 NdArray<T>& NdArray<T>::operator=(const self_type& other) {
404  m_details_ptr = make_unique<Details>(*other.m_details_ptr);
405  return *this;
406 }
407 
408 template <typename T>
409 auto NdArray<T>::slice(size_t i) -> self_type {
410  if (m_details_ptr->m_shape.size() <= 1) {
411  throw std::out_of_range("Can not slice a one dimensional array");
412  }
414  if (!m_details_ptr->m_attr_names.empty()) {
415  attrs.resize(m_details_ptr->m_attr_names.size());
416  std::copy(m_details_ptr->m_attr_names.begin(), m_details_ptr->m_attr_names.end(), attrs.begin());
417  }
418  if (i >= m_details_ptr->m_shape[0]) {
419  throw std::out_of_range("Axis 0 out of range");
420  }
421  auto offset = m_details_ptr->m_offset + i * m_details_ptr->m_stride_size[0];
422  std::vector<size_t> stride_(m_details_ptr->m_stride_size.begin() + 1, m_details_ptr->m_stride_size.end());
423  std::vector<size_t> shape_(m_details_ptr->m_shape.begin() + 1, m_details_ptr->m_shape.end());
424  return NdArray(m_details_ptr->m_container, offset, std::move(shape_), std::move(stride_), std::move(attrs));
425 }
426 
427 template <typename T>
428 auto NdArray<T>::slice(size_t i) const -> const self_type {
429  return const_cast<NdArray<T>*>(this)->slice(i);
430 }
431 
432 template <typename T>
433 auto NdArray<T>::rslice(size_t i) -> self_type {
434  if (m_details_ptr->m_shape.size() <= 1) {
435  throw std::out_of_range("Can not slice a one dimensional array");
436  }
437  if (!m_details_ptr->m_attr_names.empty()) {
438  throw std::invalid_argument("Can not slice on the last axis for arrays with attribute names");
439  }
440  if (i >= m_details_ptr->m_shape.back()) {
441  throw std::out_of_range("Axis -1 out of range");
442  }
443  auto offset = m_details_ptr->m_offset + i * m_details_ptr->m_stride_size.back();
444  std::vector<size_t> strides_(m_details_ptr->m_stride_size.begin(), m_details_ptr->m_stride_size.end() - 1);
445  std::vector<size_t> shape_(m_details_ptr->m_shape.begin(), m_details_ptr->m_shape.end() - 1);
446  return NdArray(m_details_ptr->m_container, offset, std::move(shape_), std::move(strides_),
447  m_details_ptr->m_attr_names);
448 }
449 
450 template <typename T>
451 auto NdArray<T>::rslice(size_t i) const -> const self_type {
452  return const_cast<NdArray<T>*>(this)->rslice(i);
453 }
454 
455 template <typename T>
456 void NdArray<T>::next_slice() {
457  m_details_ptr->m_offset += m_details_ptr->m_total_stride;
458 }
459 
460 template <typename T>
461 size_t NdArray<T>::get_offset(const std::vector<size_t>& coords) const {
462  if (coords.size() != m_details_ptr->m_shape.size()) {
463  throw std::out_of_range("Invalid number of coordinates, got " + std::to_string(coords.size()) + ", expected " +
464  std::to_string(m_details_ptr->m_shape.size()));
465  }
466 
467  size_t offset = m_details_ptr->m_offset;
468  for (size_t i = 0; i < coords.size(); ++i) {
469  if (coords[i] >= m_details_ptr->m_shape[i]) {
470  throw std::out_of_range(std::to_string(coords[i]) + " >= " + std::to_string(m_details_ptr->m_shape[i]) +
471  " for axis " + std::to_string(i));
472  }
473  offset += coords[i] * m_details_ptr->m_stride_size[i];
474  }
475 
476  assert(offset < m_details_ptr->m_container->nbytes());
477  return offset;
478 }
479 
480 template <typename T>
481 size_t NdArray<T>::get_attr_offset(const std::string& attr) const {
482  auto i = std::find(m_details_ptr->m_attr_names.begin(), m_details_ptr->m_attr_names.end(), attr);
483  if (i == m_details_ptr->m_attr_names.end())
484  throw std::out_of_range(attr);
485  return (i - m_details_ptr->m_attr_names.begin()) * sizeof(T);
486 }
487 
488 template <typename T>
489 void NdArray<T>::update_strides() {
490  m_details_ptr->m_stride_size.resize(m_details_ptr->m_shape.size());
491 
492  size_t acc = sizeof(T);
493  for (size_t i = m_details_ptr->m_stride_size.size(); i > 0; --i) {
494  m_details_ptr->m_stride_size[i - 1] = acc;
495  acc *= m_details_ptr->m_shape[i - 1];
496  }
497 
498  m_details_ptr->m_total_stride = m_details_ptr->m_shape.front() * m_details_ptr->m_stride_size.front();
499 }
500 
504 template <typename T>
505 template <typename... D>
506 T& NdArray<T>::at_helper(size_t offset_acc, size_t axis, size_t i, D... rest) {
507  assert(axis <= m_details_ptr->m_shape.size() && i < m_details_ptr->m_shape[axis]);
508  offset_acc += i * m_details_ptr->m_stride_size[axis];
509  return at_helper(offset_acc, ++axis, rest...);
510 }
511 
512 template <typename T>
513 T& NdArray<T>::at_helper(size_t offset_acc, ELEMENTS_UNUSED size_t axis) {
514  assert(axis == m_details_ptr->m_shape.size());
515  assert(offset_acc < m_details_ptr->m_container->nbytes());
516  return m_details_ptr->m_container->get(offset_acc);
517 }
518 
519 template <typename T>
520 T& NdArray<T>::at_helper(size_t offset_acc, ELEMENTS_UNUSED size_t axis, const std::string& attr) {
521  offset_acc += get_attr_offset(attr);
522  assert(axis == m_details_ptr->m_shape.size() - 1);
523  assert(offset_acc < m_details_ptr->m_container->nbytes());
524  return m_details_ptr->m_container->get(offset_acc);
525 }
526 
527 template <typename T>
528 template <typename... D>
529 const T& NdArray<T>::at_helper(size_t offset_acc, size_t axis, size_t i, D... rest) const {
530  assert(axis <= m_details_ptr->m_shape.size() && i < m_details_ptr->m_shape[axis]);
531  offset_acc += i * m_details_ptr->m_stride_size[axis];
532  return at_helper(offset_acc, ++axis, rest...);
533 }
534 
535 template <typename T>
536 const T& NdArray<T>::at_helper(size_t offset_acc, ELEMENTS_UNUSED size_t axis) const {
537  assert(axis == m_details_ptr->m_shape.size());
538  assert(offset_acc < m_details_ptr->m_container->nbytes());
539  // cppcheck-suppress returnTempReference
540  return m_details_ptr->m_container->get(offset_acc);
541 }
542 
543 template <typename T>
544 const T& NdArray<T>::at_helper(size_t offset_acc, ELEMENTS_UNUSED size_t axis, const std::string& attr) const {
545  offset_acc += get_attr_offset(attr);
546  assert(axis == m_details_ptr->m_shape.size() - 1);
547  assert(offset_acc < m_details_ptr->m_container->nbytes());
548  // cppcheck-suppress returnTempReference
549  return m_details_ptr->m_container->get(offset_acc);
550 }
551 
552 template <typename T>
553 template <typename... D>
554 auto NdArray<T>::reshape_helper(std::vector<size_t>& acc, size_t i, D... rest) -> self_type& {
555  acc.push_back(i);
556  return reshape_helper(acc, rest...);
557 }
558 
559 template <typename T>
560 auto NdArray<T>::reshape_helper(std::vector<size_t>& acc) -> self_type& {
561  return reshape(acc);
562 }
563 
564 template <typename T>
565 std::ostream& operator<<(std::ostream& out, const NdArray<T>& ndarray) {
566  auto shape = ndarray.shape();
567 
568  if (ndarray.size()) {
569  out << "<";
570  size_t i;
571  for (i = 0; i < shape.size() - 1; ++i) {
572  out << shape[i] << ",";
573  }
574  out << shape[i] << ">";
575 
576  auto iter = ndarray.begin(), end_iter = ndarray.end() - 1;
577  for (; iter != end_iter; ++iter) {
578  out << *iter << ",";
579  }
580  out << *iter;
581  }
582  return out;
583 }
584 
585 } // end of namespace NdArray
586 } // end of namespace Euclid
587 
588 #endif // NDARRAY_IMPL
T copy(T...args)
T front(T...args)
T to_string(T...args)
bool operator!=(const Euclid::SourceCatalog::Source::id_type &a, const Euclid::SourceCatalog::Source::id_type &b)
boost::variant specifies an equality operator (==), but, in older boost versions, not an inequality o...
Definition: Source.h:145
T end(T...args)
T resize(T...args)
STL class.
T push_back(T...args)
T move(T...args)
T find(T...args)
T size(T...args)
T is_sorted(T...args)
T begin(T...args)
std::unique_ptr< T > make_unique(Args &&...args)
Constructs an object of type T and wraps it in a std::unique_ptr using args as the parameter list for...
Definition: memory_tools.h:42
T back(T...args)
T accumulate(T...args)
STL class.
STL class.
T forward(T...args)