Alexandria  2.18
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 <algorithm>
22 
23 namespace Euclid {
24 namespace NdArray {
25 
26 template <typename T>
27 template <bool Const>
28 NdArray<T>::Iterator<Const>::Iterator(ContainerInterface* container_ptr, size_t offset)
29  : m_container_ptr{container_ptr}, m_offset{offset} {}
30 
31 template <typename T>
32 template <bool Const>
33 NdArray<T>::Iterator<Const>::Iterator(const Iterator<false>& other)
34  : m_container_ptr{other.m_container_ptr}, m_offset{other.m_offset} {}
35 
36 template <typename T>
37 template <bool Const>
38 auto NdArray<T>::Iterator<Const>::operator++() -> Iterator& {
39  ++m_offset;
40  return *this;
41 }
42 
43 template <typename T>
44 template <bool Const>
45 auto NdArray<T>::Iterator<Const>::operator++(int) -> Iterator {
46  return Iterator{m_container_ptr, m_offset + 1};
47 }
48 
49 template <typename T>
50 template <bool Const>
51 bool NdArray<T>::Iterator<Const>::operator==(const Iterator& other) const {
52  return m_container_ptr == other.m_container_ptr && m_offset == other.m_offset;
53 }
54 
55 template <typename T>
56 template <bool Const>
57 bool NdArray<T>::Iterator<Const>::operator!=(const Iterator& other) const {
58  return m_container_ptr != other.m_container_ptr || m_offset != other.m_offset;
59 }
60 
61 template <typename T>
62 template <bool Const>
63 auto NdArray<T>::Iterator<Const>::operator*() -> value_t& {
64  return m_container_ptr->at(m_offset);
65 }
66 
67 template <typename T>
68 template <bool Const>
69 auto NdArray<T>::Iterator<Const>::operator*() const -> value_t {
70  return m_container_ptr->at(m_offset);
71 }
72 
73 template <typename T>
74 template <bool Const>
75 auto NdArray<T>::Iterator<Const>::operator+=(size_t n) -> Iterator& {
76  m_offset += n;
77  return *this;
78 }
79 
80 template <typename T>
81 template <bool Const>
82 auto NdArray<T>::Iterator<Const>::operator+(size_t n) -> Iterator {
83  return Iterator{m_container_ptr, m_offset + n};
84 }
85 
86 template <typename T>
87 template <bool Const>
88 auto NdArray<T>::Iterator<Const>::operator-=(size_t n) -> Iterator& {
89  assert(n <= m_offset);
90  m_offset -= n;
91  return *this;
92 }
93 
94 template <typename T>
95 template <bool Const>
96 auto NdArray<T>::Iterator<Const>::operator-(size_t n) -> Iterator {
97  assert(n <= m_offset);
98  return Iterator{m_container_ptr, m_offset - n};
99 }
100 
101 template <typename T>
102 template <bool Const>
103 auto NdArray<T>::Iterator<Const>::operator-(const Iterator& other) -> difference_type {
104  assert(m_container_ptr == other.m_container_ptr);
105  return m_offset - other.m_offset;
106 }
107 
108 template <typename T>
109 template <bool Const>
110 auto NdArray<T>::Iterator<Const>::operator[](size_t i) -> value_t& {
111  return m_container_ptr->at(m_offset + i);
112 }
113 
114 template <typename T>
115 template <bool Const>
116 bool NdArray<T>::Iterator<Const>::operator<(const Iterator& other) {
117  assert(m_container_ptr == other.m_container_ptr);
118  return m_offset < other.m_offset;
119 }
120 
121 template <typename T>
122 template <bool Const>
123 bool NdArray<T>::Iterator<Const>::operator>(const Iterator& other) {
124  assert(m_container_ptr == other.m_container_ptr);
125  return m_offset > other.m_offset;
126 }
127 
128 template <typename T>
129 NdArray<T>::NdArray(const std::vector<size_t>& shape_)
130  : m_shape{shape_}
131  , m_size{std::accumulate(m_shape.begin(), m_shape.end(), 1u, std::multiplies<size_t>())}
132  , m_container(new ContainerWrapper<std::vector>(m_size)) {
133  update_strides();
134 }
135 
136 template <typename T>
137 template <template <class...> class Container>
138 NdArray<T>::NdArray(const std::vector<size_t>& shape_, const Container<T>& data)
139  : m_shape{shape_}
140  , m_size{std::accumulate(m_shape.begin(), m_shape.end(), 1u, std::multiplies<size_t>())}
141  , m_container{new ContainerWrapper<Container>(data)} {
142  if (m_size != m_container->size()) {
143  throw std::invalid_argument("Data size does not match the shape");
144  }
145  update_strides();
146 }
147 
148 template <typename T>
149 template <template <class...> class Container>
150 NdArray<T>::NdArray(const std::vector<size_t>& shape_, Container<T>&& data)
151  : m_shape{shape_}
152  , m_size{std::accumulate(m_shape.begin(), m_shape.end(), 1u, std::multiplies<size_t>())}
153  , m_container{new ContainerWrapper<Container>(std::move(data))} {
154  if (m_size != m_container->size()) {
155  throw std::invalid_argument("Data size does not match the shape");
156  }
157  update_strides();
158 }
159 
160 template <typename T>
161 template <typename II>
162 NdArray<T>::NdArray(const std::vector<size_t>& shape_, II ibegin, II iend)
163  : m_shape{shape_}
164  , m_size{std::accumulate(m_shape.begin(), m_shape.end(), 1u, std::multiplies<size_t>())}
165  , m_container{new ContainerWrapper<std::vector>(ibegin, iend)} {
166  if (m_size != m_container->size()) {
167  throw std::invalid_argument("Data size does not match the shape");
168  }
169  update_strides();
170 }
171 
172 template <typename T>
173 NdArray<T>::NdArray(const self_type* other)
174  : m_shape{other->m_shape}
175  , m_attr_names{other->m_attr_names}
176  , m_size{std::accumulate(m_shape.begin(), m_shape.end(), 1u, std::multiplies<size_t>())}
177  , m_container{other->m_container->copy()} {
178  update_strides();
179 }
180 
181 inline std::vector<size_t> appendAttrShape(std::vector<size_t> shape, size_t append) {
182  if (append)
183  shape.push_back(append);
184  return shape;
185 }
186 
187 template <typename T>
188 template <typename... Args>
189 NdArray<T>::NdArray(const std::vector<size_t>& shape_, const std::vector<std::string>& attr_names, Args&&... args)
190  : NdArray{appendAttrShape(shape_, attr_names.size()), std::forward<Args>(args)...} {
191  m_attr_names = attr_names;
192 }
193 
194 template <typename T>
195 auto NdArray<T>::reshape(const std::vector<size_t> new_shape) -> self_type& {
196  if (!m_attr_names.empty())
197  throw std::invalid_argument("Can not reshape arrays with attribute names");
198 
199  size_t new_size = std::accumulate(new_shape.begin(), new_shape.end(), 1, std::multiplies<size_t>());
200  if (new_size != m_size) {
201  throw std::range_error("New shape does not match the number of contained elements");
202  }
203  m_shape = new_shape;
204  update_strides();
205  return *this;
206 }
207 
208 template <typename T>
209 template <typename... D>
210 auto NdArray<T>::reshape(size_t i, D... rest) -> self_type& {
211  std::vector<size_t> acc{i};
212  return reshape_helper(acc, rest...);
213 }
214 
215 template <typename T>
216 T& NdArray<T>::at(const std::vector<size_t>& coords) {
217  auto offset = get_offset(coords);
218  return m_container->at(offset);
219 }
220 
221 template <typename T>
222 const T& NdArray<T>::at(const std::vector<size_t>& coords) const {
223  auto offset = get_offset(coords);
224  return m_container->at(offset);
225 }
226 
227 template <typename T>
228 T& NdArray<T>::at(const std::vector<size_t>& coords, const std::string& attr) {
229  auto offset = get_offset(coords, attr);
230  return m_container->at(offset);
231 }
232 
233 template <typename T>
234 const T& NdArray<T>::at(const std::vector<size_t>& coords, const std::string& attr) const {
235  auto offset = get_offset(coords, attr);
236  return m_container->at(offset);
237 }
238 
239 template <typename T>
240 template <typename... D>
241 T& NdArray<T>::at(size_t i, D... rest) {
242  std::vector<size_t> acc{i};
243  return at_helper(acc, rest...);
244 }
245 
246 template <typename T>
247 template <typename... D>
248 const T& NdArray<T>::at(size_t i, D... rest) const {
249  std::vector<size_t> acc{i};
250  return at_helper(acc, rest...);
251 }
252 
253 template <typename T>
254 auto NdArray<T>::begin() -> iterator {
255  return iterator{m_container.get(), 0};
256 }
257 
258 template <typename T>
259 auto NdArray<T>::end() -> iterator {
260  return iterator{m_container.get(), m_size};
261 }
262 
263 template <typename T>
264 auto NdArray<T>::begin() const -> const_iterator {
265  return const_iterator{m_container.get(), 0};
266 }
267 
268 template <typename T>
269 auto NdArray<T>::end() const -> const_iterator {
270  return const_iterator{m_container.get(), m_size};
271 }
272 
273 template <typename T>
274 size_t NdArray<T>::size() const {
275  return m_size;
276 }
277 
278 template <typename T>
279 bool NdArray<T>::operator==(const self_type& b) const {
280  if (shape() != b.shape())
281  return false;
282  for (auto ai = begin(), bi = b.begin(); ai != end() && bi != b.end(); ++ai, ++bi) {
283  if (*ai != *bi)
284  return false;
285  }
286  return true;
287 }
288 
289 template <typename T>
290 bool NdArray<T>::operator!=(const self_type& b) const {
291  return !(*this == b);
292 }
293 
294 template <typename T>
295 const std::vector<std::string>& NdArray<T>::attributes() const {
296  return m_attr_names;
297 }
298 
299 template <typename T>
300 auto NdArray<T>::concatenate(const self_type& other) -> self_type& {
301  // Verify dimensionality
302  if (m_shape.size() != other.m_shape.size()) {
303  throw std::length_error("Can not concatenate arrays with different dimensionality");
304  }
305  for (size_t i = 1; i < m_shape.size(); ++i) {
306  if (m_shape[i] != other.m_shape[i])
307  throw std::length_error("The size of all axis except for the first one must match");
308  }
309 
310  // New shape
311  auto old_size = m_container->size();
312  auto new_shape = m_shape;
313  new_shape[0] += other.m_shape[0];
314 
315  // Resize container
316  m_container->resize(new_shape);
317 
318  // Copy to the end
319  std::copy(std::begin(other), std::end(other), m_container->m_data_ptr + old_size);
320  // Done!
321  m_shape = new_shape;
322  return *this;
323 }
324 
325 template <typename T>
326 size_t NdArray<T>::get_offset(const std::vector<size_t>& coords) const {
327  if (coords.size() != m_shape.size()) {
328  throw std::out_of_range("Invalid number of coordinates, got " + std::to_string(coords.size()) + ", expected " +
329  std::to_string(m_shape.size()));
330  }
331 
332  size_t offset = 0;
333  for (size_t i = 0; i < coords.size(); ++i) {
334  if (coords[i] >= m_shape[i]) {
335  throw std::out_of_range(std::to_string(coords[i]) + " >= " + std::to_string(m_shape[i]) + " for axis " + std::to_string(i));
336  }
337  offset += coords[i] * m_stride_size[i];
338  }
339 
340  assert(offset < m_container->size());
341  return offset;
342 }
343 
344 template <typename T>
345 size_t NdArray<T>::get_offset(std::vector<size_t> coords, const std::string& attr) const {
346  auto i = std::find(m_attr_names.begin(), m_attr_names.end(), attr);
347  if (i == m_attr_names.end())
348  throw std::out_of_range(attr);
349  auto last = i - m_attr_names.begin();
350  coords.emplace_back(last);
351  return get_offset(coords);
352 }
353 
354 template <typename T>
355 void NdArray<T>::update_strides() {
356  m_stride_size.resize(m_shape.size());
357 
358  size_t acc = 1;
359  for (size_t i = m_stride_size.size(); i > 0; --i) {
360  m_stride_size[i - 1] = acc;
361  acc *= m_shape[i - 1];
362  }
363 }
364 
365 /**
366  * Helper to expand at with a variable number of arguments
367  */
368 template <typename T>
369 template <typename... D>
370 T& NdArray<T>::at_helper(std::vector<size_t>& acc, size_t i, D... rest) {
371  acc.push_back(i);
372  return at_helper(acc, rest...);
373 }
374 
375 template <typename T>
376 T& NdArray<T>::at_helper(std::vector<size_t>& acc) {
377  return at(acc);
378 }
379 
380 template <typename T>
381 T& NdArray<T>::at_helper(std::vector<size_t>& acc, const std::string& attr) {
382  return at(acc, attr);
383 }
384 
385 template <typename T>
386 template <typename... D>
387 const T& NdArray<T>::at_helper(std::vector<size_t>& acc, size_t i, D... rest) const {
388  acc.push_back(i);
389  return at_helper(acc, rest...);
390 }
391 
392 template <typename T>
393 const T& NdArray<T>::at_helper(std::vector<size_t>& acc) const {
394  return at(acc);
395 }
396 
397 template <typename T>
398 template <typename... D>
399 auto NdArray<T>::reshape_helper(std::vector<size_t>& acc, size_t i, D... rest) -> self_type& {
400  acc.push_back(i);
401  return reshape_helper(acc, rest...);
402 }
403 
404 template <typename T>
405 auto NdArray<T>::reshape_helper(std::vector<size_t>& acc) -> self_type& {
406  return reshape(acc);
407 }
408 
409 template <typename T>
410 std::ostream& operator<<(std::ostream& out, const NdArray<T>& ndarray) {
411  auto shape = ndarray.shape();
412 
413  if (ndarray.size()) {
414  out << "<";
415  size_t i;
416  for (i = 0; i < shape.size() - 1; ++i) {
417  out << shape[i] << ",";
418  }
419  out << shape[i] << ">";
420 
421  auto iter = ndarray.begin(), end_iter = ndarray.end() - 1;
422  for (; iter != end_iter; ++iter) {
423  out << *iter << ",";
424  }
425  out << *iter;
426  }
427  return out;
428 }
429 
430 } // end of namespace NdArray
431 } // end of namespace Euclid
432 
433 #endif // NDARRAY_IMPL