LA: storage-aware non-inplace map on vectors

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<_>) =
+    let inline mapInPlace f (v: #Vector<_>) = v.MapInplace((fun x -> f x), true)
-        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<_>) =
+    let inline mapiInPlace f (v: #Vector<_>) = v.MapIndexedInplace((fun i x -> f i x), true)
-        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<_>) =
+    let inline mapSkipZerosInPlace f (v: #Vector<_>) = v.MapInplace((fun x -> f x), false)
-        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<float>) =
+    let inline mapiSkipZerosInPlace f (v: #Vector<_>) = v.MapIndexedInplace((fun i x -> f i x), false)
-        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 inline map f (v: #Vector<_>) = v.Map((fun x -> f x), true)
-        let w = v.Clone()
-        w.MapInplace((fun x -> f x), true)
-        w
 
     /// 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 inline mapSkipZeros f (v: #Vector<_>) = v.Map((fun x -> f x), false)
-        let w = v.Clone()
-        w.MapInplace((fun x -> f x), false)
-        w
 
     /// Maps a vector to a new vector by applying a function to every element.
-    let inline mapi f (v: #Vector<_>) =
+    let inline mapi f (v: #Vector<_>) = v.MapIndexed((fun i x -> f i x), true)
-        let w = v.Clone()
-        w.MapIndexedInplace((fun i x -> f i x), true)
-        w
 
     /// 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 inline mapiSkipZeros f (v: #Vector<_>) = v.MapIndexed((fun i x -> f i x), false)
-        let w = v.Clone()
-        w.MapIndexedInplace((fun i x -> f i x), false)
-        w
 
 
 

src/Numerics/LinearAlgebra/Builder.cs

@@ -389,7 +389,8 @@ namespace MathNet.Numerics.LinearAlgebra
         /// <summary>
         /// Create a new matrix with the same kind of the provided example.
         /// </summary>
-        public Matrix<T> SameAs(Matrix<T> example, int rows, int columns, bool fullyMutable = false)
+        public Matrix<T> SameAs<TU>(Matrix<TU> example, int rows, int columns, bool fullyMutable = false)
+            where TU : struct, IEquatable<TU>, IFormattable
         {
             var storage = example.Storage;
             if (storage is DenseColumnMajorMatrixStorage<T>) return Dense(rows, columns);
@@ -401,7 +402,8 @@ namespace MathNet.Numerics.LinearAlgebra
         /// <summary>
         /// Create a new matrix with the same kind and dimensions of the provided example.
         /// </summary>
-        public Matrix<T> SameAs(Matrix<T> example)
+        public Matrix<T> SameAs<TU>(Matrix<TU> example)
+            where TU : struct, IEquatable<TU>, IFormattable
         {
             return SameAs(example, example.RowCount, example.ColumnCount);
         }
@@ -1338,7 +1340,8 @@ namespace MathNet.Numerics.LinearAlgebra
         /// <summary>
         /// Create a new vector with the same kind of the provided example.
         /// </summary>
-        public Vector<T> SameAs(Vector<T> example, int length)
+        public Vector<T> SameAs<TU>(Vector<TU> example, int length)
+            where TU : struct, IEquatable<TU>, IFormattable
         {
             return example.Storage.IsDense ? Dense(length) : Sparse(length);
         }
@@ -1346,7 +1349,8 @@ namespace MathNet.Numerics.LinearAlgebra
         /// <summary>
         /// Create a new vector with the same kind and dimension of the provided example.
         /// </summary>
-        public Vector<T> SameAs(Vector<T> example)
+        public Vector<T> SameAs<TU>(Vector<TU> example)
+            where TU : struct, IEquatable<TU>, IFormattable
         {
             return example.Storage.IsDense ? Dense(example.Count) : Sparse(example.Count);
         }

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<T>(length, data);
         }
 
@@ -212,7 +212,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             var denseTarget = target as DenseColumnMajorMatrixStorage<T>;
             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<T, T> f, bool forceMapZeros = false)
+        internal override void MapToUnchecked<TU>(VectorStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
         {
-            CommonParallel.For(0, Data.Length, 4096, (a, b) =>
+            var denseTarget = target as DenseVectorStorage<TU>;
+            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<int, T, T> f, bool forceMapZeros = false)
+        internal override void MapIndexedToUnchecked<TU>(VectorStorage<TU> target, Func<int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
         {
-            CommonParallel.For(0, Data.Length, 4096, (a, b) =>
+            var denseTarget = target as DenseVectorStorage<TU>;
+            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]));
+            }
         }
     }
 }

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<T, T> f, bool forceMapZeros = false)
+        internal override void MapToUnchecked<TU>(VectorStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
         {
-            var indices = new List<int>();
+            var sparseTarget = target as SparseVectorStorage<TU>;
-            var values = new List<T>();
+            if (sparseTarget != null)
-            if (forceMapZeros || !Zero.Equals(f(Zero)))
             {
-                int k = 0;
+                var indices = new List<int>();
-                for (int i = 0; i < Length; i++)
+                var values = new List<TU>();
+                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);
+                    for (int i = 0; i < ValueCount; i++)
-                    if (!Zero.Equals(item))
                     {
-                        values.Add(item);
+                        var item = f(Values[i]);
-                        indices.Add(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<TU>;
+            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]);
+                    int k = 0;
-                    if (!Zero.Equals(item))
+                    for (int i = 0; i < Length; i++)
                     {
-                        values.Add(item);
+                        denseTarget.Data[i] = k < ValueCount && (Indices[k]) == i
-                        indices.Add(Indices[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();
+            // FALL BACK
-            ValueCount = values.Count;
+
+            base.MapToUnchecked(target, f, forceMapZeros, skipClearing);
         }
 
-        public override void MapIndexedInplace(Func<int, T, T> f, bool forceMapZeros = false)
+        internal override void MapIndexedToUnchecked<TU>(VectorStorage<TU> target, Func<int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
         {
-            var indices = new List<int>();
+            var sparseTarget = target as SparseVectorStorage<TU>;
-            var values = new List<T>();
+            if (sparseTarget != null)
-            if (forceMapZeros || !Zero.Equals(f(0, Zero)))
             {
-                int k = 0;
+                var indices = new List<int>();
-                for (int i = 0; i < Length; i++)
+                var values = new List<TU>();
+                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);
+                    for (int i = 0; i < ValueCount; i++)
-                    if (!Zero.Equals(item))
                     {
-                        values.Add(item);
+                        var item = f(Indices[i], Values[i]);
-                        indices.Add(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<TU>;
+            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]);
+                    int k = 0;
-                    if (!Zero.Equals(item))
+                    for (int i = 0; i < Length; i++)
                     {
-                        values.Add(item);
+                        denseTarget.Data[i] = k < ValueCount && (Indices[k]) == i
-                        indices.Add(Indices[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();
+            // FALL BACK
-            ValueCount = values.Count;
+
+            base.MapIndexedToUnchecked(target, f, forceMapZeros, skipClearing);
         }
     }
 }

src/Numerics/LinearAlgebra/Storage/VectorStorage.cs

@@ -401,19 +401,53 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
         // FUNCTIONAL COMBINATORS
 
-        public virtual void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        public void MapTo<TU>(VectorStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
+            where TU : struct, IEquatable<TU>, 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<TU>(VectorStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
+            where TU : struct, IEquatable<TU>, IFormattable
         {
             for (int i = 0; i < Length; i++)
             {
-                At(i, f(At(i)));
+                target.At(i, f(At(i)));
             }
         }
 
-        public virtual void MapIndexedInplace(Func<int, T, T> f, bool forceMapZeros = false)
+        public void MapIndexedTo<TU>(VectorStorage<TU> target, Func<int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
+            where TU : struct, IEquatable<TU>, 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<TU>(VectorStorage<TU> target, Func<int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
+            where TU : struct, IEquatable<TU>, IFormattable
         {
             for (int i = 0; i < Length; i++)
             {
-                At(i, f(i, At(i)));
+                target.At(i, f(i, At(i)));
             }
         }
     }

src/Numerics/LinearAlgebra/Vector.cs

@@ -318,7 +318,7 @@ namespace MathNet.Numerics.LinearAlgebra
         /// </summary>
         public void MapInplace(Func<T, T> f, bool forceMapZeros = false)
         {
-            Storage.MapInplace(f, forceMapZeros);
+            Storage.MapToUnchecked(Storage, f, forceMapZeros, skipClearing: true);
         }
 
         /// <summary>
@@ -329,7 +329,57 @@ namespace MathNet.Numerics.LinearAlgebra
         /// </summary>
         public void MapIndexedInplace(Func<int, T, T> f, bool forceMapZeros = false)
         {
-            Storage.MapIndexedInplace(f, forceMapZeros);
+            Storage.MapIndexedToUnchecked(Storage, f, forceMapZeros, skipClearing: true);
+        }
+
+        /// <summary>
+        /// 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).
+        /// </summary>
+        public void Map<TU>(Func<T, TU> f, Vector<TU> result, bool forceMapZeros = false)
+            where TU : struct, IEquatable<TU>, IFormattable
+        {
+            Storage.MapTo(result.Storage, f, forceMapZeros, skipClearing: forceMapZeros);
+        }
+
+        /// <summary>
+        /// 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).
+        /// </summary>
+        public void MapIndexed<TU>(Func<int, T, TU> f, Vector<TU> result, bool forceMapZeros = false)
+            where TU : struct, IEquatable<TU>, IFormattable
+        {
+            Storage.MapIndexedTo(result.Storage, f, forceMapZeros, skipClearing: forceMapZeros);
+        }
+
+        /// <summary>
+        /// 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).
+        /// </summary>
+        public Vector<TU> Map<TU>(Func<T, TU> f, bool forceMapZeros = false)
+            where TU : struct, IEquatable<TU>, IFormattable
+        {
+            var result = Vector<TU>.Build.SameAs(this);
+            Storage.MapToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true);
+            return result;
+        }
+
+        /// <summary>
+        /// 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).
+        /// </summary>
+        public Vector<TU> MapIndexed<TU>(Func<int, T, TU> f, bool forceMapZeros = false)
+            where TU : struct, IEquatable<TU>, IFormattable
+        {
+            var result = Vector<TU>.Build.SameAs(this);
+            Storage.MapIndexedToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true);
+            return result;
         }
     }
 }