/*------------------------------------------------------------------------------
 *	Copyright (c) 2023 by Bai Bing (seread@163.com)
 *	See COPYING file for copying and redistribution conditions.
 *
 *  Alians IT Studio.
 *----------------------------------------------------------------------------*/
#pragma once

#include <algorithm>
#include <array>
#include <cstring>
#include <execution>
#include <iterator>
#include <list>
#include <numeric>
#include <random>
#include <stdexcept>
#include <utility>
#include <vector>

#include "_Define.h"
#include "core/ComplexOperators.h"
#include "core/Endian.h"
#include "core/Enums.h"
#include "core/Error.h"
#include "core/FileSystem.h"
#include "core/StaticAsserts.h"
#include "core/String.h"
#include "core/TypeTraits.h"
#include "ASIterators.h"
#include "ASShape.h"
#include "ASSlice.h"
#include "core/FunCheck.h"

namespace ais
{
    template <typename dtype = double, class Allocator = std::allocator<dtype>>
    class AIS_EXPORT Matrix
    {
    private:
        STATIC_ASSERT_VALID_DTYPE(dtype);
        static_assert(is_same_v<dtype, typename Allocator::value_type>,
                      "value_type and Allocator::value_type must match");

        using AllocType = typename std::allocator_traits<Allocator>::template rebind_alloc<dtype>;
        using AllocTraits = std::allocator_traits<AllocType>;
        using self_type = Matrix<dtype, AllocType>;

    public:
        using value_type = dtype;
        using allocator_type = Allocator;
        using pointer = typename AllocTraits::pointer;
        using const_pointer = typename AllocTraits::const_pointer;
        using reference = dtype &;
        using const_reference = const dtype &;
        using difference_type = typename AllocTraits::difference_type;

        using iterator = Iterator<dtype, pointer, difference_type>;
        using const_iterator = ConstIterator<dtype, const_pointer, difference_type>;
        using reverse_iterator = std::reverse_iterator<iterator>;
        using const_reverse_iterator = std::reverse_iterator<const_iterator>;

        using column_iterator = ColumnIterator<dtype, size_t, pointer, difference_type>;
        using const_column_iterator = ConstColumnIterator<dtype, size_t, const_pointer, difference_type>;
        using reverse_column_iterator = std::reverse_iterator<column_iterator>;
        using const_reverse_column_iterator = std::reverse_iterator<const_column_iterator>;

        //============================================================================
        // Method Description:
        /// Default Constructor, not very useful...
        ///
        Matrix() = default;

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param squareSize: square number of rows and columns
        ///
        explicit Matrix(size_t squareSize) : shape_(squareSize, squareSize),
                                             size_(squareSize * squareSize)
        {
            new_matrix();
        }

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param rows
        /// @param columns
        ///
        Matrix(size_t rows, size_t columns) : shape_(rows, columns),
                                              size_(rows * columns)
        {
            new_matrix();
        }

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param shape
        ///
        explicit Matrix(const Shape &shape) : shape_(shape),
                                              size_(shape_.size())
        {
            new_matrix();
        }

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param list
        ///
        Matrix(std::initializer_list<dtype> list) : shape_(1, (list.size())),
                                                    size_(shape_.size())
        {
            new_matrix();
            // copy data from shape
            if (size_ > 0)
            {
                std::copy(std::execution::par_unseq, list.begin(), list.end(), begin());
            }
        }

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param list: 2D initializer list
        ///
        Matrix(const std::initializer_list<std::initializer_list<dtype>> &list) : shape_((list.size()), 0)
        {
            for (const auto &li : list)
            {
                if (shape_.cols == 0)
                {
                    shape_.cols = (li.size());
                }
                else if (li.size() != shape_.cols)
                {
                    THROW_INVALID_ARGUMENT("All rows of the initializer list needs to have the same number of elements");
                }
            }

            size_ = shape_.size();
            new_matrix();
            size_t row = 0;
            for (const auto &li : list)
            {
                const auto ptr = begin() += row * shape_.cols;
                std::copy(std::execution::par_unseq, li.begin(), li.end(), ptr);
                ++row;
            }
        }

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param array
        /// @param copy: (optional) boolean for whether to make a copy and own the data, or
        /// act as a non-owning shell. Default true.
        ///
        template <size_t arraySize, std::enable_if_t<is_valid_dtype_v<dtype>, int> = 0>
        Matrix(std::array<dtype, arraySize> &array, bool copy = true) : shape_(1, arraySize),
                                                                        size_(shape_.size())
        {
            if (copy)
            {
                new_matrix();
                if (size_ > 0)
                {
                    std::copy(std::execution::par_unseq, array.begin(), array.end(), begin());
                }
            }
            else
            {
                array_ = array.data();
                ownsPtr_ = false;
            }
        }

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param array2D
        /// @param copy: (optional) boolean for whether to make a copy and own the data, or
        /// act as a non-owning shell. Default true.
        ///
        template <size_t d0Size, size_t d1Size>
        Matrix(std::array<std::array<dtype, d1Size>, d0Size> &array2D, bool copy = true) : shape_(d0Size, d1Size),
                                                                                           size_(shape_.size())
        {
            if (copy)
            {
                new_matrix();
                if (size_ > 0)
                {
                    const auto start = array2D.front().begin();
                    std::copy(std::execution::par_unseq, start, start + size_, begin());
                }
            }
            else
            {
                array_ = array2D.front().data();
                ownsPtr_ = false;
            }
        }

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param vector
        /// @param copy: (optional) boolean for whether to make a copy and own the data, or
        /// act as a non-owning shell. Default true.
        ///
        template <std::enable_if_t<is_valid_dtype_v<dtype>, int> = 0>
        Matrix(std::vector<dtype> &vector, bool copy = true) : shape_(1, vector.size()),
                                                               size_(shape_.size())
        {
            if (copy)
            {
                new_matrix();
                if (size_ > 0)
                {
                    std::copy(std::execution::par_unseq, vector.begin(), vector.end(), begin());
                }
            }
            else
            {
                array_ = vector.data();
                ownsPtr_ = false;
            }
        }

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param vector2D
        ///
        explicit Matrix(const std::vector<std::vector<dtype>> &vector2D) : shape_(vector2D.size(), 0)
        {
            for (const auto &row : vector2D)
            {
                if (shape_.cols == 0)
                {
                    shape_.cols = row.size();
                }
                else if (row.size() != shape_.cols)
                {
                    THROW_INVALID_ARGUMENT("All rows of the 2d vector need to have the same number of elements");
                }
            }

            size_ = shape_.size();

            new_matrix();
            auto currentPosition = begin();
            for (const auto &row : vector2D)
            {
                std::copy(std::execution::par_unseq, row.begin(), row.end(), currentPosition);
                currentPosition += shape_.cols;
            }
        }

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param list
        ///
        explicit Matrix(const std::list<dtype> &list) : shape_(1, list.size()),
                                                        size_(shape_.size())
        {
            new_matrix();
            if (size_ > 0)
            {
                std::copy(std::execution::par_unseq, list.begin(), list.end(), begin());
            }
        }

        //============================================================================
        // Method Description:
        /// Constructor
        ///
        /// @param begin
        /// @param end
        ///
        template <typename Iterator,
                  std::enable_if_t<is_same_v<typename std::iterator_traits<Iterator>::value_type, dtype>, int> = 0>
        Matrix(Iterator begin, Iterator end) : shape_(1, (std::distance(begin, end))),
                                               size_(shape_.size())
        {
            new_matrix();
            if (size_ > 0)
            {
                std::copy(std::execution::par_unseq, begin, end, begin());
            }
        }

        //============================================================================
        // Method Description:
        /// Constructor.  Copies the contents of the buffer into
        /// the array.
        ///
        /// @param ptr: const_pointer to beginning of buffer
        /// @param size: number of elements in buffer
        ///
        Matrix(const_pointer ptr, size_t size) : shape_(1, size),
                                                 size_(size)
        {
            new_matrix();
            if (ptr != nullptr && size_ > 0)
            {
                std::copy(std::execution::par_unseq, ptr, ptr + size_, begin());
            }
        }

        //============================================================================
        // Method Description:
        /// Constructor.  Copies the contents of the buffer into
        /// the array.
        ///
        /// @param ptr: const_pointer to beginning of buffer
        /// @param rows: number of rows of the buffer
        /// @param columns: number of cols of the buffer
        ///
        template <typename UIntType1,
                  typename UIntType2,
                  std::enable_if_t<!is_same_v<UIntType1, bool>, int> = 0,
                  std::enable_if_t<!is_same_v<UIntType2, bool>, int> = 0>
        Matrix(const_pointer ptr, UIntType1 rows, UIntType2 columns) : shape_(rows, columns),
                                                                       size_(shape_.size())
        {
            new_matrix();
            if (ptr != nullptr && size_ > 0)
            {
                std::copy(std::execution::par_unseq, ptr, ptr + size_, begin());
            }
        }

        //============================================================================
        // Method Description:
        /// Constructor. Operates as a shell around an already existing
        /// array of data.
        ///
        /// @param ptr: pointer to beginning of the array
        /// @param size: the number of elements in the array
        /// @param takeOwnership: whether or not to take ownership of the data
        /// and call delete[] in the destructor.
        ///
        Matrix(pointer ptr, size_t size, bool takeOwnership) noexcept : shape_(1, size),
                                                                        size_(size),
                                                                        array_(ptr),
                                                                        ownsPtr_(takeOwnership)
        {
        }

        //============================================================================
        // Method Description:
        /// Constructor. Operates as a shell around an already existing
        /// array of data.
        ///
        /// @param ptr: pointer to beginning of the array
        /// @param rows: the number of rows in the array
        /// @param columns: the number of column in the array
        /// @param takeOwnership: whether or not to take ownership of the data
        /// and call delete[] in the destructor.
        ///
        Matrix(pointer ptr, size_t rows, size_t columns, bool takeOwnership) noexcept : shape_(rows, columns),
                                                                                        size_(rows * columns),
                                                                                        array_(ptr),
                                                                                        ownsPtr_(takeOwnership)
        {
        }

        //============================================================================
        // Method Description:
        /// Copy Constructor
        ///
        /// @param otherArray
        ///
        Matrix(const Matrix<dtype> &otherArray) : shape_(otherArray.shape_),
                                                  size_(otherArray.size_),
                                                  endianess_(otherArray.endianess_)
        {
            new_matrix();
            if (size_ > 0)
            {
                std::copy(std::execution::par_unseq, otherArray.cbegin(), otherArray.cend(), begin());
            }
        }

        //============================================================================
        // Method Description:
        /// Move Constructor
        ///
        /// @param otherArray
        ///
        Matrix(Matrix<dtype> &&otherArray) noexcept : shape_(otherArray.shape_),
                                                      size_(otherArray.size_),
                                                      endianess_(otherArray.endianess_),
                                                      array_(otherArray.array_),
                                                      ownsPtr_(otherArray.ownsPtr_)
        {
            otherArray.shape_.rows = otherArray.shape_.cols = 0;
            otherArray.size_ = 0;
            otherArray.ownsPtr_ = false;
            otherArray.array_ = nullptr;
        }

        //============================================================================
        // Method Description:
        /// Destructor
        ///
        ~Matrix() noexcept
        {
            clear();
        }

        //============================================================================
        // Method Description:
        /// Assignment operator, performs a deep copy
        ///
        /// @param rhs
        /// @return Matrix<dtype>
        ///
        Matrix<dtype> &operator=(const Matrix<dtype> &rhs)
        {
            if (&rhs != this)
            {
                if (rhs.size_ > 0)
                {
                    new_matrix(rhs.shape_);
                    endianess_ = rhs.endianess_;

                    std::copy(std::execution::par_unseq, rhs.cbegin(), rhs.cend(), begin());
                }
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// Assignment operator, sets the entire array to a single
        /// scaler value.
        ///
        /// @param value
        /// @return Matrix<dtype>
        ///
        Matrix<dtype> &operator=(value_type value) noexcept
        {
            if (array_ != nullptr)
            {
                std::fill(std::execution::par_unseq, begin(), end(), value);
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// Move operator, performs a deep move
        ///
        /// @param rhs
        /// @return Matrix<dtype>
        ///
        Matrix<dtype> &operator=(Matrix<dtype> &&rhs) noexcept
        {
            if (&rhs != this)
            {
                clear();
                shape_ = rhs.shape_;
                size_ = rhs.size_;
                endianess_ = rhs.endianess_;
                array_ = rhs.array_;
                ownsPtr_ = rhs.ownsPtr_;

                rhs.shape_.rows = rhs.shape_.cols = rhs.size_ = 0;
                rhs.array_ = nullptr;
                rhs.ownsPtr_ = false;
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// 1D access operator with no bounds checking
        ///
        /// @param pos
        /// @return value
        ///
        inline reference operator[](int64_t pos)
        {
            //if ((size_t)std::abs(pos) >= size_)
            //{
            //    THROW_INVALID_ARGUMENT("pos out of data range");
            //}

            if (pos < 0)
            {
                pos += size_;
            }

            return array_[pos];
        }

        //============================================================================
        // Method Description:
        /// const 1D access operator with no bounds checking
        ///
        /// @param pos
        /// @return value
        ///
        inline const_reference operator[](int64_t pos) const
        {
            //if ((size_t)std::abs(pos) >= size_)
            //{
            //    THROW_INVALID_ARGUMENT("pos out of data range");
            //}

            if (pos < 0)
            {
                // return value from backward, just like python [-1]
                pos += size_;
            }

            return array_[pos];
        }

        //============================================================================
        // Method Description:
        /// 2D access operator with no bounds checking
        ///
        /// @param rowPos
        /// @param columnPos
        /// @return value
        ///
        inline reference operator()(int64_t rowPos, int64_t columnPos)
        {
            if (rowPos < 0)
            {
                rowPos += shape_.rows;
            }

            if (columnPos < 0)
            {
                columnPos += shape_.cols;
            }

            auto pos = rowPos * shape_.cols + columnPos;
            //if (pos >= size_)
            //{
            //    THROW_INVALID_ARGUMENT("pos out of data range");
            //}

            return array_[pos];
        }

        //============================================================================
        // Method Description:
        /// const 2D access operator with no bounds checking
        ///
        /// @param rowPos
        /// @param columnPos
        /// @return value
        ///
        inline const_reference operator()(int64_t rowPos, int64_t columnPos) const
        {
            if (rowPos < 0)
            {
                rowPos += shape_.rows;
            }

            if (columnPos < 0)
            {
                columnPos += shape_.cols;
            }

            auto pos = rowPos * shape_.cols + columnPos;
            //if (pos >= size_)
            //{
            //    THROW_INVALID_ARGUMENT("pos out of data range");
            //}

            return array_[pos];
        }

        //============================================================================
        // Method Description:
        /// 1D Slicing access operator with bounds checking.
        /// returned array is of the range [start, stop).
        ///
        /// @param slice
        /// @return Matrix
        ///
        inline Matrix<dtype> operator[](const Slice &slice) const
        {
            Slice sliceCopy(slice);

            size_t counter = 0;
            Matrix<dtype> returnArray(1, sliceCopy.element_number(size_));
            for (size_t i = sliceCopy.start; i < sliceCopy.stop; i += sliceCopy.step)
            {
                returnArray[counter++] = at(i);
            }

            return returnArray;
        }

        //============================================================================
        // Method Description:
        /// Returns the values from the input mask
        ///
        /// @param mask
        /// @return Matrix
        ///
        inline Matrix<dtype> operator[](const Matrix<bool> &mask) const
        {
            if (mask.shape() != shape_)
            {
                THROW_INVALID_ARGUMENT(
                    "input mask must have the same shape as the matrix it will be masking.");
            }

            auto indices = mask.flatnonzero();
            auto outArray = Matrix<dtype>(1, indices.size());
            for (size_t i = 0; i < indices.size(); ++i)
            {
                outArray[i] = operator[](indices[i]);
            }

            return outArray;
        }

        //============================================================================
        // Method Description:
        /// Returns the values from the input indices
        ///
        /// @param indices
        /// @return Matrix
        ///
        ///
        template <typename Indices, enable_if_t<is_same_v<Indices, Matrix<size_t>>, int> = 0>
        inline Matrix<dtype> operator[](const Indices &indices) const
        {
            if (indices.max().item() > size_ - 1)
            {
                THROW_INVALID_ARGUMENT("input indices must be less than the array size.");
            }

            auto outArray = Matrix<dtype>(1, static_cast<size_t>(indices.size()));
            size_t i = 0;
            for (auto &index : indices)
            {
                outArray[i++] = operator[](index);
            }

            return outArray;
        }

        //============================================================================
        // Method Description:
        /// 2D Slicing access operator with bounds checking.
        /// returned array is of the range [start, stop).
        ///
        /// @param rowSlice
        /// @param columnSlice
        /// @return Matrix
        ///
        inline Matrix<dtype> operator()(Slice rowSlice, Slice columnSlice) const
        {
            Matrix<dtype> returnArray(rowSlice.element_number(shape_.rows), columnSlice.element_number(shape_.cols));

            size_t rowCounter = 0;
            size_t columnCounter = 0;
            for (size_t row = rowSlice.start; row < rowSlice.stop; row += rowSlice.step)
            {
                for (size_t col = columnSlice.start; col < columnSlice.stop; col += columnSlice.step)
                {
                    returnArray(rowCounter, columnCounter++) = at(row, col);
                }
                columnCounter = 0;
                ++rowCounter;
            }

            return returnArray;
        }

        //============================================================================
        // Method Description:
        /// 2D Slicing access operator with bounds checking.
        /// returned array is of the range [start, stop).
        ///
        /// @param rowSlice
        /// @param columnPos
        /// @return Matrix
        ///
        inline Matrix<dtype> operator()(Slice rowSlice, size_t columnPos) const
        {
            Matrix<dtype> returnArray(rowSlice.element_number(shape_.rows), 1);

            size_t rowCounter = 0;
            for (size_t row = rowSlice.start; row < rowSlice.stop; row += rowSlice.step)
            {
                returnArray(rowCounter++, 0) = at(row, columnPos);
            }

            return returnArray;
        }

        //============================================================================
        // Method Description:
        /// 2D Slicing access operator with bounds checking.
        /// returned array is of the range [start, stop).
        ///
        /// @param rowPos
        /// @param columnSlice
        /// @return Matrix
        ///
        inline Matrix<dtype> operator()(size_t rowPos, Slice columnSlice) const
        {
            Matrix<dtype> returnArray(1, columnSlice.element_number(shape_.cols));

            size_t columnCounter = 0;
            for (size_t col = columnSlice.start; col < columnSlice.stop; col += columnSlice.step)
            {
                returnArray(0, columnCounter++) = at(rowPos, col);
            }

            return returnArray;
        }

        //============================================================================
        // Method Description:
        /// 2D index access operator with bounds checking.
        /// returned array is of the range.
        ///
        /// @param rowIndices
        /// @param columnPos
        /// @return Matrix
        ///
        template <typename Indices,
                  enable_if_t<is_same_v<Indices, Matrix<size_t>> || is_same_v<Indices, Matrix<size_t>>, int> = 0>
        inline Matrix<dtype> operator()(const Indices &rowIndices, size_t columnPos) const
        {
            const Matrix<size_t> columnIndices = {columnPos};
            return operator()(rowIndices, columnIndices);
        }

        //============================================================================
        // Method Description:
        /// 2D index access operator with bounds checking.
        /// returned array is of the range.
        ///
        /// @param rowIndices
        /// @param columnSlice
        /// @return Matrix
        ///
        template <typename Indices,
                  enable_if_t<is_same_v<Indices, Matrix<size_t>> || is_same_v<Indices, Matrix<size_t>>, int> = 0>
        inline Matrix<dtype> operator()(const Indices &rowIndices, Slice columnSlice) const
        {
            return operator()(rowIndices, to_indices(columnSlice, Axis::COLUMN));
        }

        //============================================================================
        // Method Description:
        /// 2D index access operator with bounds checking.
        /// returned array is of the range.
        ///
        /// @param rowIndex
        /// @param columnIndices
        /// @return Matrix
        ///
        template <typename Indices,
                  enable_if_t<is_same_v<Indices, Matrix<size_t>> || is_same_v<Indices, Matrix<size_t>>, int> = 0>
        inline Matrix<dtype> operator()(size_t rowIndex, const Indices &columnIndices) const
        {
            const Matrix<size_t> rowIndices = {rowIndex};
            return operator()(rowIndices, columnIndices);
        }

        //============================================================================
        // Method Description:
        /// 2D index access operator with bounds checking.
        /// returned array is of the range.
        ///
        /// @param rowSlice
        /// @param columnIndices
        /// @return Matrix
        ///
        template <typename Indices,
                  enable_if_t<is_same_v<Indices, Matrix<size_t>> || is_same_v<Indices, Matrix<size_t>>, int> = 0>
        inline Matrix<dtype> operator()(Slice rowSlice, const Indices &columnIndices) const
        {
            return operator()(to_indices(rowSlice, Axis::ROW), columnIndices);
        }

        //============================================================================
        // Method Description:
        /// 2D index access operator with bounds checking.
        /// returned array is of the range.
        ///
        /// @param rowIndices
        /// @param columnIndices
        /// @return Matrix
        ///
        template <typename RowIndices,
                  typename ColumnIndices,
                  enable_if_t<is_same_v<RowIndices, Matrix<size_t>> || is_same_v<RowIndices, Matrix<size_t>>, int> = 0,
                  enable_if_t<is_same_v<ColumnIndices, Matrix<size_t>> || is_same_v<ColumnIndices, Matrix<size_t>>, int> = 0>
        inline Matrix<dtype> operator()(RowIndices rowIndices, ColumnIndices columnIndices) const
        {
            rowIndices.sort();
            columnIndices.sort();

            std::vector<size_t> rowIndicesUnique(rowIndices.size());
            std::vector<size_t> columnIndicesUnique(columnIndices.size());

            const auto lastRow =
                std::unique_copy(std::execution::par_unseq, rowIndices.begin(), rowIndices.end(), rowIndicesUnique.begin());
            const auto lastColumn =
                std::unique_copy(std::execution::par_unseq, columnIndices.begin(), columnIndices.end(), columnIndicesUnique.begin());

            Matrix<dtype> returnArray((lastRow - rowIndicesUnique.begin()),
                                      (lastColumn - columnIndicesUnique.begin()));

            size_t rowCounter = 0;
            for (auto rowIter = rowIndicesUnique.begin(); rowIter != lastRow; ++rowIter)
            {
                size_t columnCounter = 0;
                for (auto columnIter = columnIndicesUnique.begin(); columnIter != lastColumn; ++columnIter)
                {
                    returnArray(rowCounter, columnCounter++) = at(*rowIter, *columnIter);
                }

                ++rowCounter;
            }

            return returnArray;
        }

        //============================================================================
        // Method Description:
        /// Returns a Slice object for slicing a row to the end of
        /// array.
        ///
        /// @param startIndex (default 0)
        /// @param step (default 1)
        /// @return Slice
        ///
        Slice column_slice(size_t startIndex = 0, size_t step = 1) const noexcept
        {
            return Slice(startIndex, shape_.cols, step);
        }

        //============================================================================
        // Method Description:
        /// Returns a Slice object for slicing a column to the end
        /// of the array.
        ///
        /// @param startIndex (default 0)
        /// @param step (default 1)
        /// @return Slice
        ///
        Slice row_slice(size_t startIndex = 0, size_t step = 1) const noexcept
        {
            return Slice(startIndex, shape_.rows, step);
        }

        //============================================================================
        // Method Description:
        /// 1D access method with bounds checking
        ///
        /// @param pos
        /// @return value
        ///
        inline reference at(int64_t pos)
        {
            if ((size_t)std::abs(pos) >= size_)
            {
                std::string errStr = "Input index " + std::to_string(pos);
                errStr += " is out of bounds for array of size " + std::to_string(size_) + ".";
                THROW_INVALID_ARGUMENT(errStr);
            }

            return operator[](pos);
        }

        //============================================================================
        // Method Description:
        /// const 1D access method with bounds checking
        ///
        /// @param pos
        /// @return value
        ///
        inline const_reference at(int64_t pos) const
        {
            if (std::abs(pos) >= size_)
            {
                std::string errStr = "Input index " + std::to_string(pos);
                errStr += " is out of bounds for array of size " + std::to_string(size_) + ".";
                THROW_INVALID_ARGUMENT(errStr);
            }

            return operator[](pos);
        }

        //============================================================================
        // Method Description:
        /// 2D access method with bounds checking
        ///
        /// @param rowPos
        /// @param columnPos
        /// @return value
        ///
        inline reference at(int64_t rowPos, int64_t columnPos)
        {
            // this doesn't allow for calling the first element as -size_...
            // but why would you really want to do that anyway?
            if ((size_t)std::abs(rowPos) >= shape_.rows)
            {
                std::string errStr = "Row index " + std::to_string(rowPos);
                errStr += " is out of bounds for array of size " + std::to_string(shape_.rows) + ".";
                THROW_INVALID_ARGUMENT(errStr);
            }

            // this doesn't allow for calling the first element as -size_...
            // but why would you really want to that anyway?
            if ((size_t)std::abs(columnPos) >= shape_.cols)
            {
                std::string errStr = "Column index " + std::to_string(columnPos);
                errStr += " is out of bounds for array of size " + std::to_string(shape_.cols) + ".";
                THROW_INVALID_ARGUMENT(errStr);
            }

            return operator()(rowPos, columnPos);
        }

        //============================================================================
        // Method Description:
        /// const 2D access method with bounds checking
        ///
        /// @param rowPos
        /// @param columnPos
        /// @return value
        ///
        inline const_reference at(int64_t rowPos, int64_t columnPos) const
        {
            // this doesn't allow for calling the first element as -size_...
            // but why would you really want to do that anyway?
            if ((size_t)std::abs(rowPos) >= shape_.rows)
            {
                std::string errStr = "Row index " + std::to_string(rowPos);
                errStr += " is out of bounds for array of size " + std::to_string(shape_.rows) + ".";
                THROW_INVALID_ARGUMENT(errStr);
            }

            // this doesn't allow for calling the first element as -size_...
            // but why would you really want to do that anyway?
            if ((size_t)std::abs(columnPos) >= shape_.cols)
            {
                std::string errStr = "Column index " + std::to_string(columnPos);
                errStr += " is out of bounds for array of size " + std::to_string(shape_.cols) + ".";
                THROW_INVALID_ARGUMENT(errStr);
            }

            return operator()(rowPos, columnPos);
        }

        //============================================================================
        // Method Description:
        /// 2D access method with bounds checking
        ///
        /// @param rowPos
        /// @param columnPos
        /// @return value
        ///
        std::pair<int32_t, int32_t> get_pos(iterator ci) const
        {
            if (size_ == 0)
                THROW_INVALID_ARGUMENT("Matrix is empty.");

            auto dis = std::distance(iterator(array_), ci);
            auto r = int32_t(dis / shape_.cols);
            auto c = int32_t(dis % shape_.cols);

            return std::make_pair(r, c);
        }

        //============================================================================
        // Method Description:
        /// get <row, col> from index
        ///
        /// @return value
        ///
        std::pair<int32_t, int32_t> get_pos(size_t i) const
        {
            if (size_ == 0)
                THROW_INVALID_ARGUMENT("Matrix is empty.");

            auto r = int32_t(i / shape_.cols);
            auto c = int32_t(i % shape_.cols);

            return std::make_pair(r, c);
        }

        //============================================================================
        // Method Description:
        /// Return index of specials element
        ///
        /// @param rowPos
        /// @param columnPos
        /// @return index
        ///
        const size_t index(int64_t rowPos, int64_t columnPos) const
        {
            if (rowPos < 0)
            {
                rowPos += shape_.rows;
            }

            if (columnPos < 0)
            {
                columnPos += shape_.cols;
            }

            auto pos = rowPos * shape_.cols + columnPos;
            if (pos >= size_)
            {
                THROW_INVALID_ARGUMENT("pos out of data range");
            }

            return pos;
        }

        //============================================================================
        // Method Description:i
        /// const 1D access method with bounds checking
        ///
        /// @param slice
        /// @return Matrix
        ///
        Matrix<dtype> at(const Slice &slice) const
        {
            // the slice operator already provides bounds checking. just including
            // the at method for completeness
            return operator[](slice);
        }

        //============================================================================
        // Method Description:
        /// const 2D access method with bounds checking
        ///
        /// @param rowSlice
        /// @param columnSlice
        /// @return Matrix
        ///
        Matrix<dtype> at(const Slice &rowSlice, const Slice &columnSlice) const
        {
            // the slice operator already provides bounds checking. just including
            // the at method for completeness
            return operator()(rowSlice, columnSlice);
        }

        //============================================================================
        // Method Description:
        /// const 2D access method with bounds checking
        ///
        /// @param rowSlice
        /// @param columnPos
        /// @return Matrix
        ///
        Matrix<dtype> at(const Slice &rowSlice, size_t columnPos) const
        {
            // the slice operator already provides bounds checking. just including
            // the at method for completeness
            return operator()(rowSlice, columnPos);
        }

        //============================================================================
        // Method Description:
        /// const 2D access method with bounds checking
        ///
        /// @param rowPos
        /// @param columnSlice
        /// @return Matrix
        ///
        Matrix<dtype> at(size_t rowPos, const Slice &columnSlice) const
        {
            // the slice operator already provides bounds checking. just including
            // the at method for completeness
            return operator()(rowPos, columnSlice);
        }

        //============================================================================
        // Method Description:
        /// const 2D access method with bounds checking
        ///
        /// @param rowIndices
        /// @param columnIndices
        /// @return Matrix
        ///
        Matrix<dtype> at(const Matrix<size_t> &rowIndices, const Matrix<size_t> &columnIndices) const
        {
            // the slice operator already provides bounds checking. just including
            // the at method for completeness
            return operator()(rowIndices, columnIndices);
        }

        //============================================================================
        // Method Description:
        /// iterator to the beginning of the flattened array
        /// @return iterator
        ///
        iterator begin() noexcept
        {
            return iterator(array_);
        }

        //============================================================================
        // Method Description:
        /// iterator to the beginning of the input row
        ///
        /// @param row
        /// @return iterator
        ///
        iterator begin(size_t row)
        {
            if (row >= shape_.rows)
            {
                THROW_INVALID_ARGUMENT("input row is greater than the number of rows in the array.");
            }

            return begin() += (row * shape_.cols);
        }

        //============================================================================
        // Method Description:
        /// const iterator to the beginning of the flattened array
        /// @return const_iterator
        ///
        const_iterator begin() const noexcept
        {
            return cbegin();
        }

        //============================================================================
        // Method Description:
        /// const iterator to the beginning of the input row
        ///
        /// @param row
        /// @return const_iterator
        ///
        const_iterator begin(size_t row) const
        {
            return cbegin(row);
        }

        //============================================================================
        // Method Description:
        /// const iterator to the beginning of the flattened array
        ///
        /// @return const_iterator
        ///
        const_iterator cbegin() const noexcept
        {
            return const_iterator(array_);
        }

        //============================================================================
        // Method Description:
        /// const iterator to the beginning of the input row
        ///
        /// @param row
        /// @return const_iterator
        ///
        const_iterator cbegin(size_t row) const
        {
            if (row >= shape_.rows)
            {
                THROW_INVALID_ARGUMENT("input row is greater than the number of rows in the array.");
            }

            return cbegin() += (row * shape_.cols);
        }

        //============================================================================
        // Method Description:
        /// column_iterator to the beginning of the flattened array
        /// @return column_iterator
        ///
        column_iterator colbegin() noexcept
        {
            return column_iterator(array_, shape_.rows, shape_.cols);
        }

        //============================================================================
        // Method Description:
        /// column_iterator to the beginning of the input column
        ///
        /// @param column
        /// @return column_iterator
        ///
        column_iterator colbegin(size_t column)
        {
            if (column >= shape_.cols)
            {
                THROW_INVALID_ARGUMENT("input col is greater than the number of cols in the array.");
            }

            return colbegin() += (column * shape_.rows);
        }

        //============================================================================
        // Method Description:
        /// const column_iterator to the beginning of the flattened array
        /// @return const_column_iterator
        ///
        const_column_iterator colbegin() const noexcept
        {
            return ccolbegin();
        }

        //============================================================================
        // Method Description:
        /// const column_iterator to the beginning of the input column
        ///
        /// @param column
        /// @return const_column_iterator
        ///
        const_column_iterator colbegin(size_t column) const
        {
            return ccolbegin(column);
        }

        //============================================================================
        // Method Description:
        /// const_column_iterator to the beginning of the flattened array
        ///
        /// @return const_column_iterator
        ///
        const_column_iterator ccolbegin() const noexcept
        {
            return const_column_iterator(array_, shape_.rows, shape_.cols);
        }

        //============================================================================
        // Method Description:
        /// const_column_iterator to the beginning of the input column
        ///
        /// @param column
        /// @return const_column_iterator
        ///
        const_column_iterator ccolbegin(size_t column) const
        {
            if (column >= shape_.cols)
            {
                THROW_INVALID_ARGUMENT("input col is greater than the number of cols in the array.");
            }

            return ccolbegin() += (column * shape_.rows);
        }

        //============================================================================
        // Method Description:
        /// reverse_iterator to the beginning of the flattened array
        /// @return reverse_iterator
        ///
        reverse_iterator rbegin() noexcept
        {
            return reverse_iterator(end());
        }

        //============================================================================
        // Method Description:
        /// reverse_iterator to the beginning of the input row
        ///
        /// @param row
        /// @return reverse_iterator
        ///
        reverse_iterator rbegin(size_t row)
        {
            if (row >= shape_.rows)
            {
                THROW_INVALID_ARGUMENT("input row is greater than the number of rows in the array.");
            }

            return rbegin() += (shape_.rows - row - 1) * shape_.cols;
        }

        //============================================================================
        // Method Description:
        /// const iterator to the beginning of the flattened array
        /// @return const_iterator
        ///
        const_reverse_iterator rbegin() const noexcept
        {
            return crbegin();
        }

        //============================================================================
        // Method Description:
        /// const iterator to the beginning of the input row
        ///
        /// @param row
        /// @return const_iterator
        ///
        const_reverse_iterator rbegin(size_t row) const
        {
            return crbegin(row);
        }

        //============================================================================
        // Method Description:
        /// const_reverse_iterator to the beginning of the flattened array
        ///
        /// @return const_reverse_iterator
        ///
        const_reverse_iterator crbegin() const noexcept
        {
            return const_reverse_iterator(cend());
        }

        //============================================================================
        // Method Description:
        /// const_reverse_iterator to the beginning of the input row
        ///
        /// @param row
        /// @return const_reverse_iterator
        ///
        const_reverse_iterator crbegin(size_t row) const
        {
            if (row >= shape_.rows)
            {
                THROW_INVALID_ARGUMENT("input row is greater than the number of rows in the array.");
            }

            return crbegin() += (shape_.rows - row - 1) * shape_.cols;
        }

        //============================================================================
        // Method Description:
        /// reverse_column_iterator to the beginning of the flattened array
        /// @return reverse_column_iterator
        ///
        reverse_column_iterator rcolbegin() noexcept
        {
            return reverse_column_iterator(colend());
        }

        //============================================================================
        // Method Description:
        /// reverse_column_iterator to the beginning of the input column
        ///
        /// @param column
        /// @return reverse_column_iterator
        ///
        reverse_column_iterator rcolbegin(size_t column)
        {
            if (column >= shape_.cols)
            {
                THROW_INVALID_ARGUMENT("input col is greater than the number of cols in the array.");
            }

            return rcolbegin() += (shape_.cols - column - 1) * shape_.rows;
        }

        //============================================================================
        // Method Description:
        /// const iterator to the beginning of the flattened array
        /// @return const_iterator
        ///
        const_reverse_column_iterator rcolbegin() const noexcept
        {
            return crcolbegin();
        }

        //============================================================================
        // Method Description:
        /// const iterator to the beginning of the input column
        ///
        /// @param column
        /// @return const_iterator
        ///
        const_reverse_column_iterator rcolbegin(size_t column) const
        {
            return crcolbegin(column);
        }

        //============================================================================
        // Method Description:
        /// const_reverse_column_iterator to the beginning of the flattened array
        ///
        /// @return const_reverse_column_iterator
        ///
        const_reverse_column_iterator crcolbegin() const noexcept
        {
            return const_reverse_column_iterator(ccolend());
        }

        //============================================================================
        // Method Description:
        /// const_reverse_column_iterator to the beginning of the input column
        ///
        /// @param column
        /// @return const_reverse_column_iterator
        ///
        const_reverse_column_iterator crcolbegin(size_t column) const
        {
            if (column >= shape_.cols)
            {
                THROW_INVALID_ARGUMENT("input col is greater than the number of cols in the array.");
            }

            return crcolbegin() += (shape_.cols - column - 1) * shape_.rows;
        }

        //============================================================================
        // Method Description:
        /// iterator to 1 past the end of the flattened array
        /// @return iterator
        ///
        iterator end() noexcept
        {
            return begin() += size_;
        }

        //============================================================================
        // Method Description:
        /// iterator to the 1 past end of the row
        ///
        /// @param row
        /// @return iterator
        ///
        iterator end(size_t row)
        {
            if (row >= shape_.rows)
            {
                THROW_INVALID_ARGUMENT("input row is greater than the number of rows in the array.");
            }

            return begin(row) += shape_.cols;
        }

        //============================================================================
        // Method Description:
        /// const iterator to 1 past the end of the flattened array
        /// @return const_iterator
        ///
        const_iterator end() const noexcept
        {
            return cend();
        }

        //============================================================================
        // Method Description:
        /// const iterator to the 1 past end of the row
        ///
        /// @param row
        /// @return const_iterator
        ///
        const_iterator end(size_t row) const
        {
            return cend(row);
        }

        //============================================================================
        // Method Description:
        /// const iterator to 1 past the end of the flattened array
        ///
        /// @return const_iterator
        ///
        const_iterator cend() const noexcept
        {
            return cbegin() += size_;
        }

        //============================================================================
        // Method Description:
        /// const iterator to 1 past the end of the input row
        ///
        /// @param row
        /// @return const_iterator
        ///
        const_iterator cend(size_t row) const
        {
            if (row >= shape_.rows)
            {
                THROW_INVALID_ARGUMENT("input row is greater than the number of rows in the array.");
            }

            return cbegin(row) += shape_.cols;
        }

        //============================================================================
        // Method Description:
        /// reverse_iterator to 1 past the end of the flattened array
        /// @return reverse_iterator
        ///
        reverse_iterator rend() noexcept
        {
            return rbegin() += size_;
        }

        //============================================================================
        // Method Description:
        /// reverse_iterator to the 1 past end of the row
        ///
        /// @param row
        /// @return reverse_iterator
        ///
        reverse_iterator rend(size_t row)
        {
            if (row >= shape_.rows)
            {
                THROW_INVALID_ARGUMENT("input row is greater than the number of rows in the array.");
            }

            return rbegin(row) += shape_.cols;
        }

        //============================================================================
        // Method Description:
        /// const_reverse_iterator to 1 past the end of the flattened array
        /// @return const_reverse_iterator
        ///
        const_reverse_iterator rend() const noexcept
        {
            return crend();
        }

        //============================================================================
        // Method Description:
        /// const_reverse_iterator to the 1 past end of the row
        ///
        /// @param row
        /// @return const_reverse_iterator
        ///
        const_reverse_iterator rend(size_t row) const
        {
            return crend(row);
        }

        //============================================================================
        // Method Description:
        /// const_reverse_iterator to 1 past the end of the flattened array
        ///
        /// @return const_reverse_iterator
        ///
        const_reverse_iterator crend() const noexcept
        {
            return crbegin() += size_;
        }

        //============================================================================
        // Method Description:
        /// const_reverse_iterator to 1 past the end of the input row
        ///
        /// @param row
        /// @return const_reverse_iterator
        ///
        const_reverse_iterator crend(size_t row) const
        {
            if (row >= shape_.rows)
            {
                THROW_INVALID_ARGUMENT("input row is greater than the number of rows in the array.");
            }

            return crbegin(row) += shape_.cols;
        }

        //============================================================================
        // Method Description:
        /// column_iterator to 1 past the end of the flattened array
        /// @return column_iterator
        ///
        column_iterator colend() noexcept
        {
            return colbegin() += size_;
        }

        //============================================================================
        // Method Description:
        /// column_iterator to the 1 past end of the column
        ///
        /// @param column
        /// @return column_iterator
        ///
        column_iterator colend(size_t column)
        {
            if (column >= shape_.cols)
            {
                THROW_INVALID_ARGUMENT("input col is greater than the number of cols in the array.");
            }

            return colbegin(column) += shape_.rows;
        }

        //============================================================================
        // Method Description:
        /// const column_iterator to 1 past the end of the flattened array
        /// @return const_column_iterator
        ///
        const_column_iterator colend() const noexcept
        {
            return ccolend();
        }

        //============================================================================
        // Method Description:
        /// const column_iterator to the 1 past end of the column
        ///
        /// @param column
        /// @return const_column_iterator
        ///
        const_column_iterator colend(size_t column) const
        {
            return ccolend(column);
        }

        //============================================================================
        // Method Description:
        /// const_column_iterator to 1 past the end of the flattened array
        ///
        /// @return const_column_iterator
        ///
        const_column_iterator ccolend() const noexcept
        {
            return ccolbegin() += size_;
        }

        //============================================================================
        // Method Description:
        /// const_column_iterator to 1 past the end of the input col
        ///
        /// @param column
        /// @return const_column_iterator
        ///
        const_column_iterator ccolend(size_t column) const
        {
            if (column >= shape_.cols)
            {
                THROW_INVALID_ARGUMENT("input col is greater than the number of cols in the array.");
            }

            return ccolbegin(column) += shape_.rows;
        }

        //============================================================================
        // Method Description:
        /// reverse_column_iterator to 1 past the end of the flattened array
        /// @return reverse_column_iterator
        ///
        reverse_column_iterator rcolend() noexcept
        {
            return rcolbegin() += size_;
        }

        //============================================================================
        // Method Description:
        /// reverse_column_iterator to the 1 past end of the column
        ///
        /// @param column
        /// @return reverse_column_iterator
        ///
        reverse_column_iterator rcolend(size_t column)
        {
            if (column >= shape_.cols)
            {
                THROW_INVALID_ARGUMENT("input col is greater than the number of cols in the array.");
            }

            return rcolbegin(column) += shape_.rows;
        }

        //============================================================================
        // Method Description:
        /// const_reverse_column_iterator to 1 past the end of the flattened array
        /// @return const_reverse_column_iterator
        ///
        const_reverse_column_iterator rcolend() const noexcept
        {
            return crcolend();
        }

        //============================================================================
        // Method Description:
        /// const_reverse_column_iterator to the 1 past end of the column
        ///
        /// @param column
        /// @return const_reverse_column_iterator
        ///
        const_reverse_column_iterator rcolend(size_t column) const
        {
            return crcolend(column);
        }

        //============================================================================
        // Method Description:
        /// const_reverse_column_iterator to 1 past the end of the flattened array
        ///
        /// @return const_reverse_column_iterator
        ///
        const_reverse_column_iterator crcolend() const noexcept
        {
            return crcolbegin() += size_;
        }

        //============================================================================
        // Method Description:
        /// const_reverse_column_iterator to 1 past the end of the input col
        ///
        /// @param column
        /// @return const_reverse_column_iterator
        ///
        const_reverse_column_iterator crcolend(size_t column) const
        {
            if (column >= shape_.cols)
            {
                THROW_INVALID_ARGUMENT("input col is greater than the number of cols in the array.");
            }

            return crcolbegin(column) += shape_.rows;
        }

        //============================================================================
        // Method Description:
        /// Returns True if all elements evaluate to True or non zero
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<bool> all(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto function = [](dtype i) -> bool
            { return i != dtype{0}; };

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<bool> returnArray = {std::all_of(std::execution::par_unseq, cbegin(), cend(), function)};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<bool> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) = std::all_of(std::execution::par_unseq, cbegin(row), cend(row), function);
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().all(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Returns True if all elements evaluate to True or non zero
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<bool> all(size_t begin_index, size_t end_index, Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto function = [](dtype i) -> bool
            { return i != dtype{0}; };

            if ((cbegin() + begin_index) > cend() || (cbegin() + end_index) > cend())
            {
                return {};
            }

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<bool> returnArray = {std::all_of(std::execution::par_unseq, cbegin() + begin_index, cend() + end_index, function)};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
            }
            case Axis::ROW:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Returns True if any elements evaluate to True or non zero
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<bool> any(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto function = [](dtype i) -> bool
            { return i != dtype{0}; };

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<bool> returnArray = {std::any_of(std::execution::par_unseq, cbegin(), cend(), function)};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<bool> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) = std::any_of(std::execution::par_unseq, cbegin(row), cend(row), function);
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().any(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Return indices of the maximum values along the given axis.
        /// Only the first index is returned.
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<size_t> argmax(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto comparator = [](dtype lhs, dtype rhs) noexcept -> bool
            { return lhs < rhs; };

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<size_t> returnArray = {(
                    std::max_element(std::execution::par_unseq, cbegin(), cend(), comparator) - cbegin())};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<size_t> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) = (std::max_element(std::execution::par_unseq, cbegin(row), cend(row), comparator) - cbegin(row));
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().argmax(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Return indices of the minimum values along the given axis.
        /// Only the first index is returned.
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<size_t> argmin(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto comparator = [](dtype lhs, dtype rhs) noexcept -> bool
            { return lhs < rhs; };

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<size_t> returnArray = {(
                    std::min_element(std::execution::par_unseq, cbegin(), cend(), comparator) - cbegin())};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<size_t> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) = (std::min_element(std::execution::par_unseq, cbegin(row), cend(row), comparator) - cbegin(row));
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().argmin(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Returns the indices that would sort this array.
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<size_t> argsort(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            switch (axis)
            {
            case Axis::NONE:
            {
                std::vector<size_t> idx(size_);
                std::iota(idx.begin(), idx.end(), 0);

                const auto function = [this](size_t i1, size_t i2) noexcept -> bool
                { return (*this)[i1] < (*this)[i2]; };

                std::stable_sort(idx.begin(), idx.end(), function);
                return Matrix<size_t>(idx);
            }
            case Axis::COLUMN:
            {
                Matrix<size_t> returnArray(shape_);
                std::vector<size_t> idx(shape_.cols);

                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    std::iota(idx.begin(), idx.end(), 0);

                    const auto function = [this, row](size_t i1, size_t i2) noexcept -> bool
                    { return operator()(row, i1) < operator()(row, i2); };

                    std::stable_sort(idx.begin(), idx.end(), function);

                    for (size_t col = 0; col < shape_.cols; ++col)
                    {
                        returnArray(row, col) = idx[col];
                    }
                }
                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().argsort(Axis::COLUMN).transpose();
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Returns a copy of the array, cast to a specified type.
        /// Arithmetic to Arithmetic
        ///
        /// @return Matrix
        ///
        template <typename dtypeOut,
                  typename dtype_ = dtype,
                  enable_if_t<is_same_v<dtype_, dtype>, int> = 0,
                  enable_if_t<is_arithmetic_v<dtype_>, int> = 0,
                  enable_if_t<is_arithmetic_v<dtypeOut>, int> = 0>
        Matrix<dtypeOut> astype() const
        {
            Matrix<dtypeOut> outArray(shape_);
            std::transform(std::execution::par_unseq, cbegin(),
                           cend(),
                           outArray.begin(),
                           [](dtype value) -> dtypeOut
                           { return static_cast<dtypeOut>(value); });

            return outArray;
        }

        //============================================================================
        // Method Description:
        /// Returns a copy of the array, cast to a specified type.
        /// Arithmetic to Complex
        ///
        /// @return Matrix
        ///
        template <typename dtypeOut,
                  typename dtype_ = dtype,
                  enable_if_t<is_same_v<dtype_, dtype>, int> = 0,
                  enable_if_t<is_arithmetic_v<dtype_>, int> = 0,
                  enable_if_t<is_complex_v<dtypeOut>, int> = 0>
        Matrix<dtypeOut> astype() const
        {
            Matrix<dtypeOut> outArray(shape_);

            const auto function = [](const_reference value) -> dtypeOut
            { return std::complex<typename dtypeOut::value_type>(value); };

            std::transform(std::execution::par_unseq, cbegin(), cend(), outArray.begin(), function);

            return outArray;
        }

        //============================================================================
        // Method Description:
        /// Returns a copy of the array, cast to a specified type.
        /// Complex to Complex
        ///
        /// @return Matrix
        ///
        template <typename dtypeOut,
                  typename dtype_ = dtype,
                  enable_if_t<is_same_v<dtype_, dtype>, int> = 0,
                  enable_if_t<is_complex_v<dtype_>, int> = 0,
                  enable_if_t<is_complex_v<dtypeOut>, int> = 0>
        Matrix<dtypeOut> astype() const
        {
            Matrix<dtypeOut> outArray(shape_);

            if (is_same_v<dtypeOut, dtype>)
            {
                std::copy(std::execution::par_unseq, cbegin(), cend(), outArray.begin());
            }
            else
            {
                const auto function = [](const_reference value) noexcept -> dtypeOut
                { return complex_cast<typename dtypeOut::value_type>(value); };

                std::transform(std::execution::par_unseq, cbegin(), cend(), outArray.begin(), function);
            }

            return outArray;
        }

        //============================================================================
        // Method Description:
        /// Returns a copy of the array, cast to a specified type.
        /// Complex to Arithmetic
        ///
        /// @return Matrix
        ///
        template <typename dtypeOut,
                  typename dtype_ = dtype,
                  enable_if_t<is_same_v<dtype_, dtype>, int> = 0,
                  enable_if_t<is_complex_v<dtype_>, int> = 0,
                  enable_if_t<is_arithmetic_v<dtypeOut>, int> = 0>
        Matrix<dtypeOut> astype() const
        {
            Matrix<dtypeOut> outArray(shape_);

            const auto function = [](const_reference value) -> dtypeOut
            { return static_cast<dtypeOut>(value.real()); };

            std::transform(std::execution::par_unseq, cbegin(), cend(), outArray.begin(), function);

            return outArray;
        }

        //============================================================================
        // Method Description:
        /// Returns a copy of the last element of the flattened array.
        ///
        /// @return dtype
        ///
        const_reference back() const noexcept
        {
            return *(cend() - 1);
        }

        //============================================================================
        // Method Description:
        /// Returns a reference the last element of the flattened array.
        ///
        /// @return dtype
        ///
        reference back() noexcept
        {
            return *(end() - 1);
        }

        //============================================================================
        // Method Description:
        /// Returns a copy of the last element of the input row.
        ///
        /// @return dtype
        ///
        const_reference back(size_t row) const
        {
            return *(cend(row) - 1);
        }

        //============================================================================
        // Method Description:
        /// Returns a reference the last element of the input row.
        ///
        /// @return dtype
        ///
        reference back(size_t row)
        {
            return *(end(row) - 1);
        }

        //============================================================================
        // Method Description:
        /// Swap the bytes of the array elements in place
        ///
        /// @return Matrix
        ///
        Matrix<dtype> &byteswap() noexcept
        {
            STATIC_ASSERT_INTEGER(dtype);

            std::for_each(std::execution::par_unseq, begin(),
                          end(),
                          [](dtype &value) noexcept -> void
                          { value = endian::byteSwap(value); });

            switch (endianess_)
            {
            case Endian::NATIVE:
            {
                endianess_ = endian::isLittleEndian() ? Endian::BIG : Endian::LITTLE;
                break;
            }
            case Endian::LITTLE:
            {
                endianess_ = Endian::BIG;
                break;
            }
            case Endian::BIG:
            {
                endianess_ = Endian::LITTLE;
                break;
            }
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// Returns an array whose values are limited to [min, max].
        ///
        /// @param min: min value to clip to
        /// @param max: max value to clip to
        /// @return clipped value
        ///
        Matrix<dtype> clip(value_type min, value_type max) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            Matrix<dtype> outArray(shape_);
            std::transform(std::execution::par_unseq, cbegin(),
                           cend(),
                           outArray.begin(),
                           [min, max](dtype value) noexcept -> dtype
                           {
#ifdef __cpp_lib_clamp
                               const auto comparator = [](dtype lhs, dtype rhs) noexcept -> bool
                               { return lhs < rhs; };

                               return std::clamp(value, min, max, comparator);
#else
                    if (value < min)
                    {
                        return min;
                    }
                    else if (value > max)
                    {
                        return max;
                    }

                    return value;
#endif
                           });

            return outArray;
        }

        //============================================================================
        // Method Description:
        /// Returns the full column of the array
        ///
        ///
        /// @return Shape
        ///
        Matrix<dtype> column(size_t column)
        {
            return operator()(row_slice(), column);
        }

        //============================================================================
        // Method Description:
        /// returns whether or not a value is included the array
        ///
        /// @param value
        /// @param axis (Optional, default NONE)
        /// @return bool
        ///
        Matrix<bool> contains(value_type value, Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<bool> returnArray = {std::find(std::execution::par_unseq, cbegin(), cend(), value) != cend()};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<bool> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) = std::find(std::execution::par_unseq, cbegin(row), cend(row), value) != cend(row);
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().contains(value, Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// returns whether or not a value is included the array
        ///
        /// @param value
        /// @param axis (Optional, default NONE)
        /// @return bool
        ///
        bool containsB(value_type value, Axis axis = Axis::NONE) const
        {
            Matrix<bool> returnArray = contains(value, axis);

            return std::any_of(std::execution::par_unseq, returnArray.begin(), returnArray.end(), [](bool v)
                               { return v; });
        }

        //============================================================================
        // Method Description:
        /// Return a copy of the array
        ///
        /// @return Matrix
        ///
        Matrix<dtype> copy() const
        {
            return Matrix<dtype>(*this);
        }

        //============================================================================
        // Method Description:
        /// Return the cumulative product of the elements along the given axis.
        ///
        /// [[1,2], [3,4]] =(Axis::NONE)=> [1,2,4,12]
        /// [[1,2], [3,4]] =(Axis::0)=> [[1,2],[3,8]]   cumprod by row
        /// [[1,2], [3,4]] =(Axis::1)=> [[1,2],[3,12]]  cumprod by column
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<dtype> cumprod(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<dtype> returnArray(1, size_);
                returnArray[0] = front();
                for (size_t i = 1; i < size_; ++i)
                {
                    returnArray[i] = returnArray[i - 1] * array_[i];
                }

                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<dtype> returnArray(shape_);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(row, 0) = operator()(row, 0);
                    for (size_t col = 1; col < shape_.cols; ++col)
                    {
                        returnArray(row, col) = returnArray(row, col - 1) * operator()(row, col);
                    }
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().cumprod(Axis::COLUMN).transpose();
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Return the cumulative sum of the elements along the given axis.
        ///
        /// [[1,2], [3,4]] =(Axis::NONE)=> [1,3,6,10]
        /// [[1,2], [3,4]] =(Axis::0)=> [[1,2],[4,6]]  cumsum by row
        /// [[1,2], [3,4]] =(Axis::1)=> [[1,3],[3,7]]  cumsum by column
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<dtype> cumsum(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<dtype> returnArray(1, size_);
                returnArray[0] = front();
                for (size_t i = 1; i < size_; ++i)
                {
                    returnArray[i] = returnArray[i - 1] + array_[i];
                }

                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<dtype> returnArray(shape_);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(row, 0) = operator()(row, 0);
                    for (size_t col = 1; col < shape_.cols; ++col)
                    {
                        returnArray(row, col) = returnArray(row, col - 1) + operator()(row, col);
                    }
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().cumsum(Axis::COLUMN).transpose();
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Returns the raw pointer to the underlying data
        /// @return pointer
        ///
        pointer data() noexcept
        {
            return array_;
        }

        //============================================================================
        // Method Description:
        /// Returns the raw pointer to the underlying data
        /// @return const_pointer
        ///
        const_pointer data() const noexcept
        {
            return array_;
        }

        //============================================================================
        // Method Description:
        /// Releases the internal data pointer so that the destructor
        /// will not call delete on it, and returns the raw pointer
        /// to the underlying data.
        /// @return pointer
        ///
        pointer release_data() noexcept
        {
            ownsPtr_ = false;
            return data();
        }

        //============================================================================
        // Method Description:
        /// Return specified diagonals.
        ///
        /// @param offset: Offset of the diagonal from the main diagonal. Can be both positive and negative. Defaults
        /// to 0.
        /// @param axis: (Optional, default ROW) axis the offset is applied to
        /// @return Matrix
        ///
        Matrix<dtype> diagonal(int64_t offset = 0, Axis axis = Axis::ROW) const
        {
            switch (axis)
            {
            case Axis::COLUMN:
            {
                std::vector<dtype> diagonalValues;
                size_t col = 0;
                for (size_t row = offset; row < shape_.rows; ++row)
                {
                    if (row < 0)
                    {
                        ++col;
                        continue;
                    }
                    if (col >= shape_.cols)
                    {
                        break;
                    }

                    diagonalValues.push_back(operator()(row, col));
                    ++col;
                }

                return Matrix<dtype>(diagonalValues);
            }
            case Axis::ROW:
            {
                return transpose().diagonal(offset, Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Dot product of two arrays.
        ///
        /// For 2-D arrays it is equivalent to matrix multiplication,
        /// and for 1-D arrays to inner product of vectors.
        ///
        ///  [1, 2], dot [1, 2], ==> [1*1 + 2*2, 1*1 + 2*2], = [5, 4],
        ///  [3, 4]      [2, 1]      [3*1 + 4*2, 3*2 + 4*1]    [11, 10]
        ///
        /// @param otherArray
        /// @return dot product
        ///
        Matrix<dtype> dot(const Matrix<dtype> &otherArray) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            if (shape_ == otherArray.shape_ && (shape_.rows == 1 || shape_.cols == 1))
            {
                dtype dotProduct = std::inner_product(cbegin(), cend(), otherArray.cbegin(), dtype{0});
                Matrix<dtype> returnArray = {dotProduct};
                return returnArray;
            }
            if (shape_.cols == otherArray.shape_.rows)
            {
                // 2D array, use matrix multiplication
                Matrix<dtype> returnArray(shape_.rows, otherArray.shape_.cols);
                auto otherArrayT = otherArray.transpose();

                for (size_t i = 0; i < shape_.rows; ++i)
                {
                    for (size_t j = 0; j < otherArrayT.shape_.rows; ++j)
                    {
                        returnArray(i, j) =
                            std::inner_product(otherArrayT.cbegin(j), otherArrayT.cend(j), cbegin(i), dtype{0});
                    }
                }

                return returnArray;
            }

            std::string errStr = "shapes of [" + std::to_string(shape_.rows) + ", " + std::to_string(shape_.cols) + "]";
            errStr += " and [" + std::to_string(otherArray.shape_.rows) + ", " +
                      std::to_string(otherArray.shape_.cols) + "]";
            errStr += " are not consistent.";
            THROW_INVALID_ARGUMENT(errStr);

            return Matrix<dtype>(); // get rid of compiler warning
        }

        //============================================================================
        // Method Description:
        /// Dump a binary file of the array to the specified file.
        /// The array can be read back with ais::load.
        ///
        /// @param filename
        ///
        void dump(const std::string &filename) const
        {
            filesystem::File f(filename);
            if (!f.hasExt())
            {
                f.withExt(".bin");
            }

            std::ofstream ofile(f.fullname().c_str(), std::ios::binary);
            if (!ofile.good())
            {
                THROW_RUNTIME_ERROR("Unable to open the input file:\n\t" + filename);
            }

            if (array_ != nullptr)
            {
                ofile.write(reinterpret_cast<const char *>(array_), size_ * sizeof(dtype));
            }
            ofile.close();
        }

        //============================================================================
        // Method Description:
        /// Return the matrix's endianess
        ///
        /// @return Endian
        ///
        Endian endianess() const noexcept
        {
            STATIC_ASSERT_ARITHMETIC(dtype);

            return endianess_;
        }

        //============================================================================
        // Method Description:
        /// Fill the array with a scaler value.
        ///
        /// @param fillValue
        /// @return *this
        ///
        Matrix<dtype> &fill(value_type fillValue) noexcept
        {
            std::fill(std::execution::par_unseq, begin(), end(), fillValue);
            return std::forward<Matrix<dtype> &>(*this);
        }

        //============================================================================
        // Method Description:
        /// Return the indices of the flattened array of the
        /// elements that are non-zero.
        ///
        /// @return Matrix
        ///
        Matrix<size_t> flatnonzero() const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            std::vector<size_t> indices;
            size_t idx = 0;
            for (auto value : *this)
            {
                if (value != dtype{0})
                {
                    indices.push_back(idx);
                }
                ++idx;
            }

            return Matrix<size_t>(indices);
        }

        //============================================================================
        // Method Description:
        /// Return a copy of the array collapsed into one dimension.
        ///
        /// @return Matrix
        ///
        Matrix<dtype> flatten() const
        {
            Matrix<dtype> outArray(1, size_);
            std::copy(std::execution::par_unseq, cbegin(), cend(), outArray.begin());
            return outArray;
        }

        //============================================================================
        // Method Description:
        /// Returns a copy of the first element of the flattened array.
        ///
        /// @return dtype
        ///
        const_reference front() const noexcept
        {
            return *cbegin();
        }

        //============================================================================
        // Method Description:
        /// Returns a reference to the first element of the flattened array.
        ///
        /// @return dtype
        ///
        reference front() noexcept
        {
            return *begin();
        }

        //============================================================================
        // Method Description:
        /// Returns a copy of the first element of the input row.
        ///
        /// @return dtype
        ///
        const_reference front(size_t row) const
        {
            return *cbegin(row);
        }

        //============================================================================
        // Method Description:
        /// Returns a reference to the first element of the input row.
        ///
        /// @return dtype
        ///
        reference front(size_t row)
        {
            return *begin(row);
        }

        //============================================================================
        // Method Description:
        /// Returns a new flat array with the givin flat input indices.
        ///
        /// @param indices
        /// @return values
        ///
        Matrix<dtype> get_by_indices(const Matrix<size_t> &indices) const
        {
            return operator[](indices);
        }

        //============================================================================
        // Method Description:
        /// Takes in a boolean mask the same size as the array
        /// and returns a flattened array with the values corresponding
        /// to the input mask.
        ///
        /// @param mask
        /// @return values
        ///
        Matrix<dtype> get_by_mask(const Matrix<bool> &mask) const
        {
            return operator[](mask);
        }

        //============================================================================
        // Method Description:
        /// Return if the Matrix is empty. ie the default constructor
        /// was used.
        ///
        /// @return boolean
        ///
        bool is_empty() const noexcept
        {
            return size_ == 0;
        }

        //============================================================================
        // Method Description:
        /// Return if the Matrix is flat. ie the number of columns or
        /// rows is equal to one.
        ///
        /// @return boolean
        ///
        bool is_flat() const noexcept
        {
            return shape_.rows == 1 || shape_.cols == 1;
        }

        //============================================================================
        // Method Description:
        /// Return if the Matrix is sorted.
        ///
        /// @param axis
        /// @return boolean
        ///
        Matrix<bool> is_sorted(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto comparator = [](dtype lhs, dtype rhs) noexcept -> bool
            { return lhs < rhs; };

            switch (axis)
            {
            case Axis::NONE:
            {
                return {std::is_sorted(std::execution::par_unseq, cbegin(), cend(), comparator)};
            }
            case Axis::COLUMN:
            {
                Matrix<bool> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) = std::is_sorted(std::execution::par_unseq, cbegin(row), cend(row), comparator);
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().is_sorted(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Return if the Matrix is square.
        ///
        /// @return boolean
        ///
        bool is_square() const noexcept
        {
            return shape_.is_square();
        }

        //============================================================================
        // Method Description:
        /// Copy an element of an array to a standard C++ scaler and return it.
        ///
        /// @return array element
        ///
        value_type item() const
        {
            if (size_ != 1)
            {
                THROW_INVALID_ARGUMENT("Can only convert an array of size 1 to a C++ scaler");
            }

            return front();
        }

        //============================================================================
        // Method Description:
        /// Return the maximum along a given axis.
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<dtype> max(Axis axis = Axis::NONE) const
        {
            // STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto comparator = [](dtype lhs, dtype rhs) noexcept -> bool
            { return lhs < rhs; };

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<dtype> returnArray = {*std::max_element(std::execution::par_unseq, cbegin(), cend(), comparator)};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<dtype> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) = *std::max_element(std::execution::par_unseq, cbegin(row), cend(row), comparator);
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().max(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Return the minimum along a given axis.
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<dtype> min(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto comparator = [](dtype lhs, dtype rhs) noexcept -> bool
            { return lhs < rhs; };

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<dtype> returnArray = {*std::min_element(std::execution::par_unseq, cbegin(), cend(), comparator)};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<dtype> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) = *std::min_element(std::execution::par_unseq, cbegin(row), cend(row), comparator);
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().min(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Return the median along a given axis.
        /// If the dtype is floating point then the middle elements will be
        /// averaged for arrays of even number of elements.
        /// If the dtype is integral then the middle elements will be integer
        /// averaged (rounded down to integer) for arrays of even number of elements.
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<dtype> median(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto comparator = [](dtype lhs, dtype rhs) noexcept -> bool
            { return lhs < rhs; };

            if (size_ == 0)
            {
                THROW_RUNTIME_ERROR("Median is undefined for an array of size = 0.");
            }

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<dtype> copyArray(*this);

                const size_t middleIdx = size_ / 2; // integer division
                std::nth_element(std::execution::par_unseq,
                                 copyArray.begin(),
                                 copyArray.begin() + middleIdx,
                                 copyArray.end(),
                                 comparator);

                dtype medianValue = copyArray.array_[middleIdx];
                if (size_ % 2 == 0)
                {
                    const size_t lhsIndex = middleIdx - 1;
                    std::nth_element(std::execution::par_unseq,
                                     copyArray.begin(),
                                     copyArray.begin() + lhsIndex,
                                     copyArray.end(),
                                     comparator);
                    medianValue =
                        (medianValue + copyArray.array_[lhsIndex]) / dtype{2}; // potentially integer division, ok
                }

                return {medianValue};
            }
            case Axis::COLUMN:
            {
                Matrix<dtype> copyArray(*this);
                Matrix<dtype> returnArray(1, shape_.rows);

                const bool isEven = shape_.cols % 2 == 0;
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    const size_t middleIdx = shape_.cols / 2; // integer division
                    std::nth_element(std::execution::par_unseq,
                                     copyArray.begin(row),
                                     copyArray.begin(row) + middleIdx,
                                     copyArray.end(row),
                                     comparator);

                    dtype medianValue = copyArray(row, middleIdx);
                    if (isEven)
                    {
                        const size_t lhsIndex = middleIdx - 1;
                        std::nth_element(std::execution::par_unseq,
                                         copyArray.begin(row),
                                         copyArray.begin(row) + lhsIndex,
                                         copyArray.end(row),
                                         comparator);
                        medianValue = (medianValue + copyArray(row, lhsIndex)) / dtype{2}; // potentially integer division, ok
                    }

                    returnArray(0, row) = medianValue;
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().median(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Fills the array with nans.
        ///
        Matrix<dtype> &nans() noexcept
        {
            STATIC_ASSERT_FLOAT(dtype);

            fill(std::nan("1"));
            return *this;
        }

        //============================================================================
        // Method Description:
        /// Fills the array with empty.
        ///
        Matrix<dtype> &emptys() noexcept
        {
            fill({});
            return *this;
        }

        //============================================================================
        // Method Description:
        /// Returns the number of bytes held by the array
        ///
        /// @return number of bytes
        ///
        uint64_t nbytes() const noexcept
        {
            return static_cast<uint64_t>(sizeof(dtype) * size_);
        }

        //============================================================================
        // Method Description:
        /// Return the array with the same data viewed with a
        /// different byte order. only works for integer types.
        ///
        /// @param endianess
        /// @return Matrix
        ///
        Matrix<dtype> new_byte_order(Endian endianess) const
        {
            STATIC_ASSERT_INTEGER(dtype);

            const bool nativeIsLittle = endian::isLittleEndian();

            switch (endianess_)
            {
            case Endian::NATIVE:
            {
                switch (endianess)
                {
                case Endian::NATIVE:
                {
                    return Matrix(*this);
                }
                case Endian::BIG:
                {
                    if (nativeIsLittle)
                    {
                        Matrix<dtype> outArray(shape_);

                        std::transform(std::execution::par_unseq, cbegin(), end(), outArray.begin(), endian::byteSwap<dtype>);

                        outArray.endianess_ = Endian::BIG;
                        return outArray;
                    }
                    else
                    {
                        auto outArray = Matrix(*this);
                        outArray.endianess_ = Endian::BIG;
                        return outArray;
                    }
                }
                case Endian::LITTLE:
                {
                    if (nativeIsLittle)
                    {
                        auto outArray = Matrix(*this);
                        outArray.endianess_ = Endian::LITTLE;
                        return outArray;
                    }
                    else
                    {
                        Matrix<dtype> outArray(shape_);

                        std::transform(std::execution::par_unseq, cbegin(), end(), outArray.begin(), endian::byteSwap<dtype>);

                        outArray.endianess_ = Endian::LITTLE;
                        return outArray;
                    }
                }
                default:
                {
                    THROW_INVALID_ARGUMENT("Unimplemented endian type.");
                    return {}; // get rid of compiler warning
                }
                }
                break;
            }
            case Endian::BIG:
            {
                switch (endianess)
                {
                case Endian::NATIVE:
                {
                    if (nativeIsLittle)
                    {
                        Matrix<dtype> outArray(shape_);

                        std::transform(std::execution::par_unseq, cbegin(), end(), outArray.begin(), endian::byteSwap<dtype>);

                        outArray.endianess_ = Endian::NATIVE;
                        return outArray;
                    }
                    else
                    {
                        auto outArray = Matrix(*this);
                        outArray.endianess_ = Endian::NATIVE;
                        return outArray;
                    }
                }
                case Endian::BIG:
                {
                    return Matrix(*this);
                }
                case Endian::LITTLE:
                {
                    Matrix<dtype> outArray(shape_);

                    std::transform(std::execution::par_unseq, cbegin(), end(), outArray.begin(), endian::byteSwap<dtype>);

                    outArray.endianess_ = Endian::LITTLE;
                    return outArray;
                }
                default:
                {
                    THROW_INVALID_ARGUMENT("Unimplemented endian type.");
                    return {}; // get rid of compiler warning
                }
                }
                break;
            }
            case Endian::LITTLE:
            {
                switch (endianess)
                {
                case Endian::NATIVE:
                {
                    if (nativeIsLittle)
                    {
                        auto outArray = Matrix(*this);
                        outArray.endianess_ = Endian::NATIVE;
                        return outArray;
                    }
                    else
                    {
                        Matrix<dtype> outArray(shape_);

                        std::transform(std::execution::par_unseq, cbegin(), end(), outArray.begin(), endian::byteSwap<dtype>);

                        outArray.endianess_ = Endian::NATIVE;
                        return outArray;
                    }
                }
                case Endian::BIG:
                {
                    Matrix<dtype> outArray(shape_);

                    std::transform(std::execution::par_unseq, cbegin(), end(), outArray.begin(), endian::byteSwap<dtype>);

                    outArray.endianess_ = Endian::BIG;
                    return outArray;
                }
                case Endian::LITTLE:
                {
                    return Matrix(*this);
                }
                default:
                {
                    THROW_INVALID_ARGUMENT("Unimplemented endian type.");
                    return {}; // get rid of compiler warning
                }
                }
                break;
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented endian type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Returns True if none elements evaluate to True or non zero
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<bool> none(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto function = [](dtype i) -> bool
            { return i != dtype{0}; };

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<bool> returnArray = {std::none_of(std::execution::par_unseq, cbegin(), cend(), function)};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<bool> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) = std::none_of(std::execution::par_unseq, cbegin(row), cend(row), function);
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().none(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Return the row/col indices of the array of the
        /// elements that are non-zero.
        ///
        /// @return std::pair<Matrix, Matrix> where first is the row indices and second is the
        /// column indices
        ///
        std::pair<Matrix<size_t>, Matrix<size_t>> nonzero() const;

        //============================================================================
        // Method Description:
        /// Return the row/col indices of the array of the
        /// elements that are non-zero.
        ///
        /// @return std::pair<Matrix, Matrix> where first is the row indices and second is the
        /// column indices
        ///
        std::pair<Matrix<size_t>, Matrix<size_t>> nonempty() const;

        //============================================================================
        // Method Description:
        /// Return the row/col indices of the array of the
        /// elements that are not nan.
        ///
        /// @return std::pair<Matrix, Matrix> where first is the row indices and second is the
        /// column indices
        ///
        std::pair<Matrix<size_t>, Matrix<size_t>> data_area() const;

        //============================================================================
        // Method Description:
        /// Returns the number of columns in the array
        ///
        ///
        /// @return size_t
        ///
        size_t columns() const noexcept
        {
            return shape_.cols;
        }

        //============================================================================
        // Method Description:
        /// Returns the number of rows in the array
        ///
        ///
        /// @return size_t
        ///
        size_t rows() const noexcept
        {
            return shape_.rows;
        }

        //============================================================================
        // Method Description:
        /// Fills the array with ones
        ///
        ///
        Matrix<dtype> &ones() noexcept
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            fill(dtype{1});
            return *this;
        }

        //============================================================================
        // Method Description:
        /// Returns whether or not the array object owns the underlying data
        ///
        /// @return bool
        ///
        bool owns_internal_data() noexcept
        {
            return ownsPtr_;
        }

        //============================================================================
        // Method Description:
        /// Rearranges the elements in the array in such a way that
        /// value of the element in kth position is in the position it
        /// would be in a sorted array. All elements smaller than the kth
        /// element are moved before this element and all equal or greater
        /// are moved behind it. The ordering of the elements in the two
        /// partitions is undefined.
        ///
        /// @param kth: kth element
        /// @param axis (Optional, default NONE)
        /// @return None
        ///
        Matrix<dtype> &partition(size_t kth, Axis axis = Axis::NONE)
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto comparator = [](dtype lhs, dtype rhs) noexcept -> bool
            { return lhs < rhs; };

            switch (axis)
            {
            case Axis::NONE:
            {
                if (kth >= size_)
                {
                    std::string errStr = "kth(=" + std::to_string(kth);
                    errStr += ") out of bounds (" + std::to_string(size_) + ")";
                    THROW_INVALID_ARGUMENT(errStr);
                }

                std::nth_element(std::execution::par_unseq, begin(), begin() + kth, end(), comparator);
                break;
            }
            case Axis::COLUMN:
            {
                if (kth >= shape_.cols)
                {
                    std::string errStr = "kth(=" + std::to_string(kth);
                    errStr += ") out of bounds (" + std::to_string(shape_.cols) + ")";
                    THROW_INVALID_ARGUMENT(errStr);
                }

                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    std::nth_element(std::execution::par_unseq, begin(row), begin(row) + kth, end(row), comparator);
                }
                break;
            }
            case Axis::ROW:
            {
                if (kth >= shape_.rows)
                {
                    std::string errStr = "kth(=" + std::to_string(kth);
                    errStr += ") out of bounds (" + std::to_string(shape_.rows) + ")";
                    THROW_INVALID_ARGUMENT(errStr);
                }

                Matrix<dtype> transposedArray = transpose();
                for (size_t row = 0; row < transposedArray.shape_.rows; ++row)
                {
                    std::nth_element(std::execution::par_unseq, transposedArray.begin(row),
                                     transposedArray.begin(row) + kth,
                                     transposedArray.end(row),
                                     comparator);
                }
                *this = transposedArray.transpose();
                break;
            }
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// Prints the array to the console.
        ///
        ///
        void print() const
        {
            std::cout << *this;
        }

        //============================================================================
        // Method Description:
        /// Return the product of the array elements over the given axis
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<dtype> prod(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            switch (axis)
            {
            case Axis::NONE:
            {
                dtype product = std::accumulate(cbegin(), cend(), dtype{1}, std::multiplies<dtype>());
                Matrix<dtype> returnArray = {product};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<dtype> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) =
                        std::accumulate(cbegin(row), cend(row), dtype{1}, std::multiplies<dtype>());
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().prod(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Peak to peak (maximum - minimum) value along a given axis.
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<dtype> ptp(Axis axis = Axis::NONE) const
        {
            STATIC_ASSERT_ARITHMETIC_OR_COMPLEX(dtype);

            const auto comparator = [](dtype lhs, dtype rhs) noexcept -> bool
            { return lhs < rhs; };

            switch (axis)
            {
            case Axis::NONE:
            {
                const auto result = std::max_element(std::execution::par_unseq, cbegin(), cend(), comparator);
                Matrix<dtype> returnArray = {*result.second - *result.first};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<dtype> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    const auto result = std::max_element(std::execution::par_unseq, cbegin(row), cend(row), comparator);
                    returnArray(0, row) = *result.second - *result.first;
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().ptp(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// set the flat index element to the value
        ///
        /// @param pos
        /// @param value
        ///
        Matrix<dtype> &put(size_t pos, value_type value)
        {
            at(pos) = value;

            return *this;
        }

        //============================================================================
        // Method Description:
        /// set the 2D row/col index element to the value
        ///
        /// @param row
        /// @param column
        /// @param value
        ///
        Matrix<dtype> &put(size_t row, size_t column, value_type value)
        {
            at(row, column) = value;

            return *this;
        }

        //============================================================================
        // Method Description:
        /// set the flat index elements to the values
        ///
        /// @param pos
        /// @param values
        ///
        template <class ct>
        Matrix<dtype> &putMany(size_t pos, const ct &list)
        {
            if ((pos + list.size()) > size_)
            {
                THROW_INVALID_ARGUMENT("pos out of data range");
            }

            const auto ptr = begin() + pos;
            std::copy(std::execution::par_unseq, list.begin(), list.end(), ptr);

            return *this;
        }

        //============================================================================
        // Method Description:
        /// set the 2D row/col index element to the value
        ///
        /// @param row
        /// @param column
        /// @param values
        ///
        template <class ct>
        Matrix<dtype> &putMany(size_t row, size_t column, const ct &list)
        {
            auto pos = row * shape_.cols + column;
            return putMany(pos, list);
        }

        //============================================================================
        // Method Description:
        /// Set a.flat[n] = values for all n in indices.
        ///
        /// @param indices
        /// @param value
        ///
        Matrix<dtype> &put(const Matrix<size_t> &indices, value_type value)
        {
            if (indices.size() > 0)
            {
                for (auto index : indices)
                {
                    put(index, value);
                }
            }
            return *this;
        }

        //============================================================================
        // Method Description:
        /// Set a.flat[n] = values[n] for all n in indices.
        ///
        /// @param indices
        /// @param values
        ///
        Matrix<dtype> &put(const Matrix<size_t> &indices, const Matrix<dtype> &values)
        {
            if (indices.size() != values.size())
            {
                THROW_INVALID_ARGUMENT("Input indices do not match values dimensions.");
            }

            size_t counter = 0;
            if (indices.size() > 0)
            {
                for (auto index : indices)
                {
                    put(index, values[counter++]);
                }
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// Set the slice indices to the input value.
        ///
        /// @param slice
        /// @param value
        ///
        Matrix<dtype> &put(const Slice &slice, value_type value)
        {
            Slice sliceCopy(slice);
            sliceCopy.make_positive_and_validate(size_);

            for (size_t i = sliceCopy.start; i < sliceCopy.stop; i += sliceCopy.step)
            {
                put(i, value);
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// Set the slice indices to the input values.
        ///
        /// @param slice
        /// @param values
        ///
        Matrix<dtype> &put(const Slice &slice, const Matrix<dtype> &values)
        {
            Slice sliceCopy(slice);
            sliceCopy.make_positive_and_validate(size_);

            std::vector<size_t> indices;
            for (size_t i = sliceCopy.start; i < sliceCopy.stop; i += sliceCopy.step)
            {
                indices.push_back(i);
            }

            return put(Matrix<size_t>(indices), values);
        }

        //============================================================================
        // Method Description:
        /// Set the slice indices to the input value.
        ///
        /// @param rowSlice
        /// @param columnSlice
        /// @param value
        ///
        Matrix<dtype> &put(const Slice &rowSlice, const Slice &columnSlice, value_type value)
        {
            Slice rowSliceCopy(rowSlice);
            Slice columnSliceCopy(columnSlice);

            rowSliceCopy.make_positive_and_validate(shape_.rows);
            columnSliceCopy.make_positive_and_validate(shape_.cols);

            std::vector<size_t> indices;
            for (size_t row = rowSliceCopy.start; row < rowSliceCopy.stop; row += rowSliceCopy.step)
            {
                for (size_t col = columnSliceCopy.start; col < columnSliceCopy.stop; col += columnSliceCopy.step)
                {
                    put(row, col, value);
                }
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// Set the slice indices to the input value.
        ///
        /// @param rowSlice
        /// @param columnPos
        /// @param value
        ///
        Matrix<dtype> &put(const Slice &rowSlice, size_t columnPos, value_type value)
        {
            Slice rowSliceCopy(rowSlice);
            rowSliceCopy.make_positive_and_validate(shape_.rows);

            std::vector<size_t> indices;
            for (size_t row = rowSliceCopy.start; row < rowSliceCopy.stop; row += rowSliceCopy.step)
            {
                put(row, columnPos, value);
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// Set the slice indices to the input value.
        ///
        /// @param rowPos
        /// @param columnSlice
        /// @param value
        ///
        Matrix<dtype> &put(size_t rowPos, const Slice &columnSlice, value_type value)
        {
            Slice columnSliceCopy(columnSlice);
            columnSliceCopy.make_positive_and_validate(shape_.cols);

            std::vector<size_t> indices;
            for (size_t col = columnSliceCopy.start; col < columnSliceCopy.stop; col += columnSliceCopy.step)
            {
                put(rowPos, col, value);
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// Set the slice indices to the input values.
        ///
        /// @param rowSlice
        /// @param columnSlice
        /// @param values
        ///
        Matrix<dtype> &put(const Slice &rowSlice, const Slice &columnSlice, const Matrix<dtype> &values)
        {
            Slice rowSliceCopy(rowSlice);
            Slice columnSliceCopy(columnSlice);

            rowSliceCopy.make_positive_and_validate(shape_.rows);
            columnSliceCopy.make_positive_and_validate(shape_.cols);

            std::vector<size_t> indices;
            for (size_t row = rowSliceCopy.start; row < rowSliceCopy.stop; row += rowSliceCopy.step)
            {
                for (size_t col = columnSliceCopy.start; col < columnSliceCopy.stop; col += columnSliceCopy.step)
                {
                    const size_t index = row * shape_.cols + col;
                    indices.push_back(index);
                }
            }

            return put(Matrix<size_t>(indices), values);
        }

        //============================================================================
        // Method Description:
        /// Set the slice indices to the input values.
        ///
        /// @param rowSlice
        /// @param columnPos
        /// @param values
        ///
        Matrix<dtype> &put(const Slice &rowSlice, size_t columnPos, const Matrix<dtype> &values)
        {
            Slice rowSliceCopy(rowSlice);
            rowSliceCopy.make_positive_and_validate(shape_.rows);

            std::vector<size_t> indices;
            for (size_t row = rowSliceCopy.start; row < rowSliceCopy.stop; row += rowSliceCopy.step)
            {
                const size_t index = row * shape_.cols + columnPos;
                indices.push_back(index);
            }

            return put(Matrix<size_t>(indices), values);
        }

        //============================================================================
        // Method Description:
        /// Set the slice indices to the input values.
        ///
        /// @param rowPos
        /// @param columnSlice
        /// @param values
        ///
        Matrix<dtype> &put(size_t rowPos, const Slice &columnSlice, const Matrix<dtype> &values)
        {
            Slice columnSliceCopy(columnSlice);
            columnSliceCopy.make_positive_and_validate(shape_.cols);

            std::vector<size_t> indices;
            for (size_t col = columnSliceCopy.start; col < columnSliceCopy.stop; col += columnSliceCopy.step)
            {
                const size_t index = rowPos * shape_.cols + col;
                indices.push_back(index);
            }

            return put(Matrix<size_t>(indices), values);
        }

        //============================================================================
        // Method Description:
        /// Set the mask indices to the input value.
        ///
        /// @param mask
        /// @param value
        ///
        Matrix<dtype> &put_mask(const Matrix<bool> &mask, value_type value)
        {
            if (mask.shape() != shape_)
            {
                THROW_INVALID_ARGUMENT("input mask must be the same shape as the array it is masking.");
            }

            return put(mask.flatnonzero(), value);
        }

        //============================================================================
        // Method Description:
        /// Set the mask indices to the input values.
        ///
        /// @param mask
        /// @param values
        ///
        Matrix<dtype> &put_mask(const Matrix<bool> &mask, const Matrix<dtype> &values)
        {
            if (mask.shape() != shape_)
            {
                THROW_INVALID_ARGUMENT("input mask must be the same shape as the array it is masking.");
            }

            return put(mask.flatnonzero(), values);
        }

        //============================================================================
        // Method Description:
        /// Flattens the array but does not make a copy.
        ///
        /// @return Matrix
        ///
        Matrix<dtype> &ravel() noexcept
        {
            reshape(size_);
            return *this;
        }

        //============================================================================
        // Method Description:
        /// Repeat elements of an array.
        ///
        /// @param rows
        /// @param columns
        /// @return Matrix
        ///
        Matrix<dtype> repeat(size_t rows, size_t columns) const
        {
            Matrix<dtype> returnArray(shape_.rows * rows, shape_.cols * columns);

            for (size_t row = 0; row < rows; ++row)
            {
                for (size_t col = 0; col < columns; ++col)
                {
                    std::vector<size_t> indices(shape_.size());

                    const size_t rowStart = row * shape_.rows;
                    const size_t colStart = col * shape_.cols;

                    const size_t rowEnd = (row + 1) * shape_.rows;
                    const size_t colEnd = (col + 1) * shape_.cols;

                    size_t counter = 0;
                    for (size_t rowIdx = rowStart; rowIdx < rowEnd; ++rowIdx)
                    {
                        for (size_t colIdx = colStart; colIdx < colEnd; ++colIdx)
                        {
                            indices[counter++] = rowIdx * returnArray.shape_.cols + colIdx;
                        }
                    }

                    returnArray.put(Matrix<size_t>(indices), *this);
                }
            }

            return returnArray;
        }

        //============================================================================
        // Method Description:
        /// Repeat elements of an array.
        ///
        /// @param repeatShape
        /// @return Matrix
        ///
        Matrix<dtype> repeat(const Shape &repeatShape) const
        {
            return repeat(repeatShape.rows, repeatShape.cols);
        }

        //============================================================================
        // Method Description:
        /// Replaces a value of the array with another value
        ///
        /// @param oldValue: the value to replace
        /// @param newValue: the value to replace with
        ///
        inline void replace(value_type oldValue, value_type newValue)
        {
            if (std::isnan(oldValue))
                std::replace_if(std::execution::par_unseq, begin(), end(), [](double v) { return std::isnan(v); }, newValue);
            else
                std::replace(std::execution::par_unseq, begin(), end(), oldValue, newValue);
        }

        //============================================================================
        // Method Description:
        /// The new shape should be compatible with the original shape. If an single integer,
        /// then the result will be a 1-D array of that length. One shape dimension
        /// can be -1. In this case, the value is inferred from the length of the
        /// array and remaining dimensions.
        ///
        /// @param size
        ///
        Matrix<dtype> &reshape(size_t size)
        {
            if (size != size_)
            {
                std::string errStr = "Cannot reshape array of size " + std::to_string(size_) + " into shape ";
                errStr += "[" + std::to_string(1) + ", " + std::to_string(size) + "]";
                THROW_RUNTIME_ERROR(errStr);
            }

            shape_.rows = 1;
            shape_.cols = size;

            return *this;
        }

        //============================================================================
        // Method Description:
        /// The new shape should be compatible with the original shape. If an single integer,
        /// then the result will be a 1-D array of that length. One shape dimension
        /// can be -1. In this case, the value is inferred from the length of the
        /// array and remaining dimensions.
        ///
        /// @param rows
        /// @param columns
        ///
        Matrix<dtype> &reshape(int64_t rows, int64_t columns)
        {
            if (rows < 0)
            {
                if (size_ % columns == 0)
                {
                    return reshape(size_ / columns, columns);
                }

                std::string errStr = "Cannot reshape array of size " + std::to_string(size_) + " into a shape ";
                errStr += "with " + std::to_string(columns) + " columns";
                THROW_INVALID_ARGUMENT(errStr);
            }

            if (columns < 0)
            {
                if (size_ % rows == 0)
                {
                    return reshape(rows, size_ / rows);
                }

                std::string errStr = "Cannot reshape array of size " + std::to_string(size_) + " into a shape ";
                errStr += "with " + std::to_string(rows) + " rows";
                THROW_INVALID_ARGUMENT(errStr);
            }

            if ((rows * columns) != size_)
            {
                std::string errStr = "Cannot reshape array of size " + std::to_string(size_) + " into shape ";
                errStr += "[" + std::to_string(rows) + ", " + std::to_string(columns) + "]";
                THROW_INVALID_ARGUMENT(errStr);
            }

            shape_.rows = (rows);
            shape_.cols = (columns);

            return *this;
        }

        //============================================================================
        // Method Description:
        /// The new shape should be compatible with the original shape. If an single integer,
        /// then the result will be a 1-D array of that length. One shape dimension
        /// can be -1. In this case, the value is inferred from the length of the
        /// array and remaining dimensions.
        ///
        /// @param shape
        ///
        Matrix<dtype> &reshape(const Shape &shape)
        {
            return reshape(shape.rows, shape.cols);
        }

        //============================================================================
        // Method Description:
        /// Change shape and size of array in-place. All previous
        /// data of the array is lost.
        ///
        /// @param rows
        /// @param columns
        ///
        Matrix<dtype> &resize_fast(size_t rows, size_t columns)
        {
            new_matrix(Shape(rows, columns));
            return *this;
        }

        //============================================================================
        // Method Description:
        /// Change shape and size of array in-place. All previous
        /// data of the array is lost.
        ///
        /// @param shape
        ///
        Matrix<dtype> &resize_fast(const Shape &shape)
        {
            return resize_fast(shape.rows, shape.cols);
        }

        //============================================================================
        // Method Description:
        /// Return a new array with the specified shape. If new shape
        /// is larger than old shape then array will be padded with zeros.
        /// If new shape is smaller than the old shape then the data will
        /// be discarded.
        ///
        /// @param rows
        /// @param columns
        ///
        Matrix<dtype> &resize_slow(size_t rows, size_t columns)
        {
            std::vector<dtype> oldData(size_);
            std::copy(std::execution::par_unseq, begin(), end(), oldData.begin());

            const Shape shape(rows, columns);
            const Shape oldShape = shape_;

            new_matrix(shape);

            for (size_t row = 0; row < shape.rows; ++row)
            {
                for (size_t col = 0; col < shape.cols; ++col)
                {
                    if (row >= oldShape.rows || col >= oldShape.cols)
                    {
                        operator()(row, col) = dtype{0}; // zero fill
                    }
                    else
                    {
                        operator()(row, col) = oldData[row * oldShape.cols + col];
                    }
                }
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// Return a new array with the specified shape. If new shape
        /// is larger than old shape then array will be padded with zeros.
        /// If new shape is smaller than the old shape then the data will
        /// be discarded.
        ///
        /// @param shape
        ///
        Matrix<dtype> &resize_slow(const Shape &shape)
        {
            return resize_slow(shape.rows, shape.cols);
        }

        //============================================================================
        // Method Description:
        /// Return a with each element rounded to the given number
        /// of decimals.
        ///
        /// @param number_decimals (default 0)
        /// @return Matrix
        ///
        Matrix<dtype> round(uint8_t number_decimals = 0) const
        {
            STATIC_ASSERT_FLOAT(dtype);

            Matrix<dtype> returnArray(shape_);
            const double multFactor = std::pow(10., number_decimals);
            const auto function = [multFactor](dtype value) noexcept -> dtype
            { return static_cast<dtype>(std::nearbyint(static_cast<double>(value) * multFactor) / multFactor); };

            std::transform(std::execution::par_unseq, cbegin(), cend(), returnArray.begin(), function);

            return returnArray;
        }

        //============================================================================
        // Method Description:
        /// Returns the full row of the array
        ///
        ///
        /// @return Shape
        ///
        Matrix<dtype> row(size_t row)
        {
            return Matrix<dtype>(cbegin(row), cend(row));
        }

        //============================================================================
        // Method Description:
        /// Return the shape of the array
        ///
        /// @return Shape
        ///
        Shape shape() const noexcept
        {
            return shape_;
        }

        //============================================================================
        // Method Description:
        /// Return the size of the array
        ///
        /// @return size
        ///
        size_t size() const noexcept
        {
            return size_;
        }

        //============================================================================
        // Method Description:
        /// Sort an array, in-place.
        ///
        /// @param axis (Optional, default NONE)
        /// @return size
        ///
        Matrix<dtype> &sort(Axis axis = Axis::NONE)
        {

            const auto comparator = [](dtype lhs, dtype rhs) noexcept -> bool
            { return lhs < rhs; };

            switch (axis)
            {
            case Axis::NONE:
            {
                std::sort(std::execution::par_unseq, begin(), end(), comparator);
                break;
            }
            case Axis::COLUMN:
            {
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    std::sort(std::execution::par_unseq, begin(row), end(row), comparator);
                }
                break;
            }
            case Axis::ROW:
            {
                Matrix<dtype> transposedArray = transpose();
                for (size_t row = 0; row < transposedArray.shape_.rows; ++row)
                {
                    std::sort(std::execution::par_unseq, transposedArray.begin(row), transposedArray.end(row), comparator);
                }

                *this = transposedArray.transpose();
                break;
            }
            }

            return *this;
        }

        //============================================================================
        // Method Description:
        /// returns the Matrix as a string representation
        ///
        /// @return string
        ///
        std::string str() const
        {

            std::string out;
            out += "[\n";
            for (size_t row = 0; row < shape_.rows; ++row)
            {
                out += "[";
                // for (size_t col = 0; col < shape_.cols; ++col)
                // {
                //     out += utils::value2str(operator()(row, col)) + ", ";
                // }

                auto ib = begin() + row * shape_.cols;
                auto ie = ib + shape_.cols;
                out += str_join(ib, ie, ", ");

                if (row == shape_.rows - 1)
                {
                    out += "]\n";
                }
                else
                {
                    out += "],\n";
                }
            }
            out += "]\n";
            return out;
        }

        //============================================================================
        // Method Description:
        /// Return the sum of the array elements over the given axis.
        ///
        /// @param axis (Optional, default NONE)
        /// @return Matrix
        ///
        Matrix<dtype> sum(Axis axis = Axis::NONE) const
        {

            switch (axis)
            {
            case Axis::NONE:
            {
                Matrix<dtype> returnArray = {std::accumulate(cbegin(), cend(), dtype{0})};
                return returnArray;
            }
            case Axis::COLUMN:
            {
                Matrix<dtype> returnArray(1, shape_.rows);
                for (size_t row = 0; row < shape_.rows; ++row)
                {
                    returnArray(0, row) = std::accumulate(cbegin(row), cend(row), dtype{0});
                }

                return returnArray;
            }
            case Axis::ROW:
            {
                return transpose().sum(Axis::COLUMN);
            }
            default:
            {
                THROW_INVALID_ARGUMENT("Unimplemented axis type.");
                return {}; // get rid of compiler warning
            }
            }
        }

        //============================================================================
        // Method Description:
        /// Interchange two axes of an array. Equivalent to transpose...
        ///
        /// @return Matrix
        ///
        Matrix<dtype> swap_axes() const
        {
            return transpose();
        }

        //============================================================================
        // Method Description:
        /// Swaps rows of the array
        ///
        /// @param colIdx1
        /// @param colIdx2
        ///
        void swap_columns(size_t colIdx1, size_t colIdx2)
        {
            for (size_t row = 0; row < static_cast<size_t>(shape_.rows); ++row)
            {
                std::swap(operator()(row, colIdx1), operator()(row, colIdx2));
            }
        }

        //============================================================================
        // Method Description:
        /// Swaps rows of the array
        ///
        /// @param rowIdx1
        /// @param rowIdx2
        ///
        void swap_rows(size_t rowIdx1, size_t rowIdx2)
        {
            for (size_t col = 0; col < static_cast<size_t>(shape_.cols); ++col)
            {
                std::swap(operator()(rowIdx1, col), operator()(rowIdx2, col));
            }
        }

        //============================================================================
        // Method Description:
        /// Write array to a file as binary.
        /// The data produced by this method can be recovered
        /// using the function fromfile().
        ///
        /// @param filename
        /// @return None
        ///
        void tofile(const std::string &filename) const
        {
            dump(filename);
        }

        //============================================================================
        // Method Description:
        /// Write array to a file as text.
        /// The data produced by this method can be recovered
        /// using the function fromfile().
        ///
        /// @param filename
        /// @param separator: Separator between array items for text output.
        /// @return None
        ///
        void tofile(const std::string &filename, const char separator) const
        {

            filesystem::File f(filename);
            if (!f.hasExt())
            {
                f.withExt("txt");
            }

            std::ofstream ofile(f.fullname().c_str());
            if (!ofile.good())
            {
                THROW_RUNTIME_ERROR("Input file could not be opened:\n\t" + filename);
            }

            size_t counter = 0;
            for (auto value : *this)
            {
                ofile << value;
                if (counter++ != size_ - 1)
                {
                    ofile << separator;
                }
            }
            ofile.close();
        }

        //============================================================================
        // Method Description:
        /// Converts the slice object to an Matrix of indices for this array
        ///
        /// @param slice: the slice object
        /// @param axis: the array axis
        ///
        /// @return Matrix<size_t>
        ///
        Matrix<size_t> to_indices(Slice slice, Axis axis = Axis::ROW) const
        {
            size_t element_number = 0;
            switch (axis)
            {
            case Axis::NONE:
            {
                element_number = slice.element_number(size_);
                break;
            }
            case Axis::ROW:
            {
                element_number = slice.element_number(shape_.rows);
                break;
            }
            case Axis::COLUMN:
            {
                element_number = slice.element_number(shape_.cols);
                break;
            }
            default:
            {
                // not actually possible, getting rid of compiler warning
                THROW_INVALID_ARGUMENT("Invalid 'axis' option");
            }
            }

            if (element_number == 0)
            {
                return {};
            }

            Matrix<size_t> indices(1, element_number);
            indices[0] = slice.start;
            for (size_t i = 1; i < indices.size(); ++i)
            {
                indices[i] = (indices[i - 1] + slice.step);
            }

            return indices;
        }

        //============================================================================
        // Method Description:
        /// Write flattened array to an STL vector
        ///
        /// @return std::vector
        ///
        std::vector<dtype> to_vector() const
        {
            return std::vector<dtype>(cbegin(), cend());
        }

        //============================================================================
        // Method Description:
        /// Return the sum along diagonals of the array.
        ///
        /// @param offset: Offset of the diagonal from the main diagonal. Can be both positive and negative. Defaults
        /// to 0.
        /// @param axis: (Optional, default ROW) Axis to offset from
        ///
        /// @return value
        ///
        value_type trace(size_t offset = 0, Axis axis = Axis::ROW) const
        {

            size_t rowStart = 0;
            size_t colStart = 0;
            switch (axis)
            {
            case Axis::ROW:
            {
                rowStart += offset;
                break;
            }
            case Axis::COLUMN:
            {
                colStart += offset;
                break;
            }
            default:
            {
                // if the user input NONE, override back to ROW
                axis = Axis::ROW;
                break;
            }
            }

            if (rowStart >= shape_.rows || colStart >= shape_.cols)
            {
                return dtype{0};
            }

            size_t col = colStart;
            dtype sum = 0;
            for (size_t row = rowStart; row < shape_.rows; ++row)
            {
                if (col >= shape_.cols)
                {
                    break;
                }
                sum += operator()(row, col++);
            }

            return sum;
        }

        //============================================================================
        // Method Description:
        /// Transpose the rows and columns of an array
        ///
        /// @return Matrix
        ///
        Matrix<dtype> transpose() const
        {
            Matrix<dtype> transArray(shape_.cols, shape_.rows);
            for (size_t row = 0; row < shape_.rows; ++row)
            {
                for (size_t col = 0; col < shape_.cols; ++col)
                {
                    transArray(col, row) = operator()(row, col);
                }
            }
            return transArray;
        }

        //============================================================================
        // Method Description:
        /// Fills the array with zeros
        ///
        ///
        Matrix<dtype> &zeros() noexcept
        {

            fill(dtype{0});
            return *this;
        }

        //============================================================================
        // Method Description:
        /// Return random integer from low (inclusive) to high (exclusive),
        /// with the given shape. If no high value is input then the range will
        /// go from [0, 100).
        ///
        /// @param begin position for filling area, default 0
        /// @param end position for filling area, default size
        /// @param low  default 0
        /// @param high default 100.
        /// @return NdArray
        ///
        Matrix<dtype> &fill_random(size_t begin, size_t end, dtype low = 0, dtype high = 100)
        {
            STATIC_ASSERT_INTEGER(dtype);

            if (low == high)
            {
                THROW_INVALID_ARGUMENT("input low value must be less than the input high value.");
            }
            else if (low > high)
            {
                std::swap(low, high);
            }

            std::uniform_int_distribution<dtype> dist(low, high - 1);
            std::mt19937_64 generator{};
            auto b = this->begin() + begin;
            auto e = this->begin() + end;
            auto f = [&dist, &generator](dtype &value) noexcept -> void
            { value = int(dist(generator)); };
            std::for_each(std::execution::par_unseq, b, e, f);

            return *this;
        }

    private:
        //====================================Attributes==============================
        allocator_type allocator_{};
        Shape shape_{0, 0};
        size_t size_{0};
        Endian endianess_{Endian::NATIVE};
        pointer array_{nullptr};
        bool ownsPtr_{false};

        //============================================================================
        // Method Description:
        /// Return iterator of specials index
        ///
        /// @param rowPos
        /// @param columnPos
        /// @return iterator
        ///
        iterator pos_iterator(int64_t rowPos, int64_t columnPos = 0)
        {
            // this doesn't allow for calling the first element as -size_...
            // but why would you really want to do that anyway?
            if (std::abs(rowPos) >= shape_.rows)
            {
                std::string errStr = "Row index " + std::to_string(rowPos);
                errStr += " is out of bounds for array of size " + std::to_string(shape_.rows) + ".";
                THROW_INVALID_ARGUMENT(errStr);
            }

            // this doesn't allow for calling the first element as -size_...
            // but why would you really want to that anyway?
            if (std::abs(columnPos) >= shape_.cols)
            {
                std::string errStr = "Column index " + std::to_string(columnPos);
                errStr += " is out of bounds for array of size " + std::to_string(shape_.cols) + ".";
                THROW_INVALID_ARGUMENT(errStr);
            }

            rowPos = rowPos >= 0 ? rowPos : rowPos + shape_.rows;
            columnPos = columnPos >= 0 ? columnPos : shape_.cols + columns();

            iterator out_iter = begin() + rowPos * columns() + columnPos;

            return out_iter;
        }

        //============================================================================
        // Method Description:
        /// Deletes the internal array
        ///
        void clear() noexcept
        {
            if (ownsPtr_ && array_ != nullptr)
            {
                allocator_.deallocate(array_, size_);
            }

            array_ = nullptr;
            shape_.rows = shape_.cols = 0;
            size_ = 0;
            ownsPtr_ = false;
            endianess_ = Endian::NATIVE;
        }

        //============================================================================
        // Method Description:
        /// Creates a new internal array
        ///
        void new_matrix()
        {
            if (size_ > 0)
            {
                array_ = allocator_.allocate(size_);
                ownsPtr_ = true;
            }
        }

        //============================================================================
        // Method Description:
        /// Creates a new internal array
        ///
        /// @param shape
        ///
        void new_matrix(const Shape &shape)
        {
            clear();

            shape_ = shape;
            size_ = shape.size();
            new_matrix();
        }

        friend std::ostream &operator<<(std::ostream &stream, const self_type &matrix)
        {
            stream << matrix.str();
            return stream;
        }
    };

    // NOTE: this needs to be defined outside of the class to get rid of a compiler
    // error in Visual Studio
    template <typename dtype, class _Alloc>
    AIS_EXPORT std::pair<Matrix<size_t>, Matrix<size_t>> Matrix<dtype, _Alloc>::nonzero() const
    {

        std::vector<size_t> rowIndices;
        std::vector<size_t> columnIndices;

        for (size_t row = 0; row < shape_.rows; ++row)
        {
            for (size_t col = 0; col < shape_.cols; ++col)
            {
                auto d = operator()(row, col);
                if (d != dtype{0})
                {
                    rowIndices.push_back(row);
                    columnIndices.push_back(col);
                }
            }
        }

        return std::make_pair(Matrix<size_t>(rowIndices), Matrix<size_t>(columnIndices));
    };

    template <typename dtype, class _Alloc>
    AIS_EXPORT std::pair<Matrix<size_t>, Matrix<size_t>> Matrix<dtype, _Alloc>::nonempty() const
    {

        std::vector<size_t> rowIndices;
        std::vector<size_t> columnIndices;

        for (size_t row = 0; row < shape_.rows; ++row)
        {
            for (size_t col = 0; col < shape_.cols; ++col)
            {
                auto d = operator()(row, col);
                if (!ais::isnan(d))
                {
                    rowIndices.push_back(row);
                    columnIndices.push_back(col);
                }
            }
        }

        return std::make_pair(Matrix<size_t>(rowIndices), Matrix<size_t>(columnIndices));
    };

    template <typename dtype, class _Alloc>
    AIS_EXPORT std::pair<Matrix<size_t>, Matrix<size_t>> Matrix<dtype, _Alloc>::data_area() const
    {

        std::vector<size_t> rowIndices;
        std::vector<size_t> columnIndices;

        for (size_t row = 0; row < shape_.rows; ++row)
        {
            for (size_t col = 0; col < shape_.cols; ++col)
            {
                auto d = operator()(row, col);
                if (!ais::isnan(d))
                {
                    rowIndices.push_back(row);
                    columnIndices.push_back(col);
                }
            }
        }

        return std::make_pair(Matrix<size_t>(rowIndices), Matrix<size_t>(columnIndices));
    }

    using FMatrix = AIS_EXPORT Matrix<double>;
} // namespace ais