LA: Functional Map/MapIndexed on matrices, with F# support

src/FSharp/LinearAlgebra.Double.Matrix.fs

@@ -36,6 +36,82 @@ open MathNet.Numerics.LinearAlgebra.Generic
 [<CompilationRepresentation(CompilationRepresentationFlags.ModuleSuffix)>]
 module Matrix =
 
+    /// Transform a vector into a 2D array.
+    let inline toArray2 (A: #Matrix<float>) = A.ToArray()
+
+    /// In-place map of every matrix element using a function.
+    let inline mapInPlace (f: float -> float) (A: #Matrix<float>) =
+        A.MapInplace((fun x -> f x), true)
+
+    /// In-place map of every matrix element using a position dependent function.
+    let inline mapiInPlace (f: int -> int -> float -> float) (A: #Matrix<float>) =
+        A.MapIndexedInplace((fun i j x -> f i j x), true)
+
+    /// In-place map of every matrix element using a function.
+    /// Zero-values may be skipped (relevant mostly for sparse matrices).
+    let inline mapnzInPlace (f: float -> float) (A: #Matrix<float>) =
+        A.MapInplace((fun x -> f x), false)
+
+    /// In-place map of every matrix element using a position dependent function.
+    /// Zero-values may be skipped (relevant mostly for sparse matrices).
+    let inline mapinzInPlace (f: int -> int -> float -> float) (A: #Matrix<float>) =
+        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<float> -> Vector<float>) (A: #Matrix<float>) =
+        for j = 0 to A.ColumnCount-1 do
+            A.SetColumn(j, f j (A.Column(j)))
+
+    /// In-place map every matrix row using the given position dependent function.
+    let inline mapRowsInPlace (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
+        for i = 0 to A.RowCount-1 do
+            A.SetRow(i, f i (A.Row(i)))
+
+    [<System.ObsoleteAttribute("Use mapiInPlace instead. Scheduled for removal in v3.0.")>]
+    let inplaceMapi = mapiInPlace
+    [<System.ObsoleteAttribute("Use mapColsInPlace instead. Scheduled for removal in v3.0.")>]
+    let inplaceMapCols = mapColsInPlace
+    [<System.ObsoleteAttribute("Use mapRowsInPlace instead. Scheduled for removal in v3.0.")>]
+    let inplaceMapRows = mapRowsInPlace
+
+    /// Map every matrix element using the given function.
+    let inline map (f: float -> float) (A: #Matrix<float>) =
+        let A = A.Clone()
+        A.MapInplace((fun x -> f x), true)
+        A
+
+    /// Map every matrix element using the given function.
+    /// Zero-values may be skipped (relevant mostly for sparse matrices).
+    let inline mapnz (f: float -> float) (A: #Matrix<float>) =
+        let A = A.Clone()
+        A.MapInplace((fun x -> f x), false)
+        A
+
+    /// Map every matrix element using the given position dependent function.
+    let inline mapi (f: int -> int -> float -> float) (A: #Matrix<float>) =
+        let A = A.Clone()
+        A.MapIndexedInplace((fun i j x -> f i j x), true)
+        A
+
+    /// Map every matrix element using the given position dependent function.
+    /// Zero-values may be skipped (relevant mostly for sparse matrices).
+    let inline mapinz (f: int -> int -> float -> float) (A: #Matrix<float>) =
+        let A = A.Clone()
+        A.MapIndexedInplace((fun i j x -> f i j x), false)
+        A
+
+    /// Map every matrix column using the given position dependent function.
+    let inline mapCols (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
+        let A = A.Clone()
+        mapColsInPlace f A
+        A
+
+    /// Map every matrix row using the given position dependent function.
+    let inline mapRows (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
+        let A = A.Clone()
+        mapRowsInPlace f A
+        A
+
     /// Fold a function over all matrix elements.
     let inline fold (f: 'a -> float -> 'a) (acc0: 'a) (A: #Matrix<float>) =
         let n = A.RowCount
@@ -66,12 +142,6 @@ module Matrix =
                 acc <- f i j acc (A.At(i,j))
         acc
 
-    /// Create a 2D array from a matrix.
-    let inline toArray2 (A: #Matrix<float>) =
-        let n = A.RowCount
-        let m = A.ColumnCount
-        Array2D.init n m (fun i j -> (A.Item(i,j)))
-
     /// Checks whether a predicate holds for all elements of a matrix.
     let inline forall (p: float -> bool) (A: #Matrix<float>) =
         let mutable b = true
@@ -116,61 +186,9 @@ module Matrix =
             if j = A.ColumnCount then i <- i+1; j <- 0
         b
 
-    /// Map every matrix element using the given function.
-    let inline map (f: float -> float) (A: #Matrix<float>) =
-        let N = A.RowCount
-        let M = A.ColumnCount
-        let C = A.Clone()
-        for i=0 to N-1 do
-            for j=0 to M-1 do
-                C.At(i, j, f (C.At(i,j)))
-        C
-
-    /// Map every matrix element using the given position dependent function.
-    let inline mapi (f: int -> int -> float -> float) (A: #Matrix<float>) =
-        let N = A.RowCount
-        let M = A.ColumnCount
-        let C = A.Clone()
-        for i=0 to N-1 do
-            for j=0 to M-1 do
-                C.At(i, j, f i j (C.At(i,j)))
-        C
-
-    /// In-place map every matrix column using the given position dependent function.
-    let inline inplaceMapCols (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
-        for j = 0 to A.ColumnCount-1 do
-            A.SetColumn(j, f j (A.Column(j)))
-        ()
-
-    /// In-place map every matrix row using the given position dependent function.
-    let inline inplaceMapRows (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
-        for i = 0 to A.RowCount-1 do
-            A.SetRow(i, f i (A.Row(i)))
-        ()
-
-    /// Map every matrix column using the given position dependent function.
-    let inline mapCols (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
-        let A = A.Clone()
-        inplaceMapCols f A
-        A
-
-    /// Map every matrix row using the given position dependent function.
-    let inline mapRows (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
-        let A = A.Clone()
-        inplaceMapRows f A
-        A
-
     /// In-place assignment.
     let inline inplaceAssign (f: int -> int -> float) (A: #Matrix<float>) =
-        for i=0 to A.RowCount-1 do
-            for j=0 to A.ColumnCount-1 do
-                A.At(i, j, f i j)
-
-    /// In-place map of every matrix element using a position dependent function.
-    let inline inplaceMapi (f: int -> int -> float -> float) (A: #Matrix<float>) =
-        for i=0 to A.RowCount-1 do
-            for j=0 to A.ColumnCount-1 do
-                A.At(i, j, f i j (A.At(i,j)))
+        A.MapIndexedInplace((fun i j x -> f i j), true)
 
     /// Creates a sequence that iterates the non-zero entries in the matrix.
     let inline nonZeroEntries (A: #Matrix<float>) =
@@ -199,14 +217,12 @@ module Matrix =
         for i=0 to A.RowCount-1 do
             for j=0 to A.ColumnCount-1 do
                 f (A.At(i,j))
-        ()
 
     /// Iterates over all elements of a matrix using the element indices.
     let inline iteri (f: int -> int -> float -> unit) (A: #Matrix<float>) =
         for i=0 to A.RowCount-1 do
             for j=0 to A.ColumnCount-1 do
                 f i j (A.At(i,j))
-        ()
 
     /// Fold one column.
     let inline foldCol (f: 'a -> float -> 'a) acc (A: #Matrix<float>) k =
@@ -298,10 +314,7 @@ module DenseMatrix =
 
     /// Create a matrix with the given entries.
     [<System.ObsoleteAttribute("Use ofSeqi instead. Scheduled for removal in v3.0.")>]
-    let inline initDense (rows: int) (cols: int) (es: #seq<int * int * float>) =
-        let A = new DenseMatrix(rows,cols)
-        Seq.iter (fun (i,j,f) -> A.At(i,j,f)) es
-        A
+    let inline initDense (rows: int) (cols: int) (es: #seq<int * int * float>) = ofSeqi rows cols es
 
     /// Create a square matrix with constant diagonal entries.
     let inline constDiag (n: int) (f: float) =

src/FSharp/LinearAlgebra.Double.Vector.fs

@@ -45,24 +45,20 @@ module Vector =
     /// In-place mutation by applying a function to every element of the vector.
     let inline mapInPlace (f: float -> float) (v: #Vector<float>) =
         v.MapInplace((fun x -> f x), true)
-        ()
 
     /// In-place mutation by applying a function to every element of the vector.
     let inline mapiInPlace (f: int -> float -> float) (v: #Vector<float>) =
         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 mapnzInPlace (f: float -> float) (v: #Vector<float>) =
         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 mapinzInPlace (f: int -> float -> float) (v: #Vector<float>) =
         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<float>) =

src/FSharpUnitTests/DenseVectorTests.fs

@@ -4,6 +4,8 @@ open NUnit.Framework
 open FsUnit
 open MathNet.Numerics.LinearAlgebra.Generic
 open MathNet.Numerics.LinearAlgebra.Double
+open MathNet.Numerics.Distributions
+open MathNet.Numerics.Statistics
 
 /// Unit tests for the dense vector type.
 module DenseVectorTests =
@@ -14,8 +16,23 @@ module DenseVectorTests =
     /// A large vector with increasingly large entries
     let largev = new DenseVector( Array.init 100 (fun i -> float i / 100.0) )
 
+    [<Test>]
+    let ``DenseVector.zeroCreate`` () =
+        (DenseVector.zeroCreate 100) + largev |> should equal largev
+        
+    [<Test>]
+    let ``DenseVector.randomCreate`` () =
+        let m = DenseVector.randomCreate 100 (Normal.WithMeanStdDev(100.0,0.1))
+        m.Values |> ArrayStatistics.Mean |> should (equalWithin 10.0) 100.0
+        m.Count |> should equal 100
+
+    [<Test>]
+    let ``DenseVector.create`` () =
+        DenseVector.create 5 0.3 |> should equal smallv
+
     [<Test>]
     let ``DenseVector.init`` () =
+        DenseVector.init 5 (fun i -> 0.3) |> should equal smallv
         DenseVector.init 100 (fun i -> float i / 100.0) |> should equal largev
 
     [<Test>]

src/FSharpUnitTests/MatrixTests.fs

@@ -12,6 +12,9 @@ module MatrixTests =
     let smallM = DenseMatrix.OfArray( Array2D.create 2 2 0.3 )
     let failingFoldBackM = DenseMatrix.init 2 3 (fun i j -> 1.0)
 
+    /// A small sparse matrix.
+    let sparseM = SparseMatrix.ofListi 2 3 [(1,0,0.3)]
+
     /// A large vector with increasingly large entries
     let largeM = DenseMatrix.OfArray( Array2D.init 100 100 (fun i j -> float i * 100.0 + float j) )
 
@@ -41,6 +44,70 @@ module MatrixTests =
             m.[..1,2..] <- DenseMatrix(2,2,[|5.;6.;7.;8.|]);
             m |> should equal (DenseMatrix(4,4,[|0.;100.;200.;300.;1.;101.;201.;301.;5.;6.;202.;302.;7.;8.;203.;303.|]))
 
+    [<Test>]
+    let ``Matrix.toArray2`` () =
+        Matrix.toArray2 smallM |> should array2_equal (Array2D.create 2 2 0.3)
+
+    [<Test>]
+    let ``Matrix.mapInPlace.Dense`` () =
+        let M = largeM.Clone()
+        M |> Matrix.mapInPlace (fun x -> 3.0 * x)
+        M |> should equal (3.0 * largeM)
+
+    [<Test>]
+    let ``Matrix.mapInPlace.Sparse`` () =
+        let M = sparseM.Clone()
+        M |> Matrix.mapInPlace (fun x -> 3.0 * x)
+        M |> should equal (3.0 * sparseM)
+
+    [<Test>]
+    let ``Matrix.mapnzInPlace.Sparse`` () =
+        let M = sparseM.Clone()
+        M |> Matrix.mapnzInPlace (fun x -> 3.0 * x)
+        M |> should equal (3.0 * sparseM)
+
+    [<Test>]
+    let ``Matrix.mapiInPlace.Dense`` () =
+        let M = largeM.Clone()
+        M |> Matrix.mapiInPlace (fun i j x -> 2.0 * (float i * 100.0 + float j) + x)
+        M |> should equal (3.0 * largeM)
+
+    [<Test>]
+    let ``Matrix.mapiInPlace.Sparse`` () =
+        let M = sparseM.Clone()
+        M |> Matrix.mapiInPlace (fun i j x -> if i=j then 2.0*x+1.0 else 2.0*x)
+        M |> should equal (2.0 * sparseM + SparseMatrix.init 2 3 (fun i j -> if i=j then 1.0 else 0.0))
+
+    [<Test>]
+    let ``Matrix.mapinzInPlace.Sparse`` () =
+        let M = sparseM.Clone()
+        M |> Matrix.mapinzInPlace (fun i j x -> 2.0*x)
+        M |> should equal (2.0 * sparseM)
+
+    [<Test>]
+    let ``Matrix.map`` () =
+        Matrix.map (fun x -> 2.0 * x) smallM |> should equal (2.0 * smallM)
+
+    [<Test>]
+    let ``Matrix.mapnz`` () =
+        Matrix.mapnz (fun x -> 2.0 * x) smallM |> should equal (2.0 * smallM)
+
+    [<Test>]
+    let ``Matrix.mapi`` () =
+        Matrix.mapi (fun i j x -> float i * 100.0 + float j + x) largeM |> should equal (2.0 * largeM)
+
+    [<Test>]
+    let ``Matrix.mapinz`` () =
+        Matrix.mapinz (fun i j x -> float i * 100.0 + float j + x) largeM |> should equal (2.0 * largeM)
+
+    [<Test>]
+    let ``Matrix.mapCols`` () =
+        Matrix.mapCols (fun j col -> col.Add(float j)) smallM |> should (approximately_matrix_equal 14) (matrix [[0.3;1.3];[0.3;1.3]])
+
+    [<Test>]
+    let ``Matrix.mapRows`` () =
+        Matrix.mapRows (fun i row -> row.Add(float i)) smallM |> should (approximately_matrix_equal 14) (matrix [[0.3;0.3];[1.3;1.3]])
+
     [<Test>]
     let ``Matrix.fold`` () =
         Matrix.fold (fun a b -> a - b) 0.0 smallM |> should equal -1.2
@@ -57,10 +124,6 @@ module MatrixTests =
     let ``Matrix.foldi`` () =
         Matrix.foldi (fun i j acc x -> acc + x + float (i+j)) 0.0 smallM |> should equal 5.2
 
-    [<Test>]
-    let ``Matrix.toArray2`` () =
-        Matrix.toArray2 smallM |> should array2_equal (Array2D.create 2 2 0.3)
-
     [<Test>]
     let ``Matrix.forall`` () =
         Matrix.forall (fun x -> x = 0.3) smallM |> should equal true
@@ -77,34 +140,12 @@ module MatrixTests =
     let ``Matrix.existsi`` () =
         Matrix.existsi (fun i j x -> x = float i * 100.0 + float j) largeM |> should equal true
 
-    [<Test>]
-    let ``Matrix.map`` () =
-        Matrix.map (fun x -> 2.0 * x) smallM |> should equal (2.0 * smallM)
-
-    [<Test>]
-    let ``Matrix.mapi`` () =
-        Matrix.mapi (fun i j x -> float i * 100.0 + float j + x) largeM |> should equal (2.0 * largeM)
-
-    [<Test>]
-    let ``Matrix.mapCols`` () =
-        Matrix.mapCols (fun j col -> col.Add(float j)) smallM |> should (approximately_matrix_equal 14) (matrix [[0.3;1.3];[0.3;1.3]])
-
-    [<Test>]
-    let ``Matrix.mapRows`` () =
-        Matrix.mapRows (fun i row -> row.Add(float i)) smallM |> should (approximately_matrix_equal 14) (matrix [[0.3;0.3];[1.3;1.3]])
-
     [<Test>]
     let ``Matrix.inplaceAssign`` () =
         let N = smallM.Clone()
-        Matrix.inplaceAssign (fun i j -> 0.0) N
+        N |> Matrix.inplaceAssign (fun i j -> 0.0)
         N |> should equal (0.0 * smallM)
 
-    [<Test>]
-    let ``Matrix.inplaceMapi`` () =
-        let N = largeM.Clone()
-        Matrix.inplaceMapi (fun i j x -> 2.0 * (float i * 100.0 + float j) + x) N
-        N |> should equal (3.0 * largeM)
-
     [<Test>]
     let ``Matrix.nonZeroEntries`` () =
         Seq.length (Matrix.nonZeroEntries smallM) |> should equal 4

src/FSharpUnitTests/VectorTests.fs

@@ -9,13 +9,13 @@ open MathNet.Numerics.LinearAlgebra.Double
 module VectorTests =
 
     /// A small uniform vector.
-    let smallv = new DenseVector([|0.3;0.3;0.3;0.3;0.3|]) :> Vector<float>
+    let smallv = DenseVector([|0.3;0.3;0.3;0.3;0.3|])
 
-    /// A small uniform vector.
-    let sparsev = SparseVector.OfIndexedEnumerable(5, [(1,0.3)]) :> Vector<float>
+    /// A small sparse vector.
+    let sparsev = SparseVector.ofListi 5 [(1,0.3)]
 
     /// A large vector with increasingly large entries
-    let largev = new DenseVector(Array.init 100 (fun i -> float i / 100.0)) :> Vector<float>
+    let largev = DenseVector(Array.init 100 (fun i -> float i / 100.0))
 
     [<Test>]
     let ``Vector.GetSlice`` () =
@@ -54,37 +54,37 @@ module VectorTests =
     [<Test>]
     let ``Vector.mapInPlace.Dense`` () =
         let w = smallv.Clone()
-        Vector.mapInPlace (fun x -> 2.0 * x) w
+        w |> Vector.mapInPlace (fun x -> 2.0 * x)
         w |> should equal (2.0 * smallv)
 
     [<Test>]
     let ``Vector.mapInPlace.Sparse`` () =
         let w = sparsev.Clone()
-        Vector.mapInPlace (fun x -> 2.0 * x) w
+        w |> Vector.mapInPlace (fun x -> 2.0 * x)
         w |> should equal (2.0 * sparsev)
 
     [<Test>]
     let ``Vector.mapnzInPlace.Sparse`` () =
         let w = sparsev.Clone()
-        Vector.mapnzInPlace (fun x -> 2.0 * x) w
+        w |> Vector.mapnzInPlace (fun x -> 2.0 * x)
         w |> should equal (2.0 * sparsev)
 
     [<Test>]
     let ``Vector.mapiInPlace.Dense`` () =
         let w = largev.Clone()
-        Vector.mapiInPlace (fun i x -> float i / 100.0) w
+        w |> Vector.mapiInPlace (fun i x -> float i / 100.0)
         w |> should equal (largev)
-        
+
     [<Test>]
     let ``Vector.mapiInPlace.Sparse`` () =
         let w = sparsev.Clone()
-        Vector.mapiInPlace (fun i x -> 2.0 * float i * x) w
+        w |> Vector.mapiInPlace (fun i x -> 2.0 * float i * x)
         w |> should equal (2.0 * sparsev)
 
     [<Test>]
     let ``Vector.mapinzInPlace.Sparse`` () =
         let w = sparsev.Clone()
-        Vector.mapinzInPlace (fun i x -> 2.0 * float i * x) w
+        w |> Vector.mapinzInPlace (fun i x -> 2.0 * float i * x)
         w |> should equal (2.0 * sparsev)
 
     [<Test>]

src/Numerics/LinearAlgebra/Generic/Matrix.cs

@@ -1438,5 +1438,26 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         {
             return Storage.ToRowMajorArray();
         }
+
+        /// <summary>
+        /// Applies a function to each value of this matrix and replaces the value with its result.
+        /// 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).
+        /// </summary>
+        public void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        {
+            Storage.MapInplace(f, forceMapZeros);
+        }
+
+        /// <summary>
+        /// Applies a function to each value of this matrix and replaces the value with its result.
+        /// The row and column indices of each value (zero-based) are passed as first arguments 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).
+        /// </summary>
+        public void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        {
+            Storage.MapIndexedInplace(f, forceMapZeros);
+        }
     }
 }

src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@@ -364,5 +364,28 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
             return ret;
         }
+
+        // FUNCTIONAL COMBINATORS
+
+        public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        {
+            for (int i = 0; i < Data.Length; i++)
+            {
+                Data[i] = f(Data[i]);
+            }
+        }
+
+        public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        {
+            int index = 0;
+            for (int j = 0; j < ColumnCount; j++)
+            {
+                for (int i = 0; i < RowCount; i++)
+                {
+                    Data[index] = f(i, j, Data[index]);
+                    index++;
+                }
+            }
+        }
     }
 }

src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs

@@ -525,5 +525,27 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
             return ret;
         }
+
+        // FUNCTIONAL COMBINATORS
+
+        public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        {
+            // we deliberately ignore forceMapZeros since we would not actually
+            // support any non-zero results outside of the diagonal anyway
+            for (int i = 0; i < Data.Length; i++)
+            {
+                Data[i] = f(Data[i]);
+            }
+        }
+
+        public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        {
+            // we deliberately ignore forceMapZeros since we would not actually
+            // support any non-zero results outside of the diagonal anyway
+            for (int i = 0; i < Data.Length; i++)
+            {
+                Data[i] = f(i, i, Data[i]);
+            }
+        }
     }
 }

src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs

@@ -411,5 +411,29 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
             return ret;
         }
+
+        // FUNCTIONAL COMBINATORS
+
+        public virtual void MapInplace(Func<T, T> 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<int, int, T, T> 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)));
+                }
+            }
+        }
     }
 }

src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@@ -456,6 +456,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 rowPointers[row] = index;
                 if (trows[row] != null)
                 {
+                    // TODO: Don't we need to sort here!?
                     foreach (var item in trows[row])
                     {
                         values.Add(item.Item2);
@@ -907,5 +908,107 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
             return ret;
         }
+
+        // FUNCTIONAL COMBINATORS
+
+        public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        {
+            var newRowPointers = new int[RowCount];
+            var newColumnIndices = new List<int>();
+            var newValues = new List<T>();
+
+            if (forceMapZeros || !Zero.Equals(f(Zero)))
+            {
+                int k = 0;
+                for (int row = 0; row < RowCount; row++)
+                {
+                    newRowPointers[row] = newValues.Count;
+                    for (int col = 0; col < ColumnCount; col++)
+                    {
+                        var item = k < (row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount) && (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 = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
+                    for (var j = startIndex; j < endIndex; j++)
+                    {
+                        var item = f(Values[j]);
+                        if (!Zero.Equals(item))
+                        {
+                            newValues.Add(item);
+                            newColumnIndices.Add(ColumnIndices[j]);
+                        }
+                    }
+                }
+            }
+
+            ColumnIndices = newColumnIndices.ToArray();
+            Values = newValues.ToArray();
+            ValueCount = newValues.Count;
+            Array.Copy(newRowPointers, RowPointers, RowCount);
+        }
+
+        public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        {
+            var newRowPointers = new int[RowCount];
+            var newColumnIndices = new List<int>();
+            var newValues = new List<T>();
+
+            if (forceMapZeros || !Zero.Equals(f(0,0,Zero)))
+            {
+                int k = 0;
+                for (int row = 0; row < RowCount; row++)
+                {
+                    newRowPointers[row] = newValues.Count;
+                    for (int col = 0; col < ColumnCount; col++)
+                    {
+                        var item = k < (row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount) && (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 = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
+                    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]);
+                        }
+                    }
+                }
+            }
+
+            ColumnIndices = newColumnIndices.ToArray();
+            Values = newValues.ToArray();
+            ValueCount = newValues.Count;
+            Array.Copy(newRowPointers, RowPointers, RowCount);
+        }
     }
 }

src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs

@@ -612,7 +612,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
             else
             {
-                for (int i = 0; i < Values.Length; i++)
+                for (int i = 0; i < ValueCount; i++)
                 {
                     var item = f(Values[i]);
                     if (!Zero.Equals(item))
@@ -646,7 +646,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
             else
             {
-                for (int i = 0; i < Values.Length; i++)
+                for (int i = 0; i < ValueCount; i++)
                 {
                     var item = f(Indices[i], Values[i]);
                     if (!Zero.Equals(item))