data.hpp

Go to the documentation of this file.

/*
 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
 * Copyright (c) 2022 Lodz University of Technology
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 */
#ifndef PLASK__DATA_H
#define PLASK__DATA_H
 
#include <iterator>
#include <algorithm>
#include <iostream>
#include <initializer_list>
#include <atomic>
#include <type_traits>
#include <memory>       // std::unique_ptr
#include <cassert>
#include <type_traits>  // std::is_trivially_destructible, std::false_type, std::true_type
 
#include "memalloc.hpp"
#include "math.hpp"
#include "exceptions.hpp"
#include "log/log.hpp"
 
#include <boost/type_traits.hpp>
 
namespace plask {
 
namespace detail {
 
    template <class T>
    inline void _construct_array(T* first, T* last, const std::false_type&) {
        while(first != last) {
            new(first) T();
            ++first;
        }
    }
 
    template <class T>
    inline void _construct_array(T*, T*, const std::true_type&) {
    }
 
    template <class T>
    inline void _destroy_array(T* first, T* last, const std::false_type&) {
        while(last != first) {
            --last;
            last->~T();
        }
    }
 
    template <class T>
    inline void _destroy_array(T*, T*, const std::true_type&) {
    }
 
    template <class T>
    inline void construct_array(T* first, T* last) {
        _construct_array(first, last, std::is_trivially_default_constructible<T>());
    }
 
    template <class T>
    inline void destroy_array(T* first, T* last) {
        _destroy_array(first, last, std::is_trivially_destructible<T>());
    }
 
    // We consider dcomplex as POD (do not value-initialize memory)
 
    template <>
    inline void construct_array(dcomplex* first, dcomplex* last) {
        _construct_array(first, last, std::true_type());
    }
 
    template <>
    inline void destroy_array(dcomplex* first, dcomplex* last) {
        _destroy_array(first, last, std::true_type());
    }
 
 
    struct DataVectorGC {
        // Count is atomic so many threads can increment and decrement it at same time.
        // If it is 0, it means that there has been only one DataVector object, so probably one thread uses it.
        std::atomic<unsigned> count;
 
        std::function<void(void*)>* deleter;
 
        explicit DataVectorGC(unsigned initial): count(initial), deleter(nullptr) {}
 
        explicit DataVectorGC(unsigned initial, const std::function<void(void*)>& deleter):
            count(initial), deleter(new std::function<void(void*)>(deleter)) {}
 
        explicit DataVectorGC(unsigned initial, std::function<void(void*)>&& deleter):
            count(initial), deleter(new std::function<void(void*)>(std::forward<std::function<void(void*)>>(deleter))) {}
 
        void free(void* data) {
            if (deleter) (*deleter)(data);
            else aligned_free(data);
        }
 
        ~DataVectorGC() { delete deleter; }
    };
 
}   // namespace detail
 
template <typename T>
struct DataVector {
 
    typedef typename std::remove_const<T>::type VT;
    typedef const T CT;
 
    typedef detail::DataVectorGC Gc;
  private:
 
    std::size_t size_;                  
    Gc* gc_;                            
    T* data_;                           
 
    void dec_ref() {    // see http://www.boost.org/doc/libs/1_53_0/doc/html/atomic/usage_examples.html "Reference counting" for optimal memory access description
        if (gc_ && gc_->count.fetch_sub(1, std::memory_order_release) == 1) {
            std::atomic_thread_fence(std::memory_order_acquire);
            detail::destroy_array(data_, data_ + size_);
            gc_->free(reinterpret_cast<void*>(const_cast<VT*>(data_)));
            delete gc_;
        }
    }
 
    void inc_ref() {
        if (gc_) gc_->count.fetch_add(1, std::memory_order_relaxed);
    }
 
    friend struct DataVector<VT>;
    friend struct DataVector<CT>;
 
  public:
 
    typedef T value_type;               
 
    typedef T* iterator;                
    typedef const T* const_iterator;    
 
    DataVector(): size_(0), gc_(nullptr), data_(nullptr) {}
 
    explicit DataVector(std::size_t size): size_(size), gc_(new Gc(1)), data_(aligned_malloc<T>(size)) {
        detail::construct_array(data_, data_ + size);
    }
 
    DataVector(std::size_t size, const T& value): size_(size) {
        std::unique_ptr<typename std::remove_const<T>::type[], aligned_deleter<T>>
            data_non_const(aligned_malloc<VT>(size));
        std::fill_n(data_non_const.get(), size, value);   // this may throw, but no memory leak than
        gc_ = new Gc(1);
        data_ = data_non_const.release();
    }
 
    DataVector(const DataVector<CT>& src): size_(src.size_), gc_(src.gc_), data_(const_cast<T*>(src.data_)) { inc_ref(); }
 
    DataVector(const DataVector<VT>& src): size_(src.size_), gc_(src.gc_), data_(const_cast<T*>(src.data_)) { inc_ref(); }
 
    template <typename TS>
    DataVector(const DataVector<TS>& src): size_(src.size()), gc_(new Gc(1)), data_(aligned_malloc<T>(src.size())) {
        write_debug("copying data in constructor");
        std::copy(src.begin(), src.end(), const_cast<VT*>(data_));
    }
 
    DataVector<T>& operator=(const DataVector<CT>& M) {
        const_cast<DataVector<CT>&>(M).inc_ref();   // must be called before dec_ref in case of self-asigment with 1 reference
        this->dec_ref();                            // release old content, this can delete the old data
        size_ = M.size_;
        data_ = const_cast<T*>(M.data_);
        gc_ = M.gc_;
        return *this;
    }
 
    DataVector<T>& operator=(const DataVector<VT>& M) {
        const_cast<DataVector<VT>&>(M).inc_ref();   // must be called before dec_ref in case of self-asigment with 1 reference
        this->dec_ref();                            // release old content, this can delete the old data
        size_ = M.size_;
        data_ = const_cast<T*>(M.data_);
        gc_ = M.gc_;
        return *this;
    }
 
    DataVector(DataVector<CT>&& src) noexcept: size_(src.size_), gc_(src.gc_), data_(const_cast<T*>(src.data_)) {
        src.gc_ = nullptr;
    }
 
    DataVector(DataVector<VT>&& src) noexcept: size_(src.size_), gc_(src.gc_), data_(const_cast<T*>(src.data_)) {
        src.gc_ = nullptr;
    }
 
    DataVector<T>& operator=(DataVector<T>&& src) noexcept {
        swap(src);
        return *this;
    }
 
    template <typename TS>
    DataVector(TS* existing_data, std::size_t size):
        size_(size), gc_(nullptr), data_(existing_data) {}
 
    template <typename TS>
    DataVector(TS* existing_data, std::size_t size, const std::function<void(void*)>& deleter):
        size_(size), gc_(new Gc(1, deleter)), data_(existing_data) {
    }
 
    template <typename TS>
    DataVector(TS* existing_data, std::size_t size, std::function<void(void*)>&& deleter):
        size_(size), gc_(new Gc(1, std::forward<std::function<void(void*)>>(deleter))), data_(existing_data) {
    }
 
    DataVector(std::initializer_list<T> init): size_(init.size()), gc_(new Gc(1)), data_(aligned_malloc<T>(size_)) {
        std::copy(init.begin(), init.end(), const_cast<VT*>(data_));
    }
 
    template <typename TS>
    DataVector(std::initializer_list<TS> init): size_(init.size()), gc_(new Gc(1)), data_(aligned_malloc<T>(size_)) {
        std::copy(init.begin(), init.end(), const_cast<VT*>(data_));
    }
 
    ~DataVector() { dec_ref(); }
 
    void reset() {
        dec_ref();
        size_ = 0;
        gc_ = nullptr;
        data_ = nullptr;
    }
 
    template <typename TS>
    void reset(TS* existing_data, std::size_t size) {
        dec_ref();
        gc_ = nullptr;
        size_ = size;
        data_ = existing_data;
    }
 
    template <typename TS>
    void reset(TS* existing_data, std::size_t size, const std::function<void(void*)>& deleter) {
        dec_ref();
        gc_ = new Gc(1, deleter);
        size_ = size;
        data_ = existing_data;
    }
 
    template <typename TS>
    void reset(TS* existing_data, std::size_t size, std::function<void(void*)>&& deleter) {
        dec_ref();
        gc_ = new Gc(1, std::forward<std::function<void(void*)>>(deleter));
        size_ = size;
        data_ = existing_data;
    }
 
    void reset(std::size_t size) {
        //TODO to consider: (when clang will be fixed, now it has no std::is_default_constructible but only non-standard std::has_trivial_default_constructor)
        //if (std::is_default_constructible<T>::value &&   //this is known at compile time and I believe that compiler optimize-out whole if when it is false
        //    size == size_ && gc_ && gc_->count == 1 && ! gc_->deleter) return;
        dec_ref();
        data_ = aligned_malloc<T>(size);
        detail::construct_array(data_, data_ + size);
        gc_ = new Gc(1);
        size_ = size;
    }
 
    void reset(std::size_t size, const T& value) {
        std::unique_ptr<VT[], aligned_deleter<T>>
            data_non_const(aligned_malloc<VT>(size));
        std::fill_n(data_non_const.get(), size, value);   // this may throw, than our data will not be changed
        dec_ref();
        gc_ = new Gc(1);    //this also may throw
        data_ = data_non_const.release();
        size_ = size;
    }
 
    template <typename InIterT>
    void reset(InIterT begin, InIterT end) {
        std::unique_ptr<VT[], aligned_deleter<T>>
            data_non_const(aligned_malloc<VT>(std::distance(begin,end)));
        std::copy(begin, end, data_non_const.get());    // this may throw, and than our vector will not be changed
        dec_ref();
        gc_ = new Gc(1);
        size_ = std::distance(begin, end);
        data_ = data_non_const.release();
    }
 
    const_iterator begin() const { return data_; }
 
    iterator begin() { return data_; }
 
    const_iterator end() const { return data_ + size_; }
 
    iterator end() { return data_ + size_; }
 
    std::size_t size() const { return size_; }
 
    const T* data() const { return data_; }
 
    T* data() { return data_; }
 
    const T& operator[](std::size_t n) const { assert(n < size_); return data_[n]; }
 
    T& operator[](std::size_t n) { assert(n < size_); return data_[n]; }
 
    operator bool() const { return data_ != nullptr; }
 
    DataVector<VT> copy() const {
        DataVector<VT> new_data(size_);
        std::copy(begin(), end(), new_data.begin());
        return new_data;
    }
 
    bool unique() const {
        return (gc_ != nullptr) && (gc_->count == 1);
    }
 
    DataVector<VT> remove_const() const {
        DataVector<VT> result(const_cast<VT*>(this->data()), this->size());
        result.gc_ = this->gc_;
        result.inc_ref();
        return result;
    }
 
    DataVector<VT> claim() const {
        return (unique() && !gc_->deleter)? remove_const() : copy();
    }
 
    void swap(DataVector<T>& other) noexcept {
        std::swap(size_, other.size_);
        std::swap(gc_, other.gc_);
        std::swap(data_, other.data_);
    }
 
    template <typename O>
    void fill(const O& value) {
        std::fill_n(data_, size_, value);
    }
 
    inline DataVector<T> getSubarrayRef(std::size_t begin_index, std::size_t subarray_size) {
        assert(begin_index + subarray_size <= size_);
        return DataVector<T>(data_ + begin_index, subarray_size);
    }
 
    inline DataVector<const T> getSubarrayRef(std::size_t begin_index, std::size_t subarray_size) const {
        assert(begin_index + subarray_size <= size_);
        return DataVector<const T>(data_ + begin_index, subarray_size);
    }
 
    inline DataVector<T> getSubarrayCopy(std::size_t begin_index, std::size_t subarray_size) const {
        return getSubarrayRef(begin_index, subarray_size).copy();
    }
 
    inline DataVector<T> getSubrangeRef(iterator begin, iterator end) {
        assert(this->begin() <= begin && begin <= end && end <= this->end());
        return DataVector<T>(begin, end - begin);
    }
 
    inline DataVector<const T> getSubrangeRef(const_iterator begin, const_iterator end) const {
        assert(this->begin() <= begin && begin <= end && end <= this->end());
        return DataVector<const T>(begin, end-begin);
    }
 
    inline DataVector<T> getSubrangeCopy(const_iterator begin, const_iterator end) {
        assert(this->begin() <= begin && begin <= end && end <= this->end());
        return getSubrangeRef(begin, end).copy();
    }
};
 
template<class T1, class T2> inline
bool operator== ( DataVector<T1> const& a, DataVector<T2> const& b)
{ return a.size() == b.size() && std::equal(a.begin(), a.end(), b.begin()); }
 
template<class T1, class T2> inline
bool operator!= ( DataVector<T1> const& a, DataVector<T2> const& b) { return !(a==b); }
 
template<class T1, class T2> inline
bool operator< ( DataVector<T1> const& a, DataVector<T2> const& b)
{ return std::lexicographical_compare(a.begin(), a.end(), b.begin(), b.end()); }
 
template<class T1, class T2> inline
bool operator> ( DataVector<T1> const& a, DataVector<T2> const& b) { return b < a; }
 
template<class T1, class T2> inline
bool operator<= ( DataVector<T1> const& a, DataVector<T2> const& b) { return !(b < a); }
 
template<class T1, class T2> inline
bool operator>= ( DataVector<T1> const& a, DataVector<T2> const& b) { return !(a < b); }
 
template<class T>
std::ostream& operator<<(std::ostream& out, DataVector<T> const& to_print) {
    out << '[';
    return print_seq(out, to_print.begin(), to_print.end()) << ']';
}
 
template<class T, class S>
DataVector<T>& operator+=(DataVector<T>& to_inc, DataVector<S> const& inc_val) {
    if (to_inc.size() != inc_val.size())
        throw DataError("data vectors sizes differ ([{0}] += [%2])", to_inc.size(), inc_val.size());
    for (std::size_t i = 0; i < to_inc.size(); ++i)
        to_inc[i] += inc_val[i];
    return to_inc;
}
 
template<class T, class S>
DataVector<T>& operator-=(DataVector<T>& to_dec, DataVector<S> const& dec_val) {
    if (to_dec.size() != dec_val.size())
        throw DataError("data vectors sizes differ ([{0}] -= [%2])", to_dec.size(), dec_val.size());
    for (std::size_t i = 0; i < to_dec.size(); ++i)
        to_dec[i] -= dec_val[i];
    return to_dec;
}
 
template<typename T, typename S>
DataVector<T>& operator*=(DataVector<T>& vec, S a) {
    for (std::size_t i = 0; i < vec.size(); ++i)
        vec[i] *= a;
    return vec;
}
 
template<typename T, typename S>
DataVector<T>& operator/=(DataVector<T>& vec, S a) {
    auto ia = 1. / a;
    for (std::size_t i = 0; i < vec.size(); ++i)
        vec[i] *= ia;
    return vec;
}
 
template<typename T1, typename T2>
DataVector<typename std::remove_cv<decltype(T1()+T2())>::type> operator+(DataVector<T1> const& vec1, DataVector<T2> const& vec2) {
    if (vec1.size() != vec2.size())
        throw DataError("data vectors sizes differ ([{0}] + [%2])", vec1.size(), vec2.size());
    DataVector<typename std::remove_cv<decltype(T1()+T2())>::type> result(vec1.size());
    for (std::size_t i = 0; i < vec1.size(); ++i)
        result[i] = vec1[i] + vec2[i];
    return result;
}
 
template<typename T1, typename T2>
DataVector<typename std::remove_cv<decltype(T1()-T2())>::type> operator-(DataVector<T1> const& vec1, DataVector<T2> const& vec2) {
    if (vec1.size() != vec2.size())
        throw DataError("data vectors sizes differ ([{0}] - [%2])", vec1.size(), vec2.size());
    DataVector<typename std::remove_cv<decltype(T1()-T2())>::type> result(vec1.size());
    for (std::size_t i = 0; i < vec1.size(); ++i)
        result [i] = vec1[i] - vec2[i];
    return result;
}
 
template<class T>
DataVector<typename std::remove_cv<T>::type> operator-(DataVector<T> const& vec) {
    DataVector<typename std::remove_cv<T>::type> result(vec.size());
    for (std::size_t i = 0; i < vec.size(); ++i)
        result [i] = -vec[i];
    return result;
}
 
template<typename T, typename S>
DataVector<typename std::remove_cv<decltype(T()*S())>::type> operator*(DataVector<T> const& vec, S a) {
    DataVector<typename std::remove_cv<decltype(T()*S())>::type> result(vec.size());
    for (std::size_t i = 0; i < vec.size(); ++i)
        result[i] = vec[i] * a;
    return result;
}
 
template<typename T, typename S>
DataVector<typename std::remove_cv<decltype(S()*T())>::type> operator*(S a, DataVector<T> const& vec) {
    DataVector<typename std::remove_cv<decltype(S()*T())>::type> result(vec.size());
    for (std::size_t i = 0; i < vec.size(); ++i)
        result[i] = a * vec[i];
    return result;
}
 
template<typename T, typename S>
DataVector<typename std::remove_cv<decltype(T()*(1./S()))>::type> operator/(DataVector<T> const& vec, S a) {
    auto ia = 1. / a;
    DataVector<typename std::remove_cv<decltype(T()*(1./S()))>::type> result(vec.size());
    for (std::size_t i = 0; i < vec.size(); ++i)
        result[i] = vec[i] * ia;
    return result;
}
 
 
 
template <class T>
typename std::remove_cv<T>::type accumulate(const DataVector<T>& to_accum, typename std::remove_cv<T>::type initial=typename std::remove_cv<T>::type()) {
    for (std::size_t i = 0; i < to_accum.size(); ++i)
        initial += to_accum[i];
    return initial;
}
 
template <class T>
T average(const DataVector<T>& v) {
    return accumulate(v) / v.size();
}
 
template<typename RT, typename T>
inline DataVector<RT> const_data_cast(const DataVector<T>& src) {
    return src.remove_const();
}
 
/*
PLASK_API_EXTERN_TEMPLATE_STRUCT(DataVector<double>)
PLASK_API_EXTERN_TEMPLATE_STRUCT(DataVector<const double>)
PLASK_API_EXTERN_TEMPLATE_STRUCT(PLASK_API DataVector<std::complex<double>>)
PLASK_API_EXTERN_TEMPLATE_STRUCT(DataVector<const std::complex<double>>)
*/
 
}   // namespace plask
 
namespace std {
    template <typename T>
    void swap(plask::DataVector<T>& s1, plask::DataVector<T>& s2) noexcept {
      s1.swap(s2);
    }
}   // namespace std
 
 
#endif // PLASK__DATA_H