Alexandria 2.25.0
SDC-CH common library for the Euclid project
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
26namespace Euclid {
27namespace NdArray {
28
29template <typename T>
30template <bool Const>
31NdArray<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
40template <typename T>
41template <bool Const>
42NdArray<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
46template <typename T>
47template <bool Const>
48NdArray<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
55template <typename T>
56template <bool Const>
57auto NdArray<T>::Iterator<Const>::operator++() -> Iterator& {
58 ++m_i;
59 return *this;
60}
61
62template <typename T>
63template <bool Const>
64auto 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
68template <typename T>
69template <bool Const>
70bool NdArray<T>::Iterator<Const>::operator==(const Iterator& other) const {
71 return std::memcmp(this, &other, sizeof(*this)) == 0;
72}
73
74template <typename T>
75template <bool Const>
76bool NdArray<T>::Iterator<Const>::operator!=(const Iterator& other) const {
77 return !(*this == other);
78}
79
80template <typename T>
81template <bool Const>
82auto NdArray<T>::Iterator<Const>::operator*() -> value_t& {
83 return operator[](0);
84}
85
86template <typename T>
87template <bool Const>
88auto NdArray<T>::Iterator<Const>::operator*() const -> value_t {
89 return operator[](0);
90}
91
92template <typename T>
93template <bool Const>
94auto NdArray<T>::Iterator<Const>::operator+=(size_t n) -> Iterator& {
95 m_i += n;
96 return *this;
97}
98
99template <typename T>
100template <bool Const>
101auto 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
105template <typename T>
106template <bool Const>
107auto NdArray<T>::Iterator<Const>::operator-=(size_t n) -> Iterator& {
108 assert(n <= m_i);
109 m_i -= n;
110 return *this;
111}
112
113template <typename T>
114template <bool Const>
115auto 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
120template <typename T>
121template <bool Const>
122auto 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
127template <typename T>
128template <bool Const>
129auto 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
133template <typename T>
134template <bool Const>
135auto 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
139template <typename T>
140template <bool Const>
141bool 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
146template <typename T>
147template <bool Const>
148bool 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
153template <typename T>
154NdArray<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
162template <typename T>
163template <template <class...> class Container>
164NdArray<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
175template <typename T>
176template <template <class...> class Container>
177NdArray<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
188template <typename T>
189template <template <class...> class Container>
190NdArray<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
205template <typename T>
206template <typename II>
207NdArray<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
218template <typename T>
219NdArray<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
228inline 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
234template <typename T>
235template <typename... Args>
236NdArray<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
241template <typename T>
242auto 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
255template <typename T>
256template <typename... D>
257auto 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
262template <typename T>
263T& 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
269template <typename T>
270const 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
276template <typename T>
277T& 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
283template <typename T>
284const 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
290template <typename T>
291template <typename... D>
292T& NdArray<T>::at(size_t i, D... rest) {
293 return at_helper(m_offset, 0, i, rest...);
294}
295
296template <typename T>
297template <typename... D>
298const T& NdArray<T>::at(size_t i, D... rest) const {
299 return at_helper(m_offset, 0, i, rest...);
300}
301
302template <typename T>
303auto NdArray<T>::begin() -> iterator {
304 return iterator{m_container.get(), m_offset, m_shape, m_stride_size, 0};
305}
306
307template <typename T>
308auto NdArray<T>::end() -> iterator {
309 return iterator{m_container.get(), m_offset, m_shape, m_stride_size, size()};
310}
311
312template <typename T>
313auto NdArray<T>::begin() const -> const_iterator {
314 return const_iterator{m_container.get(), m_offset, m_shape, m_stride_size, 0};
315}
316
317template <typename T>
318auto NdArray<T>::end() const -> const_iterator {
319 return const_iterator{m_container.get(), m_offset, m_shape, m_stride_size, size()};
320}
321
322template <typename T>
323size_t NdArray<T>::size() const {
324 return m_size;
325}
326
327template <typename T>
328bool 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
338template <typename T>
339bool NdArray<T>::operator!=(const self_type& b) const {
340 return !(*this == b);
341}
342
343template <typename T>
344const std::vector<std::string>& NdArray<T>::attributes() const {
345 return m_attr_names;
346}
347
348template <typename T>
349auto 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
375template <typename T>
376NdArray<T>::NdArray(std::shared_ptr<ContainerInterface> container, size_t offset, std::vector<size_t> shape_,
377 std::vector<size_t> stride, std::vector<std::string> attr_names)
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
386template <typename T>
387auto 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 }
391 std::vector<std::string> attrs;
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
405template <typename T>
406auto NdArray<T>::slice(size_t i) const -> const self_type {
407 return const_cast<NdArray<T>*>(this)->slice(i);
408}
409
410template <typename T>
411auto 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
427template <typename T>
428auto NdArray<T>::rslice(size_t i) const -> const self_type {
429 return const_cast<NdArray<T>*>(this)->rslice(i);
430}
431
432template <typename T>
433void NdArray<T>::next_slice() {
434 m_offset += m_total_stride;
435}
436
437template <typename T>
438size_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
457template <typename T>
458size_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
465template <typename T>
466void 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
478/**
479 * Helper to expand at with a variable number of arguments
480 */
481template <typename T>
482template <typename... D>
483T& 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
489template <typename T>
490T& 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
496template <typename T>
497T& 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
504template <typename T>
505template <typename... D>
506const 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
512template <typename T>
513const 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
520template <typename T>
521const 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
529template <typename T>
530template <typename... D>
531auto 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
536template <typename T>
537auto NdArray<T>::reshape_helper(std::vector<size_t>& acc) -> self_type& {
538 return reshape(acc);
539}
540
541template <typename T>
542std::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