MultiArray.h

Go to the documentation of this file.

// Copyright András Vukics 2022–2026. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE.txt)
#pragma once
 
#include "Traits.h"
 
#include <boost/serialization/array.hpp>
#include <boost/serialization/vector.hpp>
 
#include <array>
#include <concepts>
#include <span>
#include <stdexcept>
#include <ranges>
#include <vector>
 
 
namespace cppqedutils {
 
template <typename Type, std::size_t... sizes>
constexpr auto concatenate(const std::array<Type, sizes>&... arrays)
{
  std::array<Type, (sizes + ...)> result;
  std::size_t index{};
 
  ((std::copy_n(arrays.begin(), sizes, result.begin() + index), index += sizes), ...);
 
  return result;
}
 
 
template <size_t RANK>
using Extents=std::array<size_t,RANK>;
 
 
template <size_t RANK>
auto& incrementMultiIndex(Extents<RANK>& idx, Extents<RANK> extents)
{
  const auto increment=[&](auto n, auto inc)
  {
    using namespace hana::literals;
    if constexpr (n!=0_c) {
      if (idx[n]==extents[n]-1) {idx[n]=0; inc(n-1_c,inc);}
      else idx[n]++;
    }
    else idx[0]++; // This will eventually put the index into an illegal state, but this is how all (unchecked) iterators work.
  };
  
  increment(hana::llong_c<RANK-1>,increment);
  
  return idx;
}
 
inline auto& incrementMultiIndex(Extents<0>& idx, Extents<0>) {return idx;}
 
 
namespace multiarray {
 
inline constexpr struct Placeholder {} _{};
 
} // multiarray
 
 
 
template <typename T, size_t RANK>
class MultiArrayView
{
public:
  MultiArrayView(const MultiArrayView&) = default; MultiArrayView(MultiArrayView&&) = default; MultiArrayView() = default;
  MultiArrayView& operator=(const MultiArrayView&) = default; MultiArrayView& operator=(MultiArrayView&&) = default;
  
  MultiArrayView(Extents<RANK> extents, Extents<RANK> strides, size_t offset, auto&&... dataView)
    : extents{extents}, strides{strides}, offset{offset}, dataView{FWD(dataView)...} {}
 
  operator MultiArrayView<const T, RANK>() const requires ( !std::is_const_v<T> ) {return {extents,strides,offset,dataView}; }
 
  void checkBounds(std::convertible_to<size_t> auto ... i) const
  {
#ifndef   NDEBUG
    auto e=extents.begin();
    (... , [&] {
      if (auto eComp=e++; static_cast<size_t>(i) >= *eComp)
        throw std::range_error("Index position: "+std::to_string(eComp-extents.begin())+", index value: "+std::to_string(i)+", extent: "+std::to_string(*eComp));
    } () );
#else
    (... , (void)i); // silence unused warning
#endif // NDEBUG
  }
 
 
  T& operator() (Extents<RANK> idx) const
  {
    std::apply([this] (auto ... args) { this->checkBounds(args ...); },idx);
 
    return dataView[std::ranges::fold_left( std::views::zip(idx,strides), offset,
                                            [&] (auto init, auto ids) {return init+get<0>(ids)*get<1>(ids);} ) ];
  }
  
  T& operator() (std::convertible_to<size_t> auto ... i) const requires (sizeof...(i)==RANK)
  {
    checkBounds(i...);
    size_t idx=0;
    auto s=strides.begin();
    return dataView[ offset + ( ... , [&] { return idx+=(*s++)*i;} () ) ];
  }
  
  T& operator() (std::convertible_to<size_t> auto i) const requires (RANK==1) {checkBounds(i); return dataView[offset+strides[0]*i];}
 
  T& operator() (std::convertible_to<size_t> auto i, std::convertible_to<size_t> auto j) const requires (RANK==2) {
    checkBounds(i,j); return dataView[offset+strides[0]*i+strides[1]*j];
  }
 
  auto operator() (std::convertible_to<size_t> auto i, multiarray::Placeholder) const requires (RANK==2) {
    return MultiArrayView<T,1>{ {extents[1]}, {strides[1]}, offset+i*strides[0], dataView};
  }
 
  auto operator() (multiarray::Placeholder, std::convertible_to<size_t> auto j) const requires (RANK==2) {
    return MultiArrayView<T,1>{ {extents[0]}, {strides[0]}, offset+j*strides[1], dataView};
  }
 
  Extents<RANK> extents, strides;
 
  size_t offset;
  
  std::span<T> dataView;
  
};
 
 
template <typename T, size_t RANK>
constexpr auto passByValue_v<MultiArrayView<T,RANK>> = true;
 
 
template <typename T, size_t RANK> requires ( !std::is_const_v<T> )
using MultiArrayConstView = MultiArrayView<const T, RANK>;
 
 
// MultiArrayView is itself a reference — passing it by value is equivalent to passing by reference.
template <typename T, size_t RANK> requires ( !std::is_const_v<T> )
struct ReferenceMF<MultiArrayView<T,RANK>> : std::type_identity<MultiArrayView<T,RANK>> {};
 
template <typename T, size_t RANK> requires ( !std::is_const_v<T> )
struct ConstReferenceMF<MultiArrayView<T,RANK>> : std::type_identity<MultiArrayConstView<T,RANK>> {};
 
 
template <typename T1, typename T2, size_t RANK>
void checkExtents(MultiArrayConstView<T1,RANK> m1, MultiArrayConstView<T2,RANK> m2, std::string message)
{
#ifndef   NDEBUG
  if (m1.extents!=m2.extents) throw std::runtime_error("Extent mismatch in "+message+": "+toStringJSON(m1.extents)+" "+toStringJSON(m2.extents));
#else
  (void)m1; (void)m2; (void)message; // silence unused warnings
#endif // NDEBUG
}
 
template <typename T1, typename T2, size_t RANK>
bool isEqual(MultiArrayConstView<T1,RANK> m1, MultiArrayConstView<T2,RANK> m2)
{
  checkExtents(m1,m2,"MultiArrayView comparison");
  bool res=true;
  for (Extents<RANK> idx{}; idx[0]!=m1.extents[0]; incrementMultiIndex(idx,m1.extents))
    res &= ( m1(idx)==m2(idx) );
  return res;
}
 
 
namespace multiarray {
 
template <size_t RANK>
auto calculateStrides(Extents<RANK> extents)
{
  Extents<RANK> strides; strides[0]=1;
  for (size_t i=1; i<RANK; ++i) strides[i]=strides[i-1]*extents[i-1];
  return strides;
}
 
template <size_t RANK>
auto calculateExtent(Extents<RANK> extents)
{
  return std::ranges::fold_left(extents,1uz,std::multiplies{});
}
 
template <typename T>
inline constexpr auto noInit = [] (size_t e) -> std::vector<T> {return std::vector<T>(e);};
 
template <typename T>
inline constexpr auto zeroInit = [] (size_t e) -> std::vector<T> {return std::vector<T>(e,T(0.));};
 
 
template <typename T>
using StorageType = std::vector<T>;
 
 
template <typename T>
auto copyInit(const std::ranges::sized_range auto& input) {return [&] (size_t e) {
  if (size(input) != e) throw std::runtime_error("Extent mismatch in MultiArray copyInit: "+std::to_string(size(input))+" "+std::to_string(e));
  return StorageType<T>(begin(input),end(input));
};}
 
 
} // multiarray
 
 
template <typename T, size_t RANK>
void checkContiguous(MultiArrayView<T,RANK> mav, std::string message)
{
#ifndef   NDEBUG
  if (mav.offset!=0 || mav.strides!=multiarray::calculateStrides(mav.extents))
    throw std::invalid_argument(message+" needs a full, contiguous array; got offset "+std::to_string(mav.offset)+" and strides "+toStringJSON(mav.strides));
#else
  (void)mav; (void)message; // silence unused warnings
#endif // NDEBUG
}
 
 
 
template <typename T, size_t RANK>
requires ( !std::is_const_v<T> )
class MultiArray : public MultiArrayConstView<T,RANK>
{
public:
  using StorageType = multiarray::StorageType<T>;
  
  MultiArray(const MultiArray&) = delete; MultiArray& operator=(const MultiArray&) = delete;
  
  MultiArray(MultiArray&&) = default; MultiArray& operator=(MultiArray&&) = default;
 
  MultiArray(Extents<RANK> extents, auto initializer) requires requires (size_t e) { { initializer(e) } -> std::convertible_to<StorageType>;}
    : MultiArrayConstView<T,RANK>{extents,multiarray::calculateStrides(extents),0}, data_{initializer(multiarray::calculateExtent(extents))}
  {
    this->dataView=std::span<T>(data_);
  }
 
  explicit MultiArray(Extents<RANK> extents) : MultiArray{extents,multiarray::noInit<T>} {}
 
  // @brief Rank-1 constructor from an existing storage vector (copies the data).
  // explicit MultiArray(const StorageType& st) requires (RANK==1) : MultiArray{{st.size()},multiarray::copyInit<T>(st)} {}
  // removed: It contradicts the library's own semantics, as this constructor is the one place where a silent deep copy hides behind constructor syntax.
 
  friend MultiArray copy(const MultiArray& ma) {return MultiArray{ma.extents, [&] (size_t) {return ma.dataStorage();}};}
 
  template <typename TT>
  void assignTo(MultiArrayConstView<TT,RANK> macv)
  {
#ifndef   NDEBUG
    if (this->dataView.data()==macv.dataView.data())
      throw std::runtime_error("Self assignment attempted in MultiArray");
#endif // NDEBUG
    checkExtents(*this,macv,"MultiArrayView assignment");
 
    for (Extents<RANK> idx{}; idx[0]!=this->extents[0]; incrementMultiIndex(idx,this->extents)) (*this)(idx)=macv(idx);
  }
 
  MultiArrayConstView<T,RANK> constView() const {return *this;}
  
  auto mutableView() {return MultiArrayView<T, RANK>{this->extents,this->strides,0,data_};}
  
  operator MultiArrayView<T,RANK>() {return mutableView();}
 
  T& operator()(auto&&... i) {return const_cast<T&>(static_cast<MultiArrayConstView<T,RANK>>(*this)(FWD(i)...)) ;}
  
  const T& operator()(auto&&... i) const {return static_cast<MultiArrayConstView<T,RANK>>(*this)(FWD(i)...) ;}
 
  StorageType& dataStorage() {return data_;}
  const StorageType& dataStorage() const {return data_;}
 
  friend void tag_invoke( const json::value_from_tag&, json::value& jv, const MultiArray& ma )
  {
    jv = json::object{
      { "extents" , json::value_from(ma.extents) },
      //    { "strides" , ma.strides },
      { "data", json::value_from(ma.data_) }
    };
  }
 
  friend MultiArray tag_invoke( const json::value_to_tag< MultiArray >&, const json::value& jv )
  {
    return MultiArray{
      json::value_to<Extents<RANK>>(jv.as_object().at("extents")),
      multiarray::copyInit<T>( json::value_to<StorageType>(jv.as_object().at("data") ) ) };
  }
  
private:
  StorageType data_;
 
  friend class boost::serialization::access;
  
  template<class Archive> void save(Archive& ar, const unsigned int) const {ar & this->extents & data_;}
 
  template<class Archive> void load(Archive& ar, const unsigned int)
  {
    ar & this->extents & data_; this->strides=multiarray::calculateStrides(this->extents); this->offset=0;
    this->dataView=std::span<T>(data_);
  }
 
  BOOST_SERIALIZATION_SPLIT_MEMBER()
  
};
 
 
namespace multiarray {
 
template <typename T>
MultiArray<T,1> fromStorage(std::vector<T> v)
{
  const size_t s = v.size();
  return MultiArray<T,1>{ {s}, [&] (size_t) {return std::move(v);} };
}
 
} // multiarray
 
template <typename T, size_t RANK>
constexpr auto passByValue_v<MultiArray<T,RANK>> = false;
 
 
} // cppqedutils