DoubleMatrixDataSparse.java

package org.djunits.value.vdouble.matrix.data;

import java.util.Arrays;
import java.util.Collection;
import java.util.concurrent.atomic.AtomicInteger;
import java.util.stream.IntStream;

import org.djunits.unit.Unit;
import org.djunits.unit.scale.Scale;
import org.djunits.value.ValueRuntimeException;
import org.djunits.value.storage.StorageType;
import org.djunits.value.vdouble.function.DoubleFunction;
import org.djunits.value.vdouble.function.DoubleFunction2;
import org.djunits.value.vdouble.matrix.base.DoubleSparseValue;
import org.djunits.value.vdouble.scalar.base.DoubleScalar;
import org.djutils.exceptions.Throw;

/**
 * Stores sparse data for a DoubleMatrix and carries out basic operations. The index in the sparse matrix data is calculated as
 * <code>r * columns + c</code>, where r is the row number, cols is the total number of columns, and c is the column number.
 * <p>
 * Copyright (c) 2013-2024 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
 * BSD-style license. See <a href="https://djunits.org/docs/license.html">DJUNITS License</a>.
 * </p>
 * @author <a href="https://www.tudelft.nl/averbraeck">Alexander Verbraeck</a>
 * @author <a href="https://www.tudelft.nl/staff/p.knoppers/">Peter Knoppers</a>
 */
public class DoubleMatrixDataSparse extends DoubleMatrixData
{
    /** */
    private static final long serialVersionUID = 1L;

    /** the index values of the Matrix. The index is calculated as r * columns + c. */
    private long[] indices;

    /**
     * Create a matrix with sparse data.
     * @param matrixSI double[]; the data to store
     * @param indices long[]; the index values of the Matrix, with &lt;tt&gt;index = row * cols + col&lt;/tt&gt;
     * @param rows int; the number of rows
     * @param cols int; the number of columns
     */
    public DoubleMatrixDataSparse(final double[] matrixSI, final long[] indices, final int rows, final int cols)
    {
        super(StorageType.SPARSE);
        this.matrixSI = matrixSI;
        this.indices = indices;
        this.rows = rows;
        this.cols = cols;
    }

    /**
     * Create a matrix with sparse data.
     * @param dataSI Collection&lt;DoubleSparseValue&lt;U, S&gt;&gt;; the sparse [X, Y, SI] values to store
     * @param rows int; the number of rows of the matrix
     * @param cols int; the number of columns of the matrix
     * @throws NullPointerException when storageType is null or dataSI is null
     * @param <U> the unit
     * @param <S> the corresponding scalar type
     */
    public <U extends Unit<U>, S extends DoubleScalar<U, S>> DoubleMatrixDataSparse(
            final Collection<DoubleSparseValue<U, S>> dataSI, final int rows, final int cols) throws NullPointerException
    {
        super(StorageType.SPARSE);
        Throw.whenNull(dataSI, "matrixSI is null");

        int length = (int) dataSI.stream().parallel().filter(d -> d.getValueSI() != 0.0).count();
        this.rows = rows;
        this.cols = cols;
        this.matrixSI = new double[length];
        this.indices = new long[length];
        int index = 0;
        for (DoubleSparseValue<U, S> data : dataSI)
        {
            if (data.getValueSI() != 0.0)
            {
                this.indices[index] = data.getRow() * this.cols + data.getColumn();
                this.matrixSI[index] = data.getValueSI();
                index++;
            }
        }
    }

    @Override
    public final int cardinality()
    {
        return this.indices.length;
    }

    /**
     * Fill the sparse data structures matrixSI[] and indices[]. Note: output vectors have to be initialized at the right size.
     * Cannot be parallelized because of stateful and sequence-sensitive count.
     * @param data double[][]; the input data
     * @param matrixSI double[]; the matrix data to write
     * @param indices long[]; the indices to write
     * @throws ValueRuntimeException in case matrix is ragged
     */
    @SuppressWarnings("checkstyle:finalparameters")
    private static void fill(final double[][] data, double[] matrixSI, long[] indices) throws ValueRuntimeException
    {
        int rows = data.length;
        int cols = rows == 0 ? 0 : data[0].length;
        if (cols == 0)
        {
            rows = 0;
        }
        int count = 0;
        for (int r = 0; r < rows; r++)
        {
            double[] row = data[r];
            // Row length check has been done by checkRectangularAndNonEmpty
            for (int c = 0; c < cols; c++)
            {
                int index = r * cols + c;
                if (row[c] != 0.0)
                {
                    matrixSI[count] = row[c];
                    indices[count] = index;
                    count++;
                }
            }
        }
    }

    /**
     * Fill the sparse data structures matrixSI[] and indices[]. Note: output vectors have to be initialized at the right size.
     * Cannot be parallelized because of stateful and sequence-sensitive count.
     * @param data double[][]; the input data
     * @param matrixSI double[]; the matrix data to write
     * @param indices long[]; the indices to write
     * @param scale Scale; Scale, the scale that will convert the data to SI
     * @throws ValueRuntimeException in case matrix is ragged
     */
    @SuppressWarnings("checkstyle:finalparameters")
    private static void fill(final double[][] data, double[] matrixSI, long[] indices, final Scale scale)
            throws ValueRuntimeException
    {
        int rows = data.length;
        int cols = rows == 0 ? 0 : data[0].length;
        if (cols == 0)
        {
            rows = 0;
        }
        int count = 0;
        for (int r = 0; r < rows; r++)
        {
            double[] row = data[r];
            // Row length check has been done by checkRectangularAndNonEmpty
            for (int c = 0; c < cols; c++)
            {
                int index = r * cols + c;
                double value = scale.toStandardUnit(row[c]);
                if (value != 0.0)
                {
                    matrixSI[count] = value;
                    indices[count] = index;
                    count++;
                }
            }
        }
    }

    @Override
    public DoubleMatrixData assign(final DoubleFunction doubleFunction)
    {
        if (doubleFunction.apply(0d) != 0d)
        {
            // It is most unlikely that the result AND the left and right operands are efficiently stored in Sparse format
            DoubleMatrixDataSparse result = toDense().assign(doubleFunction).toSparse();
            this.indices = result.indices;
            this.matrixSI = result.matrixSI;
            return this;
        }
        // The code below relies on the fact that doubleFunction.apply(0d) yields 0d
        int currentSize = rows() * cols();
        if (currentSize > 16)
        {
            currentSize = 16;
        }
        long[] newIndices = new long[currentSize];
        double[] newValues = new double[currentSize];
        int nonZeroValues = 0;
        int ownIndex = 0;
        while (ownIndex < this.indices.length)
        {
            long index = this.indices[ownIndex];
            double value = doubleFunction.apply(this.matrixSI[ownIndex]);
            ownIndex++;
            if (value != 0d)
            {
                if (nonZeroValues >= currentSize)
                {
                    // increase size of arrays
                    currentSize *= 2;
                    if (currentSize > rows() * cols())
                    {
                        currentSize = rows() * cols();
                    }
                    long[] newNewIndices = new long[currentSize];
                    System.arraycopy(newIndices, 0, newNewIndices, 0, newIndices.length);
                    newIndices = newNewIndices;
                    double[] newNewValues = new double[currentSize];
                    System.arraycopy(newValues, 0, newNewValues, 0, newValues.length);
                    newValues = newNewValues;
                }
                newIndices[nonZeroValues] = index;
                newValues[nonZeroValues] = value;
                nonZeroValues++;
            }
        }
        if (nonZeroValues < currentSize)
        {
            // reduce size of arrays
            long[] newNewIndices = new long[nonZeroValues];
            System.arraycopy(newIndices, 0, newNewIndices, 0, nonZeroValues);
            newIndices = newNewIndices;
            double[] newNewValues = new double[nonZeroValues];
            System.arraycopy(newValues, 0, newNewValues, 0, nonZeroValues);
            newValues = newNewValues;
        }
        this.indices = newIndices;
        this.matrixSI = newValues;
        return this;
    }

    @Override
    public final DoubleMatrixDataSparse assign(final DoubleFunction2 doubleFunction, final DoubleMatrixData right)
    {
        if (doubleFunction.apply(0d, 0d) != 0d)
        {
            // It is most unlikely that the result AND the left and right operands are efficiently stored in Sparse format
            DoubleMatrixDataSparse result = toDense().assign(doubleFunction, right).toSparse();
            this.indices = result.indices;
            this.matrixSI = result.matrixSI;
            return this;
        }
        // The code below relies on the fact that doubleFunction.apply(0d, 0d) yields 0d
        checkSizes(right);
        int currentSize = rows() * cols();
        if (currentSize > 16)
        {
            currentSize = 16;
        }
        long[] newIndices = new long[currentSize];
        double[] newValues = new double[currentSize];
        int nonZeroValues = 0;
        int ownIndex = 0;
        int otherIndex = 0;
        if (right.isSparse())
        { // both are sparse, result must be sparse
            DoubleMatrixDataSparse other = (DoubleMatrixDataSparse) right;
            while (ownIndex < this.indices.length || otherIndex < other.indices.length)
            {
                double value;
                long index;
                if (ownIndex < this.indices.length && otherIndex < other.indices.length)
                { // neither we nor right has run out of values
                    if (this.indices[ownIndex] == other.indices[otherIndex])
                    {
                        value = doubleFunction.apply(this.matrixSI[ownIndex], other.matrixSI[otherIndex]);
                        index = this.indices[ownIndex];
                        ownIndex++;
                        otherIndex++;
                    }
                    else if (this.indices[ownIndex] < other.indices[otherIndex])
                    {
                        // we have a non-zero; right has a zero
                        value = doubleFunction.apply(this.matrixSI[ownIndex], 0d);
                        index = this.indices[ownIndex];
                        ownIndex++;
                    }
                    else
                    {
                        // we have a zero; right has a non-zero
                        value = doubleFunction.apply(0d, other.matrixSI[otherIndex]);
                        index = other.indices[otherIndex];
                        otherIndex++;
                    }
                }
                else if (ownIndex < this.indices.length)
                { // right has run out of values; we have not
                    value = doubleFunction.apply(this.matrixSI[ownIndex], 0d);
                    index = this.indices[ownIndex];
                    ownIndex++;
                }
                else
                { // we have run out of values; right has not
                    value = doubleFunction.apply(0d, other.matrixSI[otherIndex]);
                    index = other.indices[otherIndex];
                    otherIndex++;
                }
                if (value != 0d)
                {
                    if (nonZeroValues >= currentSize)
                    {
                        // increase size of arrays
                        currentSize *= 2;
                        if (currentSize > rows() * cols())
                        {
                            currentSize = rows() * cols();
                        }
                        long[] newNewIndices = new long[currentSize];
                        System.arraycopy(newIndices, 0, newNewIndices, 0, newIndices.length);
                        newIndices = newNewIndices;
                        double[] newNewValues = new double[currentSize];
                        System.arraycopy(newValues, 0, newNewValues, 0, newValues.length);
                        newValues = newNewValues;
                    }
                    newIndices[nonZeroValues] = index;
                    newValues[nonZeroValues] = value;
                    nonZeroValues++;
                }
            }
        }
        else
        { // we are sparse; other is dense, result must be sparse
            DoubleMatrixDataDense other = (DoubleMatrixDataDense) right;
            while (otherIndex < right.matrixSI.length)
            {
                double value;
                int index = otherIndex;
                if (ownIndex < this.indices.length)
                { // neither we nor right has run out of values
                    if (this.indices[ownIndex] == otherIndex)
                    {
                        value = doubleFunction.apply(this.matrixSI[ownIndex], other.matrixSI[otherIndex]);
                        ownIndex++;
                    }
                    else
                    {
                        // we have a zero; other has a value
                        value = doubleFunction.apply(0d, other.matrixSI[otherIndex]);
                    }
                    otherIndex++;
                }
                else
                { // we have run out of values; right has not
                    value = doubleFunction.apply(0d, other.matrixSI[otherIndex]);
                    otherIndex++;
                }
                if (value != 0d)
                {
                    if (nonZeroValues >= currentSize)
                    {
                        // increase size of arrays
                        currentSize *= 2;
                        if (currentSize > rows() * cols())
                        {
                            currentSize = rows() * cols();
                        }
                        long[] newNewIndices = new long[currentSize];
                        System.arraycopy(newIndices, 0, newNewIndices, 0, newIndices.length);
                        newIndices = newNewIndices;
                        double[] newNewValues = new double[currentSize];
                        System.arraycopy(newValues, 0, newNewValues, 0, newValues.length);
                        newValues = newNewValues;
                    }
                    newIndices[nonZeroValues] = index;
                    newValues[nonZeroValues] = value;
                    nonZeroValues++;
                }
            }
        }
        if (nonZeroValues < currentSize)
        {
            // reduce size of arrays
            long[] newNewIndices = new long[nonZeroValues];
            System.arraycopy(newIndices, 0, newNewIndices, 0, nonZeroValues);
            newIndices = newNewIndices;
            double[] newNewValues = new double[nonZeroValues];
            System.arraycopy(newValues, 0, newNewValues, 0, nonZeroValues);
            newValues = newNewValues;
        }
        this.indices = newIndices;
        this.matrixSI = newValues;
        return this;
    }

    @Override
    public final DoubleMatrixDataDense toDense()
    {
        double[] denseSI = new double[this.rows * this.cols];
        for (int index = 0; index < this.matrixSI.length; index++)
        {
            denseSI[(int) this.indices[index]] = this.matrixSI[index];
        }
        try
        {
            return new DoubleMatrixDataDense(denseSI, this.rows, this.cols);
        }
        catch (ValueRuntimeException exception)
        {
            throw new RuntimeException(exception); // cannot happen -- denseSI has the right size
        }
    }

    @Override
    public final DoubleMatrixDataSparse toSparse()
    {
        return this;
    }

    @Override
    public final double getSI(final int row, final int col)
    {
        long index = row * this.cols + col;
        int internalIndex = Arrays.binarySearch(this.indices, index);
        return internalIndex < 0 ? 0.0 : this.matrixSI[internalIndex];
    }

    @Override
    public final void setSI(final int row, final int col, final double valueSI)
    {
        long index = row * this.cols + col;
        int internalIndex = Arrays.binarySearch(this.indices, index);
        if (internalIndex >= 0)
        {
            this.matrixSI[internalIndex] = valueSI;
            return;
        }

        // make room. TODO increase size in chunks
        internalIndex = -internalIndex - 1;
        long[] indicesNew = new long[this.indices.length + 1];
        double[] matrixSINew = new double[this.matrixSI.length + 1];
        System.arraycopy(this.indices, 0, indicesNew, 0, internalIndex);
        System.arraycopy(this.matrixSI, 0, matrixSINew, 0, internalIndex);
        System.arraycopy(this.indices, internalIndex, indicesNew, internalIndex + 1, this.indices.length - internalIndex);
        System.arraycopy(this.matrixSI, internalIndex, matrixSINew, internalIndex + 1, this.indices.length - internalIndex);
        indicesNew[internalIndex] = index;
        matrixSINew[internalIndex] = valueSI;
        this.indices = indicesNew;
        this.matrixSI = matrixSINew;
    }

    @Override
    public final double[][] getDenseMatrixSI()
    {
        return toDense().getDenseMatrixSI();
    }

    @Override
    public final DoubleMatrixDataSparse copy()
    {
        double[] vCopy = new double[this.matrixSI.length];
        System.arraycopy(this.matrixSI, 0, vCopy, 0, this.matrixSI.length);
        long[] iCopy = new long[this.indices.length];
        System.arraycopy(this.indices, 0, iCopy, 0, this.indices.length);
        return new DoubleMatrixDataSparse(vCopy, iCopy, this.rows, this.cols);
    }

    /**
     * Instantiate a DoubleMatrixDataSparse from an array.
     * @param valuesSI double[][]; the (SI) values to store
     * @return the DoubleMatrixDataSparse
     * @throws ValueRuntimeException in case matrix is ragged
     */
    public static DoubleMatrixDataSparse instantiate(final double[][] valuesSI) throws ValueRuntimeException
    {
        checkRectangularAndNonNull(valuesSI);
        int length = nonZero(valuesSI);
        int rows = valuesSI.length;
        final int cols = rows == 0 ? 0 : valuesSI[0].length;
        if (cols == 0)
        {
            rows = 0;
        }
        double[] sparseSI = new double[length];
        long[] indices = new long[length];
        fill(valuesSI, sparseSI, indices);
        return new DoubleMatrixDataSparse(sparseSI, indices, rows, cols);
    }

    /**
     * Instantiate a DoubleMatrixDataSparse from an array.
     * @param values double[][]; the values to store
     * @param scale Scale; the scale that will convert values to SI
     * @return the DoubleMatrixDataSparse
     * @throws ValueRuntimeException in case matrix is ragged
     */
    public static DoubleMatrixDataSparse instantiate(final double[][] values, final Scale scale) throws ValueRuntimeException
    {
        checkRectangularAndNonNull(values);
        int length = nonZero(values);
        int rows = values.length;
        final int cols = rows == 0 ? 0 : values[0].length;
        if (cols == 0)
        {
            rows = 0;
        }
        double[] sparseSI = new double[length];
        long[] indices = new long[length];
        fill(values, sparseSI, indices, scale);
        return new DoubleMatrixDataSparse(sparseSI, indices, rows, cols);
    }

    /**
     * Calculate the number of non-zero values in a double[][] matrix.
     * @param valuesSI double[][]; the double[][] matrix
     * @return the number of non-zero values in the double[][] matrix
     */
    private static int nonZero(final double[][] valuesSI)
    {
        // determine number of non-null cells
        AtomicInteger atomicLength = new AtomicInteger(0);
        IntStream.range(0, valuesSI.length).parallel().forEach(r -> IntStream.range(0, valuesSI[0].length).forEach(c ->
        {
            if (valuesSI[r][c] != 0.0)
            {
                atomicLength.incrementAndGet();
            }
        }));

        return atomicLength.get();
    }

    @Override
    public DoubleMatrixData plus(final DoubleMatrixData right) throws ValueRuntimeException
    {
        if (right.isDense())
        {
            return right.copy().incrementBy(this);
        }
        return this.copy().incrementBy(right);
    }

    @Override
    public final DoubleMatrixData minus(final DoubleMatrixData right)
    {
        if (right.isDense())
        {
            return this.toDense().decrementBy(right);
        }
        return this.copy().decrementBy(right);
    }

    @Override
    public DoubleMatrixDataSparse times(final DoubleMatrixData right) throws ValueRuntimeException
    {
        checkSizes(right);
        DoubleMatrixDataSparse result = this.copy();
        result.multiplyBy(right);
        return result;
    }

    @Override
    public DoubleMatrixData divide(final DoubleMatrixData right) throws ValueRuntimeException
    {
        if (right.isSparse())
        {
            // Sparse divided by sparse makes a dense
            return this.toDense().divide(right);
        }
        // Sparse divided by dense makes a sparse
        checkSizes(right);
        return this.copy().divideBy(right);
    }

    @Override
    public int hashCode()
    {
        return super.hashCode();
    }

    @Override
    @SuppressWarnings({"checkstyle:needbraces", "checkstyle:designforextension"})
    public boolean equals(final Object obj)
    {
        if (this == obj)
            return true;
        if (obj == null)
            return false;
        if (!(obj instanceof DoubleMatrixData))
            return false;
        DoubleMatrixData other = (DoubleMatrixData) obj;
        if (this.rows != other.rows)
            return false;
        if (this.cols != other.cols)
            return false;
        if (other instanceof DoubleMatrixDataDense)
            return super.equals(other);
        if (getClass() != obj.getClass())
            return false;
        // Both are sparse
        if (!Arrays.equals(this.indices, ((DoubleMatrixDataSparse) other).indices))
            return false;
        return Arrays.equals(this.matrixSI, ((DoubleMatrixDataSparse) other).matrixSI);
    }

    @Override
    public String toString()
    {
        return "DoubleMatrixDataSparse [storageType=" + getStorageType() + ", indices=" + Arrays.toString(this.indices)
                + ", matrixSI=" + Arrays.toString(this.matrixSI) + "]";
    }

}