Program Listing for File array.hpp

↰ Return to documentation for file (source/pele/array.hpp)

#ifndef PYGMIN_ARRAY_H
#define PYGMIN_ARRAY_H

#include <algorithm>
#include <cassert>
#include <cmath>
#include <initializer_list>
#include <iostream>
#include <memory>
#include <numeric>
#include <stdexcept>
#include <vector>

namespace pele {

template <typename dtype>
class _ArrayMemory {
  std::vector<dtype> _vector;
  dtype
      *_data;
  size_t _size;

 public:
  _ArrayMemory() : _vector(), _data(_vector.data()), _size(_vector.size()) {}

  _ArrayMemory(size_t size)
      : _vector(size), _data(_vector.data()), _size(_vector.size()) {}

  _ArrayMemory(size_t size, dtype const &val)
      : _vector(size, val), _data(_vector.data()), _size(_vector.size()) {}

  _ArrayMemory(dtype *data, size_t size)
      : _vector(), _data(data), _size(size) {}

  inline size_t size() const { return _size; }

  inline dtype *data() { return _data; }
  inline dtype const *data() const { return _data; }
};

template <typename dtype>
class Array {
 protected:
  std::shared_ptr<_ArrayMemory<dtype>> _memory;
  dtype *_data;
  size_t _size;
 public:
  Array()
      : _memory(new _ArrayMemory<dtype>()),
        _data(_memory->data()),
        _size(_memory->size()) {}

  Array(size_t size)
      : _memory(new _ArrayMemory<dtype>(size)),
        _data(_memory->data()),
        _size(_memory->size()) {}

  Array(size_t size, dtype const &val)
      : _memory(new _ArrayMemory<dtype>(size, val)),
        _data(_memory->data()),
        _size(_memory->size())

  {}

  Array(std::initializer_list<dtype> data)
      : _memory(new _ArrayMemory<dtype>(data.size())),
        _data(_memory->data()),
        _size(_memory->size()) {
    std::copy(data.begin(), data.end(), _data);
  }

  Array(dtype *data, size_t size)
      : _memory(new _ArrayMemory<dtype>(data, size)),
        _data(_memory->data()),
        _size(_memory->size()) {}

  Array(dtype *data_begin, dtype *data_end)
      : _memory(new _ArrayMemory<dtype>(data_begin, data_end - data_begin)),
        _data(_memory->data()),
        _size(_memory->size()) {}

  Array(std::vector<dtype> &x)
      : _memory(new _ArrayMemory<dtype>(x.data(), x.size())),
        _data(_memory->data()),
        _size(_memory->size()) {}


  inline void wrap(Array<dtype> x) {
    _memory = x._memory;
    _data = x._data;
    _size = x._size;
  }

  inline dtype *data() { return _data; }
  inline dtype const *data() const { return _data; }

  inline size_t size() const { return _size; }

  inline dtype &operator[](const size_t i) { return data()[i]; }
  inline dtype const &operator[](const size_t i) const { return data()[i]; }

  typedef dtype *iterator;
  typedef dtype const *const_iterator;
  inline iterator begin() { return data(); }
  inline iterator end() { return data() + size(); }
  inline const_iterator begin() const { return data(); }
  inline const_iterator end() const { return data() + size(); }

  /*
   * Assignment operator: WRAP the data
   *
   * This is commented because it just duplicates the default assignment
   operator
   *
   Array<dtype> &operator=(const Array<dtype> & rhs){
   _memory = rhs._memory;
   _data = rhs._data;
   _size = rhs._size;
   }
  */

  inline bool operator==(Array<dtype> const rhs) const {
    return data() == rhs.data() and size() == rhs.size();
  }
  inline bool operator!=(Array<dtype> const rhs) const {
    return !operator==(rhs);
  }

  Array<dtype> &assign(const Array<dtype> &rhs) {
    // check for self assignment
    if ((*this) != rhs) {
      if (size() != rhs.size()) {
        throw std::runtime_error(
            "arrays must have the same size during assignment");
      }
      std::copy(rhs.begin(), rhs.end(), begin());
    }
    return *this;
  }

  Array<dtype> &assign(dtype const &d) {
    std::fill(begin(), end(), d);
    return *this;
  }

  // Inactive pointer-assignment prototype.
  //    Array<dtype> &assign(dtype const * const d) {
  //        std::copy(rhs.begin(), rhs.end(), begin());
  //        for (size_t i = 0; i < size(); ++i) {
  //
  //        }
  //        return *this;
  //    }

  Array<dtype> copy() const {
    Array<dtype> newarray(size());
    newarray.assign(*this);
    return newarray;
  }

  inline void free() {
    _memory = std::make_shared<_ArrayMemory<dtype>>();
    _data = _memory->data();
    _size = _memory->size();
  }

  inline bool empty() const { return size() == 0; }

  inline long int reference_count() const { return _memory.use_count(); }


  Array<dtype> operator+(const Array<dtype> &rhs) {
    if (size() != rhs.size()) {
      throw std::runtime_error("operator+: arrays must have the same size");
    }
    Array<dtype> newarray(size());
    const_iterator iter = rhs.begin();
    const_iterator iternew = (*this).begin();
    for (dtype &val : newarray) {
      val = *iternew + *iter;
      ++iter;
      ++iternew;
    }
    return newarray;
  }

  Array<dtype> operator-(const Array<dtype> &rhs) {
    if (size() != rhs.size()) {
      throw std::runtime_error("operator+: arrays must have the same size");
    }
    Array<dtype> newarray(size());
    const_iterator iter = rhs.begin();
    const_iterator iternew = (*this).begin();
    for (dtype &val : newarray) {
      val = *iternew - *iter;
      ++iter;
      ++iternew;
    }
    return newarray;
  }

  Array<dtype> operator-() {
    Array<dtype> newarray(size());
    const_iterator iter = (*this).begin();
    for (dtype &val : newarray) {
      val = -*iter;
      ++iter;
    }
    return newarray;
  }

  Array<dtype> &operator+=(const Array<dtype> &rhs) {
    if (size() != rhs.size()) {
      throw std::runtime_error("operator+=: arrays must have the same size");
    }
    const_iterator iter = rhs.begin();
    for (dtype &val : *this) {
      val += *iter;
      ++iter;
    }
    return *this;
  }

  Array<dtype> &operator+=(const dtype &rhs) {
    for (dtype &val : (*this)) {
      val += rhs;
    }
    return *this;
  }

  Array<dtype> &operator-=(const Array<dtype> &rhs) {
    if (size() != rhs.size()) {
      throw std::runtime_error("operator-=: arrays must have the same size");
    }
    const_iterator iter = rhs.begin();
    for (dtype &val : *this) {
      val -= *iter;
      ++iter;
    }
    return *this;
  }

  Array<dtype> &operator-=(const dtype &rhs) {
    for (dtype &val : (*this)) {
      val -= rhs;
    }
    return *this;
  }

  Array<dtype> &operator*=(const Array<dtype> &rhs) {
    if (size() != rhs.size()) {
      throw std::runtime_error("operator*=: arrays must have the same size");
    }
    const_iterator iter = rhs.begin();
    for (dtype &val : *this) {
      val *= *iter;
      ++iter;
    }
    return *this;
  }

  Array<dtype> &operator*=(const dtype &rhs) {
    for (dtype &val : (*this)) {
      val *= rhs;
    }
    return *this;
  }

  Array<dtype> operator*(const dtype rhs) {
    Array<dtype> result = this->copy();
    return (result *= rhs).copy();
  }

  Array<dtype> &operator/=(const Array<dtype> &rhs) {
    if (size() != rhs.size()) {
      throw std::runtime_error("operator/=: arrays must have the same size");
    }
    const_iterator iter = rhs.begin();
    for (dtype &val : *this) {
      val /= *iter;
      ++iter;
    }
    return *this;
  }

  Array<dtype> &operator/=(const dtype &rhs) {
    for (dtype &val : (*this)) {
      val /= rhs;
    }
    return *this;
  }

  /*SOME OTHER ARITHMETIC UTILITIES*/

  dtype sum() const {
    if (empty()) {
      throw std::runtime_error(
          "array::sum(): array is empty, can't sum array elements");
    }
    return std::accumulate(begin(), end(), dtype(0));
  }

  dtype prod() const {
    if (empty()) {
      throw std::runtime_error(
          "array::prod(): array is empty, can't take "
          "product of array elements");
    }
    return std::accumulate(begin(), end(), dtype(1), std::multiplies<dtype>());
  }

  std::vector<dtype> to_std_vector_copy() const {
    // copy the data into a new std::vector
    std::vector<dtype> std_vec(_size);
    std::copy(begin(), end(), std_vec.begin());
    return std_vec;
  }

  Array<dtype> view(size_t ibegin, size_t iend) const {
    if (iend <= ibegin) {
      throw std::invalid_argument("iend must larger than ibegin");
    }
    if (iend > size()) {
      throw std::invalid_argument("iend cannot be larger than array size");
    }
    Array<dtype> newarray(*this);
    newarray._data += ibegin;
    newarray._size = iend - ibegin;
    return newarray;
  }

  dtype get_max() const { return *std::max_element(begin(), end()); }
  dtype get_min() const { return *std::min_element(begin(), end()); }
};

// for array printing
template <class dtype>
inline std::ostream &operator<<(std::ostream &out, const Array<dtype> &a) {
  for (size_t i = 0; i < a.size(); ++i) {
    if (i > 0) out << ", ";
    out << a[i];
  }
  out << std::endl;
  return out;
}

inline double dot(Array<double> const &v1, Array<double> const &v2) {
  assert(v1.size() == v2.size());
  return std::inner_product(v1.begin(), v1.end(), v2.begin(), double(0));
}

inline double norm(Array<double> const &v) { return sqrt(dot(v, v)); }

template <class T, class U>
Array<T> operator*(const U rhs, const Array<T> &lhs) {
  Array<T> result = lhs.copy();
  return (result *= rhs).copy();
}

}  // namespace pele

#endif  // #ifndef PYGMIN_ARRAY_H