Added expression templates (#98)

- Works with NDArray
- Works with NDView
This commit is contained in:
Erik Fröjdh
2024-11-15 15:17:52 +01:00
committed by GitHub
parent 0d058274d5
commit e77b615293
12 changed files with 341 additions and 64 deletions

View File

@ -0,0 +1,99 @@
#pragma once
#include <cstdint> //int64_t
#include <cstddef> //size_t
#include <array>
#include <cassert>
namespace aare {
template <typename E, int64_t Ndim> class ArrayExpr {
public:
static constexpr bool is_leaf = false;
auto operator[](size_t i) const { return static_cast<E const &>(*this)[i]; }
auto operator()(size_t i) const { return static_cast<E const &>(*this)[i]; }
auto size() const { return static_cast<E const &>(*this).size(); }
std::array<int64_t, Ndim> shape() const { return static_cast<E const &>(*this).shape(); }
};
template <typename A, typename B, int64_t Ndim>
class ArrayAdd : public ArrayExpr<ArrayAdd<A, B, Ndim>, Ndim> {
const A &arr1_;
const B &arr2_;
public:
ArrayAdd(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
assert(arr1.size() == arr2.size());
}
auto operator[](int i) const { return arr1_[i] + arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
};
template <typename A, typename B, int64_t Ndim>
class ArraySub : public ArrayExpr<ArraySub<A, B, Ndim>, Ndim> {
const A &arr1_;
const B &arr2_;
public:
ArraySub(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
assert(arr1.size() == arr2.size());
}
auto operator[](int i) const { return arr1_[i] - arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
};
template <typename A, typename B, int64_t Ndim>
class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>,Ndim> {
const A &arr1_;
const B &arr2_;
public:
ArrayMul(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
assert(arr1.size() == arr2.size());
}
auto operator[](int i) const { return arr1_[i] * arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
};
template <typename A, typename B, int64_t Ndim>
class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
const A &arr1_;
const B &arr2_;
public:
ArrayDiv(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
assert(arr1.size() == arr2.size());
}
auto operator[](int i) const { return arr1_[i] / arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
};
template <typename A, typename B, int64_t Ndim>
auto operator+(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArrayAdd<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}
template <typename A, typename B, int64_t Ndim>
auto operator-(const ArrayExpr<A,Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArraySub<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}
template <typename A, typename B, int64_t Ndim>
auto operator*(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArrayMul<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}
template <typename A, typename B, int64_t Ndim>
auto operator/(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArrayDiv<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}
} // namespace aare

View File

@ -7,6 +7,7 @@ memory.
TODO! Add expression templates for operators
*/
#include "aare/ArrayExpr.hpp"
#include "aare/NDView.hpp"
#include <algorithm>
@ -20,7 +21,8 @@ TODO! Add expression templates for operators
namespace aare {
template <typename T, int64_t Ndim = 2> class NDArray {
template <typename T, int64_t Ndim = 2>
class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::array<int64_t, Ndim> shape_;
std::array<int64_t, Ndim> strides_;
size_t size_{};
@ -43,7 +45,7 @@ template <typename T, int64_t Ndim = 2> class NDArray {
: shape_(shape), strides_(c_strides<Ndim>(shape_)),
size_(std::accumulate(shape_.begin(), shape_.end(), 1,
std::multiplies<>())),
data_(new T[size_]) {};
data_(new T[size_]) {}
/**
* @brief Construct a new NDArray object with a shape and value.
@ -69,8 +71,7 @@ template <typename T, int64_t Ndim = 2> class NDArray {
NDArray(NDArray &&other) noexcept
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
size_(other.size_), data_(other.data_) {
other.reset();
other.reset(); // TODO! is this necessary?
}
// Copy constructor
@ -78,7 +79,14 @@ template <typename T, int64_t Ndim = 2> class NDArray {
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
size_(other.size_), data_(new T[size_]) {
std::copy(other.data_, other.data_ + size_, data_);
// fmt::print("NDArray(const NDArray &other)\n");
}
// Conversion operator from array expression to array
template <typename E>
NDArray(ArrayExpr<E, Ndim> &&expr) : NDArray(expr.shape()) {
for (int i = 0; i < size_; ++i) {
data_[i] = expr[i];
}
}
~NDArray() { delete[] data_; }
@ -90,14 +98,10 @@ template <typename T, int64_t Ndim = 2> class NDArray {
NDArray &operator=(NDArray &&other) noexcept; // Move assign
NDArray &operator=(const NDArray &other); // Copy assign
NDArray operator+(const NDArray &other);
NDArray &operator+=(const NDArray &other);
NDArray operator-(const NDArray &other);
NDArray &operator-=(const NDArray &other);
NDArray operator*(const NDArray &other);
NDArray &operator*=(const NDArray &other);
NDArray operator/(const NDArray &other);
// NDArray& operator/=(const NDArray& other);
template <typename V> NDArray &operator/=(const NDArray<V, Ndim> &other) {
@ -151,12 +155,16 @@ template <typename T, int64_t Ndim = 2> class NDArray {
return data_[element_offset(strides_, index...)];
}
// TODO! is int the right type for index?
T &operator()(int i) { return data_[i]; }
const T &operator()(int i) const { return data_[i]; }
T &operator[](int i) { return data_[i]; }
const T &operator[](int i) const { return data_[i]; }
T *data() { return data_; }
std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
uint64_t size() const { return size_; }
size_t size() const { return size_; }
size_t total_bytes() const { return size_ * sizeof(T); }
std::array<int64_t, Ndim> shape() const noexcept { return shape_; }
int64_t shape(int64_t i) const noexcept { return shape_[i]; }
@ -204,17 +212,11 @@ NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
return *this;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const NDArray &other) {
NDArray result(*this);
result += other;
return result;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
for (uint32_t i = 0; i < size_; ++i) {
for (size_t i = 0; i < size_; ++i) {
data_[i] += other.data_[i];
}
return *this;
@ -222,13 +224,6 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator-(const NDArray &other) {
NDArray result{*this};
result -= other;
return result;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
// check shape
@ -240,12 +235,6 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const NDArray &other) {
NDArray result = *this;
result *= other;
return result;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
@ -259,13 +248,6 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator/(const NDArray &other) {
NDArray result = *this;
result /= other;
return result;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
for (auto it = begin(); it != end(); ++it)
@ -276,7 +258,7 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
template <typename T, int64_t Ndim>
NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
if (shape_ == other.shape_) {
NDArray<bool> result{shape_};
NDArray<bool, Ndim> result{shape_};
for (int i = 0; i < size_; ++i) {
result(i) = (data_[i] > other.data_[i]);
}

View File

@ -1,4 +1,7 @@
#pragma once
#include "aare/ArrayExpr.hpp"
#include <algorithm>
#include <array>
#include <cassert>
@ -45,7 +48,7 @@ template <int64_t Ndim> std::array<int64_t, Ndim> make_array(const std::vector<i
return arr;
}
template <typename T, int64_t Ndim = 2> class NDView {
template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
public:
NDView() = default;
~NDView() = default;
@ -68,7 +71,7 @@ template <typename T, int64_t Ndim = 2> class NDView {
return buffer_[element_offset(strides_, index...)];
}
uint64_t size() const { return size_; }
size_t size() const { return size_; }
size_t total_bytes() const { return size_ * sizeof(T); }
std::array<int64_t, Ndim> strides() const noexcept { return strides_; }