View Javadoc
1   package org.djunits.value.vdouble.matrix.data;
2   
3   import java.util.Arrays;
4   import java.util.Collection;
5   import java.util.concurrent.atomic.AtomicInteger;
6   import java.util.stream.IntStream;
7   
8   import org.djunits.unit.Unit;
9   import org.djunits.unit.scale.Scale;
10  import org.djunits.value.ValueRuntimeException;
11  import org.djunits.value.storage.StorageType;
12  import org.djunits.value.vdouble.function.DoubleFunction;
13  import org.djunits.value.vdouble.function.DoubleFunction2;
14  import org.djunits.value.vdouble.matrix.base.DoubleSparseValue;
15  import org.djunits.value.vdouble.scalar.base.DoubleScalar;
16  import org.djutils.exceptions.Throw;
17  
18  /**
19   * Stores sparse data for a DoubleMatrix and carries out basic operations. The index in the sparse matrix data is calculated as
20   * <code>r * columns + c</code>, where r is the row number, cols is the total number of columns, and c is the column number.
21   * <p>
22   * Copyright (c) 2013-2024 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
23   * BSD-style license. See <a href="https://djunits.org/docs/license.html">DJUNITS License</a>.
24   * </p>
25   * @author <a href="https://www.tudelft.nl/averbraeck">Alexander Verbraeck</a>
26   * @author <a href="https://www.tudelft.nl/staff/p.knoppers/">Peter Knoppers</a>
27   */
28  public class DoubleMatrixDataSparse extends DoubleMatrixData
29  {
30      /** */
31      private static final long serialVersionUID = 1L;
32  
33      /** the index values of the Matrix. The index is calculated as r * columns + c. */
34      private long[] indices;
35  
36      /**
37       * Create a matrix with sparse data.
38       * @param matrixSI double[]; the data to store
39       * @param indices long[]; the index values of the Matrix, with &lt;tt&gt;index = row * cols + col&lt;/tt&gt;
40       * @param rows int; the number of rows
41       * @param cols int; the number of columns
42       */
43      public DoubleMatrixDataSparse(final double[] matrixSI, final long[] indices, final int rows, final int cols)
44      {
45          super(StorageType.SPARSE);
46          this.matrixSI = matrixSI;
47          this.indices = indices;
48          this.rows = rows;
49          this.cols = cols;
50      }
51  
52      /**
53       * Create a matrix with sparse data.
54       * @param dataSI Collection&lt;DoubleSparseValue&lt;U, S&gt;&gt;; the sparse [X, Y, SI] values to store
55       * @param rows int; the number of rows of the matrix
56       * @param cols int; the number of columns of the matrix
57       * @throws NullPointerException when storageType is null or dataSI is null
58       * @param <U> the unit
59       * @param <S> the corresponding scalar type
60       */
61      public <U extends Unit<U>, S extends DoubleScalar<U, S>> DoubleMatrixDataSparse(
62              final Collection<DoubleSparseValue<U, S>> dataSI, final int rows, final int cols) throws NullPointerException
63      {
64          super(StorageType.SPARSE);
65          Throw.whenNull(dataSI, "matrixSI is null");
66  
67          int length = (int) dataSI.stream().parallel().filter(d -> d.getValueSI() != 0.0).count();
68          this.rows = rows;
69          this.cols = cols;
70          this.matrixSI = new double[length];
71          this.indices = new long[length];
72          int index = 0;
73          for (DoubleSparseValue<U, S> data : dataSI)
74          {
75              if (data.getValueSI() != 0.0)
76              {
77                  this.indices[index] = data.getRow() * this.cols + data.getColumn();
78                  this.matrixSI[index] = data.getValueSI();
79                  index++;
80              }
81          }
82      }
83  
84      @Override
85      public final int cardinality()
86      {
87          return this.indices.length;
88      }
89  
90      /**
91       * Fill the sparse data structures matrixSI[] and indices[]. Note: output vectors have to be initialized at the right size.
92       * Cannot be parallelized because of stateful and sequence-sensitive count.
93       * @param data double[][]; the input data
94       * @param matrixSI double[]; the matrix data to write
95       * @param indices long[]; the indices to write
96       * @throws ValueRuntimeException in case matrix is ragged
97       */
98      @SuppressWarnings("checkstyle:finalparameters")
99      private static void fill(final double[][] data, double[] matrixSI, long[] indices) throws ValueRuntimeException
100     {
101         int rows = data.length;
102         int cols = rows == 0 ? 0 : data[0].length;
103         if (cols == 0)
104         {
105             rows = 0;
106         }
107         int count = 0;
108         for (int r = 0; r < rows; r++)
109         {
110             double[] row = data[r];
111             // Row length check has been done by checkRectangularAndNonEmpty
112             for (int c = 0; c < cols; c++)
113             {
114                 int index = r * cols + c;
115                 if (row[c] != 0.0)
116                 {
117                     matrixSI[count] = row[c];
118                     indices[count] = index;
119                     count++;
120                 }
121             }
122         }
123     }
124 
125     /**
126      * Fill the sparse data structures matrixSI[] and indices[]. Note: output vectors have to be initialized at the right size.
127      * Cannot be parallelized because of stateful and sequence-sensitive count.
128      * @param data double[][]; the input data
129      * @param matrixSI double[]; the matrix data to write
130      * @param indices long[]; the indices to write
131      * @param scale Scale; Scale, the scale that will convert the data to SI
132      * @throws ValueRuntimeException in case matrix is ragged
133      */
134     @SuppressWarnings("checkstyle:finalparameters")
135     private static void fill(final double[][] data, double[] matrixSI, long[] indices, final Scale scale)
136             throws ValueRuntimeException
137     {
138         int rows = data.length;
139         int cols = rows == 0 ? 0 : data[0].length;
140         if (cols == 0)
141         {
142             rows = 0;
143         }
144         int count = 0;
145         for (int r = 0; r < rows; r++)
146         {
147             double[] row = data[r];
148             // Row length check has been done by checkRectangularAndNonEmpty
149             for (int c = 0; c < cols; c++)
150             {
151                 int index = r * cols + c;
152                 double value = scale.toStandardUnit(row[c]);
153                 if (value != 0.0)
154                 {
155                     matrixSI[count] = value;
156                     indices[count] = index;
157                     count++;
158                 }
159             }
160         }
161     }
162 
163     @Override
164     public DoubleMatrixData assign(final DoubleFunction doubleFunction)
165     {
166         if (doubleFunction.apply(0d) != 0d)
167         {
168             // It is most unlikely that the result AND the left and right operands are efficiently stored in Sparse format
169             DoubleMatrixDataSparse result = toDense().assign(doubleFunction).toSparse();
170             this.indices = result.indices;
171             this.matrixSI = result.matrixSI;
172             return this;
173         }
174         // The code below relies on the fact that doubleFunction.apply(0d) yields 0d
175         int currentSize = rows() * cols();
176         if (currentSize > 16)
177         {
178             currentSize = 16;
179         }
180         long[] newIndices = new long[currentSize];
181         double[] newValues = new double[currentSize];
182         int nonZeroValues = 0;
183         int ownIndex = 0;
184         while (ownIndex < this.indices.length)
185         {
186             long index = this.indices[ownIndex];
187             double value = doubleFunction.apply(this.matrixSI[ownIndex]);
188             ownIndex++;
189             if (value != 0d)
190             {
191                 if (nonZeroValues >= currentSize)
192                 {
193                     // increase size of arrays
194                     currentSize *= 2;
195                     if (currentSize > rows() * cols())
196                     {
197                         currentSize = rows() * cols();
198                     }
199                     long[] newNewIndices = new long[currentSize];
200                     System.arraycopy(newIndices, 0, newNewIndices, 0, newIndices.length);
201                     newIndices = newNewIndices;
202                     double[] newNewValues = new double[currentSize];
203                     System.arraycopy(newValues, 0, newNewValues, 0, newValues.length);
204                     newValues = newNewValues;
205                 }
206                 newIndices[nonZeroValues] = index;
207                 newValues[nonZeroValues] = value;
208                 nonZeroValues++;
209             }
210         }
211         if (nonZeroValues < currentSize)
212         {
213             // reduce size of arrays
214             long[] newNewIndices = new long[nonZeroValues];
215             System.arraycopy(newIndices, 0, newNewIndices, 0, nonZeroValues);
216             newIndices = newNewIndices;
217             double[] newNewValues = new double[nonZeroValues];
218             System.arraycopy(newValues, 0, newNewValues, 0, nonZeroValues);
219             newValues = newNewValues;
220         }
221         this.indices = newIndices;
222         this.matrixSI = newValues;
223         return this;
224     }
225 
226     @Override
227     public final DoubleMatrixDataSparse assign(final DoubleFunction2 doubleFunction, final DoubleMatrixData right)
228     {
229         if (doubleFunction.apply(0d, 0d) != 0d)
230         {
231             // It is most unlikely that the result AND the left and right operands are efficiently stored in Sparse format
232             DoubleMatrixDataSparse result = toDense().assign(doubleFunction, right).toSparse();
233             this.indices = result.indices;
234             this.matrixSI = result.matrixSI;
235             return this;
236         }
237         // The code below relies on the fact that doubleFunction.apply(0d, 0d) yields 0d
238         checkSizes(right);
239         int currentSize = rows() * cols();
240         if (currentSize > 16)
241         {
242             currentSize = 16;
243         }
244         long[] newIndices = new long[currentSize];
245         double[] newValues = new double[currentSize];
246         int nonZeroValues = 0;
247         int ownIndex = 0;
248         int otherIndex = 0;
249         if (right.isSparse())
250         { // both are sparse, result must be sparse
251             DoubleMatrixDataSparse other = (DoubleMatrixDataSparse) right;
252             while (ownIndex < this.indices.length || otherIndex < other.indices.length)
253             {
254                 double value;
255                 long index;
256                 if (ownIndex < this.indices.length && otherIndex < other.indices.length)
257                 { // neither we nor right has run out of values
258                     if (this.indices[ownIndex] == other.indices[otherIndex])
259                     {
260                         value = doubleFunction.apply(this.matrixSI[ownIndex], other.matrixSI[otherIndex]);
261                         index = this.indices[ownIndex];
262                         ownIndex++;
263                         otherIndex++;
264                     }
265                     else if (this.indices[ownIndex] < other.indices[otherIndex])
266                     {
267                         // we have a non-zero; right has a zero
268                         value = doubleFunction.apply(this.matrixSI[ownIndex], 0d);
269                         index = this.indices[ownIndex];
270                         ownIndex++;
271                     }
272                     else
273                     {
274                         // we have a zero; right has a non-zero
275                         value = doubleFunction.apply(0d, other.matrixSI[otherIndex]);
276                         index = other.indices[otherIndex];
277                         otherIndex++;
278                     }
279                 }
280                 else if (ownIndex < this.indices.length)
281                 { // right has run out of values; we have not
282                     value = doubleFunction.apply(this.matrixSI[ownIndex], 0d);
283                     index = this.indices[ownIndex];
284                     ownIndex++;
285                 }
286                 else
287                 { // we have run out of values; right has not
288                     value = doubleFunction.apply(0d, other.matrixSI[otherIndex]);
289                     index = other.indices[otherIndex];
290                     otherIndex++;
291                 }
292                 if (value != 0d)
293                 {
294                     if (nonZeroValues >= currentSize)
295                     {
296                         // increase size of arrays
297                         currentSize *= 2;
298                         if (currentSize > rows() * cols())
299                         {
300                             currentSize = rows() * cols();
301                         }
302                         long[] newNewIndices = new long[currentSize];
303                         System.arraycopy(newIndices, 0, newNewIndices, 0, newIndices.length);
304                         newIndices = newNewIndices;
305                         double[] newNewValues = new double[currentSize];
306                         System.arraycopy(newValues, 0, newNewValues, 0, newValues.length);
307                         newValues = newNewValues;
308                     }
309                     newIndices[nonZeroValues] = index;
310                     newValues[nonZeroValues] = value;
311                     nonZeroValues++;
312                 }
313             }
314         }
315         else
316         { // we are sparse; other is dense, result must be sparse
317             DoubleMatrixDataDense other = (DoubleMatrixDataDense) right;
318             while (otherIndex < right.matrixSI.length)
319             {
320                 double value;
321                 int index = otherIndex;
322                 if (ownIndex < this.indices.length)
323                 { // neither we nor right has run out of values
324                     if (this.indices[ownIndex] == otherIndex)
325                     {
326                         value = doubleFunction.apply(this.matrixSI[ownIndex], other.matrixSI[otherIndex]);
327                         ownIndex++;
328                     }
329                     else
330                     {
331                         // we have a zero; other has a value
332                         value = doubleFunction.apply(0d, other.matrixSI[otherIndex]);
333                     }
334                     otherIndex++;
335                 }
336                 else
337                 { // we have run out of values; right has not
338                     value = doubleFunction.apply(0d, other.matrixSI[otherIndex]);
339                     otherIndex++;
340                 }
341                 if (value != 0d)
342                 {
343                     if (nonZeroValues >= currentSize)
344                     {
345                         // increase size of arrays
346                         currentSize *= 2;
347                         if (currentSize > rows() * cols())
348                         {
349                             currentSize = rows() * cols();
350                         }
351                         long[] newNewIndices = new long[currentSize];
352                         System.arraycopy(newIndices, 0, newNewIndices, 0, newIndices.length);
353                         newIndices = newNewIndices;
354                         double[] newNewValues = new double[currentSize];
355                         System.arraycopy(newValues, 0, newNewValues, 0, newValues.length);
356                         newValues = newNewValues;
357                     }
358                     newIndices[nonZeroValues] = index;
359                     newValues[nonZeroValues] = value;
360                     nonZeroValues++;
361                 }
362             }
363         }
364         if (nonZeroValues < currentSize)
365         {
366             // reduce size of arrays
367             long[] newNewIndices = new long[nonZeroValues];
368             System.arraycopy(newIndices, 0, newNewIndices, 0, nonZeroValues);
369             newIndices = newNewIndices;
370             double[] newNewValues = new double[nonZeroValues];
371             System.arraycopy(newValues, 0, newNewValues, 0, nonZeroValues);
372             newValues = newNewValues;
373         }
374         this.indices = newIndices;
375         this.matrixSI = newValues;
376         return this;
377     }
378 
379     @Override
380     public final DoubleMatrixDataDense toDense()
381     {
382         double[] denseSI = new double[this.rows * this.cols];
383         for (int index = 0; index < this.matrixSI.length; index++)
384         {
385             denseSI[(int) this.indices[index]] = this.matrixSI[index];
386         }
387         try
388         {
389             return new DoubleMatrixDataDense(denseSI, this.rows, this.cols);
390         }
391         catch (ValueRuntimeException exception)
392         {
393             throw new RuntimeException(exception); // cannot happen -- denseSI has the right size
394         }
395     }
396 
397     @Override
398     public final DoubleMatrixDataSparse toSparse()
399     {
400         return this;
401     }
402 
403     @Override
404     public final double getSI(final int row, final int col)
405     {
406         long index = row * this.cols + col;
407         int internalIndex = Arrays.binarySearch(this.indices, index);
408         return internalIndex < 0 ? 0.0 : this.matrixSI[internalIndex];
409     }
410 
411     @Override
412     public final void setSI(final int row, final int col, final double valueSI)
413     {
414         long index = row * this.cols + col;
415         int internalIndex = Arrays.binarySearch(this.indices, index);
416         if (internalIndex >= 0)
417         {
418             this.matrixSI[internalIndex] = valueSI;
419             return;
420         }
421 
422         // make room. TODO increase size in chunks
423         internalIndex = -internalIndex - 1;
424         long[] indicesNew = new long[this.indices.length + 1];
425         double[] matrixSINew = new double[this.matrixSI.length + 1];
426         System.arraycopy(this.indices, 0, indicesNew, 0, internalIndex);
427         System.arraycopy(this.matrixSI, 0, matrixSINew, 0, internalIndex);
428         System.arraycopy(this.indices, internalIndex, indicesNew, internalIndex + 1, this.indices.length - internalIndex);
429         System.arraycopy(this.matrixSI, internalIndex, matrixSINew, internalIndex + 1, this.indices.length - internalIndex);
430         indicesNew[internalIndex] = index;
431         matrixSINew[internalIndex] = valueSI;
432         this.indices = indicesNew;
433         this.matrixSI = matrixSINew;
434     }
435 
436     @Override
437     public final double[][] getDenseMatrixSI()
438     {
439         return toDense().getDenseMatrixSI();
440     }
441 
442     @Override
443     public final DoubleMatrixDataSparse copy()
444     {
445         double[] vCopy = new double[this.matrixSI.length];
446         System.arraycopy(this.matrixSI, 0, vCopy, 0, this.matrixSI.length);
447         long[] iCopy = new long[this.indices.length];
448         System.arraycopy(this.indices, 0, iCopy, 0, this.indices.length);
449         return new DoubleMatrixDataSparse(vCopy, iCopy, this.rows, this.cols);
450     }
451 
452     /**
453      * Instantiate a DoubleMatrixDataSparse from an array.
454      * @param valuesSI double[][]; the (SI) values to store
455      * @return the DoubleMatrixDataSparse
456      * @throws ValueRuntimeException in case matrix is ragged
457      */
458     public static DoubleMatrixDataSparse instantiate(final double[][] valuesSI) throws ValueRuntimeException
459     {
460         checkRectangularAndNonNull(valuesSI);
461         int length = nonZero(valuesSI);
462         int rows = valuesSI.length;
463         final int cols = rows == 0 ? 0 : valuesSI[0].length;
464         if (cols == 0)
465         {
466             rows = 0;
467         }
468         double[] sparseSI = new double[length];
469         long[] indices = new long[length];
470         fill(valuesSI, sparseSI, indices);
471         return new DoubleMatrixDataSparse(sparseSI, indices, rows, cols);
472     }
473 
474     /**
475      * Instantiate a DoubleMatrixDataSparse from an array.
476      * @param values double[][]; the values to store
477      * @param scale Scale; the scale that will convert values to SI
478      * @return the DoubleMatrixDataSparse
479      * @throws ValueRuntimeException in case matrix is ragged
480      */
481     public static DoubleMatrixDataSparse instantiate(final double[][] values, final Scale scale) throws ValueRuntimeException
482     {
483         checkRectangularAndNonNull(values);
484         int length = nonZero(values);
485         int rows = values.length;
486         final int cols = rows == 0 ? 0 : values[0].length;
487         if (cols == 0)
488         {
489             rows = 0;
490         }
491         double[] sparseSI = new double[length];
492         long[] indices = new long[length];
493         fill(values, sparseSI, indices, scale);
494         return new DoubleMatrixDataSparse(sparseSI, indices, rows, cols);
495     }
496 
497     /**
498      * Calculate the number of non-zero values in a double[][] matrix.
499      * @param valuesSI double[][]; the double[][] matrix
500      * @return the number of non-zero values in the double[][] matrix
501      */
502     private static int nonZero(final double[][] valuesSI)
503     {
504         // determine number of non-null cells
505         AtomicInteger atomicLength = new AtomicInteger(0);
506         IntStream.range(0, valuesSI.length).parallel().forEach(r -> IntStream.range(0, valuesSI[0].length).forEach(c ->
507         {
508             if (valuesSI[r][c] != 0.0)
509             {
510                 atomicLength.incrementAndGet();
511             }
512         }));
513 
514         return atomicLength.get();
515     }
516 
517     @Override
518     public DoubleMatrixData plus(final DoubleMatrixData right) throws ValueRuntimeException
519     {
520         if (right.isDense())
521         {
522             return right.copy().incrementBy(this);
523         }
524         return this.copy().incrementBy(right);
525     }
526 
527     @Override
528     public final DoubleMatrixData minus(final DoubleMatrixData right)
529     {
530         if (right.isDense())
531         {
532             return this.toDense().decrementBy(right);
533         }
534         return this.copy().decrementBy(right);
535     }
536 
537     @Override
538     public DoubleMatrixDataSparse times(final DoubleMatrixData right) throws ValueRuntimeException
539     {
540         checkSizes(right);
541         DoubleMatrixDataSparse result = this.copy();
542         result.multiplyBy(right);
543         return result;
544     }
545 
546     @Override
547     public DoubleMatrixData divide(final DoubleMatrixData right) throws ValueRuntimeException
548     {
549         if (right.isSparse())
550         {
551             // Sparse divided by sparse makes a dense
552             return this.toDense().divide(right);
553         }
554         // Sparse divided by dense makes a sparse
555         checkSizes(right);
556         return this.copy().divideBy(right);
557     }
558 
559     @Override
560     public int hashCode()
561     {
562         return super.hashCode();
563     }
564 
565     @Override
566     @SuppressWarnings({"checkstyle:needbraces", "checkstyle:designforextension"})
567     public boolean equals(final Object obj)
568     {
569         if (this == obj)
570             return true;
571         if (obj == null)
572             return false;
573         if (!(obj instanceof DoubleMatrixData))
574             return false;
575         DoubleMatrixData other = (DoubleMatrixData) obj;
576         if (this.rows != other.rows)
577             return false;
578         if (this.cols != other.cols)
579             return false;
580         if (other instanceof DoubleMatrixDataDense)
581             return super.equals(other);
582         if (getClass() != obj.getClass())
583             return false;
584         // Both are sparse
585         if (!Arrays.equals(this.indices, ((DoubleMatrixDataSparse) other).indices))
586             return false;
587         return Arrays.equals(this.matrixSI, ((DoubleMatrixDataSparse) other).matrixSI);
588     }
589 
590     @Override
591     public String toString()
592     {
593         return "DoubleMatrixDataSparse [storageType=" + getStorageType() + ", indices=" + Arrays.toString(this.indices)
594                 + ", matrixSI=" + Arrays.toString(this.matrixSI) + "]";
595     }
596 
597 }