diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index e624b098..937787ef 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -204,8 +204,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentException(Resources.ArgumentMatrixDimensions, "target"); } - // TODO this assumes that all entries matter; if "this" is a sparse matrix, - // we might be able to optimize the copying a bit. for (var i = 0; i < this.RowCount; i++) { for (var j = 0; j < this.ColumnCount; j++) diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs index 3ef77e28..8e94a045 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com +// // Copyright (c) 2009-2010 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -27,31 +31,32 @@ namespace MathNet.Numerics.LinearAlgebra.Double { using System; + using Distributions; using Properties; using Threading; - + /// - /// Sparse Matrix implementation + /// A Matrix class with sparse storage. The underlying storage scheme is 3-array CSR Format. /// - public class SparseMatrix : Matrix + public class SparseMatrix : Matrix { /// /// Object for use in "lock" /// - private readonly object lockObject = new object(); + private readonly object _lockObject = new object(); /// /// The array containing the row indices of the existing rows. Element "j" of the array gives the index of the - /// element in the array that is first non-zero element in a row "j" + /// element in the array that is first non-zero element in a row "j" /// - private readonly int[] rowIndex = new int[0]; + private readonly int[] _rowIndex = new int[0]; /// /// Array that contains the non-zero elements of matrix. Values of the non-zero elements of matrix are mapped into the values /// array using the row-major storage mapping described in a compressed sparse row (CSR) format. /// - private double[] nonZeroValues = new double[0]; - + private double[] _nonZeroValues = new double[0]; + /// /// Gets the number of non zero elements in the matrix. /// @@ -61,12 +66,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double get; private set; } - + /// /// An array containing the column indices of the non-zero values. Element "I" of the array - /// is the number of the column in matrix that contains the I-th value in the array. + /// is the number of the column in matrix that contains the I-th value in the array. /// - private int[] columnIndices = new int[0]; + private int[] _columnIndices = new int[0]; /// /// Initializes a new instance of the class. @@ -79,9 +84,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// public SparseMatrix(int rows, int columns) : base(rows, columns) { - this.rowIndex = new int[rows]; + this._rowIndex = new int[rows]; } - + /// /// Initializes a new instance of the class. This matrix is square with a given size. /// @@ -110,11 +115,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double return; } - this.NonZerosCount = rows * columns; - this.nonZeroValues = new double[this.NonZerosCount]; - this.columnIndices = new int[this.NonZerosCount]; + NonZerosCount = rows * columns; + this._nonZeroValues = new double[NonZerosCount]; + this._columnIndices = new int[NonZerosCount]; - for (int i = 0, j = 0; i < this.nonZeroValues.Length; i++, j++) + for (int i = 0, j = 0; i < this._nonZeroValues.Length; i++, j++) { // Reset column position to "0" if (j == columns) @@ -122,14 +127,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double j = 0; } - this.nonZeroValues[i] = value; - this.columnIndices[i] = j; + this._nonZeroValues[i] = value; + this._columnIndices[i] = j; } - + // Set proper row pointers - for (var i = 0; i < this.rowIndex.Length; i++) + for (var i = 0; i < this._rowIndex.Length; i++) { - this.rowIndex[i] = ((i + 1) * columns) - columns; + this._rowIndex[i] = ((i + 1) * columns) - columns; } } @@ -139,7 +144,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The number of rows. /// The number of columns. /// The one dimensional array to create this matrix from. This array should store the matrix in column-major order. - /// If lenght is less than * . + /// If length is less than * . /// public SparseMatrix(int rows, int columns, double[] array) : this(rows, columns) { @@ -204,6 +209,340 @@ namespace MathNet.Numerics.LinearAlgebra.Double return new SparseVector(size); } + /// + /// Returns a new matrix containing the lower triangle of this matrix. + /// + /// The lower triangle of this matrix. + public override Matrix LowerTriangle() + { + var result = this.CreateMatrix(this.RowCount, this.ColumnCount); + for (var row = 0; row < result.RowCount; row++) + { + var startIndex = this._rowIndex[row]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[row + 1] : this.NonZerosCount; + for (var j = startIndex; j < endIndex; j++) + { + if (row >= this._columnIndices[j]) + { + result.At(row, this._columnIndices[j], this._nonZeroValues[j]); + } + } + } + + return result; + } + + /// + /// Puts the lower triangle of this matrix into the result matrix. + /// + /// Where to store the lower triangle. + /// If is . + /// If the result matrix's dimensions are not the same as this matrix. + public override void LowerTriangle(Matrix result) + { + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (result.RowCount != this.RowCount || result.ColumnCount != this.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); + } + + if (ReferenceEquals(this, result)) + { + var tmp = result.CreateMatrix(result.RowCount, result.ColumnCount); + this.LowerTriangle(tmp); + tmp.CopyTo(result); + } + else + { + result.Clear(); + for (var row = 0; row < result.RowCount; row++) + { + var startIndex = this._rowIndex[row]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[row + 1] : this.NonZerosCount; + for (var j = startIndex; j < endIndex; j++) + { + if (row >= this._columnIndices[j]) + { + result.At(row, this._columnIndices[j], this._nonZeroValues[j]); + } + } + } + } + } + + /// + /// Returns a new matrix containing the upper triangle of this matrix. + /// + /// The upper triangle of this matrix. + public override Matrix UpperTriangle() + { + var result = this.CreateMatrix(this.RowCount, this.ColumnCount); + for (var row = 0; row < result.RowCount; row++) + { + var startIndex = this._rowIndex[row]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[row + 1] : this.NonZerosCount; + for (var j = startIndex; j < endIndex; j++) + { + if (row <= this._columnIndices[j]) + { + result.At(row, this._columnIndices[j], this._nonZeroValues[j]); + } + } + } + + return result; + } + + /// + /// Puts the upper triangle of this matrix into the result matrix. + /// + /// Where to store the lower triangle. + /// If is . + /// If the result matrix's dimensions are not the same as this matrix. + public override void UpperTriangle(Matrix result) + { + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (result.RowCount != this.RowCount || result.ColumnCount != this.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); + } + + if (ReferenceEquals(this, result)) + { + var tmp = result.CreateMatrix(result.RowCount, result.ColumnCount); + this.UpperTriangle(tmp); + tmp.CopyTo(result); + } + else + { + result.Clear(); + for (var row = 0; row < result.RowCount; row++) + { + var startIndex = this._rowIndex[row]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[row + 1] : this.NonZerosCount; + for (var j = startIndex; j < endIndex; j++) + { + if (row <= this._columnIndices[j]) + { + result.At(row, this._columnIndices[j], this._nonZeroValues[j]); + } + } + } + } + } + + /// + /// Creates a matrix that contains the values from the requested sub-matrix. + /// + /// The row to start copying from. + /// The number of rows to copy. Must be positive. + /// The column to start copying from. + /// The number of columns to copy. Must be positive. + /// The requested sub-matrix. + /// If: is + /// negative, or greater than or equal to the number of rows. + /// is negative, or greater than or equal to the number + /// of columns. + /// (columnIndex + columnLength) >= Columns + /// (rowIndex + rowLength) >= Rows + /// If or + /// is not positive. + public override Matrix SubMatrix(int rowIndex, int rowLength, int columnIndex, int columnLength) + { + if (rowIndex >= this.RowCount || rowIndex < 0) + { + throw new ArgumentOutOfRangeException("rowIndex"); + } + + if (columnIndex >= this.ColumnCount || columnIndex < 0) + { + throw new ArgumentOutOfRangeException("columnIndex"); + } + + if (rowLength < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "rowLength"); + } + + if (columnLength < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "columnLength"); + } + + var colMax = columnIndex + columnLength; + var rowMax = rowIndex + rowLength; + + if (rowMax > this.RowCount) + { + throw new ArgumentOutOfRangeException("rowLength"); + } + + if (colMax > this.ColumnCount) + { + throw new ArgumentOutOfRangeException("columnLength"); + } + + var result = this.CreateMatrix(rowLength, columnLength); + + for (int i = rowIndex, row = 0; i < rowMax; i++, row++) + { + var startIndex = this._rowIndex[i]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[i + 1] : this.NonZerosCount; + + for (int j = startIndex; j < endIndex; j++) + { + // check if the column index is in the range + if ((this._columnIndices[j] >= columnIndex) && (this._columnIndices[j] < columnIndex + columnLength)) + { + var column = this._columnIndices[j] - columnIndex; + result[row, column] = this._nonZeroValues[j]; + } + } + } + + return result; + } + + /// + /// Returns a new matrix containing the lower triangle of this matrix. The new matrix + /// does not contain the diagonal elements of this matrix. + /// + /// The lower triangle of this matrix. + public override Matrix StrictlyLowerTriangle() + { + var result = this.CreateMatrix(this.RowCount, this.ColumnCount); + for (var row = 0; row < result.RowCount; row++) + { + var startIndex = this._rowIndex[row]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[row + 1] : this.NonZerosCount; + for (var j = startIndex; j < endIndex; j++) + { + if (row > this._columnIndices[j]) + { + result.At(row, this._columnIndices[j], this._nonZeroValues[j]); + } + } + } + + return result; + } + + /// + /// Puts the strictly lower triangle of this matrix into the result matrix. + /// + /// Where to store the lower triangle. + /// If is . + /// If the result matrix's dimensions are not the same as this matrix. + public override void StrictlyLowerTriangle(Matrix result) + { + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (result.RowCount != this.RowCount || result.ColumnCount != this.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); + } + + if (ReferenceEquals(this, result)) + { + var tmp = result.CreateMatrix(result.RowCount, result.ColumnCount); + this.StrictlyLowerTriangle(tmp); + tmp.CopyTo(result); + } + else + { + result.Clear(); + for (var row = 0; row < result.RowCount; row++) + { + var startIndex = this._rowIndex[row]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[row + 1] : this.NonZerosCount; + for (var j = startIndex; j < endIndex; j++) + { + if (row > this._columnIndices[j]) + { + result.At(row, this._columnIndices[j], this._nonZeroValues[j]); + } + } + } + } + } + + /// + /// Returns a new matrix containing the upper triangle of this matrix. The new matrix + /// does not contain the diagonal elements of this matrix. + /// + /// The upper triangle of this matrix. + public override Matrix StrictlyUpperTriangle() + { + var result = this.CreateMatrix(this.RowCount, this.ColumnCount); + for (var row = 0; row < result.RowCount; row++) + { + var startIndex = this._rowIndex[row]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[row + 1] : this.NonZerosCount; + for (var j = startIndex; j < endIndex; j++) + { + if (row < this._columnIndices[j]) + { + result.At(row, this._columnIndices[j], this._nonZeroValues[j]); + } + } + } + + return result; + } + + /// + /// Puts the strictly upper triangle of this matrix into the result matrix. + /// + /// Where to store the lower triangle. + /// If is . + /// If the result matrix's dimensions are not the same as this matrix. + public override void StrictlyUpperTriangle(Matrix result) + { + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (result.RowCount != this.RowCount || result.ColumnCount != this.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); + } + + if (ReferenceEquals(this, result)) + { + var tmp = result.CreateMatrix(result.RowCount, result.ColumnCount); + this.StrictlyUpperTriangle(tmp); + tmp.CopyTo(result); + } + else + { + result.Clear(); + for (var row = 0; row < result.RowCount; row++) + { + var startIndex = this._rowIndex[row]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[row + 1] : this.NonZerosCount; + for (var j = startIndex; j < endIndex; j++) + { + if (row < this._columnIndices[j]) + { + result.At(row, this._columnIndices[j], this._nonZeroValues[j]); + } + } + } + } + } + /// /// Retrieves the requested element without range checking. /// @@ -218,13 +557,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// public override double At(int row, int column) { - lock (this.lockObject) + lock (this._lockObject) { var index = this.FindItem(row, column); - return index >= 0 ? this.nonZeroValues[index] : 0.0; + return index >= 0 ? this._nonZeroValues[index] : 0.0; } } - + /// /// Sets the value of the given element. /// @@ -239,14 +578,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// public override void At(int row, int column, double value) { - lock (this.lockObject) + lock (this._lockObject) { this.SetValueAt(row, column, value); } } #region Internal methods - CRS storage implementation - /// /// Created this method because we cannot call "virtual At" in constructor of the class, but we need to do it /// @@ -256,69 +594,69 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// WARNING: This method is not thread safe. Use "lock" with it and be sure to avoid deadlocks private void SetValueAt(int row, int column, double value) { - var index = this.FindItem(row, column); - if (index >= 0) - { - // Non-zero item found in matrix - if (value == 0.0) - { - // Delete existing item - this.DeleteItemByIndex(index, row); - } - else + var index = this.FindItem(row, column); + if (index >= 0) { - // Update item - this.nonZeroValues[index] = value; + // Non-zero item found in matrix + if (value == 0.0) + { + // Delete existing item + this.DeleteItemByIndex(index, row); + } + else + { + // Update item + this._nonZeroValues[index] = value; + } } - } - else - { - // Item not found. Add new value - if (value == 0.0) + else { - return; - } - - index = ~index; + // Item not found. Add new value + if (value == 0.0) + { + return; + } - // Check if the storage needs to be increased - if ((this.NonZerosCount == this.nonZeroValues.Length) && (this.NonZerosCount < (this.RowCount * this.ColumnCount))) - { - // Value array is completely full so we increase the size - // Determine the increase in size. We will not grow beyond the size of the matrix - var size = Math.Min(this.nonZeroValues.Length + this.GrowthSize(), this.RowCount * this.ColumnCount); - Array.Resize(ref this.nonZeroValues, size); - Array.Resize(ref this.columnIndices, size); - } + index = ~index; + + // Check if the storage needs to be increased + if ((this.NonZerosCount == this._nonZeroValues.Length) && (this.NonZerosCount < (this.RowCount * this.ColumnCount))) + { + // Value array is completely full so we increase the size + // Determine the increase in size. We will not grow beyond the size of the matrix + var size = Math.Min(this._nonZeroValues.Length + this.GrowthSize(), this.RowCount * this.ColumnCount); + Array.Resize(ref this._nonZeroValues, size); + Array.Resize(ref this._columnIndices, size); + } - // Move all values (with an position larger than index) in the value array to the next position - // move all values (with an position larger than index) in the columIndices array to the next position - for (var i = this.NonZerosCount - 1; i > index - 1; i--) - { - this.nonZeroValues[i + 1] = this.nonZeroValues[i]; - this.columnIndices[i + 1] = this.columnIndices[i]; - } + // Move all values (with an position larger than index) in the value array to the next position + // move all values (with an position larger than index) in the columIndices array to the next position + for (var i = this.NonZerosCount - 1; i > index - 1; i--) + { + this._nonZeroValues[i + 1] = this._nonZeroValues[i]; + this._columnIndices[i + 1] = this._columnIndices[i]; + } - // Add the value and the column index - this.nonZeroValues[index] = value; - this.columnIndices[index] = column; + // Add the value and the column index + this._nonZeroValues[index] = value; + this._columnIndices[index] = column; - // increase the number of non-zero numbers by one - this.NonZerosCount += 1; + // increase the number of non-zero numbers by one + this.NonZerosCount += 1; - // add 1 to all the row indices for rows bigger than rowIndex - // so that they point to the correct part of the value array again. - for (var i = row + 1; i < this.rowIndex.Length; i++) - { - this.rowIndex[i] += 1; - } + // add 1 to all the row indices for rows bigger than rowIndex + // so that they point to the correct part of the value array again. + for (var i = row + 1; i < this._rowIndex.Length; i++) + { + this._rowIndex[i] += 1; + } } } /// /// Delete value from internal storage /// - /// Index of value in nonZeroValues array + /// Index of value in nonZeroValues array /// Row number of matrix /// WARNING: This method is not thread safe. Use "lock" with it and be sure to avoid deadlocks private void DeleteItemByIndex(int itemIndex, int row) @@ -327,29 +665,29 @@ namespace MathNet.Numerics.LinearAlgebra.Double // move all values (with an position larger than index) in the columIndices array to the previous position for (var i = itemIndex + 1; i < this.NonZerosCount; i++) { - this.nonZeroValues[i - 1] = this.nonZeroValues[i]; - this.columnIndices[i - 1] = this.columnIndices[i]; + this._nonZeroValues[i - 1] = this._nonZeroValues[i]; + this._columnIndices[i - 1] = this._columnIndices[i]; } - + // Decrease value in Row - for (var i = row + 1; i < this.rowIndex.Length; i++) + for (var i = row + 1; i < this._rowIndex.Length; i++) { - this.rowIndex[i] -= 1; + this._rowIndex[i] -= 1; } this.NonZerosCount -= 1; // Check if the storage needs to be shrink. This is reasonable to do if // there are a lot of non-zero elements and storage is two times bigger - if ((this.NonZerosCount > 1024) && (this.NonZerosCount < this.nonZeroValues.Length / 2)) + if ((this.NonZerosCount > 1024) && (this.NonZerosCount < this._nonZeroValues.Length / 2)) { - Array.Resize(ref this.nonZeroValues, this.NonZerosCount); - Array.Resize(ref this.columnIndices, this.NonZerosCount); + Array.Resize(ref this._nonZeroValues, this.NonZerosCount); + Array.Resize(ref this._columnIndices, this.NonZerosCount); } } - + /// - /// Find item Index in nonZeroValues array + /// Find item Index in nonZeroValues array /// /// Matrix row index /// Matrix column index @@ -358,37 +696,37 @@ namespace MathNet.Numerics.LinearAlgebra.Double private int FindItem(int row, int column) { // Determin bounds in columnIndices array where this item should be searched (using rowIndex) - var startIndex = this.rowIndex[row]; - var endIndex = row < this.rowIndex.Length - 1 ? this.rowIndex[row + 1] : this.NonZerosCount; - return Array.BinarySearch(this.columnIndices, startIndex, endIndex - startIndex, column); + var startIndex = this._rowIndex[row]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[row + 1] : this.NonZerosCount; + return Array.BinarySearch(this._columnIndices, startIndex, endIndex - startIndex, column); } - + /// - /// Calculate grows size + /// Calculates the amount with which to grow the storage array's if they need to be + /// increased in size. /// - /// Proposed new size + /// The amount grown. private int GrowthSize() { int delta; - if (this.nonZeroValues.Length > 1024) + if (this._nonZeroValues.Length > 1024) { - delta = this.nonZeroValues.Length / 4; + delta = this._nonZeroValues.Length / 4; } else { - if (this.nonZeroValues.Length > 256) + if (this._nonZeroValues.Length > 256) { delta = 512; } else { - delta = this.nonZeroValues.Length > 64 ? 128 : 32; + delta = this._nonZeroValues.Length > 64 ? 128 : 32; } } return delta; } - #endregion /// @@ -397,7 +735,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override void Clear() { this.NonZerosCount = 0; - Array.Clear(this.rowIndex, 0, this.rowIndex.Length); + Array.Clear(this._rowIndex, 0, this._rowIndex.Length); } /// @@ -433,16 +771,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double } // Lets copy only needed data. Portion of needed data is determined by NonZerosCount value - sparseTarget.nonZeroValues = new double[this.NonZerosCount]; - sparseTarget.columnIndices = new int[this.NonZerosCount]; + sparseTarget._nonZeroValues = new double[this.NonZerosCount]; + sparseTarget._columnIndices = new int[this.NonZerosCount]; sparseTarget.NonZerosCount = this.NonZerosCount; - Buffer.BlockCopy(this.nonZeroValues, 0, sparseTarget.nonZeroValues, 0, this.NonZerosCount * Constants.SizeOfDouble); - Buffer.BlockCopy(this.columnIndices, 0, sparseTarget.columnIndices, 0, this.NonZerosCount * Constants.SizeOfInt); - Buffer.BlockCopy(this.rowIndex, 0, sparseTarget.rowIndex, 0, this.RowCount * Constants.SizeOfInt); + Buffer.BlockCopy(this._nonZeroValues, 0, sparseTarget._nonZeroValues, 0, this.NonZerosCount * Constants.SizeOfDouble); + Buffer.BlockCopy(this._columnIndices, 0, sparseTarget._columnIndices, 0, this.NonZerosCount * Constants.SizeOfInt); + Buffer.BlockCopy(this._rowIndex, 0, sparseTarget._rowIndex, 0, this.RowCount * Constants.SizeOfInt); } } - + /// /// Indicates whether the current object is equal to another object of the same type. /// @@ -475,7 +813,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double // If all else fails, perform element wise comparison. for (var index = 0; index < this.NonZerosCount; index++) { - if (!this.nonZeroValues[index].AlmostEqual(sparseMatrix.nonZeroValues[index]) || this.columnIndices[index] != sparseMatrix.columnIndices[index]) + if (!this._nonZeroValues[index].AlmostEqual(sparseMatrix._nonZeroValues[index]) || this._columnIndices[index] != sparseMatrix._columnIndices[index]) { return false; } @@ -483,7 +821,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double return true; } - + /// /// Returns a hash code for this instance. /// @@ -497,9 +835,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var i = 0; i < hashNum; i++) { #if SILVERLIGHT - hash ^= Precision.DoubleToInt64Bits(this.nonZeroValues[i]); + hash ^= Precision.DoubleToInt64Bits(this._nonZeroValues[i]); #else - hash ^= BitConverter.DoubleToInt64Bits(this.nonZeroValues[i]); + hash ^= BitConverter.DoubleToInt64Bits(this._nonZeroValues[i]); #endif } @@ -513,13 +851,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override Matrix Transpose() { var ret = new SparseMatrix(this.ColumnCount, this.RowCount); - + // Do an 'inverse' CopyTo iterate over the rows - for (var i = 0; i < this.rowIndex.Length; i++) + for (var i = 0; i < this._rowIndex.Length; i++) { // Get the begin / end index for the current row - var startIndex = this.rowIndex[i]; - var endIndex = i < this.rowIndex.Length - 1 ? this.rowIndex[i + 1] : this.NonZerosCount; + var startIndex = this._rowIndex[i]; + var endIndex = i < this._rowIndex.Length - 1 ? this._rowIndex[i + 1] : this.NonZerosCount; // Get the values for the current row if (startIndex == endIndex) @@ -530,13 +868,128 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var j = startIndex; j < endIndex; j++) { - ret[this.columnIndices[j], i] = this.nonZeroValues[j]; + ret[this._columnIndices[j], i] = this._nonZeroValues[j]; } } return ret; } + /// + /// Copies the requested row elements into a new . + /// + /// The row to copy elements from. + /// The column to start copying from. + /// The number of elements to copy. + /// The to copy the column into. + /// If the result is . + /// If is negative, + /// or greater than or equal to the number of columns. + /// If is negative, + /// or greater than or equal to the number of rows. + /// If + + /// is greater than or equal to the number of rows. + /// If is not positive. + /// If result.Count < length. + public override void Row(int rowIndex, int columnIndex, int length, Vector result) + { + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (rowIndex >= this.RowCount || rowIndex < 0) + { + throw new ArgumentOutOfRangeException("rowIndex"); + } + + if (columnIndex >= this.ColumnCount || columnIndex < 0) + { + throw new ArgumentOutOfRangeException("columnIndex"); + } + + if (columnIndex + length > this.ColumnCount) + { + throw new ArgumentOutOfRangeException("length"); + } + + if (length < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + if (result.Count < length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); + } + + // Determin bounds in columnIndices array where this item should be searched (using rowIndex) + var startIndex = this._rowIndex[rowIndex]; + var endIndex = rowIndex < this._rowIndex.Length - 1 ? this._rowIndex[rowIndex + 1] : this.NonZerosCount; + + if (startIndex == endIndex) + { + // TODO: Maybe it is reasonable to add "Clear" method in Vector class? + // Only zero elements in rowIndex row. Clear the vector + result.Multiply(0); + } + else + { + // If there are non-zero elements use base class implementation + for (int i = columnIndex, j = 0; i < columnIndex + length; i++, j++) + { + result[j] = this.At(rowIndex, i); + } + } + } + + /// + /// Diagonally stacks this matrix on top of the given matrix and places the combined matrix into the result matrix. + /// + /// The lower, right matrix. + /// The combined matrix + /// If lower is . + /// If the result matrix is . + /// If the result matrix's dimensions are not (this.Rows + lower.rows) x (this.Columns + lower.Columns). + public override void DiagonalStack(Matrix lower, Matrix result) + { + var lowerSparseMatrix = lower as SparseMatrix; + var resultSparseMatrix = result as SparseMatrix; + + if ((lowerSparseMatrix == null) || (resultSparseMatrix == null)) + { + base.DiagonalStack(lower, result); + } + else + { + if (resultSparseMatrix.RowCount != this.RowCount + lowerSparseMatrix.RowCount || resultSparseMatrix.ColumnCount != this.ColumnCount + lowerSparseMatrix.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); + } + + resultSparseMatrix.NonZerosCount = this.NonZerosCount + lowerSparseMatrix.NonZerosCount; + resultSparseMatrix._nonZeroValues = new double[resultSparseMatrix.NonZerosCount]; + resultSparseMatrix._columnIndices = new int[resultSparseMatrix.NonZerosCount]; + + Array.Copy(this._nonZeroValues, 0, resultSparseMatrix._nonZeroValues, 0, this.NonZerosCount); + Array.Copy(lowerSparseMatrix._nonZeroValues, 0, resultSparseMatrix._nonZeroValues, this.NonZerosCount, lowerSparseMatrix.NonZerosCount); + + Array.Copy(this._columnIndices, 0, resultSparseMatrix._columnIndices, 0, this.NonZerosCount); + Array.Copy(this._rowIndex, 0, resultSparseMatrix._rowIndex, 0, this.RowCount); + + // Copy and adjust lower column indices and rowIndex + for (int i = this.NonZerosCount, j = 0; i < resultSparseMatrix.NonZerosCount; i++, j++) + { + resultSparseMatrix._columnIndices[i] = lowerSparseMatrix._columnIndices[j] + this.ColumnCount; + } + + for (int i = this.RowCount, j = 0; i < resultSparseMatrix.RowCount; i++, j++) + { + resultSparseMatrix._rowIndex[i] = lowerSparseMatrix._rowIndex[j] + this.NonZerosCount; + } + } + } + #region Elementary operations /// @@ -560,7 +1013,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double } else { - this.Add(m); + Add(m); } } @@ -585,26 +1038,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var i = 0; i < other.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = other.rowIndex[i]; - var endIndex = i < other.rowIndex.Length - 1 ? other.rowIndex[i + 1] : other.NonZerosCount; + var startIndex = other._rowIndex[i]; + var endIndex = i < other._rowIndex.Length - 1 ? other._rowIndex[i + 1] : other.NonZerosCount; for (var j = startIndex; j < endIndex; j++) { - var index = this.FindItem(i, other.columnIndices[j]); + var index = this.FindItem(i, other._columnIndices[j]); if (index >= 0) { - if (this.nonZeroValues[index] + other.nonZeroValues[j] == 0.0) + if (this._nonZeroValues[index] + other._nonZeroValues[j] == 0.0) { this.DeleteItemByIndex(index, i); } else { - this.nonZeroValues[index] += other.nonZeroValues[j]; + this._nonZeroValues[index] += other._nonZeroValues[j]; } } else { - this.SetValueAt(i, other.columnIndices[j], other.nonZeroValues[j]); + this.SetValueAt(i, other._columnIndices[j], other._nonZeroValues[j]); } } } @@ -657,26 +1110,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var i = 0; i < other.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = other.rowIndex[i]; - var endIndex = i < other.rowIndex.Length - 1 ? other.rowIndex[i + 1] : other.NonZerosCount; + var startIndex = other._rowIndex[i]; + var endIndex = i < other._rowIndex.Length - 1 ? other._rowIndex[i + 1] : other.NonZerosCount; for (var j = startIndex; j < endIndex; j++) { - var index = this.FindItem(i, other.columnIndices[j]); + var index = this.FindItem(i, other._columnIndices[j]); if (index >= 0) { - if (this.nonZeroValues[index] - other.nonZeroValues[j] == 0.0) + if (this._nonZeroValues[index] - other._nonZeroValues[j] == 0.0) { this.DeleteItemByIndex(index, i); } else { - this.nonZeroValues[index] -= other.nonZeroValues[j]; + this._nonZeroValues[index] -= other._nonZeroValues[j]; } } else { - this.SetValueAt(i, other.columnIndices[j], -other.nonZeroValues[j]); + this.SetValueAt(i, other._columnIndices[j], -other._nonZeroValues[j]); } } } @@ -699,7 +1152,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double return; } - Control.LinearAlgebraProvider.ScaleArray(scalar, this.nonZeroValues); + Control.LinearAlgebraProvider.ScaleArray(scalar, this._nonZeroValues); } /// @@ -738,18 +1191,24 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var row = 0; row < this.RowCount; row++) { // Get the begin / end index for the current row - var startIndex = this.rowIndex[row]; - var endIndex = row < this.rowIndex.Length - 1 ? this.rowIndex[row + 1] : this.NonZerosCount; + var startIndex = this._rowIndex[row]; + var endIndex = row < this._rowIndex.Length - 1 ? this._rowIndex[row + 1] : this.NonZerosCount; for (var column = 0; column < otherSparseMatrix.ColumnCount; column++) { columnVector.Clear(); otherSparseMatrix.Column(column, columnVector); - + // Multiply row of matrix A on column of matrix B - var sum = CommonParallel.Aggregate( - startIndex, - endIndex, - index => this.nonZeroValues[index] * columnVector[this.columnIndices[index]]); + var sum = 0.0; + if (startIndex != endIndex) + { + // If there are elements in that row, then calculate rowA x columnB + sum = CommonParallel.Aggregate( + startIndex, + endIndex, + index => this._nonZeroValues[index] * columnVector[this._columnIndices[index]]); + } + resultSparseMatrix.SetValueAt(row, column, sum); } } @@ -777,10 +1236,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double } var result = (SparseMatrix)this.CreateMatrix(this.RowCount, matrix.ColumnCount); - Multiply(matrix, result); + this.Multiply(matrix, result); return result; } - + /// /// Multiplies two sparse matrices. /// @@ -809,10 +1268,128 @@ namespace MathNet.Numerics.LinearAlgebra.Double return (SparseMatrix)leftSide.Multiply(rightSide); } + /// + /// Pointwise multiplies this matrix with another matrix and stores the result into the result matrix. + /// + /// The matrix to pointwise multiply with this one. + /// The matrix to store the result of the pointwise multiplication. + /// If the other matrix is . + /// If the result matrix is . + /// If this matrix and are not the same size. + /// If this matrix and are not the same size. + public override void PointwiseMultiply(Matrix other, Matrix result) + { + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (this.ColumnCount != other.ColumnCount || this.RowCount != other.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); + } + + if (this.ColumnCount != result.ColumnCount || this.RowCount != result.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); + } + + result.Clear(); + for (var i = 0; i < other.RowCount; i++) + { + // Get the begin / end index for the current row + var startIndex = this._rowIndex[i]; + var endIndex = i < this._rowIndex.Length - 1 ? this._rowIndex[i + 1] : this.NonZerosCount; + + for (var j = startIndex; j < endIndex; j++) + { + var resVal = this._nonZeroValues[j] * other[i, this._columnIndices[j]]; + if (resVal != 0.0) + { + result[i, this._columnIndices[j]] = resVal; + } + } + } + } + + /// + /// Generates matrix with random elements. + /// + /// Number of rows. + /// Number of columns. + /// Continuous Random Distribution or Source + /// + /// An numberOfRows-by-numberOfColumns matrix with elements distributed according to the provided distribution. + /// + /// If the parameter is not positive. + /// If the parameter is not positive. + public override Matrix Random(int numberOfRows, int numberOfColumns, IContinuousDistribution distribution) + { + if (numberOfRows < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfRows"); + } + + if (numberOfColumns < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfColumns"); + } + + var matrix = this.CreateMatrix(numberOfRows, numberOfColumns); + for (var i = 0; i < RowCount; i++) + { + for (var j = 0; j < matrix.ColumnCount; j++) + { + matrix[i, j] = distribution.Sample(); + } + } + + return matrix; + } + + /// + /// Generates matrix with random elements. + /// + /// Number of rows. + /// Number of columns. + /// Continuous Random Distribution or Source + /// + /// An numberOfRows-by-numberOfColumns matrix with elements distributed according to the provided distribution. + /// + /// If the parameter is not positive. + /// If the parameter is not positive. + public override Matrix Random(int numberOfRows, int numberOfColumns, IDiscreteDistribution distribution) + { + if (numberOfRows < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfRows"); + } + + if (numberOfColumns < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "numberOfColumns"); + } + + var matrix = this.CreateMatrix(numberOfRows, numberOfColumns); + for (var i = 0; i < RowCount; i++) + { + for (var j = 0; j < matrix.ColumnCount; j++) + { + matrix[i, j] = distribution.Sample(); + } + } + + return matrix; + } + #endregion #region Static constructors for special matrices. - /// /// Initializes a square with all zero's except for ones on the diagonal. /// @@ -825,21 +1402,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double { var m = new SparseMatrix(order) { - NonZerosCount = order, - nonZeroValues = new double[order], - columnIndices = new int[order] + NonZerosCount = order, + _nonZeroValues = new double[order], + _columnIndices = new int[order] }; for (var i = 0; i < order; i++) { - m.nonZeroValues[i] = 1.0; - m.columnIndices[i] = i; - m.rowIndex[i] = i; + m._nonZeroValues[i] = 1.0; + m._columnIndices[i] = i; + m._rowIndex[i] = i; } return m; } - #endregion } -} \ No newline at end of file +} diff --git a/src/Numerics/LinearAlgebra/Double/SparseVector.cs b/src/Numerics/LinearAlgebra/Double/SparseVector.cs index e324f7e4..dc5dec4b 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseVector.cs @@ -1244,7 +1244,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// If the length vector is non positive. public override Vector Random(int length, IContinuousDistribution randomDistribution) { - if (length < 0) + if (length < 1) { throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); } @@ -1270,7 +1270,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// If the n vector is non positive. public override Vector Random(int length, IDiscreteDistribution randomDistribution) { - if (length < 0) + if (length < 1) { throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); } diff --git a/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs index 38d5d9e5..6c212960 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs @@ -1,9 +1,11 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com +// // Copyright (c) 2009-2010 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -26,10 +30,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double { - using System; - using System.Collections.Generic; + using System.Collections.Generic; + using MbUnit.Framework; using LinearAlgebra.Double; - using MbUnit.Framework; public class SparseMatrixTests : MatrixTests { @@ -58,16 +61,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double { var testData = new Dictionary { - { "Singular3x3", new SparseMatrix(3, 3, new double[] { 1, 1, 1, 1, 1, 1, 2, 2, 2 }) }, - { "Square3x3", new SparseMatrix(3, 3, new[] { -1.1, 0.0, -4.4, -2.2, 1.1, 5.5, -3.3, 2.2, 6.6 }) }, - { "Square4x4", new SparseMatrix(4, 4, new[] { -1.1, 0.0, 1.0, -4.4, -2.2, 1.1, 2.1, 5.5, -3.3, 2.2, 6.2, 6.6, -4.4, 3.3, 4.3, -7.7 }) }, - { "Tall3x2", new SparseMatrix(3, 2, new[] { -1.1, 0.0, -4.4, -2.2, 1.1, 5.5 }) }, + { "Singular3x3", new SparseMatrix(3, 3, new double[] { 1, 1, 1, 1, 1, 1, 2, 2, 2 }) }, + { "Square3x3", new SparseMatrix(3, 3, new[] { -1.1, 0.0, -4.4, -2.2, 1.1, 5.5, -3.3, 2.2, 6.6 }) }, + { "Square4x4", new SparseMatrix(4, 4, new[] { -1.1, 0.0, 1.0, -4.4, -2.2, 1.1, 2.1, 5.5, -3.3, 2.2, 6.2, 6.6, -4.4, 3.3, 4.3, -7.7 }) }, + { "Tall3x2", new SparseMatrix(3, 2, new[] { -1.1, 0.0, -4.4, -2.2, 1.1, 5.5 }) }, { "Wide2x3", new SparseMatrix(2, 3, new[] { -1.1, 0.0, -2.2, 1.1, -3.3, 2.2 }) } }; foreach (var name in testData.Keys) { - Assert.AreEqual(this.testMatrices[name], testData[name]); + Assert.AreEqual(testMatrices[name], testData[name]); } } @@ -84,9 +87,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void MatrixFrom2DArrayIsCopy() { - var matrix = new SparseMatrix(this.testData2D["Singular3x3"]); + var matrix = new SparseMatrix(testData2D["Singular3x3"]); matrix[0, 0] = 10.0; - Assert.AreEqual(1.0, this.testData2D["Singular3x3"][0, 0]); + Assert.AreEqual(1.0, testData2D["Singular3x3"][0, 0]); } [Test] @@ -98,12 +101,12 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Row("Wide2x3")] public void CanCreateMatrixFrom2DArray(string name) { - var matrix = new SparseMatrix(this.testData2D[name]); - for (var i = 0; i < this.testData2D[name].GetLength(0); i++) + var matrix = new SparseMatrix(testData2D[name]); + for (var i = 0; i < testData2D[name].GetLength(0); i++) { - for (var j = 0; j < this.testData2D[name].GetLength(1); j++) + for (var j = 0; j < testData2D[name].GetLength(1); j++) { - Assert.AreEqual(this.testData2D[name][i, j], matrix[i, j]); + Assert.AreEqual(testData2D[name][i, j], matrix[i, j]); } } } @@ -148,18 +151,17 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double { var matrix = new SparseMatrix(500, 1000); var nonzero = 0; - var rnd = new Random(0); + var rnd = new System.Random(); - for (var i = 0; i < matrix.RowCount; i++) + for (int i = 0; i < matrix.RowCount; i++) { - for (var j = 0; j < matrix.ColumnCount; j++) + for (int j = 0; j < matrix.ColumnCount; j++ ) { - var value = rnd.NextDouble(); + var value = rnd.Next(10) * rnd.Next(10) * rnd.Next(10) * rnd.Next(10) * rnd.Next(10); if (value != 0) { nonzero++; } - matrix[i, j] = value; } } @@ -167,4 +169,4 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double Assert.AreEqual(matrix.NonZerosCount, nonzero); } } -} \ No newline at end of file +} diff --git a/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs b/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs index cb684e4a..cea7c1b8 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs @@ -121,10 +121,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double Assert.ForAll(vector, value => value == 5); } - // TODO: Remove [Ignore] when SparseMatrix developed [Test] [MultipleAsserts] - [Ignore] public void CanCreateSparseMatrix() { var vector = new SparseVector(3); diff --git a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.cs b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.cs index 93212f72..53adeb75 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.cs @@ -217,7 +217,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double Assert.AreElementsEqual(vector, array); } - [Test, Ignore] + [Test] [MultipleAsserts] public void CanConvertVectorToColumnMatrix() { @@ -233,7 +233,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double } } - [Test, Ignore] + [Test] [MultipleAsserts] public void CanConvertVectorToRowMatrix() {