From 6242f534fd26767c346a6df033cae9c9fd610146 Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Sun, 6 Jun 2010 20:39:39 +0100 Subject: [PATCH] Added Permutation class. Added Matrix.PermuteRows and Matrix.PermuteColumns. Further implementation of LU decomposition. --- .../LinearAlgebra/Double/Factorization/LU.cs | 29 +++------- src/Numerics/LinearAlgebra/Double/Matrix.cs | 58 +++++++++++++++++++ src/Numerics/Numerics.csproj | 1 + src/Numerics/Properties/Resources.Designer.cs | 9 +++ src/Numerics/Properties/Resources.resx | 3 + src/Silverlight/Silverlight.csproj | 3 + .../Double/Factorization/LUTests.cs | 3 +- .../LinearAlgebraTests/Double/MatrixTests.cs | 50 ++++++++++++++++ src/UnitTests/TrigonometryTest.cs | 2 +- src/UnitTests/UnitTests.csproj | 1 + 10 files changed, 136 insertions(+), 23 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs index f9a6b3ba..b1c98350 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs @@ -95,6 +95,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization get { return mFactors.GetUpperTriangle(); } } + /// + /// Return the permutation applied to LU factorization. + /// + public virtual Permutation P + { + get { return Permutation.FromInversions(mPivots); } + } + /// /// The determinant of the matrix for which the LU factorization was computed. /// @@ -167,26 +175,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The right hand side vector, b. /// The left hand side , x. public abstract void Solve(Vector input, Vector result); - - /// - /// Pivot a matrix according to this LU decomposition. - /// - /// The matrix to pivot. - public void Pivot(Matrix data) - { - for (int i = 0; i < mPivots.Length; i++) - { - if (mPivots[i] != i) - { - int p = mPivots[i]; - for (int j = 0; j < data.ColumnCount; j++) - { - double temp = data.At(p, j); - data.At(p, j, data.At(i, j)); - data.At(i, j, temp); - } - } - } - } } } diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index c5da5b90..a873ca99 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -760,5 +760,63 @@ namespace MathNet.Numerics.LinearAlgebra.Double return ret; } + + /// + /// Permute the rows of a matrix according to a permutation. + /// + /// The row permutation to apply to this matrix. + public virtual void PermuteRows(Permutation p) + { + if (p.Dimension != this.RowCount) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "p"); + } + + // Get a sequence of inversions from the permutation. + int[] inv = p.ToInversions(); + + for (int i = 0; i < p.Dimension; i++) + { + if (inv[i] != i) + { + int q = inv[i]; + for (int j = 0; j < this.ColumnCount; j++) + { + double temp = At(q, j); + At(q, j, At(i, j)); + At(i, j, temp); + } + } + } + } + + /// + /// Permute the columns of a matrix according to a permutation. + /// + /// The column permutation to apply to this matrix. + public virtual void PermuteColumns(Permutation p) + { + if (p.Dimension != this.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "p"); + } + + // Get a sequence of inversions from the permutation. + int[] inv = p.ToInversions(); + + for (int i = 0; i < p.Dimension; i++) + { + if (inv[i] != i) + { + int q = inv[i]; + for (int j = 0; j < this.RowCount; j++) + { + double temp = At(j, q); + At(j, q, At(j, i)); + At(j, i, temp); + } + } + } + } } } \ No newline at end of file diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index fe5ca4f5..5179f244 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -101,6 +101,7 @@ + diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 4d595b4d..ee6eae14 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -573,6 +573,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to The integer array does not represent a valid permutation.. + /// + internal static string PermutationAsIntArrayInvalid { + get { + return ResourceManager.GetString("PermutationAsIntArrayInvalid", resourceCulture); + } + } + /// /// Looks up a localized string similar to The sampler's proposal distribution is not upper bounding the target density.. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index 0e04d29e..50742f13 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -300,4 +300,7 @@ Arguments must be different objects. + + The integer array does not represent a valid permutation. + \ No newline at end of file diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj index ab10abe5..1e59e3cf 100644 --- a/src/Silverlight/Silverlight.csproj +++ b/src/Silverlight/Silverlight.csproj @@ -245,6 +245,9 @@ NumberTheory\IntegerTheory.Euclid.cs + + Permutation.cs + Precision.cs diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs index e97c7c04..17c8f79d 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs @@ -148,7 +148,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization // Make sure the cholesky factor times it's transpose is the original matrix. var XfromLU = L * U; - lu.Pivot(XfromLU); + var Pinv = lu.P.Inverse(); + XfromLU.PermuteRows(Pinv); for (int i = 0; i < XfromLU.RowCount; i++) { for (int j = 0; j < XfromLU.ColumnCount; j++) diff --git a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs index 14d683d1..7f2a880d 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs @@ -481,5 +481,55 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double } } } + + [Test] + [Row("Singular3x3")] + [Row("Square3x3")] + [Row("Tall3x2")] + [MultipleAsserts] + public void CanPermuteMatrixRows(string name) + { + var matrix = CreateMatrix(testData2D[name]); + var matrixp = CreateMatrix(testData2D[name]); + + var permutation = new Permutation(new int[] { 2, 0, 1 }); + matrixp.PermuteRows(permutation); + + Assert.AreNotSame(matrix, matrixp); + Assert.AreEqual(matrix.RowCount, matrixp.RowCount); + Assert.AreEqual(matrix.ColumnCount, matrixp.ColumnCount); + for (var i = 0; i < matrix.RowCount; i++) + { + for (var j = 0; j < matrix.ColumnCount; j++) + { + Assert.AreEqual(matrix[i, j], matrixp[permutation[i], j]); + } + } + } + + [Test] + [Row("Singular3x3")] + [Row("Square3x3")] + [Row("Wide2x3")] + [MultipleAsserts] + public void CanPermuteMatrixColumns(string name) + { + var matrix = CreateMatrix(testData2D[name]); + var matrixp = CreateMatrix(testData2D[name]); + + var permutation = new Permutation(new int[] { 2, 0, 1 }); + matrixp.PermuteColumns(permutation); + + Assert.AreNotSame(matrix, matrixp); + Assert.AreEqual(matrix.RowCount, matrixp.RowCount); + Assert.AreEqual(matrix.ColumnCount, matrixp.ColumnCount); + for (var i = 0; i < matrix.RowCount; i++) + { + for (var j = 0; j < matrix.ColumnCount; j++) + { + Assert.AreEqual(matrix[i, j], matrixp[i, permutation[j]]); + } + } + } } } \ No newline at end of file diff --git a/src/UnitTests/TrigonometryTest.cs b/src/UnitTests/TrigonometryTest.cs index c780911d..c508cc61 100644 --- a/src/UnitTests/TrigonometryTest.cs +++ b/src/UnitTests/TrigonometryTest.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index ab2aef88..9afe6344 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -90,6 +90,7 @@ +