Alexandria 2.31.0
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
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 <utility>
24
25namespace Euclid {
26namespace NdArray {
27
28template <typename T>
29template <bool Const>
30NdArray<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
39template <typename T>
40template <bool Const>
41NdArray<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
45template <typename T>
46template <bool Const>
47NdArray<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
54template <typename T>
55template <bool Const>
56auto NdArray<T>::Iterator<Const>::operator++() -> Iterator& {
57 ++m_i;
58 return *this;
59}
60
61template <typename T>
62template <bool Const>
63auto 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
67template <typename T>
68template <bool Const>
69bool 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
73template <typename T>
74template <bool Const>
75bool NdArray<T>::Iterator<Const>::operator!=(const Iterator& other) const {
76 return !(*this == other);
77}
78
79template <typename T>
80template <bool Const>
81auto NdArray<T>::Iterator<Const>::operator*() -> value_t& {
82 return operator[](0);
83}
84
85template <typename T>
86template <bool Const>
87auto NdArray<T>::Iterator<Const>::operator*() const -> value_t {
88 return operator[](0);
89}
90
91template <typename T>
92template <bool Const>
93auto NdArray<T>::Iterator<Const>::operator+=(size_t n) -> Iterator& {
94 m_i += n;
95 return *this;
96}
97
98template <typename T>
99template <bool Const>
100auto 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
104template <typename T>
105template <bool Const>
106auto NdArray<T>::Iterator<Const>::operator-=(size_t n) -> Iterator& {
107 assert(n <= m_i);
108 m_i -= n;
109 return *this;
110}
111
112template <typename T>
113template <bool Const>
114auto 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
119template <typename T>
120template <bool Const>
121auto 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
126template <typename T>
127template <bool Const>
128auto 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
132template <typename T>
133template <bool Const>
134auto 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
138template <typename T>
139template <bool Const>
140bool 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
145template <typename T>
146template <bool Const>
147bool 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
152template <typename T>
153NdArray<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
165template <typename T>
166template <template <class...> class Container>
167NdArray<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
181template <typename T>
182template <template <class...> class Container>
183NdArray<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
197template <typename T>
198template <template <class...> class Container>
199NdArray<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
215template <typename T>
216template <typename II>
217NdArray<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
231template <typename T>
232NdArray<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
236inline 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
242template <typename T>
243template <typename... Args>
244NdArray<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
249template <typename T>
250auto 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
263template <typename T>
264template <typename... D>
265auto 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
270template <typename T>
271const T& NdArray<T>::front() const {
272 return m_details_ptr->m_container->get(m_details_ptr->m_offset);
273}
274
275template <typename T>
276T& 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
282template <typename T>
283const 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
289template <typename T>
290T& 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
296template <typename T>
297const 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
303template <typename T>
304template <typename... D>
305T& NdArray<T>::at(size_t i, D... rest) {
306 return at_helper(m_details_ptr->m_offset, 0, i, rest...);
307}
308
309template <typename T>
310template <typename... D>
311const T& NdArray<T>::at(size_t i, D... rest) const {
312 return at_helper(m_details_ptr->m_offset, 0, i, rest...);
313}
314
315template <typename T>
316auto 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
321template <typename T>
322auto 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
327template <typename T>
328auto 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
333template <typename T>
334auto 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
339template <typename T>
340size_t NdArray<T>::size() const {
341 return m_details_ptr->m_size;
342}
343
344template <typename T>
345bool 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
355template <typename T>
356bool NdArray<T>::operator!=(const self_type& b) const {
357 return !(*this == b);
358}
359
360template <typename T>
361const std::vector<std::string>& NdArray<T>::attributes() const {
362 return m_details_ptr->m_attr_names;
363}
364
365template <typename T>
366auto 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
392template <typename T>
393NdArray<T>::NdArray(std::shared_ptr<ContainerInterface> container, size_t offset, std::vector<size_t> shape_,
394 std::vector<size_t> stride, std::vector<std::string> attr_names)
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
399template <typename T>
400NdArray<T>::NdArray(const self_type& other) : m_details_ptr(make_unique<Details>(*other.m_details_ptr)) {}
401
402template <typename T>
403NdArray<T>& NdArray<T>::operator=(const self_type& other) {
404 m_details_ptr = make_unique<Details>(*other.m_details_ptr);
405 return *this;
406}
407
408template <typename T>
409auto 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 }
413 std::vector<std::string> attrs;
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
427template <typename T>
428auto NdArray<T>::slice(size_t i) const -> const self_type {
429 return const_cast<NdArray<T>*>(this)->slice(i);
430}
431
432template <typename T>
433auto 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
450template <typename T>
451auto NdArray<T>::rslice(size_t i) const -> const self_type {
452 return const_cast<NdArray<T>*>(this)->rslice(i);
453}
454
455template <typename T>
456void NdArray<T>::next_slice() {
457 m_details_ptr->m_offset += m_details_ptr->m_total_stride;
458}
459
460template <typename T>
461size_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
480template <typename T>
481size_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
488template <typename T>
489void 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
501/**
502 * Helper to expand at with a variable number of arguments
503 */
504template <typename T>
505template <typename... D>
506T& 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
512template <typename T>
513T& 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
519template <typename T>
520T& 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
527template <typename T>
528template <typename... D>
529const 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
535template <typename T>
536const 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
543template <typename T>
544const 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
552template <typename T>
553template <typename... D>
554auto 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
559template <typename T>
560auto NdArray<T>::reshape_helper(std::vector<size_t>& acc) -> self_type& {
561 return reshape(acc);
562}
563
564template <typename T>
565std::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