From 1671a9ec323590d94198c853f5a74b6956243db9 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 26 Apr 2014 22:18:44 +0200 Subject: [PATCH] LA: Matrix MapSubMatrixIndexedTo --- .../LinearAlgebra/Complex/DiagonalMatrix.cs | 23 +- .../LinearAlgebra/Complex/SparseMatrix.cs | 17 +- .../LinearAlgebra/Complex32/DiagonalMatrix.cs | 23 +- .../LinearAlgebra/Complex32/SparseMatrix.cs | 17 +- .../LinearAlgebra/Double/DiagonalMatrix.cs | 23 +- .../LinearAlgebra/Double/SparseMatrix.cs | 17 +- src/Numerics/LinearAlgebra/Matrix.cs | 4 +- .../LinearAlgebra/Single/DiagonalMatrix.cs | 23 +- .../LinearAlgebra/Single/SparseMatrix.cs | 17 +- .../Storage/DenseColumnMajorMatrixStorage.cs | 86 +++- .../Storage/DiagonalMatrixStorage.cs | 295 ++++++++++--- .../Storage/MatrixStorage.Validation.cs | 19 +- .../LinearAlgebra/Storage/MatrixStorage.cs | 87 +++- .../SparseCompressedRowMatrixStorage.cs | 392 ++++++++++++++---- .../MatrixStructureTheory.Map.cs | 260 ++++++++++++ .../MatrixStructureTheory.Reform.cs | 4 +- src/UnitTests/UnitTests.csproj | 1 + 17 files changed, 1048 insertions(+), 260 deletions(-) create mode 100644 src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Map.cs diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs index f796e6e0..a685a329 100644 --- a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs @@ -435,14 +435,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return; } - if (RowCount == ColumnCount) + if (ColumnCount == RowCount) { - // TODO: Map is rather generic, we can do better - other.MapIndexed((i, j, x) => _data[i]*x, result, false); - return; + other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false); + } + else + { + result.Clear(); + other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true); } - - base.DoMultiply(other, result); } /// @@ -577,7 +578,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return; } - base.DoTransposeThisAndMultiply(other, result); + if (ColumnCount == RowCount) + { + other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false); + } + else + { + result.Clear(); + other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true); + } } /// diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs index 0be638d2..b335c960 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs @@ -931,12 +931,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return; } - var diagonalOther = other as DiagonalMatrix; - if (diagonalOther != null && sparseResult != null && other.RowCount == other.ColumnCount) + var diagonalOther = other.Storage as DiagonalMatrixStorage; + if (diagonalOther != null && sparseResult != null) { - var diagonal = ((DiagonalMatrixStorage)other.Storage).Data; - // TODO: Map is rather generic, we can do better - MapIndexed((i, j, x) => x*diagonal[j], result, false); + var diagonal = diagonalOther.Data; + if (other.ColumnCount == other.RowCount) + { + Storage.MapIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], false, false); + } + else + { + result.Storage.Clear(); + Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], 0, 0, RowCount, 0, 0, ColumnCount, false, true); + } return; } diff --git a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs index 2340a8d3..e71f3e38 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs @@ -429,14 +429,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return; } - if (RowCount == ColumnCount) + if (ColumnCount == RowCount) { - // TODO: Map is rather generic, we can do better - other.MapIndexed((i, j, x) => _data[i]*x, result, false); - return; + other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false); + } + else + { + result.Clear(); + other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true); } - - base.DoMultiply(other, result); } /// @@ -571,7 +572,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return; } - base.DoTransposeThisAndMultiply(other, result); + if (ColumnCount == RowCount) + { + other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false); + } + else + { + result.Clear(); + other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true); + } } /// diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs index 5b7f6442..bd5cbb40 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs @@ -925,12 +925,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return; } - var diagonalOther = other as DiagonalMatrix; - if (diagonalOther != null && sparseResult != null && other.RowCount == other.ColumnCount) + var diagonalOther = other.Storage as DiagonalMatrixStorage; + if (diagonalOther != null && sparseResult != null) { - var diagonal = ((DiagonalMatrixStorage)other.Storage).Data; - // TODO: Map is rather generic, we can do better - MapIndexed((i, j, x) => x*diagonal[j], result, false); + var diagonal = diagonalOther.Data; + if (other.ColumnCount == other.RowCount) + { + Storage.MapIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], false, false); + } + else + { + result.Storage.Clear(); + Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], 0, 0, RowCount, 0, 0, ColumnCount, false, true); + } return; } diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs index dc23fc26..ac89caf5 100644 --- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs @@ -409,14 +409,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double return; } - if (RowCount == ColumnCount) + if (ColumnCount == RowCount) { - // TODO: Map is rather generic, we can do better - other.MapIndexed((i, j, x) => _data[i]*x, result, false); - return; + other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false); + } + else + { + result.Clear(); + other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true); } - - base.DoMultiply(other, result); } /// @@ -505,7 +506,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double return; } - base.DoTransposeThisAndMultiply(other, result); + if (ColumnCount == RowCount) + { + other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false); + } + else + { + result.Clear(); + other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true); + } } /// diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs index e8137235..b85afa4a 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs @@ -927,12 +927,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double return; } - var diagonalOther = other as DiagonalMatrix; - if (diagonalOther != null && sparseResult != null && other.RowCount == other.ColumnCount) + var diagonalOther = other.Storage as DiagonalMatrixStorage; + if (diagonalOther != null && sparseResult != null) { - var diagonal = ((DiagonalMatrixStorage)other.Storage).Data; - // TODO: Map is rather generic, we can do better - MapIndexed((i, j, x) => x*diagonal[j], result, false); + var diagonal = diagonalOther.Data; + if (other.ColumnCount == other.RowCount) + { + Storage.MapIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], false, false); + } + else + { + result.Storage.Clear(); + Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], 0, 0, RowCount, 0, 0, ColumnCount, false, true); + } return; } diff --git a/src/Numerics/LinearAlgebra/Matrix.cs b/src/Numerics/LinearAlgebra/Matrix.cs index 991b4664..54190a89 100644 --- a/src/Numerics/LinearAlgebra/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Matrix.cs @@ -1518,7 +1518,7 @@ namespace MathNet.Numerics.LinearAlgebra public Matrix Map(Func f, bool forceMapZeros = false) where TU : struct, IEquatable, IFormattable { - var result = Matrix.Build.SameAs(this); + var result = Matrix.Build.SameAs(this, RowCount, ColumnCount, fullyMutable: forceMapZeros); Storage.MapToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true); return result; } @@ -1532,7 +1532,7 @@ namespace MathNet.Numerics.LinearAlgebra public Matrix MapIndexed(Func f, bool forceMapZeros = false) where TU : struct, IEquatable, IFormattable { - var result = Matrix.Build.SameAs(this); + var result = Matrix.Build.SameAs(this, RowCount, ColumnCount, fullyMutable: forceMapZeros); Storage.MapIndexedToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true); return result; } diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs index 0492dd1b..30f77e45 100644 --- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs @@ -409,14 +409,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single return; } - if (RowCount == ColumnCount) + if (ColumnCount == RowCount) { - // TODO: Map is rather generic, we can do better - other.MapIndexed((i, j, x) => _data[i]*x, result, false); - return; + other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false); + } + else + { + result.Clear(); + other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true); } - - base.DoMultiply(other, result); } /// @@ -505,7 +506,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single return; } - base.DoTransposeThisAndMultiply(other, result); + if (ColumnCount == RowCount) + { + other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false); + } + else + { + result.Clear(); + other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true); + } } /// diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs index c3c80fe5..5a93f58a 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs @@ -932,12 +932,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single return; } - var diagonalOther = other as DiagonalMatrix; - if (diagonalOther != null && sparseResult != null && other.RowCount == other.ColumnCount) + var diagonalOther = other.Storage as DiagonalMatrixStorage; + if (diagonalOther != null && sparseResult != null) { - var diagonal = ((DiagonalMatrixStorage)other.Storage).Data; - // TODO: Map is rather generic, we can do better - MapIndexed((i, j, x) => x*diagonal[j], result, false); + var diagonal = diagonalOther.Data; + if (other.ColumnCount == other.RowCount) + { + Storage.MapIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], false, false); + } + else + { + result.Storage.Clear(); + Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], 0, 0, RowCount, 0, 0, ColumnCount, false, true); + } return; } diff --git a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs index 9c7c5529..d85298f4 100644 --- a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs @@ -365,9 +365,18 @@ namespace MathNet.Numerics.LinearAlgebra.Storage return; } + // TODO: Proper Sparse Implementation + // FALL BACK - base.CopySubMatrixToUnchecked(target, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, skipClearing); + for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnIndex + columnCount; j++, jj++) + { + int index = sourceRowIndex + j*RowCount; + for (int ii = targetRowIndex; ii < targetRowIndex + rowCount; ii++) + { + target.At(ii, jj, Data[index++]); + } + } } void CopySubMatrixToUnchecked(DenseColumnMajorMatrixStorage target, @@ -503,6 +512,33 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS + public override void MapInplace(Func f, bool forceMapZeros = false) + { + CommonParallel.For(0, Data.Length, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + Data[i] = f(Data[i]); + } + }); + } + + public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + { + CommonParallel.For(0, ColumnCount, Math.Max(4096/RowCount, 32), (a, b) => + { + int index = a*RowCount; + for (int j = a; j < b; j++) + { + for (int i = 0; i < RowCount; i++) + { + Data[index] = f(i, j, Data[index]); + index++; + } + } + }); + } + internal override void MapToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { var denseTarget = target as DenseColumnMajorMatrixStorage; @@ -530,17 +566,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - public override void MapInplace(Func f, bool forceMapZeros = false) - { - CommonParallel.For(0, Data.Length, 4096, (a, b) => - { - for (int i = a; i < b; i++) - { - Data[i] = f(Data[i]); - } - }); - } - internal override void MapIndexedToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { var denseTarget = target as DenseColumnMajorMatrixStorage; @@ -573,20 +598,41 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + internal override void MapSubMatrixIndexedToUnchecked(MatrixStorage target, Func f, + int sourceRowIndex, int targetRowIndex, int rowCount, + int sourceColumnIndex, int targetColumnIndex, int columnCount, + bool forceMapZeros = false, bool skipClearing = false) { - CommonParallel.For(0, ColumnCount, Math.Max(4096/RowCount, 32), (a, b) => + var denseTarget = target as DenseColumnMajorMatrixStorage; + if (denseTarget != null) { - int index = a*RowCount; - for (int j = a; j < b; j++) + CommonParallel.For(0, columnCount, Math.Max(4096/rowCount, 32), (a, b) => { - for (int i = 0; i < RowCount; i++) + for (int j = a; j < b; j++) { - Data[index] = f(i, j, Data[index]); - index++; + int sourceIndex = sourceRowIndex + (j + sourceColumnIndex)*RowCount; + int targetIndex = targetRowIndex + (j + targetColumnIndex)*target.RowCount; + for (int i = 0; i < rowCount; i++) + { + denseTarget.Data[targetIndex++] = f(targetRowIndex + i, targetColumnIndex + j, Data[sourceIndex++]); + } } + }); + return; + } + + // TODO: Proper Sparse Implementation + + // FALL BACK + + for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnIndex + columnCount; j++, jj++) + { + int index = sourceRowIndex + j*RowCount; + for (int ii = targetRowIndex; ii < targetRowIndex + rowCount; ii++) + { + target.At(ii, jj, f(ii, jj, Data[index++])); } - }); + } } } } diff --git a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs index 4869ca39..1e43ae21 100644 --- a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs @@ -349,16 +349,40 @@ namespace MathNet.Numerics.LinearAlgebra.Storage return; } - var sparseTarget = target as SparseCompressedRowMatrixStorage; - if (sparseTarget != null) - { - CopySubMatrixToUnchecked(sparseTarget, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, skipClearing); - return; - } + // TODO: Proper Sparse Implementation // FALL BACK - base.CopySubMatrixToUnchecked(target, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, skipClearing); + if (!skipClearing) + { + target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount); + } + + if (sourceRowIndex == sourceColumnIndex) + { + for (var i = 0; i < Math.Min(columnCount, rowCount); i++) + { + target.At(targetRowIndex + i, targetColumnIndex + i, Data[sourceRowIndex + i]); + } + } + else if (sourceRowIndex > sourceColumnIndex && sourceColumnIndex + columnCount > sourceRowIndex) + { + // column by column, but skip resulting zero columns at the beginning + int columnInit = sourceRowIndex - sourceColumnIndex; + for (var i = 0; i < Math.Min(columnCount - columnInit, rowCount); i++) + { + target.At(targetRowIndex + i, columnInit + targetColumnIndex + i, Data[sourceRowIndex + i]); + } + } + else if (sourceRowIndex < sourceColumnIndex && sourceRowIndex + rowCount > sourceColumnIndex) + { + // row by row, but skip resulting zero rows at the beginning + int rowInit = sourceColumnIndex - sourceRowIndex; + for (var i = 0; i < Math.Min(columnCount, rowCount - rowInit); i++) + { + target.At(rowInit + targetRowIndex + i, targetColumnIndex + i, Data[sourceColumnIndex + i]); + } + } } void CopySubMatrixToUnchecked(DiagonalMatrixStorage target, @@ -436,45 +460,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - void CopySubMatrixToUnchecked(SparseCompressedRowMatrixStorage target, - int sourceRowIndex, int targetRowIndex, int rowCount, - int sourceColumnIndex, int targetColumnIndex, int columnCount, - bool skipClearing) - { - if (!skipClearing) - { - target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount); - } - - if (sourceRowIndex == sourceColumnIndex) - { - for (var i = 0; i < Math.Min(columnCount, rowCount); i++) - { - target.At(i + targetRowIndex, i + targetColumnIndex, Data[sourceRowIndex + i]); - } - } - else if (sourceRowIndex > sourceColumnIndex && sourceColumnIndex + columnCount > sourceRowIndex) - { - // column by column, but skip resulting zero columns at the beginning - int columnInit = sourceRowIndex - sourceColumnIndex; - for (var i = 0; i < Math.Min(columnCount - columnInit, rowCount); i++) - { - target.At(i + targetRowIndex, columnInit + i + targetColumnIndex, Data[sourceRowIndex + i]); - } - } - else if (sourceRowIndex < sourceColumnIndex && sourceRowIndex + rowCount > sourceColumnIndex) - { - // row by row, but skip resulting zero rows at the beginning - int rowInit = sourceColumnIndex - sourceRowIndex; - for (var i = 0; i < Math.Min(columnCount, rowCount - rowInit); i++) - { - target.At(rowInit + i + targetRowIndex, i + targetColumnIndex, Data[sourceColumnIndex + i]); - } - } - - // else: all zero, nop - } - // ROW COPY internal override void CopySubRowToUnchecked(VectorStorage target, int rowIndex, @@ -589,6 +574,38 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS + public override void MapInplace(Func f, bool forceMapZeros = false) + { + if (forceMapZeros) + { + throw new NotSupportedException("Cannot map non-zero off-diagonal values into a diagonal matrix"); + } + + CommonParallel.For(0, Data.Length, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + Data[i] = f(Data[i]); + } + }); + } + + public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + { + if (forceMapZeros) + { + throw new NotSupportedException("Cannot map non-zero off-diagonal values into a diagonal matrix"); + } + + CommonParallel.For(0, Data.Length, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + Data[i] = f(i, i, Data[i]); + } + }); + } + internal override void MapToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { var processZeros = forceMapZeros || !Zero.Equals(f(Zero)); @@ -637,25 +654,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - public override void MapInplace(Func f, bool forceMapZeros = false) - { - if (forceMapZeros) - { - throw new NotSupportedException("Cannot map non-zero off-diagonal values into a diagonal matrix"); - } - - CommonParallel.For(0, Data.Length, 4096, (a, b) => - { - for (int i = a; i < b; i++) - { - Data[i] = f(Data[i]); - } - }); - } - internal override void MapIndexedToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - var processZeros = forceMapZeros || !Zero.Equals(f(0, 0, Zero)); + var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero)); var diagonalTarget = target as DiagonalMatrixStorage; if (diagonalTarget != null) @@ -701,20 +702,176 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + internal override void MapSubMatrixIndexedToUnchecked(MatrixStorage target, Func f, + int sourceRowIndex, int targetRowIndex, int rowCount, + int sourceColumnIndex, int targetColumnIndex, int columnCount, + bool forceMapZeros = false, bool skipClearing = false) { - if (forceMapZeros) + var diagonalTarget = target as DiagonalMatrixStorage; + if (diagonalTarget != null) + { + MapSubMatrixIndexedToUnchecked(diagonalTarget, f, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, forceMapZeros); + return; + } + + var denseTarget = target as DenseColumnMajorMatrixStorage; + if (denseTarget != null) + { + MapSubMatrixIndexedToUnchecked(denseTarget, f, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, forceMapZeros, skipClearing); + return; + } + + // TODO: Proper Sparse Implementation + + // FALL BACK + + if (!skipClearing) + { + target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount); + } + + if (sourceRowIndex == sourceColumnIndex) + { + int targetRow = targetRowIndex; + int targetColumn = targetColumnIndex; + for (var i = 0; i < Math.Min(columnCount, rowCount); i++) + { + target.At(targetRow, targetColumn, f(targetRow, targetColumn, Data[sourceRowIndex + i])); + targetRow++; + targetColumn++; + } + } + else if (sourceRowIndex > sourceColumnIndex && sourceColumnIndex + columnCount > sourceRowIndex) + { + // column by column, but skip resulting zero columns at the beginning + int columnInit = sourceRowIndex - sourceColumnIndex; + int targetRow = targetRowIndex; + int targetColumn = targetColumnIndex + columnInit; + for (var i = 0; i < Math.Min(columnCount - columnInit, rowCount); i++) + { + target.At(targetRow, targetColumn, f(targetRow, targetColumn, Data[sourceRowIndex + i])); + targetRow++; + targetColumn++; + } + } + else if (sourceRowIndex < sourceColumnIndex && sourceRowIndex + rowCount > sourceColumnIndex) + { + // row by row, but skip resulting zero rows at the beginning + int rowInit = sourceColumnIndex - sourceRowIndex; + int targetRow = targetRowIndex + rowInit; + int targetColumn = targetColumnIndex; + for (var i = 0; i < Math.Min(columnCount, rowCount - rowInit); i++) + { + target.At(targetRow, targetColumn, f(targetRow, targetColumn, Data[sourceColumnIndex + i])); + targetRow++; + targetColumn++; + } + } + } + + void MapSubMatrixIndexedToUnchecked(DiagonalMatrixStorage target, Func f, + int sourceRowIndex, int targetRowIndex, int rowCount, + int sourceColumnIndex, int targetColumnIndex, int columnCount, + bool forceMapZeros = false) + where TU : struct, IEquatable, IFormattable + { + var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero)); + if (processZeros || sourceRowIndex - sourceColumnIndex != targetRowIndex - targetColumnIndex) { throw new NotSupportedException("Cannot map non-zero off-diagonal values into a diagonal matrix"); } - CommonParallel.For(0, Data.Length, 4096, (a, b) => + var beginInclusive = Math.Max(sourceRowIndex, sourceColumnIndex); + var count = Math.Min(sourceRowIndex + rowCount, sourceColumnIndex + columnCount) - beginInclusive; + if (count > 0) { - for (int i = a; i < b; i++) + var beginTarget = Math.Max(targetRowIndex, targetColumnIndex); + CommonParallel.For(0, count, 4096, (a, b) => { - Data[i] = f(i, i, Data[i]); + int targetIndex = beginTarget + a; + for (int i = a; i < b; i++) + { + target.Data[targetIndex] = f(targetIndex, targetIndex, Data[beginInclusive + i]); + targetIndex++; + } + }); + } + } + + void MapSubMatrixIndexedToUnchecked(DenseColumnMajorMatrixStorage target, Func f, + int sourceRowIndex, int targetRowIndex, int rowCount, + int sourceColumnIndex, int targetColumnIndex, int columnCount, + bool forceMapZeros = false, bool skipClearing = false) + where TU : struct, IEquatable, IFormattable + { + var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero)); + if (!skipClearing && !processZeros) + { + target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount); + } + + if (processZeros) + { + CommonParallel.For(0, columnCount, Math.Max(4096/rowCount, 32), (a, b) => + { + int sourceColumn = sourceColumnIndex + a; + int targetColumn = targetColumnIndex + a; + for (int j = a; j < b; j++) + { + int targetIndex = targetRowIndex + (j + targetColumnIndex)*target.RowCount; + int sourceRow = sourceRowIndex; + int targetRow = targetRowIndex; + for (int i = 0; i < rowCount; i++) + { + target.Data[targetIndex++] = f(targetRow++, targetColumn, sourceRow++ == sourceColumn ? Data[sourceColumn] : Zero); + } + sourceColumn++; + targetColumn++; + } + }); + } + else + { + if (sourceRowIndex > sourceColumnIndex && sourceColumnIndex + columnCount > sourceRowIndex) + { + // column by column, but skip resulting zero columns at the beginning + + int columnInit = sourceRowIndex - sourceColumnIndex; + int offset = (columnInit + targetColumnIndex)*target.RowCount + targetRowIndex; + int step = target.RowCount + 1; + int count = Math.Min(columnCount - columnInit, rowCount); + + for (int k = 0, j = offset; k < count; j += step, k++) + { + target.Data[j] = f(targetRowIndex + k, targetColumnIndex + columnInit + k, Data[sourceRowIndex + k]); + } } - }); + else if (sourceRowIndex < sourceColumnIndex && sourceRowIndex + rowCount > sourceColumnIndex) + { + // row by row, but skip resulting zero rows at the beginning + + int rowInit = sourceColumnIndex - sourceRowIndex; + int offset = targetColumnIndex*target.RowCount + rowInit + targetRowIndex; + int step = target.RowCount + 1; + int count = Math.Min(columnCount, rowCount - rowInit); + + for (int k = 0, j = offset; k < count; j += step, k++) + { + target.Data[j] = f(targetRowIndex + rowInit + k, targetColumnIndex + k, Data[sourceColumnIndex + k]); + } + } + else + { + int offset = targetColumnIndex*target.RowCount + targetRowIndex; + int step = target.RowCount + 1; + var count = Math.Min(columnCount, rowCount); + + for (int k = 0, j = offset; k < count; j += step, k++) + { + target.Data[j] = f(targetRowIndex + k, targetColumnIndex + k, Data[sourceRowIndex + k]); + } + } + } } } } diff --git a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.Validation.cs b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.Validation.cs index 087bad08..8eb3cc63 100644 --- a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.Validation.cs +++ b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.Validation.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2014 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -34,6 +34,7 @@ using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Storage { // ReSharper disable UnusedParameter.Local + public partial class MatrixStorage { void ValidateRange(int row, int column) @@ -49,9 +50,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - void ValidateSubMatrixRange(MatrixStorage target, + void ValidateSubMatrixRange(MatrixStorage target, int sourceRowIndex, int targetRowIndex, int rowCount, int sourceColumnIndex, int targetColumnIndex, int columnCount) + where TU : struct, IEquatable, IFormattable { if (rowCount < 1) { @@ -114,7 +116,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - void ValidateRowRange(VectorStorage target, int rowIndex) + void ValidateRowRange(VectorStorage target, int rowIndex) + where TU : struct, IEquatable, IFormattable { if (rowIndex >= RowCount || rowIndex < 0) { @@ -127,7 +130,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - void ValidateColumnRange(VectorStorage target, int columnIndex) + void ValidateColumnRange(VectorStorage target, int columnIndex) + where TU : struct, IEquatable, IFormattable { if (columnIndex >= ColumnCount || columnIndex < 0) { @@ -140,8 +144,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - void ValidateSubRowRange(VectorStorage target, int rowIndex, + void ValidateSubRowRange(VectorStorage target, int rowIndex, int sourceColumnIndex, int targetColumnIndex, int columnCount) + where TU : struct, IEquatable, IFormattable { if (columnCount < 1) { @@ -178,8 +183,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - void ValidateSubColumnRange(VectorStorage target, int columnIndex, + void ValidateSubColumnRange(VectorStorage target, int columnIndex, int sourceRowIndex, int targetRowIndex, int rowCount) + where TU : struct, IEquatable, IFormattable { if (rowCount < 1) { @@ -216,5 +222,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } } + // ReSharper restore UnusedParameter.Local } diff --git a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs index 98a92bf8..b6971ecd 100644 --- a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs @@ -388,10 +388,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage public virtual T[] ToRowMajorArray() { - var ret = new T[RowCount * ColumnCount]; + var ret = new T[RowCount*ColumnCount]; for (int i = 0; i < RowCount; i++) { - var offset = i * ColumnCount; + var offset = i*ColumnCount; for (int j = 0; j < ColumnCount; j++) { ret[offset + j] = At(i, j); @@ -402,10 +402,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage public virtual T[] ToColumnMajorArray() { - var ret = new T[RowCount * ColumnCount]; + var ret = new T[RowCount*ColumnCount]; for (int j = 0; j < ColumnCount; j++) { - var offset = j * RowCount; + var offset = j*RowCount; for (int i = 0; i < RowCount; i++) { ret[offset + i] = At(i, j); @@ -416,7 +416,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage public virtual T[,] ToArray() { - var ret = new T[RowCount,ColumnCount]; + var ret = new T[RowCount, ColumnCount]; for (int i = 0; i < RowCount; i++) { for (int j = 0; j < ColumnCount; j++) @@ -483,6 +483,28 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS + public virtual void MapInplace(Func f, bool forceMapZeros = false) + { + for (int i = 0; i < RowCount; i++) + { + for (int j = 0; j < ColumnCount; j++) + { + At(i, j, f(At(i, j))); + } + } + } + + public virtual void MapIndexedInplace(Func f, bool forceMapZeros = false) + { + for (int i = 0; i < RowCount; i++) + { + for (int j = 0; j < ColumnCount; j++) + { + At(i, j, f(i, j, At(i, j))); + } + } + } + public void MapTo(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) where TU : struct, IEquatable, IFormattable { @@ -512,17 +534,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - public virtual void MapInplace(Func f, bool forceMapZeros = false) - { - for (int i = 0; i < RowCount; i++) - { - for (int j = 0; j < ColumnCount; j++) - { - At(i, j, f(At(i, j))); - } - } - } - public void MapIndexedTo(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) where TU : struct, IEquatable, IFormattable { @@ -543,22 +554,54 @@ namespace MathNet.Numerics.LinearAlgebra.Storage internal virtual void MapIndexedToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) where TU : struct, IEquatable, IFormattable { - for (int i = 0; i < RowCount; i++) + for (int j = 0; j < ColumnCount; j++) { - for (int j = 0; j < ColumnCount; j++) + for (int i = 0; i < RowCount; i++) { target.At(i, j, f(i, j, At(i, j))); } } } - public virtual void MapIndexedInplace(Func f, bool forceMapZeros = false) + public void MapSubMatrixIndexedTo(MatrixStorage target, Func f, + int sourceRowIndex, int targetRowIndex, int rowCount, + int sourceColumnIndex, int targetColumnIndex, int columnCount, + bool forceMapZeros = false, bool skipClearing = false) + where TU : struct, IEquatable, IFormattable { - for (int i = 0; i < RowCount; i++) + if (target == null) { - for (int j = 0; j < ColumnCount; j++) + throw new ArgumentNullException("target"); + } + + if (rowCount == 0 || columnCount == 0) + { + return; + } + + if (ReferenceEquals(this, target)) + { + throw new NotSupportedException(); + } + + ValidateSubMatrixRange(target, + sourceRowIndex, targetRowIndex, rowCount, + sourceColumnIndex, targetColumnIndex, columnCount); + + MapSubMatrixIndexedToUnchecked(target, f, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, forceMapZeros, skipClearing); + } + + internal virtual void MapSubMatrixIndexedToUnchecked(MatrixStorage target, Func f, + int sourceRowIndex, int targetRowIndex, int rowCount, + int sourceColumnIndex, int targetColumnIndex, int columnCount, + bool forceMapZeros = false, bool skipClearing = false) + where TU : struct, IEquatable, IFormattable + { + for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnIndex + columnCount; j++, jj++) + { + for (int i = sourceRowIndex, ii = targetRowIndex; i < sourceRowIndex + rowCount; i++, ii++) { - At(i, j, f(i, j, At(i, j))); + target.At(ii, jj, f(ii, jj, At(i, j))); } } } diff --git a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs index f7035e95..cd2aa145 100644 --- a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs @@ -926,7 +926,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage var columnIndices = new List(ValueCount); var rowPointers = target.RowPointers; - for (int i = sourceRowIndex, row = 0; i < sourceRowIndex + rowCount; i++, row++) + for (int i = sourceRowIndex; i < sourceRowIndex + rowCount; i++) { rowPointers[i + rowOffset] = values.Count; @@ -934,13 +934,13 @@ namespace MathNet.Numerics.LinearAlgebra.Storage var endIndex = RowPointers[i + 1]; // note: we might be able to replace this loop with Array.Copy (perf) - for (int j = startIndex; j < endIndex; j++) + for (int k = startIndex; k < endIndex; k++) { // check if the column index is in the range - if ((ColumnIndices[j] >= sourceColumnIndex) && (ColumnIndices[j] < sourceColumnIndex + columnCount)) + if ((ColumnIndices[k] >= sourceColumnIndex) && (ColumnIndices[k] < sourceColumnIndex + columnCount)) { - values.Add(Values[j]); - columnIndices.Add(ColumnIndices[j] + columnOffset); + values.Add(Values[k]); + columnIndices.Add(ColumnIndices[k] + columnOffset); } } } @@ -963,7 +963,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } // NOTE: potential for more efficient implementation - for (int i = sourceRowIndex, row = 0; i < sourceRowIndex + rowCount; i++, row++) + for (int i = sourceRowIndex, row = 0; row < rowCount; i++, row++) { var startIndex = RowPointers[i]; var endIndex = RowPointers[i + 1]; @@ -1118,6 +1118,112 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS + public override void MapInplace(Func f, bool forceMapZeros = false) + { + if (forceMapZeros || !Zero.Equals(f(Zero))) + { + var newRowPointers = RowPointers; + var newColumnIndices = new List(ColumnIndices.Length); + var newValues = new List(Values.Length); + + int k = 0; + for (int row = 0; row < RowCount; row++) + { + newRowPointers[row] = newValues.Count; + for (int col = 0; col < ColumnCount; col++) + { + var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(Values[k++]) : f(Zero); + if (!Zero.Equals(item)) + { + newValues.Add(item); + newColumnIndices.Add(col); + } + } + } + + ColumnIndices = newColumnIndices.ToArray(); + Values = newValues.ToArray(); + newRowPointers[RowCount] = newValues.Count; + } + else + { + // we can safely do this in-place: + int nonZero = 0; + for (int row = 0; row < RowCount; row++) + { + var startIndex = RowPointers[row]; + var endIndex = RowPointers[row + 1]; + RowPointers[row] = nonZero; + for (var j = startIndex; j < endIndex; j++) + { + var item = f(Values[j]); + if (!Zero.Equals(item)) + { + Values[nonZero] = item; + ColumnIndices[nonZero] = ColumnIndices[j]; + nonZero++; + } + } + } + Array.Resize(ref ColumnIndices, nonZero); + Array.Resize(ref Values, nonZero); + RowPointers[RowCount] = nonZero; + } + } + + public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + { + if (forceMapZeros || !Zero.Equals(f(0, 1, Zero))) + { + var newRowPointers = RowPointers; + var newColumnIndices = new List(ColumnIndices.Length); + var newValues = new List(Values.Length); + + int k = 0; + for (int row = 0; row < RowCount; row++) + { + newRowPointers[row] = newValues.Count; + for (int col = 0; col < ColumnCount; col++) + { + var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(row, col, Values[k++]) : f(row, col, Zero); + if (!Zero.Equals(item)) + { + newValues.Add(item); + newColumnIndices.Add(col); + } + } + } + + ColumnIndices = newColumnIndices.ToArray(); + Values = newValues.ToArray(); + newRowPointers[RowCount] = newValues.Count; + } + else + { + // we can safely do this in-place: + int nonZero = 0; + for (int row = 0; row < RowCount; row++) + { + var startIndex = RowPointers[row]; + var endIndex = RowPointers[row + 1]; + RowPointers[row] = nonZero; + for (var j = startIndex; j < endIndex; j++) + { + var item = f(row, ColumnIndices[j], Values[j]); + if (!Zero.Equals(item)) + { + Values[nonZero] = item; + ColumnIndices[nonZero] = ColumnIndices[j]; + nonZero++; + } + } + } + Array.Resize(ref ColumnIndices, nonZero); + Array.Resize(ref Values, nonZero); + RowPointers[RowCount] = nonZero; + } + } + internal override void MapToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { var processZeros = forceMapZeros || !Zero.Equals(f(Zero)); @@ -1212,62 +1318,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - public override void MapInplace(Func f, bool forceMapZeros = false) - { - if (forceMapZeros || !Zero.Equals(f(Zero))) - { - var newRowPointers = RowPointers; - var newColumnIndices = new List(ColumnIndices.Length); - var newValues = new List(Values.Length); - - int k = 0; - for (int row = 0; row < RowCount; row++) - { - newRowPointers[row] = newValues.Count; - for (int col = 0; col < ColumnCount; col++) - { - var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(Values[k++]) : f(Zero); - if (!Zero.Equals(item)) - { - newValues.Add(item); - newColumnIndices.Add(col); - } - } - } - - ColumnIndices = newColumnIndices.ToArray(); - Values = newValues.ToArray(); - newRowPointers[RowCount] = newValues.Count; - } - else - { - // we can safely do this in-place: - int nonZero = 0; - for (int row = 0; row < RowCount; row++) - { - var startIndex = RowPointers[row]; - var endIndex = RowPointers[row + 1]; - RowPointers[row] = nonZero; - for (var j = startIndex; j < endIndex; j++) - { - var item = f(Values[j]); - if (!Zero.Equals(item)) - { - Values[nonZero] = item; - ColumnIndices[nonZero] = ColumnIndices[j]; - nonZero++; - } - } - } - Array.Resize(ref ColumnIndices, nonZero); - Array.Resize(ref Values, nonZero); - RowPointers[RowCount] = nonZero; - } - } - internal override void MapIndexedToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - var processZeros = forceMapZeros || !Zero.Equals(f(0, 0, Zero)); + var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero)); var sparseTarget = target as SparseCompressedRowMatrixStorage; if (sparseTarget != null) @@ -1359,56 +1412,213 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } - public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + internal override void MapSubMatrixIndexedToUnchecked(MatrixStorage target, Func f, + int sourceRowIndex, int targetRowIndex, int rowCount, + int sourceColumnIndex, int targetColumnIndex, int columnCount, + bool forceMapZeros = false, bool skipClearing = false) { - if (forceMapZeros || !Zero.Equals(f(0, 0, Zero))) + var sparseTarget = target as SparseCompressedRowMatrixStorage; + if (sparseTarget != null) { - var newRowPointers = RowPointers; - var newColumnIndices = new List(ColumnIndices.Length); - var newValues = new List(Values.Length); + MapSubMatrixIndexedToUnchecked(sparseTarget, f, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, forceMapZeros, skipClearing); + return; + } - int k = 0; - for (int row = 0; row < RowCount; row++) + // FALL BACK + + var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero)); + if (!skipClearing && !processZeros) + { + target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount); + } + + if (processZeros) + { + for (int sr = sourceRowIndex, tr = targetRowIndex; sr < sourceRowIndex + rowCount; sr++, tr++) { - newRowPointers[row] = newValues.Count; - for (int col = 0; col < ColumnCount; col++) + var index = RowPointers[sr]; + var endIndex = RowPointers[sr + 1]; + + // move forward to our sub-range + for (; ColumnIndices[index] < sourceColumnIndex && index < endIndex; index++) { - var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(row, col, Values[k++]) : f(row, col, Zero); - if (!Zero.Equals(item)) + } + for (int sc = sourceColumnIndex, tc = targetColumnIndex; sc < sourceColumnIndex + columnCount; sc++, tc++) + { + if (index < endIndex && sc == ColumnIndices[index]) { - newValues.Add(item); - newColumnIndices.Add(col); + target.At(tr, tc, f(tr, tc, Values[index])); + index = Math.Min(index + 1, endIndex); + } + else + { + target.At(tr, tc, f(tr, tc, Zero)); } } } + } + else + { + int columnOffset = targetColumnIndex - sourceColumnIndex; + for (int sr = sourceRowIndex, tr = targetRowIndex; sr < sourceRowIndex + rowCount; sr++, tr++) + { + var startIndex = RowPointers[sr]; + var endIndex = RowPointers[sr + 1]; + for (int k = startIndex; k < endIndex; k++) + { + // check if the column index is in the range + if ((ColumnIndices[k] >= sourceColumnIndex) && (ColumnIndices[k] < sourceColumnIndex + columnCount)) + { + int tc = ColumnIndices[k] + columnOffset; + target.At(tr, tc, f(tr, tc, Values[k])); + } + } + } + } + } - ColumnIndices = newColumnIndices.ToArray(); - Values = newValues.ToArray(); - newRowPointers[RowCount] = newValues.Count; + void MapSubMatrixIndexedToUnchecked(SparseCompressedRowMatrixStorage target, Func f, + int sourceRowIndex, int targetRowIndex, int rowCount, + int sourceColumnIndex, int targetColumnIndex, int columnCount, + bool forceMapZeros = false, bool skipClearing = false) + where TU : struct, IEquatable, IFormattable + { + var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero)); + if (!skipClearing && !processZeros) + { + target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount); + } + + var rowOffset = targetRowIndex - sourceRowIndex; + var columnOffset = targetColumnIndex - sourceColumnIndex; + var zero = Matrix.Zero; + + // special case for empty target - much faster + if (target.ValueCount == 0) + { + var values = new List(ValueCount); + var columnIndices = new List(ValueCount); + var rowPointers = target.RowPointers; + + if (processZeros) + { + for (int sr = sourceRowIndex; sr < sourceRowIndex + rowCount; sr++) + { + int tr = sr + rowOffset; + rowPointers[tr] = values.Count; + + var index = RowPointers[sr]; + var endIndex = RowPointers[sr + 1]; + + // move forward to our sub-range + for (; ColumnIndices[index] < sourceColumnIndex && index < endIndex; index++) + { + } + for (int sc = sourceColumnIndex, tc = targetColumnIndex; sc < sourceColumnIndex + columnCount; sc++, tc++) + { + if (index < endIndex && sc == ColumnIndices[index]) + { + TU item = f(tr, tc, Values[index]); + if (!zero.Equals(item)) + { + values.Add(item); + columnIndices.Add(tc); + } + index = Math.Min(index + 1, endIndex); + } + else + { + TU item = f(tr, tc, Zero); + if (!zero.Equals(item)) + { + values.Add(item); + columnIndices.Add(tc); + } + } + } + } + } + else + { + for (int sr = sourceRowIndex; sr < sourceRowIndex + rowCount; sr++) + { + int tr = sr + rowOffset; + rowPointers[tr] = values.Count; + + var startIndex = RowPointers[sr]; + var endIndex = RowPointers[sr + 1]; + + for (int k = startIndex; k < endIndex; k++) + { + // check if the column index is in the range + if ((ColumnIndices[k] >= sourceColumnIndex) && (ColumnIndices[k] < sourceColumnIndex + columnCount)) + { + int tc = ColumnIndices[k] + columnOffset; + TU item = f(tr, tc, Values[k]); + if (!zero.Equals(item)) + { + values.Add(item); + columnIndices.Add(tc); + } + } + } + } + } + + for (int i = targetRowIndex + rowCount; i < rowPointers.Length; i++) + { + rowPointers[i] = values.Count; + } + + target.RowPointers[target.RowCount] = values.Count; + target.Values = values.ToArray(); + target.ColumnIndices = columnIndices.ToArray(); + return; + } + + // TODO: proper general sparse case - the following is essentially a fall back, not leveraging the target data structure + + if (processZeros) + { + for (int sr = sourceRowIndex, tr = targetRowIndex; sr < sourceRowIndex + rowCount; sr++, tr++) + { + var index = RowPointers[sr]; + var endIndex = RowPointers[sr + 1]; + + // move forward to our sub-range + for (; ColumnIndices[index] < sourceColumnIndex && index < endIndex; index++) + { + } + for (int sc = sourceColumnIndex, tc = targetColumnIndex; sc < sourceColumnIndex + columnCount; sc++, tc++) + { + if (index < endIndex && sc == ColumnIndices[index]) + { + target.At(tr, tc, f(tr, tc, Values[index])); + index = Math.Min(index + 1, endIndex); + } + else + { + target.At(tr, tc, f(tr, tc, Zero)); + } + } + } } else { - // we can safely do this in-place: - int nonZero = 0; - for (int row = 0; row < RowCount; row++) + for (int sr = sourceRowIndex, tr = targetRowIndex; sr < sourceRowIndex + rowCount; sr++, tr++) { - var startIndex = RowPointers[row]; - var endIndex = RowPointers[row + 1]; - RowPointers[row] = nonZero; - for (var j = startIndex; j < endIndex; j++) + var startIndex = RowPointers[sr]; + var endIndex = RowPointers[sr + 1]; + for (int k = startIndex; k < endIndex; k++) { - var item = f(row, ColumnIndices[j], Values[j]); - if (!Zero.Equals(item)) + // check if the column index is in the range + if ((ColumnIndices[k] >= sourceColumnIndex) && (ColumnIndices[k] < sourceColumnIndex + columnCount)) { - Values[nonZero] = item; - ColumnIndices[nonZero] = ColumnIndices[j]; - nonZero++; + int tc = ColumnIndices[k] + columnOffset; + target.At(tr, tc, f(tr, tc, Values[k])); } } } - Array.Resize(ref ColumnIndices, nonZero); - Array.Resize(ref Values, nonZero); - RowPointers[RowCount] = nonZero; } } } diff --git a/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Map.cs b/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Map.cs new file mode 100644 index 00000000..89dbaac1 --- /dev/null +++ b/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Map.cs @@ -0,0 +1,260 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2014 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System.Linq; +using MathNet.Numerics.LinearAlgebra; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.LinearAlgebraTests +{ + partial class MatrixStructureTheory + { + [Theory] + public void CanMap(Matrix matrix) + { + Matrix a = matrix.Map(x => x, false); + Assert.That(a, Is.EqualTo(matrix)); + Assert.That(a.Storage.IsDense, Is.EqualTo(matrix.Storage.IsDense)); + Assert.That(a.Storage.IsFullyMutable, Is.EqualTo(matrix.Storage.IsFullyMutable)); + + T one = Matrix.Build.One; + Assert.That(matrix.Map(x => x, true), Is.EqualTo(matrix)); + Assert.That(matrix.Map(x => one, true), Is.EqualTo(Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, one))); + + // Map into existing - we skip zeros, but existing values must still be reset to zero + var dense = Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, one); + matrix.Map(x => x, dense, false); + Assert.That(dense, Is.EqualTo(matrix)); + + // Map into self, without using the proper MapInplace method: + var copy = matrix.Clone(); + copy.Map(x => x, copy, false); + Assert.That(copy, Is.EqualTo(matrix)); + if (matrix.Storage.IsFullyMutable) + { + copy.Map(x => one, copy, false); + Assert.That(copy, Is.EqualTo(Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, one))); + } + } + + [Theory] + public void CanMapIndexed(Matrix matrix) + { + Matrix a = matrix.MapIndexed((i, j, x) => + { + if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x)); + return x; + }, false); + Assert.That(a, Is.EqualTo(matrix)); + Assert.That(a.Storage.IsDense, Is.EqualTo(matrix.Storage.IsDense)); + Assert.That(a.Storage.IsFullyMutable, Is.EqualTo(matrix.Storage.IsFullyMutable)); + + T one = Matrix.Build.One; + var d = matrix.MapIndexed((i, j, x) => i == j ? one : x, false); + Assert.That(d.Diagonal().All(x => one.Equals(x) || Zero.Equals(x)), Is.True); + Assert.That(d.EnumerateIndexed().All(z => (z.Item1 == z.Item2) || (matrix.At(z.Item1, z.Item2).Equals(z.Item3)))); + + Assert.That(matrix.MapIndexed((i, j, x) => x, true), Is.EqualTo(matrix)); + Assert.That(matrix.MapIndexed((i, j, x) => one, true), Is.EqualTo(Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, one))); + + // Map into existing - we skip zeros, but existing values must still be reset to zero + var dense = Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, one); + matrix.MapIndexed((i, j, x) => x, dense, false); + Assert.That(dense, Is.EqualTo(matrix)); + + // Map into self, without using the proper MapInplace method: + var copy = matrix.Clone(); + copy.MapIndexed((i, j, x) => x, copy, false); + Assert.That(copy, Is.EqualTo(matrix)); + if (matrix.Storage.IsFullyMutable) + { + copy.MapIndexed((i, j, x) => one, copy, false); + Assert.That(copy, Is.EqualTo(Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, one))); + } + } + + [Theory] + public void CanMapInplace(Matrix matrix) + { + var a = matrix.Clone(); + a.MapInplace(x => x, false); + Assert.That(a, Is.EqualTo(matrix)); + + if (matrix.Storage.IsFullyMutable) + { + a.MapInplace(x => x, true); + Assert.That(a, Is.EqualTo(matrix)); + + T one = Matrix.Build.One; + a.MapInplace(x => one, true); + Assert.That(a, Is.EqualTo(Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, one))); + } + } + + [Theory] + public void CanMapIndexedInplace(Matrix matrix) + { + var a = matrix.Clone(); + a.MapIndexedInplace((i, j, x) => + { + if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x)); + return x; + + }, false); + Assert.That(a, Is.EqualTo(matrix)); + + if (matrix.Storage.IsFullyMutable) + { + a.MapIndexedInplace((i, j, x) => x, true); + Assert.That(a, Is.EqualTo(matrix)); + + T one = Matrix.Build.One; + a.MapIndexedInplace((i, j, x) => i == j ? one : x, false); + Assert.That(a.Diagonal().All(x => one.Equals(x) || Zero.Equals(x)), Is.True); + Assert.That(a.EnumerateIndexed().All(z => (z.Item1 == z.Item2) || (matrix.At(z.Item1, z.Item2).Equals(z.Item3)))); + + a.MapIndexedInplace((i, j, x) => one, true); + Assert.That(a, Is.EqualTo(Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, one))); + } + } + + [Theory] + public void CanMapSubMatrixToSame(Matrix matrix) + { + T one = Matrix.Build.One; + + // Full Range - not forced + Matrix target = Matrix.Build.SameAs(matrix); + matrix.Storage.MapSubMatrixIndexedTo(target.Storage, (i, j, x) => x, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, false, false); + Assert.That(target, Is.EqualTo(matrix), "Full Range - not forced"); + matrix.Storage.MapSubMatrixIndexedTo(target.Storage, (i, j, x) => Zero.Equals(x) ? Zero : one, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, false, false); + Assert.That(target.Enumerate().All(x => Zero.Equals(x) || one.Equals(x)), Is.True); + + if (matrix.Storage.IsFullyMutable) + { + // Full Range - forced + target = Matrix.Build.SameAs(matrix); + matrix.Storage.MapSubMatrixIndexedTo(target.Storage, (i, j, x) => x, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, true, false); + Assert.That(target, Is.EqualTo(matrix), "Full Range - forced"); + matrix.Storage.MapSubMatrixIndexedTo(target.Storage, (i, j, x) => Zero.Equals(x) ? Zero : one, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, true, false); + Assert.That(target.Enumerate().All(x => Zero.Equals(x) || one.Equals(x)), Is.True); + } + } + + [Theory] + public void CanMapSubMatrixToDense(Matrix matrix) + { + T one = Matrix.Build.One; + + // Full Range - not forced + Matrix dense = Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, one); + matrix.Storage.MapSubMatrixIndexedTo(dense.Storage, (i, j, x) => + { + if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x)); + return x; + }, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, false, false); + Assert.That(dense, Is.EqualTo(matrix), "Full Range - not forced"); + + // Full Range - forced + dense = Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, one); + matrix.Storage.MapSubMatrixIndexedTo(dense.Storage, (i, j, x) => + { + if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x)); + return x; + }, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, true, false); + Assert.That(dense, Is.EqualTo(matrix), "Full Range - forced"); + + // Sub Range - not forced - all except first column padded into 1-border + dense = Matrix.Build.Dense(matrix.RowCount + 2, matrix.ColumnCount + 1, one); + matrix.Storage.MapSubMatrixIndexedTo(dense.Storage, (i, j, x) => x, 0, 1, matrix.RowCount, 1, 1, matrix.ColumnCount - 1, false, false); + Assert.That(dense.SubMatrix(1, dense.RowCount - 2, 1, dense.ColumnCount - 2), + Is.EqualTo(matrix.SubMatrix(0, matrix.RowCount, 1, matrix.ColumnCount - 1)), "Sub Range - not forced - range"); + dense.SetSubMatrix(1, 0, matrix.RowCount, 1, 0, matrix.ColumnCount - 1, Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount - 1, one)); + Assert.That(dense.Enumerate().All(one.Equals), Is.True); + + // Sub Range - forced - all except first row padded into 1-border + dense = Matrix.Build.Dense(matrix.RowCount + 1, matrix.ColumnCount + 2, one); + matrix.Storage.MapSubMatrixIndexedTo(dense.Storage, (i, j, x) => x, 1, 1, matrix.RowCount - 1, 0, 1, matrix.ColumnCount, true, false); + Assert.That(dense.SubMatrix(1, dense.RowCount - 2, 1, dense.ColumnCount - 2), + Is.EqualTo(matrix.SubMatrix(1, matrix.RowCount - 1, 0, matrix.ColumnCount)), "Sub Range - forced - range"); + dense.SetSubMatrix(1, 0, matrix.RowCount - 1, 1, 0, matrix.ColumnCount, Matrix.Build.Dense(matrix.RowCount - 1, matrix.ColumnCount, one)); + Assert.That(dense.Enumerate().All(one.Equals), Is.True); + } + + [Theory] + public void CanMapSubMatrixToSparse(Matrix matrix) + { + T one = Matrix.Build.One; + + // Full Range - filled, not forced + Matrix sparse = Matrix.Build.Sparse(matrix.RowCount, matrix.ColumnCount, one); + matrix.Storage.MapSubMatrixIndexedTo(sparse.Storage, (i, j, x) => + { + if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x)); + return x; + }, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, false, false); + Assert.That(sparse, Is.EqualTo(matrix), "Full Range - filled, not forced"); + + // Full Range - empty, not forced + sparse = Matrix.Build.Sparse(matrix.RowCount, matrix.ColumnCount); + matrix.Storage.MapSubMatrixIndexedTo(sparse.Storage, (i, j, x) => + { + if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x)); + return x; + }, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, false, false); + Assert.That(sparse, Is.EqualTo(matrix), "Full Range - empty, not forced"); + + // Full Range - filled, forced + sparse = Matrix.Build.Sparse(matrix.RowCount, matrix.ColumnCount, one); + matrix.Storage.MapSubMatrixIndexedTo(sparse.Storage, (i, j, x) => + { + if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x)); + return x; + }, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, true, false); + Assert.That(sparse, Is.EqualTo(matrix), "Full Range - filled, forced"); + + // Sub Range - filled, not forced - all except first column padded into 1-border + sparse = Matrix.Build.Sparse(matrix.RowCount + 2, matrix.ColumnCount + 1, one); + matrix.Storage.MapSubMatrixIndexedTo(sparse.Storage, (i, j, x) => x, 0, 1, matrix.RowCount, 1, 1, matrix.ColumnCount - 1, false, false); + Assert.That(sparse.SubMatrix(1, sparse.RowCount - 2, 1, sparse.ColumnCount - 2), + Is.EqualTo(matrix.SubMatrix(0, matrix.RowCount, 1, matrix.ColumnCount - 1)), "Sub Range - filled, not forced - range"); + sparse.SetSubMatrix(1, 0, matrix.RowCount, 1, 0, matrix.ColumnCount - 1, Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount - 1, one)); + Assert.That(sparse.Enumerate().All(one.Equals), Is.True); + + // Sub Range - filled, forced - all except first row padded into 1-border + sparse = Matrix.Build.Sparse(matrix.RowCount + 1, matrix.ColumnCount + 2, one); + matrix.Storage.MapSubMatrixIndexedTo(sparse.Storage, (i, j, x) => x, 1, 1, matrix.RowCount - 1, 0, 1, matrix.ColumnCount, true, false); + Assert.That(sparse.SubMatrix(1, sparse.RowCount - 2, 1, sparse.ColumnCount - 2), + Is.EqualTo(matrix.SubMatrix(1, matrix.RowCount - 1, 0, matrix.ColumnCount)), "Sub Range - filled, forced - range"); + sparse.SetSubMatrix(1, 0, matrix.RowCount - 1, 1, 0, matrix.ColumnCount, Matrix.Build.Dense(matrix.RowCount - 1, matrix.ColumnCount, one)); + Assert.That(sparse.Enumerate().All(one.Equals), Is.True); + } + } +} diff --git a/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Reform.cs b/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Reform.cs index dd123a96..7efb3095 100644 --- a/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Reform.cs +++ b/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Reform.cs @@ -28,10 +28,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.LinearAlgebra; -using NUnit.Framework; using System; using System.Linq; +using MathNet.Numerics.LinearAlgebra; +using NUnit.Framework; namespace MathNet.Numerics.UnitTests.LinearAlgebraTests { diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 5ab21410..fc3e9f17 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -305,6 +305,7 @@ +