2 * Copyright (C) 2012-2021 Euclid Science Ground Segment
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)
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
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
21 #include <ElementsKernel/Unused.h>
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,
33 : m_container_ptr(container_ptr)
35 , m_row_size(std::accumulate(shape.begin(), shape.end(), 1, std::multiplies<size_t>()))
36 , m_stride{strides.back()}
41 NdArray<T>::Iterator<Const>::Iterator(ContainerInterface* container_ptr, size_t offset, size_t row_size, size_t stride,
43 : m_container_ptr(container_ptr), m_offset(offset), m_row_size(row_size), m_stride(stride), m_i(start) {}
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}
56 auto NdArray<T>::Iterator<Const>::operator++() -> Iterator& {
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};
69 bool NdArray<T>::Iterator<Const>::operator==(const Iterator& other) const {
70 return std::memcmp(this, &other, sizeof(*this)) == 0;
75 bool NdArray<T>::Iterator<Const>::operator!=(const Iterator& other) const {
76 return !(*this == other);
81 auto NdArray<T>::Iterator<Const>::operator*() -> value_t& {
87 auto NdArray<T>::Iterator<Const>::operator*() const -> value_t {
93 auto NdArray<T>::Iterator<Const>::operator+=(size_t n) -> Iterator& {
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};
104 template <typename T>
105 template <bool Const>
106 auto NdArray<T>::Iterator<Const>::operator-=(size_t n) -> Iterator& {
112 template <typename T>
113 template <bool Const>
114 auto NdArray<T>::Iterator<Const>::operator-(size_t n) -> Iterator {
116 return Iterator{m_container_ptr, m_offset, m_row_size, m_stride, m_i - n};
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;
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);
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);
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;
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;
152 template <typename T>
153 NdArray<T>::NdArray(const std::vector<size_t>& 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)) {
161 template <typename T>
162 template <template <class...> class Container>
163 NdArray<T>::NdArray(const std::vector<size_t>& shape_, const Container<T>& data)
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");
174 template <typename T>
175 template <template <class...> class Container>
176 NdArray<T>::NdArray(const std::vector<size_t>& shape_, Container<T>&& data)
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");
187 template <typename T>
188 template <typename II>
189 NdArray<T>::NdArray(const std::vector<size_t>& shape_, II ibegin, II iend)
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");
200 template <typename T>
201 NdArray<T>::NdArray(const self_type* other)
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()) {
210 inline std::vector<size_t> appendAttrShape(std::vector<size_t> shape, size_t append) {
212 shape.push_back(append);
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;
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");
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");
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...);
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);
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);
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);
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);
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...);
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...);
284 template <typename T>
285 auto NdArray<T>::begin() -> iterator {
286 return iterator{m_container.get(), m_offset, m_shape, m_stride_size, 0};
289 template <typename T>
290 auto NdArray<T>::end() -> iterator {
291 return iterator{m_container.get(), m_offset, m_shape, m_stride_size, size()};
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};
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()};
304 template <typename T>
305 size_t NdArray<T>::size() const {
309 template <typename T>
310 bool NdArray<T>::operator==(const self_type& b) const {
311 if (shape() != b.shape())
313 for (auto ai = begin(), bi = b.begin(); ai != end() && bi != b.end(); ++ai, ++bi) {
320 template <typename T>
321 bool NdArray<T>::operator!=(const self_type& b) const {
322 return !(*this == b);
325 template <typename T>
326 const std::vector<std::string>& NdArray<T>::attributes() const {
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");
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");
342 auto old_size = m_container->size();
343 auto new_shape = m_shape;
344 new_shape[0] += other.m_shape[0];
347 m_container->resize(new_shape);
350 std::copy(std::begin(other), std::end(other), m_container->m_data_ptr + old_size);
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)
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) {}
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");
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());
376 if (i >= m_shape[0]) {
377 throw std::out_of_range("Axis 0 out of range");
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));
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);
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");
395 if (!m_attr_names.empty()) {
396 throw std::invalid_argument("Can not slice on the last axis for arrays with attribute names");
398 if (i >= m_shape.back()) {
399 throw std::out_of_range("Axis -1 out of range");
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);
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);
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()));
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 " +
425 offset += coords[i] * m_stride_size[i];
428 assert(offset < m_container->size());
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();
440 template <typename T>
441 void NdArray<T>::update_strides() {
442 m_stride_size.resize(m_shape.size());
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];
452 * Helper to expand at with a variable number of arguments
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...);
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);
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);
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...);
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);
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);
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& {
504 return reshape_helper(acc, rest...);
507 template <typename T>
508 auto NdArray<T>::reshape_helper(std::vector<size_t>& acc) -> self_type& {
512 template <typename T>
513 std::ostream& operator<<(std::ostream& out, const NdArray<T>& ndarray) {
514 auto shape = ndarray.shape();
516 if (ndarray.size()) {
519 for (i = 0; i < shape.size() - 1; ++i) {
520 out << shape[i] << ",";
522 out << shape[i] << ">";
524 auto iter = ndarray.begin(), end_iter = ndarray.end() - 1;
525 for (; iter != end_iter; ++iter) {
533 } // end of namespace NdArray
534 } // end of namespace Euclid
536 #endif // NDARRAY_IMPL