arraynd.h

Go to the documentation of this file.

#pragma once
#ifndef SYDEVS_ARRAYND_H_
#define SYDEVS_ARRAYND_H_

#include <sydevs/core/arraynd_base.h>
#include <ostream>

namespace sydevs {

    
template<typename T, int64 ndims>
class arraynd : public arraynd_base<T, ndims>
{
friend class arraynd<T, ndims+1>;
public:
    arraynd();

    arraynd(const std::array<int64, ndims>& dims, const T& value);

    arraynd(const std::array<int64, ndims>& dims, const std::vector<T>& data);

    arraynd(const std::array<int64, ndims>& dims, std::function<T(const std::array<int64, ndims>& indices)> func);
    
    arraynd(const arraynd<T, ndims>&)                      = default;  
    arraynd<T, ndims>& operator=(const arraynd<T, ndims>&) = default;  
    arraynd(arraynd<T, ndims>&&)                           = default;  
    arraynd<T, ndims>& operator=(arraynd<T, ndims>&&)      = default;  
    ~arraynd()                                             = default;  

    const arraynd<T, ndims-1> operator[](int64 index) const;  
    arraynd<T, ndims-1> operator[](int64 index);              

    const arraynd<T, ndims> operator[](range r) const;  
    arraynd<T, ndims> operator[](range r);              

    const arraynd<T, ndims> view() const;  
    arraynd<T, ndims> view();              
    arraynd<T, ndims> copy() const;        

    const arraynd<T, ndims> view_transposed() const;  
    arraynd<T, ndims> view_transposed();              
    arraynd<T, ndims> copy_transposed() const;        

    const arraynd<T, ndims> view_swapped_axes(int64 idim0, int64 idim1) const;  
    arraynd<T, ndims> view_swapped_axes(int64 idim0, int64 idim1);              
    arraynd<T, ndims> copy_swapped_axes(int64 idim0, int64 idim1) const;        

    const arraynd<T, ndims + 1> view_subdivided_axis(int64 idim, const std::array<int64, 2>& dims) const;  
    arraynd<T, ndims + 1> view_subdivided_axis(int64 idim, const std::array<int64, 2>& dims);              
    arraynd<T, ndims + 1> copy_subdivided_axis(int64 idim, const std::array<int64, 2>& dims) const;        

    const arraynd<T, ndims - 1> view_absorbed_axis(int64 idim) const;  
    arraynd<T, ndims - 1> view_absorbed_axis(int64 idim);              
    arraynd<T, ndims - 1> copy_absorbed_axis(int64 idim) const;        

private:
    arraynd(const arraynd<T, ndims+1>& rhs, int64 index, bool readonly);
    arraynd(const arraynd<T, ndims>& rhs, range r, bool readonly);
    arraynd(const arraynd<T, ndims>& rhs, bool is_view, bool readonly);
};


template<typename T>
class arraynd<T, 1> : public arraynd_base<T, 1>
{
friend class arraynd<T, 2>;
public:
    arraynd();

    arraynd(const std::array<int64, 1>& dims, const T& value);

    arraynd(const std::array<int64, 1>& dims, const std::vector<T>& data);

    arraynd(const std::array<int64, 1>& dims, std::function<T(const std::array<int64, 1>& indices)> func);

    arraynd<T, 1>(const arraynd<T, 1>&)            = default;  
    arraynd<T, 1>& operator=(const arraynd<T, 1>&) = default;  
    arraynd<T, 1>(arraynd<T, 1>&&)                 = default;  
    arraynd<T, 1>& operator=(arraynd<T, 1>&&)      = default;  
    ~arraynd<T, 1>()                               = default;  

    const T& operator[](int64 index) const;  
    T& operator[](int64 index);              

    const arraynd<T, 1> operator[](range r) const;  
    arraynd<T, 1> operator[](range r);              

    const arraynd<T, 1> view() const;  
    arraynd<T, 1> view();              
    arraynd<T, 1> copy() const;        

    const arraynd<T, 2> view_subdivided_axis(int64 idim, const std::array<int64, 2>& dims) const;  
    arraynd<T, 2> view_subdivided_axis(int64 idim, const std::array<int64, 2>& dims);              
    arraynd<T, 2> copy_subdivided_axis(int64 idim, const std::array<int64, 2>& dims) const;        

private:
    arraynd(const arraynd<T, 2>& rhs, int64 index, bool readonly);
    arraynd(const arraynd<T, 1>& rhs, range r, bool readonly);
    arraynd(const arraynd<T, 1>& rhs, bool is_view, bool readonly);
};


// n-dimensional array aliases

template<typename T> using array1d = arraynd<T, 1>;
template<typename T> using array2d = arraynd<T, 2>;
template<typename T> using array3d = arraynd<T, 3>;
template<typename T> using array4d = arraynd<T, 4>;
template<typename T> using array5d = arraynd<T, 5>;
template<typename T> using array6d = arraynd<T, 6>;
template<typename T> using array7d = arraynd<T, 7>;
template<typename T> using array8d = arraynd<T, 8>;
template<typename T> using array9d = arraynd<T, 9>;


// n-dimensional array (ndims > 1) implementation

template<typename T, int64 ndims>
arraynd<T, ndims>::arraynd()
    : arraynd_base<T, ndims>()
{
}


template<typename T, int64 ndims>
arraynd<T, ndims>::arraynd(const std::array<int64, ndims>& dims, const T& value)
    : arraynd_base<T, ndims>(dims, value)
{
}


template<typename T, int64 ndims>
arraynd<T, ndims>::arraynd(const std::array<int64, ndims>& dims, const std::vector<T>& data)
    : arraynd_base<T, ndims>(dims, data)
{
}


template<typename T, int64 ndims>
arraynd<T, ndims>::arraynd(const std::array<int64, ndims>& dims, std::function<T(const std::array<int64, ndims>& indices)> func)
    : arraynd_base<T, ndims>(dims, func)
{
}

    
template<typename T, int64 ndims>
const arraynd<T, ndims-1> arraynd<T, ndims>::operator[](int64 index) const
{
    return arraynd<T, ndims-1>(*this, index, true);
}


template<typename T, int64 ndims>
arraynd<T, ndims-1> arraynd<T, ndims>::operator[](int64 index)
{
    return arraynd<T, ndims-1>(*this, index, this->is_readonly());
}


template<typename T, int64 ndims>
const arraynd<T, ndims> arraynd<T, ndims>::operator[](range r) const
{
    return arraynd<T, ndims>(*this, r, true);
}


template<typename T, int64 ndims>
arraynd<T, ndims> arraynd<T, ndims>::operator[](range r)
{
    return arraynd<T, ndims>(*this, r, this->is_readonly());
}


template<typename T, int64 ndims>
const arraynd<T, ndims> arraynd<T, ndims>::view() const
{
    return arraynd<T, ndims>(*this, true, true);
}


template<typename T, int64 ndims>
arraynd<T, ndims> arraynd<T, ndims>::view()
{
    return arraynd<T, ndims>(*this, true, this->is_readonly());
}


template<typename T, int64 ndims>
arraynd<T, ndims> arraynd<T, ndims>::copy() const
{
    return arraynd<T, ndims>(*this, false, false);
}


template<typename T, int64 ndims>
const arraynd<T, ndims> arraynd<T, ndims>::view_transposed() const
{
    auto arr = view();
    this->transpose(arr);
    return arr;
}


template<typename T, int64 ndims>
arraynd<T, ndims> arraynd<T, ndims>::view_transposed()
{
    auto arr = view();
    this->transpose(arr);
    return arr;
}


template<typename T, int64 ndims>
arraynd<T, ndims> arraynd<T, ndims>::copy_transposed() const
{
    auto arr = view();
    this->transpose(arr);
    return arr.copy();
}


template<typename T, int64 ndims>
const arraynd<T, ndims> arraynd<T, ndims>::view_swapped_axes(int64 idim0, int64 idim1) const
{
    auto arr = view();
    this->swap_axes(arr, idim0, idim1);
    return arr;
}


template<typename T, int64 ndims>
arraynd<T, ndims> arraynd<T, ndims>::view_swapped_axes(int64 idim0, int64 idim1)
{
    auto arr = view();
    this->swap_axes(arr, idim0, idim1);
    return arr;
}


template<typename T, int64 ndims>
arraynd<T, ndims> arraynd<T, ndims>::copy_swapped_axes(int64 idim0, int64 idim1) const
{
    auto arr = view();
    this->swap_axes(arr, idim0, idim1);
    return arr.copy();
}

    
template<typename T, int64 ndims>
const arraynd<T, ndims + 1> arraynd<T, ndims>::view_subdivided_axis(int64 idim, const std::array<int64, 2>& dims) const
{
    auto arr = arraynd<T, ndims + 1>();
    this->subdivide_axis(*this, arr, idim, dims);
    return arr;
}


template<typename T, int64 ndims>
arraynd<T, ndims + 1> arraynd<T, ndims>::view_subdivided_axis(int64 idim, const std::array<int64, 2>& dims)
{
    auto arr = arraynd<T, ndims + 1>();
    this->subdivide_axis(*this, arr, idim, dims);
    return arr;
}


template<typename T, int64 ndims>
arraynd<T, ndims + 1> arraynd<T, ndims>::copy_subdivided_axis(int64 idim, const std::array<int64, 2>& dims) const
{
    auto arr0 = copy();
    return arr0.view_subdivided_axis(idim, dims);
}


template<typename T, int64 ndims>
const arraynd<T, ndims - 1> arraynd<T, ndims>::view_absorbed_axis(int64 idim) const
{
    auto arr = arraynd<T, ndims - 1>(std::array<int64, ndims - 1>(), std::vector<T>());
    this->absorb_axis(*this, arr, idim);
    return arr;
}


template<typename T, int64 ndims>
arraynd<T, ndims - 1> arraynd<T, ndims>::view_absorbed_axis(int64 idim)
{
    auto arr = arraynd<T, ndims - 1>(std::array<int64, ndims - 1>(), std::vector<T>());
    this->absorb_axis(*this, arr, idim);
    return arr;
}


template<typename T, int64 ndims>
arraynd<T, ndims - 1> arraynd<T, ndims>::copy_absorbed_axis(int64 idim) const
{
    auto arr0 = copy();
    return arr0.view_absorbed_axis(idim);
}


template<typename T, int64 ndims>
arraynd<T, ndims>::arraynd(const arraynd<T, ndims+1>& rhs, int64 index, bool readonly)
    : arraynd_base<T, ndims>(rhs, index, readonly)
{
}


template<typename T, int64 ndims>
arraynd<T, ndims>::arraynd(const arraynd<T, ndims>& rhs, range r, bool readonly)
    : arraynd_base<T, ndims>(rhs, r, readonly)
{
}

    
template<typename T, int64 ndims>
arraynd<T, ndims>::arraynd(const arraynd<T, ndims>& rhs, bool is_view, bool readonly)
    : arraynd_base<T, ndims>(rhs, is_view, readonly)
{
}

    
// 1-dimensional array implementation

template<typename T>
arraynd<T, 1>::arraynd()
    : arraynd_base<T, 1>()
{
}


template<typename T>
arraynd<T, 1>::arraynd(const std::array<int64, 1>& dims, const T& value)
    : arraynd_base<T, 1>(dims, value)
{
}


template<typename T>
arraynd<T, 1>::arraynd(const std::array<int64, 1>& dims, const std::vector<T>& data)
    : arraynd_base<T, 1>(dims, data)
{
}


template<typename T>
arraynd<T, 1>::arraynd(const std::array<int64, 1>& dims, std::function<T(const std::array<int64, 1>& indices)> func)
    : arraynd_base<T, 1>(dims, func)
{
}


template<typename T>
const T& arraynd<T, 1>::operator[](int64 index) const
{
    return arraynd_base<T, 1>::data()[arraynd_base<T, 1>::offset() + index*arraynd_base<T, 1>::strides()[0]];
}


template<typename T>
T& arraynd<T, 1>::operator[](int64 index)
{
    if (this->is_readonly()) throw std::logic_error("Attempt to obtain a non-const reference to readonly multidimensional array data");
    return arraynd_base<T, 1>::data()[arraynd_base<T, 1>::offset() + index*arraynd_base<T, 1>::strides()[0]];
}


template<typename T>
const arraynd<T, 1> arraynd<T, 1>::operator[](range r) const
{
    return arraynd<T, 1>(*this, r, true);
}


template<typename T>
arraynd<T, 1> arraynd<T, 1>::operator[](range r)
{
    return arraynd<T, 1>(*this, r, this->is_readonly());
}


template<typename T>
const arraynd<T, 1> arraynd<T, 1>::view() const
{
    return arraynd<T, 1>(*this, true, true);
}


template<typename T>
arraynd<T, 1> arraynd<T, 1>::view()
{
    return arraynd<T, 1>(*this, true, this->is_readonly());
}


template<typename T>
arraynd<T, 1> arraynd<T, 1>::copy() const
{
    return arraynd<T, 1>(*this, false, false);
}


template<typename T>
const arraynd<T, 2> arraynd<T, 1>::view_subdivided_axis(int64 idim, const std::array<int64, 2>& dims) const
{
    auto arr = arraynd<T, 2>();
    this->subdivide_axis(*this, arr, idim, dims);
    return arr;
}


template<typename T>
arraynd<T, 2> arraynd<T, 1>::view_subdivided_axis(int64 idim, const std::array<int64, 2>& dims)
{
    auto arr = arraynd<T, 2>();
    this->subdivide_axis(*this, arr, idim, dims);
    return arr;
}


template<typename T>
arraynd<T, 2> arraynd<T, 1>::copy_subdivided_axis(int64 idim, const std::array<int64, 2>& dims) const
{
    auto arr0 = copy();
    return arr0.view_subdivided_axis(idim, dims);
}


template<typename T>
arraynd<T, 1>::arraynd(const arraynd<T, 2>& rhs, int64 index, bool readonly)
    : arraynd_base<T, 1>(rhs, index, readonly)
{
}


template<typename T>
arraynd<T, 1>::arraynd(const arraynd<T, 1>& rhs, range r, bool readonly)
    : arraynd_base<T, 1>(rhs, r, readonly)
{
}

    
template<typename T>
arraynd<T, 1>::arraynd(const arraynd<T, 1>& rhs, bool is_view, bool readonly)
    : arraynd_base<T, 1>(rhs, is_view, readonly)
{
}


// multi-dimensional array operations

template<typename T, int64 ndims>
bool aligned(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    auto okay = true;
    for (int64 idim = 0; okay && idim < ndims; ++idim) {
        okay = (lhs.dims()[idim] == rhs.dims()[idim]);
    }
    return okay;
}


template<typename T, int64 ndims>
arraynd<T, ndims> operator+(const arraynd<T, ndims>& rhs)
{
    return rhs.copy();
}


template<typename T, int64 ndims>
arraynd<T, ndims> operator-(const arraynd<T, ndims>& rhs)
{
    return arraynd<T, ndims>(rhs.dims(), [&rhs](const std::array<int64, ndims>& indices) {
        return -rhs(indices);
    });
}


template<typename T, int64 ndims>
arraynd<T, ndims> operator+(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply + operator to multidimensional arrays with inconsistent dimensions");    
    return arraynd<T, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) + rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> operator+(const arraynd<T, ndims>& lhs, const U& rhs)
{
    return arraynd<T, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) + rhs;
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> operator+(const U& lhs, const arraynd<T, ndims>& rhs)
{
    return arraynd<T, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs + rhs(indices);
    });
}


template<typename T, int64 ndims>
arraynd<T, ndims> operator-(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply - operator to multidimensional arrays with inconsistent dimensions");    
    return arraynd<T, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) - rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> operator-(const arraynd<T, ndims>& lhs, const U& rhs)
{
    return arraynd<T, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) - rhs;
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> operator-(const U& lhs, const arraynd<T, ndims>& rhs)
{
    return arraynd<T, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs - rhs(indices);
    });
}


template<typename T, int64 ndims>
arraynd<T, ndims> operator*(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply * operator to multidimensional arrays with inconsistent dimensions");    
    return arraynd<T, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices)*rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> operator*(const arraynd<T, ndims>& lhs, const U& rhs)
{
    return arraynd<T, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices)*rhs;
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> operator*(const U& lhs, const arraynd<T, ndims>& rhs)
{
    return arraynd<T, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs*rhs(indices);
    });
}


template<typename T, int64 ndims>
arraynd<T, ndims> operator/(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply / operator to multidimensional arrays with inconsistent dimensions");    
    return arraynd<T, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices)/rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> operator/(const arraynd<T, ndims>& lhs, const U& rhs)
{
    return arraynd<T, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices)/rhs;
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> operator/(const U& lhs, const arraynd<T, ndims>& rhs)
{
    return arraynd<T, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs/rhs(indices);
    });
}


template<int64 ndims>
arraynd<bool, ndims> operator!(const arraynd<bool, ndims>& arr)
{
    return arraynd<bool, ndims>(arr.dims(), [&arr](const std::array<int64, ndims>& indices) {
        return !arr(indices);
    });
}


template<int64 ndims>
bool all(const arraynd<bool, ndims>& arr)
{
    bool conjunction = true;
    arr.traverse([&conjunction](const std::array<int64, ndims>& indices, const bool& value) {
        conjunction = (conjunction && value);
        return conjunction;
    });
    return conjunction;
}


template<int64 ndims>
bool any(const arraynd<bool, ndims>& arr)
{
    bool disjunction = false;
    arr.traverse([&disjunction](const std::array<int64, ndims>& indices, const bool& value) {
        disjunction = (disjunction || value);
        return !disjunction;
    });
    return disjunction;
}


template<typename T, int64 ndims>
arraynd<bool, ndims> operator<(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply < operator to multidimensional arrays with inconsistent dimensions");
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) < rhs(indices);
    });
}


template<typename T, int64 ndims>
arraynd<bool, ndims> operator>(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply > operator to multidimensional arrays with inconsistent dimensions");
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) > rhs(indices);
    });
}


template<typename T, int64 ndims>
arraynd<bool, ndims> operator<=(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply <= operator to multidimensional arrays with inconsistent dimensions");
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) <= rhs(indices);
    });
}


template<typename T, int64 ndims>
arraynd<bool, ndims> operator>=(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply >= operator to multidimensional arrays with inconsistent dimensions");
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) >= rhs(indices);
    });
}


template<typename T, int64 ndims>
arraynd<bool, ndims> operator==(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply == operator to multidimensional arrays with inconsistent dimensions");
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) == rhs(indices);
    });
}


template<typename T, int64 ndims>
arraynd<bool, ndims> operator!=(const arraynd<T, ndims>& lhs, const arraynd<T, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply != operator to multidimensional arrays with inconsistent dimensions");
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) != rhs(indices);
    });
}


template<int64 ndims>
arraynd<bool, ndims> operator&&(const arraynd<bool, ndims>& lhs, const arraynd<bool, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply && operator to multidimensional arrays with inconsistent dimensions");
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) && rhs(indices);
    });
}


template<int64 ndims>
arraynd<bool, ndims> operator||(const arraynd<bool, ndims>& lhs, const arraynd<bool, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to apply || operator to multidimensional arrays with inconsistent dimensions");
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) || rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator<(const arraynd<T, ndims>& lhs, const U& rhs)
{
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) < rhs;
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator>(const arraynd<T, ndims>& lhs, const U& rhs)
{
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) > rhs;
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator<=(const arraynd<T, ndims>& lhs, const U& rhs)
{
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) <= rhs;
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator>=(const arraynd<T, ndims>& lhs, const U& rhs)
{
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) >= rhs;
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator==(const arraynd<T, ndims>& lhs, const U& rhs)
{
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) == rhs;
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator!=(const arraynd<T, ndims>& lhs, const U& rhs)
{
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) != rhs;
    });
}


template<int64 ndims>
arraynd<bool, ndims> operator&&(const arraynd<bool, ndims>& lhs, const bool& rhs)
{
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) && rhs;
    });
}


template<int64 ndims>
arraynd<bool, ndims> operator||(const arraynd<bool, ndims>& lhs, const bool& rhs)
{
    return arraynd<bool, ndims>(lhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs(indices) || rhs;
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator<(const U& lhs, const arraynd<T, ndims>& rhs)
{
    return arraynd<bool, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs < rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator>(const U& lhs, const arraynd<T, ndims>& rhs)
{
    return arraynd<bool, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs > rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator<=(const U& lhs, const arraynd<T, ndims>& rhs)
{
    return arraynd<bool, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs <= rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator>=(const U& lhs, const arraynd<T, ndims>& rhs)
{
    return arraynd<bool, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs >= rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator==(const U& lhs, const arraynd<T, ndims>& rhs)
{
    return arraynd<bool, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs == rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<bool, ndims> operator!=(const U& lhs, const arraynd<T, ndims>& rhs)
{
    return arraynd<bool, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs != rhs(indices);
    });
}


template<int64 ndims>
arraynd<bool, ndims> operator&&(const bool& lhs, const arraynd<bool, ndims>& rhs)
{
    return arraynd<bool, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs && rhs(indices);
    });
}


template<int64 ndims>
arraynd<bool, ndims> operator||(const bool& lhs, const arraynd<bool, ndims>& rhs)
{
    return arraynd<bool, ndims>(rhs.dims(), [&lhs, &rhs](const std::array<int64, ndims>& indices) {
        return lhs || rhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> to(const arraynd<U, ndims>& arr)
{
    return arraynd<T, ndims>(arr.dims(), [&arr](const std::array<int64, ndims>& indices) {
        return static_cast<T>(arr(indices));
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> replace(const arraynd<T, ndims>& lhs, const arraynd<bool, ndims>& selection, const U& rhs)
{
    return arraynd<T, ndims>(lhs.dims(), [&lhs, &selection, &rhs](const std::array<int64, ndims>& indices) {
        return selection(indices) ? rhs : lhs(indices);
    });
}


template<typename T, int64 ndims, typename U>
arraynd<T, ndims> replace(const arraynd<T, ndims>& lhs, const arraynd<bool, ndims>& selection, const arraynd<U, ndims>& rhs)
{
    if (!aligned(lhs, rhs)) throw std::domain_error("Attempt to replace values using multidimensional arrays with inconsistent dimensions");
    return arraynd<T, ndims>(lhs.dims(), [&lhs, &selection, &rhs](const std::array<int64, ndims>& indices) {
        return selection(indices) ? rhs(indices) : lhs(indices);
    });
}


template<typename T, int64 ndims>
std::ostream& operator<<(std::ostream& os, const arraynd<T, ndims>& rhs)
{
    auto indices = std::array<int64, ndims>();
    for (int64 idim = 0; idim < ndims; ++idim) {
        indices[idim] = 0;
    }
    int64 offset = rhs.offset();
    int64 idim = 0;
    os << "{";
    while (idim >= 0) {
        if (indices[idim] < rhs.dims()[idim]) {
            // The index is not beyond the last element on the current axis.
            // Continue outputting the current sub-array.
            if (indices[idim] > 0) {
                // This is not the first element on the current axis.
                // Output a delimiter.
                os << ", ";
            }
            if (idim == ndims - 1) {
                // The current axis is the last one.
                // Output the current element and advance the index.
                os << rhs.data()[offset];
                ++indices[idim];
                offset += rhs.strides()[idim];
            }
            else {
                // The current axis is not the last one.
                // Start outputting a new sub-array.
                ++idim;
                os << "{";
            }
        }
        else {
            // The index is beyond the last element on the current axis.
            // Finish outputting the current sub-array.
            indices[idim] = 0;
            offset -= rhs.dims()[idim]*rhs.strides()[idim];
            --idim;
            if (idim >= 0) {
                ++indices[idim];
                offset += rhs.strides()[idim];
                os << "}";
            }
        }
    }
    os << "}";
    return os;
}


}  // namespace


namespace std {
template<typename T>
struct less<sydevs::array1d<T>>
{
    bool operator()(const sydevs::array1d<T>& lhs, const sydevs::array1d<T>& rhs) const {
        return std::lexicographical_compare(lhs.data(), lhs.data() + lhs.size(), 
                                            rhs.data(), rhs.data() + rhs.size());
    }
};
}


#endif