diff --git a/src/Numerics/LinearAlgebra/Double/SparseVector.cs b/src/Numerics/LinearAlgebra/Double/SparseVector.cs index fa2d00dc..1ff9cdf4 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseVector.cs @@ -33,6 +33,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double using System; using System.Collections.Generic; using System.Globalization; + using Distributions; using NumberTheory; using Properties; using Threading; @@ -264,6 +265,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double return new SparseVector(size); } + public void Clear() + { + NonZerosCount = 0; + } + /// /// Copies the values of this vector into the target vector. /// @@ -425,13 +431,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double } else if (alpha == -1.0) { - NonZerosCount = 0; // Vector is subtracted from itself + Clear(); // Vector is subtracted from itself + return; } else { for (int i = 0; i < this.NonZerosCount; i++) { - this[other.NonZeroIndices[i]] += alpha * this.NonZeroValues[i]; + this.NonZeroValues[i] += alpha * this.NonZeroValues[i]; } } } @@ -727,7 +734,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double } if (scalar == 0) { - NonZerosCount = 0; // Set array empty + Clear(); // Set array empty + return; } Control.LinearAlgebraProvider.ScaleArray(scalar, this.NonZeroValues); } @@ -760,7 +768,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double } return result; } - /// /// Multiplies a vector with a scalar. /// @@ -844,58 +851,410 @@ namespace MathNet.Numerics.LinearAlgebra.Double ret.Multiply(1.0 / rightSide); return ret; } - #endregion - #region Vector Norms + /// + /// Returns the index of the absolute minimum element. + /// + /// The index of absolute minimum element. + public override int AbsoluteMinimumIndex() + { + if (this.NonZerosCount == 0) // No non-zero elements. Return 0 + return 0; + + var index = 0; + var min = Math.Abs(this.NonZeroValues[index]); + for (var i = 1; i < this.NonZerosCount; i++) + { + var test = Math.Abs(this.NonZeroValues[i]); + if (test < min) + { + index = i; + min = test; + } + } + + return this.NonZeroIndices[index]; + } + + /// + /// Creates a vector containing specified elements. + /// + /// The first element to begin copying from. + /// The number of elements to copy. + /// A vector containing a copy of the specified elements. + /// If is not positive or + /// greater than or equal to the size of the vector. + /// If + is greater than or equal to the size of the vector. + /// + /// If is not positive. + public override Vector SubVector(int index, int length) + { + if (index < 0 || index >= this.Count) + { + throw new ArgumentOutOfRangeException("index"); + } + + if (length <= 0) + { + throw new ArgumentOutOfRangeException("length"); + } + + if (index + length > this.Count) + { + throw new ArgumentOutOfRangeException("length"); + } + + var result = new SparseVector(length); + for (int i = index; i < index + length; i++) + result[i - index] = this[i]; + + return result; + } + + /// + /// Set the values of this vector to the given values. + /// + /// The array containing the values to use. + /// If is . + /// If is not the same size as this vector. + public override void SetValues(double[] values) + { + if (values == null) + { + throw new ArgumentNullException("values"); + } + + if (values.Length != this.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "values"); + } + + for (int i = 0; i < values.Length; i++ ) + this[i] = values[i]; + } + /// + /// Returns the index of the absolute maximum element. + /// + /// The index of absolute maximum element. + public override int MaximumIndex() + { + if (this.NonZerosCount == 0) + return 0; + + var index = 0; + var max = this.NonZeroValues[0]; + for (var i = 1; i < this.NonZerosCount; i++) + { + if (max < this.NonZeroValues[i]) + { + index = i; + max = this.NonZeroValues[i]; + } + } + + return this.NonZeroIndices[index]; + } + /// + /// Returns the index of the minimum element. + /// + /// The index of minimum element. + public override int MinimumIndex() + { + if (this.NonZerosCount == 0) + return 0; + + var index = 0; + var min = this.NonZeroValues[0]; + for (var i = 1; i < this.NonZerosCount; i++) + { + if (min > this.NonZeroValues[i]) + { + index = i; + min = this.NonZeroValues[i]; + } + } + + return this.NonZeroIndices[index]; + } /// - /// Euclidean Norm also known as 2-Norm. + /// Computes the sum of the vector's elements. /// - /// Scalar ret = sqrt(sum(this[i]^2)) - public override double Norm() + /// The sum of the vector's elements. + public override double Sum() { - var sum = 0.0; + double result = 0; + for (var i = 0; i < this.NonZerosCount; i++) + { + result += this.NonZeroValues[i]; + } + + return result; + } + /// + /// Computes the sum of the absolute value of the vector's elements. + /// + /// The sum of the absolute value of the vector's elements. + public override double SumMagnitudes() + { + double result = 0; + for (var i = 0; i < this.NonZerosCount; i++) + { + result += Math.Abs(this.NonZeroValues[i]); + } + + return result; + } + /// + /// Pointwise multiplies this vector with another vector. + /// + /// The vector to pointwise multiply with this one. + /// If the other vector is . + /// If this vector and are not the same size. + public override void PointWiseMultiply(Vector other) + { + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (this.Count != other.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + } + + // We cannot iterate using NonZeroCount because the value may be changed (if multiply by 0) for (var i = 0; i < this.Count; i++) { - sum = SpecialFunctions.Hypotenuse(sum, this[i]); + this[i] *= other[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. + /// 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 PointWiseMultiply(Vector other, Vector result) + { + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (other == null) + { + throw new ArgumentNullException("other"); } - return sum; + if (this.Count != other.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + } + + if (this.Count != result.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); + } + + if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) + { + var tmp = result.CreateVector(result.Count); + this.PointWiseMultiply(other, tmp); + tmp.CopyTo(result); + } + else + { + this.CopyTo(result); + result.PointWiseMultiply(other); + } } /// - /// 1-Norm also known as Manhattan Norm or Taxicab Norm. + /// Pointwise divide this vector with another vector. /// - /// Scalar ret = sum(abs(this[i])) - public override double Norm1() + /// The vector to pointwise divide this one by. + /// If the other vector is . + /// If this vector and are not the same size. + public override void PointWiseDivide(Vector other) { - return CommonParallel.Aggregate( - 0, - this.NonZerosCount, - index => Math.Abs(this.NonZeroValues[index])); + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (this.Count != other.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + } + + // base implementation iterates though all elements, but we need only take non-zeros + for (var i = 0; i < this.NonZerosCount; i++) + { + this[this.NonZeroIndices[i]] /= other[this.NonZeroIndices[i]]; + } } /// - /// Computes the p-Norm. + /// Pointwise divide this vector with another vector and stores the result into the result vector. /// - /// The p value. - /// Scalar ret = (sum(abs(this[i])^p))^(1/p) - public override double NormP(int p) + /// The vector to pointwise divide this one by. + /// The vector to store the result of the pointwise division. + /// 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 PointWiseDivide(Vector other, Vector result) { - if (1 > p) + if (result == null) { - throw new ArgumentOutOfRangeException("p"); + throw new ArgumentNullException("result"); + } + + if (other == null) + { + throw new ArgumentNullException("other"); } - if (1 == p) + if (this.Count != other.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); + } + + if (this.Count != result.Count) { - return this.Norm1(); + throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } - if (2 == p) + if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) { - return this.Norm(); + var tmp = result.CreateVector(result.Count); + this.PointWiseDivide(other, tmp); + tmp.CopyTo(result); + } + else + { + this.CopyTo(result); + result.PointWiseDivide(other); + } + } + + /// + /// Outer product of two vectors + /// + /// First vector + /// Second vector + /// 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) + { + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (v == null) + { + throw new ArgumentNullException("v"); + } + + throw new NotImplementedException(); + //var matrix = new DenseMatrix(u.Count, v.Count); + //CommonParallel.For( + // 0, + // u.Count, + // i => + // { + // for (int j = 0; j < v.Count; j++) + // { + // matrix.At(i, j, u.Data[i] * v.Data[j]); + // } + // }); + //return matrix; + } + + /// + /// Generates a vector with random elements + /// + /// Number of elements in the vector. + /// Continuous Random Distribution or Source + /// + /// A vector with n-random elements distributed according + /// to the specified random distribution. + /// + /// If the length vector is non poisitive. + public override Vector Random(int length, IContinuousDistribution randomDistribution) + { + if (length < 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + var v = (SparseVector)this.CreateVector(length); + for (var index = 0; index < v.Count; index++) + { + v[index] = randomDistribution.Sample(); + } + + return v; + } + /// + /// Generates a vector with random elements + /// + /// Number of elements in the vector. + /// Continuous Random Distribution or Source + /// + /// A vector with n-random elements distributed according + /// to the specified random distribution. + /// + /// If the n vector is non poisitive. + public override Vector Random(int length, IDiscreteDistribution randomDistribution) + { + if (length < 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + var v = (SparseVector)this.CreateVector(length); + for (var index = 0; index < v.Count; index++) + { + v[index] = randomDistribution.Sample(); + } + + return v; + } + + /// + /// Tensor Product (Outer) of this and another vector. + /// + /// The vector to operate on. + /// + /// Matrix M[i,j] = this[i] * v[j]. + /// + /// + public Matrix TensorMultiply(SparseVector v) + { + return OuterProduct(this, v); + } + #endregion + + #region Vector Norms + /// + /// Computes the p-Norm. + /// + /// The p value. + /// Scalar ret = (sum(abs(this[i])^p))^(1/p) + public override double NormP(int p) + { + if (1 > p) + { + throw new ArgumentOutOfRangeException("p"); } var sum = CommonParallel.Aggregate( @@ -1122,7 +1481,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double itemIndex = ~itemIndex; //Index where to put new value // Check if the storage needs to be increased - if (NonZerosCount == NonZeroValues.Length) + if ((NonZerosCount == NonZeroValues.Length) && (NonZerosCount < Count)) { // Value and Indices arrays are completely full so we increase the size int size = Math.Min(NonZeroValues.Length + GrowthSize(), Count); diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj index 1c4772bc..ce5b8d66 100644 --- a/src/Silverlight/Silverlight.csproj +++ b/src/Silverlight/Silverlight.csproj @@ -242,6 +242,9 @@ LinearAlgebra\Double\Matrix.cs + + LinearAlgebra\Double\SparseVector.cs + LinearAlgebra\Double\Vector.cs diff --git a/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs b/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs index 75a96424..e35a064c 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs @@ -180,7 +180,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double } [Test] - public void CanCallUnaryNegationOperatorOnDenseVector() + public void CanCallUnaryNegationOperatorOnSparseVector() { var vector = new SparseVector(_data); var other = -vector; @@ -257,6 +257,39 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double Assert.AreEqual(_data[i] / 2.0, vector[i]); } } + + [Test] + public void CanCalculateOuterProductForSparseVector() + { + var vector1 = this.CreateVector(this._data); + var vector2 = this.CreateVector(this._data); + Matrix m = Vector.OuterProduct(vector1, vector2); + for (var i = 0; i < vector1.Count; i++) + { + for (var j = 0; j < vector2.Count; j++) + { + Assert.AreEqual(m[i, j], vector1[i] * vector2[j]); + } + } + } + + [Test] + [ExpectedArgumentNullException] + public void OuterProducForSparseVectortWithFirstParameterNullShouldThrowException() + { + SparseVector vector1 = null; + var vector2 = this.CreateVector(this._data); + Vector.OuterProduct(vector1, vector2); + } + + [Test] + [ExpectedArgumentNullException] + public void OuterProductForSparseVectorWithSecondParameterNullShouldThrowException() + { + var vector1 = this.CreateVector(this._data); + SparseVector vector2 = null; + Vector.OuterProduct(vector1, vector2); + } #endregion [Test] @@ -365,5 +398,22 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double var vector = new SparseVector(data); } + + [Test] + public void PointWiseMultiplySparseVector() + { + var zeroArray = new[] { 0.0, 1.0, 0.0, 1.0, 0.0 }; + var vector1 = new SparseVector(this._data); + var vector2 = new SparseVector(zeroArray); + var result = new SparseVector(vector1.Count); + + vector1.PointWiseMultiply(vector2, result); + + for (var i = 0; i < vector1.Count; i++) + { + Assert.AreEqual(this._data[i] * zeroArray[i], result[i]); + } + Assert.AreEqual(2, result.NonZerosCount); + } } }