BioDynaMo  v1.03.58-27764645
math_array.h
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 CERN & University of Surrey for the benefit of the
4 // BioDynaMo collaboration. All Rights Reserved.
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 //
9 // See the LICENSE file distributed with this work for details.
10 // See the NOTICE file distributed with this work for additional information
11 // regarding copyright ownership.
12 //
13 // -----------------------------------------------------------------------------
14 
15 #ifndef CORE_CONTAINER_MATH_ARRAY_H_
16 #define CORE_CONTAINER_MATH_ARRAY_H_
17 
18 #include <algorithm>
19 #include <cassert>
20 #include <cmath>
21 #include <numeric>
22 #include <ostream>
23 #include <stdexcept>
24 #include <utility>
25 
26 #include "core/real_t.h"
27 #include "core/util/log.h"
28 #include "core/util/root.h"
29 
30 namespace bdm {
31 
35 template <class T, std::size_t N>
36 class MathArray { // NOLINT
37  public:
40 #pragma omp simd
41  for (size_t i = 0; i < N; i++) {
42  data_[i] = T();
43  }
44  }
45 
49  constexpr MathArray(std::initializer_list<T> l) {
50  assert(l.size() <= N);
51  auto it = l.begin();
52  for (uint64_t i = 0; i < N; i++) {
53  data_[i] = *(it++);
54  }
55  for (uint64_t i = l.size(); i < N; i++) {
56  data_[i] = T();
57  }
58  }
59 
62  inline const T* data() const { return &data_[0]; } // NOLINT
63 
66  inline const size_t size() const { return N; } // NOLINT
67 
70  inline const bool empty() const { return N == 0; } // NOLINT
71 
76  T& operator[](size_t idx) { return data_[idx]; }
77 
81  const T& operator[](size_t idx) const { return data_[idx]; }
82 
88  T& at(size_t idx) noexcept(false) { // NOLINT
89  if (idx > size() || idx < 0) {
90  throw std::out_of_range("The index is out of range");
91  }
92  return data_[idx];
93  }
94 
95  const T* begin() const { return &(data_[0]); } // NOLINT
96 
97  const T* end() const { return &(data_[N]); } // NOLINT
98 
99  T* begin() { return &(data_[0]); } // NOLINT
100 
101  T* end() { return &(data_[N]); } // NOLINT
102 
105  T& front() { return *(this->begin()); } // NOLINT
106 
109  T& back() { // NOLINT
110  auto tmp = this->end();
111  tmp--;
112  return *tmp;
113  }
114 
118  MathArray& operator=(const MathArray& other) {
119  if (this != &other) {
120  assert(other.size() == N);
121  std::copy(other.data_, other.data_ + other.size(), data_);
122  }
123  return *this;
124  }
125 
129  bool operator==(const MathArray& other) const {
130  if (other.size() != N) {
131  return false;
132  }
133  for (size_t i = 0; i < N; i++) {
134  if (other[i] != data_[i]) {
135  return false;
136  }
137  }
138  return true;
139  }
140 
141  bool operator!=(const MathArray& other) const { return !operator==(other); }
142 
144 #pragma omp simd
145  for (size_t i = 0; i < N; i++) {
146  ++data_[i];
147  }
148  return *this;
149  }
150 
152  MathArray tmp(*this);
153  operator++();
154  return tmp;
155  }
156 
158 #pragma omp simd
159  for (size_t i = 0; i < N; i++) {
160  --data_[i];
161  }
162  return *this;
163  }
164 
166  MathArray tmp(*this);
167  operator--();
168  return tmp;
169  }
170 
172  assert(N == rhs.size());
173 #pragma omp simd
174  for (size_t i = 0; i < N; i++) {
175  data_[i] += rhs[i];
176  }
177  return *this;
178  }
179 
181  assert(size() == rhs.size());
182  MathArray tmp;
183 #pragma omp simd
184  for (size_t i = 0; i < N; i++) {
185  tmp[i] = data_[i] + rhs[i];
186  }
187  return tmp;
188  }
189 
190  const MathArray operator+(const MathArray& rhs) const {
191  assert(size() == rhs.size());
192  MathArray tmp;
193 #pragma omp simd
194  for (size_t i = 0; i < N; i++) {
195  tmp[i] = data_[i] + rhs[i];
196  }
197  return tmp;
198  }
199 
200  MathArray& operator+=(const T& rhs) {
201 #pragma omp simd
202  for (size_t i = 0; i < N; i++) {
203  data_[i] += rhs;
204  }
205  return *this;
206  }
207 
208  MathArray operator+(const T& rhs) {
209  MathArray tmp;
210 #pragma omp simd
211  for (size_t i = 0; i < N; i++) {
212  tmp[i] = data_[i] + rhs;
213  }
214  return tmp;
215  }
216 
218  assert(size() == rhs.size());
219 #pragma omp simd
220  for (size_t i = 0; i < N; i++) {
221  data_[i] -= rhs[i];
222  }
223  return *this;
224  }
225 
227  assert(size() == rhs.size());
228  MathArray tmp;
229 #pragma omp simd
230  for (size_t i = 0; i < N; i++) {
231  tmp[i] = data_[i] - rhs[i];
232  }
233  return tmp;
234  }
235 
236  const MathArray operator-(const MathArray& rhs) const {
237  assert(size() == rhs.size());
238  MathArray tmp;
239 #pragma omp simd
240  for (size_t i = 0; i < N; i++) {
241  tmp[i] = data_[i] - rhs[i];
242  }
243  return tmp;
244  }
245 
246  MathArray& operator-=(const T& rhs) {
247 #pragma omp simd
248  for (size_t i = 0; i < N; i++) {
249  data_[i] -= rhs;
250  }
251  return *this;
252  }
253 
254  MathArray operator-(const T& rhs) {
255  MathArray tmp;
256 #pragma omp simd
257  for (size_t i = 0; i < N; i++) {
258  tmp[i] = data_[i] - rhs;
259  }
260  return tmp;
261  }
262 
263  T& operator*=(const MathArray& rhs) = delete;
264 
265  T operator*(const MathArray& rhs) {
266  assert(size() == rhs.size());
267  T result = 0;
268 #pragma omp simd
269  for (size_t i = 0; i < N; i++) {
270  result += data_[i] * rhs[i];
271  }
272  return result;
273  }
274 
275  const T operator*(const MathArray& rhs) const {
276  assert(size() == rhs.size());
277  T result = 0;
278 #pragma omp simd
279  for (size_t i = 0; i < N; i++) {
280  result += data_[i] * rhs[i];
281  }
282  return result;
283  }
284 
285  MathArray& operator*=(const T& k) {
286 #pragma omp simd
287  for (size_t i = 0; i < N; i++) {
288  data_[i] *= k;
289  }
290  return *this;
291  }
292 
293  MathArray operator*(const T& k) {
294  MathArray tmp;
295 #pragma omp simd
296  for (size_t i = 0; i < N; i++) {
297  tmp[i] = data_[i] * k;
298  }
299  return tmp;
300  }
301 
302  const MathArray operator*(const T& k) const {
303  MathArray tmp;
304 #pragma omp simd
305  for (size_t i = 0; i < N; i++) {
306  tmp[i] = data_[i] * k;
307  }
308  return tmp;
309  }
310 
311  MathArray& operator/=(const T& k) {
312 #pragma omp simd
313  for (size_t i = 0; i < N; i++) {
314  data_[i] /= k;
315  }
316  return *this;
317  }
318 
319  MathArray operator/(const T& k) {
320  MathArray tmp;
321 #pragma omp simd
322  for (size_t i = 0; i < N; i++) {
323  tmp[i] = data_[i] / k;
324  }
325  return tmp;
326  }
327 
331  MathArray& fill(const T& k) { // NOLINT
332  std::fill(std::begin(data_), std::end(data_), k);
333  return *this;
334  }
335 
338  T Sum() const { return std::accumulate(begin(), end(), 0); }
339 
341  bool IsZero() const {
342  for (size_t i = 0; i < N; i++) {
343  if (data_[i] != 0) {
344  return false;
345  }
346  }
347  return true;
348  }
349 
352  T Norm() const {
353  T result = 0;
354 #pragma omp simd
355  for (size_t i = 0; i < N; i++) {
356  result += data_[i] * data_[i];
357  }
358  result = std::sqrt(result);
359 
360  return result;
361  }
362 
364  void Normalize() {
365  T norm = Norm();
366  Normalize(norm);
367  }
368 
373  void Normalize(T norm) {
374  if (norm == 0) {
375  Log::Fatal("MathArray::Normalize",
376  "You tried to normalize a zero vector. "
377  "This cannot be done. Exiting.");
378  }
379 #pragma omp simd
380  for (size_t i = 0; i < N; i++) {
381  data_[i] /= norm;
382  }
383  }
384 
387  MathArray normalized_array(*this);
388  normalized_array.Normalize();
389  return normalized_array;
390  }
391 
397  assert(rhs.size() == N);
398  MathArray tmp;
399 #pragma omp simd
400  for (size_t i = 0; i < N; ++i) {
401  tmp[i] = data_[i] * rhs[i];
402  }
403  return tmp;
404  }
405 
406  private:
407  T data_[N];
408  BDM_CLASS_DEF_NV(MathArray, 1); // NOLINT
409 };
410 
411 template <class T, std::size_t N>
412 std::ostream& operator<<(std::ostream& o, const MathArray<T, N>& arr) {
413  for (size_t i = 0; i < N; i++) {
414  o << arr[i];
415  if (i != N - 1) {
416  o << ", ";
417  }
418  }
419  return o;
420 }
421 
422 // Note: 1) We pass by value to allow for copy-elision and move semantics
423 // optimization.
424 // 2) We do return MathArray<T, N> because a const MathArray<T, N>
425 // prevents move semantics in C++11.
426 // see https://tinyurl.com/left-multiply
428 template <class T, std::size_t N>
429 MathArray<T, N> operator*(T const& scalar, MathArray<T, N> array) {
430  return array *= scalar;
431 }
432 
437 
442 
443 } // namespace bdm
444 
445 #endif // CORE_CONTAINER_MATH_ARRAY_H_
bdm::MathArray::operator[]
const T & operator[](size_t idx) const
Definition: math_array.h:81
bdm::MathArray::front
T & front()
Definition: math_array.h:105
bdm::MathArray::operator+
const MathArray operator+(const MathArray &rhs) const
Definition: math_array.h:190
bdm::MathArray::Sum
T Sum() const
Definition: math_array.h:338
bdm::MathArray::operator/
MathArray operator/(const T &k)
Definition: math_array.h:319
bdm::MathArray::operator*
const MathArray operator*(const T &k) const
Definition: math_array.h:302
bdm::MathArray::operator=
MathArray & operator=(const MathArray &other)
Definition: math_array.h:118
bdm::MathArray::data_
T data_[N]
Definition: math_array.h:407
bdm
Definition: agent.cc:39
bdm::MathArray::begin
T * begin()
Definition: math_array.h:99
bdm::MathArray::end
T * end()
Definition: math_array.h:101
bdm::MathArray::operator--
MathArray & operator--()
Definition: math_array.h:157
bdm::MathArray::operator-=
MathArray & operator-=(const T &rhs)
Definition: math_array.h:246
bdm::MathArray::operator*=
MathArray & operator*=(const T &k)
Definition: math_array.h:285
bdm::operator*
MathArray< T, N > operator*(T const &scalar, MathArray< T, N > array)
Template function to multiply array with scalar from the left.
Definition: math_array.h:429
bdm::MathArray::operator+
MathArray operator+(const MathArray &rhs)
Definition: math_array.h:180
bdm::MathArray::operator*
const T operator*(const MathArray &rhs) const
Definition: math_array.h:275
bdm::MathArray::data
const T * data() const
Definition: math_array.h:62
bdm::MathArray::operator==
bool operator==(const MathArray &other) const
Definition: math_array.h:129
bdm::MathArray::fill
MathArray & fill(const T &k)
Definition: math_array.h:331
bdm::MathArray::operator+=
MathArray & operator+=(const T &rhs)
Definition: math_array.h:200
bdm::MathArray::Normalize
void Normalize(T norm)
Definition: math_array.h:373
bdm::MathArray::operator-=
MathArray & operator-=(const MathArray &rhs)
Definition: math_array.h:217
bdm::MathArray::operator/=
MathArray & operator/=(const T &k)
Definition: math_array.h:311
bdm::MathArray::Norm
T Norm() const
Definition: math_array.h:352
bdm::MathArray::operator++
MathArray operator++(int)
Definition: math_array.h:151
bdm::MathArray::operator!=
bool operator!=(const MathArray &other) const
Definition: math_array.h:141
bdm::MathArray::operator-
MathArray operator-(const MathArray &rhs)
Definition: math_array.h:226
bdm::MathArray::begin
const T * begin() const
Definition: math_array.h:95
root.h
bdm::MathArray::operator[]
T & operator[](size_t idx)
Definition: math_array.h:76
bdm::MathArray::operator*
MathArray operator*(const T &k)
Definition: math_array.h:293
bdm::MathArray::BDM_CLASS_DEF_NV
BDM_CLASS_DEF_NV(MathArray, 1)
bdm::Log::Fatal
static void Fatal(const std::string &location, const Args &... parts)
Prints fatal error message.
Definition: log.h:115
log.h
bdm::MathArray::MathArray
constexpr MathArray(std::initializer_list< T > l)
Definition: math_array.h:49
bdm::MathArray::empty
const bool empty() const
Definition: math_array.h:70
bdm::MathArray::end
const T * end() const
Definition: math_array.h:97
bdm::operator<<
std::ostream & operator<<(std::ostream &o, const MathArray< T, N > &arr)
Definition: math_array.h:412
bdm::MathArray::operator--
MathArray operator--(int)
Definition: math_array.h:165
real_t.h
bdm::MathArray::MathArray
MathArray()
Default constructor.
Definition: math_array.h:39
bdm::MathArray::operator+
MathArray operator+(const T &rhs)
Definition: math_array.h:208
bdm::MathArray::IsZero
bool IsZero() const
Checks if vector is a zero vector, e.g. if all entries are zero.
Definition: math_array.h:341
bdm::MathArray::EntryWiseProduct
MathArray EntryWiseProduct(const MathArray &rhs) const
Definition: math_array.h:396
bdm::MathArray
Definition: math_array.h:36
bdm::MathArray::back
T & back()
Definition: math_array.h:109
bdm::MathArray::size
const size_t size() const
Definition: math_array.h:66
bdm::MathArray::operator-
MathArray operator-(const T &rhs)
Definition: math_array.h:254
bdm::MathArray::operator+=
MathArray & operator+=(const MathArray &rhs)
Definition: math_array.h:171
bdm::MathArray::operator-
const MathArray operator-(const MathArray &rhs) const
Definition: math_array.h:236
bdm::MathArray::at
T & at(size_t idx) noexcept(false)
Definition: math_array.h:88
bdm::MathArray::operator++
MathArray & operator++()
Definition: math_array.h:143
bdm::MathArray::Normalize
void Normalize()
Normalize the array in-place.
Definition: math_array.h:364
bdm::MathArray::GetNormalizedArray
MathArray GetNormalizedArray() const
Get a nomalized copy of the MathArray.
Definition: math_array.h:386
bdm::MathArray::operator*=
T & operator*=(const MathArray &rhs)=delete
bdm::MathArray::operator*
T operator*(const MathArray &rhs)
Definition: math_array.h:265