math_array.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// Copyright (C) 2021 CERN & University of Surrey for the benefit of the
// BioDynaMo collaboration. All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
//
// See the LICENSE file distributed with this work for details.
// See the NOTICE file distributed with this work for additional information
// regarding copyright ownership.
//
// -----------------------------------------------------------------------------
 
#ifndef CORE_CONTAINER_MATH_ARRAY_H_
#define CORE_CONTAINER_MATH_ARRAY_H_
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <numeric>
#include <ostream>
#include <stdexcept>
#include <utility>
 
#include "core/real_t.h"
#include "core/util/log.h"
#include "core/util/root.h"
 
namespace bdm {
 
template <class T, std::size_t N>
class MathArray {  // NOLINT
 public:
  MathArray() {
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      data_[i] = T();
    }
  }
 
  constexpr MathArray(std::initializer_list<T> l) {
    assert(l.size() <= N);
    auto it = l.begin();
    for (uint64_t i = 0; i < N; i++) {
      data_[i] = *(it++);
    }
    for (uint64_t i = l.size(); i < N; i++) {
      data_[i] = T();
    }
  }
 
  inline const T* data() const { return &data_[0]; }  // NOLINT
 
  inline const size_t size() const { return N; }  // NOLINT
 
  inline const bool empty() const { return N == 0; }  // NOLINT
 
  T& operator[](size_t idx) { return data_[idx]; }
 
  const T& operator[](size_t idx) const { return data_[idx]; }
 
  T& at(size_t idx) noexcept(false) {  // NOLINT
    if (idx > size() || idx < 0) {
      throw std::out_of_range("The index is out of range");
    }
    return data_[idx];
  }
 
  const T* begin() const { return &(data_[0]); }  // NOLINT
 
  const T* end() const { return &(data_[N]); }  // NOLINT
 
  T* begin() { return &(data_[0]); }  // NOLINT
 
  T* end() { return &(data_[N]); }  // NOLINT
 
  T& front() { return *(this->begin()); }  // NOLINT
 
  T& back() {  // NOLINT
    auto tmp = this->end();
    tmp--;
    return *tmp;
  }
 
  MathArray& operator=(const MathArray& other) {
    if (this != &other) {
      assert(other.size() == N);
      std::copy(other.data_, other.data_ + other.size(), data_);
    }
    return *this;
  }
 
  bool operator==(const MathArray& other) const {
    if (other.size() != N) {
      return false;
    }
    for (size_t i = 0; i < N; i++) {
      if (other[i] != data_[i]) {
        return false;
      }
    }
    return true;
  }
 
  bool operator!=(const MathArray& other) const { return !operator==(other); }
 
  MathArray& operator++() {
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      ++data_[i];
    }
    return *this;
  }
 
  MathArray operator++(int) {
    MathArray tmp(*this);
    operator++();
    return tmp;
  }
 
  MathArray& operator--() {
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      --data_[i];
    }
    return *this;
  }
 
  MathArray operator--(int) {
    MathArray tmp(*this);
    operator--();
    return tmp;
  }
 
  MathArray& operator+=(const MathArray& rhs) {
    assert(N == rhs.size());
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      data_[i] += rhs[i];
    }
    return *this;
  }
 
  MathArray operator+(const MathArray& rhs) {
    assert(size() == rhs.size());
    MathArray tmp;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      tmp[i] = data_[i] + rhs[i];
    }
    return tmp;
  }
 
  const MathArray operator+(const MathArray& rhs) const {
    assert(size() == rhs.size());
    MathArray tmp;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      tmp[i] = data_[i] + rhs[i];
    }
    return tmp;
  }
 
  MathArray& operator+=(const T& rhs) {
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      data_[i] += rhs;
    }
    return *this;
  }
 
  MathArray operator+(const T& rhs) {
    MathArray tmp;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      tmp[i] = data_[i] + rhs;
    }
    return tmp;
  }
 
  MathArray& operator-=(const MathArray& rhs) {
    assert(size() == rhs.size());
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      data_[i] -= rhs[i];
    }
    return *this;
  }
 
  MathArray operator-(const MathArray& rhs) {
    assert(size() == rhs.size());
    MathArray tmp;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      tmp[i] = data_[i] - rhs[i];
    }
    return tmp;
  }
 
  const MathArray operator-(const MathArray& rhs) const {
    assert(size() == rhs.size());
    MathArray tmp;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      tmp[i] = data_[i] - rhs[i];
    }
    return tmp;
  }
 
  MathArray& operator-=(const T& rhs) {
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      data_[i] -= rhs;
    }
    return *this;
  }
 
  MathArray operator-(const T& rhs) {
    MathArray tmp;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      tmp[i] = data_[i] - rhs;
    }
    return tmp;
  }
 
  T& operator*=(const MathArray& rhs) = delete;
 
  T operator*(const MathArray& rhs) {
    assert(size() == rhs.size());
    T result = 0;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      result += data_[i] * rhs[i];
    }
    return result;
  }
 
  const T operator*(const MathArray& rhs) const {
    assert(size() == rhs.size());
    T result = 0;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      result += data_[i] * rhs[i];
    }
    return result;
  }
 
  MathArray& operator*=(const T& k) {
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      data_[i] *= k;
    }
    return *this;
  }
 
  MathArray operator*(const T& k) {
    MathArray tmp;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      tmp[i] = data_[i] * k;
    }
    return tmp;
  }
 
  const MathArray operator*(const T& k) const {
    MathArray tmp;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      tmp[i] = data_[i] * k;
    }
    return tmp;
  }
 
  MathArray& operator/=(const T& k) {
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      data_[i] /= k;
    }
    return *this;
  }
 
  MathArray operator/(const T& k) {
    MathArray tmp;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      tmp[i] = data_[i] / k;
    }
    return tmp;
  }
 
  MathArray& fill(const T& k) {  // NOLINT
    std::fill(std::begin(data_), std::end(data_), k);
    return *this;
  }
 
  T Sum() const { return std::accumulate(begin(), end(), 0); }
 
  bool IsZero() const {
    for (size_t i = 0; i < N; i++) {
      if (data_[i] != 0) {
        return false;
      }
    }
    return true;
  }
 
  T Norm() const {
    T result = 0;
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      result += data_[i] * data_[i];
    }
    result = std::sqrt(result);
 
    return result;
  }
 
  void Normalize() {
    T norm = Norm();
    Normalize(norm);
  }
 
  void Normalize(T norm) {
    if (norm == 0) {
      Log::Fatal("MathArray::Normalize",
                 "You tried to normalize a zero vector. "
                 "This cannot be done. Exiting.");
    }
#pragma omp simd
    for (size_t i = 0; i < N; i++) {
      data_[i] /= norm;
    }
  }
 
  MathArray GetNormalizedArray() const {
    MathArray normalized_array(*this);
    normalized_array.Normalize();
    return normalized_array;
  }
 
  MathArray EntryWiseProduct(const MathArray& rhs) const {
    assert(rhs.size() == N);
    MathArray tmp;
#pragma omp simd
    for (size_t i = 0; i < N; ++i) {
      tmp[i] = data_[i] * rhs[i];
    }
    return tmp;
  }
 
 private:
  T data_[N];
  BDM_CLASS_DEF_NV(MathArray, 1);  // NOLINT
};
 
template <class T, std::size_t N>
std::ostream& operator<<(std::ostream& o, const MathArray<T, N>& arr) {
  for (size_t i = 0; i < N; i++) {
    o << arr[i];
    if (i != N - 1) {
      o << ", ";
    }
  }
  return o;
}
 
// Note: 1) We pass by value to allow for copy-elision and move semantics
//          optimization.
//       2) We do return MathArray<T, N> because a const MathArray<T, N>
//          prevents move semantics in C++11.
//       see https://tinyurl.com/left-multiply
template <class T, std::size_t N>
MathArray<T, N> operator*(T const& scalar, MathArray<T, N> array) {
  return array *= scalar;
}
 
using Real3 = MathArray<real_t, 3>;
using Float3 = MathArray<float, 3>;
using Double3 = MathArray<double, 3>;
 
using Real4 = MathArray<real_t, 4>;
using Float4 = MathArray<float, 4>;
using Double4 = MathArray<double, 4>;
 
}  // namespace bdm
 
#endif  // CORE_CONTAINER_MATH_ARRAY_H_