LA: Vector.Map2, storage-optimized

src/Numerics/LinearAlgebra/Double/SparseVector.cs

@@ -153,7 +153,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
             if (ReferenceEquals(this, result))
             {
-                //populate a new vector with the scalar
+                // populate a new vector with the scalar
                 var vnonZeroValues = new double[Count];
                 var vnonZeroIndices = new int[Count];
                 for (int index = 0; index < Count; index++)
@@ -162,7 +162,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                     vnonZeroValues[index] = scalar;
                 }
 
-                //populate the non zero values from this
+                // populate the non zero values from this
                 var indices = _storage.Indices;
                 var values = _storage.Values;
                 for (int j = 0; j < _storage.ValueCount; j++)
@@ -170,7 +170,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                     vnonZeroValues[indices[j]] = values[j] + scalar;
                 }
 
-                //assign this vectors arrary to the new arrays.
+                // assign this vectors array to the new arrays.
                 _storage.Values = vnonZeroValues;
                 _storage.Indices = vnonZeroIndices;
                 _storage.ValueCount = Count;

src/Numerics/LinearAlgebra/Storage/DenseVectorStorage.cs

@@ -172,6 +172,29 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 {
                     Array.Copy(Data, 0, denseTarget.Data, 0, Data.Length);
                 }
+
+                return;
+            }
+
+            var sparseTarget = target as SparseVectorStorage<T>;
+            if (sparseTarget != null)
+            {
+                var indices = new List<int>();
+                var values = new List<T>();
+
+                for (int i = 0; i < Data.Length; i++)
+                {
+                    var item = Data[i];
+                    if (!Zero.Equals(item))
+                    {
+                        values.Add(item);
+                        indices.Add(i);
+                    }
+                }
+
+                sparseTarget.Indices = indices.ToArray();
+                sparseTarget.Values = values.ToArray();
+                sparseTarget.ValueCount = values.Count;
                 return;
             }
 
@@ -364,5 +387,60 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 target.At(i, f(i, Data[i]));
             }
         }
+
+        internal override void Map2ToUnchecked(VectorStorage<T> target, VectorStorage<T> other, Func<T, T, T> f, Zeros zeros = Zeros.AllowSkip, ExistingData existingData = ExistingData.Clear)
+        {
+            if (target is SparseVectorStorage<T>)
+            {
+                // Recursive to dense target at first, since the operation is
+                // effectively dense anyway because at least one operand is dense
+                var intermediate = new DenseVectorStorage<T>(target.Length);
+                Map2ToUnchecked(intermediate, other, f, zeros, ExistingData.AssumeZeros);
+                intermediate.CopyTo(target, existingData);
+                return;
+            }
+
+            var denseTarget = target as DenseVectorStorage<T>;
+            var denseOther = other as DenseVectorStorage<T>;
+            if (denseTarget != null && denseOther != null)
+            {
+                CommonParallel.For(0, Data.Length, 4096, (a, b) =>
+                {
+                    for (int i = a; i < b; i++)
+                    {
+                        denseTarget.Data[i] = f(Data[i], denseOther.Data[i]);
+                    }
+                });
+
+                return;
+            }
+
+            var sparseOther = other as SparseVectorStorage<T>;
+            if (denseTarget != null && sparseOther != null)
+            {
+                T[] targetData = denseTarget.Data;
+                int[] otherIndices = sparseOther.Indices;
+                T[] otherValues = sparseOther.Values;
+                int otherValueCount = sparseOther.ValueCount;
+
+                int k = 0;
+                for (int i = 0; i < Data.Length; i++)
+                {
+                    if (k < otherValueCount && otherIndices[k] == i)
+                    {
+                        targetData[i] = f(Data[i], otherValues[k]);
+                        k++;
+                    }
+                    else
+                    {
+                        targetData[i] = f(Data[i], Zero);
+                    }
+                }
+
+                return;
+            }
+
+            base.Map2ToUnchecked(target, other, f, zeros, existingData);
+        }
     }
 }

src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs

@@ -30,6 +30,7 @@
 
 using System;
 using System.Collections.Generic;
+using System.Diagnostics;
 using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Threading;
@@ -768,5 +769,156 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
             base.MapIndexedToUnchecked(target, f, zeros, existingData);
         }
+
+        internal override void Map2ToUnchecked(VectorStorage<T> target, VectorStorage<T> other, Func<T, T, T> f, Zeros zeros = Zeros.AllowSkip, ExistingData existingData = ExistingData.Clear)
+        {
+            var processZeros = zeros == Zeros.Include || !Zero.Equals(f(Zero, Zero));
+
+            var denseTarget = target as DenseVectorStorage<T>;
+            var denseOther = other as DenseVectorStorage<T>;
+
+            if (denseTarget == null && (denseOther != null || processZeros))
+            {
+                // The handling is effectively dense but we're supposed to push
+                // to a sparse target. Let's use a dense target instead,
+                // then copy it normalized back to the sparse target.
+                var intermediate = new DenseVectorStorage<T>(target.Length);
+                Map2ToUnchecked(intermediate, other, f, zeros, ExistingData.AssumeZeros);
+                intermediate.CopyTo(target, existingData);
+                return;
+            }
+
+            if (denseOther != null)
+            {
+                T[] targetData = denseTarget.Data;
+                T[] otherData = denseOther.Data;
+
+                int k = 0;
+                for (int i = 0; i < otherData.Length; i++)
+                {
+                    if (k < ValueCount && Indices[k] == i)
+                    {
+                        targetData[i] = f(Values[k], otherData[i]);
+                        k++;
+                    }
+                    else
+                    {
+                        targetData[i] = f(Zero, otherData[i]);
+                    }
+                }
+
+                return;
+            }
+
+            var sparseOther = other as SparseVectorStorage<T>;
+            if (sparseOther != null && denseTarget != null)
+            {
+                T[] targetData = denseTarget.Data;
+                int[] otherIndices = sparseOther.Indices;
+                T[] otherValues = sparseOther.Values;
+                int otherValueCount = sparseOther.ValueCount;
+
+                if (processZeros)
+                {
+                    int p = 0, q = 0;
+                    for (int i = 0; i < targetData.Length; i++)
+                    {
+                        var left = p < ValueCount && Indices[p] == i ? Values[p++] : Zero;
+                        var right = q < otherValueCount && otherIndices[q] == i ? otherValues[q++] : Zero;
+                        targetData[i] = f(left, right);
+                    }
+                }
+                else
+                {
+                    if (existingData == ExistingData.Clear)
+                    {
+                        denseTarget.Clear();
+                    }
+
+                    int p = 0, q = 0;
+                    while (p < ValueCount || q < otherValueCount)
+                    {
+                        if (q >= otherValueCount || p < ValueCount && Indices[p] < otherIndices[q])
+                        {
+                            targetData[Indices[p]] = f(Values[p], Zero);
+                            p++;
+                        }
+                        else if (p >= ValueCount || q < otherValueCount && Indices[p] > otherIndices[q])
+                        {
+                            targetData[otherIndices[q]] = f(Zero, otherValues[q]);
+                            q++;
+                        }
+                        else
+                        {
+                            Debug.Assert(Indices[p] == otherIndices[q]);
+                            targetData[Indices[p]] = f(Values[p], otherValues[q]);
+                            p++;
+                            q++;
+                        }
+                    }
+                }
+
+                return;
+            }
+
+            var sparseTarget = target as SparseVectorStorage<T>;
+            if (sparseOther != null && sparseTarget != null)
+            {
+                var indices = new List<int>();
+                var values = new List<T>();
+                int[] otherIndices = sparseOther.Indices;
+                T[] otherValues = sparseOther.Values;
+                int otherValueCount = sparseOther.ValueCount;
+
+                int p = 0, q = 0;
+                while (p < ValueCount || q < otherValueCount)
+                {
+                    if (q >= otherValueCount || p < ValueCount && Indices[p] < otherIndices[q])
+                    {
+                        var value = f(Values[p], Zero);
+                        if (!Zero.Equals(value))
+                        {
+                            indices.Add(Indices[p]);
+                            values.Add(value);
+                        }
+
+                        p++;
+                    }
+                    else if (p >= ValueCount || q < otherValueCount && Indices[p] > otherIndices[q])
+                    {
+                        var value = f(Zero, otherValues[q]);
+                        if (!Zero.Equals(value))
+                        {
+                            indices.Add(otherIndices[q]);
+                            values.Add(value);
+                        }
+
+                        q++;
+                    }
+                    else
+                    {
+                        Debug.Assert(Indices[p] == otherIndices[q]);
+                        var value = f(Values[p], otherValues[q]);
+                        if (!Zero.Equals(value))
+                        {
+                            indices.Add(Indices[p]);
+                            values.Add(value);
+                        }
+
+                        p++;
+                        q++;
+                    }
+                }
+
+                sparseTarget.Indices = indices.ToArray();
+                sparseTarget.Values = values.ToArray();
+                sparseTarget.ValueCount = values.Count;
+                return;
+            }
+
+            // FALL BACK
+
+            base.Map2ToUnchecked(target, other, f, zeros, existingData);
+        }
     }
 }

src/Numerics/LinearAlgebra/Storage/VectorStorage.cs

@@ -469,5 +469,40 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 target.At(i, f(i, At(i)));
             }
         }
+
+        public void Map2To(VectorStorage<T> target, VectorStorage<T> other, Func<T, T, T> f,
+            Zeros zeros = Zeros.AllowSkip, ExistingData existingData = ExistingData.Clear)
+        {
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            if (other == null)
+            {
+                throw new ArgumentNullException("other");
+            }
+
+            if (Length != target.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength, "target");
+            }
+
+            if (Length != other.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
+            }
+
+            Map2ToUnchecked(target, other, f, zeros, existingData);
+        }
+
+        internal virtual void Map2ToUnchecked(VectorStorage<T> target, VectorStorage<T> other, Func<T, T, T> f,
+            Zeros zeros = Zeros.AllowSkip, ExistingData existingData = ExistingData.Clear)
+        {
+            for (int i = 0; i < Length; i++)
+            {
+                target.At(i, f(At(i), other.At(i)));
+            }
+        }
     }
 }

src/Numerics/LinearAlgebra/Vector.cs

@@ -457,5 +457,23 @@ namespace MathNet.Numerics.LinearAlgebra
             Storage.MapIndexedToUnchecked(result.Storage, f, zeros, ExistingData.AssumeZeros);
             return result;
         }
+
+        /// <summary>
+        /// Applies a function to each value pair of two vectors and replaces the value in the result vector.
+        /// </summary>
+        public void Map2(Func<T, T, T> f, Vector<T> other, Vector<T> result, Zeros zeros = Zeros.AllowSkip)
+        {
+            Storage.Map2To(result.Storage, other.Storage, f, zeros, ExistingData.Clear);
+        }
+
+        /// <summary>
+        /// Applies a function to each value pair of two vectors and returns the results as a new vector.
+        /// </summary>
+        public Vector<T> Map2(Func<T, T, T> f, Vector<T> other, Zeros zeros = Zeros.AllowSkip)
+        {
+            var result = Build.SameAs(this);
+            Storage.Map2To(result.Storage, other.Storage, f, zeros, ExistingData.AssumeZeros);
+            return result;
+        }
     }
 }

src/UnitTests/LinearAlgebraTests/VectorStorageTests.cs

@@ -0,0 +1,184 @@
+// <copyright file="VectorStorageTests.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Storage;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.LinearAlgebraTests
+{
+    [TestFixture, Category("LA")]
+    public class VectorStorageTests
+    {
+        [Test]
+        public void Map2SkipZeros()
+        {
+            double[] a = { 1.0, 2.0, 0.0, 4.0, 0.0, 6.0 };
+            double[] b = { 11.0, 12.0, 13.0, 0.0, 0.0, 16.0 };
+            var adense = DenseVectorStorage<double>.OfEnumerable(a);
+            var asparse = SparseVectorStorage<double>.OfEnumerable(a);
+            var bdense = DenseVectorStorage<double>.OfEnumerable(b);
+            var bsparse = SparseVectorStorage<double>.OfEnumerable(b);
+
+            var rdense = new DenseVectorStorage<double>(a.Length);
+            var rsparse = new SparseVectorStorage<double>(a.Length);
+
+            var expected = new DenseVectorStorage<double>(6, new[] { 12.0, 14.0, 13.0, 4.0, 0.0, 22.0 });
+
+            rdense.Clear();
+            adense.Map2To(rdense, bdense, (u, v) => u+v, Zeros.AllowSkip);
+            Assert.That(rdense.Equals(expected), "dense*dense->dense");
+
+            rsparse.Clear();
+            adense.Map2To(rsparse, bdense, (u, v) => u+v, Zeros.AllowSkip);
+            Assert.That(rsparse.Equals(expected), "dense*dense->sparse");
+
+            rdense.Clear();
+            adense.Map2To(rdense, bsparse, (u, v) => u+v, Zeros.AllowSkip);
+            Assert.That(rdense.Equals(expected), "dense*sparse->dense");
+
+            rsparse.Clear();
+            adense.Map2To(rsparse, bsparse, (u, v) => u+v, Zeros.AllowSkip);
+            Assert.That(rsparse.Equals(expected), "dense*sparse->sparse");
+
+            rdense.Clear();
+            asparse.Map2To(rdense, bdense, (u, v) => u+v, Zeros.AllowSkip);
+            Assert.That(rdense.Equals(expected), "sparse*dense->dense");
+
+            rsparse.Clear();
+            asparse.Map2To(rsparse, bdense, (u, v) => u+v, Zeros.AllowSkip);
+            Assert.That(rsparse.Equals(expected), "sparse*dense->sparse");
+
+            rdense.Clear();
+            asparse.Map2To(rdense, bsparse, (u, v) => u+v, Zeros.AllowSkip);
+            Assert.That(rdense.Equals(expected), "sparse*sparse->dense");
+
+            rsparse.Clear();
+            asparse.Map2To(rsparse, bsparse, (u, v) => u+v, Zeros.AllowSkip);
+            Assert.That(rsparse.Equals(expected), "sparse*sparse->sparse");
+        }
+
+        [Test]
+        public void Map2ForceIncludeZeros()
+        {
+            double[] a = { 1.0, 2.0, 0.0, 4.0, 0.0, 6.0 };
+            double[] b = { 11.0, 12.0, 13.0, 0.0, 0.0, 16.0 };
+            var adense = DenseVectorStorage<double>.OfEnumerable(a);
+            var asparse = SparseVectorStorage<double>.OfEnumerable(a);
+            var bdense = DenseVectorStorage<double>.OfEnumerable(b);
+            var bsparse = SparseVectorStorage<double>.OfEnumerable(b);
+
+            var rdense = new DenseVectorStorage<double>(a.Length);
+            var rsparse = new SparseVectorStorage<double>(a.Length);
+
+            var expected = new DenseVectorStorage<double>(6, new[] { 13.0, 15.0, 14.0, 5.0, 1.0, 23.0 });
+
+            rdense.Clear();
+            adense.Map2To(rdense, bdense, (u, v) => u+v+1.0, Zeros.Include);
+            Assert.That(rdense.Equals(expected), "dense*dense->dense");
+
+            rsparse.Clear();
+            adense.Map2To(rsparse, bdense, (u, v) => u+v+1.0, Zeros.Include);
+            Assert.That(rsparse.Equals(expected), "dense*dense->sparse");
+
+            rdense.Clear();
+            adense.Map2To(rdense, bsparse, (u, v) => u+v+1.0, Zeros.Include);
+            Assert.That(rdense.Equals(expected), "dense*sparse->dense");
+
+            rsparse.Clear();
+            adense.Map2To(rsparse, bsparse, (u, v) => u+v+1.0, Zeros.Include);
+            Assert.That(rsparse.Equals(expected), "dense*sparse->sparse");
+
+            rdense.Clear();
+            asparse.Map2To(rdense, bdense, (u, v) => u+v+1.0, Zeros.Include);
+            Assert.That(rdense.Equals(expected), "sparse*dense->dense");
+
+            rsparse.Clear();
+            asparse.Map2To(rsparse, bdense, (u, v) => u+v+1.0, Zeros.Include);
+            Assert.That(rsparse.Equals(expected), "sparse*dense->sparse");
+
+            rdense.Clear();
+            asparse.Map2To(rdense, bsparse, (u, v) => u+v+1.0, Zeros.Include);
+            Assert.That(rdense.Equals(expected), "sparse*sparse->dense");
+
+            rsparse.Clear();
+            asparse.Map2To(rsparse, bsparse, (u, v) => u+v+1.0, Zeros.Include);
+            Assert.That(rsparse.Equals(expected), "sparse*sparse->sparse");
+        }
+
+        [Test]
+        public void Map2AutoIncludeZeros()
+        {
+            double[] a = { 1.0, 2.0, 0.0, 4.0, 0.0, 6.0 };
+            double[] b = { 11.0, 12.0, 13.0, 0.0, 0.0, 16.0 };
+            var adense = DenseVectorStorage<double>.OfEnumerable(a);
+            var asparse = SparseVectorStorage<double>.OfEnumerable(a);
+            var bdense = DenseVectorStorage<double>.OfEnumerable(b);
+            var bsparse = SparseVectorStorage<double>.OfEnumerable(b);
+
+            var rdense = new DenseVectorStorage<double>(a.Length);
+            var rsparse = new SparseVectorStorage<double>(a.Length);
+
+            var expected = new DenseVectorStorage<double>(6, new[] { 13.0, 15.0, 14.0, 5.0, 1.0, 23.0 });
+
+            rdense.Clear();
+            adense.Map2To(rdense, bdense, (u, v) => u+v+1.0, Zeros.AllowSkip);
+            Assert.That(rdense.Equals(expected), "dense*dense->dense");
+
+            rsparse.Clear();
+            adense.Map2To(rsparse, bdense, (u, v) => u+v+1.0, Zeros.AllowSkip);
+            Assert.That(rsparse.Equals(expected), "dense*dense->sparse");
+
+            rdense.Clear();
+            adense.Map2To(rdense, bsparse, (u, v) => u+v+1.0, Zeros.AllowSkip);
+            Assert.That(rdense.Equals(expected), "dense*sparse->dense");
+
+            rsparse.Clear();
+            adense.Map2To(rsparse, bsparse, (u, v) => u+v+1.0, Zeros.AllowSkip);
+            Assert.That(rsparse.Equals(expected), "dense*sparse->sparse");
+
+            rdense.Clear();
+            asparse.Map2To(rdense, bdense, (u, v) => u+v+1.0, Zeros.AllowSkip);
+            Assert.That(rdense.Equals(expected), "sparse*dense->dense");
+
+            rsparse.Clear();
+            asparse.Map2To(rsparse, bdense, (u, v) => u+v+1.0, Zeros.AllowSkip);
+            Assert.That(rsparse.Equals(expected), "sparse*dense->sparse");
+
+            rdense.Clear();
+            asparse.Map2To(rdense, bsparse, (u, v) => u+v+1.0, Zeros.AllowSkip);
+            Assert.That(rdense.Equals(expected), "sparse*sparse->dense");
+
+            rsparse.Clear();
+            asparse.Map2To(rsparse, bsparse, (u, v) => u+v+1.0, Zeros.AllowSkip);
+            Assert.That(rsparse.Equals(expected), "sparse*sparse->sparse");
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -350,6 +350,7 @@
     <Compile Include="EuclidTests\GcdRelatedTest.cs" />
     <Compile Include="EuclidTests\GcdRelatedTestBigInteger.cs" />
     <Compile Include="EuclidTests\IntegerTheoryTest.cs" />
+    <Compile Include="LinearAlgebraTests\VectorStorageTests.cs" />
     <Compile Include="Random\SystemRandomSourceTests.cs" />
     <Compile Include="RootFindingTests\BisectionTest.cs" />
     <Compile Include="PermutationTest.cs" />