From fbfe077552ad4eb39ac9343f36a9320532355174 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 21 Mar 2015 19:57:22 +0100 Subject: [PATCH] LA: migrate mixed-storage fallback implementations to use higher order functions --- src/Numerics/LinearAlgebra/Complex/Matrix.cs | 486 +++++++--------- src/Numerics/LinearAlgebra/Complex/Vector.cs | 85 +-- .../LinearAlgebra/Complex32/Matrix.cs | 484 +++++++--------- .../LinearAlgebra/Complex32/Vector.cs | 85 +-- src/Numerics/LinearAlgebra/Double/Matrix.cs | 526 ++++++++---------- src/Numerics/LinearAlgebra/Double/Vector.cs | 122 ++-- src/Numerics/LinearAlgebra/Matrix.cs | 18 + src/Numerics/LinearAlgebra/Single/Matrix.cs | 525 ++++++++--------- src/Numerics/LinearAlgebra/Single/Vector.cs | 122 ++-- src/Numerics/LinearAlgebra/Vector.cs | 4 + 10 files changed, 1022 insertions(+), 1435 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs index 8a8594ea..0b111ca4 100644 --- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -65,201 +65,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex MapInplace(x => x.Magnitude < threshold ? Complex.Zero : x, Zeros.AllowSkip); } - /// Calculates the induced L1 norm of this matrix. - /// The maximum absolute column sum of the matrix. - public override double L1Norm() - { - var norm = 0d; - for (var j = 0; j < ColumnCount; j++) - { - var s = 0d; - for (var i = 0; i < RowCount; i++) - { - s += At(i, j).Magnitude; - } - norm = Math.Max(norm, s); - } - return norm; - } - - /// Calculates the induced infinity norm of this matrix. - /// The maximum absolute row sum of the matrix. - public override double InfinityNorm() - { - var norm = 0d; - for (var i = 0; i < RowCount; i++) - { - var s = 0d; - for (var j = 0; j < ColumnCount; j++) - { - s += At(i, j).Magnitude; - } - norm = Math.Max(norm, s); - } - return norm; - } - - /// Calculates the entry-wise Frobenius norm of this matrix. - /// The square root of the sum of the squared values. - public override double FrobeniusNorm() - { - var transpose = ConjugateTranspose(); - var aat = this*transpose; - var norm = 0d; - for (var i = 0; i < RowCount; i++) - { - norm += aat.At(i, i).Magnitude; - } - return Math.Sqrt(norm); - } - /// - /// Calculates the p-norms of all row vectors. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override Vector RowNorms(double norm) - { - if (norm <= 0.0) - { - throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); - } - - var ret = new double[RowCount]; - if (norm == 2.0) - { - Storage.FoldByRowUnchecked(ret, (s, x) => s + x.MagnitudeSquared(), (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); - } - else if (norm == 1.0) - { - Storage.FoldByRowUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); - } - else if (double.IsPositiveInfinity(norm)) - { - Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); - } - else - { - double invnorm = 1.0/norm; - Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); - } - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the p-norms of all column vectors. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override Vector ColumnNorms(double norm) - { - if (norm <= 0.0) - { - throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); - } - - var ret = new double[ColumnCount]; - if (norm == 2.0) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.MagnitudeSquared(), (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); - } - else if (norm == 1.0) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); - } - else if (double.IsPositiveInfinity(norm)) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); - } - else - { - double invnorm = 1.0/norm; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); - } - return Vector.Build.Dense(ret); - } - - /// - /// Normalizes all row vectors to a unit p-norm. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override sealed Matrix NormalizeRows(double norm) - { - var norminv = ((DenseVectorStorage)RowNorms(norm).Storage).Data; - for (int i = 0; i < norminv.Length; i++) - { - norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; - } - - var result = Build.SameAs(this, RowCount, ColumnCount); - Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[i]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); - return result; - } - - /// - /// Normalizes all column vectors to a unit p-norm. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override sealed Matrix NormalizeColumns(double norm) - { - var norminv = ((DenseVectorStorage)ColumnNorms(norm).Storage).Data; - for (int i = 0; i < norminv.Length; i++) - { - norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; - } - - var result = Build.SameAs(this, RowCount, ColumnCount); - Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[j]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); - return result; - } - - /// - /// Calculates the value sum of each row vector. - /// - public override Vector RowSums() - { - var ret = new Complex[RowCount]; - Storage.FoldByRowUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the absolute value sum of each row vector. - /// - public override Vector RowAbsoluteSums() - { - var ret = new Complex[RowCount]; - Storage.FoldByRowUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the value sum of each column vector. + /// Returns the conjugate transpose of this matrix. /// - public override Vector ColumnSums() + /// The conjugate transpose of this matrix. + public override sealed Matrix ConjugateTranspose() { - var ret = new Complex[ColumnCount]; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); + var ret = Transpose(); + ret.MapInplace(c => c.Conjugate(), Zeros.AllowSkip); + return ret; } /// - /// Calculates the absolute value sum of each column vector. + /// Complex conjugates each element of this matrix and place the results into the result matrix. /// - public override Vector ColumnAbsoluteSums() + /// The result of the conjugation. + protected override void DoConjugate(Matrix result) { - var ret = new Complex[ColumnCount]; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); + Map(Complex.Conjugate, result, Zeros.AllowSkip); } /// - /// Returns the conjugate transpose of this matrix. + /// Negate each element of this matrix and place the results into the result matrix. /// - /// The conjugate transpose of this matrix. - public override sealed Matrix ConjugateTranspose() + /// The result of the negation. + protected override void DoNegate(Matrix result) { - var ret = Transpose(); - ret.MapInplace(c => c.Conjugate(), Zeros.AllowSkip); - return ret; + Map(Complex.Negate, result, Zeros.AllowSkip); } /// @@ -269,13 +101,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The matrix to store the result of the addition. protected override void DoAdd(Complex scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) + scalar); - } - } + Map(x => x + scalar, result, Zeros.Include); } /// @@ -287,13 +113,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// If the two matrices don't have the same dimensions. protected override void DoAdd(Matrix other, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) + other.At(i, j)); - } - } + Map2(Complex.Add, other, result, Zeros.AllowSkip); } /// @@ -303,13 +123,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The matrix to store the result of the subtraction. protected override void DoSubtract(Complex scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) - scalar); - } - } + Map(x => x - scalar, result, Zeros.Include); } /// @@ -321,13 +135,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// If the two matrices don't have the same dimensions. protected override void DoSubtract(Matrix other, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) - other.At(i, j)); - } - } + Map2(Complex.Subtract, other, result, Zeros.AllowSkip); } /// @@ -337,13 +145,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The matrix to store the result of the multiplication. protected override void DoMultiply(Complex scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j)*scalar); - } - } + Map(x => x*scalar, result, Zeros.AllowSkip); } /// @@ -392,23 +194,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The matrix to store the result of the division. protected override void DoDivide(Complex divisor, Matrix result) { - DoMultiply(1.0/divisor, result); + Map(x => x/divisor, result, divisor.IsZero() ? Zeros.Include : Zeros.AllowSkip); } /// /// Divides a scalar by each element of the matrix and stores the result in the result matrix. /// - /// The scalar to add. + /// The scalar to divide by each element of the matrix. /// The matrix to store the result of the division. protected override void DoDivideByThis(Complex dividend, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, dividend/At(i, j)); - } - } + Map(x => dividend/x, result, Zeros.Include); } /// @@ -531,36 +327,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } } - /// - /// Negate each element of this matrix and place the results into the result matrix. - /// - /// The result of the negation. - protected override void DoNegate(Matrix result) - { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, -At(i, j)); - } - } - } - - /// - /// Complex conjugates each element of this matrix and place the results into the result matrix. - /// - /// The result of the conjugation. - protected override void DoConjugate(Matrix result) - { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j).Conjugate()); - } - } - } - /// /// Pointwise multiplies this matrix with another matrix and stores the result into the result matrix. /// @@ -568,13 +334,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The matrix to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Matrix other, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, At(i, j)*other.At(i, j)); - } - } + Map2(Complex.Multiply, other, result, Zeros.AllowSkip); } /// @@ -584,13 +344,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The matrix to store the result of the pointwise division. protected override void DoPointwiseDivide(Matrix divisor, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, At(i, j)/divisor.At(i, j)); - } - } + Map2(Complex.Divide, divisor, result, Zeros.Include); } /// @@ -600,7 +354,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Complex exponent, Matrix result) { - Map(x => x.Power(exponent), result, Zeros.AllowSkip); + Map(x => x.Power(exponent), result, Zeros.Include); } /// @@ -687,6 +441,194 @@ namespace MathNet.Numerics.LinearAlgebra.Complex Map(Complex.Log, result, Zeros.Include); } + + + /// Calculates the induced L1 norm of this matrix. + /// The maximum absolute column sum of the matrix. + public override double L1Norm() + { + var norm = 0d; + for (var j = 0; j < ColumnCount; j++) + { + var s = 0d; + for (var i = 0; i < RowCount; i++) + { + s += At(i, j).Magnitude; + } + norm = Math.Max(norm, s); + } + return norm; + } + + /// Calculates the induced infinity norm of this matrix. + /// The maximum absolute row sum of the matrix. + public override double InfinityNorm() + { + var norm = 0d; + for (var i = 0; i < RowCount; i++) + { + var s = 0d; + for (var j = 0; j < ColumnCount; j++) + { + s += At(i, j).Magnitude; + } + norm = Math.Max(norm, s); + } + return norm; + } + + /// Calculates the entry-wise Frobenius norm of this matrix. + /// The square root of the sum of the squared values. + public override double FrobeniusNorm() + { + var transpose = ConjugateTranspose(); + var aat = this*transpose; + var norm = 0d; + for (var i = 0; i < RowCount; i++) + { + norm += aat.At(i, i).Magnitude; + } + return Math.Sqrt(norm); + } + + /// + /// Calculates the p-norms of all row vectors. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override Vector RowNorms(double norm) + { + if (norm <= 0.0) + { + throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); + } + + var ret = new double[RowCount]; + if (norm == 2.0) + { + Storage.FoldByRowUnchecked(ret, (s, x) => s + x.MagnitudeSquared(), (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); + } + else if (norm == 1.0) + { + Storage.FoldByRowUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); + } + else if (double.IsPositiveInfinity(norm)) + { + Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); + } + else + { + double invnorm = 1.0/norm; + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); + } + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the p-norms of all column vectors. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override Vector ColumnNorms(double norm) + { + if (norm <= 0.0) + { + throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); + } + + var ret = new double[ColumnCount]; + if (norm == 2.0) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.MagnitudeSquared(), (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); + } + else if (norm == 1.0) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); + } + else if (double.IsPositiveInfinity(norm)) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); + } + else + { + double invnorm = 1.0/norm; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); + } + return Vector.Build.Dense(ret); + } + + /// + /// Normalizes all row vectors to a unit p-norm. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override sealed Matrix NormalizeRows(double norm) + { + var norminv = ((DenseVectorStorage)RowNorms(norm).Storage).Data; + for (int i = 0; i < norminv.Length; i++) + { + norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; + } + + var result = Build.SameAs(this, RowCount, ColumnCount); + Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[i]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); + return result; + } + + /// + /// Normalizes all column vectors to a unit p-norm. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override sealed Matrix NormalizeColumns(double norm) + { + var norminv = ((DenseVectorStorage)ColumnNorms(norm).Storage).Data; + for (int i = 0; i < norminv.Length; i++) + { + norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; + } + + var result = Build.SameAs(this, RowCount, ColumnCount); + Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[j]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); + return result; + } + + /// + /// Calculates the value sum of each row vector. + /// + public override Vector RowSums() + { + var ret = new Complex[RowCount]; + Storage.FoldByRowUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the absolute value sum of each row vector. + /// + public override Vector RowAbsoluteSums() + { + var ret = new Complex[RowCount]; + Storage.FoldByRowUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the value sum of each column vector. + /// + public override Vector ColumnSums() + { + var ret = new Complex[ColumnCount]; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the absolute value sum of each column vector. + /// + public override Vector ColumnAbsoluteSums() + { + var ret = new Complex[ColumnCount]; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + /// /// Computes the trace of this matrix. /// diff --git a/src/Numerics/LinearAlgebra/Complex/Vector.cs b/src/Numerics/LinearAlgebra/Complex/Vector.cs index 6cda7da6..a057e15f 100644 --- a/src/Numerics/LinearAlgebra/Complex/Vector.cs +++ b/src/Numerics/LinearAlgebra/Complex/Vector.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -63,6 +63,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex MapInplace(x => x.Magnitude < threshold ? Complex.Zero : x, Zeros.AllowSkip); } + /// + /// Conjugates vector and save result to + /// + /// Target vector + protected override void DoConjugate(Vector result) + { + Map(Complex.Conjugate, result, Zeros.AllowSkip); + } + + /// + /// Negates vector and saves result to + /// + /// Target vector + protected override void DoNegate(Vector result) + { + Map(Complex.Negate, result, Zeros.AllowSkip); + } + /// /// Adds a scalar to each element of the vector and stores the result in the result vector. /// @@ -74,10 +92,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoAdd(Complex scalar, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) + scalar); - } + Map(x => x + scalar, result, Zeros.Include); } /// @@ -91,10 +106,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoAdd(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) + other.At(index)); - } + Map2(Complex.Add, other, result, Zeros.AllowSkip); } /// @@ -108,7 +120,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoSubtract(Complex scalar, Vector result) { - DoAdd(-scalar, result); + Map(x => x - scalar, result, Zeros.Include); } /// @@ -122,10 +134,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoSubtract(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) - other.At(index)); - } + Map2(Complex.Subtract, other, result, Zeros.AllowSkip); } /// @@ -139,10 +148,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoMultiply(Complex scalar, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) * scalar); - } + Map(x => x*scalar, result, Zeros.AllowSkip); } /// @@ -156,7 +162,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoDivide(Complex divisor, Vector result) { - DoMultiply(1 / divisor, result); + Map(x => x/divisor, result, divisor.IsZero() ? Zeros.Include : Zeros.AllowSkip); } /// @@ -166,10 +172,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The vector to store the result of the division. protected override void DoDivideByThis(Complex dividend, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, dividend / At(index)); - } + Map(x => dividend/x, result, Zeros.Include); } /// @@ -179,10 +182,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The vector to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) * other.At(index)); - } + Map2(Complex.Multiply, other, result, Zeros.AllowSkip); } /// @@ -192,10 +192,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The vector to store the result of the pointwise division. protected override void DoPointwiseDivide(Vector divisor, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) / divisor.At(index)); - } + Map2(Complex.Divide, divisor, result, Zeros.Include); } /// @@ -205,7 +202,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Complex exponent, Vector result) { - Map(x => x.Power(exponent), result, Zeros.AllowSkip); + Map(x => x.Power(exponent), result, Zeros.Include); } /// @@ -453,30 +450,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return Math.Pow(sum, 1.0/p); } - /// - /// Conjugates vector and save result to - /// - /// Target vector - protected override void DoConjugate(Vector result) - { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index).Conjugate()); - } - } - - /// - /// Negates vector and saves result to - /// - /// Target vector - protected override void DoNegate(Vector result) - { - for (var index = 0; index < Count; index++) - { - result.At(index, -At(index)); - } - } - /// /// Returns the index of the absolute maximum element. /// diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs index cc9d5c95..a4501582 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -60,201 +60,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 MapInplace(x => x.Magnitude < threshold ? Complex32.Zero : x, Zeros.AllowSkip); } - /// Calculates the induced L1 norm of this matrix. - /// The maximum absolute column sum of the matrix. - public override double L1Norm() - { - var norm = 0d; - for (var j = 0; j < ColumnCount; j++) - { - var s = 0d; - for (var i = 0; i < RowCount; i++) - { - s += At(i, j).Magnitude; - } - norm = Math.Max(norm, s); - } - return norm; - } - - /// Calculates the induced infinity norm of this matrix. - /// The maximum absolute row sum of the matrix. - public override double InfinityNorm() - { - var norm = 0d; - for (var i = 0; i < RowCount; i++) - { - var s = 0d; - for (var j = 0; j < ColumnCount; j++) - { - s += At(i, j).Magnitude; - } - norm = Math.Max(norm, s); - } - return norm; - } - - /// Calculates the entry-wise Frobenius norm of this matrix. - /// The square root of the sum of the squared values. - public override double FrobeniusNorm() - { - var transpose = ConjugateTranspose(); - var aat = this*transpose; - var norm = 0d; - for (var i = 0; i < RowCount; i++) - { - norm += aat.At(i, i).Magnitude; - } - return Math.Sqrt(norm); - } - - /// - /// Calculates the p-norms of all row vectors. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override Vector RowNorms(double norm) - { - if (norm <= 0.0) - { - throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); - } - - var ret = new double[RowCount]; - if (norm == 2.0) - { - Storage.FoldByRowUnchecked(ret, (s, x) => s + x.MagnitudeSquared, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); - } - else if (norm == 1.0) - { - Storage.FoldByRowUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); - } - else if (double.IsPositiveInfinity(norm)) - { - Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); - } - else - { - double invnorm = 1.0/norm; - Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); - } - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the p-norms of all column vectors. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override Vector ColumnNorms(double norm) - { - if (norm <= 0.0) - { - throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); - } - - var ret = new double[ColumnCount]; - if (norm == 2.0) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.MagnitudeSquared, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); - } - else if (norm == 1.0) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); - } - else if (double.IsPositiveInfinity(norm)) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); - } - else - { - double invnorm = 1.0/norm; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); - } - return Vector.Build.Dense(ret); - } - - /// - /// Normalizes all row vectors to a unit p-norm. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override sealed Matrix NormalizeRows(double norm) - { - var norminv = ((DenseVectorStorage)RowNorms(norm).Storage).Data; - for (int i = 0; i < norminv.Length; i++) - { - norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; - } - - var result = Build.SameAs(this, RowCount, ColumnCount); - Storage.MapIndexedTo(result.Storage, (i, j, x) => ((float)norminv[i])*x, Zeros.AllowSkip, ExistingData.AssumeZeros); - return result; - } - /// - /// Normalizes all column vectors to a unit p-norm. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override sealed Matrix NormalizeColumns(double norm) - { - var norminv = ((DenseVectorStorage)ColumnNorms(norm).Storage).Data; - for (int i = 0; i < norminv.Length; i++) - { - norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; - } - - var result = Build.SameAs(this, RowCount, ColumnCount); - Storage.MapIndexedTo(result.Storage, (i, j, x) => ((float)norminv[j])*x, Zeros.AllowSkip, ExistingData.AssumeZeros); - return result; - } - - /// - /// Calculates the value sum of each row vector. - /// - public override Vector RowSums() - { - var ret = new Complex32[RowCount]; - Storage.FoldByRowUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the absolute value sum of each row vector. - /// - public override Vector RowAbsoluteSums() - { - var ret = new Complex32[RowCount]; - Storage.FoldByRowUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the value sum of each column vector. + /// Returns the conjugate transpose of this matrix. /// - public override Vector ColumnSums() + /// The conjugate transpose of this matrix. + public override sealed Matrix ConjugateTranspose() { - var ret = new Complex32[ColumnCount]; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); + var ret = Transpose(); + ret.MapInplace(c => c.Conjugate(), Zeros.AllowSkip); + return ret; } /// - /// Calculates the absolute value sum of each column vector. + /// Complex conjugates each element of this matrix and place the results into the result matrix. /// - public override Vector ColumnAbsoluteSums() + /// The result of the conjugation. + protected override void DoConjugate(Matrix result) { - var ret = new Complex32[ColumnCount]; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); + Map(Complex32.Conjugate, result, Zeros.AllowSkip); } /// - /// Returns the conjugate transpose of this matrix. + /// Negate each element of this matrix and place the results into the result matrix. /// - /// The conjugate transpose of this matrix. - public override sealed Matrix ConjugateTranspose() + /// The result of the negation. + protected override void DoNegate(Matrix result) { - var ret = Transpose(); - ret.MapInplace(c => c.Conjugate(), Zeros.AllowSkip); - return ret; + Map(Complex32.Negate, result, Zeros.AllowSkip); } /// @@ -264,13 +96,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The matrix to store the result of the addition. protected override void DoAdd(Complex32 scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) + scalar); - } - } + Map(x => x + scalar, result, Zeros.Include); } /// @@ -282,13 +108,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// If the two matrices don't have the same dimensions. protected override void DoAdd(Matrix other, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) + other.At(i, j)); - } - } + Map2(Complex32.Add, other, result, Zeros.AllowSkip); } /// @@ -298,13 +118,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The matrix to store the result of the subtraction. protected override void DoSubtract(Complex32 scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) - scalar); - } - } + Map(x => x - scalar, result, Zeros.Include); } /// @@ -316,13 +130,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// If the two matrices don't have the same dimensions. protected override void DoSubtract(Matrix other, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) - other.At(i, j)); - } - } + Map2(Complex32.Subtract, other, result, Zeros.AllowSkip); } /// @@ -332,13 +140,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The matrix to store the result of the multiplication. protected override void DoMultiply(Complex32 scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j)*scalar); - } - } + Map(x => x*scalar, result, Zeros.AllowSkip); } /// @@ -366,23 +168,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The matrix to store the result of the division. protected override void DoDivide(Complex32 divisor, Matrix result) { - DoMultiply(1.0f/divisor, result); + Map(x => x/divisor, result, divisor.IsZero() ? Zeros.Include : Zeros.AllowSkip); } /// /// Divides a scalar by each element of the matrix and stores the result in the result matrix. /// - /// The scalar to add. + /// The scalar to divide by each element of the matrix. /// The matrix to store the result of the division. protected override void DoDivideByThis(Complex32 dividend, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, dividend/At(i, j)); - } - } + Map(x => dividend/x, result, Zeros.Include); } /// @@ -526,36 +322,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } } - /// - /// Negate each element of this matrix and place the results into the result matrix. - /// - /// The result of the negation. - protected override void DoNegate(Matrix result) - { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, -At(i, j)); - } - } - } - - /// - /// Complex conjugates each element of this matrix and place the results into the result matrix. - /// - /// The result of the conjugation. - protected override void DoConjugate(Matrix result) - { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j).Conjugate()); - } - } - } - /// /// Pointwise multiplies this matrix with another matrix and stores the result into the result matrix. /// @@ -563,13 +329,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The matrix to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Matrix other, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, At(i, j)*other.At(i, j)); - } - } + Map2(Complex32.Multiply, other, result, Zeros.AllowSkip); } /// @@ -579,13 +339,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The matrix to store the result of the pointwise division. protected override void DoPointwiseDivide(Matrix divisor, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, At(i, j)/divisor.At(i, j)); - } - } + Map2(Complex32.Divide, divisor, result, Zeros.Include); } /// @@ -595,7 +349,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Complex32 exponent, Matrix result) { - Map(x => x.Power(exponent), result, Zeros.AllowSkip); + Map(x => x.Power(exponent), result, Zeros.Include); } /// @@ -682,6 +436,192 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 Map(Complex32.Log, result, Zeros.Include); } + /// Calculates the induced L1 norm of this matrix. + /// The maximum absolute column sum of the matrix. + public override double L1Norm() + { + var norm = 0d; + for (var j = 0; j < ColumnCount; j++) + { + var s = 0d; + for (var i = 0; i < RowCount; i++) + { + s += At(i, j).Magnitude; + } + norm = Math.Max(norm, s); + } + return norm; + } + + /// Calculates the induced infinity norm of this matrix. + /// The maximum absolute row sum of the matrix. + public override double InfinityNorm() + { + var norm = 0d; + for (var i = 0; i < RowCount; i++) + { + var s = 0d; + for (var j = 0; j < ColumnCount; j++) + { + s += At(i, j).Magnitude; + } + norm = Math.Max(norm, s); + } + return norm; + } + + /// Calculates the entry-wise Frobenius norm of this matrix. + /// The square root of the sum of the squared values. + public override double FrobeniusNorm() + { + var transpose = ConjugateTranspose(); + var aat = this*transpose; + var norm = 0d; + for (var i = 0; i < RowCount; i++) + { + norm += aat.At(i, i).Magnitude; + } + return Math.Sqrt(norm); + } + + /// + /// Calculates the p-norms of all row vectors. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override Vector RowNorms(double norm) + { + if (norm <= 0.0) + { + throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); + } + + var ret = new double[RowCount]; + if (norm == 2.0) + { + Storage.FoldByRowUnchecked(ret, (s, x) => s + x.MagnitudeSquared, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); + } + else if (norm == 1.0) + { + Storage.FoldByRowUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); + } + else if (double.IsPositiveInfinity(norm)) + { + Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); + } + else + { + double invnorm = 1.0/norm; + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); + } + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the p-norms of all column vectors. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override Vector ColumnNorms(double norm) + { + if (norm <= 0.0) + { + throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); + } + + var ret = new double[ColumnCount]; + if (norm == 2.0) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.MagnitudeSquared, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); + } + else if (norm == 1.0) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); + } + else if (double.IsPositiveInfinity(norm)) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); + } + else + { + double invnorm = 1.0/norm; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); + } + return Vector.Build.Dense(ret); + } + + /// + /// Normalizes all row vectors to a unit p-norm. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override sealed Matrix NormalizeRows(double norm) + { + var norminv = ((DenseVectorStorage)RowNorms(norm).Storage).Data; + for (int i = 0; i < norminv.Length; i++) + { + norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; + } + + var result = Build.SameAs(this, RowCount, ColumnCount); + Storage.MapIndexedTo(result.Storage, (i, j, x) => ((float)norminv[i])*x, Zeros.AllowSkip, ExistingData.AssumeZeros); + return result; + } + + /// + /// Normalizes all column vectors to a unit p-norm. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override sealed Matrix NormalizeColumns(double norm) + { + var norminv = ((DenseVectorStorage)ColumnNorms(norm).Storage).Data; + for (int i = 0; i < norminv.Length; i++) + { + norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; + } + + var result = Build.SameAs(this, RowCount, ColumnCount); + Storage.MapIndexedTo(result.Storage, (i, j, x) => ((float)norminv[j])*x, Zeros.AllowSkip, ExistingData.AssumeZeros); + return result; + } + + /// + /// Calculates the value sum of each row vector. + /// + public override Vector RowSums() + { + var ret = new Complex32[RowCount]; + Storage.FoldByRowUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the absolute value sum of each row vector. + /// + public override Vector RowAbsoluteSums() + { + var ret = new Complex32[RowCount]; + Storage.FoldByRowUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the value sum of each column vector. + /// + public override Vector ColumnSums() + { + var ret = new Complex32[ColumnCount]; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the absolute value sum of each column vector. + /// + public override Vector ColumnAbsoluteSums() + { + var ret = new Complex32[ColumnCount]; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + /// /// Computes the trace of this matrix. /// diff --git a/src/Numerics/LinearAlgebra/Complex32/Vector.cs b/src/Numerics/LinearAlgebra/Complex32/Vector.cs index 865296d5..41d93489 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Vector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Vector.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -58,6 +58,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 MapInplace(x => x.Magnitude < threshold ? Complex32.Zero : x, Zeros.AllowSkip); } + /// + /// Conjugates vector and save result to + /// + /// Target vector + protected override void DoConjugate(Vector result) + { + Map(Complex32.Conjugate, result, Zeros.AllowSkip); + } + + /// + /// Negates vector and saves result to + /// + /// Target vector + protected override void DoNegate(Vector result) + { + Map(Complex32.Negate, result, Zeros.AllowSkip); + } + /// /// Adds a scalar to each element of the vector and stores the result in the result vector. /// @@ -69,10 +87,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoAdd(Complex32 scalar, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) + scalar); - } + Map(x => x + scalar, result, Zeros.Include); } /// @@ -86,10 +101,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoAdd(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) + other.At(index)); - } + Map2(Complex32.Add, other, result, Zeros.AllowSkip); } /// @@ -103,7 +115,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoSubtract(Complex32 scalar, Vector result) { - DoAdd(-scalar, result); + Map(x => x - scalar, result, Zeros.Include); } /// @@ -117,10 +129,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoSubtract(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) - other.At(index)); - } + Map2(Complex32.Subtract, other, result, Zeros.AllowSkip); } /// @@ -134,10 +143,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoMultiply(Complex32 scalar, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) * scalar); - } + Map(x => x*scalar, result, Zeros.AllowSkip); } /// @@ -151,7 +157,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoDivide(Complex32 divisor, Vector result) { - DoMultiply(1 / divisor, result); + Map(x => x/divisor, result, divisor.IsZero() ? Zeros.Include : Zeros.AllowSkip); } /// @@ -161,10 +167,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The vector to store the result of the division. protected override void DoDivideByThis(Complex32 dividend, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, dividend / At(index)); - } + Map(x => dividend/x, result, Zeros.Include); } /// @@ -174,10 +177,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The vector to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) * other.At(index)); - } + Map2(Complex32.Multiply, other, result, Zeros.AllowSkip); } /// @@ -187,10 +187,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The vector to store the result of the pointwise division. protected override void DoPointwiseDivide(Vector divisor, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) / divisor.At(index)); - } + Map2(Complex32.Divide, divisor, result, Zeros.Include); } /// @@ -200,7 +197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Complex32 exponent, Vector result) { - Map(x => x.Power(exponent), result, Zeros.AllowSkip); + Map(x => x.Power(exponent), result, Zeros.Include); } /// @@ -448,30 +445,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return Math.Pow(sum, 1.0/p); } - /// - /// Conjugates vector and save result to - /// - /// Target vector - protected override void DoConjugate(Vector result) - { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index).Conjugate()); - } - } - - /// - /// Negates vector and saves result to - /// - /// Target vector - protected override void DoNegate(Vector result) - { - for (var index = 0; index < Count; index++) - { - result.At(index, -At(index)); - } - } - /// /// Returns the index of the absolute maximum element. /// diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index 75dd70a0..edf8c603 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -58,199 +58,36 @@ namespace MathNet.Numerics.LinearAlgebra.Double MapInplace(x => Math.Abs(x) < threshold ? 0d : x, Zeros.AllowSkip); } - /// Calculates the induced L1 norm of this matrix. - /// The maximum absolute column sum of the matrix. - public override double L1Norm() - { - var norm = 0d; - for (var j = 0; j < ColumnCount; j++) - { - var s = 0d; - for (var i = 0; i < RowCount; i++) - { - s += Math.Abs(At(i, j)); - } - norm = Math.Max(norm, s); - } - return norm; - } - - /// Calculates the induced infinity norm of this matrix. - /// The maximum absolute row sum of the matrix. - public override double InfinityNorm() - { - var norm = 0d; - for (var i = 0; i < RowCount; i++) - { - var s = 0d; - for (var j = 0; j < ColumnCount; j++) - { - s += Math.Abs(At(i, j)); - } - norm = Math.Max(norm, s); - } - return norm; - } - - /// Calculates the entry-wise Frobenius norm of this matrix. - /// The square root of the sum of the squared values. - public override double FrobeniusNorm() - { - var transpose = Transpose(); - var aat = this*transpose; - var norm = 0d; - for (var i = 0; i < RowCount; i++) - { - norm += aat.At(i, i); - } - return Math.Sqrt(norm); - } - - /// - /// Calculates the p-norms of all row vectors. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override Vector RowNorms(double norm) - { - if (norm <= 0.0) - { - throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); - } - - var ret = new double[RowCount]; - if (norm == 2.0) - { - Storage.FoldByRowUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); - } - else if (norm == 1.0) - { - Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); - } - else if (double.IsPositiveInfinity(norm)) - { - Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); - } - else - { - double invnorm = 1.0/norm; - Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); - } - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the p-norms of all column vectors. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override Vector ColumnNorms(double norm) - { - if (norm <= 0.0) - { - throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); - } - - var ret = new double[ColumnCount]; - if (norm == 2.0) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); - } - else if (norm == 1.0) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); - } - else if (double.IsPositiveInfinity(norm)) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); - } - else - { - double invnorm = 1.0/norm; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); - } - return Vector.Build.Dense(ret); - } - /// - /// Normalizes all row vectors to a unit p-norm. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// Returns the conjugate transpose of this matrix. /// - public override sealed Matrix NormalizeRows(double norm) + /// The conjugate transpose of this matrix. + public override sealed Matrix ConjugateTranspose() { - var norminv = ((DenseVectorStorage)RowNorms(norm).Storage).Data; - for (int i = 0; i < norminv.Length; i++) - { - norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; - } - - var result = Build.SameAs(this, RowCount, ColumnCount); - Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[i]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); - return result; + return Transpose(); } /// - /// Normalizes all column vectors to a unit p-norm. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// Complex conjugates each element of this matrix and place the results into the result matrix. /// - public override sealed Matrix NormalizeColumns(double norm) + /// The result of the conjugation. + protected override sealed void DoConjugate(Matrix result) { - var norminv = ((DenseVectorStorage)ColumnNorms(norm).Storage).Data; - for (int i = 0; i < norminv.Length; i++) + if (ReferenceEquals(this, result)) { - norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; + return; } - var result = Build.SameAs(this, RowCount, ColumnCount); - Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[j]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); - return result; - } - - /// - /// Calculates the value sum of each row vector. - /// - public override Vector RowSums() - { - var ret = new double[RowCount]; - Storage.FoldByRowUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the absolute value sum of each row vector. - /// - public override Vector RowAbsoluteSums() - { - var ret = new double[RowCount]; - Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the value sum of each column vector. - /// - public override Vector ColumnSums() - { - var ret = new double[ColumnCount]; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the absolute value sum of each column vector. - /// - public override Vector ColumnAbsoluteSums() - { - var ret = new double[ColumnCount]; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); + CopyTo(result); } /// - /// Returns the conjugate transpose of this matrix. + /// Negate each element of this matrix and place the results into the result matrix. /// - /// The conjugate transpose of this matrix. - public override sealed Matrix ConjugateTranspose() + /// The result of the negation. + protected override void DoNegate(Matrix result) { - return Transpose(); + Map(x => -x, result, Zeros.AllowSkip); } /// @@ -260,13 +97,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The matrix to store the result of the addition. protected override void DoAdd(double scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) + scalar); - } - } + Map(x => x + scalar, result, Zeros.Include); } /// @@ -278,13 +109,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// If the two matrices don't have the same dimensions. protected override void DoAdd(Matrix other, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) + other.At(i, j)); - } - } + Map2((x, y) => x + y, other, result, Zeros.AllowSkip); } /// @@ -294,13 +119,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The matrix to store the result of the subtraction. protected override void DoSubtract(double scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) - scalar); - } - } + Map(x => x - scalar, result, Zeros.Include); } /// @@ -312,13 +131,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// If the two matrices don't have the same dimensions. protected override void DoSubtract(Matrix other, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) - other.At(i, j)); - } - } + Map2((x, y) => x - y, other, result, Zeros.AllowSkip); } /// @@ -328,13 +141,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The matrix to store the result of the multiplication. protected override void DoMultiply(double scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j)*scalar); - } - } + Map(x => x*scalar, result, Zeros.AllowSkip); } /// @@ -362,23 +169,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The matrix to store the result of the division. protected override void DoDivide(double divisor, Matrix result) { - DoMultiply(1.0/divisor, result); + Map(x => x/divisor, result, divisor == 0.0 ? Zeros.Include : Zeros.AllowSkip); } /// /// Divides a scalar by each element of the matrix and stores the result in the result matrix. /// - /// The scalar to add. + /// The scalar to divide by each element of the matrix. /// The matrix to store the result of the division. protected override void DoDivideByThis(double dividend, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, dividend/At(i, j)); - } - } + Map(x => dividend/x, result, Zeros.Include); } /// @@ -492,35 +293,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double DoTransposeThisAndMultiply(rightSide, result); } - /// - /// Negate each element of this matrix and place the results into the result matrix. - /// - /// The result of the negation. - protected override void DoNegate(Matrix result) - { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, -At(i, j)); - } - } - } - - /// - /// Complex conjugates each element of this matrix and place the results into the result matrix. - /// - /// The result of the conjugation. - protected override sealed void DoConjugate(Matrix result) - { - if (ReferenceEquals(this, result)) - { - return; - } - - CopyTo(result); - } - /// /// Computes the canonical modulus, where the result has the sign of the divisor, /// for the given divisor each element of the matrix. @@ -529,13 +301,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Matrix to store the results in. protected override void DoModulus(double divisor, Matrix result) { - for (var row = 0; row < RowCount; row++) - { - for (var column = 0; column < ColumnCount; column++) - { - result.At(row, column, Euclid.Modulus(At(row, column), divisor)); - } - } + Map(x => Euclid.Modulus(x, divisor), result, Zeros.Include); } /// @@ -546,13 +312,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoModulusByThis(double dividend, Matrix result) { - for (var row = 0; row < RowCount; row++) - { - for (var column = 0; column < ColumnCount; column++) - { - result.At(row, column, Euclid.Modulus(dividend, At(row, column))); - } - } + Map(x => Euclid.Modulus(dividend, x), result, Zeros.Include); } /// @@ -563,13 +323,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Matrix to store the results in. protected override void DoRemainder(double divisor, Matrix result) { - for (var row = 0; row < RowCount; row++) - { - for (var column = 0; column < ColumnCount; column++) - { - result.At(row, column, At(row, column)%divisor); - } - } + Map(x => Euclid.Remainder(x, divisor), result, Zeros.Include); } /// @@ -580,13 +334,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoRemainderByThis(double dividend, Matrix result) { - for (var row = 0; row < RowCount; row++) - { - for (var column = 0; column < ColumnCount; column++) - { - result.At(row, column, dividend%At(row, column)); - } - } + Map(x => Euclid.Remainder(dividend, x), result, Zeros.Include); } /// @@ -596,13 +344,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The matrix to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Matrix other, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, At(i, j)*other.At(i, j)); - } - } + Map2((x, y) => x*y, other, result, Zeros.AllowSkip); } /// @@ -612,13 +354,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The matrix to store the result of the pointwise division. protected override void DoPointwiseDivide(Matrix divisor, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, At(i, j)/divisor.At(i, j)); - } - } + Map2((x, y) => x/y, divisor, result, Zeros.Include); } /// @@ -628,7 +364,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(double exponent, Matrix result) { - Map(x => Math.Pow(x, exponent), result, Zeros.AllowSkip); + Map(x => Math.Pow(x, exponent), result, exponent > 0.0 ? Zeros.AllowSkip : Zeros.Include); } /// @@ -639,13 +375,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The result of the modulus. protected override void DoPointwiseModulus(Matrix divisor, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, Euclid.Modulus(At(i, j), divisor.At(i, j))); - } - } + Map2(Euclid.Modulus, divisor, result, Zeros.Include); } /// @@ -656,13 +386,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The result of the modulus. protected override void DoPointwiseRemainder(Matrix divisor, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, At(i, j)%divisor.At(i, j)); - } - } + Map2(Euclid.Remainder, divisor, result, Zeros.Include); } /// @@ -683,6 +407,192 @@ namespace MathNet.Numerics.LinearAlgebra.Double Map(Math.Log, result, Zeros.Include); } + /// Calculates the induced L1 norm of this matrix. + /// The maximum absolute column sum of the matrix. + public override double L1Norm() + { + var norm = 0d; + for (var j = 0; j < ColumnCount; j++) + { + var s = 0d; + for (var i = 0; i < RowCount; i++) + { + s += Math.Abs(At(i, j)); + } + norm = Math.Max(norm, s); + } + return norm; + } + + /// Calculates the induced infinity norm of this matrix. + /// The maximum absolute row sum of the matrix. + public override double InfinityNorm() + { + var norm = 0d; + for (var i = 0; i < RowCount; i++) + { + var s = 0d; + for (var j = 0; j < ColumnCount; j++) + { + s += Math.Abs(At(i, j)); + } + norm = Math.Max(norm, s); + } + return norm; + } + + /// Calculates the entry-wise Frobenius norm of this matrix. + /// The square root of the sum of the squared values. + public override double FrobeniusNorm() + { + var transpose = Transpose(); + var aat = this*transpose; + var norm = 0d; + for (var i = 0; i < RowCount; i++) + { + norm += aat.At(i, i); + } + return Math.Sqrt(norm); + } + + /// + /// Calculates the p-norms of all row vectors. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override Vector RowNorms(double norm) + { + if (norm <= 0.0) + { + throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); + } + + var ret = new double[RowCount]; + if (norm == 2.0) + { + Storage.FoldByRowUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); + } + else if (norm == 1.0) + { + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); + } + else if (double.IsPositiveInfinity(norm)) + { + Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); + } + else + { + double invnorm = 1.0/norm; + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); + } + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the p-norms of all column vectors. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override Vector ColumnNorms(double norm) + { + if (norm <= 0.0) + { + throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); + } + + var ret = new double[ColumnCount]; + if (norm == 2.0) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); + } + else if (norm == 1.0) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); + } + else if (double.IsPositiveInfinity(norm)) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); + } + else + { + double invnorm = 1.0/norm; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); + } + return Vector.Build.Dense(ret); + } + + /// + /// Normalizes all row vectors to a unit p-norm. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override sealed Matrix NormalizeRows(double norm) + { + var norminv = ((DenseVectorStorage)RowNorms(norm).Storage).Data; + for (int i = 0; i < norminv.Length; i++) + { + norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; + } + + var result = Build.SameAs(this, RowCount, ColumnCount); + Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[i]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); + return result; + } + + /// + /// Normalizes all column vectors to a unit p-norm. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override sealed Matrix NormalizeColumns(double norm) + { + var norminv = ((DenseVectorStorage)ColumnNorms(norm).Storage).Data; + for (int i = 0; i < norminv.Length; i++) + { + norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; + } + + var result = Build.SameAs(this, RowCount, ColumnCount); + Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[j]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); + return result; + } + + /// + /// Calculates the value sum of each row vector. + /// + public override Vector RowSums() + { + var ret = new double[RowCount]; + Storage.FoldByRowUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the absolute value sum of each row vector. + /// + public override Vector RowAbsoluteSums() + { + var ret = new double[RowCount]; + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the value sum of each column vector. + /// + public override Vector ColumnSums() + { + var ret = new double[ColumnCount]; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the absolute value sum of each column vector. + /// + public override Vector ColumnAbsoluteSums() + { + var ret = new double[ColumnCount]; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + /// /// Computes the trace of this matrix. /// diff --git a/src/Numerics/LinearAlgebra/Double/Vector.cs b/src/Numerics/LinearAlgebra/Double/Vector.cs index 01e80270..d72d12f8 100644 --- a/src/Numerics/LinearAlgebra/Double/Vector.cs +++ b/src/Numerics/LinearAlgebra/Double/Vector.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -56,6 +56,29 @@ namespace MathNet.Numerics.LinearAlgebra.Double MapInplace(x => Math.Abs(x) < threshold ? 0d : x, Zeros.AllowSkip); } + /// + /// Conjugates vector and save result to + /// + /// Target vector + protected override sealed void DoConjugate(Vector result) + { + if (ReferenceEquals(this, result)) + { + return; + } + + CopyTo(result); + } + + /// + /// Negates vector and saves result to + /// + /// Target vector + protected override void DoNegate(Vector result) + { + Map(x => -x, result, Zeros.AllowSkip); + } + /// /// Adds a scalar to each element of the vector and stores the result in the result vector. /// @@ -67,10 +90,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoAdd(double scalar, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) + scalar); - } + Map(x => x + scalar, result, Zeros.Include); } /// @@ -84,10 +104,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoAdd(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) + other.At(index)); - } + Map2((x, y) => x + y, other, result, Zeros.AllowSkip); } /// @@ -101,7 +118,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoSubtract(double scalar, Vector result) { - DoAdd(-scalar, result); + Map(x => x - scalar, result, Zeros.Include); } /// @@ -115,10 +132,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoSubtract(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) - other.At(index)); - } + Map2((x, y) => x - y, other, result, Zeros.AllowSkip); } /// @@ -132,10 +146,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoMultiply(double scalar, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) * scalar); - } + Map(x => x*scalar, result, Zeros.AllowSkip); } /// @@ -149,7 +160,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoDivide(double divisor, Vector result) { - DoMultiply(1 / divisor, result); + Map(x => x/divisor, result, divisor == 0.0 ? Zeros.Include : Zeros.AllowSkip); } /// @@ -159,10 +170,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The vector to store the result of the division. protected override void DoDivideByThis(double dividend, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, dividend / At(index)); - } + Map(x => dividend/x, result, Zeros.Include); } /// @@ -172,10 +180,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The vector to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) * other.At(index)); - } + Map2((x, y) => x*y, other, result, Zeros.AllowSkip); } /// @@ -185,10 +190,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The vector to store the result of the pointwise division. protected override void DoPointwiseDivide(Vector divisor, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) / divisor.At(index)); - } + Map2((x, y) => x/y, divisor, result, Zeros.Include); } /// @@ -198,7 +200,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(double exponent, Vector result) { - Map(x => Math.Pow(x, exponent), result, Zeros.AllowSkip); + Map(x => Math.Pow(x, exponent), result, exponent > 0.0 ? Zeros.AllowSkip : Zeros.Include); } /// @@ -209,10 +211,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The result of the modulus. protected override void DoPointwiseModulus(Vector divisor, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, Euclid.Modulus(At(index), divisor.At(index))); - } + Map2(Euclid.Modulus, divisor, result, Zeros.Include); } /// @@ -223,10 +222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The result of the modulus. protected override void DoPointwiseRemainder(Vector divisor, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index)%divisor.At(index)); - } + Map2(Euclid.Remainder, divisor, result, Zeros.Include); } /// @@ -280,10 +276,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoModulus(double divisor, Vector result) { - for (int i = 0; i < Count; i++) - { - result.At(i, Euclid.Modulus(At(i), divisor)); - } + Map(x => Euclid.Modulus(x, divisor), result, Zeros.Include); } /// @@ -294,10 +287,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoModulusByThis(double dividend, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, Euclid.Modulus(dividend, At(index))); - } + Map(x => Euclid.Modulus(dividend, x), result, Zeros.Include); } /// @@ -308,10 +298,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoRemainder(double divisor, Vector result) { - for (int i = 0; i < Count; i++) - { - result.At(i, At(i)%divisor); - } + Map(x => Euclid.Remainder(x, divisor), result, Zeros.Include); } /// @@ -322,10 +309,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoRemainderByThis(double dividend, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, dividend%At(index)); - } + Map(x => Euclid.Remainder(dividend, x), result, Zeros.Include); } /// @@ -459,32 +443,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double return Math.Pow(sum, 1.0/p); } - /// - /// Conjugates vector and save result to - /// - /// Target vector - protected override void DoConjugate(Vector result) - { - if (ReferenceEquals(this, result)) - { - return; - } - - CopyTo(result); - } - - /// - /// Negates vector and saves result to - /// - /// Target vector - protected override void DoNegate(Vector result) - { - for (var index = 0; index < Count; index++) - { - result.At(index, -At(index)); - } - } - /// /// Returns the index of the absolute maximum element. /// diff --git a/src/Numerics/LinearAlgebra/Matrix.cs b/src/Numerics/LinearAlgebra/Matrix.cs index 4c67a260..10ac4171 100644 --- a/src/Numerics/LinearAlgebra/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Matrix.cs @@ -1736,6 +1736,24 @@ namespace MathNet.Numerics.LinearAlgebra return EnumerateColumns().Aggregate(f); } + /// + /// Applies a function to each value pair of two matrices and replaces the value in the result vector. + /// + public void Map2(Func f, Matrix other, Matrix result, Zeros zeros = Zeros.AllowSkip) + { + Storage.Map2To(result.Storage, other.Storage, f, zeros, ExistingData.Clear); + } + + /// + /// Applies a function to each value pair of two matrices and returns the results as a new vector. + /// + public Matrix Map2(Func f, Matrix other, Zeros zeros = Zeros.AllowSkip) + { + var result = Build.SameAs(this); + Storage.Map2To(result.Storage, other.Storage, f, zeros, ExistingData.AssumeZeros); + return result; + } + /// /// Applies a function to update the status with each value pair of two matrices and returns the resulting status. /// diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs index dc020b5a..599fec1c 100644 --- a/src/Numerics/LinearAlgebra/Single/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -58,198 +58,36 @@ namespace MathNet.Numerics.LinearAlgebra.Single MapInplace(x => Math.Abs(x) < threshold ? 0f : x, Zeros.AllowSkip); } - /// Calculates the induced L1 norm of this matrix. - /// The maximum absolute column sum of the matrix. - public override double L1Norm() - { - var norm = 0d; - for (var j = 0; j < ColumnCount; j++) - { - var s = 0d; - for (var i = 0; i < RowCount; i++) - { - s += Math.Abs(At(i, j)); - } - norm = Math.Max(norm, s); - } - return norm; - } - - /// Calculates the induced infinity norm of this matrix. - /// The maximum absolute row sum of the matrix. - public override double InfinityNorm() - { - var norm = 0d; - for (var i = 0; i < RowCount; i++) - { - var s = 0d; - for (var j = 0; j < ColumnCount; j++) - { - s += Math.Abs(At(i, j)); - } - norm = Math.Max(norm, s); - } - return norm; - } - - /// Calculates the entry-wise Frobenius norm of this matrix. - /// The square root of the sum of the squared values. - public override double FrobeniusNorm() - { - var transpose = Transpose(); - var aat = this*transpose; - var norm = 0d; - for (var i = 0; i < RowCount; i++) - { - norm += aat.At(i, i); - } - return Math.Sqrt(norm); - } - - /// - /// Calculates the p-norms of all row vectors. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override Vector RowNorms(double norm) - { - if (norm <= 0.0) - { - throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); - } - - var ret = new double[RowCount]; - if (norm == 2.0) - { - Storage.FoldByRowUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); - } - else if (norm == 1.0) - { - Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); - } - else if (double.IsPositiveInfinity(norm)) - { - Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); - } - else - { - double invnorm = 1.0/norm; - Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); - } - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the p-norms of all column vectors. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) - /// - public override Vector ColumnNorms(double norm) - { - if (norm <= 0.0) - { - throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); - } - - var ret = new double[ColumnCount]; - if (norm == 2.0) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); - } - else if (norm == 1.0) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); - } - else if (double.IsPositiveInfinity(norm)) - { - Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); - } - else - { - double invnorm = 1.0/norm; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); - } - return Vector.Build.Dense(ret); - } - /// - /// Normalizes all row vectors to a unit p-norm. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// Returns the conjugate transpose of this matrix. /// - public override sealed Matrix NormalizeRows(double norm) + /// The conjugate transpose of this matrix. + public override sealed Matrix ConjugateTranspose() { - var norminv = ((DenseVectorStorage)RowNorms(norm).Storage).Data; - for (int i = 0; i < norminv.Length; i++) - { - norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; - } - - var result = Build.SameAs(this, RowCount, ColumnCount); - Storage.MapIndexedTo(result.Storage, (i, j, x) => (float)norminv[i]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); - return result; + return Transpose(); } /// - /// Normalizes all column vectors to a unit p-norm. - /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// Complex conjugates each element of this matrix and place the results into the result matrix. /// - public override sealed Matrix NormalizeColumns(double norm) + /// The result of the conjugation. + protected override sealed void DoConjugate(Matrix result) { - var norminv = ((DenseVectorStorage)ColumnNorms(norm).Storage).Data; - for (int i = 0; i < norminv.Length; i++) + if (ReferenceEquals(this, result)) { - norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; + return; } - var result = Build.SameAs(this, RowCount, ColumnCount); - Storage.MapIndexedTo(result.Storage, (i, j, x) => (float)norminv[j]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); - return result; - } - /// - /// Calculates the value sum of each row vector. - /// - public override Vector RowSums() - { - var ret = new float[RowCount]; - Storage.FoldByRowUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the absolute value sum of each row vector. - /// - public override Vector RowAbsoluteSums() - { - var ret = new float[RowCount]; - Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the value sum of each column vector. - /// - public override Vector ColumnSums() - { - var ret = new float[ColumnCount]; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); - } - - /// - /// Calculates the absolute value sum of each column vector. - /// - public override Vector ColumnAbsoluteSums() - { - var ret = new float[ColumnCount]; - Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); - return Vector.Build.Dense(ret); + CopyTo(result); } /// - /// Returns the conjugate transpose of this matrix. + /// Negate each element of this matrix and place the results into the result matrix. /// - /// The conjugate transpose of this matrix. - public override sealed Matrix ConjugateTranspose() + /// The result of the negation. + protected override void DoNegate(Matrix result) { - return Transpose(); + Map(x => -x, result, Zeros.AllowSkip); } /// @@ -259,13 +97,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The matrix to store the result of the addition. protected override void DoAdd(float scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) + scalar); - } - } + Map(x => x + scalar, result, Zeros.Include); } /// @@ -277,13 +109,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// If the two matrices don't have the same dimensions. protected override void DoAdd(Matrix other, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) + other.At(i, j)); - } - } + Map2((x, y) => x + y, other, result, Zeros.AllowSkip); } /// @@ -293,13 +119,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The matrix to store the result of the subtraction. protected override void DoSubtract(float scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) - scalar); - } - } + Map(x => x - scalar, result, Zeros.Include); } /// @@ -311,13 +131,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// If the two matrices don't have the same dimensions. protected override void DoSubtract(Matrix other, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j) - other.At(i, j)); - } - } + Map2((x, y) => x - y, other, result, Zeros.AllowSkip); } /// @@ -327,13 +141,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The matrix to store the result of the multiplication. protected override void DoMultiply(float scalar, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, At(i, j)*scalar); - } - } + Map(x => x*scalar, result, Zeros.AllowSkip); } /// @@ -382,23 +190,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The matrix to store the result of the division. protected override void DoDivide(float divisor, Matrix result) { - DoMultiply(1.0f/divisor, result); + Map(x => x/divisor, result, divisor == 0.0f ? Zeros.Include : Zeros.AllowSkip); } /// /// Divides a scalar by each element of the matrix and stores the result in the result matrix. /// - /// The scalar to add. + /// The scalar to divide by each element of the matrix. /// The matrix to store the result of the division. protected override void DoDivideByThis(float dividend, Matrix result) { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, dividend/At(i, j)); - } - } + Map(x => dividend/x, result, Zeros.Include); } /// @@ -499,13 +301,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Matrix to store the results in. protected override void DoModulus(float divisor, Matrix result) { - for (var row = 0; row < RowCount; row++) - { - for (var column = 0; column < ColumnCount; column++) - { - result.At(row, column, Euclid.Modulus(At(row, column), divisor)); - } - } + Map(x => Euclid.Modulus(x, divisor), result, Zeros.Include); } /// @@ -516,13 +312,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoModulusByThis(float dividend, Matrix result) { - for (var row = 0; row < RowCount; row++) - { - for (var column = 0; column < ColumnCount; column++) - { - result.At(row, column, Euclid.Modulus(dividend, At(row, column))); - } - } + Map(x => Euclid.Modulus(dividend, x), result, Zeros.Include); } /// @@ -533,13 +323,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Matrix to store the results in. protected override void DoRemainder(float divisor, Matrix result) { - for (var row = 0; row < RowCount; row++) - { - for (var column = 0; column < ColumnCount; column++) - { - result.At(row, column, At(row, column)%divisor); - } - } + Map(x => Euclid.Remainder(x, divisor), result, Zeros.Include); } /// @@ -550,42 +334,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoRemainderByThis(float dividend, Matrix result) { - for (var row = 0; row < RowCount; row++) - { - for (var column = 0; column < ColumnCount; column++) - { - result.At(row, column, dividend%At(row, column)); - } - } - } - - /// - /// Negate each element of this matrix and place the results into the result matrix. - /// - /// The result of the negation. - protected override void DoNegate(Matrix result) - { - for (var i = 0; i < RowCount; i++) - { - for (var j = 0; j < ColumnCount; j++) - { - result.At(i, j, -At(i, j)); - } - } - } - - /// - /// Complex conjugates each element of this matrix and place the results into the result matrix. - /// - /// The result of the conjugation. - protected override sealed void DoConjugate(Matrix result) - { - if (ReferenceEquals(this, result)) - { - return; - } - - CopyTo(result); + Map(x => Euclid.Remainder(dividend, x), result, Zeros.Include); } /// @@ -595,13 +344,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The matrix to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Matrix other, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, At(i, j)*other.At(i, j)); - } - } + Map2((x, y) => x*y, other, result, Zeros.AllowSkip); } /// @@ -611,13 +354,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The matrix to store the result of the pointwise division. protected override void DoPointwiseDivide(Matrix divisor, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, At(i, j)/divisor.At(i, j)); - } - } + Map2((x, y) => x/y, divisor, result, Zeros.Include); } /// @@ -627,7 +364,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(float exponent, Matrix result) { - Map(x => (float)Math.Pow(x, exponent), result, Zeros.AllowSkip); + Map(x => (float)Math.Pow(x, exponent), result, exponent > 0.0f ? Zeros.AllowSkip : Zeros.Include); } /// @@ -638,13 +375,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The result of the modulus. protected override void DoPointwiseModulus(Matrix divisor, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, Euclid.Modulus(At(i, j), divisor.At(i, j))); - } - } + Map2(Euclid.Modulus, divisor, result, Zeros.Include); } /// @@ -655,13 +386,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The result of the modulus. protected override void DoPointwiseRemainder(Matrix divisor, Matrix result) { - for (var j = 0; j < ColumnCount; j++) - { - for (var i = 0; i < RowCount; i++) - { - result.At(i, j, At(i, j)%divisor.At(i, j)); - } - } + Map2(Euclid.Remainder, divisor, result, Zeros.Include); } /// @@ -703,6 +428,192 @@ namespace MathNet.Numerics.LinearAlgebra.Single return sum; } + /// Calculates the induced L1 norm of this matrix. + /// The maximum absolute column sum of the matrix. + public override double L1Norm() + { + var norm = 0d; + for (var j = 0; j < ColumnCount; j++) + { + var s = 0d; + for (var i = 0; i < RowCount; i++) + { + s += Math.Abs(At(i, j)); + } + norm = Math.Max(norm, s); + } + return norm; + } + + /// Calculates the induced infinity norm of this matrix. + /// The maximum absolute row sum of the matrix. + public override double InfinityNorm() + { + var norm = 0d; + for (var i = 0; i < RowCount; i++) + { + var s = 0d; + for (var j = 0; j < ColumnCount; j++) + { + s += Math.Abs(At(i, j)); + } + norm = Math.Max(norm, s); + } + return norm; + } + + /// Calculates the entry-wise Frobenius norm of this matrix. + /// The square root of the sum of the squared values. + public override double FrobeniusNorm() + { + var transpose = Transpose(); + var aat = this*transpose; + var norm = 0d; + for (var i = 0; i < RowCount; i++) + { + norm += aat.At(i, i); + } + return Math.Sqrt(norm); + } + + /// + /// Calculates the p-norms of all row vectors. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override Vector RowNorms(double norm) + { + if (norm <= 0.0) + { + throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); + } + + var ret = new double[RowCount]; + if (norm == 2.0) + { + Storage.FoldByRowUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); + } + else if (norm == 1.0) + { + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); + } + else if (double.IsPositiveInfinity(norm)) + { + Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); + } + else + { + double invnorm = 1.0/norm; + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); + } + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the p-norms of all column vectors. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override Vector ColumnNorms(double norm) + { + if (norm <= 0.0) + { + throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); + } + + var ret = new double[ColumnCount]; + if (norm == 2.0) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); + } + else if (norm == 1.0) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); + } + else if (double.IsPositiveInfinity(norm)) + { + Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); + } + else + { + double invnorm = 1.0/norm; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); + } + return Vector.Build.Dense(ret); + } + + /// + /// Normalizes all row vectors to a unit p-norm. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override sealed Matrix NormalizeRows(double norm) + { + var norminv = ((DenseVectorStorage)RowNorms(norm).Storage).Data; + for (int i = 0; i < norminv.Length; i++) + { + norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; + } + + var result = Build.SameAs(this, RowCount, ColumnCount); + Storage.MapIndexedTo(result.Storage, (i, j, x) => (float)norminv[i]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); + return result; + } + + /// + /// Normalizes all column vectors to a unit p-norm. + /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm) + /// + public override sealed Matrix NormalizeColumns(double norm) + { + var norminv = ((DenseVectorStorage)ColumnNorms(norm).Storage).Data; + for (int i = 0; i < norminv.Length; i++) + { + norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i]; + } + + var result = Build.SameAs(this, RowCount, ColumnCount); + Storage.MapIndexedTo(result.Storage, (i, j, x) => (float)norminv[j]*x, Zeros.AllowSkip, ExistingData.AssumeZeros); + return result; + } + + /// + /// Calculates the value sum of each row vector. + /// + public override Vector RowSums() + { + var ret = new float[RowCount]; + Storage.FoldByRowUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the absolute value sum of each row vector. + /// + public override Vector RowAbsoluteSums() + { + var ret = new float[RowCount]; + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the value sum of each column vector. + /// + public override Vector ColumnSums() + { + var ret = new float[ColumnCount]; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x, (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + + /// + /// Calculates the absolute value sum of each column vector. + /// + public override Vector ColumnAbsoluteSums() + { + var ret = new float[ColumnCount]; + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); + return Vector.Build.Dense(ret); + } + /// /// Evaluates whether this matrix is hermitian (conjugate symmetric). /// diff --git a/src/Numerics/LinearAlgebra/Single/Vector.cs b/src/Numerics/LinearAlgebra/Single/Vector.cs index 0ff911fa..1a8de510 100644 --- a/src/Numerics/LinearAlgebra/Single/Vector.cs +++ b/src/Numerics/LinearAlgebra/Single/Vector.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -56,6 +56,29 @@ namespace MathNet.Numerics.LinearAlgebra.Single MapInplace(x => Math.Abs(x) < threshold ? 0f : x, Zeros.AllowSkip); } + /// + /// Conjugates vector and save result to + /// + /// Target vector + protected override sealed void DoConjugate(Vector result) + { + if (ReferenceEquals(this, result)) + { + return; + } + + CopyTo(result); + } + + /// + /// Negates vector and saves result to + /// + /// Target vector + protected override void DoNegate(Vector result) + { + Map(x => -x, result, Zeros.AllowSkip); + } + /// /// Adds a scalar to each element of the vector and stores the result in the result vector. /// @@ -67,10 +90,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoAdd(float scalar, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) + scalar); - } + Map(x => x + scalar, result, Zeros.Include); } /// @@ -84,10 +104,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoAdd(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) + other.At(index)); - } + Map2((x, y) => x + y, other, result, Zeros.AllowSkip); } /// @@ -101,7 +118,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoSubtract(float scalar, Vector result) { - DoAdd(-scalar, result); + Map(x => x - scalar, result, Zeros.Include); } /// @@ -115,10 +132,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoSubtract(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) - other.At(index)); - } + Map2((x, y) => x - y, other, result, Zeros.AllowSkip); } /// @@ -132,10 +146,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoMultiply(float scalar, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) * scalar); - } + Map(x => x*scalar, result, Zeros.AllowSkip); } /// @@ -149,7 +160,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoDivide(float divisor, Vector result) { - DoMultiply(1 / divisor, result); + Map(x => x/divisor, result, divisor == 0.0f ? Zeros.Include : Zeros.AllowSkip); } /// @@ -159,10 +170,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The vector to store the result of the division. protected override void DoDivideByThis(float dividend, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, dividend / At(index)); - } + Map(x => dividend/x, result, Zeros.Include); } /// @@ -172,10 +180,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The vector to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Vector other, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) * other.At(index)); - } + Map2((x, y) => x*y, other, result, Zeros.AllowSkip); } /// @@ -185,10 +190,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The vector to store the result of the pointwise division. protected override void DoPointwiseDivide(Vector divisor, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index) / divisor.At(index)); - } + Map2((x, y) => x/y, divisor, result, Zeros.Include); } /// @@ -198,7 +200,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(float exponent, Vector result) { - Map(x => (float)Math.Pow(x, exponent), result, Zeros.AllowSkip); + Map(x => (float)Math.Pow(x, exponent), result, exponent > 0.0f ? Zeros.AllowSkip : Zeros.Include); } /// @@ -209,10 +211,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The result of the modulus. protected override void DoPointwiseModulus(Vector divisor, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, Euclid.Modulus(At(index), divisor.At(index))); - } + Map2(Euclid.Modulus, divisor, result, Zeros.Include); } /// @@ -223,10 +222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The result of the modulus. protected override void DoPointwiseRemainder(Vector divisor, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, At(index)%divisor.At(index)); - } + Map2(Euclid.Remainder, divisor, result, Zeros.Include); } /// @@ -280,10 +276,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoModulus(float divisor, Vector result) { - for (int i = 0; i < Count; i++) - { - result.At(i, Euclid.Modulus(At(i), divisor)); - } + Map(x => Euclid.Modulus(x, divisor), result, Zeros.Include); } /// @@ -294,10 +287,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoModulusByThis(float dividend, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, Euclid.Modulus(dividend, At(index))); - } + Map(x => Euclid.Modulus(dividend, x), result, Zeros.Include); } /// @@ -308,10 +298,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoRemainder(float divisor, Vector result) { - for (int i = 0; i < Count; i++) - { - result.At(i, At(i)%divisor); - } + Map(x => Euclid.Remainder(x, divisor), result, Zeros.Include); } /// @@ -322,10 +309,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoRemainderByThis(float dividend, Vector result) { - for (var index = 0; index < Count; index++) - { - result.At(index, dividend%At(index)); - } + Map(x => Euclid.Remainder(dividend, x), result, Zeros.Include); } /// @@ -459,32 +443,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single return Math.Pow(sum, 1.0/p); } - /// - /// Conjugates vector and save result to - /// - /// Target vector - protected override void DoConjugate(Vector result) - { - if (ReferenceEquals(this, result)) - { - return; - } - - CopyTo(result); - } - - /// - /// Negates vector and saves result to - /// - /// Target vector - protected override void DoNegate(Vector result) - { - for (var index = 0; index < Count; index++) - { - result.At(index, -At(index)); - } - } - /// /// Returns the index of the absolute maximum element. /// diff --git a/src/Numerics/LinearAlgebra/Vector.cs b/src/Numerics/LinearAlgebra/Vector.cs index 1988103a..e02c06cd 100644 --- a/src/Numerics/LinearAlgebra/Vector.cs +++ b/src/Numerics/LinearAlgebra/Vector.cs @@ -369,6 +369,7 @@ namespace MathNet.Numerics.LinearAlgebra /// public void MapInplace(Func f, Zeros zeros = Zeros.AllowSkip) { + // TODO: actual in-place Storage.MapToUnchecked(Storage, f, zeros, ExistingData.AssumeZeros); } @@ -380,6 +381,7 @@ namespace MathNet.Numerics.LinearAlgebra /// public void MapIndexedInplace(Func f, Zeros zeros = Zeros.AllowSkip) { + // TODO: actual in-place Storage.MapIndexedToUnchecked(Storage, f, zeros, ExistingData.AssumeZeros); } @@ -392,6 +394,7 @@ namespace MathNet.Numerics.LinearAlgebra where TU : struct, IEquatable, IFormattable { // TODO: in v4 update this method to replace TU with T (consistent with Matrix, see MapConvert) + // then automatically do in-place if possible. Storage.MapTo(result.Storage, f, zeros, zeros == Zeros.Include ? ExistingData.AssumeZeros : ExistingData.Clear); } @@ -405,6 +408,7 @@ namespace MathNet.Numerics.LinearAlgebra where TU : struct, IEquatable, IFormattable { // TODO: in v4 update this method to replace TU with T (consistent with Matrix, see MapIndexedConvert) + // then automatically do in-place if possible. Storage.MapIndexedTo(result.Storage, f, zeros, zeros == Zeros.Include ? ExistingData.AssumeZeros : ExistingData.Clear); }