Matrix.java

/*
 * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package com.irurueta.algebra;

import com.irurueta.statistics.GaussianRandomizer;
import com.irurueta.statistics.UniformRandomizer;

import java.io.Serializable;
import java.util.Arrays;
import java.util.Objects;
import java.util.Random;

/**
 * Defines a matrix of numerical data.
 * Values of a matrix are stored inside an internal array of data.
 */
@SuppressWarnings("DuplicatedCode")
public class Matrix implements Serializable, Cloneable {
    /**
     * Constant defining the default order in which values are stored in a
     * matrix. This can be useful when copying data to and from an array.
     */
    public static final boolean DEFAULT_USE_COLUMN_ORDER = true;

    /**
     * Number of matrix rows.
     */
    private int rows;

    /**
     * Number of matrix columns.
     */
    private int columns;

    /**
     * Array containing data of matrix. Data is stored linearly in memory
     * using column order.
     */
    private double[] buffer;

    /**
     * Array used for indexing column start positions within buffer. This
     * is used for faster access to matrix elements.
     */
    private int[] columnIndex;

    /**
     * Constructor of this class.
     * By default, when instantiating a matrix its values are undefined.
     *
     * @param rows    Defines number of rows in matrix.
     * @param columns Defines number of columns in matrix.
     * @throws WrongSizeException Exception thrown if matrix must be
     *                            empty with provided size. In other words, matrices must have at least
     *                            one row and one column. Leaving number of rows or columns to zero will
     *                            raise this exception.
     */
    public Matrix(final int rows, final int columns) throws WrongSizeException {
        internalResize(rows, columns);
    }

    /**
     * Copy constructor.
     *
     * @param m matrix to copy from.
     */
    public Matrix(final Matrix m) {
        try {
            internalResize(m.getRows(), m.getColumns());
        } catch (final WrongSizeException ignore) {
            // never happens
        }

        System.arraycopy(m.getBuffer(), 0, buffer, 0, buffer.length);
    }

    /**
     * Returns number of rows in matrix.
     *
     * @return Number of rows in matrix.
     */
    public int getRows() {
        return rows;
    }

    /**
     * Returns number of columns in matrix.
     *
     * @return Number of columns in matrix.
     */
    public int getColumns() {
        return columns;
    }

    /**
     * Obtains element in matrix located at position (row and column).
     * Notice that row and column position are zero-indexed.
     *
     * @param row    Row to be used for element location.
     * @param column Column to be used for element location.
     * @return Value of element in matrix located at provided position.
     * @throws ArrayIndexOutOfBoundsException Exception raised if attempting
     *                                        to access a location that lies outside the boundaries of the internal
     *                                        array containing matrix data. Notice that internal data is stored in
     *                                        column order, hence, if row position exceeds the number of rows in the
     *                                        matrix, this exception might not be raised depending on column position,
     *                                        however, if the column value exceeds the number of columns, it will
     *                                        always raise the exception.
     */
    public double getElementAt(final int row, final int column) {
        return buffer[columnIndex[column] + row];
    }

    /**
     * Returns index within internal buffer corresponding to provided row
     * and column positions.
     *
     * @param row    row to be used for element location.
     * @param column column to be used for element location.
     * @return index within internal buffer.
     */
    public int getIndex(final int row, final int column) {
        return columnIndex[column] + row;
    }

    /**
     * Obtains element in matrix located at provided index value using the
     * order specified by DEFAULT_USE_COLUMN_ORDER.
     * Note: an index indicates linear position within matrix, in layman terms,
     * if using column order, then position (i, j), where i stands for row
     * index and j for column index, becomes index = j * rows + 1. On the other
     * hand, if using row order then index becomes index = i * columns + j.
     *
     * @param index Linear position.
     * @return Value of element contained at provided position.
     * @throws ArrayIndexOutOfBoundsException Exception raised if index lays
     *                                        outside valid values, which range from zero (inclusive) to
     *                                        rows * columns (exclusive)
     */
    public double getElementAtIndex(final int index) {
        return getElementAtIndex(index, DEFAULT_USE_COLUMN_ORDER);
    }

    /**
     * Obtain element in matrix located at provided index value using provided
     * order (either column or row order).
     *
     * @param index         Linear position.
     * @param isColumnOrder if true indicates that values are retrieved
     *                      assuming that are stored in column order.
     * @return Value of element contained at provided position.
     * @throws ArrayIndexOutOfBoundsException Exception raised if index lays
     *                                        outside valid values, which range from zero (inclusive) to
     *                                        rows * columns (exclusive)
     */
    public double getElementAtIndex(final int index, final boolean isColumnOrder) {

        if (isColumnOrder) {
            return buffer[index];
        } else {
            final var row = index / columns;
            final var column = index % columns;
            return buffer[columnIndex[column] + row];
        }
    }

    /**
     * Sets element in matrix located at provided position (row and column).
     *
     * @param row    Row to be used for element location to be set.
     * @param column Column to be used for element location to be set.
     * @param value  Value to be set at provided position.
     * @throws ArrayIndexOutOfBoundsException Exception raised if attempting
     *                                        to access a location that lies outside the boundaries of the internal
     *                                        array containing matrix data. Notice that internal data is stored in
     *                                        column order, hence, if row position exceeds the number of rows in the
     *                                        matrix, this exception might not be raised depending on column position,
     *                                        however, if the column value exceeds the number of columns, it will
     *                                        always raise the exception.
     */
    public void setElementAt(final int row, final int column, final double value) {
        buffer[columnIndex[column] + row] = value;
    }

    /**
     * Sets element in matrix located at provided index using the order
     * specified by DEFAULT_USE_COLUMN_ORDER.
     * Note: an index indicates linear position within matrix, in layman terms,
     * if using column order, then position (i, j), where i stands for row
     * index and j for column index, becomes index = j * rows + 1. On the other
     * hand, if using row order then index becomes index = i * columns + j.
     *
     * @param index Linear position.
     * @param value Value of element contained at provided position.
     * @throws ArrayIndexOutOfBoundsException Exception raised if index lays
     *                                        outside valid values, which range from zero (inclusive) to
     *                                        rows * columns (exclusive)
     */
    public void setElementAtIndex(final int index, final double value) {
        setElementAtIndex(index, value, DEFAULT_USE_COLUMN_ORDER);
    }

    /**
     * Sets element in matrix located at provided index using provided order
     * (either column or row order).
     *
     * @param index         Linear position.
     * @param value         Value of element to be set at provided position.
     * @param isColumnOrder if true indicates that values are retrieved
     *                      assuming that are stored in column order.
     * @throws ArrayIndexOutOfBoundsException Exception raised if index lays
     *                                        outside valid values, which range from zero (inclusive) to
     *                                        rows * columns (exclusive).
     */
    public void setElementAtIndex(final int index, final double value, final boolean isColumnOrder) {
        if (isColumnOrder) {
            buffer[index] = value;
        } else {
            final var row = index / columns;
            final var column = index % columns;
            buffer[columnIndex[column] + row] = value;
        }
    }

    /**
     * Returns a new matrix instance containing the same data as this instance.
     *
     * @return A copy of this matrix instance.
     * @throws CloneNotSupportedException if clone fails.
     */
    @Override
    public Matrix clone() throws CloneNotSupportedException {
        var out = (Matrix) super.clone();
        out.copyFrom(this);
        return out;
    }

    /**
     * Copies this matrix data into provided matrix. Provided output matrix will
     * be resized if needed.
     *
     * @param output Destination matrix where data will be copied to.
     * @throws NullPointerException Exception raised if provided output matrix
     *                              is null.
     */
    public void copyTo(final Matrix output) {
        // reset output size if not equal, otherwise reuse buffer and column
        // index
        if (output.getRows() != rows || output.getColumns() != columns) {
            // resets size and column index
            try {
                output.resize(rows, columns);
            } catch (final WrongSizeException ignore) {
                // never happens
            }
        }
        // copies content
        System.arraycopy(buffer, 0, output.buffer, 0, buffer.length);
    }

    /**
     * Copies the contents of provided matrix into this instance. This instance
     * will be resized if needed.
     *
     * @param input Input matrix where data will be copied from.
     * @throws NullPointerException Exception raised if provided input matrix is
     *                              null.
     */
    public void copyFrom(final Matrix input) {
        // reset size if not equal, otherwise reuse buffer and column index
        if (input.getRows() != rows || input.getColumns() != columns) {
            // resets size and column index
            try {
                resize(input.getRows(), input.getColumns());
            } catch (final WrongSizeException ignore) {
                // never happens
            }
        }
        // copies content
        System.arraycopy(input.buffer, 0, buffer, 0, buffer.length);
    }

    /**
     * Adds another matrix to this matrix instance and stores the result in
     * provided result matrix. If provided result matrix doesn't have proper
     * size, it will be resized.
     *
     * @param other  Matrix to be added to current instance.
     * @param result Matrix where result of summation is stored.
     * @throws WrongSizeException   Exception thrown if provided matrix to be
     *                              added (i.e. other) does not have the same size as this matrix.
     * @throws NullPointerException Exception raised if provided matrices are
     *                              null.
     */
    public void add(final Matrix other, final Matrix result) throws WrongSizeException {
        if (other.getRows() != rows || other.getColumns() != columns) {
            throw new WrongSizeException();
        }

        // resize result if needed
        if (result.getRows() != rows || result.getColumns() != columns) {
            result.resize(rows, columns);
        }

        internalAdd(other, result);
    }

    /**
     * Adds provided matrix to this instance and returns the result as a new
     * matrix instance.
     *
     * @param other Matrix to be added.
     * @return Returns a new matrix containing the sum of this matrix with
     * provided matrix.
     * @throws WrongSizeException   Exception raised if provided matrix does
     *                              not have the same size as this matrix.
     * @throws NullPointerException Exception raised if provided matrix is null.
     */
    public Matrix addAndReturnNew(final Matrix other) throws WrongSizeException {
        if (other.getRows() != rows || other.getColumns() != columns) {
            throw new WrongSizeException();
        }

        final var out = new Matrix(rows, columns);
        internalAdd(other, out);
        return out;
    }

    /**
     * Adds provided matrix to this instance.
     *
     * @param other Matrix to be added.
     * @throws WrongSizeException   Exception raised if provided matrix does
     *                              not have the same size as this matrix.
     * @throws NullPointerException Exception raised if provided matrix is null.
     */
    public void add(final Matrix other) throws WrongSizeException {
        if (other.getRows() != rows || other.getColumns() != columns) {
            throw new WrongSizeException();
        }

        internalAdd(other, this);
    }

    /**
     * Subtracts another matrix from this matrix instance and stores the result
     * in provided result matrix. If provided result matrix doesn't have proper
     * size, it will be resized.
     *
     * @param other  Matrix to be added to current instance.
     * @param result Matrix where result of subtraction is stored.
     * @throws WrongSizeException   Exception thrown if provided matrix to be
     *                              subtracted (i.e. other) does not have the same size as this matrix.
     * @throws NullPointerException Exception raised if provided matrices are
     *                              null.
     */
    public void subtract(final Matrix other, final Matrix result) throws WrongSizeException {
        if (other.getRows() != rows || other.getColumns() != columns) {
            throw new WrongSizeException();
        }

        // resize result if needed
        if (result.getRows() != rows || result.getColumns() != columns) {
            result.resize(rows, columns);
        }

        internalSubtract(other, result);
    }

    /**
     * Subtracts provided matrix from this instance and returns the result as a
     * new matrix instance.
     *
     * @param other Matrix to be subtracted from.
     * @return Returns a new matrix containing the subtraction of provided
     * matrix from this matrix.
     * @throws WrongSizeException   Exception raised if provided matrix does
     *                              not have the same size as this matrix.
     * @throws NullPointerException Exception raised if provided matrix is null.
     */
    public Matrix subtractAndReturnNew(final Matrix other) throws WrongSizeException {
        if (other.getRows() != rows || other.getColumns() != columns) {
            throw new WrongSizeException();
        }

        final var out = new Matrix(rows, columns);
        internalSubtract(other, out);
        return out;
    }

    /**
     * Subtracts provided matrix from this instance.
     *
     * @param other Matrix to be subtracted from.
     * @throws WrongSizeException   Exception raised if provided matrix does
     *                              not have the same size as this matrix.
     * @throws NullPointerException Exception raised if provided matrix is null.
     */
    public void subtract(final Matrix other) throws WrongSizeException {
        if (other.getRows() != rows || other.getColumns() != columns) {
            throw new WrongSizeException();
        }

        internalSubtract(other, this);
    }

    /**
     * Multiplies another matrix to this matrix instance and stores the result
     * in provided result matrix. If provided result matrix doesn't have proper
     * size, it will be resized.
     *
     * @param other  Matrix to be multiplied to current instance.
     * @param result Matrix where result of product is stored.
     * @throws WrongSizeException   Exception thrown when current and provided
     *                              matrix (i.e. other) has incompatible size for product computation.
     * @throws NullPointerException Exception raised if provided matrices are
     *                              null.
     */
    public void multiply(final Matrix other, final Matrix result) throws WrongSizeException {
        if (columns != other.rows) {
            throw new WrongSizeException();
        }

        // resize result if needed
        if (result.rows != rows || result.columns != other.columns) {
            result.resize(rows, other.columns);
        }

        internalMultiply(other, result);
    }

    /**
     * Multiplies this matrix with provided matrix and returns the result as
     * a new instance.
     * If this matrix m1 has size m x n and provided matrix m2 has size p x q,
     * then n must be equal to p so that product m1 * m2 can be correctly
     * computed obtaining a matrix of size m x q, otherwise an
     * IllegalArgumentException will be raised.
     *
     * @param other Right operand of matrix product
     * @return Matrix containing result of multiplication
     * @throws WrongSizeException   Exception thrown when current and
     *                              provided matrices have incompatible sizes for product computation.
     * @throws NullPointerException Exception thrown if provided matrix is null
     */
    public Matrix multiplyAndReturnNew(final Matrix other) throws WrongSizeException {

        if (columns != other.rows) {
            throw new WrongSizeException();
        }

        final var out = new Matrix(rows, other.columns);
        internalMultiply(other, out);
        return out;
    }

    /**
     * Multiplies this matrix with provided matrix.
     * If this matrix m1 has size m x n and provided matrix m2 has size p x q,
     * then n must be equal to p so that product m1 * m2 can be correctly
     * computed resizing this matrix to a new one having size m x q, otherwise
     * an IllegalArgumentException will be raised.
     *
     * @param other Right operand of matrix product
     * @throws WrongSizeException   Exception thrown when current and
     *                              provided matrices have incompatible sizes for product computation.
     * @throws NullPointerException Exception thrown if provided matrix is null
     */
    public void multiply(final Matrix other) throws WrongSizeException {

        if (columns != other.rows) {
            throw new WrongSizeException();
        }

        // instantiate new buffer and column index
        final var resultBuffer = new double[rows * other.columns];
        final var resultColumnIndex = new int[other.columns];
        var counter = 0;
        for (int i = 0; i < other.columns; i++) {
            resultColumnIndex[i] = counter;
            counter += rows;
        }

        internalMultiply(other, resultBuffer, resultColumnIndex);
        // update matrix data
        columns = other.columns;
        columnIndex = resultColumnIndex;
        buffer = resultBuffer;
    }

    /**
     * Computes the Kronecker product with provided matrix and stores the
     * result in provided result matrix. If provided result matrix doesn't
     * have proper size, it will be resized.
     *
     * @param other  other matrix to be Kronecker multiplied to current matrix.
     * @param result matrix where result will be stored.
     */
    public void multiplyKronecker(final Matrix other, final Matrix result) {
        // resize result if needed
        final var resultRows = rows * other.rows;
        final var resultCols = columns * other.columns;
        if (result.rows != resultRows || result.columns != resultCols) {
            try {
                result.resize(resultRows, resultCols);
            } catch (final WrongSizeException ignore) {
                // never thrown
            }
        }

        internalMultiplyKronecker(other, result);
    }

    /**
     * Computes the Kronecker product with provided matrix and returns the
     * result as a new instance.
     * If this matrix m1 has size mxn and provided matrix m2 has size pxq,
     * the resulting matrix will be m*pxn*q.
     *
     * @param other other matrix to be Kronecker multiplied to current matrix.
     * @return matrix containing result of Kronecker multiplication.
     */
    public Matrix multiplyKroneckerAndReturnNew(final Matrix other) {
        Matrix out = null;
        try {
            out = new Matrix(rows * other.rows, columns * other.columns);
            internalMultiplyKronecker(other, out);
        } catch (final WrongSizeException ignore) {
            // never thrown
        }
        return out;
    }

    /**
     * Computes the Kronecker product of this matrix with provided matrix and
     * updates this matrix with the result of the multiplication.
     * If this matrix m1 has size mxn and provided matrix m2 has size pxq,
     * the resulting matrix will be resized to m*pxn*q.
     *
     * @param other other matrix to be Kronecker multiplied to current matrix.
     */
    public void multiplyKronecker(final Matrix other) {

        // instantiate new buffer and column index
        final var resultRows = rows * other.rows;
        final var resultCols = columns * other.columns;
        final var resultBuffer = new double[resultRows * resultCols];
        final var resultColumnIndex = new int[resultCols];
        var counter = 0;
        for (int i = 0; i < resultCols; i++) {
            resultColumnIndex[i] = counter;
            counter += resultRows;
        }

        internalMultiplyKronecker(other, resultBuffer, resultColumnIndex);
        // update matrix data
        rows = resultRows;
        columns = resultCols;
        columnIndex = resultColumnIndex;
        buffer = resultBuffer;
    }

    /**
     * Computes product by scalar of this instance multiplying all its elements
     * by provided scalar value and storing the results in provided result
     * matrix. If provided result matrix doesn't have proper size, it will be
     * automatically resized.
     *
     * @param scalar Scalar amount that current matrix will be multiplied by.
     * @param result Matrix where result of operation is stored.
     */
    private void multiplyByScalar(final double scalar, final Matrix result) {
        final var length = rows * columns;
        for (int i = 0; i < length; i++) {
            result.buffer[i] = scalar * buffer[i];
        }
    }

    /**
     * Computes product by scalar of this instance multiplying all its elements
     * by provided scalar value and returning the result as a new instance.
     *
     * @param scalar Scalar amount that current matrix will be multiplied by.
     * @return Returns a new matrix instance that contains result of product by
     * scalar.
     */
    public Matrix multiplyByScalarAndReturnNew(final double scalar) {
        Matrix out = null;
        try {
            out = new Matrix(rows, columns);
            multiplyByScalar(scalar, out);
        } catch (final WrongSizeException ignore) {
            // never happens
        }
        return out;
    }

    /**
     * Computes product by scalar of this instance multiplying all its elements
     * by provided scalar value and returning the result as a new instance.
     *
     * @param scalar Scalar amount that current matrix will be multiplied by.
     */
    public void multiplyByScalar(final double scalar) {
        multiplyByScalar(scalar, this);
    }

    /**
     * Checks if provided object is a Matrix instance having exactly the same
     * contents as this matrix instance.
     *
     * @param obj Object to be compared
     * @return Returns true if both objects are considered to be equal.
     */
    @Override
    public boolean equals(final Object obj) {
        if (obj == null) {
            return false;
        }
        if (obj == this) {
            return true;
        }
        if (!(obj instanceof Matrix other)) {
            return false;
        }

        return equals(other);
    }

    /**
     * Computes and returns hash code for this instance. Hash codes are almost
     * unique values that are useful for fast classification and storage of
     * objects in collections.
     *
     * @return Hash code
     */
    @Override
    public int hashCode() {
        return Objects.hash(this.rows, this.columns, Arrays.hashCode(this.buffer), Arrays.hashCode(this.columnIndex));
    }

    /**
     * Checks if provided matrix has exactly the same contents as this matrix
     * instance.
     *
     * @param other Matrix to be compared.
     * @return Returns true if both objects are considered to be equal (same
     * content and size)
     */
    public boolean equals(final Matrix other) {
        return equals(other, 0.0);
    }

    /**
     * Checks if provided matrix has contents similar to this matrix by checking
     * that all values have a maximum difference equal to provided threshold and
     * same size.
     *
     * @param other     Matrix to be compared
     * @param threshold Maximum difference allowed between values on same
     *                  position to determine that matrices are equal
     * @return True if matrices are considered to be equal (almost equal content
     * and same size)
     */
    public boolean equals(final Matrix other, final double threshold) {
        if (other == null) {
            return false;
        }
        if (other.getRows() != rows) {
            return false;
        }
        if (other.getColumns() != columns) {
            return false;
        }

        // check contents
        final var length = rows * columns;
        for (int i = 0; i < length; i++) {
            if (Math.abs(buffer[i] - other.buffer[i]) > threshold) {
                return false;
            }
        }
        return true;
    }

    /**
     * Computes element by element product (i.e. Hadamard product) between
     * current and provided (i.e. other) matrix, and stores the result in
     * provided result matrix.
     *
     * @param other  Instance that will be used for element by element product
     *               with current instance.
     * @param result Matrix where result of operation is stored.
     * @throws WrongSizeException   Exception raised if attempting to perform
     *                              element by element product on matrices of different size.
     * @throws NullPointerException Exception raised if provided matrix is null
     */
    public void elementByElementProduct(final Matrix other, final Matrix result) throws WrongSizeException {
        if (other.getRows() != rows || other.getColumns() != columns) {
            throw new WrongSizeException();
        }

        // resize result if needed
        if (result.getRows() != rows || result.getColumns() != columns) {
            result.resize(rows, columns);
        }

        internalElementByElementProduct(other, result);
    }

    /**
     * Computes element by element product (i.e. Hadamard product) between
     * current and provided instances and returns the result as a new instance.
     * Example:
     * Having matrices
     * [1, 2, 3]		 [10, 11, 12]
     * m1 =	[4, 5, 6]	m2 = [13, 14, 15]
     * [7, 8, 9]		 [16, 17, 18]
     * Then their element by element product will be
     * [10,  22,  36 ]
     * m3 =    [52,  70,  90 ]
     * [112, 136, 162]
     * <p>
     * Note: Attempting to perform element product on matrices of different size
     * will raise an IllegalArgumentException
     *
     * @param other Instance that will be used for element by element product
     *              with current instance.
     * @return A new Matrix containing element by element product result
     * @throws WrongSizeException   Exception raised if attempting to
     *                              perform element by element product on matrices of different size.
     * @throws NullPointerException Exception raised if provided matrix is null
     */
    public Matrix elementByElementProductAndReturnNew(final Matrix other) throws WrongSizeException {
        if (other.getRows() != rows || other.getColumns() != columns) {
            throw new WrongSizeException();
        }

        final var out = new Matrix(rows, columns);
        internalElementByElementProduct(other, out);
        return out;
    }

    /**
     * Computes element by element product (i.e. Hadamard product) between
     * current and provided instances.
     * Example:
     * Having matrices
     * [1, 2, 3]		 [10, 11, 12]
     * m1 =	[4, 5, 6]	m2 = [13, 14, 15]
     * [7, 8, 9]		 [16, 17, 18]
     * Then their element by element product will be
     * [10,  22,  36 ]
     * m3 =    [52,  70,  90 ]
     * [112, 136, 162]
     * <p>
     * Note: Attempting to perform element product on matrices of different size
     * will raise an IllegalArgumentException
     *
     * @param other Instance that will be used for element by element product
     *              with current instance.
     * @throws WrongSizeException   Exception raised if attempting to
     *                              perform element by element product on matrices of different size.
     * @throws NullPointerException Exception raised if provided matrix is null
     */
    public void elementByElementProduct(final Matrix other) throws WrongSizeException {
        if (other.getRows() != rows || other.getColumns() != columns) {
            throw new WrongSizeException();
        }

        internalElementByElementProduct(other, this);
    }

    /**
     * Transposes current matrix and stores result in provided matrix. If
     * provided matrix doesn't have proper size, it will be resized.
     *
     * @param result Instance where transposed matrix is stored.
     */
    public void transpose(final Matrix result) {
        // resize result if needed
        if (result.getRows() != columns || result.getColumns() != rows) {
            try {
                result.resize(columns, rows);
            } catch (final WrongSizeException ignore) {
                // never happens
            }
        }
        internalTranspose(result);
    }

    /**
     * Transposes current matrix and returns result as a new instance.
     * Transposition of a matrix is done by exchanging rows and columns, in
     * Layman terms, given a matrix m1 with elements located at m1(i,j), where i
     * is the row index and j is the column index, then it follows that its
     * transposed matrix m2 has the following property m2(i,j) = m1(j, i).
     *
     * @return A new Matrix instance containing transposed matrix.
     */
    public Matrix transposeAndReturnNew() {
        Matrix out = null;
        try {
            out = new Matrix(columns, rows);
            internalTranspose(out);
        } catch (final WrongSizeException ignore) {
            // never happens
        }
        return out;
    }

    /**
     * Transposes current matrix.
     * Transposition of a matrix is done by exchanging rows and columns, in
     * Layman terms, given a matrix m1 with elements located at m1(i,j), where i
     * is the row index and j is the column index, then it follows that its
     * transposed matrix m2 has the following property m2(i,j) = m1(j, i).
     */
    public void transpose() {

        final var newBuffer = new double[rows * columns];
        final var newColumnIndex = new int[rows];
        var counter = 0;
        for (var i = 0; i < rows; i++) {
            newColumnIndex[i] = counter;
            counter += columns;
        }
        internalTranspose(newBuffer, newColumnIndex);

        // update matrix data

        // swap rows and columns
        final var tmp = rows;
        rows = columns;
        columns = tmp;

        columnIndex = newColumnIndex;
        buffer = newBuffer;
    }

    /**
     * Sets the contents of this matrix to provided value in all of its elements
     *
     * @param initValue Value to be set on all the elements of this matrix.
     */
    public void initialize(final double initValue) {
        // initialize buffer array to provided value
        Arrays.fill(buffer, initValue);
    }

    /**
     * Resizes current instance by removing its contents and resizing it to
     * provided size.
     *
     * @param rows    Number of rows to be set
     * @param columns Number of columns to be set
     * @throws WrongSizeException Exception raised if either rows or
     *                            columns is zero.
     */
    public void resize(final int rows, final int columns) throws WrongSizeException {
        internalResize(rows, columns);
    }

    /**
     * Resets current instance by removing its contents, resizing it to provided
     * size and setting all its elements to provided value.
     *
     * @param rows      Number of rows to be set
     * @param columns   Number of columns to be set
     * @param initValue Value to be set in all of its elements
     * @throws WrongSizeException Exception raised if either rows or
     *                            columns is zero.
     */
    public void reset(final int rows, final int columns, final double initValue) throws WrongSizeException {
        internalResize(rows, columns);
        initialize(initValue);
    }

    /**
     * Returns the contents of the matrix as an array of values using
     * DEFAULT_USE_COLUMN_ORDER to pick elements.
     *
     * @return Contents of matrix as an array
     */
    public double[] toArray() {
        return toArray(DEFAULT_USE_COLUMN_ORDER);
    }

    /**
     * Returns the contents of the matrix as an array of values using provided
     * order to pick elements.
     *
     * @param isColumnOrder If true, picks elements from matrix using column
     *                      order, otherwise row order is used.
     * @return Contents of matrix as an array,
     */
    public double[] toArray(final boolean isColumnOrder) {
        final var length = rows * columns;

        if (isColumnOrder) {
            return Arrays.copyOf(buffer, length);
        } else {
            final var out = new double[length];
            double value;
            var counter = 0;

            for (var j = 0; j < rows; j++) {
                for (var i = 0; i < columns; i++) {
                    value = buffer[columnIndex[i] + j];
                    out[counter] = value;
                    counter++;
                }
            }

            return out;
        }
    }

    /**
     * Copies the contents of the matrix to an array of values using column
     * order.
     *
     * @param result array where values will be copied to.
     * @throws WrongSizeException if provided result array does not have the
     *                            same number of elements as the matrix (i.e. rows x columns).
     */
    public void toArray(final double[] result) throws WrongSizeException {
        toArray(result, DEFAULT_USE_COLUMN_ORDER);
    }

    /**
     * Copies the contents of the matrix to an array of values using provided
     * order to pick elements.
     *
     * @param result        array where values will be copied to.
     * @param isColumnOrder if true, picks elements from matrix using column
     *                      order, otherwise row order is used.
     * @throws WrongSizeException if provided result array does not have the
     *                            same number of elements as the matrix (i.e. rows x columns).
     */
    public void toArray(final double[] result, final boolean isColumnOrder) throws WrongSizeException {
        if (result.length != buffer.length) {
            throw new WrongSizeException("result array must be equal to rows x columns");
        }

        if (isColumnOrder) {
            System.arraycopy(buffer, 0, result, 0, buffer.length);
        } else {
            double value;
            var counter = 0;

            for (var j = 0; j < rows; j++) {
                for (var i = 0; i < columns; i++) {
                    value = buffer[columnIndex[i] + j];
                    result[counter] = value;
                    counter++;
                }
            }
        }
    }

    /**
     * Returns current matrix internal buffer of data.
     *
     * @return Internal buffer of data.
     */
    public double[] getBuffer() {
        return buffer;
    }

    /**
     * Obtains a sub-matrix of current matrix instance. Sub-matrix is obtained by
     * copying all elements contained within provided coordinates (both top-left
     * and bottom-right points are included within sub-matrix).
     *
     * @param topLeftRow        Top-left row index where sub-matrix starts.
     * @param topLeftColumn     Top-left column index where sub-matrix starts.
     * @param bottomRightRow    Bottom-right row index where sub-matrix ends.
     * @param bottomRightColumn Bottom-right column index where sub-matrix ends.
     * @param result            Instance where sub-matrix data is stored.
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located below or at right side of bottom-right corner.
     * @throws NullPointerException     If provided result matrix is null.
     */
    public void getSubmatrix(final int topLeftRow, final int topLeftColumn, final int bottomRightRow,
                             final int bottomRightColumn, final Matrix result) {
        if (topLeftRow < 0 || topLeftRow >= rows || topLeftColumn < 0 || topLeftColumn >= columns
                || bottomRightRow < 0 || bottomRightRow >= rows || bottomRightColumn < 0 || bottomRightColumn >= columns
                || topLeftRow > bottomRightRow || topLeftColumn > bottomRightColumn) {
            throw new IllegalArgumentException();
        }

        final var subRows = bottomRightRow - topLeftRow + 1;
        final var subCols = bottomRightColumn - topLeftColumn + 1;
        if (result.getRows() != subRows || result.getColumns() != subCols) {
            // resize result
            try {
                result.resize(subRows, subCols);
            } catch (final WrongSizeException ignore) {
                // never happens
            }
        }
        internalGetSubmatrix(topLeftRow, topLeftColumn, bottomRightRow, bottomRightColumn, result);
    }

    /**
     * Obtains a sub-matrix of current matrix instance. Sub-matrix is obtained by
     * copying all elements contained within provided coordinates (both top-left
     * and bottom-right points are included within sub-matrix).
     *
     * @param topLeftRow        Top-left row index where sub-matrix starts.
     * @param topLeftColumn     Top-left column index where sub-matrix starts.
     * @param bottomRightRow    Bottom-right row index where sub-matrix ends.
     * @param bottomRightColumn Bottom-right column index where sub-matrix ends.
     * @return A new instance containing selected sub-matrix.
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located belo or at right side of bottom-right corner.
     */
    public Matrix getSubmatrix(final int topLeftRow, final int topLeftColumn,
                               final int bottomRightRow, final int bottomRightColumn) {
        if (topLeftRow < 0 || topLeftRow >= rows || topLeftColumn < 0 || topLeftColumn >= columns
                || bottomRightRow < 0 || bottomRightRow >= rows || bottomRightColumn < 0 || bottomRightColumn >= columns
                || topLeftRow > bottomRightRow || topLeftColumn > bottomRightColumn) {
            throw new IllegalArgumentException();
        }

        Matrix out;
        try {
            out = new Matrix(bottomRightRow - topLeftRow + 1, bottomRightColumn - topLeftColumn + 1);
        } catch (final WrongSizeException e) {
            throw new IllegalArgumentException(e);
        }
        internalGetSubmatrix(topLeftRow, topLeftColumn, bottomRightRow, bottomRightColumn, out);
        return out;
    }

    /**
     * Retrieves a sub-matrix of current matrix instance as an array of values
     * using column order and storing the result in provided array.
     * Sub-matrix is obtained by copying all elements contained within provided
     * coordinates (both top-left and bottom-right points are included within
     * sub-matrix).
     *
     * @param topLeftRow        Top-left row index where sub-matrix starts
     * @param topLeftColumn     Top-left column index where sub-matrix starts
     * @param bottomRightRow    Bottom-right row index where sub-matrix ends
     * @param bottomRightColumn Bottom-right column index where sub-matrix ends.
     * @param array             Array where sub-matrix data is stored.
     * @throws IllegalArgumentException Thrown if provided coordinates lie
     *                                  outside of matrix boundaries or if top-left corner is at the bottom or
     *                                  right side of bottom-right corner
     * @throws WrongSizeException       thrown if length of provided array does not
     *                                  match the number of elements to be extracted from this matrix
     */
    public void getSubmatrixAsArray(final int topLeftRow, final int topLeftColumn,
                                    final int bottomRightRow, final int bottomRightColumn,
                                    final double[] array) throws WrongSizeException {
        if (topLeftRow < 0 || topLeftRow >= rows || topLeftColumn < 0 || topLeftColumn >= columns || bottomRightRow < 0
                || bottomRightRow >= rows || bottomRightColumn < 0 || bottomRightColumn >= columns
                || topLeftRow > bottomRightRow || topLeftColumn > bottomRightColumn) {
            throw new IllegalArgumentException();
        }

        final var length = (bottomRightRow - topLeftRow + 1) * (bottomRightColumn - topLeftColumn + 1);
        if (array.length != length) {
            throw new WrongSizeException();
        }

        getSubmatrixAsArray(topLeftRow, topLeftColumn, bottomRightRow, bottomRightColumn, DEFAULT_USE_COLUMN_ORDER,
                array);
    }

    /**
     * Retrieves a sub-matrix of current matrix instance as an array of values
     * using provided column order and storing the result in provided array.
     * Sub-matrix is obtained by copying all elements contained within provided
     * coordinates (both top-left and bottom-right points are included within
     * sub-matrix).
     *
     * @param topLeftRow        Top-left row index where sub-matrix starts
     * @param topLeftColumn     Top-left column index where sub-matrix starts
     * @param bottomRightRow    Bottom-right row index where sub-matrix ends
     * @param bottomRightColumn Bottom-right column index where sub-matrix ends.
     * @param isColumnOrder     If true, picks elements from matrix using column
     *                          order, otherwise row order is used.
     * @param array             Array where sub-matrix data is stored.
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located below or at right side of bottom-right corner.
     * @throws WrongSizeException       If provided array doesn't have proper length,
     *                                  which must be equal to the amount of elements in desired sub-matrix.
     */
    public void getSubmatrixAsArray(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow, final int bottomRightColumn,
            final boolean isColumnOrder, final double[] array) throws WrongSizeException {
        if (topLeftRow < 0 || topLeftRow >= rows || topLeftColumn < 0 || topLeftColumn >= columns || bottomRightRow < 0
                || bottomRightRow >= rows || bottomRightColumn < 0 || bottomRightColumn >= columns
                || topLeftRow > bottomRightRow || topLeftColumn > bottomRightColumn) {
            throw new IllegalArgumentException();
        }

        final var length = (bottomRightRow - topLeftRow + 1) * (bottomRightColumn - topLeftColumn + 1);
        if (array.length != length) {
            throw new WrongSizeException();
        }

        internalGetSubmatrixAsArray(topLeftRow, topLeftColumn, bottomRightRow, bottomRightColumn, isColumnOrder, array);
    }

    /**
     * Obtains sub-matrix of current matrix instance as an array of values using
     * DEFAULT_USE_COLUMN_ORDER. Array is obtained by copying all elements
     * contained within provided coordinates (both top-left
     * and bottom-right points are included within sub-matrix).
     *
     * @param topLeftRow        Top-left row index where sub-matrix starts
     * @param topLeftColumn     Top-left column index where sub-matrix starts
     * @param bottomRightRow    Bottom-right row index where sub-matrix ends
     * @param bottomRightColumn Bottom-right column index where sub-matrix ends.
     * @return An array containing sub-matrix elements
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located belo or at right side of bottom-right corner.
     */
    public double[] getSubmatrixAsArray(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow, final int bottomRightColumn) {
        return getSubmatrixAsArray(topLeftRow, topLeftColumn, bottomRightRow, bottomRightColumn,
                DEFAULT_USE_COLUMN_ORDER);
    }

    /**
     * Obtains sub-matrix of current matrix instance as an array of values using
     * provided order. Array is obtained by copying all elements
     * contained within provided coordinates (both top-left
     * and bottom-right points are included within sub-matrix).
     *
     * @param topLeftRow        Top-left row index where sub-matrix starts
     * @param topLeftColumn     Top-left column index where sub-matrix starts
     * @param bottomRightRow    Bottom-right row index where sub-matrix ends
     * @param bottomRightColumn Bottom-right column index where sub-matrix ends.
     * @param isColumnOrder     If true, picks elements from matrix using column
     *                          order, otherwise row order is used.
     * @return An array containing sub-matrix elements
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located below or at right side of bottom-right corner.
     */
    public double[] getSubmatrixAsArray(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow, final int bottomRightColumn,
            final boolean isColumnOrder) {
        if (topLeftRow < 0 || topLeftRow >= rows || topLeftColumn < 0 || topLeftColumn >= columns || bottomRightRow < 0
                || bottomRightRow >= rows || bottomRightColumn < 0 || bottomRightColumn >= columns
                || topLeftRow > bottomRightRow || topLeftColumn > bottomRightColumn) {
            throw new IllegalArgumentException();
        }

        final var length = (bottomRightRow - topLeftRow + 1) * (bottomRightColumn - topLeftColumn + 1);

        final var out = new double[length];
        internalGetSubmatrixAsArray(topLeftRow, topLeftColumn, bottomRightRow, bottomRightColumn, isColumnOrder, out);
        return out;
    }

    /**
     * Copies elements from provided sub-matrix into this matrix at provided
     * location.
     *
     * @param topLeftRow        Top-left row index where sub-matrix copy starts
     * @param topLeftColumn     Top-left column index where sub-matrix copy starts
     * @param bottomRightRow    Bottom-right row index where sub-matrix copy ends
     * @param bottomRightColumn Bottom-right column index where sub-matrix copy
     *                          ends.
     * @param submatrix         Sub-matrix to be copied
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located below or at right side of bottom-right corner.
     */
    public void setSubmatrix(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow, final int bottomRightColumn,
            final Matrix submatrix) {
        setSubmatrix(topLeftRow, topLeftColumn, bottomRightRow, bottomRightColumn, submatrix, 0,
                0, submatrix.getRows() - 1,
                submatrix.getColumns() - 1);
    }

    /**
     * Copies elements from provided sub-matrix into this matrix at provided
     * location
     *
     * @param topLeftRow                 Top-left row index where sub-matrix copy starts
     * @param topLeftColumn              Top-left column index where sub-matrix copy starts
     * @param bottomRightRow             Bottom-right row index where sub-matrix copy ends
     * @param bottomRightColumn          Bottom-right column index where sub-matrix copy
     *                                   ends.
     * @param submatrix                  Sub-matrix to be copied
     * @param submatrixTopLeftRow        Top-left row index of sub-matrix where copy
     *                                   starts
     * @param submatrixTopLeftColumn     Top-left column index of sub-matrix where
     *                                   copy starts
     * @param submatrixBottomRightRow    Bottom-right row index of sub-matrix where
     *                                   copy ends
     * @param submatrixBottomRightColumn Bottom-right column index of sub-matrix
     *                                   where copy ends
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current or provided matrices, or if
     *                                  top-left corners are indeed located below or at right side of
     *                                  bottom-right corners.
     */
    public void setSubmatrix(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow, final int bottomRightColumn,
            final Matrix submatrix, final int submatrixTopLeftRow, final int submatrixTopLeftColumn,
            final int submatrixBottomRightRow, final int submatrixBottomRightColumn) {
        if (topLeftRow < 0 || topLeftRow >= rows || topLeftColumn < 0 || topLeftColumn >= columns || bottomRightRow < 0
                || bottomRightRow >= rows || bottomRightColumn < 0 || bottomRightColumn >= columns
                || topLeftRow > bottomRightRow || topLeftColumn > bottomRightColumn) {
            throw new IllegalArgumentException();
        }

        if (submatrixTopLeftRow < 0 || submatrixTopLeftColumn < 0 || submatrixBottomRightRow < 0
                || submatrixBottomRightRow >= submatrix.rows || submatrixBottomRightColumn < 0
                || submatrixBottomRightColumn >= submatrix.columns || submatrixTopLeftRow > submatrixBottomRightRow
                || submatrixTopLeftColumn > submatrixBottomRightColumn) {
            throw new IllegalArgumentException();
        }

        final var matrixRows = bottomRightRow - topLeftRow + 1;
        final var matrixColumns = bottomRightColumn - topLeftColumn + 1;
        final var submatrixRows = submatrixBottomRightRow - submatrixTopLeftRow + 1;
        final var submatrixColumns = submatrixBottomRightColumn - submatrixTopLeftColumn + 1;
        if (matrixRows != submatrixRows || matrixColumns != submatrixColumns) {
            throw new IllegalArgumentException();
        }

        var j2 = submatrixTopLeftColumn;
        int destPos;
        int sourcePos;
        for (var j = topLeftColumn; j <= bottomRightColumn; j++) {
            destPos = columnIndex[j] + topLeftRow;
            sourcePos = submatrix.columnIndex[j2] + submatrixTopLeftRow;
            for (var i = topLeftRow; i <= bottomRightRow; i++) {
                // Lines below are equivalent to commented code
                buffer[destPos] = submatrix.buffer[sourcePos];
                destPos++;
                sourcePos++;
            }
            j2++;
        }
    }

    /**
     * Sets elements in provided region to provided value.
     *
     * @param topLeftRow        Top-left row index of region (inclusive).
     * @param topLeftColumn     Top-left column index of region (inclusive).
     * @param bottomRightRow    Bottom-right row index of region (inclusive).
     * @param bottomRightColumn Bottom-right column index of region (inclusive).
     * @param value             Value to be set.
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located below or at right side of bottom-right corner.
     */
    public void setSubmatrix(
            int topLeftRow, int topLeftColumn, int bottomRightRow, int bottomRightColumn, double value) {
        if (topLeftRow < 0 || topLeftRow >= rows || topLeftColumn < 0 || topLeftColumn >= columns || bottomRightRow < 0
                || bottomRightRow >= rows || bottomRightColumn < 0 || bottomRightColumn >= columns
                || topLeftRow > bottomRightRow || topLeftColumn > bottomRightColumn) {
            throw new IllegalArgumentException();
        }

        for (var j = topLeftColumn; j <= bottomRightColumn; j++) {
            for (var i = topLeftRow; i <= bottomRightRow; i++) {
                buffer[columnIndex[j] + i] = value;
            }
        }
    }

    /**
     * Copies elements from provided array into this matrix at provided
     * location. Elements in array are copied into this matrix considering
     * DEFAULT_USE_COLUMN_ORDER.
     *
     * @param topLeftRow        Top-left row index where copy starts.
     * @param topLeftColumn     Top-left column index where copy starts.
     * @param bottomRightRow    Bottom-right row index where copy ends.
     * @param bottomRightColumn Bottom-right column index where copy
     *                          ends.
     * @param values            Array to be copied.
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located below or at right side of bottom-right corner.
     */
    public void setSubmatrix(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow, final int bottomRightColumn,
            final double[] values) {
        setSubmatrix(topLeftRow, topLeftColumn, bottomRightRow, bottomRightColumn, values, DEFAULT_USE_COLUMN_ORDER);
    }

    /**
     * Copies elements from provided array into this matrix at provided
     * location. Elements in array are copied into this matrix following
     * provided order (either row or column order).
     *
     * @param topLeftRow        Top-left row index where copy starts.
     * @param topLeftColumn     Top-left column index where copy ends.
     * @param bottomRightRow    Bottom-right row index where copy ends.
     * @param bottomRightColumn Bottom-right column index where copy ends.
     * @param values            Array to be copied.
     * @param isColumnOrder     If true values are copied consecutively from array
     *                          following column order on the destination matrix, otherwise row order is
     *                          used.
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located below or at right side of bottom-right corner.
     */
    public void setSubmatrix(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow, final int bottomRightColumn,
            final double[] values, final boolean isColumnOrder) {
        setSubmatrix(topLeftRow, topLeftColumn, bottomRightRow, bottomRightColumn, values, 0,
                values.length - 1, isColumnOrder);
    }

    /**
     * Copies elements from provided array into this matrix at provided
     * location. Elements in array are copied into this matrix considering
     * DEFAULT_USE_COLUMN_ORDER starting at provided location until end provided
     * position.
     *
     * @param topLeftRow        Top-left row index where copy starts (inclusive).
     * @param topLeftColumn     Top-left column index where copy starts (inclusive).
     * @param bottomRightRow    Bottom-right row index where copy ends (inclusive).
     * @param bottomRightColumn Bottom-right column index where copy
     *                          ends. (inclusive)
     * @param values            Array to be copied.
     * @param valuesStart       Position where copy from array will start (inclusive).
     * @param valuesEnd         Position where copy from array will finish (inclusive).
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located below or at right side of bottom-right corner,
     *                                  or if valuesStart and valuesEnd lie outside valid array positions or
     *                                  start is greater than end position.
     */
    public void setSubmatrix(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow, final int bottomRightColumn,
            final double[] values, final int valuesStart, final int valuesEnd) {
        setSubmatrix(topLeftRow, topLeftColumn, bottomRightRow, bottomRightColumn, values, valuesStart, valuesEnd,
                DEFAULT_USE_COLUMN_ORDER);
    }

    /**
     * Copies elements from provided array into this matrix at provided
     * location. Elements in array are copied into this matrix following
     * provided order (either row or column order) starting at provided location
     * until end provided position.
     *
     * @param topLeftRow        Top-left row index where copy starts (inclusive)
     * @param topLeftColumn     Top-left column index where copy starts (inclusive)
     * @param bottomRightRow    Bottom-right row index where copy ends (inclusive)
     * @param bottomRightColumn Bottom-right column index where copy
     *                          ends. (inclusive)
     * @param values            Array to be copied.
     * @param valuesStart       Position where copy from array will start (inclusive).
     * @param valuesEnd         Position where copy from array will finish (inclusive).
     * @param isColumnOrder     If true values are copied consecutively from array
     *                          following column order on the destination matrix, otherwise row order is
     *                          used.
     * @throws IllegalArgumentException Exception raised whenever top-left or
     *                                  bottom-right corners lie outside current matrix instance, or if top-left
     *                                  corner is indeed located below or at right side of bottom-right corner,
     *                                  or if valuesStart and valuesEnd lie outside valid array positions or
     *                                  start is greater than end position.
     */
    public void setSubmatrix(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow, final int bottomRightColumn,
            final double[] values, final int valuesStart, final int valuesEnd, final boolean isColumnOrder) {

        if (topLeftRow < 0 || topLeftRow >= rows || topLeftColumn < 0 || topLeftColumn >= columns || bottomRightRow < 0
                || bottomRightRow >= rows || bottomRightColumn < 0 || bottomRightColumn >= columns
                || topLeftRow > bottomRightRow || topLeftColumn > bottomRightColumn) {
            throw new IllegalArgumentException();
        }

        if (valuesStart < 0 || valuesEnd < 0 || valuesEnd >= values.length || valuesStart > valuesEnd) {
            throw new IllegalArgumentException();
        }

        final var matrixRows = bottomRightRow - topLeftRow + 1;
        final var matrixColumns = bottomRightColumn - topLeftColumn + 1;
        final var matrixLength = matrixRows * matrixColumns;
        final var valuesLength = valuesEnd - valuesStart + 1;

        if (matrixLength != valuesLength) {
            throw new IllegalArgumentException();
        }

        var counter = valuesStart;
        if (isColumnOrder) {
            int destPos;
            for (var j = topLeftColumn; j <= bottomRightColumn; j++) {
                destPos = columnIndex[j] + topLeftRow;
                for (var i = topLeftRow; i <= bottomRightRow; i++) {
                    // Two Lines below are equivalent to:
                    // buffer[columnIndex[j] + i] = values[counter]
                    buffer[destPos] = values[counter];
                    destPos++;
                    counter++;
                }
            }

        } else {
            for (var i = topLeftRow; i <= bottomRightRow; i++) {
                for (var j = topLeftColumn; j <= bottomRightColumn; j++) {
                    buffer[columnIndex[j] + i] = values[counter];
                    counter++;
                }
            }
        }
    }

    /**
     * Sets values into provided matrix to make it an identity matrix (all
     * elements in the diagonal equal to one, and remaining elements to zero).
     *
     * @param m Matrix where identity values are set.
     */
    public static void identity(final Matrix m) {
        m.initialize(0.0);

        final var minSize = Math.min(m.rows, m.columns);
        for (var i = 0; i < minSize; i++) {
            m.buffer[m.columnIndex[i] + i] = 1.0;
        }
    }

    /**
     * Creates and returns a new matrix instance having all the elements on
     * the diagonal equal to one and the remaining ones equal to zero.
     *
     * @param rows    Number of rows of created instance
     * @param columns Number of columns of created instance
     * @return An identity matrix
     * @throws WrongSizeException Raised if either rows or columns is
     *                            equal to zero
     */
    public static Matrix identity(final int rows, final int columns) throws WrongSizeException {
        final var out = new Matrix(rows, columns);
        identity(out);
        return out;
    }

    /**
     * Fills provided matrix with random uniform values ranging from minValue to
     * maxValue.
     *
     * @param minValue Minimum value of uniform random generated values
     * @param maxValue Maximum value of uniform random generated values
     * @param random   Random generator
     * @param result   Matrix where random values are stored.
     * @throws IllegalArgumentException if minValue &lt;= maxValue
     * @throws NullPointerException     if either provided random or result
     *                                  instances are null
     */
    public static void fillWithUniformRandomValues(
            final double minValue, final double maxValue, final Random random, final Matrix result) {
        final var randomizer = new UniformRandomizer(random);

        final var length = result.rows * result.columns;

        for (int i = 0; i < length; i++) {
            result.buffer[i] = randomizer.nextDouble(minValue, maxValue);
        }
    }

    /**
     * Fills provided matrix with random uniform values ranging from minValue to
     * maxValue.
     *
     * @param minValue Minimum value of uniform random generated values
     * @param maxValue Maximum value of uniform random generated values
     * @param result   Matrix where random values are stored
     * @throws IllegalArgumentException if minValue &lt;= maxValue
     * @throws NullPointerException     if provided result matrix is null
     */
    public static void fillWithUniformRandomValues(final double minValue, final double maxValue, final Matrix result) {
        fillWithUniformRandomValues(minValue, maxValue, new Random(), result);
    }

    /**
     * Creates new matrix instance using provided size and containing uniformly
     * distributed random values with provided range.
     *
     * @param rows     Number of rows of instantiated matrix
     * @param columns  Number of columns of instantiated matrix
     * @param minValue Minimum value of uniform random generated values
     * @param maxValue Maximum value of uniform random generated values
     * @return A new matrix containing uniform random values
     * @throws WrongSizeException       Exception thrown if either rows or
     *                                  columns is zero, or if the minimum random value is greater or equal
     *                                  than the maximum random value
     * @throws IllegalArgumentException if minValue &lt;= maxValue
     */
    public static Matrix createWithUniformRandomValues(
            final int rows, final int columns, final double minValue, final double maxValue) throws WrongSizeException {
        return createWithUniformRandomValues(rows, columns, minValue, maxValue, new Random());
    }

    /**
     * Creates new matrix instance using provided size and containing uniformly
     * distributed random values with provided range and using provided random
     * generator
     *
     * @param rows     Number of rows of instantiated matrix
     * @param columns  Number of columns of instantiated matrix
     * @param minValue Minimum value of uniform random generated values
     * @param maxValue Maximum value of uniform random generated values
     * @param random   A random generator.
     * @return A new matrix containing uniform random values
     * @throws WrongSizeException       Exception thrown if either rows or
     *                                  columns is zero, or if the minimum random value is greater or equal
     *                                  than the maximum random value
     * @throws IllegalArgumentException if minValue &lt;= maxValue
     */
    public static Matrix createWithUniformRandomValues(
            final int rows, final int columns, final double minValue, final double maxValue, final Random random)
            throws WrongSizeException {

        final var out = new Matrix(rows, columns);
        fillWithUniformRandomValues(minValue, maxValue, random, out);
        return out;
    }

    /**
     * Fills provided matrix with random Gaussian values with provided mean and
     * standard deviation.
     *
     * @param mean              Mean value of generated random values
     * @param standardDeviation Standard deviation of generated random values
     * @param random            Random generator
     * @param result            Matrix where random values are stored
     * @throws IllegalArgumentException if standard deviation is negative or
     *                                  zero
     * @throws NullPointerException     if provided result matrix is null
     */
    public static void fillWithGaussianRandomValues(
            final double mean, final double standardDeviation, final Random random, final Matrix result) {

        final var randomizer = new GaussianRandomizer(random, mean, standardDeviation);

        final var length = result.rows * result.columns;

        for (var i = 0; i < length; i++) {
            result.buffer[i] = randomizer.nextDouble();
        }
    }

    /**
     * Fills provided matrix with random Gaussian values with provided mean and
     * standard deviation
     *
     * @param mean              Mean value of generated random values
     * @param standardDeviation Standard deviation of generated random values
     * @param result            Matrix where random values are stored
     * @throws IllegalArgumentException if standard deviation is negative or
     *                                  zero
     * @throws NullPointerException     if provided result matrix is null
     */
    public static void fillWithGaussianRandomValues(
            final double mean, final double standardDeviation, final Matrix result) {
        fillWithGaussianRandomValues(mean, standardDeviation, new Random(), result);
    }

    /**
     * Creates new matrix instance using provided size and containing
     * gaussian/normal distributed random values with provided median and
     * standard deviation.
     *
     * @param rows              Number of rows of instantiated matrix
     * @param columns           Number of columns of instantiated matrix
     * @param mean              Mean value of gaussian random generated values
     * @param standardDeviation Standard deviation of gaussian random generated
     *                          values
     * @return A new matrix containing gaussian random values
     * @throws WrongSizeException       Exception thrown if either rows or
     *                                  columns is zero, or if the standard deviation is negative or zero.
     * @throws IllegalArgumentException thrown if provided standard deviation
     *                                  is negative or zero.
     */
    public static Matrix createWithGaussianRandomValues(
            final int rows, final int columns, final double mean, final double standardDeviation)
            throws WrongSizeException {
        return createWithGaussianRandomValues(rows, columns, mean, standardDeviation, new Random());
    }

    /**
     * Creates new matrix instance using provided size and containing
     * gaussian/normal distributed random values with provided median and
     * standard deviation and using provided random generator
     *
     * @param rows              Number of rows of instantiated matrix
     * @param columns           Number of columns of instantiated matrix
     * @param mean              Mean value of gaussian random generated values
     * @param standardDeviation Standard deviation of gaussian random generated
     *                          values
     * @param random            A random generator.
     * @return A new matrix containing gaussian random values
     * @throws WrongSizeException       Exception thrown if either rows or
     *                                  columns is zero, or if the standard deviation is negative or zero.
     * @throws IllegalArgumentException thrown if provided standard deviation
     *                                  is negative or zero.
     */
    public static Matrix createWithGaussianRandomValues(
            final int rows, final int columns, final double mean, final double standardDeviation, final Random random)
            throws WrongSizeException {

        final var out = new Matrix(rows, columns);
        fillWithGaussianRandomValues(mean, standardDeviation, random, out);
        return out;
    }

    /**
     * Makes provided result matrix a diagonal matrix containing provided
     * elements in the diagonal. Elements outside the diagonal will be set to
     * zero
     *
     * @param diagonal Array containing the elements to be set on the diagonal
     * @param result   Matrix where values are stored
     * @throws NullPointerException Exception thrown if provided parameters are
     *                              null
     */
    public static void diagonal(final double[] diagonal, final Matrix result) {

        result.initialize(0.0);
        // set diagonal elements
        //noinspection ManualArrayCopy
        for (var i = 0; i < diagonal.length; i++) {
            result.buffer[result.columnIndex[i] + i] = diagonal[i];
        }
    }

    /**
     * Creates a diagonal matrix having all the elements in provided array in
     * its diagonal and the remaining elements equal to zero.
     * Returned matrix will have size n x n, where n is the length of the array
     *
     * @param diagonal Array containing the elements to be set on the diagonal
     * @return A diagonal matrix
     * @throws NullPointerException Raised if provided diagonal array is null
     */
    public static Matrix diagonal(final double[] diagonal) {

        Matrix out = null;
        try {
            out = new Matrix(diagonal.length, diagonal.length);
            diagonal(diagonal, out);
        } catch (final WrongSizeException ignore) {
            // never happens
        }

        return out;
    }

    /**
     * Instantiates new matrix from array using DEFAULT_USE_COLUMN_ORDER
     *
     * @param array Array used as source to copy values from
     * @return Returns matrix created from array
     */
    public static Matrix newFromArray(final double[] array) {
        return newFromArray(array, DEFAULT_USE_COLUMN_ORDER);
    }

    /**
     * Instantiates new matrix from array using either column or row order
     *
     * @param array         Array used as source to copy values from
     * @param isColumnOrder True if column order must be used, false otherwise
     * @return Returns matrix created from array
     */
    public static Matrix newFromArray(final double[] array, final boolean isColumnOrder) {
        Matrix m = null;
        try {
            if (isColumnOrder) {
                m = new Matrix(array.length, 1);
                m.setSubmatrix(0, 0, array.length - 1, 0, array);
            } else {
                m = new Matrix(1, array.length);
                m.setSubmatrix(0, 0, 0, array.length - 1, array);
            }
        } catch (final WrongSizeException ignore) {
            // never happens
        }
        return m;
    }

    /**
     * Copies elements of array into this instance using column order.
     *
     * @param array array to copy values from.
     * @throws WrongSizeException if provided array length is not equal to the
     *                            number of rows multiplied per the number of columns of this instance.
     */
    public void fromArray(final double[] array) throws WrongSizeException {
        fromArray(array, DEFAULT_USE_COLUMN_ORDER);
    }

    /**
     * Copies elements of array into this instance using provided order.
     *
     * @param array         array to copy values from.
     * @param isColumnOrder true to use column order, false otherwise.
     * @throws WrongSizeException if provided array length is not equal to the
     *                            number of rows multiplied per the number of columns of this instance.
     */
    public void fromArray(final double[] array, final boolean isColumnOrder) throws WrongSizeException {
        if (array.length != buffer.length) {
            throw new WrongSizeException("array length must be equal to rows x columns");
        }

        if (isColumnOrder) {
            System.arraycopy(array, 0, buffer, 0, array.length);
        } else {
            var counter = 0;

            for (var j = 0; j < rows; j++) {
                for (var i = 0; i < columns; i++) {
                    buffer[columnIndex[i] + j] = array[counter];
                    counter++;
                }
            }
        }
    }

    /**
     * Symmetrizes this instance and stores the result into provided instance.
     * Symmetrization is done by averaging this instance with its transpose
     * (i.e. S = (M+M')/2
     *
     * @param result instance where symmetrized version of this instance will be
     *               stored.
     * @throws WrongSizeException if this instance is not square or provided
     *                            result instance doesn't have the same size as this instance.
     */
    public void symmetrize(final Matrix result) throws WrongSizeException {
        if (rows != columns) {
            throw new WrongSizeException("matrix must be square");
        }
        if (result.getRows() != rows || result.getColumns() != columns) {
            throw new WrongSizeException("result matrix must have the size of this instance");
        }

        // S = (M+M')/2
        double value1;
        double value2;
        double avg;
        int pos1;
        int pos2;
        for (var i = 0; i < columns; i++) {
            for (var j = i; j < rows; j++) {
                // value at (i, j)
                pos1 = columnIndex[i] + j;
                value1 = buffer[pos1];
                // transposed value (i.e. value at (j,i))
                pos2 = columnIndex[j] + i;
                value2 = buffer[pos2];

                avg = 0.5 * (value1 + value2);
                result.buffer[pos1] = avg;
                result.buffer[pos2] = avg;
            }
        }
    }

    /**
     * Symmetrizes this instance and returns the result as a new matrix
     * instance.
     * Symmetrization is done by averaging this instance with its transpose
     * (i.e. S = (M+M')/2
     *
     * @return a new symmetrizes version of this instance.
     * @throws WrongSizeException if this instance is not square.
     */
    public Matrix symmetrizeAndReturnNew() throws WrongSizeException {
        final var m = new Matrix(rows, columns);
        symmetrize(m);
        return m;
    }

    /**
     * Symmetrizes this instance and updates it with computed value.
     * Symmetrization is done by averaging this instance with its transpose
     * (i.e. S = (M+M')/2
     *
     * @throws WrongSizeException if this instance is not square.
     */
    public void symmetrize() throws WrongSizeException {
        symmetrize(this);
    }

    /**
     * Method to internally add two matrices.
     *
     * @param other  Matrix to be added to current matrix
     * @param result Matrix where result will be stored.
     */
    private void internalAdd(final Matrix other, final Matrix result) {
        final var length = rows * columns;
        for (var i = 0; i < length; i++) {
            result.buffer[i] = buffer[i] + other.buffer[i];
        }
    }

    /**
     * Method to internally subtract two matrices.
     *
     * @param other  Matrix to be subtracted from current matrix.
     * @param result Matrix where result will be stored.
     */
    private void internalSubtract(final Matrix other, final Matrix result) {
        final var length = rows * columns;
        for (var i = 0; i < length; i++) {
            result.buffer[i] = buffer[i] - other.buffer[i];
        }
    }

    /**
     * Method to internally multiply two matrices.
     *
     * @param other             Matrix to be multiplied to current matrix
     * @param resultBuffer      Matrix buffer of data where result will be stored.
     * @param resultColumnIndex Array of matrix column indices where result will
     *                          be stored.
     */
    private void internalMultiply(
            final Matrix other, final double[] resultBuffer, final int[] resultColumnIndex) {
        final var columns2 = other.columns;
        double value;
        for (var k = 0; k < columns2; k++) {
            for (var j = 0; j < rows; j++) {
                value = 0.0;
                for (var i = 0; i < columns; i++) {
                    value += buffer[columnIndex[i] + j] *
                            other.buffer[other.columnIndex[k] + i];
                }
                resultBuffer[resultColumnIndex[k] + j] = value;
            }
        }
    }

    /**
     * Method to internally multiply two matrices.
     *
     * @param other  Matrix to be multiplied to current matrix.
     * @param result Matrix where result will be stored.
     */
    private void internalMultiply(final Matrix other, final Matrix result) {
        internalMultiply(other, result.buffer, result.columnIndex);
    }

    /**
     * Method to internally compute the Kronecker product between two matrices.
     *
     * @param other             other matrix to be Kronecker multiplied to current matrix.
     * @param resultBuffer      matrix buffer of data where result will be stored.
     * @param resultColumnIndex array of matrix column indices where result will
     *                          be stored.
     */
    private void internalMultiplyKronecker(
            final Matrix other, final double[] resultBuffer, final int[] resultColumnIndex) {
        final var rows2 = other.rows;
        final var columns2 = other.columns;

        for (var j1 = 0; j1 < rows; j1++) {
            final var startJ3 = j1 * other.rows;
            for (var i1 = 0; i1 < columns; i1++) {
                final var startI3 = i1 * other.columns;
                final var value1 = buffer[columnIndex[i1] + j1];

                for (var j2 = 0; j2 < rows2; j2++) {
                    final var j3 = startJ3 + j2;
                    for (var i2 = 0; i2 < columns2; i2++) {
                        final var i3 = startI3 + i2;
                        final var value2 = other.buffer[other.columnIndex[i2] + j2];

                        final var value3 = value1 * value2;
                        resultBuffer[resultColumnIndex[i3] + j3] = value3;
                    }
                }
            }
        }
    }

    /**
     * Method to internally compute the Kronecker product between two matrices.
     *
     * @param other  other matrix to be Kronecker multiplied to current matrix.
     * @param result matrix where result will be stored.
     */
    private void internalMultiplyKronecker(final Matrix other, final Matrix result) {
        internalMultiplyKronecker(other, result.buffer, result.columnIndex);
    }

    /**
     * Method to internally compute element by element product of two matrices.
     *
     * @param other  Matrix to be element by element multiplied to current matrix
     * @param result Matrix where result will be stored.
     */
    private void internalElementByElementProduct(final Matrix other, final Matrix result) {
        final var length = rows * columns;
        for (var i = 0; i < length; i++) {
            result.buffer[i] = buffer[i] * other.buffer[i];
        }
    }

    /**
     * Method to internally compute matrix transposition.
     *
     * @param resultBuffer      Buffer where transposed matrix data is stored.
     * @param resultColumnIndex Buffer where indices of transposed matrix data
     *                          is stored.
     */
    private void internalTranspose(final double[] resultBuffer, final int[] resultColumnIndex) {
        for (var j = 0; j < rows; j++) {
            for (var i = 0; i < columns; i++) {
                resultBuffer[resultColumnIndex[j] + i] = buffer[columnIndex[i] + j];
            }
        }
    }

    /**
     * Method to internally compute matrix transposition.
     *
     * @param result Matrix where transposed data is stored.
     */
    private void internalTranspose(final Matrix result) {
        internalTranspose(result.buffer, result.columnIndex);
    }

    /**
     * Method used internally to remove matrix contents and resizing it.
     *
     * @param rows    Number of rows to be set
     * @param columns Number of columns to be set.
     * @throws WrongSizeException Exception raised if either rows or
     *                            columns is zero.
     */
    private void internalResize(final int rows, final int columns) throws WrongSizeException {
        if (rows == 0 || columns == 0) {
            throw new WrongSizeException();
        }

        this.rows = rows;
        this.columns = columns;

        // instantiate buffers of data
        buffer = new double[rows * columns];
        columnIndex = new int[columns];

        // initialize column index
        var counter = 0;
        for (var i = 0; i < columns; i++) {
            columnIndex[i] = counter;
            counter += rows;
        }
    }

    /**
     * Internal method to retrieve a sub-matrix of current matrix instance.
     * Sub-matrix is obtained by copying all elements contained within provided
     * coordinates (both top-left and bottom-right points are included within
     * sub-matrix).
     *
     * @param topLeftRow        Top-left row index where sub-matrix starts.
     * @param topLeftColumn     Top-left column index where sub-matrix starts.
     * @param bottomRightRow    Bottom-right row index where sub-matrix ends.
     * @param bottomRightColumn Bottom-right column index where sub-matrix ends.
     * @param result            Instance where sub-matrix data is stored.
     */
    private void internalGetSubmatrix(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow, final int bottomRightColumn,
            final Matrix result) {
        var i2 = 0;
        var j2 = 0;
        for (var j = topLeftColumn; j <= bottomRightColumn; j++) {
            for (var i = topLeftRow; i <= bottomRightRow; i++) {
                result.buffer[result.columnIndex[j2] + i2] = buffer[columnIndex[j] + i];
                i2++;
            }
            i2 = 0;
            j2++;
        }
    }

    /**
     * Internal method to retrieve a sub-matrix of current matrix instance as an
     * array of values using provided column order and storing the result in
     * provided array.
     * Sub-matrix is obtained by copying all elements contained within provided
     * coordinates (both top-left and bottom-right points are included within
     * sub-matrix).
     *
     * @param topLeftRow        Top-left row index where sub-matrix starts
     * @param topLeftColumn     Top-left column index where sub-matrix starts
     * @param bottomRightRow    Bottom-right row index where sub-matrix ends
     * @param bottomRightColumn Bottom-right column index where sub-matrix ends.
     * @param isColumnOrder     If true, picks elements from matrix using column
     *                          order, otherwise row order is used.
     * @param result            Array where sub-matrix data is stored.
     */
    private void internalGetSubmatrixAsArray(
            final int topLeftRow, final int topLeftColumn, final int bottomRightRow,
            final int bottomRightColumn, final boolean isColumnOrder, final double[] result) {
        var counter = 0;
        if (isColumnOrder) {
            for (var j = topLeftColumn; j <= bottomRightColumn; j++) {
                for (var i = topLeftRow; i <= bottomRightRow; i++) {
                    result[counter] = buffer[columnIndex[j] + i];
                    counter++;
                }
            }
        } else {
            for (var i = topLeftRow; i <= bottomRightRow; i++) {
                for (var j = topLeftColumn; j <= bottomRightColumn; j++) {
                    result[counter] = buffer[columnIndex[j] + i];
                    counter++;
                }
            }
        }
    }
}