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