diff --git a/src/FSharp/LinearAlgebra.Matrix.fs b/src/FSharp/LinearAlgebra.Matrix.fs index 7ad6436c..672263d5 100644 --- a/src/FSharp/LinearAlgebra.Matrix.fs +++ b/src/FSharp/LinearAlgebra.Matrix.fs @@ -207,22 +207,18 @@ module Matrix = /// In-place map of every matrix element using a function. - let inline mapInPlace f (A: #Matrix<_>) = - A.MapInplace((fun x -> f x), true) - - /// In-place map of every matrix element using a position dependent function. - let inline mapiInPlace f (A: #Matrix<_>) = - A.MapIndexedInplace((fun i j x -> f i j x), true) + let inline mapInPlace f (A: #Matrix<_>) = A.MapInplace((fun x -> f x), true) /// In-place map of every matrix element using a function. /// Zero-values may be skipped (relevant mostly for sparse matrices). - let inline mapSkipZerosInPlace f (A: #Matrix<_>) = - A.MapInplace((fun x -> f x), false) + let inline mapSkipZerosInPlace f (A: #Matrix<_>) = A.MapInplace((fun x -> f x), false) + + /// In-place map of every matrix element using a position dependent function. + let inline mapiInPlace f (A: #Matrix<_>) = A.MapIndexedInplace((fun i j x -> f i j x), true) /// In-place map of every matrix element using a position dependent function. /// Zero-values may be skipped (relevant mostly for sparse matrices). - let inline mapiSkipZerosInPlace f (A: #Matrix<_>) = - A.MapIndexedInplace((fun i j x -> f i j x), false) + let inline mapiSkipZerosInPlace f (A: #Matrix<_>) = A.MapIndexedInplace((fun i j x -> f i j x), false) /// In-place map every matrix column using the given position dependent function. let inline mapColsInPlace (f: int -> Vector<'a> -> Vector<'a>) (A: #Matrix<_>) = @@ -236,30 +232,18 @@ module Matrix = /// Map every matrix element using the given function. - let inline map f (A: #Matrix<_>) = - let A = A.Clone() - A.MapInplace((fun x -> f x), true) - A + let inline map f (A: #Matrix<_>) = A.Map((fun x -> f x), true) /// Map every matrix element using the given function. /// Zero-values may be skipped (relevant mostly for sparse matrices). - let inline mapSkipZeros f (A: #Matrix<_>) = - let A = A.Clone() - A.MapInplace((fun x -> f x), false) - A + let inline mapSkipZeros f (A: #Matrix<_>) = A.Map((fun x -> f x), false) /// Map every matrix element using the given position dependent function. - let inline mapi f (A: #Matrix<_>) = - let A = A.Clone() - A.MapIndexedInplace((fun i j x -> f i j x), true) - A + let inline mapi f (A: #Matrix<_>) = A.MapIndexed((fun i j x -> f i j x), true) /// Map every matrix element using the given position dependent function. /// Zero-values may be skipped (relevant mostly for sparse matrices). - let inline mapiSkipZeros f (A: #Matrix<_>) = - let A = A.Clone() - A.MapIndexedInplace((fun i j x -> f i j x), false) - A + let inline mapiSkipZeros f (A: #Matrix<_>) = A.MapIndexed((fun i j x -> f i j x), false) /// Map every matrix column using the given position dependent function. let inline mapCols (f: int -> Vector<'a> -> Vector<'a>) (A: #Matrix<_>) = diff --git a/src/Numerics/LinearAlgebra/Matrix.cs b/src/Numerics/LinearAlgebra/Matrix.cs index b349f593..d420c968 100644 --- a/src/Numerics/LinearAlgebra/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Matrix.cs @@ -1455,7 +1455,7 @@ namespace MathNet.Numerics.LinearAlgebra /// public void MapInplace(Func f, bool forceMapZeros = false) { - Storage.MapInplace(f, forceMapZeros); + Storage.MapToUnchecked(Storage, f, forceMapZeros, skipClearing: true); } /// @@ -1466,7 +1466,57 @@ namespace MathNet.Numerics.LinearAlgebra /// public void MapIndexedInplace(Func f, bool forceMapZeros = false) { - Storage.MapIndexedInplace(f, forceMapZeros); + Storage.MapIndexedToUnchecked(Storage, f, forceMapZeros, skipClearing: true); + } + + /// + /// Applies a function to each value of this matrix and replaces the value in the result matrix. + /// If forceMapZero is not set to true, zero values may or may not be skipped depending + /// on the actual data storage implementation (relevant mostly for sparse matrices). + /// + public void Map(Func f, Matrix result, bool forceMapZeros = false) + where TU : struct, IEquatable, IFormattable + { + Storage.MapTo(result.Storage, f, forceMapZeros, skipClearing: forceMapZeros); + } + + /// + /// Applies a function to each value of this matrix and replaces the value in the result matrix. + /// The index of each value (zero-based) is passed as first argument to the function. + /// If forceMapZero is not set to true, zero values may or may not be skipped depending + /// on the actual data storage implementation (relevant mostly for sparse matrices). + /// + public void MapIndexed(Func f, Matrix result, bool forceMapZeros = false) + where TU : struct, IEquatable, IFormattable + { + Storage.MapIndexedTo(result.Storage, f, forceMapZeros, skipClearing: forceMapZeros); + } + + /// + /// Applies a function to each value of this matrix and returns the results as a new matrix. + /// If forceMapZero is not set to true, zero values may or may not be skipped depending + /// on the actual data storage implementation (relevant mostly for sparse matrices). + /// + public Matrix Map(Func f, bool forceMapZeros = false) + where TU : struct, IEquatable, IFormattable + { + var result = Matrix.Build.SameAs(this); + Storage.MapToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true); + return result; + } + + /// + /// Applies a function to each value of this matrix and returns the results as a new matrix. + /// The index of each value (zero-based) is passed as first argument to the function. + /// If forceMapZero is not set to true, zero values may or may not be skipped depending + /// on the actual data storage implementation (relevant mostly for sparse matrices). + /// + public Matrix MapIndexed(Func f, bool forceMapZeros = false) + where TU : struct, IEquatable, IFormattable + { + var result = Matrix.Build.SameAs(this); + Storage.MapIndexedToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true); + return result; } } } diff --git a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs index 77c21efe..35947b20 100644 --- a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs @@ -503,26 +503,58 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS - public override void MapInplace(Func f, bool forceMapZeros = false) + internal override void MapToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - CommonParallel.For(0, Data.Length, 4096, (a, b) => + var denseTarget = target as DenseColumnMajorMatrixStorage; + if (denseTarget != null) + { + CommonParallel.For(0, Data.Length, 4096, (a, b) => { for (int i = a; i < b; i++) { - Data[i] = f(Data[i]); + denseTarget.Data[i] = f(Data[i]); } }); + return; + } + + // FALL BACK + + int index = 0; + for (int j = 0; j < ColumnCount; j++) + { + for (int i = 0; i < RowCount; i++) + { + target.At(i, j, f(Data[index++])); + } + } } - public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + internal override void MapIndexedToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - int index = 0; + var denseTarget = target as DenseColumnMajorMatrixStorage; + if (denseTarget != null) + { + int index = 0; + for (int j = 0; j < ColumnCount; j++) + { + for (int i = 0; i < RowCount; i++) + { + denseTarget.Data[index] = f(i, j, Data[index]); + index++; + } + } + return; + } + + // FALL BACK + + int index2 = 0; for (int j = 0; j < ColumnCount; j++) { for (int i = 0; i < RowCount; i++) { - Data[index] = f(i, j, Data[index]); - index++; + target.At(i, j, f(i, j, Data[index2++])); } } } diff --git a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs index f986ff2c..b863ec8e 100644 --- a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs @@ -589,30 +589,100 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS - public override void MapInplace(Func f, bool forceMapZeros = false) + internal override void MapToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - // we deliberately ignore forceMapZeros since we would not actually - // support any non-zero results outside of the diagonal anyway - CommonParallel.For(0, Data.Length, 4096, (a, b) => + var processZeros = forceMapZeros || !Zero.Equals(f(Zero)); + + var diagonalTarget = target as DiagonalMatrixStorage; + if (diagonalTarget != null) + { + if (processZeros) + { + 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]); + diagonalTarget.Data[i] = f(Data[i]); } }); + return; + } + + // FALL BACK + + if (!skipClearing && !processZeros) + { + target.Clear(); + } + + if (processZeros) + { + for (int j = 0; j < ColumnCount; j++) + { + for (int i = 0; i < RowCount; i++) + { + target.At(i, j, f(i == j ? Data[i] : Zero)); + } + } + } + else + { + for (int i = 0; i < Data.Length; i++) + { + target.At(i, i, f(Data[i])); + } + } } - public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + internal override void MapIndexedToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - // we deliberately ignore forceMapZeros since we would not actually - // support any non-zero results outside of the diagonal anyway - CommonParallel.For(0, Data.Length, 4096, (a, b) => + var processZeros = forceMapZeros || !Zero.Equals(f(0, 0, Zero)); + + var diagonalTarget = target as DiagonalMatrixStorage; + if (diagonalTarget != null) + { + if (processZeros) + { + 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]); + diagonalTarget.Data[i] = f(i, i, Data[i]); } }); + return; + } + + // FALL BACK + + if (!skipClearing && !processZeros) + { + target.Clear(); + } + + if (processZeros) + { + for (int j = 0; j < ColumnCount; j++) + { + for (int i = 0; i < RowCount; i++) + { + target.At(i, j, f(i, j, i == j ? Data[i] : Zero)); + } + } + } + else + { + for (int i = 0; i < Data.Length; i++) + { + target.At(i, i, f(i, i, Data[i])); + } + } } } } diff --git a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs index d3d4fab9..45cdeb85 100644 --- a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs @@ -468,24 +468,60 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS - public virtual void MapInplace(Func f, bool forceMapZeros = false) + public void MapTo(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) + where TU : struct, IEquatable, IFormattable + { + if (target == null) + { + throw new ArgumentNullException("target"); + } + + if (RowCount != target.RowCount || ColumnCount != target.ColumnCount) + { + var message = string.Format(Resources.ArgumentMatrixDimensions2, RowCount + "x" + ColumnCount, target.RowCount + "x" + target.ColumnCount); + throw new ArgumentException(message, "target"); + } + + MapToUnchecked(target, f, forceMapZeros, skipClearing); + } + + internal virtual void MapToUnchecked(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++) { - At(i, j, f(At(i, j))); + target.At(i, j, f(At(i, j))); } } } - public virtual void MapIndexedInplace(Func f, bool forceMapZeros = false) + public void MapIndexedTo(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) + where TU : struct, IEquatable, IFormattable + { + if (target == null) + { + throw new ArgumentNullException("target"); + } + + if (RowCount != target.RowCount || ColumnCount != target.ColumnCount) + { + var message = string.Format(Resources.ArgumentMatrixDimensions2, RowCount + "x" + ColumnCount, target.RowCount + "x" + target.ColumnCount); + throw new ArgumentException(message, "target"); + } + + MapIndexedToUnchecked(target, f, forceMapZeros, skipClearing); + } + + 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++) { - At(i, j, f(i, j, At(i, j))); + target.At(i, j, f(i, j, At(i, j))); } } } diff --git a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs index 6d45c9e7..3d8c7e5a 100644 --- a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs @@ -30,6 +30,7 @@ using System; using System.Collections.Generic; +using System.Diagnostics; using System.Linq; using MathNet.Numerics.Properties; @@ -42,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // [ruegg] public fields are OK here /// - /// The array containing the row indices of the existing rows. Element "i" of the array gives the index of the + /// The array containing the row indices of the existing rows. Element "i" of the array gives the index of the /// element in the array that is first non-zero element in a row "i". /// The last value is equal to ValueCount, so that the number of non-zero entries in row "i" is always /// given by RowPointers[i+i] - RowPointers[i]. This array thus has length RowCount+1. @@ -50,13 +51,13 @@ namespace MathNet.Numerics.LinearAlgebra.Storage public readonly int[] RowPointers; /// - /// An array containing the column indices of the non-zero values. Element "j" of the array + /// An array containing the column indices of the non-zero values. Element "j" of the array /// is the number of the column in matrix that contains the j-th value in the array. /// public int[] ColumnIndices; /// - /// Array that contains the non-zero elements of matrix. Values of the non-zero elements of matrix are mapped into the values + /// Array that contains the non-zero elements of matrix. Values of the non-zero elements of matrix are mapped into the values /// array using the row-major storage mapping described in a compressed sparse row (CSR) format. /// public T[] Values; @@ -214,7 +215,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage valueCount -= 1; - // Check whether we need to shrink the arrays. This is reasonable to do if + // Check whether we need to shrink the arrays. This is reasonable to do if // there are a lot of non-zero elements and storage is two times bigger if ((valueCount > 1024) && (valueCount < Values.Length / 2)) { @@ -374,7 +375,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage /// Returns a hash code for this instance. /// /// - /// A hash code for this instance, suitable for use in hashing algorithms and data structures like a hash table. + /// A hash code for this instance, suitable for use in hashing algorithms and data structures like a hash table. /// public override int GetHashCode() { @@ -1104,25 +1105,82 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS - public override void MapInplace(Func f, bool forceMapZeros = false) + internal override void MapToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - var newRowPointers = new int[RowCount+1]; - var newColumnIndices = new List(); - var newValues = new List(); + var processZeros = forceMapZeros || !Zero.Equals(f(Zero)); + + var sparseTarget = target as SparseCompressedRowMatrixStorage; + if (sparseTarget != null) + { + var newRowPointers = sparseTarget.RowPointers; + var newColumnIndices = new List(); + var newValues = new List(); + + if (processZeros) + { + 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); + } + } + } + } + else + { + for (int row = 0; row < RowCount; row++) + { + newRowPointers[row] = newValues.Count; + var startIndex = RowPointers[row]; + var endIndex = RowPointers[row + 1]; + for (var j = startIndex; j < endIndex; j++) + { + var item = f(Values[j]); + if (!Zero.Equals(item)) + { + newValues.Add(item); + newColumnIndices.Add(ColumnIndices[j]); + } + } + } + } + + sparseTarget.ColumnIndices = newColumnIndices.ToArray(); + sparseTarget.Values = newValues.ToArray(); + newRowPointers[RowCount] = newValues.Count; + return; + } + + // FALL BACK - if (forceMapZeros || !Zero.Equals(f(Zero))) + if (!skipClearing && !processZeros) + { + target.Clear(); + } + + if (processZeros) { - int k = 0; for (int row = 0; row < RowCount; row++) { - newRowPointers[row] = newValues.Count; - for (int col = 0; col < ColumnCount; col++) + var index = RowPointers[row]; + var endIndex = RowPointers[row + 1]; + for (int j = 0; j < ColumnCount; j++) { - var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(Values[k++]) : f(Zero); - if (!Zero.Equals(item)) + if (j == ColumnIndices[index]) { - newValues.Add(item); - newColumnIndices.Add(col); + target.At(row, j, f(Values[j])); + index = Math.Min(index + 1, endIndex); + } + else + { + target.At(row, j, f(Zero)); } } } @@ -1131,46 +1189,92 @@ namespace MathNet.Numerics.LinearAlgebra.Storage { for (int row = 0; row < RowCount; row++) { - newRowPointers[row] = newValues.Count; var startIndex = RowPointers[row]; var endIndex = RowPointers[row + 1]; for (var j = startIndex; j < endIndex; j++) { - var item = f(Values[j]); - if (!Zero.Equals(item)) + target.At(row, ColumnIndices[j], f(Values[j])); + } + } + } + } + + internal override void MapIndexedToUnchecked(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) + { + var processZeros = forceMapZeros || !Zero.Equals(f(0, 0, Zero)); + + var sparseTarget = target as SparseCompressedRowMatrixStorage; + if (sparseTarget != null) + { + var newRowPointers = sparseTarget.RowPointers; + var newColumnIndices = new List(); + var newValues = new List(); + + if (processZeros) + { + int k = 0; + for (int row = 0; row < RowCount; row++) + { + newRowPointers[row] = newValues.Count; + for (int col = 0; col < ColumnCount; col++) { - newValues.Add(item); - newColumnIndices.Add(ColumnIndices[j]); + 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); + } } } } + else + { + for (int row = 0; row < RowCount; row++) + { + newRowPointers[row] = newValues.Count; + var startIndex = RowPointers[row]; + var endIndex = RowPointers[row + 1]; + for (var j = startIndex; j < endIndex; j++) + { + var item = f(row, ColumnIndices[j], Values[j]); + if (!Zero.Equals(item)) + { + newValues.Add(item); + newColumnIndices.Add(ColumnIndices[j]); + } + } + } + } + + sparseTarget.ColumnIndices = newColumnIndices.ToArray(); + sparseTarget.Values = newValues.ToArray(); + newRowPointers[RowCount] = newValues.Count; + return; } - ColumnIndices = newColumnIndices.ToArray(); - Values = newValues.ToArray(); - newRowPointers[RowCount] = newValues.Count; - Array.Copy(newRowPointers, RowPointers, newRowPointers.Length); - } + // FALL BACK - public override void MapIndexedInplace(Func f, bool forceMapZeros = false) - { - var newRowPointers = new int[RowCount+1]; - var newColumnIndices = new List(); - var newValues = new List(); + if (!skipClearing && !processZeros) + { + target.Clear(); + } - if (forceMapZeros || !Zero.Equals(f(0, 0, Zero))) + if (processZeros) { - int k = 0; for (int row = 0; row < RowCount; row++) { - newRowPointers[row] = newValues.Count; - for (int col = 0; col < ColumnCount; col++) + var index = RowPointers[row]; + var endIndex = RowPointers[row + 1]; + for (int j = 0; j < ColumnCount; j++) { - var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(row, col, Values[k++]) : f(row, col, Zero); - if (!Zero.Equals(item)) + if (j == ColumnIndices[index]) + { + target.At(row, j, f(row, j, Values[j])); + index = Math.Min(index + 1, endIndex); + } + else { - newValues.Add(item); - newColumnIndices.Add(col); + target.At(row, j, f(row, j, Zero)); } } } @@ -1179,25 +1283,14 @@ namespace MathNet.Numerics.LinearAlgebra.Storage { for (int row = 0; row < RowCount; row++) { - newRowPointers[row] = newValues.Count; var startIndex = RowPointers[row]; var endIndex = RowPointers[row + 1]; for (var j = startIndex; j < endIndex; j++) { - var item = f(row, ColumnIndices[j], Values[j]); - if (!Zero.Equals(item)) - { - newValues.Add(item); - newColumnIndices.Add(ColumnIndices[j]); - } + target.At(row, ColumnIndices[j], f(row, ColumnIndices[j], Values[j])); } } } - - ColumnIndices = newColumnIndices.ToArray(); - Values = newValues.ToArray(); - newRowPointers[RowCount] = newValues.Count; - Array.Copy(newRowPointers, RowPointers, newRowPointers.Length); } } }