#ifndef included_ArrayClass_hpp
#define included_ArrayClass_hpp

#include "common/Array.h"
#include "common/FunctionTable.h"
#include "common/FunctionTable.hpp"
#include "common/Utilities.h"

#include <algorithm>
#include <cmath>
#include <complex>
#include <cstring>
#include <limits>
#include <memory>


/********************************************************
 *  External instantiations                              *
 ********************************************************/
extern template class Array<char>;
extern template class Array<uint8_t>;
extern template class Array<uint16_t>;
extern template class Array<uint32_t>;
extern template class Array<uint64_t>;
extern template class Array<int8_t>;
extern template class Array<int16_t>;
extern template class Array<int32_t>;
extern template class Array<int64_t>;
extern template class Array<double>;
extern template class Array<float>;


/********************************************************
 *  Macros to help instantiate functions                 *
 ********************************************************/
// clang-format off
#define instantiateArrayConstructors( TYPE )                                    \
    template Array<TYPE>::Array();                                              \
    template Array<TYPE>::~Array();                                             \
    template Array<TYPE>::Array( const ArraySize & );                           \
    template Array<TYPE>::Array( size_t );                                      \
    template Array<TYPE>::Array( size_t, size_t );                              \
    template Array<TYPE>::Array( size_t, size_t, size_t );                      \
    template Array<TYPE>::Array( size_t, size_t, size_t, size_t );              \
    template Array<TYPE>::Array( size_t, size_t, size_t, size_t, size_t );      \
    template Array<TYPE>::Array( const std::vector<size_t> &, const TYPE * );   \
    template Array<TYPE>::Array( std::initializer_list<TYPE> );                 \
    template Array<TYPE>::Array( std::initializer_list<std::initializer_list<TYPE>> ); \
    template Array<TYPE>::Array( const Array<TYPE> & );                         \
    template Array<TYPE>::Array( Array<TYPE> && );                              \
    template void Array<TYPE>::reshape( ArraySize const& );                     \
    template void Array<TYPE>::squeeze();                                       \
    template std::unique_ptr<const Array<TYPE>>                                 \
        Array<TYPE>::constView(ArraySize const&, std::shared_ptr<TYPE const> const&); \
    template void Array<TYPE>::viewRaw( ArraySize const&, TYPE*, bool, bool );  \
    template void Array<TYPE>::view2(ArraySize const&, std::shared_ptr<TYPE> ); \
    template Array<TYPE> &Array<TYPE>::operator=( const Array<TYPE> & );        \
    template Array<TYPE> &Array<TYPE>::operator=( Array<TYPE> && );
// clang-format on


/********************************************************
 *  Constructors                                         *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array()
    : d_isCopyable( true ), d_isFixedSize( false ), d_data( nullptr )
{
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( const ArraySize &N )
    : d_isCopyable( true ), d_isFixedSize( false )
{
    allocate( N );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( size_t N ) : d_isCopyable( true ), d_isFixedSize( false )
{
    allocate( ArraySize( N ) );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( size_t N_rows, size_t N_cols )
    : d_isCopyable( true ), d_isFixedSize( false )
{
    allocate( ArraySize( N_rows, N_cols ) );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( size_t N1, size_t N2, size_t N3 )
    : d_isCopyable( true ), d_isFixedSize( false )
{
    allocate( ArraySize( N1, N2, N3 ) );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( size_t N1, size_t N2, size_t N3, size_t N4 )
    : d_isCopyable( true ), d_isFixedSize( false )
{
    allocate( ArraySize( N1, N2, N3, N4 ) );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( size_t N1, size_t N2, size_t N3, size_t N4, size_t N5 )
    : d_isCopyable( true ), d_isFixedSize( false )
{
    allocate( ArraySize( N1, N2, N3, N4, N5 ) );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( const std::vector<size_t> &N, const TYPE *data )
    : d_isCopyable( true ), d_isFixedSize( false )
{
    allocate( N );
    if ( data )
        copy( data );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( std::string str ) : d_isCopyable( true ), d_isFixedSize( false )
{
    allocate( 0 );
    if ( (int) std::count( str.begin(), str.end(), ' ' ) == (int) str.length() ) {
        // Empty string
        return;
    }
    // Remove unnecessary whitespace
    while ( str.front() == ' ' )
        str.erase( 0, 1 );
    while ( str.back() == ' ' )
        str.resize( str.length() - 1 );
    while ( str.find( ',' ) != std::string::npos )
        str[str.find( ',' )] = ' ';
    while ( str.find( "  " ) != std::string::npos )
        str.replace( str.find( "  " ), 2, " " );
    // Check if the string is of the format [...]
    if ( str.front() == '[' && str.back() == ']' ) {
        str.erase( 0, 1 );
        str.resize( str.length() - 1 );
        *this = Array( str );
        return;
    }
    // Check if we are dealing with a 2D array
    if ( str.find( ';' ) != std::string::npos ) {
        size_t i1 = 0;
        size_t i2 = str.find( ';' );
        std::vector<Array> x;
        while ( i2 > i1 ) {
            auto tmp = str.substr( i1, i2 - i1 );
            x.emplace_back( Array( tmp ) );
            i1 = i2 + 1;
            i2 = str.find( ';', i1 + 1 );
            if ( i2 == std::string::npos )
                i2 = str.length();
        }
        for ( auto &y : x )
            y.reshape( { 1, y.length() } );
        *this = cat( x, 0 );
        return;
    }
    // Begin parsing the array constructor
    size_t i1 = 0;
    size_t i2 = str.find( ' ' );
    std::vector<TYPE> data;
    while ( i2 > i1 ) {
        auto tmp = str.substr( i1, i2 - i1 );
        int type = std::count( tmp.begin(), tmp.end(), ':' );
        if ( type == 0 ) {
            data.push_back( std::stod( tmp ) );
        } else if ( type == 1 ) {
            size_t k   = tmp.find( ':' );
            TYPE x1    = std::stod( tmp.substr( 0, k ) );
            TYPE x2    = std::stod( tmp.substr( k + 1 ) );
            double tol = 1e-8 * ( x2 - x1 );
            for ( TYPE x = x1; x <= x2 + tol; x += 1 )
                data.push_back( x );
        } else if ( type == 2 ) {
            size_t k1  = tmp.find( ':' );
            size_t k2  = tmp.find( ':', k1 + 1 );
            TYPE x1    = std::stod( tmp.substr( 0, k1 ) );
            TYPE dx    = std::stod( tmp.substr( k1 + 1, k2 - k1 - 1 ) );
            TYPE x2    = std::stod( tmp.substr( k2 + 1 ) );
            double tol = 1e-8 * ( x2 - x1 );
            for ( TYPE x = x1; x <= x2 + tol; x += dx )
                data.push_back( x );
        } else {
            throw std::logic_error( "Failed to parse string constructor: " + str );
        }
        i1 = i2;
        i2 = str.find( ' ', i1 + 1 );
        if ( i2 == std::string::npos )
            i2 = str.length();
    }
    allocate( data.size() );
    if ( std::is_same<TYPE, bool>::value ) {
        for ( size_t i = 0; i < data.size(); i++ )
            d_data[i] = data[i];
    } else {
        copy( data.data() );
    }
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( std::initializer_list<TYPE> x )
    : d_isCopyable( true ), d_isFixedSize( false )
{
    allocate( { x.size() } );
    auto it = x.begin();
    for ( size_t i = 0; i < x.size(); ++i, ++it )
        d_data[i] = *it;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( std::initializer_list<std::initializer_list<TYPE>> x )
    : d_isCopyable( true ), d_isFixedSize( false )
{
    size_t Nx = x.size();
    size_t Ny = 0;
    for ( const auto y : x )
        Ny = std::max<size_t>( Ny, y.size() );
    allocate( { Nx, Ny } );
    auto itx = x.begin();
    for ( size_t i = 0; i < x.size(); ++i, ++itx ) {
        auto ity = itx->begin();
        for ( size_t j = 0; j < itx->size(); ++j, ++ity ) {
            d_data[i + j * Nx] = *ity;
        }
    }
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::allocate( const ArraySize &N )
{
    if ( d_isFixedSize )
        throw std::logic_error( "Array cannot be resized" );
    d_size      = N;
    auto length = d_size.length();
    d_data      = nullptr;
    if ( length > 0 ) {
        try {
            d_data = new TYPE[length];
        } catch ( ... ) {
            throw std::logic_error( "Failed to allocate array" );
        }
    }
    d_ptr.reset( d_data, []( TYPE *p ) { delete[] p; } );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( const Array &rhs )
    : d_isCopyable( true ), d_isFixedSize( false )
{
    if ( !rhs.d_isCopyable )
        throw std::logic_error( "Array cannot be copied" );
    allocate( rhs.size() );
    copy( rhs.d_data );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::Array( Array &&rhs )
    : d_isCopyable( rhs.d_isCopyable ),
      d_isFixedSize( rhs.d_isFixedSize ),
      d_size( rhs.d_size ),
      d_data( rhs.d_data ),
      d_ptr( std::move( rhs.d_ptr ) )
{
    rhs.d_size = ArraySize();
    rhs.d_data = nullptr;
    rhs.d_ptr  = nullptr;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> &Array<TYPE, FUN, Allocator>::operator=( const Array &rhs )
{
    if ( this == &rhs )
        return *this;
    if ( !rhs.d_isCopyable )
        throw std::logic_error( "Array cannot be copied" );
    allocate( rhs.size() );
    copy( rhs.d_data );
    return *this;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> &Array<TYPE, FUN, Allocator>::operator=( Array &&rhs )
{
    if ( this == &rhs )
        return *this;
    d_isCopyable  = rhs.d_isCopyable;
    d_isFixedSize = rhs.d_isFixedSize;
    d_size        = rhs.d_size;
    d_data        = rhs.d_data;
    d_ptr         = std::move( rhs.d_ptr );
    rhs.d_size    = ArraySize();
    rhs.d_data    = nullptr;
    rhs.d_ptr     = nullptr;
    return *this;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> &Array<TYPE, FUN, Allocator>::operator=( const std::vector<TYPE> &rhs )
{
    allocate( ArraySize( rhs.size() ) );
    for ( size_t i = 0; i < rhs.size(); i++ )
        d_data[i] = rhs[i];
    return *this;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>::~Array()
{
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::clear()
{
    d_isCopyable  = true;
    d_isFixedSize = false;
    d_size        = ArraySize();
    d_ptr.reset();
    d_data = nullptr;
}


/********************************************************
 *  Copy/move values from one array to another (resize)  *
 ********************************************************/
template<class TYPE>
static inline void moveValues( const ArraySize &N1, const ArraySize &N2, TYPE *data1, TYPE *data2 )
{
    for ( size_t i5 = 0; i5 < std::min( N1[4], N2[4] ); i5++ ) {
        for ( size_t i4 = 0; i4 < std::min( N1[3], N2[3] ); i4++ ) {
            for ( size_t i3 = 0; i3 < std::min( N1[2], N2[2] ); i3++ ) {
                for ( size_t i2 = 0; i2 < std::min( N1[1], N2[1] ); i2++ ) {
                    for ( size_t i1 = 0; i1 < std::min( N1[0], N2[0] ); i1++ ) {
                        size_t index1 = N1.index( i1, i2, i3, i4, i5 );
                        size_t index2 = N2.index( i1, i2, i3, i4, i5 );
                        data2[index2] = std::move( data1[index1] );
                    }
                }
            }
        }
    }
}
template<class TYPE>
static inline void
copyValues( const ArraySize &N1, const ArraySize &N2, const TYPE *data1, TYPE *data2 )
{
    for ( size_t i5 = 0; i5 < std::min( N1[4], N2[4] ); i5++ ) {
        for ( size_t i4 = 0; i4 < std::min( N1[3], N2[3] ); i4++ ) {
            for ( size_t i3 = 0; i3 < std::min( N1[2], N2[2] ); i3++ ) {
                for ( size_t i2 = 0; i2 < std::min( N1[1], N2[1] ); i2++ ) {
                    for ( size_t i1 = 0; i1 < std::min( N1[0], N2[0] ); i1++ ) {
                        size_t index1 = N1.index( i1, i2, i3, i4, i5 );
                        size_t index2 = N2.index( i1, i2, i3, i4, i5 );
                        data2[index2] = data1[index1];
                    }
                }
            }
        }
    }
}


/********************************************************
 *  Resize the array                                     *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::resize( const ArraySize &N )
{
    // Check if the array actually changed size
    bool equal = true;
    for ( size_t i = 0; i < ArraySize::maxDim(); i++ )
        equal = equal && N[i] == d_size[i];
    if ( equal ) {
        d_size = N;
        return;
    }
    // Store the old data
    auto N0    = d_size;
    auto data0 = d_ptr;
    // Allocate new data
    allocate( N );
    // Copy the old values
    if ( N.length() > 0 && d_size.length() > 0 ) {
        if ( data0.use_count() <= 1 ) {
            // We own the data, use std:move
            moveValues( N0, N, data0.get(), d_data );
        } else if ( std::is_copy_constructible<TYPE>::value ) {
            // We do not own the data, copy
            copyValues( N0, N, data0.get(), d_data );
        } else {
            throw std::logic_error( "No copy constructor" );
        }
    }
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::resizeDim( int dim, size_t N, const TYPE &value )
{
    if ( dim < 0 || dim > d_size.ndim() )
        throw std::out_of_range( "Invalid dimension" );
    size_t N0 = d_size[dim];
    auto size = d_size;
    size.resize( dim, N );
    resize( size );
    size_t n1 = 1, n2 = 1;
    for ( int d = 0; d < dim; d++ )
        n1 *= size[d];
    for ( size_t d = dim + 1; d < size.ndim(); d++ )
        n2 *= size[d];
    for ( size_t k = 0; k < n2; k++ ) {
        for ( size_t j = N0; j < N; j++ ) {
            for ( size_t i = 0; i < n1; i++ ) {
                d_data[i + j * n1 + k * n1 * N] = value;
            }
        }
    }
}


/********************************************************
 *  Reshape/squeeze the array                            *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::reshape( const ArraySize &N )
{
    if ( N.length() != d_size.length() )
        throw std::logic_error( "reshape is not allowed to change the array size" );
    d_size = N;
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::squeeze()
{
    d_size.squeeze();
}


/********************************************************
 *  Shift/permute the array                              *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::shiftDim( int N ) const
{
    if ( N > 0 )
        N = N % d_size.ndim();
    if ( N == 0 ) {
        // No shift required
        return *this;
    } else if ( N > 0 ) {
        // Shift to the left and wrap
        std::vector<uint8_t> index( d_size.ndim() );
        size_t i = 0;
        for ( size_t j=N; j<index.size(); j++, i++)
            index[i] = j;
        for ( size_t j=0; i<index.size(); j++, i++)
            index[i] = j;
        return permute( index );
    } else  {
        // Shift to the right (padding with singletons)
        N = -N;
        ASSERT( d_size.ndim() + N < (int) ArraySize::maxDim() );
        size_t dims[10] = { 1 };
        for ( int i = 0; i < ndim(); i++ )
            dims[N+i] = d_size[i];
        auto y = *this;
        y.reshape( ArraySize( ndim() + N, dims ) );
        return y;
    }
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::permute( const std::vector<uint8_t> &index ) const
{
    // Check the permutation
    ASSERT( (int) index.size() == ndim() );
    for ( int i=0; i < ndim(); i++) {
        ASSERT( index[i] < ndim() );
        for ( int j=0; j < i; j++)
            ASSERT( index[i] != index[j] );
    }
    // Create the new Array
    size_t dims[5] = { 1u, 1u, 1u, 1u, 1u };
    for ( size_t i=0; i<index.size(); i++)
        dims[i] = d_size[index[i]];
    Array y( ArraySize( ndim(), dims ) );
    y.fill( -1 );
    ASSERT( y.length() == this->length() );
    // Fill the data
    size_t N[5] = { 1u, 1u, 1u, 1u, 1u };
    for ( int i=0; i < ndim(); i++) {
        std::array<size_t, 5> ijk = { 0, 0, 0, 0, 0 };
        ijk[index[i]] = 1;
        N[i] = d_size.index( ijk );
    }
    size_t tmp = ( dims[0] - 1 ) * N[0] + ( dims[1] - 1 ) * N[1] + ( dims[2] - 1 ) * N[2] + ( dims[3] - 1 ) * N[3] + ( dims[4] - 1 ) * N[4] + 1;
    ASSERT( tmp == length() );
    for ( size_t i4 = 0; i4 < dims[4]; i4++ ) {
        for ( size_t i3 = 0; i3 < dims[3]; i3++ ) {
            for ( size_t i2 = 0; i2 < dims[2]; i2++ ) {
                for ( size_t i1 = 0; i1 < dims[1]; i1++ ) {
                    for ( size_t i0 = 0; i0 < dims[0]; i0++ ) {
                        size_t index2 = i0 * N[0] + i1 * N[1] + i2 * N[2] + i3 * N[3] + i4 * N[4];
                        y( i0, i1, i2, i3, i4 ) = d_data[index2];
                    }
                }
            }
        }
    }
    return y;
}


/********************************************************
 *  Subset the array                                     *
 ********************************************************/
// Helper function to check subset indices
template<class TYPE, class FUN, class Allocator>
inline void
Array<TYPE, FUN, Allocator>::checkSubsetIndex( const std::vector<Range<size_t>> &range ) const
{
    bool test = (int) range.size() == d_size.ndim();
    for ( size_t d = 0; d < range.size(); d++ )
        test = test && range[d].j <= d_size[d];
    if ( !test )
        throw std::logic_error( "indices for subset are invalid" );
}
template<class TYPE, class FUN, class Allocator>
std::vector<Range<size_t>>
Array<TYPE, FUN, Allocator>::convert( const std::vector<size_t> &index ) const
{
    std::vector<Range<size_t>> range( d_size.ndim() );
    if ( index.size() % 2 != 0 || static_cast<int>( index.size() / 2 ) < d_size.ndim() )
        throw std::logic_error( "indices for subset are invalid" );
    for ( int d = 0; d < d_size.ndim(); d++ )
        range[d] = Range<size_t>( index[2 * d + 0], index[2 * d + 1] );
    return range;
}
// Helper function to return dimensions for the subset array
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::getSubsetArrays( const std::vector<Range<size_t>> &index,
                                                   std::array<size_t, 5> &first,
                                                   std::array<size_t, 5> &last,
                                                   std::array<size_t, 5> &inc,
                                                   std::array<size_t, 5> &N )
{
    first.fill( 0 );
    last.fill( 0 );
    inc.fill( 1 );
    N.fill( 1 );
    size_t ndim = index.size();
    for ( size_t d = 0; d < ndim; d++ ) {
        first[d] = index[d].i;
        last[d]  = index[d].j;
        inc[d]   = index[d].k;
        N[d]     = ( last[d] - first[d] + inc[d] ) / inc[d];
    }
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>
Array<TYPE, FUN, Allocator>::subset( const std::vector<Range<size_t>> &index ) const
{
    // Get the subset indicies
    checkSubsetIndex( index );
    std::array<size_t, 5> first, last, inc, N1;
    getSubsetArrays( index, first, last, inc, N1 );
    ArraySize S1( d_size.ndim(), N1.data() );
    // Create the new array
    Array<TYPE, FUN, Allocator> subset_array( S1 );
    // Fill the new array
    TYPE *subset_data = subset_array.data();
    for ( size_t i4 = first[4], k1 = 0; i4 <= last[4]; i4 += inc[4] ) {
        for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) {
            for ( size_t i2 = first[2]; i2 <= last[2]; i2 += inc[2] ) {
                for ( size_t i1 = first[1]; i1 <= last[1]; i1 += inc[1] ) {
                    for ( size_t i0 = first[0]; i0 <= last[0]; i0 += inc[0], k1++ ) {
                        size_t k2       = d_size.index( i0, i1, i2, i3, i4 );
                        subset_data[k1] = d_data[k2];
                    }
                }
            }
        }
    }
    return subset_array;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>
Array<TYPE, FUN, Allocator>::subset( const std::vector<size_t> &index ) const
{
    auto range = convert( index );
    return subset( range );
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::copySubset( const std::vector<Range<size_t>> &index,
                                              const Array<TYPE, FUN, Allocator> &subset )
{
    // Get the subset indices
    checkSubsetIndex( index );
    std::array<size_t, 5> first, last, inc, N1;
    getSubsetArrays( index, first, last, inc, N1 );
    // Copy the sub-array
    const TYPE *src_data = subset.data();
    for ( size_t i4 = first[4], k1 = 0; i4 <= last[4]; i4 += inc[4] ) {
        for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) {
            for ( size_t i2 = first[2]; i2 <= last[2]; i2 += inc[2] ) {
                for ( size_t i1 = first[1]; i1 <= last[1]; i1 += inc[1] ) {
                    for ( size_t i0 = first[0]; i0 <= last[0]; i0 += inc[0], k1++ ) {
                        size_t k2  = d_size.index( i0, i1, i2, i3, i4 );
                        d_data[k2] = src_data[k1];
                    }
                }
            }
        }
    }
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::addSubset( const std::vector<Range<size_t>> &index,
                                             const Array<TYPE, FUN, Allocator> &subset )
{
    // Get the subset indices
    checkSubsetIndex( index );
    std::array<size_t, 5> first, last, inc, N1;
    getSubsetArrays( index, first, last, inc, N1 );
    // add the sub-array
    for ( size_t i4 = first[4], k1 = 0; i4 <= last[4]; i4 += inc[4] ) {
        for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) {
            for ( size_t i2 = first[2]; i2 <= last[2]; i2 += inc[2] ) {
                for ( size_t i1 = first[1]; i1 <= last[1]; i1 += inc[1] ) {
                    for ( size_t i0 = first[0]; i0 <= last[0]; i0 += inc[0], k1++ ) {
                        size_t k2 = d_size.index( i0, i1, i2, i3, i4 );
                        d_data[k2] += subset.d_data[k1];
                    }
                }
            }
        }
    }
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::copySubset( const std::vector<size_t> &index,
                                              const Array<TYPE, FUN, Allocator> &subset )
{
    auto range = convert( index );
    copySubset( range, subset );
}

template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::addSubset( const std::vector<size_t> &index,
                                             const Array<TYPE, FUN, Allocator> &subset )
{
    auto range = convert( index );
    addSubset( range, subset );
}


/********************************************************
 *  Operator overloading                                 *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
bool Array<TYPE, FUN, Allocator>::operator==( const Array &rhs ) const
{
    if ( this == &rhs )
        return true;
    if ( d_size != rhs.d_size )
        return false;
    bool match = true;
    for ( size_t i = 0; i < d_size.length(); i++ )
        match = match && d_data[i] == rhs.d_data[i];
    return match;
}


/********************************************************
 *  Get a view of an C array                             *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
std::unique_ptr<Array<TYPE, FUN, Allocator>>
Array<TYPE, FUN, Allocator>::view( const ArraySize &N, std::shared_ptr<TYPE> data )
{
    auto array    = std::make_unique<Array<TYPE, FUN, Allocator>>();
    array->d_size = N;
    array->d_ptr  = data;
    array->d_data = array->d_ptr.get();
    return array;
}
template<class TYPE, class FUN, class Allocator>
std::unique_ptr<const Array<TYPE, FUN, Allocator>>
Array<TYPE, FUN, Allocator>::constView( const ArraySize &N,
                                        std::shared_ptr<const TYPE> const &data )
{
    auto array    = std::make_unique<Array<TYPE, FUN, Allocator>>();
    array->d_size = N;
    array->d_ptr  = std::const_pointer_cast<TYPE>( data );
    array->d_data = array->d_ptr.get();
    return array;
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::view2( Array<TYPE, FUN, Allocator> &src )
{
    view2( src.size(), src.getPtr() );
    d_data = src.d_data;
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::view2( const ArraySize &N, std::shared_ptr<TYPE> data )
{
    d_size = N;
    d_ptr  = data;
    d_data = d_ptr.get();
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::viewRaw( const ArraySize &N,
                                           TYPE *data,
                                           bool isCopyable,
                                           bool isFixedSize )
{
    d_isCopyable  = isCopyable;
    d_isFixedSize = isFixedSize;
    d_ptr.reset();
    d_size = N;
    d_data = data;
}


/********************************************************
 *  Basic functions                                      *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::swap( Array &other )
{
    // check that dimensions match
    if ( d_size != other.d_size )
        throw std::logic_error( "length of arrays do not match" );
    // swap the data
    std::swap( d_data, other.d_data );
    std::swap( d_ptr, other.d_ptr );
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::pow( const Array<TYPE, FUN, Allocator> &baseArray,
                                       const TYPE &exp )
{
    // not insisting on the shapes being the same
    // but insisting on the total size being the same
    if ( d_size.length() != baseArray.length() )
        throw std::logic_error( "length of arrays do not match" );

    const auto base_data = baseArray.data();
    for ( size_t i = 0; i < d_size.length(); i++ )
        d_data[i] = std::pow( base_data[i], exp );
}


/********************************************************
 *  Replicate the array                                  *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>
Array<TYPE, FUN, Allocator>::repmat( const std::vector<size_t> &N_rep ) const
{
    std::vector<size_t> N2( d_size.begin(), d_size.end() );
    if ( N2.size() < N_rep.size() )
        N2.resize( N_rep.size(), 1 );
    std::array<size_t, 5> N1, Nr;
    N1.fill( 1 );
    Nr.fill( 1 );
    for ( size_t d = 0; d < N_rep.size(); d++ ) {
        N1[d] = d_size[d];
        Nr[d] = N_rep[d];
        N2[d] *= N_rep[d];
    }
    Array<TYPE, FUN, Allocator> y( N2 );
    TYPE *y2 = y.data();
    for ( size_t i4 = 0, index = 0; i4 < N1[4]; i4++ ) {
        for ( size_t j4 = 0; j4 < Nr[4]; j4++ ) {
            for ( size_t i3 = 0; i3 < N1[3]; i3++ ) {
                for ( size_t j3 = 0; j3 < Nr[3]; j3++ ) {
                    for ( size_t i2 = 0; i2 < N1[2]; i2++ ) {
                        for ( size_t j2 = 0; j2 < Nr[2]; j2++ ) {
                            for ( size_t i1 = 0; i1 < N1[1]; i1++ ) {
                                for ( size_t j1 = 0; j1 < Nr[1]; j1++ ) {
                                    for ( size_t i0 = 0; i0 < N1[0]; i0++ ) {
                                        size_t k = d_size.index( i0, i1, i2, i3, i4 );
                                        TYPE x   = d_data[k];
                                        for ( size_t j0 = 0; j0 < Nr[0]; j0++, index++ )
                                            y2[index] = x;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    return y;
}


/********************************************************
 *  Simple math operations                               *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
bool Array<TYPE, FUN, Allocator>::NaNs() const
{
    bool test = false;
    for ( size_t i = 0; i < d_size.length(); i++ )
        test = test || d_data[i] != d_data[i];
    return test;
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::mean( void ) const
{
    TYPE x = sum() / d_size.length();
    return x;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::min( int dir ) const
{
    auto size_ans = d_size;
    size_ans.resize( dir, 1 );
    Array<TYPE, FUN, Allocator> ans( size_ans );
    size_t N1 = 1, N2 = 1, N3 = 1;
    for ( int d = 0; d < std::min<int>( dir, d_size.ndim() ); d++ )
        N1 *= d_size[d];
    N2 = d_size[dir];
    for ( size_t d = dir + 1; d < d_size.ndim(); d++ )
        N3 *= d_size[d];
    TYPE *data2 = ans.d_data;
    for ( size_t i3 = 0; i3 < N3; i3++ ) {
        for ( size_t i1 = 0; i1 < N1; i1++ ) {
            TYPE x = d_data[i1 + i3 * N1 * N2];
            for ( size_t i2 = 0; i2 < N2; i2++ )
                x = std::min( x, d_data[i1 + i2 * N1 + i3 * N1 * N2] );
            data2[i1 + i3 * N1] = x;
        }
    }
    return ans;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::max( int dir ) const
{
    auto size_ans = d_size;
    size_ans.resize( dir, 1 );
    Array<TYPE, FUN, Allocator> ans( size_ans );
    size_t N1 = 1, N2 = 1, N3 = 1;
    for ( int d = 0; d < std::min<int>( dir, d_size.ndim() ); d++ )
        N1 *= d_size[d];
    N2 = d_size[dir];
    DISABLE_WARNINGS // Suppress false array subscript is above array bounds
        for ( size_t d = dir + 1; d < d_size.ndim(); d++ ) N3 *= d_size[d];
    ENABLE_WARNINGS // Enable warnings
        TYPE *data2 = ans.d_data;
    for ( size_t i3 = 0; i3 < N3; i3++ ) {
        for ( size_t i1 = 0; i1 < N1; i1++ ) {
            TYPE x = d_data[i1 + i3 * N1 * N2];
            for ( size_t i2 = 0; i2 < N2; i2++ )
                x = std::max( x, d_data[i1 + i2 * N1 + i3 * N1 * N2] );
            data2[i1 + i3 * N1] = x;
        }
    }
    return ans;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::sum( int dir ) const
{
    auto size_ans = d_size;
    size_ans.resize( dir, 1 );
    Array<TYPE, FUN, Allocator> ans( size_ans );
    size_t N1 = 1, N2 = 1, N3 = 1;
    for ( int d = 0; d < std::min<int>( dir, d_size.ndim() ); d++ )
        N1 *= d_size[d];
    N2 = d_size[dir];
    DISABLE_WARNINGS
    for ( size_t d = dir + 1; d < d_size.ndim(); d++ )
        N3 *= d_size[d];
    ENABLE_WARNINGS
    TYPE *data2 = ans.d_data;
    for ( size_t i3 = 0; i3 < N3; i3++ ) {
        for ( size_t i1 = 0; i1 < N1; i1++ ) {
            TYPE x = 0;
            for ( size_t i2 = 0; i2 < N2; i2++ )
                x += d_data[i1 + i2 * N1 + i3 * N1 * N2];
            data2[i1 + i3 * N1] = x;
        }
    }
    return ans;
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::min( const std::vector<Range<size_t>> &range ) const
{
    // Get the subset indicies
    checkSubsetIndex( range );
    std::array<size_t, 5> first, last, inc, N1;
    getSubsetArrays( range, first, last, inc, N1 );
    TYPE x = std::numeric_limits<TYPE>::max();
    for ( size_t i4 = first[4]; i4 <= last[4]; i4 += inc[4] ) {
        for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) {
            for ( size_t i2 = first[2]; i2 <= last[2]; i2 += inc[2] ) {
                for ( size_t i1 = first[1]; i1 <= last[1]; i1 += inc[1] ) {
                    for ( size_t i0 = first[0]; i0 <= last[0]; i0 += inc[0] ) {
                        size_t k1 = d_size.index( i0, i1, i2, i3, i4 );
                        x         = std::min( x, d_data[k1] );
                    }
                }
            }
        }
    }
    return x;
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::max( const std::vector<Range<size_t>> &range ) const
{
    // Get the subset indicies
    checkSubsetIndex( range );
    std::array<size_t, 5> first, last, inc, N1;
    getSubsetArrays( range, first, last, inc, N1 );
    TYPE x = std::numeric_limits<TYPE>::min();
    for ( size_t i4 = first[4]; i4 <= last[4]; i4 += inc[4] ) {
        for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) {
            for ( size_t i2 = first[2]; i2 <= last[2]; i2 += inc[2] ) {
                for ( size_t i1 = first[1]; i1 <= last[1]; i1 += inc[1] ) {
                    for ( size_t i0 = first[0]; i0 <= last[0]; i0 += inc[0] ) {
                        size_t k1 = d_size.index( i0, i1, i2, i3, i4 );
                        x         = std::max( x, d_data[k1] );
                    }
                }
            }
        }
    }
    return x;
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::sum( const std::vector<Range<size_t>> &range ) const
{
    // Get the subset indicies
    checkSubsetIndex( range );
    std::array<size_t, 5> first, last, inc, N1;
    getSubsetArrays( range, first, last, inc, N1 );
    TYPE x = 0;
    for ( size_t i4 = first[4]; i4 <= last[4]; i4 += inc[4] ) {
        for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) {
            for ( size_t i2 = first[2]; i2 <= last[2]; i2 += inc[2] ) {
                for ( size_t i1 = first[1]; i1 <= last[1]; i1 += inc[1] ) {
                    for ( size_t i0 = first[0]; i0 <= last[0]; i0 += inc[0] ) {
                        size_t k1 = d_size.index( i0, i1, i2, i3, i4 );
                        x += d_data[k1];
                    }
                }
            }
        }
    }
    return x;
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::mean( const std::vector<Range<size_t>> &range ) const
{
    // Get the subset indicies
    checkSubsetIndex( range );
    std::array<size_t, 5> first, last, inc, N1;
    getSubsetArrays( range, first, last, inc, N1 );
    size_t n = 1;
    for ( auto &d : N1 )
        n *= d;
    TYPE x = sum( range ) / n;
    return x;
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::min( const std::vector<size_t> &index ) const
{
    auto range = convert( index );
    return min( range );
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::max( const std::vector<size_t> &index ) const
{
    auto range = convert( index );
    return max( range );
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::sum( const std::vector<size_t> &index ) const
{
    auto range = convert( index );
    return sum( range );
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::mean( const std::vector<size_t> &index ) const
{
    auto range = convert( index );
    return mean( range );
}


/********************************************************
 *  Find all elements that match the given operation     *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
std::vector<size_t>
Array<TYPE, FUN, Allocator>::find( const TYPE &value,
                                   std::function<bool( const TYPE &, const TYPE & )> compare ) const
{
    std::vector<size_t> result;
    result.reserve( d_size.length() );
    for ( size_t i = 0; i < d_size.length(); i++ ) {
        if ( compare( d_data[i], value ) )
            result.push_back( i );
    }
    return result;
}


/********************************************************
 *  Print an array to an output stream                   *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::print( std::ostream &os,
                                         const std::string &name,
                                         const std::string &prefix ) const
{
    if ( d_size.ndim() == 1 ) {
        for ( size_t i = 0; i < d_size[0]; i++ )
            os << prefix << name << "[" << i << "] = " << d_data[i] << std::endl;
    } else if ( d_size.ndim() == 2 ) {
        os << prefix << name << ":" << std::endl;
        for ( size_t i = 0; i < d_size[0]; i++ ) {
            for ( size_t j = 0; j < d_size[1]; j++ )
                os << prefix << "  " << operator()( i, j );
            os << std::endl;
        }
    } else {
        throw std::logic_error( "Not programmed for this dimension" );
    }
}


/********************************************************
 *  Reverse dimensions (transpose)                       *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::reverseDim() const
{
    size_t N2[5];
    for ( int d = 0; d < ArraySize::maxDim(); d++ )
        N2[d] = d_size[ArraySize::maxDim() - d - 1];
    ArraySize S2( ArraySize::maxDim(), N2 );
    Array<TYPE, FUN, Allocator> y( S2 );
    TYPE *y2 = y.data();
    for ( size_t i0 = 0; i0 < d_size[0]; i0++ ) {
        for ( size_t i1 = 0; i1 < d_size[1]; i1++ ) {
            for ( size_t i2 = 0; i2 < d_size[2]; i2++ ) {
                for ( size_t i3 = 0; i3 < d_size[3]; i3++ ) {
                    for ( size_t i4 = 0; i4 < d_size[4]; i4++ ) {
                        y2[S2.index( i4, i3, i2, i1, i0 )] =
                            d_data[d_size.index( i0, i1, i2, i3, i4 )];
                    }
                }
            }
        }
    }
    for ( int d = 0; d < d_size.ndim(); d++ )
        N2[d] = d_size[d_size.ndim() - d - 1];
    y.reshape( ArraySize( d_size.ndim(), N2 ) );
    return y;
}


/********************************************************
 *  Coarsen the array                                    *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>
Array<TYPE, FUN, Allocator>::coarsen( const Array<TYPE, FUN, Allocator> &filter ) const
{
    auto S2 = size();
    for ( size_t i = 0; i < S2.size(); i++ ) {
        size_t s = S2[i] / filter.size( i );
        S2.resize( i, s );
        if ( S2[i] * filter.size( i ) != size( i ) )
            throw std::invalid_argument( "Array must be multiple of filter size" );
    }
    Array<TYPE, FUN, Allocator> y( S2 );
    if ( d_size.ndim() > 3 )
        throw std::logic_error( "Function not programmed for more than 3 dimensions" );
    const auto &Nh = filter.d_size;
    for ( size_t k1 = 0; k1 < y.d_size[2]; k1++ ) {
        for ( size_t j1 = 0; j1 < y.d_size[1]; j1++ ) {
            for ( size_t i1 = 0; i1 < y.d_size[0]; i1++ ) {
                TYPE tmp = 0;
                for ( size_t k2 = 0; k2 < Nh[2]; k2++ ) {
                    for ( size_t j2 = 0; j2 < Nh[1]; j2++ ) {
                        for ( size_t i2 = 0; i2 < Nh[0]; i2++ ) {
                            tmp += filter( i2, j2, k2 ) * operator()( i1 *Nh[0] + i2,
                                                                      j1 * Nh[1] + j2,
                                                                      k1 * Nh[2] + k2 );
                        }
                    }
                }
                y( i1, j1, k1 ) = tmp;
            }
        }
    }
    return y;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::coarsen(
    const std::vector<size_t> &ratio,
    std::function<TYPE( const Array<TYPE, FUN, Allocator> & )> filter ) const
{
    if ( ratio.size() != d_size.ndim() )
        throw std::logic_error( "ratio size does not match ndim" );
    auto S2 = size();
    for ( size_t i = 0; i < S2.size(); i++ ) {
        S2.resize( i, S2[i] / ratio[i] );
        if ( S2[i] * ratio[i] != size( i ) )
            throw std::invalid_argument( "Array must be multiple of filter size" );
    }
    Array<TYPE, FUN, Allocator> tmp( ratio );
    Array<TYPE, FUN, Allocator> y( S2 );
    if ( d_size.ndim() > 3 )
        throw std::logic_error( "Function not programmed for more than 3 dimensions" );
    for ( size_t k1 = 0; k1 < y.d_size[2]; k1++ ) {
        for ( size_t j1 = 0; j1 < y.d_size[1]; j1++ ) {
            for ( size_t i1 = 0; i1 < y.d_size[0]; i1++ ) {
                for ( size_t k2 = 0; k2 < ratio[2]; k2++ ) {
                    for ( size_t j2 = 0; j2 < ratio[1]; j2++ ) {
                        for ( size_t i2 = 0; i2 < ratio[0]; i2++ ) {
                            tmp( i2, j2, k2 ) = operator()(
                                i1 *ratio[0] + i2, j1 * ratio[1] + j2, k1 * ratio[2] + k2 );
                        }
                    }
                }
                y( i1, j1, k1 ) = filter( tmp );
            }
        }
    }
    return y;
}


/********************************************************
 *  Concatenates the arrays                              *
 ********************************************************/
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::cat( const Array<TYPE, FUN, Allocator> &x, int dim )
{
    std::vector<Array<TYPE, FUN, Allocator>> tmp( 2 );
    tmp[0].view2( *this );
    tmp[1].view2( const_cast<Array<TYPE, FUN, Allocator> &>( x ) );
    *this = cat( tmp, dim );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::cat( const std::initializer_list<Array> &x,
                                                              int dim )
{
    return cat( x.size(), x.begin(), dim );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::cat( const std::vector<Array> &x, int dim )
{
    return cat( x.size(), x.data(), dim );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>
Array<TYPE, FUN, Allocator>::cat( size_t N_array, const Array *x, int dim )
{
    if ( N_array == 0 )
        return Array<TYPE, FUN, Allocator>();
    // Check that the dimensions match
    bool check = true;
    for ( size_t i = 1; i < N_array; i++ ) {
        check = check && x[i].ndim() == x[0].ndim();
        for ( int d = 0; d < x[0].ndim(); d++ )
            if ( d != dim )
                check = check && x[i].size( d ) == x[0].size( d );
    }
    if ( !check )
        throw std::logic_error( "Array dimensions do not match for concatenation" );
    // Create the output array
    auto size = x[0].d_size;
    for ( size_t i = 1; i < N_array; i++ )
        size.resize( dim, size[dim] + x[i].size( dim ) );
    Array<TYPE, FUN, Allocator> out( size );
    size_t N1 = 1;
    size_t N2 = size[dim];
    size_t N3 = 1;
    for ( int d = 0; d < dim; d++ )
        N1 *= size[d];
    for ( size_t d = dim + 1; d < size.ndim(); d++ )
        N3 *= size[d];
    TYPE *data = out.data();
    for ( size_t i = 0, i0 = 0; i < N_array; i++ ) {
        const TYPE *src = x[i].data();
        size_t N22      = x[i].size( dim );
        for ( size_t j2 = 0; j2 < N3; j2++ ) {
            for ( size_t i1 = 0; i1 < N22; i1++ ) {
                for ( size_t j1 = 0; j1 < N1; j1++ ) {
                    data[j1 + ( i1 + i0 ) * N1 + j2 * N1 * N2] = src[j1 + i1 * N1 + j2 * N1 * N22];
                }
            }
        }
        i0 += N22;
    }
    return out;
}


/********************************************************
 *  Interpolate                                          *
 ********************************************************/
template<class T>
constexpr bool is_compatible_double()
{
    return std::is_floating_point<T>::value || std::is_integral<T>::value;
}
template<class TYPE>
inline TYPE Array_interp_1D( double x, int N, const TYPE *data )
{
    if ( is_compatible_double<TYPE>() ) {
        int i = floor( x );
        i     = std::max( i, 0 );
        i     = std::min( i, N - 2 );
        return ( i + 1 - x ) * data[i] + ( x - i ) * data[i + 1];
    } else {
        throw std::logic_error( "Invalid conversion" );
    }
}
template<class TYPE>
inline TYPE Array_interp_2D( double x, double y, int Nx, int Ny, const TYPE *data )
{
    if ( is_compatible_double<TYPE>() ) {
        int i       = floor( x );
        i           = std::max( i, 0 );
        i           = std::min( i, Nx - 2 );
        double dx   = x - i;
        double dx2  = 1.0 - dx;
        int j       = floor( y );
        j           = std::max( j, 0 );
        j           = std::min( j, Ny - 2 );
        double dy   = y - j;
        double dy2  = 1.0 - dy;
        double f[4] = { (double) data[i + j * Nx],
                        (double) data[i + 1 + j * Nx],
                        (double) data[i + ( j + 1 ) * Nx],
                        (double) data[i + 1 + ( j + 1 ) * Nx] };
        return ( dx * f[1] + dx2 * f[0] ) * dy2 + ( dx * f[3] + dx2 * f[2] ) * dy;
    } else {
        throw std::logic_error( "Invalid conversion" );
    }
}
template<class TYPE>
inline TYPE
Array_interp_3D( double x, double y, double z, int Nx, int Ny, int Nz, const TYPE *data )
{
    if ( is_compatible_double<TYPE>() ) {
        int i       = floor( x );
        i           = std::max( i, 0 );
        i           = std::min( i, Nx - 2 );
        double dx   = x - i;
        double dx2  = 1.0 - dx;
        int j       = floor( y );
        j           = std::max( j, 0 );
        j           = std::min( j, Ny - 2 );
        double dy   = y - j;
        double dy2  = 1.0 - dy;
        int k       = floor( z );
        k           = std::max( k, 0 );
        k           = std::min( k, Nz - 2 );
        double dz   = z - k;
        double dz2  = 1.0 - dz;
        double f[8] = { (double) data[i + j * Nx + k * Nx * Ny],
                        (double) data[i + 1 + j * Nx + k * Nx * Ny],
                        (double) data[i + ( j + 1 ) * Nx + k * Nx * Ny],
                        (double) data[i + 1 + ( j + 1 ) * Nx + k * Nx * Ny],
                        (double) data[i + j * Nx + ( k + 1 ) * Nx * Ny],
                        (double) data[i + 1 + j * Nx + ( k + 1 ) * Nx * Ny],
                        (double) data[i + ( j + 1 ) * Nx + ( k + 1 ) * Nx * Ny],
                        (double) data[i + 1 + ( j + 1 ) * Nx + ( k + 1 ) * Nx * Ny] };
        double h0   = ( dx * f[1] + dx2 * f[0] ) * dy2 + ( dx * f[3] + dx2 * f[2] ) * dy;
        double h1   = ( dx * f[5] + dx2 * f[4] ) * dy2 + ( dx * f[7] + dx2 * f[6] ) * dy;
        return h0 * dz2 + h1 * dz;
    } else {
        throw std::logic_error( "Invalid conversion" );
    }
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::interp( const double *x ) const
{
    int ndim = 0, dim[5];
    double x2[5];
    for ( int d = 0; d < d_size.ndim(); d++ ) {
        if ( d_size[d] > 1 ) {
            x2[ndim]  = x[d];
            dim[ndim] = d_size[d];
            ndim++;
        }
    }
    TYPE f = 0;
    if ( ndim == 0 ) {
        // No data, do nothing
    } else if ( ndim == 1 ) {
        f = Array_interp_1D( x2[0], dim[0], d_data );
    } else if ( ndim == 2 ) {
        f = Array_interp_2D( x2[0], x2[1], dim[0], dim[1], d_data );
    } else if ( ndim == 3 ) {
        f = Array_interp_3D( x2[0], x2[1], x2[2], dim[0], dim[1], dim[2], d_data );
    } else {
        throw std::logic_error( "Not finished" );
    }
    return f;
}


/********************************************************
 *  Math operations (should call the Math class)         *
 ********************************************************/
/*template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::rand()
{
    FUN::rand( *this );
}
*/
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> &
Array<TYPE, FUN, Allocator>::operator+=( const Array<TYPE, FUN, Allocator> &rhs )
{
    auto op = []( const TYPE &a, const TYPE &b ) { return a + b; };
    FUN::transform( op, *this, rhs, *this );
    return *this;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> &
Array<TYPE, FUN, Allocator>::operator-=( const Array<TYPE, FUN, Allocator> &rhs )
{
    auto op = []( const TYPE &a, const TYPE &b ) { return a - b; };
    FUN::transform( op, *this, rhs, *this );
    return *this;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> &Array<TYPE, FUN, Allocator>::operator+=( const TYPE &rhs )
{
    auto op = [rhs]( const TYPE &x ) { return x + rhs; };
    FUN::transform( op, *this, *this );
    return *this;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> &Array<TYPE, FUN, Allocator>::operator-=( const TYPE &rhs )
{
    auto op = [rhs]( const TYPE &x ) { return x - rhs; };
    FUN::transform( op, *this, *this );
    return *this;
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::min() const
{
    const auto &op = []( const TYPE &a, const TYPE &b ) { return a < b ? a : b; };
    return FUN::reduce( op, *this, d_data[0] );
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::max() const
{
    const auto &op = []( const TYPE &a, const TYPE &b ) { return a > b ? a : b; };
    return FUN::reduce( op, *this, d_data[0] );
}
template<class TYPE, class FUN, class Allocator>
TYPE Array<TYPE, FUN, Allocator>::sum() const
{
    const auto &op = []( const TYPE &a, const TYPE &b ) { return a + b; };
    return FUN::reduce( op, *this, static_cast<TYPE>( 0 ) );
}
template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::axpby( const TYPE &alpha,
                                         const Array<TYPE, FUN, Allocator> &x,
                                         const TYPE &beta )
{
    const auto &op = [alpha, beta]( const TYPE &x, const TYPE &y ) { return alpha * x + beta * y; };
    return FUN::transform( op, x, *this, *this );
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>
Array<TYPE, FUN, Allocator>::transform( std::function<TYPE( const TYPE & )> fun,
                                        const Array<TYPE, FUN, Allocator> &x )
{
    Array<TYPE, FUN, Allocator> y;
    FUN::transform( fun, x, y );
    return y;
}
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator>
Array<TYPE, FUN, Allocator>::transform( std::function<TYPE( const TYPE &, const TYPE & )> fun,
                                        const Array<TYPE, FUN, Allocator> &x,
                                        const Array<TYPE, FUN, Allocator> &y )
{
    Array<TYPE, FUN, Allocator> z;
    FUN::transform( fun, x, y, z );
    return z;
}
template<class TYPE, class FUN, class Allocator>
bool Array<TYPE, FUN, Allocator>::equals( const Array &rhs, TYPE tol ) const
{
    return FUN::equals( *this, rhs, tol );
}


#endif