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 <tt>index = row * cols + col</tt>
* @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<DoubleSparseValue<U, S>>; 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) + "]";
}
}