diff --git a/src/MathNet.Numerics.5.1.ReSharper b/src/MathNet.Numerics.5.1.ReSharper index 46415860..af81e577 100644 --- a/src/MathNet.Numerics.5.1.ReSharper +++ b/src/MathNet.Numerics.5.1.ReSharper @@ -15,7 +15,10 @@ cancelled pre preconditioner Preconditioners -indetermined +indetermined +indices +Pointwise +Frobenius @@ -561,7 +564,7 @@ indetermined - false + true Any @@ -593,7 +596,7 @@ indetermined - false + true Any @@ -657,7 +660,7 @@ indetermined - false + true Any @@ -686,7 +689,7 @@ indetermined - false + true Any diff --git a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs index b743b1a4..630c778c 100644 --- a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs @@ -1334,5 +1334,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex index => otherVector.Data[index] = Data[index].Conjugate()); } } + + /// Gets the value at the given . + /// The index of the value to get or set. + /// The value of the vector at the given . + internal protected override Complex At(int index) + { + return Data[index]; + } + + /// Sets the at the given . + /// The index of the value to get or set. + /// The value to set. + internal protected override void At(int index, Complex value) + { + Data[index] = value; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs index ddcd222e..88b945d2 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs @@ -194,7 +194,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// The matrix to copy. public SparseMatrix(Matrix matrix) - : base(matrix.RowCount, matrix.ColumnCount) + : this(matrix.RowCount, matrix.ColumnCount) { var sparseMatrix = matrix as SparseMatrix; diff --git a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs index 551ca0b4..64146344 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs @@ -38,13 +38,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// A vector with sparse storage. /// + /// The sparse vector is not thread safe. public class SparseVector : Vector { - /// - /// Lock object for the indexer. - /// - private readonly object _lockObject = new object(); - /// /// Gets the vector's internal data. The array containing the actual values; only the non-zero values are stored. /// @@ -130,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex { for (var i = 0; i < other.Count; i++) { - this[i] = other[i]; + this[i] = other.At(i); } } else @@ -174,9 +170,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// The array to create this vector from. /// The vector copy the array. Any changes to the vector will NOT change the array. - public SparseVector(Complex[] array) : this(array.Length) + public SparseVector(IList array) : this(array.Count) { - for (var i = 0; i < array.Length; i++) + for (var i = 0; i < array.Count; i++) { this[i] = array[i]; } @@ -193,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var matrix = new SparseMatrix(Count, 1); for (var i = 0; i < NonZerosCount; i++) { - matrix[_nonZeroIndices[i], 0] = _nonZeroValues[i]; + matrix.At(_nonZeroIndices[i], 0, _nonZeroValues[i]); } return matrix; @@ -208,55 +204,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var matrix = new SparseMatrix(1, Count); for (var i = 0; i < NonZerosCount; i++) { - matrix[0, _nonZeroIndices[i]] = _nonZeroValues[i]; + matrix.At(0, _nonZeroIndices[i], _nonZeroValues[i]); } return matrix; } - /// Gets or sets the value at the given . - /// The index of the value to get or set. - /// The value of the vector at the given . - /// If is negative or - /// greater than the size of the vector. - public override Complex this[int index] - { - get - { - // If index is out of bounds - if ((index < 0) || (index >= Count)) - { - throw new IndexOutOfRangeException(); - } - - lock (_lockObject) - { - // Search if item idex exists in NonZeroIndices array in range "0 - complex nonzero values count" - var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); - if (itemIndex >= 0) - { - return _nonZeroValues[itemIndex]; - } - } - - return Complex.Zero; - } - - set - { - // If index is out of bounds - if ((index < 0) || (index >= Count)) - { - throw new IndexOutOfRangeException(); - } - - lock (_lockObject) - { - SetValue(index, value); - } - } - } - /// /// Creates a matrix with the given dimensions using the same storage type /// as this vector. @@ -330,10 +283,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var otherVector = target as SparseVector; if (otherVector == null) { - CommonParallel.For( - 0, - Count, - index => target[index] = this[index]); + target.Clear(); + + for (var index = 0; index < NonZerosCount; index++) + { + target.At(_nonZeroIndices[index], _nonZeroValues[index]); + } } else { @@ -392,147 +347,42 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } } } + #region Operators and supplementary functions /// - /// Adds a complex to each element of the vector. - /// - /// The complex to add. - /// A copy of the vector with the complex added. - public override Vector Add(Complex complex) - { - if (complex == Complex.Zero) - { - return Clone(); - } - - var copy = (SparseVector)Clone(); - for (var i = 0; i < Count; i++) - { - copy[i] += complex; - } - - return copy; - } - - /// - /// Adds a complex to each element of the vector and stores the result in the result vector. - /// - /// The complex to add. - /// The vector to store the result of the addition. - /// If the result vector is . - /// If this vector and are not the same size. - public override void Add(Complex complex, Vector result) - { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - base.Add(complex, result); - } - - /// - /// Adds another vector to this vector. + /// Adds a scalar to each element of the vector and stores the result in the result vector. /// - /// The vector to add to this one. - /// A new vector containing the sum of both vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector Add(Vector other) - { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - var sparseVector = other as SparseVector; - - if (sparseVector == null) - { - return base.Add(other); - } - - var copy = (SparseVector)Clone(); - copy.AddScaledSparseVector(Complex.One, sparseVector); - return copy; - } - - /// - /// Adds the scaled sparse vector. - /// - /// The alpha. - /// The other. - private void AddScaledSparseVector(Complex alpha, SparseVector other) + /// + /// The scalar to add. + /// + /// + /// The vector to store the result of the addition. + /// + protected override void DoAdd(Complex scalar, Vector result) { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) + if (scalar == Complex.Zero) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } + if (!ReferenceEquals(this, result)) + { + CopyTo(result); + } - if (alpha == Complex.Zero) - { return; } - // I don't use ILinearAlgebraProvider because we will get no benefit due to "lock" in this[index] - // Possible fucniton in ILinearAlgebraProvider may be AddSparseVectorToScaledSparseVector(T[] y, int[] yIndices, T alpha, T[] x, int[] xIndices); - // But it require to develop value setting algorithm and due to "lock" it will be even more greedy then implemented below - if (ReferenceEquals(this, other)) + if (ReferenceEquals(this, result)) { - // Adding the same instance of sparse vector. That means if we modify "this" then "other" will be modified too. - // To avoid such problem lets change values in internal storage of "this" - if (alpha == Complex.One) - { - for (var i = 0; i < NonZerosCount; i++) - { - _nonZeroValues[i] += _nonZeroValues[i]; - } - } - else if (alpha == -Complex.One) - { - Clear(); // Vector is subtracted from itself - return; - } - else - { - for (var i = 0; i < NonZerosCount; i++) - { - _nonZeroValues[i] += alpha * _nonZeroValues[i]; - } - } + CommonParallel.For( + 0, + NonZerosCount, + index => _nonZeroValues[index] += scalar); } else { - // "this" and "other" are different objects, so by modifying "this" the "other" object will not be changed - if (alpha == Complex.One) + for (var index = 0; index < Count; index++) { - for (var i = 0; i < other.NonZerosCount; i++) - { - this[other._nonZeroIndices[i]] += other._nonZeroValues[i]; - } - } - else - { - for (var i = 0; i < other.NonZerosCount; i++) - { - this[other._nonZeroIndices[i]] += alpha * other._nonZeroValues[i]; - } + result.At(index, At(index) + scalar); } } } @@ -540,53 +390,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// Adds another vector to this vector and stores the result into the result vector. /// - /// The vector to add to this one. - /// The vector to store the result of the addition. - /// If the other vector is . - /// If the result vector is . - /// If this vector and are not the same size. - /// If this vector and are not the same size. - public override void Add(Vector other, Vector result) + /// + /// The vector to add to this one. + /// + /// + /// The vector to store the result of the addition. + /// + protected override void DoAdd(Vector other, Vector result) { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) + if (ReferenceEquals(this, result)) { - var tmp = Add(other); - tmp.CopyTo(result); + CommonParallel.For( + 0, + NonZerosCount, + index => _nonZeroValues[index] += _nonZeroValues[index]); } else { - var sparse = result as SparseVector; - if (sparse == null) + for (var index = 0; index < Count; index++) { - base.Add(other, result); - } - else - { - var sparseother = other as SparseVector; - if (sparseother == null) - { - base.Add(other, result); - } - else - { - CopyTo(result); - sparse.AddScaledSparseVector(Complex.One, sparseother); - } + result.At(index, At(index) + other.At(index)); } } } @@ -637,130 +460,39 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } /// - /// Subtracts a complex from each element of the vector. - /// - /// The complex to subtract. - /// A new vector containing the subtraction of this vector and the complex. - public override Vector Subtract(Complex complex) - { - if (complex == Complex.Zero) - { - return Clone(); - } - - var copy = (SparseVector)Clone(); - for (var i = 0; i < Count; i++) - { - copy[i] -= complex; - } - - return copy; - } - - /// - /// Subtracts a complex from each element of the vector and stores the result in the result vector. - /// - /// The complex to subtract. - /// The vector to store the result of the subtraction. - /// If the result vector is . - /// If this vector and are not the same size. - public override void Subtract(Complex complex, Vector result) - { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - base.Subtract(complex, result); - } - - /// - /// Subtracts another vector from this vector. + /// Subtracts a scalar from each element of the vector and stores the result in the result vector. /// - /// The vector to subtract from this one. - /// A new vector containing the subtraction of the the two vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector Subtract(Vector other) + /// + /// The scalar to subtract. + /// + /// + /// The vector to store the result of the subtraction. + /// + protected override void DoSubtract(Complex scalar, Vector result) { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - var sparseVector = other as SparseVector; - - if (sparseVector == null) - { - return base.Subtract(other); - } - - var copy = (SparseVector)Clone(); - copy.AddScaledSparseVector(-Complex.One, sparseVector); - return copy; + DoAdd(-scalar, result); } /// /// Subtracts another vector to this vector and stores the result into the result vector. /// - /// The vector to subtract from this one. - /// The vector to store the result of the subtraction. - /// If the other vector is . - /// If the result vector is . - /// If this vector and are not the same size. - /// If this vector and are not the same size. - public override void Subtract(Vector other, Vector result) + /// + /// The vector to subtract from this one. + /// + /// + /// The vector to store the result of the subtraction. + /// + protected override void DoSubtract(Vector other, Vector result) { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - if (Count != result.Count) + if (ReferenceEquals(this, other)) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); + result.Clear(); + return; } - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) - { - var tmp = Subtract(other); - tmp.CopyTo(result); - } - else + for (var index = 0; index < Count; index++) { - var sparse = result as SparseVector; - if (sparse == null) - { - base.Subtract(other, result); - } - else - { - var sparseother = other as SparseVector; - if (sparseother == null) - { - base.Subtract(other, result); - } - else - { - CopyTo(result); - sparse.AddScaledSparseVector(-Complex.One, sparseother); - } - } + result.At(index, At(index) - other.At(index)); } } @@ -835,52 +567,81 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } /// - /// Multiplies a complex to each element of the vector. + /// Multiplies a scalar to each element of the vector and stores the result in the result vector. /// - /// The complex to multiply. - /// A new vector that is the multiplication of the vector and the complex. - public override Vector Multiply(Complex complex) + /// + /// The scalar to multiply. + /// + /// + /// The vector to store the result of the multiplication. + /// + protected override void DoMultiply(Complex scalar, Vector result) { - if (complex == Complex.One) + if (scalar == Complex.One) { - return Clone(); + if (!ReferenceEquals(this, result)) + { + CopyTo(result); + } + + return; } - if (complex == Complex.Zero) + if (scalar == Complex.Zero) { - return new SparseVector(Count); + result.Clear(); + return; } - var copy = new SparseVector(this); - Control.LinearAlgebraProvider.ScaleArray(complex, copy._nonZeroValues, copy._nonZeroValues); - return copy; + var sparseResult = result as SparseVector; + if (sparseResult == null) + { + result.Clear(); + for (var index = 0; index < NonZerosCount; index++) + { + result.At(_nonZeroIndices[index], scalar * _nonZeroValues[index]); + } + } + else + { + if (!ReferenceEquals(this, result)) + { + sparseResult.NonZerosCount = NonZerosCount; + sparseResult._nonZeroIndices = new int[NonZerosCount]; + Buffer.BlockCopy(_nonZeroIndices, 0, sparseResult._nonZeroIndices, 0, _nonZeroIndices.Length * Constants.SizeOfInt); + sparseResult._nonZeroValues = new Complex[_nonZeroValues.Length]; + } + + Control.LinearAlgebraProvider.ScaleArray(scalar, _nonZeroValues, sparseResult._nonZeroValues); + } } /// /// Computes the dot product between this vector and another vector. /// - /// The other vector to add. - /// The result of the addition. - /// If is not of the same size. - /// If is . - public override Complex DotProduct(Vector other) + /// + /// The other vector to add. + /// + /// s + /// The result of the addition. + /// + protected override Complex DoDotProduct(Vector other) { - if (other == null) - { - throw new ArgumentNullException("other"); - } + var result = Complex.Zero; - if (Count != other.Count) + if (ReferenceEquals(this, other)) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + for (var i = 0; i < NonZerosCount; i++) + { + result += _nonZeroValues[i] * _nonZeroValues[i]; + } } - - var result = Complex.Zero; - - // base implementation iterates though all elements, but we need only take non-zeros - for (var i = 0; i < NonZerosCount; i++) + else { - result += _nonZeroValues[i] * other[_nonZeroIndices[i]]; + for (var i = 0; i < NonZerosCount; i++) + { + result += _nonZeroValues[i] * other.At(_nonZeroIndices[i]); + } } return result; @@ -1023,7 +784,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var result = new SparseVector(length); for (var i = index; i < index + length; i++) { - result[i - index] = this[i]; + result.At(i - index, At(i)); } return result; @@ -1049,7 +810,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex for (var i = 0; i < values.Length; i++) { - this[i] = values[i]; + At(i, values[i]); } } @@ -1084,35 +845,51 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } /// - /// Pointwise multiplies this vector with another vector. + /// Pointwise multiplies this vector with another vector and stores the result into the result vector. /// /// The vector to pointwise multiply with this one. - /// A new vector which is the pointwise multiplication of the two vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector PointwiseMultiply(Vector other) + /// The vector to store the result of the pointwise multiplication. + protected override void DoPointwiseMultiply(Vector other, Vector result) { - if (other == null) + if (ReferenceEquals(this, other)) { - throw new ArgumentNullException("other"); + for (var i = 0; i < NonZerosCount; i++) + { + _nonZeroValues[i] *= _nonZeroValues[i]; + } } - - if (Count != other.Count) + else { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + for (var i = 0; i < NonZerosCount; i++) + { + var index = _nonZeroIndices[i]; + result.At(index, other.At(index) * _nonZeroValues[i]); + } } + } - var copy = new SparseVector(Count); - for (var i = 0; i < _nonZeroIndices.Length; i++) + /// + /// Pointwise multiplies this vector with another vector and stores the result into the result vector. + /// + /// The vector to pointwise multiply with this one. + /// The vector to store the result of the pointwise multiplication. + protected override void DoPointwiseDivide(Vector other, Vector result) + { + if (ReferenceEquals(this, other)) { - var d = _nonZeroValues[i] * other[_nonZeroIndices[i]]; - if (d != Complex.Zero) + for (var i = 0; i < NonZerosCount; i++) { - copy[_nonZeroIndices[i]] = d; + _nonZeroValues[i] /= _nonZeroValues[i]; + } + } + else + { + for (var i = 0; i < NonZerosCount; i++) + { + var index = _nonZeroIndices[i]; + result.At(index, _nonZeroValues[i] / other.At(index)); } } - - return copy; } /// @@ -1123,7 +900,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// Matrix M[i,j] = u[i]*v[j] /// If the u vector is . /// If the v vector is . - public static Matrix /*SparseMatrix*/ OuterProduct(SparseVector u, SparseVector v) + public static Matrix OuterProduct(SparseVector u, SparseVector v) { if (u == null) { @@ -1361,13 +1138,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } #endregion + /// + /// Gets the value at the given index. + /// + /// Value real index in array + /// The value at the given index. + internal protected override Complex At(int index) + { + // Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count" + var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); + return itemIndex >= 0 ? _nonZeroValues[itemIndex] : Complex.Zero; + } + /// /// Delete, Add or Update the value in NonZeroValues and NonZeroIndices /// /// Value real index in array /// The value to set. /// This method assume that index is between 0 and Array Size - private void SetValue(int index, Complex value) + internal protected override void At(int index, Complex value) { // Search if "index" already exists in range "0 - complex nonzero values count" var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); @@ -1549,5 +1338,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex yield return new KeyValuePair(_nonZeroIndices[i], _nonZeroValues[i]); } } + + /// + /// Returns the data contained in the vector as an array. + /// + /// + /// The vector's data as an array. + /// + public override Complex[] ToArray() + { + var ret = new Complex[Count]; + for (var i = 0; i < NonZerosCount; i++) + { + ret[_nonZeroIndices[i]] = _nonZeroValues[i]; + } + + return ret; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Vector.cs b/src/Numerics/LinearAlgebra/Complex/Vector.cs index e0440e25..39c3ca9f 100644 --- a/src/Numerics/LinearAlgebra/Complex/Vector.cs +++ b/src/Numerics/LinearAlgebra/Complex/Vector.cs @@ -63,10 +63,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoAdd(Complex scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] + scalar); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) + scalar); + } } /// @@ -80,10 +80,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoAdd(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] + other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) + other.At(index)); + } } /// @@ -111,10 +111,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoSubtract(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] - other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) - other.At(index)); + } } /// @@ -128,10 +128,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoMultiply(Complex scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] * scalar); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) * scalar); + } } /// @@ -145,10 +145,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoDivide(Complex scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] / scalar); + DoMultiply(1 / scalar, result); } /// @@ -158,10 +155,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The vector to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] * other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) * other.At(index)); + } } /// @@ -171,10 +168,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The vector to store the result of the pointwise division. protected override void DoPointwiseDivide(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] / other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) / other.At(index)); + } } /// @@ -183,7 +180,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// The other vector to add. /// - /// s + /// /// The result of the addition. /// protected override Complex DoDotProduct(Vector other) @@ -191,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return CommonParallel.Aggregate( 0, Count, - i => this[i] * other[i]); + i => At(i) * other.At(i)); } /// @@ -200,7 +197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The value of the absolute minimum element. public override Complex AbsoluteMinimum() { - return this[AbsoluteMinimumIndex()].Magnitude; + return At(AbsoluteMinimumIndex()).Magnitude; } /// @@ -210,10 +207,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public override int AbsoluteMinimumIndex() { var index = 0; - var min = this[index].Magnitude; + var min = At(index).Magnitude; for (var i = 1; i < Count; i++) { - var test = this[i].Magnitude; + var test = At(i).Magnitude; if (test < min) { index = i; @@ -230,7 +227,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The value of the absolute maximum element. public override Complex AbsoluteMaximum() { - return this[AbsoluteMaximumIndex()].Magnitude; + return At(AbsoluteMaximumIndex()).Magnitude; } /// @@ -240,10 +237,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public override int AbsoluteMaximumIndex() { var index = 0; - var max = this[index].Magnitude; + var max = At(index).Magnitude; for (var i = 1; i < Count; i++) { - var test = this[i].Magnitude; + var test = At(i).Magnitude; if (test > max) { index = i; @@ -263,7 +260,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return CommonParallel.Aggregate( 0, Count, - i => this[i]); + i => At(i)); } /// @@ -275,7 +272,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return CommonParallel.Aggregate( 0, Count, - i => this[i].Magnitude); + i => At(i).Magnitude); } /// @@ -285,7 +282,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The p value. /// /// - /// Scalar ret = (sum(abs(this[i])^p))^(1/p) + /// Scalar ret = (sum(abs(At(i))^p))^(1/p) /// public override Complex Norm(double p) { @@ -299,14 +296,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return CommonParallel.Select( 0, Count, - (index, localData) => Math.Max(localData, this[index].Magnitude), + (index, localData) => Math.Max(localData, At(index).Magnitude), Math.Max); } var sum = CommonParallel.Aggregate( 0, Count, - index => Math.Pow(this[index].Magnitude, p)); + index => Math.Pow(At(index).Magnitude, p)); return Math.Pow(sum, 1.0 / p); } @@ -317,11 +314,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// Target vector protected override void DoConjugate(Vector target) { - CopyTo(target); - CommonParallel.For( - 0, - Count, - index => target[index] = this[index].Conjugate()); + for (var index = 0; index < Count; index++) + { + target.At(index, At(index).Conjugate()); + } } /// @@ -336,10 +332,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public override Vector Negate() { var result = CreateVector(Count); - CommonParallel.For( - 0, - Count, - index => result[index] = -this[index]); + + for (var index = 0; index < Count; index++) + { + result.At(index, -At(index)); + } return result; } @@ -364,7 +361,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var vector = CreateVector(length); for (var index = 0; index < length; index++) { - vector[index] = new Complex(randomDistribution.Sample(), randomDistribution.Sample()); + vector.At(index, new Complex(randomDistribution.Sample(), randomDistribution.Sample())); } return vector; @@ -390,7 +387,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var vector = CreateVector(length); for (var index = 0; index < length; index++) { - vector[index] = new Complex(randomDistribution.Sample(), randomDistribution.Sample()); + vector.At(index, new Complex(randomDistribution.Sample(), randomDistribution.Sample())); } return vector; diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs index 2214e708..158f6888 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs @@ -1387,5 +1387,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 index => otherVector.Data[index] = Data[index].Conjugate()); } } + + /// Gets the value at the given . + /// The index of the value to get or set. + /// The value of the vector at the given . + internal protected override Complex32 At(int index) + { + return Data[index]; + } + + /// Sets the at the given . + /// The index of the value to get or set. + /// The value to set. + internal protected override void At(int index, Complex32 value) + { + Data[index] = value; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs index 584a13e4..ee35a380 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs @@ -189,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// The matrix to copy. public SparseMatrix(Matrix matrix) - : base(matrix.RowCount, matrix.ColumnCount) + : this(matrix.RowCount, matrix.ColumnCount) { var sparseMatrix = matrix as SparseMatrix; diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs index 3ceedb9f..991f32c8 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs @@ -38,13 +38,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// A vector with sparse storage. /// + /// The sparse vector is not thread safe. public class SparseVector : Vector { - /// - /// Lock object for the indexer. - /// - private readonly object _lockObject = new object(); - /// /// Gets the vector's internal data. The array containing the actual values; only the non-zero values are stored. /// @@ -130,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 { for (var i = 0; i < other.Count; i++) { - this[i] = other[i]; + this[i] = other.At(i); } } else @@ -174,9 +170,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// The array to create this vector from. /// The vector copy the array. Any changes to the vector will NOT change the array. - public SparseVector(Complex32[] array) : this(array.Length) + public SparseVector(IList array) : this(array.Count) { - for (var i = 0; i < array.Length; i++) + for (var i = 0; i < array.Count; i++) { this[i] = array[i]; } @@ -193,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var matrix = new SparseMatrix(Count, 1); for (var i = 0; i < NonZerosCount; i++) { - matrix[_nonZeroIndices[i], 0] = _nonZeroValues[i]; + matrix.At(_nonZeroIndices[i], 0, _nonZeroValues[i]); } return matrix; @@ -208,7 +204,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var matrix = new SparseMatrix(1, Count); for (var i = 0; i < NonZerosCount; i++) { - matrix[0, _nonZeroIndices[i]] = _nonZeroValues[i]; + matrix.At(0, _nonZeroIndices[i], _nonZeroValues[i]); } return matrix; @@ -229,17 +225,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 throw new IndexOutOfRangeException(); } - lock (_lockObject) - { - // Search if item idex exists in NonZeroIndices array in range "0 - complex nonzero values count" - var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); - if (itemIndex >= 0) - { - return _nonZeroValues[itemIndex]; - } - } - - return Complex32.Zero; + return At(index); } set @@ -250,10 +236,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 throw new IndexOutOfRangeException(); } - lock (_lockObject) - { - SetValue(index, value); - } + At(index, value); } } @@ -330,10 +313,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var otherVector = target as SparseVector; if (otherVector == null) { - CommonParallel.For( - 0, - Count, - index => target[index] = this[index]); + target.Clear(); + + for (var index = 0; index < NonZerosCount; index++) + { + target.At(_nonZeroIndices[index], _nonZeroValues[index]); + } } else { @@ -392,147 +377,42 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } } } + #region Operators and supplementary functions /// - /// Adds a complex to each element of the vector. + /// Adds a scalar to each element of the vector and stores the result in the result vector. /// - /// The complex to add. - /// A copy of the vector with the complex added. - public override Vector Add(Complex32 complex) - { - if (complex == Complex32.Zero) - { - return Clone(); - } - - var copy = (SparseVector)Clone(); - for (var i = 0; i < Count; i++) - { - copy[i] += complex; - } - - return copy; - } - - /// - /// Adds a complex to each element of the vector and stores the result in the result vector. - /// - /// The complex to add. - /// The vector to store the result of the addition. - /// If the result vector is . - /// If this vector and are not the same size. - public override void Add(Complex32 complex, Vector result) - { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - base.Add(complex, result); - } - - /// - /// Adds another vector to this vector. - /// - /// The vector to add to this one. - /// A new vector containing the sum of both vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector Add(Vector other) - { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - var sparseVector = other as SparseVector; - - if (sparseVector == null) - { - return base.Add(other); - } - - var copy = (SparseVector)Clone(); - copy.AddScaledSparseVector(Complex32.One, sparseVector); - return copy; - } - - /// - /// Adds the scaled sparse vector. - /// - /// The alpha. - /// The other. - private void AddScaledSparseVector(Complex32 alpha, SparseVector other) + /// + /// The scalar to add. + /// + /// + /// The vector to store the result of the addition. + /// + protected override void DoAdd(Complex32 scalar, Vector result) { - if (other == null) + if (scalar == Complex32.Zero) { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } + if (!ReferenceEquals(this, result)) + { + CopyTo(result); + } - if (alpha == Complex32.Zero) - { return; } - // I don't use ILinearAlgebraProvider because we will get no benefit due to "lock" in this[index] - // Possible fucniton in ILinearAlgebraProvider may be AddSparseVectorToScaledSparseVector(T[] y, int[] yIndices, T alpha, T[] x, int[] xIndices); - // But it require to develop value setting algorithm and due to "lock" it will be even more greedy then implemented below - if (ReferenceEquals(this, other)) + if (ReferenceEquals(this, result)) { - // Adding the same instance of sparse vector. That means if we modify "this" then "other" will be modified too. - // To avoid such problem lets change values in internal storage of "this" - if (alpha == Complex32.One) - { - for (var i = 0; i < NonZerosCount; i++) - { - _nonZeroValues[i] += _nonZeroValues[i]; - } - } - else if (alpha == -Complex32.One) - { - Clear(); // Vector is subtracted from itself - return; - } - else - { - for (var i = 0; i < NonZerosCount; i++) - { - _nonZeroValues[i] += alpha * _nonZeroValues[i]; - } - } + CommonParallel.For( + 0, + NonZerosCount, + index => _nonZeroValues[index] += scalar); } else { - // "this" and "other" are different objects, so by modifying "this" the "other" object will not be changed - if (alpha == Complex32.One) + for (var index = 0; index < Count; index++) { - for (var i = 0; i < other.NonZerosCount; i++) - { - this[other._nonZeroIndices[i]] += other._nonZeroValues[i]; - } - } - else - { - for (var i = 0; i < other.NonZerosCount; i++) - { - this[other._nonZeroIndices[i]] += alpha * other._nonZeroValues[i]; - } + result.At(index, At(index) + scalar); } } } @@ -540,53 +420,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// Adds another vector to this vector and stores the result into the result vector. /// - /// The vector to add to this one. - /// The vector to store the result of the addition. - /// If the other vector is . - /// If the result vector is . - /// If this vector and are not the same size. - /// If this vector and are not the same size. - public override void Add(Vector other, Vector result) + /// + /// The vector to add to this one. + /// + /// + /// The vector to store the result of the addition. + /// + protected override void DoAdd(Vector other, Vector result) { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) + if (ReferenceEquals(this, result)) { - var tmp = Add(other); - tmp.CopyTo(result); + CommonParallel.For( + 0, + NonZerosCount, + index => _nonZeroValues[index] += _nonZeroValues[index]); } else { - var sparse = result as SparseVector; - if (sparse == null) + for (var index = 0; index < Count; index++) { - base.Add(other, result); - } - else - { - var sparseother = other as SparseVector; - if (sparseother == null) - { - base.Add(other, result); - } - else - { - CopyTo(result); - sparse.AddScaledSparseVector(Complex32.One, sparseother); - } + result.At(index, At(index) + other.At(index)); } } } @@ -637,130 +490,39 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } /// - /// Subtracts a complex from each element of the vector. - /// - /// The complex to subtract. - /// A new vector containing the subtraction of this vector and the complex. - public override Vector Subtract(Complex32 complex) - { - if (complex == Complex32.Zero) - { - return Clone(); - } - - var copy = (SparseVector)Clone(); - for (var i = 0; i < Count; i++) - { - copy[i] -= complex; - } - - return copy; - } - - /// - /// Subtracts a complex from each element of the vector and stores the result in the result vector. - /// - /// The complex to subtract. - /// The vector to store the result of the subtraction. - /// If the result vector is . - /// If this vector and are not the same size. - public override void Subtract(Complex32 complex, Vector result) - { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - base.Subtract(complex, result); - } - - /// - /// Subtracts another vector from this vector. + /// Subtracts a scalar from each element of the vector and stores the result in the result vector. /// - /// The vector to subtract from this one. - /// A new vector containing the subtraction of the the two vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector Subtract(Vector other) + /// + /// The scalar to subtract. + /// + /// + /// The vector to store the result of the subtraction. + /// + protected override void DoSubtract(Complex32 scalar, Vector result) { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - var sparseVector = other as SparseVector; - - if (sparseVector == null) - { - return base.Subtract(other); - } - - var copy = (SparseVector)Clone(); - copy.AddScaledSparseVector(-Complex32.One, sparseVector); - return copy; + DoAdd(-scalar, result); } /// /// Subtracts another vector to this vector and stores the result into the result vector. /// - /// The vector to subtract from this one. - /// The vector to store the result of the subtraction. - /// If the other vector is . - /// If the result vector is . - /// If this vector and are not the same size. - /// If this vector and are not the same size. - public override void Subtract(Vector other, Vector result) + /// + /// The vector to subtract from this one. + /// + /// + /// The vector to store the result of the subtraction. + /// + protected override void DoSubtract(Vector other, Vector result) { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - if (Count != result.Count) + if (ReferenceEquals(this, other)) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); + result.Clear(); + return; } - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) - { - var tmp = Subtract(other); - tmp.CopyTo(result); - } - else + for (var index = 0; index < Count; index++) { - var sparse = result as SparseVector; - if (sparse == null) - { - base.Subtract(other, result); - } - else - { - var sparseother = other as SparseVector; - if (sparseother == null) - { - base.Subtract(other, result); - } - else - { - CopyTo(result); - sparse.AddScaledSparseVector(-Complex32.One, sparseother); - } - } + result.At(index, At(index) - other.At(index)); } } @@ -835,52 +597,81 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } /// - /// Multiplies a complex to each element of the vector. + /// Multiplies a scalar to each element of the vector and stores the result in the result vector. /// - /// The complex to multiply. - /// A new vector that is the multiplication of the vector and the complex. - public override Vector Multiply(Complex32 complex) + /// + /// The scalar to multiply. + /// + /// + /// The vector to store the result of the multiplication. + /// + protected override void DoMultiply(Complex32 scalar, Vector result) { - if (complex == Complex32.One) + if (scalar == Complex32.One) { - return Clone(); + if (!ReferenceEquals(this, result)) + { + CopyTo(result); + } + + return; } - if (complex == Complex32.Zero) + if (scalar == Complex32.Zero) { - return new SparseVector(Count); + result.Clear(); + return; } - var copy = new SparseVector(this); - Control.LinearAlgebraProvider.ScaleArray(complex, copy._nonZeroValues, copy._nonZeroValues); - return copy; + var sparseResult = result as SparseVector; + if (sparseResult == null) + { + result.Clear(); + for (var index = 0; index < NonZerosCount; index++) + { + result.At(_nonZeroIndices[index], scalar * _nonZeroValues[index]); + } + } + else + { + if (!ReferenceEquals(this, result)) + { + sparseResult.NonZerosCount = NonZerosCount; + sparseResult._nonZeroIndices = new int[NonZerosCount]; + Buffer.BlockCopy(_nonZeroIndices, 0, sparseResult._nonZeroIndices, 0, _nonZeroIndices.Length * Constants.SizeOfInt); + sparseResult._nonZeroValues = new Complex32[_nonZeroValues.Length]; + } + + Control.LinearAlgebraProvider.ScaleArray(scalar, _nonZeroValues, sparseResult._nonZeroValues); + } } /// /// Computes the dot product between this vector and another vector. /// - /// The other vector to add. - /// The result of the addition. - /// If is not of the same size. - /// If is . - public override Complex32 DotProduct(Vector other) + /// + /// The other vector to add. + /// + /// s + /// The result of the addition. + /// + protected override Complex32 DoDotProduct(Vector other) { - if (other == null) - { - throw new ArgumentNullException("other"); - } + var result = Complex32.Zero; - if (Count != other.Count) + if (ReferenceEquals(this, other)) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + for (var i = 0; i < NonZerosCount; i++) + { + result += _nonZeroValues[i] * _nonZeroValues[i]; + } } - - var result = Complex32.Zero; - - // base implementation iterates though all elements, but we need only take non-zeros - for (var i = 0; i < NonZerosCount; i++) + else { - result += _nonZeroValues[i] * other[_nonZeroIndices[i]]; + for (var i = 0; i < NonZerosCount; i++) + { + result += _nonZeroValues[i] * other.At(_nonZeroIndices[i]); + } } return result; @@ -1023,7 +814,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var result = new SparseVector(length); for (var i = index; i < index + length; i++) { - result[i - index] = this[i]; + result.At(i - index, At(i)); } return result; @@ -1049,7 +840,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 for (var i = 0; i < values.Length; i++) { - this[i] = values[i]; + At(i, values[i]); } } @@ -1084,35 +875,51 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } /// - /// Pointwise multiplies this vector with another vector. + /// Pointwise multiplies this vector with another vector and stores the result into the result vector. /// /// The vector to pointwise multiply with this one. - /// A new vector which is the pointwise multiplication of the two vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector PointwiseMultiply(Vector other) + /// The vector to store the result of the pointwise multiplication. + protected override void DoPointwiseMultiply(Vector other, Vector result) { - if (other == null) + if (ReferenceEquals(this, other)) { - throw new ArgumentNullException("other"); + for (var i = 0; i < NonZerosCount; i++) + { + _nonZeroValues[i] *= _nonZeroValues[i]; + } } - - if (Count != other.Count) + else { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + for (var i = 0; i < NonZerosCount; i++) + { + var index = _nonZeroIndices[i]; + result.At(index, other.At(index) * _nonZeroValues[i]); + } } + } - var copy = new SparseVector(Count); - for (var i = 0; i < _nonZeroIndices.Length; i++) + /// + /// Pointwise multiplies this vector with another vector and stores the result into the result vector. + /// + /// The vector to pointwise multiply with this one. + /// The vector to store the result of the pointwise multiplication. + protected override void DoPointwiseDivide(Vector other, Vector result) + { + if (ReferenceEquals(this, other)) { - var d = _nonZeroValues[i] * other[_nonZeroIndices[i]]; - if (d != Complex32.Zero) + for (var i = 0; i < NonZerosCount; i++) { - copy[_nonZeroIndices[i]] = d; + _nonZeroValues[i] /= _nonZeroValues[i]; + } + } + else + { + for (var i = 0; i < NonZerosCount; i++) + { + var index = _nonZeroIndices[i]; + result.At(index, _nonZeroValues[i] / other.At(index)); } } - - return copy; } /// @@ -1361,13 +1168,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } #endregion + /// + /// Gets the value at the given index. + /// + /// Value real index in array + /// The value at the given index. + internal protected override Complex32 At(int index) + { + // Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count" + var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); + return itemIndex >= 0 ? _nonZeroValues[itemIndex] : Complex32.Zero; + } + /// /// Delete, Add or Update the value in NonZeroValues and NonZeroIndices /// /// Value real index in array /// The value to set. /// This method assume that index is between 0 and Array Size - private void SetValue(int index, Complex32 value) + internal protected override void At(int index, Complex32 value) { // Search if "index" already exists in range "0 - complex nonzero values count" var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); @@ -1549,5 +1368,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 yield return new KeyValuePair(_nonZeroIndices[i], _nonZeroValues[i]); } } + + /// + /// Returns the data contained in the vector as an array. + /// + /// + /// The vector's data as an array. + /// + public override Complex32[] ToArray() + { + var ret = new Complex32[Count]; + for (var i = 0; i < NonZerosCount; i++) + { + ret[_nonZeroIndices[i]] = _nonZeroValues[i]; + } + + return ret; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Vector.cs b/src/Numerics/LinearAlgebra/Complex32/Vector.cs index ea0f3dc7..9db2acf2 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Vector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Vector.cs @@ -63,10 +63,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoAdd(Complex32 scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] + scalar); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) + scalar); + } } /// @@ -80,10 +80,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoAdd(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] + other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) + other.At(index)); + } } /// @@ -111,10 +111,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoSubtract(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] - other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) - other.At(index)); + } } /// @@ -128,10 +128,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoMultiply(Complex32 scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] * scalar); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) * scalar); + } } /// @@ -145,10 +145,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoDivide(Complex32 scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] / scalar); + DoMultiply(1 / scalar, result); } /// @@ -158,10 +155,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The vector to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] * other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) * other.At(index)); + } } /// @@ -171,10 +168,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The vector to store the result of the pointwise division. protected override void DoPointwiseDivide(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] / other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) / other.At(index)); + } } /// @@ -191,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return CommonParallel.Aggregate( 0, Count, - i => this[i] * other[i]); + i => At(i) * other.At(i)); } /// @@ -200,7 +197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The value of the absolute minimum element. public override Complex32 AbsoluteMinimum() { - return this[AbsoluteMinimumIndex()].Magnitude; + return At(AbsoluteMinimumIndex()).Magnitude; } /// @@ -210,10 +207,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public override int AbsoluteMinimumIndex() { var index = 0; - var min = this[index].Magnitude; + var min = At(index).Magnitude; for (var i = 1; i < Count; i++) { - var test = this[i].Magnitude; + var test = At(i).Magnitude; if (test < min) { index = i; @@ -230,7 +227,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The value of the absolute maximum element. public override Complex32 AbsoluteMaximum() { - return this[AbsoluteMaximumIndex()].Magnitude; + return At(AbsoluteMaximumIndex()).Magnitude; } /// @@ -240,10 +237,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public override int AbsoluteMaximumIndex() { var index = 0; - var max = this[index].Magnitude; + var max = At(index).Magnitude; for (var i = 1; i < Count; i++) { - var test = this[i].Magnitude; + var test = At(i).Magnitude; if (test > max) { index = i; @@ -263,7 +260,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return CommonParallel.Aggregate( 0, Count, - i => this[i]); + i => At(i)); } /// @@ -275,7 +272,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return CommonParallel.Aggregate( 0, Count, - i => this[i].Magnitude); + i => At(i).Magnitude); } /// @@ -285,7 +282,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The p value. /// /// - /// Scalar ret = (sum(abs(this[i])^p))^(1/p) + /// Scalar ret = (sum(abs(At(i))^p))^(1/p) /// public override Complex32 Norm(double p) { @@ -299,14 +296,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return CommonParallel.Select( 0, Count, - (index, localData) => Math.Max(localData, this[index].Magnitude), + (index, localData) => Math.Max(localData, At(index).Magnitude), Common.Max); } var sum = CommonParallel.Aggregate( 0, Count, - index => Math.Pow(this[index].Magnitude, p)); + index => Math.Pow(At(index).Magnitude, p)); return (float)Math.Pow(sum, 1.0 / p); } @@ -317,11 +314,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// Target vector protected override void DoConjugate(Vector target) { - CopyTo(target); - CommonParallel.For( - 0, - Count, - index => target[index] = this[index].Conjugate()); + for (var index = 0; index < Count; index++) + { + target.At(index, At(index).Conjugate()); + } } /// @@ -336,10 +332,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public override Vector Negate() { var result = CreateVector(Count); - CommonParallel.For( - 0, - Count, - index => result[index] = -this[index]); + + for (var index = 0; index < Count; index++) + { + result.At(index, -At(index)); + } return result; } @@ -410,7 +407,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var vector = CreateVector(length); for (var index = 0; index < length; index++) { - vector[index] = new Complex32(Convert.ToSingle(randomDistribution.Sample()), Convert.ToSingle(randomDistribution.Sample())); + vector.At(index, new Complex32(Convert.ToSingle(randomDistribution.Sample()), Convert.ToSingle(randomDistribution.Sample()))); } return vector; @@ -436,7 +433,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var vector = CreateVector(length); for (var index = 0; index < length; index++) { - vector[index] = new Complex32(Convert.ToSingle(randomDistribution.Sample()), Convert.ToSingle(randomDistribution.Sample())); + vector.At(index, new Complex32(Convert.ToSingle(randomDistribution.Sample()), Convert.ToSingle(randomDistribution.Sample()))); } return vector; diff --git a/src/Numerics/LinearAlgebra/Double/DenseVector.cs b/src/Numerics/LinearAlgebra/Double/DenseVector.cs index 1a39a081..d3003ff0 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseVector.cs @@ -1393,5 +1393,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double { Array.Clear(Data, 0, Data.Length); } + + /// Gets the value at the given . + /// The index of the value to get or set. + /// The value of the vector at the given . + internal protected override double At(int index) + { + return Data[index]; + } + + /// Sets the at the given . + /// The index of the value to get or set. + /// The value to set. + internal protected override void At(int index, double value) + { + Data[index] = value; + } } } diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs index 0aa76e6d..52e15e65 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs @@ -187,7 +187,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// the values from the given matrix. /// /// The matrix to copy. - public SparseMatrix(Matrix matrix) : base(matrix.RowCount, matrix.ColumnCount) + public SparseMatrix(Matrix matrix) : this(matrix.RowCount, matrix.ColumnCount) { var sparseMatrix = matrix as SparseMatrix; @@ -608,8 +608,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { lock (_lockObject) { - var index = FindItem(row, column); - return index >= 0 ? _nonZeroValues[index] : 0.0; + return GetValueAt(row, column); } } @@ -634,6 +633,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double } #region Internal methods - CRS storage implementation + + /// + /// Retrieves the requested element without range checking. + /// + /// + /// The row of the element. + /// + /// + /// The column of the element. + /// + /// + /// The requested element. + /// + private double GetValueAt(int row, int column) + { + var index = FindItem(row, column); + return index >= 0 ? _nonZeroValues[index] : 0.0; + } + /// /// Created this method because we cannot call "virtual At" in constructor of the class, but we need to do it /// diff --git a/src/Numerics/LinearAlgebra/Double/SparseVector.cs b/src/Numerics/LinearAlgebra/Double/SparseVector.cs index f80d0277..b8be1d2a 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseVector.cs @@ -38,13 +38,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// /// A vector with sparse storage. /// + /// The sparse vector is not thread safe. public class SparseVector : Vector { - /// - /// Lock object for the indexer. - /// - private readonly object _lockObject = new object(); - /// /// Gets the vector's internal data. The array containing the actual values; only the non-zero values are stored. /// @@ -130,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { for (var i = 0; i < other.Count; i++) { - this[i] = other[i]; + this[i] = other.At(i); } } else @@ -168,9 +164,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// /// The array to create this vector from. /// The vector copy the array. Any changes to the vector will NOT change the array. - public SparseVector(double[] array) : this(array.Length) + public SparseVector(IList array) : this(array.Count) { - for (var i = 0; i < array.Length; i++) + for (var i = 0; i < array.Count; i++) { this[i] = array[i]; } @@ -187,7 +183,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double var matrix = new SparseMatrix(Count, 1); for (var i = 0; i < NonZerosCount; i++) { - matrix[_nonZeroIndices[i], 0] = _nonZeroValues[i]; + matrix.At(_nonZeroIndices[i], 0, _nonZeroValues[i]); } return matrix; @@ -202,55 +198,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double var matrix = new SparseMatrix(1, Count); for (var i = 0; i < NonZerosCount; i++) { - matrix[0, _nonZeroIndices[i]] = _nonZeroValues[i]; + matrix.At(0, _nonZeroIndices[i], _nonZeroValues[i]); } return matrix; } - /// Gets or sets the value at the given . - /// The index of the value to get or set. - /// The value of the vector at the given . - /// If is negative or - /// greater than the size of the vector. - public override double this[int index] - { - get - { - // If index is out of bounds - if ((index < 0) || (index >= Count)) - { - throw new IndexOutOfRangeException(); - } - - lock (_lockObject) - { - // Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count" - var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); - if (itemIndex >= 0) - { - return _nonZeroValues[itemIndex]; - } - } - - return 0.0; - } - - set - { - // If index is out of bounds - if ((index < 0) || (index >= Count)) - { - throw new IndexOutOfRangeException(); - } - - lock (_lockObject) - { - SetValue(index, value); - } - } - } - /// /// Creates a matrix with the given dimensions using the same storage type /// as this vector. @@ -324,10 +277,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double var otherVector = target as SparseVector; if (otherVector == null) { - CommonParallel.For( - 0, - Count, - index => target[index] = this[index]); + target.Clear(); + + for (var index = 0; index < NonZerosCount; index++) + { + target.At(_nonZeroIndices[index], _nonZeroValues[index]); + } } else { @@ -343,145 +298,39 @@ namespace MathNet.Numerics.LinearAlgebra.Double #region Operators and supplementary functions - /// - /// Adds a scalar to each element of the vector. - /// - /// The scalar to add. - /// A copy of the vector with the scalar added. - public override Vector Add(double scalar) - { - if (scalar == 0.0) - { - return Clone(); - } - - var copy = (SparseVector)Clone(); - for (var i = 0; i < Count; i++) - { - copy[i] += scalar; - } - - return copy; - } - /// /// Adds a scalar to each element of the vector and stores the result in the result vector. /// - /// The scalar to add. - /// The vector to store the result of the addition. - /// If the result vector is . - /// If this vector and are not the same size. - public override void Add(double scalar, Vector result) - { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - base.Add(scalar, result); - } - - /// - /// Adds another vector to this vector. - /// - /// The vector to add to this one. - /// A new vector containing the sum of both vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector Add(Vector other) - { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - var sparseVector = other as SparseVector; - - if (sparseVector == null) - { - return base.Add(other); - } - - var copy = (SparseVector)Clone(); - copy.AddScaledSparseVector(1.0, sparseVector); - return copy; - } - - /// - /// Adds the scaled sparse vector. - /// - /// The alpha. - /// The other. - private void AddScaledSparseVector(double alpha, SparseVector other) + /// + /// The scalar to add. + /// + /// + /// The vector to store the result of the addition. + /// + protected override void DoAdd(double scalar, Vector result) { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) + if (scalar == 0.0) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } + if (!ReferenceEquals(this, result)) + { + CopyTo(result); + } - if (alpha == 0.0) - { return; } - // I don't use ILinearAlgebraProvider because we will get no benefit due to "lock" in this[index] - // Possible fucniton in ILinearAlgebraProvider may be AddSparseVectorToScaledSparseVector(T[] y, int[] yIndices, T alpha, T[] x, int[] xIndices); - // But it require to develop value setting algorithm and due to "lock" it will be even more greedy then implemented below - if (ReferenceEquals(this, other)) + if (ReferenceEquals(this, result)) { - // Adding the same instance of sparse vector. That means if we modify "this" then "other" will be modified too. - // To avoid such problem lets change values in internal storage of "this" - if (alpha == 1.0) - { - for (var i = 0; i < NonZerosCount; i++) - { - _nonZeroValues[i] += _nonZeroValues[i]; - } - } - else if (alpha == -1.0) - { - Clear(); // Vector is subtracted from itself - return; - } - else - { - for (var i = 0; i < NonZerosCount; i++) - { - _nonZeroValues[i] += alpha * _nonZeroValues[i]; - } - } + CommonParallel.For( + 0, + NonZerosCount, + index => _nonZeroValues[index] += scalar); } else { - // "this" and "other" are different objects, so by modifying "this" the "other" object will not be changed - if (alpha == 1.0) + for (var index = 0; index < Count; index++) { - for (var i = 0; i < other.NonZerosCount; i++) - { - this[other._nonZeroIndices[i]] += other._nonZeroValues[i]; - } - } - else - { - for (var i = 0; i < other.NonZerosCount; i++) - { - this[other._nonZeroIndices[i]] += alpha * other._nonZeroValues[i]; - } + result.At(index, At(index) + scalar); } } } @@ -489,53 +338,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// /// Adds another vector to this vector and stores the result into the result vector. /// - /// The vector to add to this one. - /// The vector to store the result of the addition. - /// If the other vector is . - /// If the result vector is . - /// If this vector and are not the same size. - /// If this vector and are not the same size. - public override void Add(Vector other, Vector result) + /// + /// The vector to add to this one. + /// + /// + /// The vector to store the result of the addition. + /// + protected override void DoAdd(Vector other, Vector result) { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - if (Count != result.Count) + if (ReferenceEquals(this, result)) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) - { - var tmp = Add(other); - tmp.CopyTo(result); + CommonParallel.For( + 0, + NonZerosCount, + index => _nonZeroValues[index] += _nonZeroValues[index]); } else { - var sparse = result as SparseVector; - if (sparse == null) - { - base.Add(other, result); - } - else + for (var index = 0; index < Count; index++) { - var sparseother = other as SparseVector; - if (sparseother == null) - { - base.Add(other, result); - } - else - { - CopyTo(result); - sparse.AddScaledSparseVector(1.0, sparseother); - } + result.At(index, At(index) + other.At(index)); } } } @@ -585,131 +407,40 @@ namespace MathNet.Numerics.LinearAlgebra.Double return leftSide.Add(rightSide); } - /// - /// Subtracts a scalar from each element of the vector. - /// - /// The scalar to subtract. - /// A new vector containing the subtraction of this vector and the scalar. - public override Vector Subtract(double scalar) - { - if (scalar == 0.0) - { - return Clone(); - } - - var copy = (SparseVector)Clone(); - for (var i = 0; i < Count; i++) - { - copy[i] -= scalar; - } - - return copy; - } - /// /// Subtracts a scalar from each element of the vector and stores the result in the result vector. /// - /// The scalar to subtract. - /// The vector to store the result of the subtraction. - /// If the result vector is . - /// If this vector and are not the same size. - public override void Subtract(double scalar, Vector result) - { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - base.Subtract(scalar, result); - } - - /// - /// Subtracts another vector from this vector. - /// - /// The vector to subtract from this one. - /// A new vector containing the subtraction of the the two vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector Subtract(Vector other) + /// + /// The scalar to subtract. + /// + /// + /// The vector to store the result of the subtraction. + /// + protected override void DoSubtract(double scalar, Vector result) { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - var sparseVector = other as SparseVector; - - if (sparseVector == null) - { - return base.Subtract(other); - } - - var copy = (SparseVector)Clone(); - copy.AddScaledSparseVector(-1.0, sparseVector); - return copy; + DoAdd(-scalar, result); } /// /// Subtracts another vector to this vector and stores the result into the result vector. /// - /// The vector to subtract from this one. - /// The vector to store the result of the subtraction. - /// If the other vector is . - /// If the result vector is . - /// If this vector and are not the same size. - /// If this vector and are not the same size. - public override void Subtract(Vector other, Vector result) + /// + /// The vector to subtract from this one. + /// + /// + /// The vector to store the result of the subtraction. + /// + protected override void DoSubtract(Vector other, Vector result) { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - if (Count != result.Count) + if (ReferenceEquals(this, other)) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); + result.Clear(); + return; } - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) + for (var index = 0; index < Count; index++) { - var tmp = Subtract(other); - tmp.CopyTo(result); - } - else - { - var sparse = result as SparseVector; - if (sparse == null) - { - base.Subtract(other, result); - } - else - { - var sparseother = other as SparseVector; - if (sparseother == null) - { - base.Subtract(other, result); - } - else - { - CopyTo(result); - sparse.AddScaledSparseVector(-1.0, sparseother); - } - } + result.At(index, At(index) - other.At(index)); } } @@ -777,6 +508,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double 0, NonZerosCount, index => result._nonZeroValues[index] = -_nonZeroValues[index]); + Buffer.BlockCopy(_nonZeroIndices, 0, result._nonZeroIndices, 0, NonZerosCount * Constants.SizeOfInt); } @@ -806,30 +538,81 @@ namespace MathNet.Numerics.LinearAlgebra.Double } /// - /// Computes the dot product between this vector and another vector. + /// Multiplies a scalar to each element of the vector and stores the result in the result vector. /// - /// The other vector to add. - /// The result of the addition. - /// If is not of the same size. - /// If is . - public override double DotProduct(Vector other) + /// + /// The scalar to multiply. + /// + /// + /// The vector to store the result of the multiplication. + /// + protected override void DoMultiply(double scalar, Vector result) { - if (other == null) + if (scalar == 1.0) { - throw new ArgumentNullException("other"); + if (!ReferenceEquals(this, result)) + { + CopyTo(result); + } + + return; } - if (Count != other.Count) + if (scalar == 0) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + result.Clear(); + return; } - double result = 0; + var sparseResult = result as SparseVector; + if (sparseResult == null) + { + result.Clear(); + for (var index = 0; index < NonZerosCount; index++) + { + result.At(_nonZeroIndices[index], scalar * _nonZeroValues[index]); + } + } + else + { + if (!ReferenceEquals(this, result)) + { + sparseResult.NonZerosCount = NonZerosCount; + sparseResult._nonZeroIndices = new int[NonZerosCount]; + Buffer.BlockCopy(_nonZeroIndices, 0, sparseResult._nonZeroIndices, 0, _nonZeroIndices.Length * Constants.SizeOfInt); + sparseResult._nonZeroValues = new double[_nonZeroValues.Length]; + } - // base implementation iterates though all elements, but we need only take non-zeros - for (var i = 0; i < NonZerosCount; i++) + Control.LinearAlgebraProvider.ScaleArray(scalar, _nonZeroValues, sparseResult._nonZeroValues); + } + } + + /// + /// Computes the dot product between this vector and another vector. + /// + /// + /// The other vector to add. + /// + /// s + /// The result of the addition. + /// + protected override double DoDotProduct(Vector other) + { + var result = 0.0; + + if (ReferenceEquals(this, other)) + { + for (var i = 0; i < NonZerosCount; i++) + { + result += _nonZeroValues[i] * _nonZeroValues[i]; + } + } + else { - result += _nonZeroValues[i] * other[_nonZeroIndices[i]]; + for (var i = 0; i < NonZerosCount; i++) + { + result += _nonZeroValues[i] * other.At(_nonZeroIndices[i]); + } } return result; @@ -972,7 +755,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double var result = new SparseVector(length); for (var i = index; i < index + length; i++) { - result[i - index] = this[i]; + result.At(i - index, At(i)); } return result; @@ -998,7 +781,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var i = 0; i < values.Length; i++) { - this[i] = values[i]; + At(i, values[i]); } } @@ -1083,35 +866,51 @@ namespace MathNet.Numerics.LinearAlgebra.Double } /// - /// Pointwise multiplies this vector with another vector. + /// Pointwise multiplies this vector with another vector and stores the result into the result vector. /// /// The vector to pointwise multiply with this one. - /// A new vector which is the pointwise multiplication of the two vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector PointwiseMultiply(Vector other) + /// The vector to store the result of the pointwise multiplication. + protected override void DoPointwiseMultiply(Vector other, Vector result) { - if (other == null) + if (ReferenceEquals(this, other)) { - throw new ArgumentNullException("other"); + for (var i = 0; i < NonZerosCount; i++) + { + _nonZeroValues[i] *= _nonZeroValues[i]; + } } - - if (Count != other.Count) + else { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + for (var i = 0; i < NonZerosCount; i++) + { + var index = _nonZeroIndices[i]; + result.At(index, other.At(index) * _nonZeroValues[i]); + } } + } - var copy = new SparseVector(Count); - for (var i = 0; i < _nonZeroIndices.Length; i++) + /// + /// Pointwise multiplies this vector with another vector and stores the result into the result vector. + /// + /// The vector to pointwise multiply with this one. + /// The vector to store the result of the pointwise multiplication. + protected override void DoPointwiseDivide(Vector other, Vector result) + { + if (ReferenceEquals(this, other)) { - var d = _nonZeroValues[i] * other[_nonZeroIndices[i]]; - if (d != 0.0) + for (var i = 0; i < NonZerosCount; i++) { - copy[_nonZeroIndices[i]] = d; + _nonZeroValues[i] /= _nonZeroValues[i]; + } + } + else + { + for (var i = 0; i < NonZerosCount; i++) + { + var index = _nonZeroIndices[i]; + result.At(index, _nonZeroValues[i] / other.At(index)); } } - - return copy; } /// @@ -1122,7 +921,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Matrix M[i,j] = u[i]*v[j] /// If the u vector is . /// If the v vector is . - public static Matrix /*SparseMatrix*/ OuterProduct(SparseVector u, SparseVector v) + public static Matrix OuterProduct(SparseVector u, SparseVector v) { if (u == null) { @@ -1361,13 +1160,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double #endregion + /// + /// Gets the value at the given index. + /// + /// Value real index in array + /// The value at the given index. + internal protected override double At(int index) + { + // Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count" + var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); + return itemIndex >= 0 ? _nonZeroValues[itemIndex] : 0.0; + } + /// /// Delete, Add or Update the value in NonZeroValues and NonZeroIndices /// /// Value real index in array /// The value to set. /// This method assume that index is between 0 and Array Size - private void SetValue(int index, double value) + internal protected override void At(int index, double value) { // Search if "index" already exists in range "0 - real nonzero values count" var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); @@ -1562,5 +1373,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double yield return new KeyValuePair(_nonZeroIndices[i], _nonZeroValues[i]); } } + + /// + /// Returns the data contained in the vector as an array. + /// + /// + /// The vector's data as an array. + /// + public override double[] ToArray() + { + var ret = new double[Count]; + for (var i = 0; i < NonZerosCount; i++) + { + ret[_nonZeroIndices[i]] = _nonZeroValues[i]; + } + + return ret; + } } } diff --git a/src/Numerics/LinearAlgebra/Double/Vector.cs b/src/Numerics/LinearAlgebra/Double/Vector.cs index 7baefb6a..dff69988 100644 --- a/src/Numerics/LinearAlgebra/Double/Vector.cs +++ b/src/Numerics/LinearAlgebra/Double/Vector.cs @@ -62,10 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoAdd(double scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] + scalar); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) + scalar); + } } /// @@ -79,10 +79,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoAdd(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] + other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) + other.At(index)); + } } /// @@ -110,10 +110,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoSubtract(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] - other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) - other.At(index)); + } } /// @@ -127,10 +127,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoMultiply(double scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] * scalar); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) * scalar); + } } /// @@ -144,10 +144,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoDivide(double scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] / scalar); + DoMultiply(1 / scalar, result); } /// @@ -157,10 +154,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The vector to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] * other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) * other.At(index)); + } } /// @@ -170,10 +167,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The vector to store the result of the pointwise division. protected override void DoPointwiseDivide(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] / other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) / other.At(index)); + } } /// @@ -190,7 +187,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double return CommonParallel.Aggregate( 0, Count, - i => this[i] * other[i]); + i => At(i) * other.At(i)); } /// @@ -199,7 +196,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The value of the absolute minimum element. public override double AbsoluteMinimum() { - return Math.Abs(this[AbsoluteMinimumIndex()]); + return Math.Abs(At(AbsoluteMinimumIndex())); } /// @@ -209,10 +206,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override int AbsoluteMinimumIndex() { var index = 0; - var min = Math.Abs(this[index]); + var min = Math.Abs(At(index)); for (var i = 1; i < Count; i++) { - var test = Math.Abs(this[i]); + var test = Math.Abs(At(i)); if (test < min) { index = i; @@ -229,7 +226,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The value of the absolute maximum element. public override double AbsoluteMaximum() { - return Math.Abs(this[AbsoluteMaximumIndex()]); + return Math.Abs(At(AbsoluteMaximumIndex())); } /// @@ -239,10 +236,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override int AbsoluteMaximumIndex() { var index = 0; - var max = Math.Abs(this[index]); + var max = Math.Abs(At(index)); for (var i = 1; i < Count; i++) { - var test = Math.Abs(this[i]); + var test = Math.Abs(At(i)); if (test > max) { index = i; @@ -262,7 +259,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double return CommonParallel.Aggregate( 0, Count, - i => this[i]); + i => At(i)); } /// @@ -274,7 +271,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double return CommonParallel.Aggregate( 0, Count, - i => Math.Abs(this[i])); + i => Math.Abs(At(i))); } /// @@ -298,14 +295,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double return CommonParallel.Select( 0, Count, - (index, localData) => Math.Max(localData, Math.Abs(this[index])), + (index, localData) => Math.Max(localData, Math.Abs(At(index))), Math.Max); } var sum = CommonParallel.Aggregate( 0, Count, - index => Math.Pow(Math.Abs(this[index]), p)); + index => Math.Pow(Math.Abs(At(index)), p)); return Math.Pow(sum, 1.0 / p); } @@ -336,10 +333,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override Vector Negate() { var result = CreateVector(Count); - CommonParallel.For( - 0, - Count, - index => result[index] = -this[index]); + + for (var index = 0; index < Count; index++) + { + result.At(index, -At(index)); + } return result; } @@ -351,10 +349,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override int MaximumIndex() { var index = 0; - var max = this[index]; + var max = At(index); for (var i = 1; i < Count; i++) { - var test = this[i]; + var test = At(i); if (test > max) { index = i; @@ -372,10 +370,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override int MinimumIndex() { var index = 0; - var min = this[index]; + var min = At(index); for (var i = 1; i < Count; i++) { - var test = this[i]; + var test = At(i); if (test < min) { index = i; @@ -434,7 +432,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double var v = CreateVector(length); for (var index = 0; index < Count; index++) { - v[index] = randomDistribution.Sample(); + v.At(index, randomDistribution.Sample()); } return v; @@ -460,7 +458,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double var v = CreateVector(length); for (var index = 0; index < Count; index++) { - this[index] = randomDistribution.Sample(); + At(index, randomDistribution.Sample()); } return v; diff --git a/src/Numerics/LinearAlgebra/Generic/Vector.cs b/src/Numerics/LinearAlgebra/Generic/Vector.cs index a1c72a5b..7d9fcd76 100644 --- a/src/Numerics/LinearAlgebra/Generic/Vector.cs +++ b/src/Numerics/LinearAlgebra/Generic/Vector.cs @@ -93,10 +93,27 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// The value of the vector at the given . /// If is negative or /// greater than the size of the vector. - public abstract T this[int index] + public virtual T this[int index] { - get; - set; + get + { + if (index < 0 || index >= Count) + { + throw new ArgumentOutOfRangeException("index"); + } + + return At(index); + } + + set + { + if (index < 0 || index >= Count) + { + throw new ArgumentOutOfRangeException("index"); + } + + At(index, value); + } } /// @@ -715,7 +732,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// Pointwise divide this vector with another vector and stores the result into the result vector. /// /// The vector to pointwise divide this one by. - /// The vector to store the result of the pointwise division. protected abstract void DoPointwiseDivide(Vector other, Vector result); /// @@ -1207,7 +1223,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic var ret = new T[Count]; for (var i = 0; i < ret.Length; i++) { - ret[i] = this[i]; + ret[i] = At(i); } return ret; @@ -1534,5 +1550,15 @@ namespace MathNet.Numerics.LinearAlgebra.Generic { CommonParallel.For(0, Count, index => this[index] = default(T)); } + + /// Gets the value at the given . + /// The index of the value to get or set. + /// The value of the vector at the given . + internal protected abstract T At(int index); + + /// Sets the at the given . + /// The index of the value to get or set. + /// The value to set. + internal protected abstract void At(int index, T value); } } diff --git a/src/Numerics/LinearAlgebra/Single/DenseVector.cs b/src/Numerics/LinearAlgebra/Single/DenseVector.cs index 3f068bd5..b3a8f7fc 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseVector.cs @@ -1393,5 +1393,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single { Array.Clear(Data, 0, Data.Length); } + + /// Gets the value at the given . + /// The index of the value to get or set. + /// The value of the vector at the given . + internal protected override float At(int index) + { + return Data[index]; + } + + /// Sets the at the given . + /// The index of the value to get or set. + /// The value to set. + internal protected override void At(int index, float value) + { + Data[index] = value; + } } } diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs index 51690d84..037dd6b8 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs @@ -188,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// /// The matrix to copy. public SparseMatrix(Matrix matrix) - : base(matrix.RowCount, matrix.ColumnCount) + : this(matrix.RowCount, matrix.ColumnCount) { var sparseMatrix = matrix as SparseMatrix; diff --git a/src/Numerics/LinearAlgebra/Single/SparseVector.cs b/src/Numerics/LinearAlgebra/Single/SparseVector.cs index 4dbfd83d..142971aa 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseVector.cs @@ -38,13 +38,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// /// A vector with sparse storage. /// + /// The sparse vector is not thread safe. public class SparseVector : Vector { - /// - /// Lock object for the indexer. - /// - private readonly object _lockObject = new object(); - /// /// Gets the vector's internal data. The array containing the actual values; only the non-zero values are stored. /// @@ -130,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single { for (var i = 0; i < other.Count; i++) { - this[i] = other[i]; + this[i] = other.At(i); } } else @@ -168,9 +164,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// /// The array to create this vector from. /// The vector copy the array. Any changes to the vector will NOT change the array. - public SparseVector(float[] array) : this(array.Length) + public SparseVector(IList array) : this(array.Count) { - for (var i = 0; i < array.Length; i++) + for (var i = 0; i < array.Count; i++) { this[i] = array[i]; } @@ -187,7 +183,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single var matrix = new SparseMatrix(Count, 1); for (var i = 0; i < NonZerosCount; i++) { - matrix[_nonZeroIndices[i], 0] = _nonZeroValues[i]; + matrix.At(_nonZeroIndices[i], 0, _nonZeroValues[i]); } return matrix; @@ -202,7 +198,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single var matrix = new SparseMatrix(1, Count); for (var i = 0; i < NonZerosCount; i++) { - matrix[0, _nonZeroIndices[i]] = _nonZeroValues[i]; + matrix.At(0, _nonZeroIndices[i], _nonZeroValues[i]); } return matrix; @@ -223,17 +219,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single throw new IndexOutOfRangeException(); } - lock (_lockObject) - { - // Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count" - var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); - if (itemIndex >= 0) - { - return _nonZeroValues[itemIndex]; - } - } - - return 0.0f; + return At(index); } set @@ -244,10 +230,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single throw new IndexOutOfRangeException(); } - lock (_lockObject) - { - SetValue(index, value); - } + At(index, value); } } @@ -324,10 +307,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single var otherVector = target as SparseVector; if (otherVector == null) { - CommonParallel.For( - 0, - Count, - index => target[index] = this[index]); + target.Clear(); + + for (var index = 0; index < NonZerosCount; index++) + { + target.At(_nonZeroIndices[index], _nonZeroValues[index]); + } } else { @@ -343,145 +328,39 @@ namespace MathNet.Numerics.LinearAlgebra.Single #region Operators and supplementary functions - /// - /// Adds a scalar to each element of the vector. - /// - /// The scalar to add. - /// A copy of the vector with the scalar added. - public override Vector Add(float scalar) - { - if (scalar == 0.0) - { - return Clone(); - } - - var copy = (SparseVector)Clone(); - for (var i = 0; i < Count; i++) - { - copy[i] += scalar; - } - - return copy; - } - /// /// Adds a scalar to each element of the vector and stores the result in the result vector. /// - /// The scalar to add. - /// The vector to store the result of the addition. - /// If the result vector is . - /// If this vector and are not the same size. - public override void Add(float scalar, Vector result) - { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - base.Add(scalar, result); - } - - /// - /// Adds another vector to this vector. - /// - /// The vector to add to this one. - /// A new vector containing the sum of both vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector Add(Vector other) - { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - var sparseVector = other as SparseVector; - - if (sparseVector == null) - { - return base.Add(other); - } - - var copy = (SparseVector)Clone(); - copy.AddScaledSparseVector(1.0f, sparseVector); - return copy; - } - - /// - /// Adds the scaled sparse vector. - /// - /// The alpha. - /// The other. - private void AddScaledSparseVector(float alpha, SparseVector other) + /// + /// The scalar to add. + /// + /// + /// The vector to store the result of the addition. + /// + protected override void DoAdd(float scalar, Vector result) { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) + if (scalar == 0.0f) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } + if (!ReferenceEquals(this, result)) + { + CopyTo(result); + } - if (alpha == 0.0) - { return; } - // I don't use ILinearAlgebraProvider because we will get no benefit due to "lock" in this[index] - // Possible fucniton in ILinearAlgebraProvider may be AddSparseVectorToScaledSparseVector(T[] y, int[] yIndices, T alpha, T[] x, int[] xIndices); - // But it require to develop value setting algorithm and due to "lock" it will be even more greedy then implemented below - if (ReferenceEquals(this, other)) + if (ReferenceEquals(this, result)) { - // Adding the same instance of sparse vector. That means if we modify "this" then "other" will be modified too. - // To avoid such problem lets change values in internal storage of "this" - if (alpha == 1.0) - { - for (var i = 0; i < NonZerosCount; i++) - { - _nonZeroValues[i] += _nonZeroValues[i]; - } - } - else if (alpha == -1.0) - { - Clear(); // Vector is subtracted from itself - return; - } - else - { - for (var i = 0; i < NonZerosCount; i++) - { - _nonZeroValues[i] += alpha * _nonZeroValues[i]; - } - } + CommonParallel.For( + 0, + NonZerosCount, + index => _nonZeroValues[index] += scalar); } else { - // "this" and "other" are different objects, so by modifying "this" the "other" object will not be changed - if (alpha == 1.0) - { - for (var i = 0; i < other.NonZerosCount; i++) - { - this[other._nonZeroIndices[i]] += other._nonZeroValues[i]; - } - } - else + for (var index = 0; index < Count; index++) { - for (var i = 0; i < other.NonZerosCount; i++) - { - this[other._nonZeroIndices[i]] += alpha * other._nonZeroValues[i]; - } + result.At(index, At(index) + scalar); } } } @@ -489,53 +368,26 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// /// Adds another vector to this vector and stores the result into the result vector. /// - /// The vector to add to this one. - /// The vector to store the result of the addition. - /// If the other vector is . - /// If the result vector is . - /// If this vector and are not the same size. - /// If this vector and are not the same size. - public override void Add(Vector other, Vector result) + /// + /// The vector to add to this one. + /// + /// + /// The vector to store the result of the addition. + /// + protected override void DoAdd(Vector other, Vector result) { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) + if (ReferenceEquals(this, result)) { - var tmp = Add(other); - tmp.CopyTo(result); + CommonParallel.For( + 0, + NonZerosCount, + index => _nonZeroValues[index] += _nonZeroValues[index]); } else { - var sparse = result as SparseVector; - if (sparse == null) + for (var index = 0; index < Count; index++) { - base.Add(other, result); - } - else - { - var sparseother = other as SparseVector; - if (sparseother == null) - { - base.Add(other, result); - } - else - { - CopyTo(result); - sparse.AddScaledSparseVector(1.0f, sparseother); - } + result.At(index, At(index) + other.At(index)); } } } @@ -585,131 +437,40 @@ namespace MathNet.Numerics.LinearAlgebra.Single return leftSide.Add(rightSide); } - /// - /// Subtracts a scalar from each element of the vector. - /// - /// The scalar to subtract. - /// A new vector containing the subtraction of this vector and the scalar. - public override Vector Subtract(float scalar) - { - if (scalar == 0.0) - { - return Clone(); - } - - var copy = (SparseVector)Clone(); - for (var i = 0; i < Count; i++) - { - copy[i] -= scalar; - } - - return copy; - } - /// /// Subtracts a scalar from each element of the vector and stores the result in the result vector. /// - /// The scalar to subtract. - /// The vector to store the result of the subtraction. - /// If the result vector is . - /// If this vector and are not the same size. - public override void Subtract(float scalar, Vector result) - { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != result.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); - } - - base.Subtract(scalar, result); - } - - /// - /// Subtracts another vector from this vector. - /// - /// The vector to subtract from this one. - /// A new vector containing the subtraction of the the two vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector Subtract(Vector other) + /// + /// The scalar to subtract. + /// + /// + /// The vector to store the result of the subtraction. + /// + protected override void DoSubtract(float scalar, Vector result) { - if (other == null) - { - throw new ArgumentNullException("other"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - var sparseVector = other as SparseVector; - - if (sparseVector == null) - { - return base.Subtract(other); - } - - var copy = (SparseVector)Clone(); - copy.AddScaledSparseVector(-1.0f, sparseVector); - return copy; + DoAdd(-scalar, result); } /// /// Subtracts another vector to this vector and stores the result into the result vector. /// - /// The vector to subtract from this one. - /// The vector to store the result of the subtraction. - /// If the other vector is . - /// If the result vector is . - /// If this vector and are not the same size. - /// If this vector and are not the same size. - public override void Subtract(Vector other, Vector result) + /// + /// The vector to subtract from this one. + /// + /// + /// The vector to store the result of the subtraction. + /// + protected override void DoSubtract(Vector other, Vector result) { - if (result == null) - { - throw new ArgumentNullException("result"); - } - - if (Count != other.Count) - { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); - } - - if (Count != result.Count) + if (ReferenceEquals(this, other)) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); + result.Clear(); + return; } - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) - { - var tmp = Subtract(other); - tmp.CopyTo(result); - } - else + for (var index = 0; index < Count; index++) { - var sparse = result as SparseVector; - if (sparse == null) - { - base.Subtract(other, result); - } - else - { - var sparseother = other as SparseVector; - if (sparseother == null) - { - base.Subtract(other, result); - } - else - { - CopyTo(result); - sparse.AddScaledSparseVector(-1.0f, sparseother); - } - } + result.At(index, At(index) - other.At(index)); } } @@ -784,52 +545,81 @@ namespace MathNet.Numerics.LinearAlgebra.Single } /// - /// Multiplies a scalar to each element of the vector. + /// Multiplies a scalar to each element of the vector and stores the result in the result vector. /// - /// The scalar to multiply. - /// A new vector that is the multiplication of the vector and the scalar. - public override Vector Multiply(float scalar) + /// + /// The scalar to multiply. + /// + /// + /// The vector to store the result of the multiplication. + /// + protected override void DoMultiply(float scalar, Vector result) { if (scalar == 1.0) { - return Clone(); + if (!ReferenceEquals(this, result)) + { + CopyTo(result); + } + + return; } if (scalar == 0) { - return new SparseVector(Count); + result.Clear(); + return; } - var copy = new SparseVector(this); - Control.LinearAlgebraProvider.ScaleArray(scalar, copy._nonZeroValues, copy._nonZeroValues); - return copy; + var sparseResult = result as SparseVector; + if (sparseResult == null) + { + result.Clear(); + for (var index = 0; index < NonZerosCount; index++) + { + result.At(_nonZeroIndices[index], scalar * _nonZeroValues[index]); + } + } + else + { + if (!ReferenceEquals(this, result)) + { + sparseResult.NonZerosCount = NonZerosCount; + sparseResult._nonZeroIndices = new int[NonZerosCount]; + Buffer.BlockCopy(_nonZeroIndices, 0, sparseResult._nonZeroIndices, 0, _nonZeroIndices.Length * Constants.SizeOfInt); + sparseResult._nonZeroValues = new float[_nonZeroValues.Length]; + } + + Control.LinearAlgebraProvider.ScaleArray(scalar, _nonZeroValues, sparseResult._nonZeroValues); + } } /// /// Computes the dot product between this vector and another vector. /// - /// The other vector to add. - /// The result of the addition. - /// If is not of the same size. - /// If is . - public override float DotProduct(Vector other) + /// + /// The other vector to add. + /// + /// s + /// The result of the addition. + /// + protected override float DoDotProduct(Vector other) { - if (other == null) - { - throw new ArgumentNullException("other"); - } + var result = 0.0f; - if (Count != other.Count) + if (ReferenceEquals(this, other)) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + for (var i = 0; i < NonZerosCount; i++) + { + result += _nonZeroValues[i] * _nonZeroValues[i]; + } } - - float result = 0; - - // base implementation iterates though all elements, but we need only take non-zeros - for (var i = 0; i < NonZerosCount; i++) + else { - result += _nonZeroValues[i] * other[_nonZeroIndices[i]]; + for (var i = 0; i < NonZerosCount; i++) + { + result += _nonZeroValues[i] * other.At(_nonZeroIndices[i]); + } } return result; @@ -998,7 +788,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single for (var i = 0; i < values.Length; i++) { - this[i] = values[i]; + At(i, values[i]); } } @@ -1083,35 +873,51 @@ namespace MathNet.Numerics.LinearAlgebra.Single } /// - /// Pointwise multiplies this vector with another vector. + /// Pointwise multiplies this vector with another vector and stores the result into the result vector. /// /// The vector to pointwise multiply with this one. - /// A new vector which is the pointwise multiplication of the two vectors. - /// If the other vector is . - /// If this vector and are not the same size. - public override Vector PointwiseMultiply(Vector other) + /// The vector to store the result of the pointwise multiplication. + protected override void DoPointwiseMultiply(Vector other, Vector result) { - if (other == null) + if (ReferenceEquals(this, other)) { - throw new ArgumentNullException("other"); + for (var i = 0; i < NonZerosCount; i++) + { + _nonZeroValues[i] *= _nonZeroValues[i]; + } } - - if (Count != other.Count) + else { - throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + for (var i = 0; i < NonZerosCount; i++) + { + var index = _nonZeroIndices[i]; + result.At(index, other.At(index) * _nonZeroValues[i]); + } } + } - var copy = new SparseVector(Count); - for (var i = 0; i < _nonZeroIndices.Length; i++) + /// + /// Pointwise multiplies this vector with another vector and stores the result into the result vector. + /// + /// The vector to pointwise multiply with this one. + /// The vector to store the result of the pointwise multiplication. + protected override void DoPointwiseDivide(Vector other, Vector result) + { + if (ReferenceEquals(this, other)) { - var d = _nonZeroValues[i] * other[_nonZeroIndices[i]]; - if (d != 0.0) + for (var i = 0; i < NonZerosCount; i++) { - copy[_nonZeroIndices[i]] = d; + _nonZeroValues[i] /= _nonZeroValues[i]; + } + } + else + { + for (var i = 0; i < NonZerosCount; i++) + { + var index = _nonZeroIndices[i]; + result.At(index, _nonZeroValues[i] / other.At(index)); } } - - return copy; } /// @@ -1122,7 +928,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Matrix M[i,j] = u[i]*v[j] /// If the u vector is . /// If the v vector is . - public static Matrix /*SparseMatrix*/ OuterProduct(SparseVector u, SparseVector v) + public static Matrix OuterProduct(SparseVector u, SparseVector v) { if (u == null) { @@ -1361,13 +1167,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single #endregion + /// + /// Gets the value at the given index. + /// + /// Value real index in array + /// The value at the given index. + internal protected override float At(int index) + { + // Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count" + var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); + return itemIndex >= 0 ? _nonZeroValues[itemIndex] : 0.0f; + } + /// /// Delete, Add or Update the value in NonZeroValues and NonZeroIndices /// /// Value real index in array /// The value to set. /// This method assume that index is between 0 and Array Size - private void SetValue(int index, float value) + internal protected override void At(int index, float value) { // Search if "index" already exists in range "0 - real nonzero values count" var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index); @@ -1562,5 +1380,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single yield return new KeyValuePair(_nonZeroIndices[i], _nonZeroValues[i]); } } + + /// + /// Returns the data contained in the vector as an array. + /// + /// + /// The vector's data as an array. + /// + public override float[] ToArray() + { + var ret = new float[Count]; + for (var i = 0; i < NonZerosCount; i++) + { + ret[_nonZeroIndices[i]] = _nonZeroValues[i]; + } + + return ret; + } } } diff --git a/src/Numerics/LinearAlgebra/Single/Vector.cs b/src/Numerics/LinearAlgebra/Single/Vector.cs index 0248c1ac..e3ab2be0 100644 --- a/src/Numerics/LinearAlgebra/Single/Vector.cs +++ b/src/Numerics/LinearAlgebra/Single/Vector.cs @@ -63,10 +63,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoAdd(float scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] + scalar); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) + scalar); + } } /// @@ -80,10 +80,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoAdd(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] + other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) + other.At(index)); + } } /// @@ -111,10 +111,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoSubtract(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] - other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) - other.At(index)); + } } /// @@ -128,10 +128,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoMultiply(float scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] * scalar); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) * scalar); + } } /// @@ -145,10 +145,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoDivide(float scalar, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] / scalar); + DoMultiply(1 / scalar, result); } /// @@ -158,10 +155,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The vector to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] * other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) * other.At(index)); + } } /// @@ -171,10 +168,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The vector to store the result of the pointwise division. protected override void DoPointwiseDivide(Vector other, Vector result) { - CommonParallel.For( - 0, - Count, - index => result[index] = this[index] / other[index]); + for (var index = 0; index < Count; index++) + { + result.At(index, At(index) / other.At(index)); + } } /// @@ -191,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single return CommonParallel.Aggregate( 0, Count, - i => this[i] * other[i]); + i => At(i) * other.At(i)); } /// @@ -200,7 +197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The value of the absolute minimum element. public override float AbsoluteMinimum() { - return Math.Abs(this[AbsoluteMinimumIndex()]); + return Math.Abs(At(AbsoluteMinimumIndex())); } /// @@ -210,10 +207,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override int AbsoluteMinimumIndex() { var index = 0; - var min = Math.Abs(this[index]); + var min = Math.Abs(At(index)); for (var i = 1; i < Count; i++) { - var test = Math.Abs(this[i]); + var test = Math.Abs(At(i)); if (test < min) { index = i; @@ -230,7 +227,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The value of the absolute maximum element. public override float AbsoluteMaximum() { - return Math.Abs(this[AbsoluteMaximumIndex()]); + return Math.Abs(At(AbsoluteMaximumIndex())); } /// @@ -240,10 +237,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override int AbsoluteMaximumIndex() { var index = 0; - var max = Math.Abs(this[index]); + var max = Math.Abs(At(index)); for (var i = 1; i < Count; i++) { - var test = Math.Abs(this[i]); + var test = Math.Abs(At(i)); if (test > max) { index = i; @@ -263,7 +260,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single return CommonParallel.Aggregate( 0, Count, - i => this[i]); + i => At(i)); } /// @@ -275,7 +272,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single return CommonParallel.Aggregate( 0, Count, - i => Math.Abs(this[i])); + i => Math.Abs(At(i))); } /// @@ -285,7 +282,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The p value. /// /// - /// Scalar ret = (sum(abs(this[i])^p))^(1/p) + /// Scalar ret = (sum(abs(At(i))^p))^(1/p) /// public override float Norm(double p) { @@ -299,14 +296,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single return CommonParallel.Select( 0, Count, - (index, localData) => Math.Max(localData, Math.Abs(this[index])), + (index, localData) => Math.Max(localData, Math.Abs(At(index))), Common.Max); } var sum = CommonParallel.Aggregate( 0, Count, - index => Math.Pow(Math.Abs(this[index]), p)); + index => Math.Pow(Math.Abs(At(index)), p)); return (float)Math.Pow(sum, 1.0 / p); } @@ -337,10 +334,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override Vector Negate() { var result = CreateVector(Count); - CommonParallel.For( - 0, - Count, - index => result[index] = -this[index]); + + for (var index = 0; index < Count; index++) + { + result.At(index, -At(index)); + } return result; } @@ -352,10 +350,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override int MaximumIndex() { var index = 0; - var max = this[index]; + var max = At(index); for (var i = 1; i < Count; i++) { - var test = this[i]; + var test = At(i); if (test > max) { index = i; @@ -373,10 +371,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override int MinimumIndex() { var index = 0; - var min = this[index]; + var min = At(index); for (var i = 1; i < Count; i++) { - var test = this[i]; + var test = At(i); if (test < min) { index = i; @@ -435,7 +433,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single var v = CreateVector(length); for (var index = 0; index < Count; index++) { - v[index] = Convert.ToSingle(randomDistribution.Sample()); + v.At(index, Convert.ToSingle(randomDistribution.Sample())); } return v; @@ -461,7 +459,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single var v = CreateVector(length); for (var index = 0; index < Count; index++) { - v[index] = Convert.ToSingle(randomDistribution.Sample()); + v.At(index, Convert.ToSingle(randomDistribution.Sample())); } return v; diff --git a/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedVectorTests.cs index b8ada682..282c4d52 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedVectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedVectorTests.cs @@ -54,19 +54,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex _data = (Complex[])data.Clone(); } - public override Complex this[int index] - { - get - { - return _data[index]; - } - - set - { - _data[index] = value; - } - } - public override Matrix CreateMatrix(int rows, int columns) { return new UserDefinedMatrix(rows, columns); @@ -77,6 +64,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex return new UserDefinedVector(size); } + protected internal override Complex At(int index) + { + return _data[index]; + } + + protected internal override void At(int index, Complex value) + { + _data[index] = value; + } } public class UserDefinedVectorTests : VectorTests diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedVectorTests.cs index 822223f1..aaf815ec 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedVectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedVectorTests.cs @@ -47,19 +47,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 _data = (Complex32[])data.Clone(); } - public override Complex32 this[int index] - { - get - { - return _data[index]; - } - - set - { - _data[index] = value; - } - } - public override Matrix CreateMatrix(int rows, int columns) { return new UserDefinedMatrix(rows, columns); @@ -69,6 +56,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 { return new UserDefinedVector(size); } + + protected internal override Complex32 At(int index) + { + return _data[index]; + } + + protected internal override void At(int index, Complex32 value) + { + _data[index] = value; + } } public class UserDefinedVectorTests : VectorTests diff --git a/src/UnitTests/LinearAlgebraTests/Double/UserDefinedVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Double/UserDefinedVectorTests.cs index 37c51450..1188561b 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/UserDefinedVectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/UserDefinedVectorTests.cs @@ -46,19 +46,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double _data = (double[])data.Clone(); } - public override double this[int index] - { - get - { - return _data[index]; - } - - set - { - _data[index] = value; - } - } - public override Matrix CreateMatrix(int rows, int columns) { return new UserDefinedMatrix(rows, columns); @@ -68,6 +55,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double { return new UserDefinedVector(size); } + + protected internal override double At(int index) + { + return _data[index]; + } + + protected internal override void At(int index, double value) + { + _data[index] = value; + } } public class UserDefinedVectorTests : VectorTests diff --git a/src/UnitTests/LinearAlgebraTests/Single/UserDefinedVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Single/UserDefinedVectorTests.cs index 7b484eef..7bb6ceaa 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/UserDefinedVectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/UserDefinedVectorTests.cs @@ -55,19 +55,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single _data = (float[])data.Clone(); } - public override float this[int index] - { - get - { - return _data[index]; - } - - set - { - _data[index] = value; - } - } - public override Matrix CreateMatrix(int rows, int columns) { return new UserDefinedMatrix(rows, columns); @@ -77,6 +64,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single { return new UserDefinedVector(size); } + + protected internal override float At(int index) + { + return _data[index]; + } + + protected internal override void At(int index, float value) + { + _data[index] = value; + } } public class UserDefinedVectorTests : VectorTests