From f73bc67ee06ab83fbf387a48825477cdcaaa0537 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Thu, 3 Apr 2014 14:24:17 +0200 Subject: [PATCH] LA: storage-aware non-inplace map on vectors --- src/FSharp/LinearAlgebra.Vector.fs | 32 +--- src/Numerics/LinearAlgebra/Builder.cs | 12 +- .../Storage/DenseVectorStorage.cs | 50 ++++-- .../Storage/SparseVectorStorage.cs | 157 +++++++++++++----- .../LinearAlgebra/Storage/VectorStorage.cs | 42 ++++- src/Numerics/LinearAlgebra/Vector.cs | 54 +++++- 6 files changed, 262 insertions(+), 85 deletions(-) diff --git a/src/FSharp/LinearAlgebra.Vector.fs b/src/FSharp/LinearAlgebra.Vector.fs index e013ce26..df2187c1 100644 --- a/src/FSharp/LinearAlgebra.Vector.fs +++ b/src/FSharp/LinearAlgebra.Vector.fs @@ -129,49 +129,33 @@ module Vector = /// In-place mutation by applying a function to every element of the vector. - let inline mapInPlace f (v: #Vector<_>) = - v.MapInplace((fun x -> f x), true) + let inline mapInPlace f (v: #Vector<_>) = v.MapInplace((fun x -> f x), true) /// In-place mutation by applying a function to every element of the vector. - let inline mapiInPlace f (v: #Vector<_>) = - v.MapIndexedInplace((fun i x -> f i x), true) + let inline mapiInPlace f (v: #Vector<_>) = v.MapIndexedInplace((fun i x -> f i x), true) /// In-place mutation by applying a function to every element of the vector. /// Zero-values may be skipped (relevant mostly for sparse vectors). - let inline mapSkipZerosInPlace f (v: #Vector<_>) = - v.MapInplace((fun x -> f x), false) + let inline mapSkipZerosInPlace f (v: #Vector<_>) = v.MapInplace((fun x -> f x), false) /// In-place mutation by applying a function to every element of the vector. /// Zero-values may be skipped (relevant mostly for sparse vectors). - let inline mapiSkipZerosInPlace (f: int -> float -> float) (v: #Vector) = - v.MapIndexedInplace((fun i x -> f i x), false) + let inline mapiSkipZerosInPlace f (v: #Vector<_>) = v.MapIndexedInplace((fun i x -> f i x), false) /// Maps a vector to a new vector by applying a function to every element. - let inline map f (v: #Vector<_>) = - let w = v.Clone() - w.MapInplace((fun x -> f x), true) - w + let inline map f (v: #Vector<_>) = v.Map((fun x -> f x), true) /// Maps a vector to a new vector by applying a function to every element. /// Zero-values may be skipped (relevant mostly for sparse vectors). - let inline mapSkipZeros f (v: #Vector<_>) = - let w = v.Clone() - w.MapInplace((fun x -> f x), false) - w + let inline mapSkipZeros f (v: #Vector<_>) = v.Map((fun x -> f x), false) /// Maps a vector to a new vector by applying a function to every element. - let inline mapi f (v: #Vector<_>) = - let w = v.Clone() - w.MapIndexedInplace((fun i x -> f i x), true) - w + let inline mapi f (v: #Vector<_>) = v.MapIndexed((fun i x -> f i x), true) /// Maps a vector to a new vector by applying a function to every element. /// Zero-values may be skipped (relevant mostly for sparse vectors). - let inline mapiSkipZeros f (v: #Vector<_>) = - let w = v.Clone() - w.MapIndexedInplace((fun i x -> f i x), false) - w + let inline mapiSkipZeros f (v: #Vector<_>) = v.MapIndexed((fun i x -> f i x), false) diff --git a/src/Numerics/LinearAlgebra/Builder.cs b/src/Numerics/LinearAlgebra/Builder.cs index 9d5251bf..5201d4a6 100644 --- a/src/Numerics/LinearAlgebra/Builder.cs +++ b/src/Numerics/LinearAlgebra/Builder.cs @@ -389,7 +389,8 @@ namespace MathNet.Numerics.LinearAlgebra /// /// Create a new matrix with the same kind of the provided example. /// - public Matrix SameAs(Matrix example, int rows, int columns, bool fullyMutable = false) + public Matrix SameAs(Matrix example, int rows, int columns, bool fullyMutable = false) + where TU : struct, IEquatable, IFormattable { var storage = example.Storage; if (storage is DenseColumnMajorMatrixStorage) return Dense(rows, columns); @@ -401,7 +402,8 @@ namespace MathNet.Numerics.LinearAlgebra /// /// Create a new matrix with the same kind and dimensions of the provided example. /// - public Matrix SameAs(Matrix example) + public Matrix SameAs(Matrix example) + where TU : struct, IEquatable, IFormattable { return SameAs(example, example.RowCount, example.ColumnCount); } @@ -1338,7 +1340,8 @@ namespace MathNet.Numerics.LinearAlgebra /// /// Create a new vector with the same kind of the provided example. /// - public Vector SameAs(Vector example, int length) + public Vector SameAs(Vector example, int length) + where TU : struct, IEquatable, IFormattable { return example.Storage.IsDense ? Dense(length) : Sparse(length); } @@ -1346,7 +1349,8 @@ namespace MathNet.Numerics.LinearAlgebra /// /// Create a new vector with the same kind and dimension of the provided example. /// - public Vector SameAs(Vector example) + public Vector SameAs(Vector example) + where TU : struct, IEquatable, IFormattable { return example.Storage.IsDense ? Dense(example.Count) : Sparse(example.Count); } diff --git a/src/Numerics/LinearAlgebra/Storage/DenseVectorStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseVectorStorage.cs index e79486c0..e99c3225 100644 --- a/src/Numerics/LinearAlgebra/Storage/DenseVectorStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/DenseVectorStorage.cs @@ -118,12 +118,12 @@ namespace MathNet.Numerics.LinearAlgebra.Storage var data = new T[length]; CommonParallel.For(0, data.Length, 4096, (a, b) => + { + for (int i = a; i < b; i++) { - for (int i = a; i < b; i++) - { - data[i] = init(i); - } - }); + data[i] = init(i); + } + }); return new DenseVectorStorage(length, data); } @@ -212,7 +212,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage var denseTarget = target as DenseColumnMajorMatrixStorage; if (denseTarget != null) { - Array.Copy(Data, 0, denseTarget.Data, columnIndex * denseTarget.RowCount, Data.Length); + Array.Copy(Data, 0, denseTarget.Data, columnIndex*denseTarget.RowCount, Data.Length); return; } @@ -253,7 +253,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage { for (int j = 0; j < Data.Length; j++) { - denseTarget.Data[(j + targetColumnIndex) * target.RowCount + rowIndex] = Data[j + sourceColumnIndex]; + denseTarget.Data[(j + targetColumnIndex)*target.RowCount + rowIndex] = Data[j + sourceColumnIndex]; } return; } @@ -317,26 +317,50 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS - public override void MapInplace(Func f, bool forceMapZeros = false) + internal override void MapToUnchecked(VectorStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - CommonParallel.For(0, Data.Length, 4096, (a, b) => + var denseTarget = target as DenseVectorStorage; + 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 + + for (int i = 0; i < Length; i++) + { + target.At(i, f(Data[i])); + } } - public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + internal override void MapIndexedToUnchecked(VectorStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - CommonParallel.For(0, Data.Length, 4096, (a, b) => + var denseTarget = target as DenseVectorStorage; + if (denseTarget != null) + { + CommonParallel.For(0, Data.Length, 4096, (a, b) => { for (int i = a; i < b; i++) { - Data[i] = f(i, Data[i]); + denseTarget.Data[i] = f(i, Data[i]); } }); + return; + } + + // FALL BACK + + for (int i = 0; i < Length; i++) + { + target.At(i, f(i, Data[i])); + } } } } diff --git a/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs index 281957fe..9340dee1 100644 --- a/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs @@ -32,6 +32,7 @@ using System; using System.Collections.Generic; using System.Linq; using MathNet.Numerics.Properties; +using MathNet.Numerics.Threading; namespace MathNet.Numerics.LinearAlgebra.Storage { @@ -618,72 +619,152 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS - public override void MapInplace(Func f, bool forceMapZeros = false) + internal override void MapToUnchecked(VectorStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - var indices = new List(); - var values = new List(); - if (forceMapZeros || !Zero.Equals(f(Zero))) + var sparseTarget = target as SparseVectorStorage; + if (sparseTarget != null) { - int k = 0; - for (int i = 0; i < Length; i++) + var indices = new List(); + var values = new List(); + if (forceMapZeros || !Zero.Equals(f(Zero))) + { + int k = 0; + for (int i = 0; i < Length; i++) + { + var item = k < ValueCount && (Indices[k]) == i ? f(Values[k++]) : f(Zero); + if (!Zero.Equals(item)) + { + values.Add(item); + indices.Add(i); + } + } + } + else { - var item = k < ValueCount && (Indices[k]) == i ? f(Values[k++]) : f(Zero); - if (!Zero.Equals(item)) + for (int i = 0; i < ValueCount; i++) { - values.Add(item); - indices.Add(i); + var item = f(Values[i]); + if (!Zero.Equals(item)) + { + values.Add(item); + indices.Add(Indices[i]); + } } } + sparseTarget.Indices = indices.ToArray(); + sparseTarget.Values = values.ToArray(); + sparseTarget.ValueCount = values.Count; + return; } - else + + var denseTarget = target as DenseVectorStorage; + if (denseTarget != null) { - for (int i = 0; i < ValueCount; i++) + if (!skipClearing) + { + denseTarget.Clear(); + } + + if (forceMapZeros || !Zero.Equals(f(Zero))) { - var item = f(Values[i]); - if (!Zero.Equals(item)) + int k = 0; + for (int i = 0; i < Length; i++) { - values.Add(item); - indices.Add(Indices[i]); + denseTarget.Data[i] = k < ValueCount && (Indices[k]) == i + ? f(Values[k++]) + : f(Zero); } } + else + { + CommonParallel.For(0, ValueCount, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + denseTarget.Data[Indices[i]] = f(Values[i]); + } + }); + } + return; } - Indices = indices.ToArray(); - Values = values.ToArray(); - ValueCount = values.Count; + + // FALL BACK + + base.MapToUnchecked(target, f, forceMapZeros, skipClearing); } - public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + internal override void MapIndexedToUnchecked(VectorStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { - var indices = new List(); - var values = new List(); - if (forceMapZeros || !Zero.Equals(f(0, Zero))) + var sparseTarget = target as SparseVectorStorage; + if (sparseTarget != null) { - int k = 0; - for (int i = 0; i < Length; i++) + var indices = new List(); + var values = new List(); + if (forceMapZeros || !Zero.Equals(f(0, Zero))) + { + int k = 0; + for (int i = 0; i < Length; i++) + { + var item = k < ValueCount && (Indices[k]) == i ? f(i, Values[k++]) : f(i, Zero); + if (!Zero.Equals(item)) + { + values.Add(item); + indices.Add(i); + } + } + } + else { - var item = k < ValueCount && (Indices[k]) == i ? f(i, Values[k++]) : f(i, Zero); - if (!Zero.Equals(item)) + for (int i = 0; i < ValueCount; i++) { - values.Add(item); - indices.Add(i); + var item = f(Indices[i], Values[i]); + if (!Zero.Equals(item)) + { + values.Add(item); + indices.Add(Indices[i]); + } } } + sparseTarget.Indices = indices.ToArray(); + sparseTarget.Values = values.ToArray(); + sparseTarget.ValueCount = values.Count; + return; } - else + + var denseTarget = target as DenseVectorStorage; + if (denseTarget != null) { - for (int i = 0; i < ValueCount; i++) + if (!skipClearing) + { + denseTarget.Clear(); + } + + if (forceMapZeros || !Zero.Equals(f(0, Zero))) { - var item = f(Indices[i], Values[i]); - if (!Zero.Equals(item)) + int k = 0; + for (int i = 0; i < Length; i++) { - values.Add(item); - indices.Add(Indices[i]); + denseTarget.Data[i] = k < ValueCount && (Indices[k]) == i + ? f(i, Values[k++]) + : f(i, Zero); } } + else + { + CommonParallel.For(0, ValueCount, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + denseTarget.Data[Indices[i]] = f(Indices[i], Values[i]); + } + }); + } + return; } - Indices = indices.ToArray(); - Values = values.ToArray(); - ValueCount = values.Count; + + // FALL BACK + + base.MapIndexedToUnchecked(target, f, forceMapZeros, skipClearing); } } } diff --git a/src/Numerics/LinearAlgebra/Storage/VectorStorage.cs b/src/Numerics/LinearAlgebra/Storage/VectorStorage.cs index b1bbccb3..ff924a63 100644 --- a/src/Numerics/LinearAlgebra/Storage/VectorStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/VectorStorage.cs @@ -401,19 +401,53 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS - public virtual void MapInplace(Func f, bool forceMapZeros = false) + public void MapTo(VectorStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) + where TU : struct, IEquatable, IFormattable + { + if (target == null) + { + throw new ArgumentNullException("target"); + } + + if (Length != target.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "target"); + } + + MapToUnchecked(target, f, forceMapZeros, skipClearing); + } + + internal virtual void MapToUnchecked(VectorStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) + where TU : struct, IEquatable, IFormattable { for (int i = 0; i < Length; i++) { - At(i, f(At(i))); + target.At(i, f(At(i))); } } - public virtual void MapIndexedInplace(Func f, bool forceMapZeros = false) + public void MapIndexedTo(VectorStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) + where TU : struct, IEquatable, IFormattable + { + if (target == null) + { + throw new ArgumentNullException("target"); + } + + if (Length != target.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "target"); + } + + MapIndexedToUnchecked(target, f, forceMapZeros, skipClearing); + } + + internal virtual void MapIndexedToUnchecked(VectorStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) + where TU : struct, IEquatable, IFormattable { for (int i = 0; i < Length; i++) { - At(i, f(i, At(i))); + target.At(i, f(i, At(i))); } } } diff --git a/src/Numerics/LinearAlgebra/Vector.cs b/src/Numerics/LinearAlgebra/Vector.cs index ba8199b1..1089e4c9 100644 --- a/src/Numerics/LinearAlgebra/Vector.cs +++ b/src/Numerics/LinearAlgebra/Vector.cs @@ -318,7 +318,7 @@ namespace MathNet.Numerics.LinearAlgebra /// public void MapInplace(Func f, bool forceMapZeros = false) { - Storage.MapInplace(f, forceMapZeros); + Storage.MapToUnchecked(Storage, f, forceMapZeros, skipClearing: true); } /// @@ -329,7 +329,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 vector and replaces the value in the result vector. + /// 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 vectors). + /// + public void Map(Func f, Vector 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 vector and replaces the value in the result vector. + /// 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 vectors). + /// + public void MapIndexed(Func f, Vector 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 vector and returns the results as a new vector. + /// 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 vectors). + /// + public Vector Map(Func f, bool forceMapZeros = false) + where TU : struct, IEquatable, IFormattable + { + var result = Vector.Build.SameAs(this); + Storage.MapToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true); + return result; + } + + /// + /// Applies a function to each value of this vector and returns the results as a new vector. + /// 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 vectors). + /// + public Vector MapIndexed(Func f, bool forceMapZeros = false) + where TU : struct, IEquatable, IFormattable + { + var result = Vector.Build.SameAs(this); + Storage.MapIndexedToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true); + return result; } } }