Array.h

Go to the documentation of this file.

/* Copyright 2017 Kristofer Björnson
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
 
#ifndef COM_DAFER45_TBTK_ARRAY
#define COM_DAFER45_TBTK_ARRAY
 
#include "TBTK/CArray.h"
#include "TBTK/Index.h"
#include "TBTK/MultiCounter.h"
#include "TBTK/Range.h"
#include "TBTK/Serializable.h"
#include "TBTK/TBTKMacros.h"
 
#include <vector>
 
namespace TBTK{
 
//Forward declation
namespace Math{
    template<typename DataType>
    class ArrayAlgorithms;
}; //End of namespace Math
 
template<typename DataType>
class Array : public Serializable{
public:
    //TBTKFeature Utilities.Array.construction.1 2019-10-31
    Array();
 
    //TBTKFeature Utilities.Array.construction.2 2019-10-31
    explicit Array(const std::initializer_list<unsigned int> &ranges);
 
    //TBTKFeature Utilities.Array.construction.3 2019-10-31
    Array(
        const std::initializer_list<unsigned int> &ranges,
        const DataType &fillValue
    );
 
    Array(const std::vector<DataType> &vector);
 
    Array(const std::vector<std::vector<DataType>> &vector);
 
    Array(const Range &range);
 
    Array(const std::string &serialization, Mode mode);
 
    static Array create(const std::vector<unsigned int> &ranges);
 
    static Array create(
        const std::vector<unsigned int> &ranges,
        const DataType &fillValue
    );
 
    template<typename CastType>
    operator Array<CastType>() const;
 
    //TBTKFeature Utilities.Array.operatorArraySubscript.1 2019-10-31
    DataType& operator[](const std::vector<unsigned int> &index);
 
    //TBTKFeature Utilities.Array.operatorArraySubscript.2 2019-10-31
    const DataType& operator[](
        const std::vector<unsigned int> &index
    ) const;
 
    //TBTKFeature Utilities.Array.operatorArraySubscript.3 2019-10-31
    DataType& operator[](unsigned int n);
 
    //TBTKFeature Utilities.Array.opreatorArraySubscript.4 2019-10-31
    const DataType& operator[](unsigned int n) const;
 
    //TBTKFeature Utilities.Array.operatorAddition.1 2019-10-31
    //TBTKFeature Utilities.Array.operatorAddition.2 2019-10-31
    //TBTKFeature Utilities.Array.operatorAddition.3 2019-10-31
    Array operator+(const Array &rhs) const;
 
    //TBTKFeature Utilities.Array.operatorSubtraction.1 2019-10-31
    //TBTKFeature Utilities.Array.operatorSubtraction.2 2019-10-31
    //TBTKFeature Utilities.Array.operatorSubtraction.3 2019-10-31
    Array operator-(const Array &rhs) const;
 
    //TBTKFeature Utilities.Array.operatorMultiplication.1 2019-10-31
    Array operator*(const DataType &rhs) const;
 
    //TBTKFeature Utilities.Array.operatorMultiplication.2 2019-10-31
    friend Array operator*(const DataType &lhs, const Array &rhs){
        Array<DataType> result = Array<DataType>::create(rhs.ranges);
        for(unsigned int n = 0; n < rhs.getSize(); n++)
            result.data[n] = lhs*rhs.data[n];
 
        return result;
    }
 
    Array operator*(const Array &rhs) const;
 
    //TBTKFeature Utilities.Array.operatorDivision.1 2019-10-31
    Array operator/(const DataType &rhs) const;
 
    bool operator==(const Array &rhs) const;
 
    static Array contract(
        const Array &array0,
        const std::vector<Subindex> &pattern0,
        const Array &array1,
        const std::vector<Subindex> &pattern1
    );
 
    //TBTKFeature Utilities.Array.getSlice.1
    Array<DataType> getSlice(const std::vector<Subindex> &index) const;
 
    Array<DataType> getArrayWithPermutedIndices(
        const std::vector<unsigned int> &permutation
    ) const;
 
    Array<DataType> getArrayWithReversedIndices() const;
 
    const std::vector<unsigned int>& getRanges() const;
 
    template<typename DT>
    friend std::ostream& operator<<(
        std::ostream &stream,
        const Array<DT> &array
    );
 
    //TBTKFeature Utilities.Array.getData.1 2019-10-31
    CArray<DataType>& getData();
 
    //TBTKFeature Utilities.Array.getData.2.C++ 2019-10-31
    const CArray<DataType>& getData() const;
 
    //TBTKFeature Utilities.Array.getSize.1 2019-10-31
    unsigned int getSize() const;
 
    std::string serialize(Mode mode) const;
private:
    CArray<DataType> data;
 
    std::vector<unsigned int> ranges;
 
    void fillSlice(
        Array &array,
        const std::vector<Subindex> &index,
        unsigned int subindex,
        unsigned int offsetSlice,
        unsigned int offsetOriginal
    ) const;
 
    void assertCompatibleRanges(
        const Array &array,
        std::string functionName
    ) const;
 
    friend class Math::ArrayAlgorithms<DataType>;
};
 
template<typename DataType>
Array<DataType>::Array(){
}
 
template<typename DataType>
Array<DataType>::Array(const std::initializer_list<unsigned int> &ranges){
    this->ranges = ranges;
    unsigned int size = 1;
    for(unsigned int n = 0; n < this->ranges.size(); n++){
        TBTKAssert(
            this->ranges[n] > 0,
            "Array::Array()",
            "Invalid ranges.",
            "'ranges' must only contain positive numbers."
        );
        size *= this->ranges[n];
    }
 
    data = CArray<DataType>(size);
}
 
template<typename DataType>
Array<DataType>::Array(
    const std::initializer_list<unsigned int> &ranges,
    const DataType &fillValue
){
    this->ranges = ranges;
    unsigned int size = 1;
    for(unsigned int n = 0; n < this->ranges.size(); n++){
        TBTKAssert(
            this->ranges[n] > 0,
            "Array::Array()",
            "Invalid ranges.",
            "'ranges' must only contain positive numbers."
        );
        size *= this->ranges[n];
    }
 
    data = CArray<DataType>(size);
 
    for(unsigned int n = 0; n < size; n++)
        data[n] = fillValue;
}
 
template<typename DataType>
Array<DataType>::Array(const std::vector<DataType> &vector){
    TBTKAssert(
        vector.size() != 0,
        "Array::Array()",
        "Invalid input.",
        "Unable to create an Array from an empty vector."
    );
 
    ranges.push_back(vector.size());
    data = CArray<DataType>(ranges[0]);
    for(unsigned int n = 0; n < ranges[0]; n++)
        data[n] = vector[n];
}
 
template<typename DataType>
Array<DataType>::Array(const std::vector<std::vector<DataType>> &vector){
    TBTKAssert(
        vector.size() != 0,
        "Array::Array()",
        "Invalid input. Unable to create an Array from an empty"
        << " vector.",
        ""
    );
    TBTKAssert(
        vector[0].size() != 0,
        "Array::Array()",
        "Invalid input. Unable to create an Array from a vector with"
        << " an empty row.",
        ""
    );
    ranges.push_back(vector.size());
    ranges.push_back(vector[0].size());
 
    for(unsigned int n = 1; n < vector.size(); n++){
        TBTKAssert(
            vector[n].size() == vector[0].size(),
            "Array::Array()",
            "Invalid input. vector[" << n << "] has a different"
            << " size than vector[0]",
            ""
        );
    }
 
    data = CArray<DataType>(ranges[0]*ranges[1]);
    for(unsigned int x = 0; x < ranges[0]; x++)
        for(unsigned int y = 0; y < ranges[1]; y++)
            data[ranges[1]*x + y] = vector[x][y];
}
 
template<typename DataType>
Array<DataType>::Array(const Range &range){
    ranges.push_back(range.getResolution());
    data = CArray<DataType>(range.getResolution());
    for(unsigned int n = 0; n < range.getResolution(); n++)
        data[n] = (DataType)range[n];
}
 
template<typename DataType>
Array<DataType>::Array(const std::string &serialization, Mode mode){
    TBTKAssert(
        validate(serialization, "Array", mode),
        "Array::Array()",
        "Unable to parse string as Array '" << serialization << "'.",
        ""
    );
 
    switch(mode){
    case Mode::JSON:
    {
        try{
            nlohmann::json j
                = nlohmann::json::parse(serialization);
            data = CArray<DataType>(j.at("data"), mode);
            ranges = j.at("ranges").get<std::vector<unsigned int>>();
        }
        catch(nlohmann::json::exception &e){
            TBTKExit(
                "Array::Array()",
                "Unable to parse string as Array '"
                << serialization << "'.",
                ""
            );
        }
 
        break;
    }
    default:
        TBTKExit(
            "Array::Array()",
            "Unable to parse string as Array '" << serialization
            << "'.",
            ""
        );
    }
}
 
template<typename DataType>
Array<DataType> Array<DataType>::create(
    const std::vector<unsigned int> &ranges
){
    Array<DataType> array;
    array.ranges = ranges;
    unsigned int size = 1;
    for(unsigned int n = 0; n < array.ranges.size(); n++){
        TBTKAssert(
            array.ranges[n] > 0,
            "Array::Array()",
            "Invalid ranges.",
            "'ranges' must only contain positive numbers."
        );
        size *= array.ranges[n];
    }
 
    array.data = CArray<DataType>(size);
 
    return array;
}
 
template<typename DataType>
Array<DataType> Array<DataType>::create(
    const std::vector<unsigned int> &ranges,
    const DataType &fillValue
){
    Array<DataType> array;
    array.ranges = ranges;
    unsigned int size = 1;
    for(unsigned int n = 0; n < array.ranges.size(); n++){
        TBTKAssert(
            array.ranges[n] > 0,
            "Array::Array()",
            "Invalid ranges.",
            "'ranges' must only contain positive numbers."
        );
        size *= array.ranges[n];
    }
 
    array.data = CArray<DataType>(size);
 
    for(unsigned int n = 0; n < size; n++)
        array.data[n] = fillValue;
 
    return array;
}
 
template<typename DataType>
template<typename CastType>
inline Array<DataType>::operator Array<CastType>() const{
    Array<CastType> result = Array<CastType>::create(ranges);
    CArray<CastType> &resultData = result.getData();
    for(unsigned int n = 0; n < data.getSize(); n++)
        resultData[n] = (CastType)data[n];
 
    return result;
}
 
template<typename DataType>
inline DataType& Array<DataType>::operator[](
    const std::vector<unsigned int> &index
){
    unsigned int idx = 0;
    for(unsigned int n = 0; n < index.size(); n++){
        if(n != 0)
            idx *= ranges[n];
        idx += index[n];
    }
 
    return data[idx];
}
 
template<typename DataType>
inline const DataType& Array<DataType>::operator[](
    const std::vector<unsigned int> &index
) const{
    unsigned int idx = 0;
    for(unsigned int n = 0; n < index.size(); n++){
        if(n != 0)
            idx *= ranges[n];
        idx += index[n];
    }
 
    return data[idx];
}
 
template<typename DataType>
inline DataType& Array<DataType>::operator[](unsigned int n){
    return data[n];
}
 
template<typename DataType>
inline const DataType& Array<DataType>::operator[](unsigned int n) const{
    return data[n];
}
 
template<typename DataType>
inline Array<DataType> Array<DataType>::operator+(
    const Array<DataType> &rhs
) const{
    assertCompatibleRanges(rhs, "operator+()");
 
    Array<DataType> result = Array<DataType>::create(ranges);
    for(unsigned int n = 0; n < data.getSize(); n++)
        result.data[n] = data[n] + rhs.data[n];
 
    return result;
}
 
template<typename DataType>
inline Array<DataType> Array<DataType>::operator-(
    const Array<DataType> &rhs
) const{
    assertCompatibleRanges(rhs, "operator+()");
 
    Array<DataType> result = Array<DataType>::create(ranges);
    for(unsigned int n = 0; n < data.getSize(); n++)
        result.data[n] = data[n] - rhs.data[n];
 
    return result;
}
 
template<typename DataType>
inline Array<DataType> Array<DataType>::operator*(
    const DataType &rhs
) const{
    Array<DataType> result = Array<DataType>::create(ranges);
    for(unsigned int n = 0; n < data.getSize(); n++)
        result.data[n] = data[n]*rhs;
 
    return result;
}
 
template<typename DataType>
Array<DataType> Array<DataType>::operator*(const Array<DataType> &rhs) const{
    switch(ranges.size()){
    case 1:
        switch(rhs.ranges.size()){
        case 1:
            return contract(
                *this,
                {IDX_ALL_(0)},
                rhs,
                {IDX_ALL_(0)}
            );
        case 2:
            return contract(
                *this,
                {IDX_ALL_(0)},
                rhs,
                {IDX_ALL_(0), IDX_ALL}
            );
        default:
            TBTKExit(
                "Array::operator*()",
                "Unsupported Array rank. Multiplication is"
                << " only possible for Arrays of rank one or"
                << " two, but the left hand side has rank '"
                << ranges.size() << "'.",
                ""
            );
        }
    case 2:
        switch(rhs.ranges.size()){
        case 1:
            return contract(
                *this,
                {IDX_ALL, IDX_ALL_(0)},
                rhs,
                {IDX_ALL_(0)}
            );
        case 2:
            return contract(
                *this,
                {IDX_ALL, IDX_ALL_(0)},
                rhs,
                {IDX_ALL_(0), IDX_ALL}
            );
        default:
            TBTKExit(
                "Array::operator*()",
                "Unsupported Array rank. Multiplication is"
                << " only possible for Arrays of rank one or"
                << " two, but the left hand side has rank '"
                << ranges.size() << "'.",
                ""
            );
        }
    default:
        TBTKExit(
            "Array::operator*()",
            "Unsupported Array rank. Multiplication is only"
            << " possible for Arrays of rank one or two, but the"
            << " left hand side has rank '" << ranges.size()
            << "'.",
            ""
        );
    }
}
 
template<typename DataType>
inline Array<DataType> Array<DataType>::operator/(const DataType &rhs) const{
    Array<DataType> result = Array<DataType>::create(ranges);
    for(unsigned int n = 0; n < data.getSize(); n++)
        result.data[n] = data[n]/rhs;
 
    return result;
}
 
template<typename DataType>
inline bool Array<DataType>::operator==(const Array<DataType> &rhs) const{
    if(ranges.size() != rhs.ranges.size())
        return false;
    for(unsigned int n = 0; n < ranges.size(); n++)
        if(ranges[n] != rhs.ranges[n])
            return false;
 
    for(unsigned int n = 0; n < getSize(); n++)
        if(data[n] != rhs.data[n])
            return false;
 
    return true;
}
 
template<typename DataType>
Array<DataType> Array<DataType>::contract(
    const Array &array0,
    const std::vector<Subindex> &pattern0,
    const Array &array1,
    const std::vector<Subindex> &pattern1
){
    const std::vector<unsigned int> &ranges0 = array0.getRanges();
    const std::vector<unsigned int> &ranges1 = array1.getRanges();
    TBTKAssert(
        ranges0.size() == pattern0.size(),
        "Array::contract()",
        "Incompatible pattern size. The number of elements in"
        << " 'pattern0' must be the same as the number of ranges in"
        << " 'array0', but 'pattern0' has '" << pattern0.size() << "'"
        << " elements while 'array0' has '"
        << ranges0.size() << "' ranges.",
        ""
    );
    TBTKAssert(
        ranges1.size() == pattern1.size(),
        "Array::contract()",
        "Incompatible pattern size. The number of elements in"
        << " 'pattern1' must be the same as the number of ranges in"
        << " 'array1', but 'pattern1' has '" << pattern1.size() << "'"
        << " elements while 'array1' has '"
        << ranges1.size() << "' ranges.",
        ""
    );
    for(unsigned int n = 0; n < pattern0.size(); n++){
        TBTKAssert(
            pattern0[n].isWildcard()
            || pattern0[n].isLabeledWildcard(),
            "Array::contract()",
            "Invalid pattern. The patterns must only contain"
            " wildcards or labeled wildcards, but found '"
            << pattern0[n] << "' in position '" << n << "' of"
            << " 'pattern0'.",
            ""
        );
    }
    for(unsigned int n = 0; n < pattern1.size(); n++){
        TBTKAssert(
            pattern1[n].isWildcard()
            || pattern1[n].isLabeledWildcard(),
            "Array::contract()",
            "Invalid pattern. The patterns must only contain"
            " wildcards or labeled wildcards, but found '"
            << pattern1[n] << "' in position '" << n << "' of"
            << " 'pattern1'.",
            ""
        );
    }
 
    std::vector<unsigned int> summationIndices0;
    std::vector<Subindex> wildcards0;
    for(unsigned int n = 0; n < pattern0.size(); n++){
        if(pattern0[n].isLabeledWildcard()){
            summationIndices0.push_back(n);
            wildcards0.push_back(pattern0[n]);
        }
    }
    std::vector<unsigned int> summationIndices1;
    std::vector<Subindex> wildcards1;
    for(unsigned int n = 0; n < pattern1.size(); n++){
        if(pattern1[n].isLabeledWildcard()){
            summationIndices1.push_back(n);
            wildcards1.push_back(pattern1[n]);
        }
    }
 
    TBTKAssert(
        wildcards0.size() == wildcards1.size(),
        "Array::contract()",
        "Incompatible patterns. The number of labeled wildcards are"
        << " different in 'pattern0' and 'pattern1'.",
        ""
    );
    for(unsigned int n = 0; n < wildcards0.size(); n++){
        for(unsigned int c = n+1; c < wildcards0.size(); c++){
            TBTKAssert(
                wildcards0[n] != wildcards0[c],
                "Array::contract()",
                "Repeated labeled wildcard in 'pattern0'",
                ""
            );
        }
    }
    for(unsigned int n = 0; n < wildcards1.size(); n++){
        for(unsigned int c = n+1; c < wildcards1.size(); c++){
            TBTKAssert(
                wildcards1[n] != wildcards1[c],
                "Array::contract()",
                "Repeated labeled wildcard in 'pattern1'",
                ""
            );
        }
    }
 
    std::vector<unsigned int> summationIndicesMap;
    for(unsigned int n = 0; n < wildcards0.size(); n++){
        for(unsigned int c = 0; c < wildcards1.size(); c++){
            if(wildcards0[n] == wildcards1[c]){
                summationIndicesMap.push_back(
                    summationIndices1[c]
                );
                break;
            }
        }
        TBTKAssert(
            summationIndicesMap.size() == n + 1,
            "Array::contract()",
            "Incompatible patterns. The labeled wildcard at"
            << " position '" << summationIndices0[n] << " in"
            << " 'pattern0' is missing in 'pattern1'.",
            ""
        );
    }
 
    std::vector<unsigned int> summationRanges;
    for(unsigned int n = 0; n < summationIndices0.size(); n++){
        TBTKAssert(
            ranges0[summationIndices0[n]]
                == ranges1[summationIndicesMap[n]],
            "Array::contract()",
            "Incompatible summation indices. Unable to contract"
            << " the Subindex at position '"
            << summationIndices0[n] << "' in 'array0' with the"
            << " Subindex at position '" << summationIndicesMap[n]
            << "' in 'array1' since they have different range.",
            ""
        );
        summationRanges.push_back(ranges0[summationIndices0[n]]);
    }
 
    std::vector<unsigned int> resultRanges;
    for(unsigned int n = 0; n < ranges0.size(); n++)
        if(pattern0[n].isWildcard())
            resultRanges.push_back(ranges0[n]);
    for(unsigned int n = 0; n < ranges1.size(); n++)
        if(pattern1[n].isWildcard())
            resultRanges.push_back(ranges1[n]);
 
    Array result;
    if(resultRanges.size() == 0){
        result = Array({1}, 0);
 
        std::vector<unsigned int> summationInitialValues(
            summationRanges.size(),
            0
        );
        std::vector<unsigned int> summationIncrements(
            summationRanges.size(),
            1
        );
        MultiCounter<unsigned int> summationCounter(
            summationInitialValues,
            summationRanges,
            summationIncrements
        );
        for(
            summationCounter.reset();
            !summationCounter.done();
            ++summationCounter
        ){
            std::vector<unsigned int> index0(
                summationCounter.getSize()
            );
            std::vector<unsigned int> index1(
                summationCounter.getSize()
            );
            for(
                unsigned int n = 0;
                n < summationCounter.getSize();
                n++
            ){
                index0[summationIndices0[n]]
                    = summationCounter[n];
                index1[summationIndicesMap[n]]
                    = summationCounter[n];
            }
            result[{0}] += array0[index0]*array1[index1];
        }
    }
    else{
        result = Array::create(resultRanges, 0);
 
        std::vector<unsigned int> resultInitialValues(
            resultRanges.size(),
            0
        );
        std::vector<unsigned int> resultIncrements(
            resultRanges.size(),
            1
        );
        MultiCounter<unsigned int> resultCounter(
            resultInitialValues,
            resultRanges,
            resultIncrements
        );
 
        std::vector<unsigned int> summationInitialValues(
        summationRanges.size(),
        0
        );
        std::vector<unsigned int> summationIncrements(
            summationRanges.size(),
            1
        );
        MultiCounter<unsigned int> summationCounter(
            summationInitialValues,
            summationRanges,
            summationIncrements
        );
        for(resultCounter.reset(); !resultCounter.done(); ++resultCounter){
            std::vector<unsigned int> resultIndex
                = (std::vector<unsigned int>)resultCounter;
            for(
                summationCounter.reset();
                !summationCounter.done();
                ++summationCounter
            ){
                std::vector<unsigned int> index0(
                    resultIndex.begin(),
                    resultIndex.begin() + ranges0.size()
                        - summationIndices0.size()
                );
                std::vector<unsigned int> index1(
                    resultIndex.begin() + ranges0.size()
                        - summationIndices0.size(),
                    resultIndex.end()
                );
                for(
                    unsigned int n = 0;
                    n < summationCounter.getSize();
                    n++
                ){
                    index1.insert(
                        index1.begin() + summationIndices1[n],
                        0
                    );
                }
                for(
                    unsigned int n = 0;
                    n < summationCounter.getSize();
                    n++
                ){
                    index0.insert(
                        index0.begin() + summationIndices0[n],
                        summationCounter[n]
                    );
                    index1[summationIndicesMap[n]]
                        = summationCounter[n];
                }
                result[resultCounter] += array0[index0]*array1[index1];
            }
        }
    }
 
    return result;
}
 
template<typename DataType>
Array<DataType> Array<DataType>::getSlice(const std::vector<Subindex> &index) const{
    TBTKAssert(
        ranges.size() == index.size(),
        "Array::getSlice()",
        "Incompatible ranges.",
        "'index' must have the same number of dimensions as 'ranges'."
    );
 
    std::vector<unsigned int> newRanges;
    for(unsigned int n = 0; n < ranges.size(); n++){
        TBTKAssert(
            index[n] < (int)ranges[n],
            "Array::getSlice()",
            "'index' out of range.",
            ""
        );
        if(index[n] < 0){
            TBTKAssert(
                index[n].isWildcard(),
                "Array::getSlice()",
                "Invalid symbol.",
                "'index' can only contain positive numbers or"
                << " 'IDX_ALL'."
            );
            newRanges.push_back(ranges[n]);
        }
    }
 
    Array array = Array::create(newRanges);
 
    fillSlice(array, index, 0, 0, 0);
 
    return array;
}
 
template<typename DataType>
Array<DataType> Array<DataType>::getArrayWithPermutedIndices(
    const std::vector<unsigned int> &permutation
) const{
    TBTKAssert(
        permutation.size() == ranges.size(),
        "Array::getArrayWithPermutedIndices()",
        "The number of permutation indices '" << permutation.size()
        << "' must be the same a the number of ranges '"
        << ranges.size() << "'.",
        ""
    );
 
    std::vector<bool> indexIncluded(permutation.size(), false);
    for(unsigned int n = 0; n < permutation.size(); n++){
        TBTKAssert(
            permutation[n] >= 0
            && permutation[n] < permutation.size(),
            "Array::getArrayWithPermutedIndices()",
            "Invalid permutation values 'permutation[" << n << "]"
            << " = " << permutation[n] << "'. Must be a number"
            << " between 0 and N-1, where N is the number of"
            << " ranges.",
            ""
        );
        indexIncluded[permutation[n]] = true;
    }
    for(unsigned int n = 0; n < indexIncluded.size(); n++){
        TBTKAssert(
            indexIncluded[n],
            "Array::getArrayWithPermutedIndices()",
            "Invalid permutation. Missing permutation index '" << n
            << "'.",
            ""
        );
    }
 
    std::vector<unsigned int> newRanges;
    for(unsigned int n = 0; n < ranges.size(); n++)
        newRanges.push_back(ranges[permutation[n]]);
    Array<DataType> newArray = Array<DataType>::create(newRanges);
 
    std::vector<unsigned int> begin = newRanges;
    std::vector<unsigned int> end = newRanges;
    std::vector<unsigned int> increment = newRanges;
    for(unsigned int n = 0; n < newRanges.size(); n++){
        begin[n] = 0;
        increment[n] = 1;
    }
    MultiCounter<unsigned int> counter(begin, end, increment);
    for(counter.reset(); !counter.done(); ++counter){
        std::vector<unsigned int> newArrayIndex = counter;
        std::vector<unsigned int> arrayIndex(newArrayIndex.size());
        for(unsigned int c = 0; c < newArrayIndex.size(); c++)
            arrayIndex[permutation[c]] = newArrayIndex[c];
        newArray[newArrayIndex] = (*this)[arrayIndex];
    }
 
    return newArray;
}
 
template<typename DataType>
Array<DataType> Array<DataType>::getArrayWithReversedIndices() const{
    std::vector<unsigned int> permutation;
    for(unsigned int n = 0; n < ranges.size(); n++)
        permutation.push_back(ranges.size() - n - 1);
 
    return getArrayWithPermutedIndices(permutation);
}
 
template<typename DataType>
void Array<DataType>::fillSlice(
    Array &array,
    const std::vector<Subindex> &index,
    unsigned int subindex,
    unsigned int offsetSlice,
    unsigned int offsetOriginal
) const{
    if(subindex == index.size()-1){
        if(index[subindex].isWildcard()){
            for(unsigned int n = 0; n < ranges[subindex]; n++){
                array.data[offsetSlice*ranges[subindex] + n]
                    = data[
                        offsetOriginal*ranges[subindex]
                        + n
                    ];
            }
        }
        else{
            array.data[offsetSlice] = data[
                offsetOriginal*ranges[subindex]
                + index[subindex]
            ];
        }
    }
    else{
        if(index[subindex].isWildcard()){
            for(unsigned int n = 0; n < ranges[subindex]; n++){
                fillSlice(
                    array,
                    index,
                    subindex+1,
                    offsetSlice*ranges[subindex] + n,
                    offsetOriginal*ranges[subindex] + n
                );
            }
        }
        else{
            fillSlice(
                array,
                index,
                subindex+1,
                offsetSlice,
                offsetOriginal*ranges[subindex]
                    + index[subindex]
            );
        }
    }
}
 
template<typename DataType>
inline const std::vector<unsigned int>& Array<DataType>::getRanges() const{
    return ranges;
}
 
template<typename DataType>
inline std::ostream& operator<<(
    std::ostream &stream,
    const Array<DataType> &array
){
    switch(array.ranges.size()){
    case 1:
        stream << "[";
        for(unsigned int n = 0; n < array.ranges[0]; n++){
            if(n != 0)
                stream << ", ";
            stream << array[{n}];
        }
        stream << "]";
 
        break;
    case 2:
        stream << "[";
        for(unsigned int row = 0; row < array.ranges[0]; row++){
            if(row != 0)
                stream << "\n";
            stream << "[";
            for(
                unsigned int column = 0;
                column < array.ranges[1];
                column++
            ){
                if(column != 0)
                    stream << ", ";
                stream << array[{row, column}];
            }
            stream << "]";
        }
        stream << "]";
 
        break;
    default:
        TBTKExit(
            "Array::operator<<()",
            "Unable to print Array of rank '"
            << array.ranges.size() << "'.",
            ""
        );
    }
 
    return stream;
}
 
template<typename DataType>
inline CArray<DataType>& Array<DataType>::getData(){
    return data;
}
 
template<typename DataType>
inline const CArray<DataType>& Array<DataType>::getData() const{
    return data;
}
 
template<typename DataType>
inline unsigned int Array<DataType>::getSize() const{
    return data.getSize();
}
 
template<typename DataType>
inline std::string Array<DataType>::serialize(Mode mode) const{
    switch(mode){
    case Mode::JSON:
    {
        nlohmann::json j;
        j["id"] = "Array";
        j["data"] = data.serialize(mode);
        j["ranges"] = ranges;
 
        return j.dump();
    }
    default:
        TBTKExit(
            "Array::serialize()",
            "Only Serializable::Mode::JSON is supported yet.",
            ""
        );
    }
}
 
template<typename DataType>
inline void Array<DataType>::assertCompatibleRanges(
    const Array<DataType> &array,
    std::string functionName
) const{
    TBTKAssert(
        ranges.size() == array.ranges.size(),
        "Array::" + functionName,
        "Incompatible ranges.",
        "Left and right hand sides must have the same number of"
        << " dimensions."
    );
    for(unsigned int n = 0; n < ranges.size(); n++){
        TBTKAssert(
            ranges[n] == array.ranges[n],
            "Array::" + functionName,
            "Incompatible ranges.",
            "Left and right hand sides must have the same ranges."
        );
    }
}
 
}; //End of namesapce TBTK
 
#endif