Alexandria  2.22.0
Please provide a description of the project.
NdArray.icpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012-2021 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 
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) -> 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 std::memcmp(this, &other, sizeof(*this)) == 0;
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) -> 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) -> 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->at(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->at(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(const std::vector<size_t>& shape_)
154  : m_offset(0)
155  , m_shape(shape_)
156  , m_size(std::accumulate(m_shape.begin(), m_shape.end(), 1u, std::multiplies<size_t>()))
157  , m_container(new ContainerWrapper<std::vector>(m_size)) {
158  update_strides();
159 }
160 
161 template <typename T>
162 template <template <class...> class Container>
163 NdArray<T>::NdArray(const std::vector<size_t>& shape_, const Container<T>& data)
164  : m_offset(0)
165  , m_shape(shape_)
166  , m_size(std::accumulate(m_shape.begin(), m_shape.end(), 1u, std::multiplies<size_t>()))
167  , m_container(new ContainerWrapper<Container>(data)) {
168  if (m_size != m_container->size()) {
169  throw std::invalid_argument("Data size does not match the shape");
170  }
171  update_strides();
172 }
173 
174 template <typename T>
175 template <template <class...> class Container>
176 NdArray<T>::NdArray(const std::vector<size_t>& shape_, Container<T>&& data)
177  : m_offset(0)
178  , m_shape(shape_)
179  , m_size(std::accumulate(m_shape.begin(), m_shape.end(), 1u, std::multiplies<size_t>()))
180  , m_container(new ContainerWrapper<Container>(std::move(data))) {
181  if (m_size != m_container->size()) {
182  throw std::invalid_argument("Data size does not match the shape");
183  }
184  update_strides();
185 }
186 
187 template <typename T>
188 template <typename II>
189 NdArray<T>::NdArray(const std::vector<size_t>& shape_, II ibegin, II iend)
190  : m_offset(0)
191  , m_shape(shape_)
192  , m_size(std::accumulate(m_shape.begin(), m_shape.end(), 1u, std::multiplies<size_t>()))
193  , m_container(new ContainerWrapper<std::vector>(ibegin, iend)) {
194  if (m_size != m_container->size()) {
195  throw std::invalid_argument("Data size does not match the shape");
196  }
197  update_strides();
198 }
199 
200 template <typename T>
201 NdArray<T>::NdArray(const self_type* other)
202  : m_offset(0)
203  , m_shape(other->m_shape)
204  , m_attr_names(other->m_attr_names)
205  , m_size(std::accumulate(m_shape.begin(), m_shape.end(), 1u, std::multiplies<size_t>()))
206  , m_container(other->m_container->copy()) {
207  update_strides();
208 }
209 
210 inline std::vector<size_t> appendAttrShape(std::vector<size_t> shape, size_t append) {
211  if (append)
212  shape.push_back(append);
213  return shape;
214 }
215 
216 template <typename T>
217 template <typename... Args>
218 NdArray<T>::NdArray(const std::vector<size_t>& shape_, const std::vector<std::string>& attr_names, Args&&... args)
219  : NdArray(appendAttrShape(shape_, attr_names.size()), std::forward<Args>(args)...) {
220  m_attr_names = attr_names;
221 }
222 
223 template <typename T>
224 auto NdArray<T>::reshape(const std::vector<size_t>& new_shape) -> self_type& {
225  if (!m_attr_names.empty())
226  throw std::invalid_argument("Can not reshape arrays with attribute names");
227 
228  size_t new_size = std::accumulate(new_shape.begin(), new_shape.end(), 1, std::multiplies<size_t>());
229  if (new_size != m_size) {
230  throw std::range_error("New shape does not match the number of contained elements");
231  }
232  m_shape = new_shape;
233  update_strides();
234  return *this;
235 }
236 
237 template <typename T>
238 template <typename... D>
239 auto NdArray<T>::reshape(size_t i, D... rest) -> self_type& {
240  std::vector<size_t> acc{i};
241  return reshape_helper(acc, rest...);
242 }
243 
244 template <typename T>
245 T& NdArray<T>::at(const std::vector<size_t>& coords) {
246  auto offset = get_offset(coords);
247  // cppcheck-suppress "returnTempReference"
248  return m_container->at(offset);
249 }
250 
251 template <typename T>
252 const T& NdArray<T>::at(const std::vector<size_t>& coords) const {
253  auto offset = get_offset(coords);
254  // cppcheck-suppress "returnTempReference"
255  return m_container->at(offset);
256 }
257 
258 template <typename T>
259 T& NdArray<T>::at(const std::vector<size_t>& coords, const std::string& attr) {
260  auto offset = get_offset(coords, attr);
261  // cppcheck-suppress "returnTempReference"
262  return m_container->at(offset);
263 }
264 
265 template <typename T>
266 const T& NdArray<T>::at(const std::vector<size_t>& coords, const std::string& attr) const {
267  auto offset = get_offset(coords, attr);
268  // cppcheck-suppress "returnTempReference"
269  return m_container->at(offset);
270 }
271 
272 template <typename T>
273 template <typename... D>
274 T& NdArray<T>::at(size_t i, D... rest) {
275  return at_helper(m_offset, 0, i, rest...);
276 }
277 
278 template <typename T>
279 template <typename... D>
280 const T& NdArray<T>::at(size_t i, D... rest) const {
281  return at_helper(m_offset, 0, i, rest...);
282 }
283 
284 template <typename T>
285 auto NdArray<T>::begin() -> iterator {
286  return iterator{m_container.get(), m_offset, m_shape, m_stride_size, 0};
287 }
288 
289 template <typename T>
290 auto NdArray<T>::end() -> iterator {
291  return iterator{m_container.get(), m_offset, m_shape, m_stride_size, size()};
292 }
293 
294 template <typename T>
295 auto NdArray<T>::begin() const -> const_iterator {
296  return const_iterator{m_container.get(), m_offset, m_shape, m_stride_size, 0};
297 }
298 
299 template <typename T>
300 auto NdArray<T>::end() const -> const_iterator {
301  return const_iterator{m_container.get(), m_offset, m_shape, m_stride_size, size()};
302 }
303 
304 template <typename T>
305 size_t NdArray<T>::size() const {
306  return m_size;
307 }
308 
309 template <typename T>
310 bool NdArray<T>::operator==(const self_type& b) const {
311  if (shape() != b.shape())
312  return false;
313  for (auto ai = begin(), bi = b.begin(); ai != end() && bi != b.end(); ++ai, ++bi) {
314  if (*ai != *bi)
315  return false;
316  }
317  return true;
318 }
319 
320 template <typename T>
321 bool NdArray<T>::operator!=(const self_type& b) const {
322  return !(*this == b);
323 }
324 
325 template <typename T>
326 const std::vector<std::string>& NdArray<T>::attributes() const {
327  return m_attr_names;
328 }
329 
330 template <typename T>
331 auto NdArray<T>::concatenate(const self_type& other) -> self_type& {
332  // Verify dimensionality
333  if (m_shape.size() != other.m_shape.size()) {
334  throw std::length_error("Can not concatenate arrays with different dimensionality");
335  }
336  for (size_t i = 1; i < m_shape.size(); ++i) {
337  if (m_shape[i] != other.m_shape[i])
338  throw std::length_error("The size of all axis except for the first one must match");
339  }
340 
341  // New shape
342  auto old_size = m_container->size();
343  auto new_shape = m_shape;
344  new_shape[0] += other.m_shape[0];
345 
346  // Resize container
347  m_container->resize(new_shape);
348 
349  // Copy to the end
350  std::copy(std::begin(other), std::end(other), m_container->m_data_ptr + old_size);
351  // Done!
352  m_shape = new_shape;
353  return *this;
354 }
355 
356 template <typename T>
357 NdArray<T>::NdArray(std::shared_ptr<ContainerInterface> container, size_t offset, const std::vector<size_t>& shape_,
358  const std::vector<size_t>& stride, const std::vector<std::string>& attr_names)
359  : m_offset(offset)
360  , m_shape(shape_)
361  , m_stride_size(stride)
362  , m_attr_names(attr_names)
363  , m_size(std::accumulate(m_shape.begin(), m_shape.end(), 1, std::multiplies<size_t>()))
364  , m_container(container) {}
365 
366 template <typename T>
367 auto NdArray<T>::slice(size_t i) -> self_type {
368  if (m_shape.size() <= 1) {
369  throw std::out_of_range("Can not slice a one dimensional array");
370  }
371  std::vector<std::string> attrs;
372  if (!m_attr_names.empty()) {
373  attrs.resize(m_attr_names.size());
374  std::copy(m_attr_names.begin(), m_attr_names.end(), attrs.begin());
375  }
376  if (i >= m_shape[0]) {
377  throw std::out_of_range("Axis 0 out of range");
378  }
379  auto offset = m_offset + i * m_stride_size[0];
380  std::vector<size_t> stride_(m_stride_size.begin() + 1, m_stride_size.end());
381  std::vector<size_t> shape_(m_shape.begin() + 1, m_shape.end());
382  return NdArray(m_container, offset, std::move(shape_), std::move(stride_), std::move(attrs));
383 }
384 
385 template <typename T>
386 auto NdArray<T>::slice(size_t i) const -> const self_type {
387  return const_cast<NdArray<T>*>(this)->slice(i);
388 }
389 
390 template <typename T>
391 auto NdArray<T>::rslice(size_t i) -> self_type {
392  if (m_shape.size() <= 1) {
393  throw std::out_of_range("Can not slice a one dimensional array");
394  }
395  if (!m_attr_names.empty()) {
396  throw std::invalid_argument("Can not slice on the last axis for arrays with attribute names");
397  }
398  if (i >= m_shape.back()) {
399  throw std::out_of_range("Axis -1 out of range");
400  }
401  auto offset = m_offset + i * m_stride_size.back();
402  std::vector<size_t> strides_(m_stride_size.begin(), m_stride_size.end() - 1);
403  std::vector<size_t> shape_(m_shape.begin(), m_shape.end() - 1);
404  return NdArray(m_container, offset, std::move(shape_), std::move(strides_), m_attr_names);
405 }
406 
407 template <typename T>
408 auto NdArray<T>::rslice(size_t i) const -> const self_type {
409  return const_cast<NdArray<T>*>(this)->rslice(i);
410 }
411 
412 template <typename T>
413 size_t NdArray<T>::get_offset(const std::vector<size_t>& coords) const {
414  if (coords.size() != m_shape.size()) {
415  throw std::out_of_range("Invalid number of coordinates, got " + std::to_string(coords.size()) + ", expected " +
416  std::to_string(m_shape.size()));
417  }
418 
419  size_t offset = m_offset;
420  for (size_t i = 0; i < coords.size(); ++i) {
421  if (coords[i] >= m_shape[i]) {
422  throw std::out_of_range(std::to_string(coords[i]) + " >= " + std::to_string(m_shape[i]) + " for axis " +
423  std::to_string(i));
424  }
425  offset += coords[i] * m_stride_size[i];
426  }
427 
428  assert(offset < m_container->size());
429  return offset;
430 }
431 
432 template <typename T>
433 size_t NdArray<T>::get_attr_offset(const std::string& attr) const {
434  auto i = std::find(m_attr_names.begin(), m_attr_names.end(), attr);
435  if (i == m_attr_names.end())
436  throw std::out_of_range(attr);
437  return i - m_attr_names.begin();
438 }
439 
440 template <typename T>
441 void NdArray<T>::update_strides() {
442  m_stride_size.resize(m_shape.size());
443 
444  size_t acc = 1;
445  for (size_t i = m_stride_size.size(); i > 0; --i) {
446  m_stride_size[i - 1] = acc;
447  acc *= m_shape[i - 1];
448  }
449 }
450 
451 /**
452  * Helper to expand at with a variable number of arguments
453  */
454 template <typename T>
455 template <typename... D>
456 T& NdArray<T>::at_helper(size_t offset_acc, size_t axis, size_t i, D... rest) {
457  assert(axis <= m_shape.size() && i < m_shape[axis]);
458  offset_acc += i * m_stride_size[axis];
459  return at_helper(offset_acc, ++axis, rest...);
460 }
461 
462 template <typename T>
463 T& NdArray<T>::at_helper(size_t offset_acc, ELEMENTS_UNUSED size_t axis) {
464  assert(axis == m_shape.size());
465  assert(offset_acc < m_container->size());
466  return m_container->at(offset_acc);
467 }
468 
469 template <typename T>
470 T& NdArray<T>::at_helper(size_t offset_acc, ELEMENTS_UNUSED size_t axis, const std::string& attr) {
471  offset_acc += get_attr_offset(attr);
472  assert(axis == m_shape.size() - 1);
473  assert(offset_acc < m_container->size());
474  return m_container->at(offset_acc);
475 }
476 
477 template <typename T>
478 template <typename... D>
479 const T& NdArray<T>::at_helper(size_t offset_acc, size_t axis, size_t i, D... rest) const {
480  assert(axis <= m_shape.size() && i < m_shape[axis]);
481  offset_acc += i * m_stride_size[axis];
482  return at_helper(offset_acc, ++axis, rest...);
483 }
484 
485 template <typename T>
486 const T& NdArray<T>::at_helper(size_t offset_acc, ELEMENTS_UNUSED size_t axis) const {
487  assert(axis == m_shape.size());
488  assert(offset_acc < m_container->size());
489  return m_container->at(offset_acc);
490 }
491 
492 template <typename T>
493 const T& NdArray<T>::at_helper(size_t offset_acc, ELEMENTS_UNUSED size_t axis, const std::string& attr) const {
494  offset_acc += get_attr_offset(attr);
495  assert(axis == m_shape.size() - 1);
496  assert(offset_acc < m_container->size());
497  return m_container->at(offset_acc);
498 }
499 
500 template <typename T>
501 template <typename... D>
502 auto NdArray<T>::reshape_helper(std::vector<size_t>& acc, size_t i, D... rest) -> self_type& {
503  acc.push_back(i);
504  return reshape_helper(acc, rest...);
505 }
506 
507 template <typename T>
508 auto NdArray<T>::reshape_helper(std::vector<size_t>& acc) -> self_type& {
509  return reshape(acc);
510 }
511 
512 template <typename T>
513 std::ostream& operator<<(std::ostream& out, const NdArray<T>& ndarray) {
514  auto shape = ndarray.shape();
515 
516  if (ndarray.size()) {
517  out << "<";
518  size_t i;
519  for (i = 0; i < shape.size() - 1; ++i) {
520  out << shape[i] << ",";
521  }
522  out << shape[i] << ">";
523 
524  auto iter = ndarray.begin(), end_iter = ndarray.end() - 1;
525  for (; iter != end_iter; ++iter) {
526  out << *iter << ",";
527  }
528  out << *iter;
529  }
530  return out;
531 }
532 
533 } // end of namespace NdArray
534 } // end of namespace Euclid
535 
536 #endif // NDARRAY_IMPL