diff --git a/src/MathNet.Numerics.5.1.ReSharper b/src/MathNet.Numerics.5.1.ReSharper index 1b3673de..37a21618 100644 --- a/src/MathNet.Numerics.5.1.ReSharper +++ b/src/MathNet.Numerics.5.1.ReSharper @@ -19,7 +19,11 @@ Matlab Matlab Matlab Endian -indices +indices +&lt +&gt +Frobenius +Pointwise diff --git a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs index 983c7cb2..ba9d78d1 100644 --- a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs @@ -39,7 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// A vector using dense storage. /// - public class DenseVector : Vector + public class DenseVector : Vector { /// /// Initializes a new instance of the class with a given size. @@ -767,7 +767,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// Returns the value of the absolute minimum element. /// /// The value of the absolute minimum element. - public override double AbsoluteMinimum() + public override Complex AbsoluteMinimum() { return Data[AbsoluteMinimumIndex()].Magnitude; } @@ -776,7 +776,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// Returns the value of the absolute maximum element. /// /// The value of the absolute maximum element. - public override double AbsoluteMaximum() + public override Complex AbsoluteMaximum() { return Data[AbsoluteMaximumIndex()].Magnitude; } @@ -869,28 +869,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The sum of the vector's elements. public override Complex Sum() { - var result = Complex.Zero; - for (var i = 0; i < Count; i++) - { - result += Data[i]; - } - - return result; + return CommonParallel.Aggregate( + 0, + Count, + i => Data[i]); } /// /// 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() + public override Complex SumMagnitudes() { - double result = 0; - for (var i = 0; i < Count; i++) - { - result += Data[i].Magnitude; - } - - return result; + return CommonParallel.Aggregate( + 0, + Count, + i => Data[i].Magnitude); } /// @@ -1101,58 +1095,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex 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 n vector is non positive. - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (DenseVector)CreateVector(length); - for (var index = 0; index < v.Data.Length; index++) - { - v.Data[index] = new Complex(randomDistribution.Sample(), 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 positive. - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (DenseVector)CreateVector(length); - for (var index = 0; index < v.Data.Length; index++) - { - v.Data[index] = new Complex(randomDistribution.Sample(), randomDistribution.Sample()); - } - - return v; - } - /// /// Outer product of this and another vector. /// @@ -1171,7 +1113,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// The p value. /// Scalar ret = (sum(abs(this[i])^p))^(1/p) - public override double Norm(double p) + public override Complex Norm(double p) { if (p < 0.0) { @@ -1208,34 +1150,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return Math.Pow(sum, 1.0 / p); } - /// - /// Normalizes this vector to a unit vector with respect to the p-norm. - /// - /// - /// The p value. - /// - /// - /// This vector normalized to a unit vector with respect to the p-norm. - /// - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(1.0 / norm, clone); - - return clone; - } - #region Parse Functions /// @@ -1395,24 +1309,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex #endregion - /// - /// Returns the index of the absolute maximum element. - /// - /// The index of absolute maximum element. - public override int MaximumIndex() - { - throw new NotSupportedException(); - } - - /// - /// Returns the index of the minimum element. - /// - /// The index of minimum element. - public override int MinimumIndex() - { - throw new NotSupportedException(); - } - /// /// Resets all values to zero. /// @@ -1457,61 +1353,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex index => otherVector.Data[index] = Data[index].Conjugate()); } } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override Complex AddT(Complex val1, Complex val2) - { - return val1 + val2; - } - - /// - /// Subtract two values T-T - /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected sealed override Complex SubtractT(Complex val1, Complex val2) - { - return val1 - val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - /// - /// Divide two values T/T - /// - /// Left operand value - /// Right operand value - /// Result of divide - protected sealed override Complex DivideT(Complex val1, Complex val2) - { - return val1 / val2; - } - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected sealed override double AbsoluteT(Complex val1) - { - return val1.Magnitude; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs index 7dd708db..6cd59fde 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs @@ -286,7 +286,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative Complex eta = 0; double theta = 0; - var tau = startNorm; + var tau = startNorm.Real; Complex rho = tau * tau; // Calculate the initial values for v @@ -342,7 +342,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative yinternal.Add(temp, d); // theta = ||pseudoResiduals||_2 / tau - theta = pseudoResiduals.Norm(2) / tau; + theta = pseudoResiduals.Norm(2).Real / tau; var c = 1 / Math.Sqrt(1 + (theta * theta)); // tau = tau * theta * c diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Preconditioners/Ilutp.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Preconditioners/Ilutp.cs index 32bdb7aa..a2bd3141 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Preconditioners/Ilutp.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Preconditioners/Ilutp.cs @@ -442,7 +442,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners // { // w(j) = 0 // } - if (workVector[j].Magnitude <= _dropTolerance * vectorNorm) + if (workVector[j].Magnitude <= _dropTolerance * vectorNorm.Real) { workVector[j] = 0.0; } diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/DivergenceStopCriterium.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/DivergenceStopCriterium.cs index ec8515d4..d34600b0 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/DivergenceStopCriterium.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/DivergenceStopCriterium.cs @@ -252,7 +252,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium // Store the infinity norms of both the solution and residual vectors // These values will be used to calculate the relative drop in residuals later on. - _residualHistory[_residualHistory.Length - 1] = residualVector.Norm(Double.PositiveInfinity); + _residualHistory[_residualHistory.Length - 1] = residualVector.Norm(Double.PositiveInfinity).Real; // Check if we have NaN's. If so we've gone way beyond normal divergence. // Stop the iteration. diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/FailureStopCriterium.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/FailureStopCriterium.cs index 9c61ac38..f6f6fc39 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/FailureStopCriterium.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/FailureStopCriterium.cs @@ -110,7 +110,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium var residualNorm = residualVector.Norm(Double.PositiveInfinity); var solutionNorm = solutionVector.Norm(Double.PositiveInfinity); - if (Double.IsNaN(solutionNorm) || Double.IsNaN(residualNorm)) + if (Double.IsNaN(solutionNorm.Real) || Double.IsNaN(residualNorm.Real)) { SetStatusToFailed(); } diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/ResidualStopCriterium.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/ResidualStopCriterium.cs index 82e8e2e4..153e1806 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/ResidualStopCriterium.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/ResidualStopCriterium.cs @@ -263,12 +263,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium // Check the residuals by calculating: // ||r_i|| <= stop_tol * ||b|| - var stopCriterium = ComputeStopCriterium(sourceVector.Norm(Double.PositiveInfinity)); + var stopCriterium = ComputeStopCriterium(sourceVector.Norm(Double.PositiveInfinity).Real); // First check that we have real numbers not NaN's. // NaN's can occur when the iterative process diverges so we // stop if that is the case. - if (double.IsNaN(stopCriterium) || double.IsNaN(residualNorm)) + if (double.IsNaN(stopCriterium) || double.IsNaN(residualNorm.Real)) { _iterationCount = 0; SetStatusToDiverged(); @@ -278,7 +278,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium // ||r_i|| <= stop_tol * ||b|| // Stop the calculation if it's clearly smaller than the tolerance var decimalMagnitude = Math.Abs(stopCriterium.Magnitude()) + 1; - if (residualNorm.IsSmallerWithDecimalPlaces(stopCriterium, decimalMagnitude)) + if (residualNorm.Real.IsSmallerWithDecimalPlaces(stopCriterium, decimalMagnitude)) { if (_lastIteration <= iterationNumber) { diff --git a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs index aa18a3a1..cb4fe7ca 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs @@ -30,7 +30,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex using System.Collections.Generic; using System.Linq; using System.Numerics; - using Distributions; using Generic; using NumberTheory; using Properties; @@ -39,7 +38,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// A vector with sparse storage. /// - public class SparseVector : Vector + public class SparseVector : Vector { /// /// Lock object for the indexer. @@ -581,11 +580,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var sparseother = other as SparseVector; if (sparseother == null) { - sparse.AddScaledSparseVector(Complex.One, sparseother); + base.Add(other, result); } else { - base.Add(other, result); + CopyTo(result); + sparse.AddScaledSparseVector(Complex.One, sparseother); } } } @@ -753,11 +753,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var sparseother = other as SparseVector; if (sparseother == null) { - sparse.AddScaledSparseVector(-Complex.One, sparseother); + base.Subtract(other, result); } else { - base.Subtract(other, result); + CopyTo(result); + sparse.AddScaledSparseVector(-Complex.One, sparseother); } } } @@ -1075,7 +1076,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// 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() + public override Complex SumMagnitudes() { double result = 0; for (var i = 0; i < NonZerosCount; i++) @@ -1153,58 +1154,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex 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 positive. - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (SparseVector)CreateVector(length); - for (var index = 0; index < v.Count; index++) - { - v[index] = new Complex(randomDistribution.Sample(), 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 positive. - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (SparseVector)CreateVector(length); - for (var index = 0; index < v.Count; index++) - { - v[index] = new Complex(randomDistribution.Sample(), randomDistribution.Sample()); - } - - return v; - } - /// /// Outer product of this and another vector. /// @@ -1227,7 +1176,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// The p value. /// Scalar ret = (sum(abs(this[i])^p))^(1/p) - public override double Norm(double p) + public override Complex Norm(double p) { if (1 > p) { @@ -1257,33 +1206,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return Math.Pow(sum, 1.0 / p); } - /// - /// Normalizes this vector to a unit vector with respect to the p-norm. - /// - /// - /// The p value. - /// - /// - /// This vector normalized to a unit vector with respect to the p-norm. - /// - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(1.0 / norm, clone); - - return clone; - } #endregion #region Parse Functions @@ -1610,79 +1532,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } #endregion - - /// - /// Returns the index of the absolute maximum element. - /// - /// The index of absolute maximum element. - public override int MaximumIndex() - { - throw new NotSupportedException(); - } - - /// - /// Returns the index of the minimum element. - /// - /// The index of minimum element. - public override int MinimumIndex() - { - throw new NotSupportedException(); - } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override Complex AddT(Complex val1, Complex val2) - { - return val1 + val2; - } - - /// - /// Subtract two values T-T - /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected sealed override Complex SubtractT(Complex val1, Complex val2) - { - return val1 - val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - /// - /// Divide two values T/T - /// - /// Left operand value - /// Right operand value - /// Result of divide - protected sealed override Complex DivideT(Complex val1, Complex val2) - { - return val1 / val2; - } - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected sealed override double AbsoluteT(Complex val1) - { - return val1.Magnitude; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex/Vector.cs b/src/Numerics/LinearAlgebra/Complex/Vector.cs new file mode 100644 index 00000000..809fe1ed --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex/Vector.cs @@ -0,0 +1,446 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// Copyright (c) 2009-2010 Math.NET +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.LinearAlgebra.Complex +{ + using System; + using System.Numerics; + using Distributions; + using Generic; + using Properties; + using Threading; + + /// + /// Complex version of the class. + /// + public abstract class Vector : Vector + { + /// + /// Initializes a new instance of the Vector class. + /// Constructs a Vector with the given size. + /// + /// + /// The size of the Vector to construct. + /// + /// + /// If is less than one. + /// + protected Vector(int size) : base(size) + { + } + + /// + /// 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. + /// + protected override void DoAdd(Complex scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] + scalar); + } + + /// + /// 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. + /// + protected override void DoAdd(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] + other[index]); + } + + /// + /// 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. + /// + protected override void DoSubtract(Complex scalar, Vector result) + { + 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. + /// + protected override void DoSubtract(Vector other, Vector result) + { + CopyTo(result); + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] - other[index]); + } + + /// + /// Multiplies a scalar to each element of the vector and stores the result in the result vector. + /// + /// + /// The scalar to multiply. + /// + /// + /// The vector to store the result of the multiplication. + /// + protected override void DoMultiply(Complex scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] * scalar); + } + + /// + /// Divides each element of the vector by a scalar and stores the result in the result vector. + /// + /// + /// The scalar to divide with. + /// + /// + /// The vector to store the result of the division. + /// + protected override void DoDivide(Complex scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] / scalar); + } + + /// + /// 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 DoPointwiseMultiply(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] * other[index]); + } + + /// + /// 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 override void DoPointwiseDivide(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] / other[index]); + } + + /// + /// Computes the dot product between this vector and another vector. + /// + /// + /// The other vector to add. + /// + /// s + /// The result of the addition. + /// + protected override Complex DoDotProduct(Vector other) + { + return CommonParallel.Aggregate( + 0, + Count, + i => this[i] * other[i]); + } + + /// + /// Returns the value of the absolute minimum element. + /// + /// The value of the absolute minimum element. + public override Complex AbsoluteMinimum() + { + return this[AbsoluteMinimumIndex()].Magnitude; + } + + /// + /// Returns the index of the absolute minimum element. + /// + /// The index of absolute minimum element. + public override int AbsoluteMinimumIndex() + { + var index = 0; + var min = this[index].Magnitude; + for (var i = 1; i < Count; i++) + { + var test = this[i].Magnitude; + if (test < min) + { + index = i; + min = test; + } + } + + return index; + } + + /// + /// Returns the value of the absolute maximum element. + /// + /// The value of the absolute maximum element. + public override Complex AbsoluteMaximum() + { + return this[AbsoluteMaximumIndex()].Magnitude; + } + + /// + /// Returns the index of the absolute maximum element. + /// + /// The index of absolute maximum element. + public override int AbsoluteMaximumIndex() + { + var index = 0; + var max = this[index].Magnitude; + for (var i = 1; i < Count; i++) + { + var test = this[i].Magnitude; + if (test > max) + { + index = i; + max = test; + } + } + + return index; + } + + /// + /// Computes the sum of the vector's elements. + /// + /// The sum of the vector's elements. + public override Complex Sum() + { + return CommonParallel.Aggregate( + 0, + Count, + i => this[i]); + } + + /// + /// 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 Complex SumMagnitudes() + { + return CommonParallel.Aggregate( + 0, + Count, + i => this[i].Magnitude); + } + + /// + /// Computes the p-Norm. + /// + /// + /// The p value. + /// + /// + /// Scalar ret = (sum(abs(this[i])^p))^(1/p) + /// + public override Complex Norm(double p) + { + if (p < 0.0) + { + throw new ArgumentOutOfRangeException("p"); + } + + if (double.IsPositiveInfinity(p)) + { + return CommonParallel.Select( + 0, + Count, + (index, localData) => Math.Max(localData, this[index].Magnitude), + Math.Max); + } + + var sum = CommonParallel.Aggregate( + 0, + Count, + index => Math.Pow(this[index].Magnitude, p)); + + return Math.Pow(sum, 1.0 / p); + } + + /// + /// Conjugates vector and save result to + /// + /// Target vector + protected override void DoConjugate(Vector target) + { + CopyTo(target); + CommonParallel.For( + 0, + Count, + index => target[index] = this[index].Conjugate()); + } + + /// + /// Returns a negated vector. + /// + /// + /// The negated vector. + /// + /// + /// Added as an alternative to the unary negation operator. + /// + public override Vector Negate() + { + var result = CreateVector(Count); + CommonParallel.For( + 0, + Count, + index => result[index] = -this[index]); + + return result; + } + + /// + /// 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-positive. + public override Vector Random(int length, IContinuousDistribution randomDistribution) + { + if (length < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + var vector = CreateVector(length); + for (var index = 0; index < length; index++) + { + vector[index] = new Complex(randomDistribution.Sample(), randomDistribution.Sample()); + } + + return vector; + } + + /// + /// 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 not positive. + public override Vector Random(int length, IDiscreteDistribution randomDistribution) + { + if (length < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + var vector = CreateVector(length); + for (var index = 0; index < length; index++) + { + vector[index] = new Complex(randomDistribution.Sample(), randomDistribution.Sample()); + } + + return vector; + } + + /// + /// Returns the index of the absolute maximum element. + /// + /// The index of absolute maximum element. + public override int MaximumIndex() + { + throw new NotSupportedException(); + } + + /// + /// Returns the index of the minimum element. + /// + /// The index of minimum element. + public override int MinimumIndex() + { + throw new NotSupportedException(); + } + + /// + /// Normalizes this vector to a unit vector with respect to the p-norm. + /// + /// + /// The p value. + /// + /// + /// This vector normalized to a unit vector with respect to the p-norm. + /// + public override Vector Normalize(double p) + { + if (p < 0.0) + { + throw new ArgumentOutOfRangeException("p"); + } + + var norm = Norm(p); + var clone = Clone(); + if (norm.Real == 0.0) + { + return clone; + } + + clone.Multiply(1.0 / norm, clone); + + return clone; + } + } +} diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs index 6c778ba4..35c8f2a7 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs @@ -39,7 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// A vector using dense storage. /// - public class DenseVector : Vector + public class DenseVector : Vector { /// /// Initializes a new instance of the class with a given size. @@ -767,7 +767,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// Returns the value of the absolute minimum element. /// /// The value of the absolute minimum element. - public override double AbsoluteMinimum() + public override Complex32 AbsoluteMinimum() { return Data[AbsoluteMinimumIndex()].Magnitude; } @@ -776,7 +776,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// Returns the value of the absolute maximum element. /// /// The value of the absolute maximum element. - public override double AbsoluteMaximum() + public override Complex32 AbsoluteMaximum() { return Data[AbsoluteMaximumIndex()].Magnitude; } @@ -869,28 +869,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The sum of the vector's elements. public override Complex32 Sum() { - var result = Complex32.Zero; - for (var i = 0; i < Count; i++) - { - result += Data[i]; - } - - return result; + return CommonParallel.Aggregate( + 0, + Count, + i => Data[i]); } /// /// 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() + public override Complex32 SumMagnitudes() { - double result = 0; - for (var i = 0; i < Count; i++) - { - result += Data[i].Magnitude; - } - - return result; + return CommonParallel.Aggregate( + 0, + Count, + i => Data[i].Magnitude); } /// @@ -1171,7 +1165,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// The p value. /// Scalar ret = (sum(abs(this[i])^p))^(1/p) - public override double Norm(double p) + public override Complex32 Norm(double p) { if (p < 0.0) { @@ -1197,7 +1191,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 0, Count, (index, localData) => Math.Max(localData, Data[index].Magnitude), - Math.Max); + Common.Max); } var sum = CommonParallel.Aggregate( @@ -1205,35 +1199,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 Count, index => Math.Pow(Data[index].Magnitude, p)); - return Math.Pow(sum, 1.0 / p); - } - - /// - /// Normalizes this vector to a unit vector with respect to the p-norm. - /// - /// - /// The p value. - /// - /// - /// This vector normalized to a unit vector with respect to the p-norm. - /// - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(Complex32.One / (float)norm, clone); - - return clone; + return (float)Math.Pow(sum, 1.0 / p); } #region Parse Functions @@ -1395,24 +1361,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 #endregion - /// - /// Returns the index of the absolute maximum element. - /// - /// The index of absolute maximum element. - public override int MaximumIndex() - { - throw new NotSupportedException(); - } - - /// - /// Returns the index of the minimum element. - /// - /// The index of minimum element. - public override int MinimumIndex() - { - throw new NotSupportedException(); - } - /// /// Resets all values to zero. /// @@ -1457,61 +1405,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 index => otherVector.Data[index] = Data[index].Conjugate()); } } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override Complex32 AddT(Complex32 val1, Complex32 val2) - { - return val1 + val2; - } - - /// - /// Subtract two values T-T - /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected sealed override Complex32 SubtractT(Complex32 val1, Complex32 val2) - { - return val1 - val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - /// - /// Divide two values T/T - /// - /// Left operand value - /// Right operand value - /// Result of divide - protected sealed override Complex32 DivideT(Complex32 val1, Complex32 val2) - { - return val1 / val2; - } - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected sealed override double AbsoluteT(Complex32 val1) - { - return val1.Magnitude; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserGramSchmidt.cs index 385e35e9..662a0c7b 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserGramSchmidt.cs @@ -70,7 +70,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization for (var k = 0; k < MatrixQ.ColumnCount; k++) { - var norm = (float)MatrixQ.Column(k).Norm(2); + var norm = MatrixQ.Column(k).Norm(2).Real; if (norm == 0.0f) { throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient); diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs index f6e671a6..0b6fbf19 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs @@ -652,7 +652,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative result.Add(orthogonalMatrix.Column(i)); // Normalize the result vector - result[i].Multiply(1 / (float)result[i].Norm(2), result[i]); + result[i].Multiply(1 / result[i].Norm(2).Real, result[i]); } return result; diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs index e7c8f1df..50bf4fd2 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs @@ -286,7 +286,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative Complex32 eta = 0; float theta = 0; - var tau = (float)startNorm; + var tau = startNorm.Real; Complex32 rho = tau * tau; // Calculate the initial values for v @@ -342,7 +342,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative yinternal.Add(temp, d); // theta = ||pseudoResiduals||_2 / tau - theta = (float)pseudoResiduals.Norm(2) / tau; + theta = pseudoResiduals.Norm(2).Real / tau; var c = 1 / (float)Math.Sqrt(1 + (theta * theta)); // tau = tau * theta * c diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Preconditioners/Ilutp.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Preconditioners/Ilutp.cs index f95d7a4e..f7769a3e 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Preconditioners/Ilutp.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Preconditioners/Ilutp.cs @@ -442,7 +442,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners // { // w(j) = 0 // } - if (workVector[j].Magnitude <= _dropTolerance * vectorNorm) + if (workVector[j].Magnitude <= _dropTolerance * vectorNorm.Real) { workVector[j] = 0.0f; } diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/DivergenceStopCriterium.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/DivergenceStopCriterium.cs index c4537d04..165e9575 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/DivergenceStopCriterium.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/DivergenceStopCriterium.cs @@ -252,7 +252,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium // Store the infinity norms of both the solution and residual vectors // These values will be used to calculate the relative drop in residuals later on. - _residualHistory[_residualHistory.Length - 1] = residualVector.Norm(Double.PositiveInfinity); + _residualHistory[_residualHistory.Length - 1] = residualVector.Norm(Double.PositiveInfinity).Real; // Check if we have NaN's. If so we've gone way beyond normal divergence. // Stop the iteration. diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/FailureStopCriterium.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/FailureStopCriterium.cs index 00be0d44..d3d1d0a7 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/FailureStopCriterium.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/FailureStopCriterium.cs @@ -110,7 +110,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium var residualNorm = residualVector.Norm(Double.PositiveInfinity); var solutionNorm = solutionVector.Norm(Double.PositiveInfinity); - if (Double.IsNaN(solutionNorm) || Double.IsNaN(residualNorm)) + if (Single.IsNaN(solutionNorm.Real) || Single.IsNaN(residualNorm.Real)) { SetStatusToFailed(); } diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/ResidualStopCriterium.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/ResidualStopCriterium.cs index 284f505b..906e6046 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/ResidualStopCriterium.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/ResidualStopCriterium.cs @@ -46,7 +46,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium /// /// The default value for the maximum value of the residual. /// - public const double DefaultMaximumResidual = 1e-6; + public const float DefaultMaximumResidual = 1e-6f; /// /// The default value for the minimum number of iterations. @@ -66,7 +66,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium /// /// The maximum value for the residual below which the calculation is considered converged. /// - private double _maximum; + private float _maximum; /// /// The minimum number of iterations for which the residual has to be below the maximum before @@ -102,7 +102,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium /// maximum residual and the default minimum number of iterations. /// /// The maximum value for the residual below which the calculation is considered converged. - public ResidualStopCriterium(double maximum) : this(maximum, DefaultMinimumIterationsBelowMaximum) + public ResidualStopCriterium(float maximum) : this(maximum, DefaultMinimumIterationsBelowMaximum) { } @@ -129,7 +129,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium /// The minimum number of iterations for which the residual has to be below the maximum before /// the calculation is considered converged. /// - public ResidualStopCriterium(double maximum, int minimumIterationsBelowMaximum) + public ResidualStopCriterium(float maximum, int minimumIterationsBelowMaximum) { if (maximum < 0) { @@ -150,7 +150,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium /// converged. /// /// Thrown if the Maximum is set to a negative value. - public double Maximum + public float Maximum { [DebuggerStepThrough] get @@ -259,16 +259,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium // Store the infinity norms of both the solution and residual vectors // These values will be used to calculate the relative drop in residuals // later on. - var residualNorm = residualVector.Norm(Double.PositiveInfinity); + var residualNorm = residualVector.Norm(float.PositiveInfinity); // Check the residuals by calculating: // ||r_i|| <= stop_tol * ||b|| - var stopCriterium = ComputeStopCriterium(sourceVector.Norm(Double.PositiveInfinity)); + var stopCriterium = ComputeStopCriterium(sourceVector.Norm(float.PositiveInfinity).Real); // First check that we have real numbers not NaN's. // NaN's can occur when the iterative process diverges so we // stop if that is the case. - if (double.IsNaN(stopCriterium) || double.IsNaN(residualNorm)) + if (float.IsNaN(stopCriterium) || float.IsNaN(residualNorm.Real)) { _iterationCount = 0; SetStatusToDiverged(); @@ -278,7 +278,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium // ||r_i|| <= stop_tol * ||b|| // Stop the calculation if it's clearly smaller than the tolerance var decimalMagnitude = Math.Abs(stopCriterium.Magnitude()) + 1; - if (residualNorm.IsSmallerWithDecimalPlaces(stopCriterium, decimalMagnitude)) + if (residualNorm.Real.IsSmallerWithDecimalPlaces(stopCriterium, decimalMagnitude)) { if (_lastIteration <= iterationNumber) { @@ -307,7 +307,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium /// /// Solution vector norm /// Criterium value - private double ComputeStopCriterium(double solutionNorm) + private float ComputeStopCriterium(float solutionNorm) { // This is criterium 1 from Templates for the solution of linear systems. // The problem with this criterium is that it's not limiting enough. For now diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs index 59253675..bc361a6e 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs @@ -29,7 +29,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 using System; using System.Collections.Generic; using System.Linq; - using Distributions; using Generic; using NumberTheory; using Numerics; @@ -39,7 +38,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// A vector with sparse storage. /// - public class SparseVector : Vector + public class SparseVector : Vector { /// /// Lock object for the indexer. @@ -581,11 +580,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var sparseother = other as SparseVector; if (sparseother == null) { - sparse.AddScaledSparseVector(Complex32.One, sparseother); + base.Add(other, result); } else { - base.Add(other, result); + CopyTo(result); + sparse.AddScaledSparseVector(Complex32.One, sparseother); } } } @@ -753,11 +753,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var sparseother = other as SparseVector; if (sparseother == null) { - sparse.AddScaledSparseVector(-Complex32.One, sparseother); + base.Subtract(other, result); } else { - base.Subtract(other, result); + CopyTo(result); + sparse.AddScaledSparseVector(-Complex32.One, sparseother); } } } @@ -1075,9 +1076,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// 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() + public override Complex32 SumMagnitudes() { - double result = 0; + var result = 0.0f; for (var i = 0; i < NonZerosCount; i++) { result += _nonZeroValues[i].Magnitude; @@ -1153,58 +1154,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 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 positive. - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (SparseVector)CreateVector(length); - for (var index = 0; index < v.Count; index++) - { - v[index] = new Complex32((float)randomDistribution.Sample(), (float)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 positive. - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (SparseVector)CreateVector(length); - for (var index = 0; index < v.Count; index++) - { - v[index] = new Complex32(randomDistribution.Sample(), randomDistribution.Sample()); - } - - return v; - } - /// /// Outer product of this and another vector. /// @@ -1227,7 +1176,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// The p value. /// Scalar ret = (sum(abs(this[i])^p))^(1/p) - public override double Norm(double p) + public override Complex32 Norm(double p) { if (1 > p) { @@ -1236,7 +1185,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 if (NonZerosCount == 0) { - return 0.0; + return Complex32.Zero; } if (2.0 == p) @@ -1246,7 +1195,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 if (Double.IsPositiveInfinity(p)) { - return CommonParallel.Select(0, NonZerosCount, (index, localData) => Math.Max(localData, _nonZeroValues[index].Magnitude), Math.Max); + return CommonParallel.Select(0, NonZerosCount, (index, localData) => Math.Max(localData, _nonZeroValues[index].Magnitude), Common.Max); } var sum = CommonParallel.Aggregate( @@ -1254,36 +1203,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 NonZerosCount, index => Math.Pow(_nonZeroValues[index].Magnitude, p)); - return Math.Pow(sum, 1.0 / p); + return (float)Math.Pow(sum, 1.0 / p); } - /// - /// Normalizes this vector to a unit vector with respect to the p-norm. - /// - /// - /// The p value. - /// - /// - /// This vector normalized to a unit vector with respect to the p-norm. - /// - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(Complex32.One / (float)norm, clone); - - return clone; - } #endregion #region Parse Functions @@ -1610,79 +1532,5 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } #endregion - - /// - /// Returns the index of the absolute maximum element. - /// - /// The index of absolute maximum element. - public override int MaximumIndex() - { - throw new NotSupportedException(); - } - - /// - /// Returns the index of the minimum element. - /// - /// The index of minimum element. - public override int MinimumIndex() - { - throw new NotSupportedException(); - } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override Complex32 AddT(Complex32 val1, Complex32 val2) - { - return val1 + val2; - } - - /// - /// Subtract two values T-T - /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected sealed override Complex32 SubtractT(Complex32 val1, Complex32 val2) - { - return val1 - val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - /// - /// Divide two values T/T - /// - /// Left operand value - /// Right operand value - /// Result of divide - protected sealed override Complex32 DivideT(Complex32 val1, Complex32 val2) - { - return val1 / val2; - } - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected sealed override double AbsoluteT(Complex32 val1) - { - return val1.Magnitude; - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Vector.cs b/src/Numerics/LinearAlgebra/Complex32/Vector.cs new file mode 100644 index 00000000..b6610cbc --- /dev/null +++ b/src/Numerics/LinearAlgebra/Complex32/Vector.cs @@ -0,0 +1,446 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// Copyright (c) 2009-2010 Math.NET +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.LinearAlgebra.Complex32 +{ + using System; + using Distributions; + using Generic; + using Properties; + using Threading; + using Complex32 = Numerics.Complex32; + + /// + /// Complex32 version of the class. + /// + public abstract class Vector : Vector + { + /// + /// Initializes a new instance of the Vector class. + /// Constructs a Vector with the given size. + /// + /// + /// The size of the Vector to construct. + /// + /// + /// If is less than one. + /// + protected Vector(int size) : base(size) + { + } + + /// + /// 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. + /// + protected override void DoAdd(Complex32 scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] + scalar); + } + + /// + /// 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. + /// + protected override void DoAdd(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] + other[index]); + } + + /// + /// 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. + /// + protected override void DoSubtract(Complex32 scalar, Vector result) + { + 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. + /// + protected override void DoSubtract(Vector other, Vector result) + { + CopyTo(result); + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] - other[index]); + } + + /// + /// Multiplies a scalar to each element of the vector and stores the result in the result vector. + /// + /// + /// The scalar to multiply. + /// + /// + /// The vector to store the result of the multiplication. + /// + protected override void DoMultiply(Complex32 scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] * scalar); + } + + /// + /// Divides each element of the vector by a scalar and stores the result in the result vector. + /// + /// + /// The scalar to divide with. + /// + /// + /// The vector to store the result of the division. + /// + protected override void DoDivide(Complex32 scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] / scalar); + } + + /// + /// 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 DoPointwiseMultiply(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] * other[index]); + } + + /// + /// 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 override void DoPointwiseDivide(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] / other[index]); + } + + /// + /// Computes the dot product between this vector and another vector. + /// + /// + /// The other vector to add. + /// + /// s + /// The result of the addition. + /// + protected override Complex32 DoDotProduct(Vector other) + { + return CommonParallel.Aggregate( + 0, + Count, + i => this[i] * other[i]); + } + + /// + /// Returns the value of the absolute minimum element. + /// + /// The value of the absolute minimum element. + public override Complex32 AbsoluteMinimum() + { + return this[AbsoluteMinimumIndex()].Magnitude; + } + + /// + /// Returns the index of the absolute minimum element. + /// + /// The index of absolute minimum element. + public override int AbsoluteMinimumIndex() + { + var index = 0; + var min = this[index].Magnitude; + for (var i = 1; i < Count; i++) + { + var test = this[i].Magnitude; + if (test < min) + { + index = i; + min = test; + } + } + + return index; + } + + /// + /// Returns the value of the absolute maximum element. + /// + /// The value of the absolute maximum element. + public override Complex32 AbsoluteMaximum() + { + return this[AbsoluteMaximumIndex()].Magnitude; + } + + /// + /// Returns the index of the absolute maximum element. + /// + /// The index of absolute maximum element. + public override int AbsoluteMaximumIndex() + { + var index = 0; + var max = this[index].Magnitude; + for (var i = 1; i < Count; i++) + { + var test = this[i].Magnitude; + if (test > max) + { + index = i; + max = test; + } + } + + return index; + } + + /// + /// Computes the sum of the vector's elements. + /// + /// The sum of the vector's elements. + public override Complex32 Sum() + { + return CommonParallel.Aggregate( + 0, + Count, + i => this[i]); + } + + /// + /// 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 Complex32 SumMagnitudes() + { + return CommonParallel.Aggregate( + 0, + Count, + i => this[i].Magnitude); + } + + /// + /// Computes the p-Norm. + /// + /// + /// The p value. + /// + /// + /// Scalar ret = (sum(abs(this[i])^p))^(1/p) + /// + public override Complex32 Norm(double p) + { + if (p < 0.0) + { + throw new ArgumentOutOfRangeException("p"); + } + + if (double.IsPositiveInfinity(p)) + { + return CommonParallel.Select( + 0, + Count, + (index, localData) => Math.Max(localData, this[index].Magnitude), + Common.Max); + } + + var sum = CommonParallel.Aggregate( + 0, + Count, + index => Math.Pow(this[index].Magnitude, p)); + + return (float)Math.Pow(sum, 1.0 / p); + } + + /// + /// Conjugates vector and save result to + /// + /// Target vector + protected override void DoConjugate(Vector target) + { + CopyTo(target); + CommonParallel.For( + 0, + Count, + index => target[index] = this[index].Conjugate()); + } + + /// + /// Returns a negated vector. + /// + /// + /// The negated vector. + /// + /// + /// Added as an alternative to the unary negation operator. + /// + public override Vector Negate() + { + var result = CreateVector(Count); + CommonParallel.For( + 0, + Count, + index => result[index] = -this[index]); + + return result; + } + + /// + /// Returns the index of the absolute maximum element. + /// + /// The index of absolute maximum element. + public override int MaximumIndex() + { + throw new NotSupportedException(); + } + + /// + /// Returns the index of the minimum element. + /// + /// The index of minimum element. + public override int MinimumIndex() + { + throw new NotSupportedException(); + } + + /// + /// Normalizes this vector to a unit vector with respect to the p-norm. + /// + /// + /// The p value. + /// + /// + /// This vector normalized to a unit vector with respect to the p-norm. + /// + public override Vector Normalize(double p) + { + if (p < 0.0) + { + throw new ArgumentOutOfRangeException("p"); + } + + var norm = Norm(p); + var clone = Clone(); + if (norm.Real == 0.0f) + { + return clone; + } + + clone.Multiply(1.0f / norm, clone); + + return clone; + } + + /// + /// 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-positive. + public override Vector Random(int length, IContinuousDistribution randomDistribution) + { + if (length < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + var vector = CreateVector(length); + for (var index = 0; index < length; index++) + { + vector[index] = new Complex32(Convert.ToSingle(randomDistribution.Sample()), Convert.ToSingle(randomDistribution.Sample())); + } + + return vector; + } + + /// + /// 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 not positive. + public override Vector Random(int length, IDiscreteDistribution randomDistribution) + { + if (length < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + var vector = CreateVector(length); + for (var index = 0; index < length; index++) + { + vector[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 98f49638..f561b356 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseVector.cs @@ -39,7 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// /// A vector using dense storage. /// - public class DenseVector : Vector + public class DenseVector : Vector { /// /// Initializes a new instance of the class with a given size. @@ -961,13 +961,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The sum of the vector's elements. public override double Sum() { - double result = 0; - for (var i = 0; i < Count; i++) - { - result += Data[i]; - } - - return result; + return CommonParallel.Aggregate( + 0, + Count, + i => Data[i]); } /// @@ -976,13 +973,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The sum of the absolute value of the vector's elements. public override double SumMagnitudes() { - double result = 0; - for (var i = 0; i < Count; i++) - { - result += Math.Abs(Data[i]); - } - - return result; + return CommonParallel.Aggregate( + 0, + Count, + i => Math.Abs(Data[i])); } /// @@ -1193,58 +1187,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double 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 n vector is non positive. - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (DenseVector)CreateVector(length); - for (var index = 0; index < v.Data.Length; index++) - { - v.Data[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 positive. - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (DenseVector)CreateVector(length); - for (var index = 0; index < v.Data.Length; index++) - { - v.Data[index] = randomDistribution.Sample(); - } - - return v; - } - /// /// Outer product of this and another vector. /// @@ -1301,34 +1243,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double return Math.Pow(sum, 1.0 / p); } - - /// - /// Normalizes this vector to a unit vector with respect to the p-norm. - /// - /// - /// The p value. - /// - /// - /// This vector normalized to a unit vector with respect to the p-norm. - /// - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(1.0 / norm, clone); - - return clone; - } + #endregion #region Parse Functions @@ -1497,61 +1412,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double { Array.Clear(Data, 0, Data.Length); } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override double AddT(double val1, double val2) - { - return val1 + val2; - } - - /// - /// Subtract two values T-T - /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected sealed override double SubtractT(double val1, double val2) - { - return val1 - val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Divide two values T/T - /// - /// Left operand value - /// Right operand value - /// Result of divide - protected sealed override double DivideT(double val1, double val2) - { - return val1 / val2; - } - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/SparseVector.cs b/src/Numerics/LinearAlgebra/Double/SparseVector.cs index c006e501..f5a44801 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseVector.cs @@ -30,7 +30,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double using System.Collections.Generic; using System.Globalization; using System.Linq; - using Distributions; using Generic; using NumberTheory; using Properties; @@ -39,7 +38,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// /// A vector with sparse storage. /// - public class SparseVector : Vector + public class SparseVector : Vector { /// /// Lock object for the indexer. @@ -530,11 +529,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double var sparseother = other as SparseVector; if (sparseother == null) { - sparse.AddScaledSparseVector(1.0, sparseother); + base.Add(other, result); } else { - base.Add(other, result); + CopyTo(result); + sparse.AddScaledSparseVector(1.0, sparseother); } } } @@ -702,11 +702,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double var sparseother = other as SparseVector; if (sparseother == null) { - sparse.AddScaledSparseVector(-1.0, sparseother); + base.Subtract(other, result); } else { - base.Subtract(other, result); + CopyTo(result); + sparse.AddScaledSparseVector(-1.0, sparseother); } } } @@ -1152,58 +1153,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double 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 positive. - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (SparseVector)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 positive. - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (SparseVector)CreateVector(length); - for (var index = 0; index < v.Count; index++) - { - v[index] = randomDistribution.Sample(); - } - - return v; - } - /// /// Outer product of this and another vector. /// @@ -1256,33 +1205,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double return Math.Pow(sum, 1.0 / p); } - /// - /// Normalizes this vector to a unit vector with respect to the p-norm. - /// - /// - /// The p value. - /// - /// - /// This vector normalized to a unit vector with respect to the p-norm. - /// - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(1.0 / norm, clone); - - return clone; - } #endregion #region Parse Functions @@ -1623,61 +1545,5 @@ namespace MathNet.Numerics.LinearAlgebra.Double return true; } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override double AddT(double val1, double val2) - { - return val1 + val2; - } - - /// - /// Subtract two values T-T - /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected sealed override double SubtractT(double val1, double val2) - { - return val1 - val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - /// - /// Divide two values T/T - /// - /// Left operand value - /// Right operand value - /// Result of divide - protected sealed override double DivideT(double val1, double val2) - { - return val1 / val2; - } - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Double/Vector.cs b/src/Numerics/LinearAlgebra/Double/Vector.cs new file mode 100644 index 00000000..e45abccf --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Vector.cs @@ -0,0 +1,470 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// Copyright (c) 2009-2010 Math.NET +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.LinearAlgebra.Double +{ + using System; + using Distributions; + using Generic; + using Properties; + using Threading; + + /// + /// double version of the class. + /// + public abstract class Vector : Vector + { + /// + /// Initializes a new instance of the Vector class. + /// Constructs a Vector with the given size. + /// + /// + /// The size of the Vector to construct. + /// + /// + /// If is less than one. + /// + protected Vector(int size) : base(size) + { + } + + /// + /// 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. + /// + protected override void DoAdd(double scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] + scalar); + } + + /// + /// 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. + /// + protected override void DoAdd(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] + other[index]); + } + + /// + /// 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. + /// + protected override void DoSubtract(double scalar, Vector result) + { + 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. + /// + protected override void DoSubtract(Vector other, Vector result) + { + CopyTo(result); + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] - other[index]); + } + + /// + /// Multiplies a scalar to each element of the vector and stores the result in the result vector. + /// + /// + /// The scalar to multiply. + /// + /// + /// The vector to store the result of the multiplication. + /// + protected override void DoMultiply(double scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] * scalar); + } + + /// + /// Divides each element of the vector by a scalar and stores the result in the result vector. + /// + /// + /// The scalar to divide with. + /// + /// + /// The vector to store the result of the division. + /// + protected override void DoDivide(double scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] / scalar); + } + + /// + /// 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 DoPointwiseMultiply(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] * other[index]); + } + + /// + /// 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 override void DoPointwiseDivide(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] / other[index]); + } + + /// + /// 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) + { + return CommonParallel.Aggregate( + 0, + Count, + i => this[i] * other[i]); + } + + /// + /// Returns the value of the absolute minimum element. + /// + /// The value of the absolute minimum element. + public override double AbsoluteMinimum() + { + return Math.Abs(this[AbsoluteMinimumIndex()]); + } + + /// + /// Returns the index of the absolute minimum element. + /// + /// The index of absolute minimum element. + public override int AbsoluteMinimumIndex() + { + var index = 0; + var min = Math.Abs(this[index]); + for (var i = 1; i < Count; i++) + { + var test = Math.Abs(this[i]); + if (test < min) + { + index = i; + min = test; + } + } + + return index; + } + + /// + /// Returns the value of the absolute maximum element. + /// + /// The value of the absolute maximum element. + public override double AbsoluteMaximum() + { + return Math.Abs(this[AbsoluteMaximumIndex()]); + } + + /// + /// Returns the index of the absolute maximum element. + /// + /// The index of absolute maximum element. + public override int AbsoluteMaximumIndex() + { + var index = 0; + var max = Math.Abs(this[index]); + for (var i = 1; i < Count; i++) + { + var test = Math.Abs(this[i]); + if (test > max) + { + index = i; + max = test; + } + } + + return index; + } + + /// + /// Computes the sum of the vector's elements. + /// + /// The sum of the vector's elements. + public override double Sum() + { + return CommonParallel.Aggregate( + 0, + Count, + i => this[i]); + } + + /// + /// 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() + { + return CommonParallel.Aggregate( + 0, + Count, + i => Math.Abs(this[i])); + } + + /// + /// Computes the p-Norm. + /// + /// + /// The p value. + /// + /// + /// Scalar ret = (sum(abs(this[i])^p))^(1/p) + /// + public override double Norm(double p) + { + if (p < 0.0) + { + throw new ArgumentOutOfRangeException("p"); + } + + if (Double.IsPositiveInfinity(p)) + { + return CommonParallel.Select( + 0, + Count, + (index, localData) => Math.Max(localData, Math.Abs(this[index])), + Math.Max); + } + + var sum = CommonParallel.Aggregate( + 0, + Count, + index => Math.Pow(Math.Abs(this[index]), p)); + + return Math.Pow(sum, 1.0 / p); + } + + /// + /// Conjugates vector and save result to + /// + /// Target vector + protected override void DoConjugate(Vector target) + { + if (ReferenceEquals(this, target)) + { + return; + } + + CopyTo(target); + } + + /// + /// Returns a negated vector. + /// + /// + /// The negated vector. + /// + /// + /// Added as an alternative to the unary negation operator. + /// + public override Vector Negate() + { + var result = CreateVector(Count); + CommonParallel.For( + 0, + Count, + index => result[index] = -this[index]); + + return result; + } + + /// + /// Returns the index of the absolute maximum element. + /// + /// The index of absolute maximum element. + public override int MaximumIndex() + { + var index = 0; + var max = this[index]; + for (var i = 1; i < Count; i++) + { + var test = this[i]; + if (test > max) + { + index = i; + max = test; + } + } + + return index; + } + + /// + /// Returns the index of the minimum element. + /// + /// The index of minimum element. + public override int MinimumIndex() + { + var index = 0; + var min = this[index]; + for (var i = 1; i < Count; i++) + { + var test = this[i]; + if (test < min) + { + index = i; + min = test; + } + } + + return index; + } + + /// + /// Normalizes this vector to a unit vector with respect to the p-norm. + /// + /// + /// The p value. + /// + /// + /// This vector normalized to a unit vector with respect to the p-norm. + /// + public override Vector Normalize(double p) + { + if (p < 0.0) + { + throw new ArgumentOutOfRangeException("p"); + } + + var norm = Norm(p); + var clone = Clone(); + if (norm == 0.0) + { + return clone; + } + + clone.Multiply(1.0 / norm, clone); + + return clone; + } + + /// + /// 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 positive. + public override Vector Random(int length, IContinuousDistribution randomDistribution) + { + if (length < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + var v = CreateVector(length); + for (var index = 0; index < 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 positive. + public override Vector Random(int length, IDiscreteDistribution randomDistribution) + { + if (length < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + var v = CreateVector(length); + for (var index = 0; index < Count; index++) + { + this[index] = randomDistribution.Sample(); + } + + return v; + } + } +} diff --git a/src/Numerics/LinearAlgebra/Generic/Common.cs b/src/Numerics/LinearAlgebra/Generic/Common.cs new file mode 100644 index 00000000..ae443530 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Generic/Common.cs @@ -0,0 +1,21 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.LinearAlgebra.Generic +{ + internal static class Common + { + /// + /// Returns the maximum value. + /// + /// The first value. + /// The second value. + /// The maximum value. + public static float Max(float a, float b) + { + return Math.Max(a, b); + } + } +} diff --git a/src/Numerics/LinearAlgebra/Generic/Matrix.cs b/src/Numerics/LinearAlgebra/Generic/Matrix.cs index 2583e0d3..0bd87fe3 100644 --- a/src/Numerics/LinearAlgebra/Generic/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Generic/Matrix.cs @@ -1888,15 +1888,15 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// /// Take absolute value /// - /// Source alue - /// True if one; otherwise false + /// Source value + /// True if one; otherwise false protected abstract double AbsoluteT(T val1); /// /// Is equal to one? /// /// Value to check - /// True if one; otherwise false + /// True if one; otherwise false private static bool IsOneT(T val1) { if (typeof(T) == typeof(Complex)) diff --git a/src/Numerics/LinearAlgebra/Generic/Vector.cs b/src/Numerics/LinearAlgebra/Generic/Vector.cs index c664b7b8..ef8fb22b 100644 --- a/src/Numerics/LinearAlgebra/Generic/Vector.cs +++ b/src/Numerics/LinearAlgebra/Generic/Vector.cs @@ -49,6 +49,16 @@ namespace MathNet.Numerics.LinearAlgebra.Generic #endif where T : struct, IEquatable, IFormattable { + /// + /// The zero value for type T. + /// + private static readonly T Zero = default(T); + + /// + /// The value on 1.0 for type T. + /// + private static readonly T One = SetOne(); + /// /// Initializes a new instance of the Vector class. /// Constructs a Vector with the given size. @@ -127,14 +137,14 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// A copy of the vector with the scalar added. public virtual Vector Add(T scalar) { - if (scalar.Equals(default(T))) + if (scalar.Equals(Zero)) { return Clone(); } - var copy = Clone(); - Add(scalar, copy); - return copy; + var result = CreateVector(Count); + Add(scalar, result); + return result; } /// @@ -169,12 +179,20 @@ namespace MathNet.Numerics.LinearAlgebra.Generic CopyTo(result); } - CommonParallel.For( - 0, - Count, - index => result[index] = AddT(result[index], scalar)); + DoAdd(scalar, result); } + /// + /// 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. + /// + protected abstract void DoAdd(T scalar, Vector result); + /// /// Returns a copy of this vector. /// @@ -214,9 +232,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); } - var copy = Clone(); - Add(other, copy); - return copy; + var result = CreateVector(Count); + Add(other, result); + return result; } /// @@ -259,13 +277,21 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } else { - CommonParallel.For( - 0, - Count, - index => result[index] = AddT(this[index], other[index])); + DoAdd(other, result); } } + /// + /// 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. + /// + protected abstract void DoAdd(Vector other, Vector result); + /// /// Subtracts a scalar from each element of the vector. /// @@ -280,9 +306,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic return Clone(); } - var copy = Clone(); - Subtract(scalar, copy); - return copy; + var result = CreateVector(Count); + Subtract(scalar, result); + return result; } /// @@ -317,12 +343,20 @@ namespace MathNet.Numerics.LinearAlgebra.Generic CopyTo(result); } - CommonParallel.For( - 0, - Count, - index => result[index] = SubtractT(result[index], scalar)); + DoSubtract(scalar, result); } + /// + /// 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. + /// + protected abstract void DoSubtract(T scalar, Vector result); + /// /// Returns a negated vector. /// @@ -359,9 +393,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); } - var copy = Clone(); - Subtract(other, copy); - return copy; + var result = CreateVector(Count); + Subtract(other, result); + return result; } /// @@ -404,14 +438,21 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } else { - CopyTo(result); - CommonParallel.For( - 0, - Count, - index => result[index] = SubtractT(this[index], other[index])); + DoSubtract(other, 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. + /// + protected abstract void DoSubtract(Vector other, Vector result); + /// /// Multiplies a scalar to each element of the vector. /// @@ -421,14 +462,14 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// A new vector that is the multiplication of the vector and the scalar. public virtual Vector Multiply(T scalar) { - if (IsOneT(scalar)) + if (scalar.Equals(One)) { return Clone(); } - var copy = Clone(); - Multiply(scalar, copy); - return copy; + var result = CreateVector(Count); + Multiply(scalar, result); + return result; } /// @@ -463,19 +504,27 @@ namespace MathNet.Numerics.LinearAlgebra.Generic CopyTo(result); } - CommonParallel.For( - 0, - Count, - index => result[index] = MultiplyT(result[index], scalar)); + DoMultiply(scalar, result); } + /// + /// Multiplies a scalar to each element of the vector and stores the result in the result vector. + /// + /// + /// The scalar to multiply. + /// + /// + /// The vector to store the result of the multiplication. + /// + protected abstract void DoMultiply(T scalar, Vector result); + /// /// Computes the dot product between this vector and another vector. /// /// /// The other vector to add. /// - /// + /// s /// The result of the addition. /// /// @@ -496,15 +545,20 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); } - var dot = default(T); - for (var i = 0; i < Count; i++) - { - dot = AddT(dot, MultiplyT(this[i], other[i])); - } - - return dot; + return DoDotProduct(other); } + /// + /// Computes the dot product between this vector and another vector. + /// + /// + /// The other vector to add. + /// + /// s + /// The result of the addition. + /// + protected abstract T DoDotProduct(Vector other); + /// /// Divides each element of the vector by a scalar. /// @@ -514,14 +568,14 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// A new vector that is the division of the vector and the scalar. public virtual Vector Divide(T scalar) { - if (IsOneT(scalar)) + if (scalar.Equals(One)) { return Clone(); } - var copy = Clone(); - Divide(scalar, copy); - return copy; + var result = CreateVector(Count); + Divide(scalar, result); + return result; } /// @@ -556,12 +610,20 @@ namespace MathNet.Numerics.LinearAlgebra.Generic CopyTo(result); } - CommonParallel.For( - 0, - Count, - index => result[index] = DivideT(result[index], scalar)); + DoDivide(scalar, result); } + /// + /// Divides each element of the vector by a scalar and stores the result in the result vector. + /// + /// + /// The scalar to divide with. + /// + /// + /// The vector to store the result of the division. + /// + protected abstract void DoDivide(T scalar, Vector result); + /// /// Pointwise multiplies this vector with another vector. /// @@ -581,9 +643,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); } - var copy = Clone(); - PointwiseMultiply(other, copy); - return copy; + var result = CreateVector(Count); + PointwiseMultiply(other, result); + return result; } /// @@ -624,13 +686,17 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } else { - CommonParallel.For( - 0, - Count, - index => result[index] = MultiplyT(this[index], other[index])); + DoPointwiseMultiply(other, result); } } + /// + /// 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 abstract void DoPointwiseMultiply(Vector other, Vector result); + /// /// Pointwise divide this vector with another vector. /// @@ -650,9 +716,9 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); } - var copy = Clone(); - PointwiseDivide(other, copy); - return copy; + var result = CreateVector(Count); + PointwiseDivide(other, result); + return result; } /// @@ -693,13 +759,17 @@ namespace MathNet.Numerics.LinearAlgebra.Generic } else { - CommonParallel.For( - 0, - Count, - index => result[index] = DivideT(this[index], other[index])); + DoPointwiseDivide(other, result); } } + /// + /// 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); + /// /// Outer product of two vectors /// @@ -769,61 +839,25 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// Returns the value of the absolute minimum element. /// /// The value of the absolute minimum element. - public virtual double AbsoluteMinimum() - { - return AbsoluteT(this[AbsoluteMinimumIndex()]); - } + public abstract T AbsoluteMinimum(); /// /// Returns the index of the absolute minimum element. /// /// The index of absolute minimum element. - public virtual int AbsoluteMinimumIndex() - { - var index = 0; - var min = AbsoluteT(this[index]); - for (var i = 1; i < Count; i++) - { - var test = AbsoluteT(this[i]); - if (test < min) - { - index = i; - min = test; - } - } + public abstract int AbsoluteMinimumIndex(); - return index; - } - /// /// Returns the value of the absolute maximum element. /// /// The value of the absolute maximum element. - public virtual double AbsoluteMaximum() - { - return AbsoluteT(this[AbsoluteMaximumIndex()]); - } + public abstract T AbsoluteMaximum(); /// /// Returns the index of the absolute maximum element. /// /// The index of absolute maximum element. - public virtual int AbsoluteMaximumIndex() - { - var index = 0; - var max = AbsoluteT(this[index]); - for (var i = 1; i < Count; i++) - { - var test = AbsoluteT(this[i]); - if (test > max) - { - index = i; - max = test; - } - } - - return index; - } + public abstract int AbsoluteMaximumIndex(); /// /// Returns the value of maximum element. @@ -859,31 +893,13 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// Computes the sum of the vector's elements. /// /// The sum of the vector's elements. - public virtual T Sum() - { - var result = default(T); - for (var i = 0; i < Count; i++) - { - result = AddT(result, this[i]); - } - - return result; - } + public abstract T Sum(); /// /// Computes the sum of the absolute value of the vector's elements. /// /// The sum of the absolute value of the vector's elements. - public virtual double SumMagnitudes() - { - double result = 0; - for (var i = 0; i < Count; i++) - { - result += AbsoluteT(this[i]); - } - - return result; - } + public abstract T SumMagnitudes(); #endregion @@ -1070,29 +1086,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// /// Scalar ret = (sum(abs(this[i])^p))^(1/p) /// - public virtual double Norm(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - if (Double.IsPositiveInfinity(p)) - { - return CommonParallel.Select( - 0, - Count, - (index, localData) => Math.Max(localData, AbsoluteT(this[index])), - Math.Max); - } - - var sum = CommonParallel.Aggregate( - 0, - Count, - index => Math.Pow(AbsoluteT(this[index]), p)); - - return Math.Pow(sum, 1.0 / p); - } + public abstract T Norm(double p); /// /// Normalizes this vector to a unit vector with respect to the p-norm. @@ -1124,13 +1118,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// Target vector public virtual void Conjugate(Vector target) { - // In case of real return copy of vector - if (typeof(T) == typeof(double) || (typeof(T) == typeof(float))) - { - CopyTo(target); - return; - } - if (target == null) { throw new ArgumentNullException("target"); @@ -1141,19 +1128,15 @@ namespace MathNet.Numerics.LinearAlgebra.Generic throw new ArgumentException(Resources.ArgumentVectorsSameLength, "target"); } - if (ReferenceEquals(this, target)) - { - var tmp = CreateVector(Count); - Conjugate(tmp); - tmp.CopyTo(target); - } - - CommonParallel.For( - 0, - Count, - index => target[index] = ConjugateT(this[index])); + DoConjugate(target); } + /// + /// Conjugates vector and save result to + /// + /// Target vector + protected abstract void DoConjugate(Vector target); + #region Copying and Conversion /// @@ -1604,114 +1587,33 @@ namespace MathNet.Numerics.LinearAlgebra.Generic CommonParallel.For(0, Count, index => this[index] = default(T)); } - #region Simple arithmetic of type T - - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected abstract T AddT(T val1, T val2); - - /// - /// Subtract two values T-T - /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected abstract T SubtractT(T val1, T val2); - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected abstract T MultiplyT(T val1, T val2); - - /// - /// Divide two values T/T - /// - /// Left operand value - /// Right operand value - /// Result of divide - protected abstract T DivideT(T val1, T val2); - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected abstract double AbsoluteT(T val1); - - /// - /// Is equal to one? - /// - /// Value to check - /// True if one; otherwise false - private static bool IsOneT(T val1) - { - if (typeof(T) == typeof(Complex)) - { - object obj1 = val1; - return Complex.One.AlmostEqual((Complex)obj1); - } - - if (typeof(T) == typeof(Complex32)) - { - object obj1 = val1; - return Complex32.One.AlmostEqual((Complex32)obj1); - } - - if (typeof(T) == typeof(double)) - { - object obj1 = val1; - return 1.0.AlmostEqualInDecimalPlaces((double)obj1, 15); - } - - if (typeof(T) == typeof(float)) - { - object obj1 = val1; - return 1.0f.AlmostEqualInDecimalPlaces((float)obj1, 7); - } - - throw new NotSupportedException(); - } - /// - /// Conjugate complex value. In real case the same value is returned + /// Sets the value of 1.0 for type T. /// - /// Value to conjugate - /// Conjugated value (complex) or the same (real) - private static T ConjugateT(T val1) + /// The value of 1.0 for type T. + private static T SetOne() { if (typeof(T) == typeof(Complex)) { - object obj = val1; - object conj = Complex.Conjugate((Complex)obj); - return (T)conj; + return (T)(object)Complex.One; } if (typeof(T) == typeof(Complex32)) { - object obj = val1; - object conj = ((Complex32)obj).Conjugate(); - return (T)conj; + return (T)(object)Complex32.One; } if (typeof(T) == typeof(double)) { - return val1; + return (T)(object)1.0; } if (typeof(T) == typeof(float)) { - return val1; + return (T)(object)1.0f; } throw new NotSupportedException(); } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/DenseVector.cs b/src/Numerics/LinearAlgebra/Single/DenseVector.cs index c43bf0b8..eec92b19 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseVector.cs @@ -39,7 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// /// A vector using dense storage. /// - public class DenseVector : Vector + public class DenseVector : Vector { /// /// Initializes a new instance of the class with a given size. @@ -819,7 +819,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Returns the value of the absolute minimum element. /// /// The value of the absolute minimum element. - public override double AbsoluteMinimum() + public override float AbsoluteMinimum() { return Math.Abs(Data[AbsoluteMinimumIndex()]); } @@ -828,7 +828,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Returns the value of the absolute maximum element. /// /// The value of the absolute maximum element. - public override double AbsoluteMaximum() + public override float AbsoluteMaximum() { return Math.Abs(Data[AbsoluteMaximumIndex()]); } @@ -961,28 +961,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The sum of the vector's elements. public override float Sum() { - float result = 0; - for (var i = 0; i < Count; i++) - { - result += Data[i]; - } - - return result; + return CommonParallel.Aggregate( + 0, + Count, + i => Data[i]); } /// /// 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() + public override float SumMagnitudes() { - float result = 0; - for (var i = 0; i < Count; i++) - { - result += Math.Abs(Data[i]); - } - - return result; + return CommonParallel.Aggregate( + 0, + Count, + i => Math.Abs(Data[i])); } /// @@ -1193,58 +1187,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single 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 n vector is non positive. - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (DenseVector)CreateVector(length); - for (var index = 0; index < v.Data.Length; index++) - { - v.Data[index] = (float)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 positive. - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (DenseVector)CreateVector(length); - for (var index = 0; index < v.Data.Length; index++) - { - v.Data[index] = randomDistribution.Sample(); - } - - return v; - } - /// /// Outer product of this and another vector. /// @@ -1265,7 +1207,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// /// The p value. /// Scalar ret = (sum(abs(this[i])^p))^(1/p) - public override double Norm(double p) + public override float Norm(double p) { if (p < 0.0) { @@ -1291,7 +1233,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single 0, Count, (index, localData) => Math.Max(localData, Math.Abs(Data[index])), - Math.Max); + Common.Max); } var sum = CommonParallel.Aggregate( @@ -1299,36 +1241,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single Count, index => Math.Pow(Math.Abs(Data[index]), p)); - return Math.Pow(sum, 1.0 / p); + return (float)Math.Pow(sum, 1.0 / p); } - /// - /// Normalizes this vector to a unit vector with respect to the p-norm. - /// - /// - /// The p value. - /// - /// - /// This vector normalized to a unit vector with respect to the p-norm. - /// - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(1.0f / (float)norm, clone); - - return clone; - } #endregion #region Parse Functions @@ -1497,61 +1412,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single { Array.Clear(Data, 0, Data.Length); } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override float AddT(float val1, float val2) - { - return val1 + val2; - } - - /// - /// Subtract two values T-T - /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected sealed override float SubtractT(float val1, float val2) - { - return val1 - val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - /// - /// Divide two values T/T - /// - /// Left operand value - /// Right operand value - /// Result of divide - protected sealed override float DivideT(float val1, float val2) - { - return val1 / val2; - } - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected sealed override double AbsoluteT(float val1) - { - return Math.Abs(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/StopCriterium/ResidualStopCriterium.cs b/src/Numerics/LinearAlgebra/Single/Solvers/StopCriterium/ResidualStopCriterium.cs index 4f828c6b..25f0790e 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/StopCriterium/ResidualStopCriterium.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/StopCriterium/ResidualStopCriterium.cs @@ -45,7 +45,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium /// /// The default value for the maximum value of the residual. /// - public const double DefaultMaximumResidual = 1e-6; + public const float DefaultMaximumResidual = 1e-6f; /// /// The default value for the minimum number of iterations. @@ -65,7 +65,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium /// /// The maximum value for the residual below which the calculation is considered converged. /// - private double _maximum; + private float _maximum; /// /// The minimum number of iterations for which the residual has to be below the maximum before @@ -101,7 +101,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium /// maximum residual and the default minimum number of iterations. /// /// The maximum value for the residual below which the calculation is considered converged. - public ResidualStopCriterium(double maximum) : this(maximum, DefaultMinimumIterationsBelowMaximum) + public ResidualStopCriterium(float maximum) : this(maximum, DefaultMinimumIterationsBelowMaximum) { } @@ -128,7 +128,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium /// The minimum number of iterations for which the residual has to be below the maximum before /// the calculation is considered converged. /// - public ResidualStopCriterium(double maximum, int minimumIterationsBelowMaximum) + public ResidualStopCriterium(float maximum, int minimumIterationsBelowMaximum) { if (maximum < 0) { @@ -149,7 +149,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium /// converged. /// /// Thrown if the Maximum is set to a negative value. - public double Maximum + public float Maximum { [DebuggerStepThrough] get @@ -258,16 +258,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium // Store the infinity norms of both the solution and residual vectors // These values will be used to calculate the relative drop in residuals // later on. - var residualNorm = residualVector.Norm(Double.PositiveInfinity); + var residualNorm = residualVector.Norm(float.PositiveInfinity); // Check the residuals by calculating: // ||r_i|| <= stop_tol * ||b|| - var stopCriterium = ComputeStopCriterium(sourceVector.Norm(Double.PositiveInfinity)); + var stopCriterium = ComputeStopCriterium(sourceVector.Norm(float.PositiveInfinity)); // First check that we have real numbers not NaN's. // NaN's can occur when the iterative process diverges so we // stop if that is the case. - if (double.IsNaN(stopCriterium) || double.IsNaN(residualNorm)) + if (float.IsNaN(stopCriterium) || float.IsNaN(residualNorm)) { _iterationCount = 0; SetStatusToDiverged(); @@ -306,7 +306,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium /// /// Solution vector norm /// Criterium value - private double ComputeStopCriterium(double solutionNorm) + private float ComputeStopCriterium(float solutionNorm) { // This is criterium 1 from Templates for the solution of linear systems. // The problem with this criterium is that it's not limiting enough. For now diff --git a/src/Numerics/LinearAlgebra/Single/SparseVector.cs b/src/Numerics/LinearAlgebra/Single/SparseVector.cs index 46937e42..8b0942d3 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseVector.cs @@ -30,7 +30,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single using System.Collections.Generic; using System.Globalization; using System.Linq; - using Distributions; using Generic; using NumberTheory; using Properties; @@ -39,7 +38,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// /// A vector with sparse storage. /// - public class SparseVector : Vector + public class SparseVector : Vector { /// /// Lock object for the indexer. @@ -530,11 +529,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single var sparseother = other as SparseVector; if (sparseother == null) { - sparse.AddScaledSparseVector(1.0f, sparseother); + base.Add(other, result); } else { - base.Add(other, result); + CopyTo(result); + sparse.AddScaledSparseVector(1.0f, sparseother); } } } @@ -702,11 +702,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single var sparseother = other as SparseVector; if (sparseother == null) { - sparse.AddScaledSparseVector(-1.0f, sparseother); + base.Subtract(other, result); } else { - base.Subtract(other, result); + CopyTo(result); + sparse.AddScaledSparseVector(-1.0f, sparseother); } } } @@ -1061,7 +1062,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The sum of the vector's elements. public override float Sum() { - float result = 0; + var result = 0.0f; for (var i = 0; i < NonZerosCount; i++) { result += _nonZeroValues[i]; @@ -1074,9 +1075,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// 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() + public override float SumMagnitudes() { - double result = 0; + var result = 0.0f; for (var i = 0; i < NonZerosCount; i++) { result += Math.Abs(_nonZeroValues[i]); @@ -1152,58 +1153,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single 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 positive. - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (SparseVector)CreateVector(length); - for (var index = 0; index < v.Count; index++) - { - v[index] = (float)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 positive. - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (SparseVector)CreateVector(length); - for (var index = 0; index < v.Count; index++) - { - v[index] = randomDistribution.Sample(); - } - - return v; - } - /// /// Outer product of this and another vector. /// @@ -1226,7 +1175,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// /// The p value. /// Scalar ret = (sum(abs(this[i])^p))^(1/p) - public override double Norm(double p) + public override float Norm(double p) { if (1 > p) { @@ -1235,7 +1184,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single if (NonZerosCount == 0) { - return 0.0; + return 0.0f; } if (2.0 == p) @@ -1245,7 +1194,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single if (Double.IsPositiveInfinity(p)) { - return CommonParallel.Select(0, NonZerosCount, (index, localData) => Math.Max(localData, Math.Abs(_nonZeroValues[index])), Math.Max); + return CommonParallel.Select(0, NonZerosCount, (index, localData) => Math.Max(localData, Math.Abs(_nonZeroValues[index])), Common.Max); } var sum = CommonParallel.Aggregate( @@ -1253,36 +1202,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single NonZerosCount, index => Math.Pow(Math.Abs(_nonZeroValues[index]), p)); - return Math.Pow(sum, 1.0 / p); + return (float)Math.Pow(sum, 1.0 / p); } - /// - /// Normalizes this vector to a unit vector with respect to the p-norm. - /// - /// - /// The p value. - /// - /// - /// This vector normalized to a unit vector with respect to the p-norm. - /// - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(1.0f / (float)norm, clone); - - return clone; - } #endregion #region Parse Functions @@ -1623,61 +1545,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single return true; } - - #region Simple arithmetic of type T - /// - /// Add two values T+T - /// - /// Left operand value - /// Right operand value - /// Result of addition - protected sealed override float AddT(float val1, float val2) - { - return val1 + val2; - } - - /// - /// Subtract two values T-T - /// - /// Left operand value - /// Right operand value - /// Result of subtract - protected sealed override float SubtractT(float val1, float val2) - { - return val1 - val2; - } - - /// - /// Multiply two values T*T - /// - /// Left operand value - /// Right operand value - /// Result of multiplication - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - /// - /// Divide two values T/T - /// - /// Left operand value - /// Right operand value - /// Result of divide - protected sealed override float DivideT(float val1, float val2) - { - return val1 / val2; - } - - /// - /// Take absolute value - /// - /// Source alue - /// True if one; otherwise false - protected sealed override double AbsoluteT(float val1) - { - return Math.Abs(val1); - } - #endregion } } diff --git a/src/Numerics/LinearAlgebra/Single/Vector.cs b/src/Numerics/LinearAlgebra/Single/Vector.cs new file mode 100644 index 00000000..56b1201d --- /dev/null +++ b/src/Numerics/LinearAlgebra/Single/Vector.cs @@ -0,0 +1,473 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// Copyright (c) 2009-2010 Math.NET +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.LinearAlgebra.Single +{ + using System; + using Distributions; + using Generic; + using Properties; + using Threading; + + /// + /// float version of the class. + /// + public abstract class Vector : Vector + { + /// + /// Initializes a new instance of the Vector class. + /// Constructs a Vector with the given size. + /// + /// + /// The size of the Vector to construct. + /// + /// + /// If is less than one. + /// + protected Vector(int size) + : base(size) + { + } + + /// + /// 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. + /// + protected override void DoAdd(float scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] + scalar); + } + + /// + /// 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. + /// + protected override void DoAdd(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] + other[index]); + } + + /// + /// 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. + /// + protected override void DoSubtract(float scalar, Vector result) + { + 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. + /// + protected override void DoSubtract(Vector other, Vector result) + { + CopyTo(result); + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] - other[index]); + } + + /// + /// Multiplies a scalar to each element of the vector and stores the result in the result vector. + /// + /// + /// The scalar to multiply. + /// + /// + /// The vector to store the result of the multiplication. + /// + protected override void DoMultiply(float scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] * scalar); + } + + /// + /// Divides each element of the vector by a scalar and stores the result in the result vector. + /// + /// + /// The scalar to divide with. + /// + /// + /// The vector to store the result of the division. + /// + protected override void DoDivide(float scalar, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = result[index] / scalar); + } + + /// + /// 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 DoPointwiseMultiply(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] * other[index]); + } + + /// + /// 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 override void DoPointwiseDivide(Vector other, Vector result) + { + CommonParallel.For( + 0, + Count, + index => result[index] = this[index] / other[index]); + } + + /// + /// Computes the dot product between this vector and another vector. + /// + /// + /// The other vector to add. + /// + /// s + /// The result of the addition. + /// + protected override float DoDotProduct(Vector other) + { + return CommonParallel.Aggregate( + 0, + Count, + i => this[i] * other[i]); + } + + /// + /// Returns the value of the absolute minimum element. + /// + /// The value of the absolute minimum element. + public override float AbsoluteMinimum() + { + return Math.Abs(this[AbsoluteMinimumIndex()]); + } + + /// + /// Returns the index of the absolute minimum element. + /// + /// The index of absolute minimum element. + public override int AbsoluteMinimumIndex() + { + var index = 0; + var min = Math.Abs(this[index]); + for (var i = 1; i < Count; i++) + { + var test = Math.Abs(this[i]); + if (test < min) + { + index = i; + min = test; + } + } + + return index; + } + + /// + /// Returns the value of the absolute maximum element. + /// + /// The value of the absolute maximum element. + public override float AbsoluteMaximum() + { + return Math.Abs(this[AbsoluteMaximumIndex()]); + } + + /// + /// Returns the index of the absolute maximum element. + /// + /// The index of absolute maximum element. + public override int AbsoluteMaximumIndex() + { + var index = 0; + var max = Math.Abs(this[index]); + for (var i = 1; i < Count; i++) + { + var test = Math.Abs(this[i]); + if (test > max) + { + index = i; + max = test; + } + } + + return index; + } + + /// + /// Computes the sum of the vector's elements. + /// + /// The sum of the vector's elements. + public override float Sum() + { + return CommonParallel.Aggregate( + 0, + Count, + i => this[i]); + } + + /// + /// 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 float SumMagnitudes() + { + return CommonParallel.Aggregate( + 0, + Count, + i => Math.Abs(this[i])); + } + + /// + /// Computes the p-Norm. + /// + /// + /// The p value. + /// + /// + /// Scalar ret = (sum(abs(this[i])^p))^(1/p) + /// + public override float Norm(double p) + { + if (p < 0.0) + { + throw new ArgumentOutOfRangeException("p"); + } + + if (float.IsPositiveInfinity((float)p)) + { + return CommonParallel.Select( + 0, + Count, + (index, localData) => Math.Max(localData, Math.Abs(this[index])), + Common.Max); + } + + var sum = CommonParallel.Aggregate( + 0, + Count, + index => Math.Pow(Math.Abs(this[index]), p)); + + return (float)Math.Pow(sum, 1.0 / p); + } + + /// + /// Conjugates vector and save result to + /// + /// Target vector + protected override void DoConjugate(Vector target) + { + if (ReferenceEquals(this, target)) + { + return; + } + + CopyTo(target); + } + + /// + /// Returns a negated vector. + /// + /// + /// The negated vector. + /// + /// + /// Added as an alternative to the unary negation operator. + /// + public override Vector Negate() + { + var result = CreateVector(Count); + CommonParallel.For( + 0, + Count, + index => result[index] = -this[index]); + + return result; + } + + /// + /// Returns the index of the absolute maximum element. + /// + /// The index of absolute maximum element. + public override int MaximumIndex() + { + var index = 0; + var max = this[index]; + for (var i = 1; i < Count; i++) + { + var test = this[i]; + if (test > max) + { + index = i; + max = test; + } + } + + return index; + } + + /// + /// Returns the index of the minimum element. + /// + /// The index of minimum element. + public override int MinimumIndex() + { + var index = 0; + var min = this[index]; + for (var i = 1; i < Count; i++) + { + var test = this[i]; + if (test < min) + { + index = i; + min = test; + } + } + + return index; + } + + /// + /// Normalizes this vector to a unit vector with respect to the p-norm. + /// + /// + /// The p value. + /// + /// + /// This vector normalized to a unit vector with respect to the p-norm. + /// + public override Vector Normalize(double p) + { + if (p < 0.0) + { + throw new ArgumentOutOfRangeException("p"); + } + + var norm = Norm(p); + var clone = Clone(); + if (norm == 0.0) + { + return clone; + } + + clone.Multiply(1.0f / norm, clone); + + return clone; + } + + /// + /// 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 positive. + public override Vector Random(int length, IContinuousDistribution randomDistribution) + { + if (length < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + var v = CreateVector(length); + for (var index = 0; index < Count; index++) + { + v[index] = Convert.ToSingle(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 positive. + public override Vector Random(int length, IDiscreteDistribution randomDistribution) + { + if (length < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); + } + + var v = CreateVector(length); + for (var index = 0; index < Count; index++) + { + v[index] = Convert.ToSingle(randomDistribution.Sample()); + } + + return v; + } + + + } +} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index d76d3fb4..148d38b1 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -126,9 +126,13 @@ + + + + @@ -305,6 +309,7 @@ + diff --git a/src/Numerics/Precision.cs b/src/Numerics/Precision.cs index 53042a3d..16af911b 100644 --- a/src/Numerics/Precision.cs +++ b/src/Numerics/Precision.cs @@ -200,6 +200,42 @@ namespace MathNet.Numerics #endif } + + /// + /// Returns the magnitude of the number. + /// + /// The value. + /// The magnitude of the number. + public static int Magnitude(this float value) + { + // Can't do this with zero because the 10-log of zero doesn't exist. + if (value.Equals(0.0f)) + { + return 0; + } + + // Note that we need the absolute value of the input because Log10 doesn't + // work for negative numbers (obviously). + var magnitude = Convert.ToSingle(Math.Log10(Math.Abs(value))); + + // To get the right number we need to know if the value is negative or positive + // truncating a positive number will always give use the correct magnitude + // truncating a negative number will give us a magnitude that is off by 1 + if (magnitude < 0) + { +#if SILVERLIGHT + return (int)Truncate(magnitude - 1); +#else + return (int)Math.Truncate(magnitude - 1); +#endif + } + +#if SILVERLIGHT + return (int)Truncate(magnitude); +#else + return (int)Math.Truncate(magnitude); +#endif + } /// /// Returns the number divided by it's magnitude, effectively returning a number between -10 and 10. /// @@ -1582,6 +1618,34 @@ namespace MathNet.Numerics return CompareToInDecimalPlaces(a, b, decimalPlaces) < 0; } + + /// + /// Compares two floats and determines if the first value is smaller than the second + /// value to within the specified number of decimal places or not. + /// + /// + /// + /// The values are equal if the difference between the two numbers is smaller than 10^(-numberOfDecimalPlaces). We divide by + /// two so that we have half the range on each side of thg. if == 2, then 0.01 will equal between + /// 0.005 and 0.015, but not 0.02 and not 0.00 + /// + /// + /// The first value. + /// The second value. + /// The number of decimal places. + /// true if the first value is smaller than the second value; otherwise false. + public static bool IsSmallerWithDecimalPlaces(this float a, float b, int decimalPlaces) + { + // If A or B are a NAN, return false. NANs are equal to nothing, + // not even themselves, and thus they're not bigger or + // smaller than anything either + if (double.IsNaN(a) || double.IsNaN(b)) + { + return false; + } + + return CompareToInDecimalPlaces(a, b, decimalPlaces) < 0; + } /// /// Compares two doubles and determines which double is bigger. diff --git a/src/Numerics/Threading/CommonParallel.cs b/src/Numerics/Threading/CommonParallel.cs index 6bc6dafa..7f730354 100644 --- a/src/Numerics/Threading/CommonParallel.cs +++ b/src/Numerics/Threading/CommonParallel.cs @@ -385,5 +385,57 @@ namespace MathNet.Numerics.Threading return ret; } + + /// + /// Selects an item (such as Max or Min). + /// + /// Starting index of the loop. + /// Ending index of the loop + /// The function to select items over a subset. + /// The function to select the item of selection from the subsets. + /// The selected value. + public static float Select(int fromInclusive, int toExclusive, Func body, Func localFinally) + { + float ret = 0; + var syncLock = new object(); + +#if SILVERLIGHT + Parallel.For( + fromInclusive, + toExclusive, + () => 0.0f, + (i, localData) => localData += body(i, localData), + localResult => + { + lock (syncLock) + { + ret = localFinally(ret, localResult); + } + }); +#else + Parallel.ForEach( + Partitioner.Create(fromInclusive, toExclusive), + new ParallelOptions { MaxDegreeOfParallelism = Control.NumberOfParallelWorkerThreads }, + () => 0.0f, + (range, loop, localData) => + { + for (var i = range.Item1; i < range.Item2; i++) + { + localData = body(i, localData); + } + + return localData; + }, + localResult => + { + lock (syncLock) + { + ret = localFinally(ret, localResult); + } + }); +#endif + + return ret; + } } } \ No newline at end of file diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj index d4b446b9..64bc4a4f 100644 --- a/src/Silverlight/Silverlight.csproj +++ b/src/Silverlight/Silverlight.csproj @@ -368,6 +368,9 @@ LinearAlgebra\Complex32\SparseVector.cs + + LinearAlgebra\Complex32\Vector.cs + LinearAlgebra\Complex\DenseMatrix.cs @@ -464,6 +467,9 @@ LinearAlgebra\Complex\SparseVector.cs + + LinearAlgebra\Complex\Vector.cs + LinearAlgebra\Double\DenseMatrix.cs @@ -572,6 +578,12 @@ LinearAlgebra\Double\SparseVector.cs + + LinearAlgebra\Double\Vector.cs + + + LinearAlgebra\Generic\Common.cs + LinearAlgebra\Generic\Factorization\Evd.cs @@ -743,6 +755,9 @@ LinearAlgebra\Generic\Vector.cs + + LinearAlgebra\Single\Vector.cs + NumberTheory\IntegerTheory.cs diff --git a/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedVectorTests.cs index fa1633e5..b8ada682 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedVectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/UserDefinedVectorTests.cs @@ -34,10 +34,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex using System.Numerics; using System.Collections.Generic; using Distributions; + using LinearAlgebra.Complex; using LinearAlgebra.Generic; using Threading; - internal class UserDefinedVector : Vector + internal class UserDefinedVector : Vector { private readonly Complex[] _data; @@ -76,128 +77,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex return new UserDefinedVector(size); } - public override Vector Negate() - { - var result = new UserDefinedVector(Count); - CommonParallel.For( - 0, - _data.Length, - index => result[index] = -_data[index]); - - return result; - } - - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException("length"); - } - - var v = (UserDefinedVector)CreateVector(length); - for (var index = 0; index < v._data.Length; index++) - { - v._data[index] = new Complex(randomDistribution.Sample(), randomDistribution.Sample()); - } - - return v; - } - - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException("length"); - } - - var v = (UserDefinedVector)CreateVector(length); - for (var index = 0; index < v._data.Length; index++) - { - v._data[index] = new Complex(randomDistribution.Sample(), randomDistribution.Sample()); - } - - return v; - } - - public override int MinimumIndex() - { - throw new NotSupportedException(); - } - - public override int MaximumIndex() - { - throw new NotSupportedException(); - } - - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(1.0 / norm, clone); - - return clone; - } - - protected sealed override Complex AddT(Complex val1, Complex val2) - { - return val1 + val2; - } - - protected sealed override Complex SubtractT(Complex val1, Complex val2) - { - return val1 - val2; - } - - protected sealed override Complex MultiplyT(Complex val1, Complex val2) - { - return val1 * val2; - } - - protected sealed override Complex DivideT(Complex val1, Complex val2) - { - return val1 / val2; - } - - protected sealed override double AbsoluteT(Complex val1) - { - return val1.Magnitude; - } - - public override void Conjugate(Vector target) - { - if (target == null) - { - throw new ArgumentNullException("target"); - } - - if (Count != target.Count) - { - throw new ArgumentException("target"); - } - - if (ReferenceEquals(this, target)) - { - var tmp = CreateVector(Count); - Conjugate(tmp); - tmp.CopyTo(target); - } - - CommonParallel.For( - 0, - Count, - index => target[index] = this[index].Conjugate()); - } - } + } public class UserDefinedVectorTests : VectorTests { diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs index d81c5ee2..156b6d67 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs @@ -707,7 +707,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 for (var j = 0; j < result.ColumnCount; j++) { var col = result.Column(j); - AssertHelpers.AlmostEqual(Complex32.One, (float)col.Norm(pValue), 6); + AssertHelpers.AlmostEqual(Complex32.One, col.Norm(pValue).Real, 6); } } @@ -721,7 +721,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 for (var i = 0; i < matrix.RowCount; i++) { var row = matrix.Row(i); - AssertHelpers.AlmostEqual(Complex32.One, (float)row.Norm(pValue), 6); + AssertHelpers.AlmostEqual(Complex32.One, row.Norm(pValue).Real, 6); } } diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/StopCriterium/ResidualStopCriteriumTest.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/StopCriterium/ResidualStopCriteriumTest.cs index c5fe4d6d..13a5e5d9 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/StopCriterium/ResidualStopCriteriumTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/StopCriterium/ResidualStopCriteriumTest.cs @@ -13,7 +13,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [ExpectedArgumentOutOfRangeException] public void CreateWithNegativeMaximum() { - new ResidualStopCriterium(-0.1); + new ResidualStopCriterium(-0.1f); Assert.Fail(); } @@ -29,10 +29,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [MultipleAsserts] public void Create() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); - Assert.AreEqual(1e-8, criterium.Maximum, "Incorrect maximum"); + Assert.AreEqual(1e-6f, criterium.Maximum, "Incorrect maximum"); Assert.AreEqual(50, criterium.MinimumIterationsBelowMaximum, "Incorrect iteration count"); } @@ -40,7 +40,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [MultipleAsserts] public void ResetMaximum() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.ResetMaximumResidualToDefault(); @@ -51,7 +51,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [MultipleAsserts] public void ResetMinimumIterationsBelowMaximum() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.ResetMinimumIterationsBelowMaximumToDefault(); @@ -62,7 +62,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [ExpectedArgumentOutOfRangeException] public void DetermineStatusWithIllegalIterationNumber() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(-1, @@ -76,7 +76,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [ExpectedArgumentNullException] public void DetermineStatusWithNullSolutionVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -90,7 +90,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [ExpectedArgumentNullException] public void DetermineStatusWithNullSourceVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -104,7 +104,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [ExpectedArgumentNullException] public void DetermineStatusWithNullResidualVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -117,7 +117,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [ExpectedArgumentException] public void DetermineStatusWithNonMatchingSolutionVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -130,7 +130,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [ExpectedArgumentException] public void DetermineStatusWithNonMatchingSourceVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -143,7 +143,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [ExpectedArgumentException] public void DetermineStatusWithNonMatchingResidualVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -156,7 +156,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [MultipleAsserts] public void DetermineStatusWithSourceNaN() { - var criterium = new ResidualStopCriterium(1e-3, 10); + var criterium = new ResidualStopCriterium(1e-3f, 10); Assert.IsNotNull(criterium, "There should be a criterium"); var solution = new DenseVector(new[] { new Complex32(1.0f, 1), new Complex32(1.0f, 1), new Complex32(2.0f, 1) }); @@ -171,7 +171,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [MultipleAsserts] public void DetermineStatusWithResidualNaN() { - var criterium = new ResidualStopCriterium(1e-3, 10); + var criterium = new ResidualStopCriterium(1e-3f, 10); Assert.IsNotNull(criterium, "There should be a criterium"); var solution = new DenseVector(new[] { new Complex32(1.0f, 1), new Complex32(1.0f, 1), new Complex32(2.0f, 1) }); @@ -204,7 +204,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [MultipleAsserts] public void DetermineStatus() { - var criterium = new ResidualStopCriterium(1e-3, 10); + var criterium = new ResidualStopCriterium(1e-3f, 10); Assert.IsNotNull(criterium, "There should be a criterium"); // Note that the solution vector isn't actually being used so ... @@ -227,7 +227,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [MultipleAsserts] public void ResetCalculationState() { - var criterium = new ResidualStopCriterium(1e-3, 10); + var criterium = new ResidualStopCriterium(1e-3f, 10); Assert.IsNotNull(criterium, "There should be a criterium"); var solution = new DenseVector(new[] { new Complex32(0.001f, 1), new Complex32(0.001f, 1), new Complex32(0.002f, 1) }); @@ -245,7 +245,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr [MultipleAsserts] public void Clone() { - var criterium = new ResidualStopCriterium(1e-3, 10); + var criterium = new ResidualStopCriterium(1e-3f, 10); Assert.IsNotNull(criterium, "There should be a criterium"); var clone = criterium.Clone(); diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedVectorTests.cs index bc7e8883..822223f1 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedVectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedVectorTests.cs @@ -3,9 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// // Copyright (c) 2009-2010 Math.NET -// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +12,8 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -30,14 +26,12 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 { - using System; - using Numerics; using System.Collections.Generic; - using Distributions; + using LinearAlgebra.Complex32; using LinearAlgebra.Generic; - using Threading; + using Complex32 = Numerics.Complex32; - internal class UserDefinedVector : Vector + internal class UserDefinedVector : Vector { private readonly Complex32[] _data; @@ -75,128 +69,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 { return new UserDefinedVector(size); } - - public override Vector Negate() - { - var result = new UserDefinedVector(Count); - CommonParallel.For( - 0, - _data.Length, - index => result[index] = -_data[index]); - - return result; - } - - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException("length"); - } - - var v = (UserDefinedVector)CreateVector(length); - for (var index = 0; index < v._data.Length; index++) - { - v._data[index] = new Complex32((float)randomDistribution.Sample(), (float)randomDistribution.Sample()); - } - - return v; - } - - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException("length"); - } - - var v = (UserDefinedVector)CreateVector(length); - for (var index = 0; index < v._data.Length; index++) - { - v._data[index] = new Complex32(randomDistribution.Sample(), randomDistribution.Sample()); - } - - return v; - } - - public override int MinimumIndex() - { - throw new NotSupportedException(); - } - - public override int MaximumIndex() - { - throw new NotSupportedException(); - } - - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(1.0f / (float)norm, clone); - - return clone; - } - - protected sealed override Complex32 AddT(Complex32 val1, Complex32 val2) - { - return val1 + val2; - } - - protected sealed override Complex32 SubtractT(Complex32 val1, Complex32 val2) - { - return val1 - val2; - } - - protected sealed override Complex32 MultiplyT(Complex32 val1, Complex32 val2) - { - return val1 * val2; - } - - protected sealed override Complex32 DivideT(Complex32 val1, Complex32 val2) - { - return val1 / val2; - } - - protected sealed override double AbsoluteT(Complex32 val1) - { - return val1.Magnitude; - } - - public override void Conjugate(Vector target) - { - if (target == null) - { - throw new ArgumentNullException("target"); - } - - if (Count != target.Count) - { - throw new ArgumentException("target"); - } - - if (ReferenceEquals(this, target)) - { - var tmp = CreateVector(Count); - Conjugate(tmp); - tmp.CopyTo(target); - } - - CommonParallel.For( - 0, - Count, - index => target[index] = this[index].Conjugate()); - } } public class UserDefinedVectorTests : VectorTests @@ -217,4 +89,4 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 return vector; } } -} \ No newline at end of file +} diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/VectorTests.Norm.cs b/src/UnitTests/LinearAlgebraTests/Complex32/VectorTests.Norm.cs index eb000bd1..4f7085a9 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/VectorTests.Norm.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/VectorTests.Norm.cs @@ -40,21 +40,21 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 public void CanComputeNorm() { var vector = CreateVector(Data); - AssertHelpers.AlmostEqual(7.7459666f, (float)vector.Norm(2), 7); + AssertHelpers.AlmostEqual(7.7459666f, vector.Norm(2).Real, 7); } [Test] public void CanComputeNorm1() { var vector = CreateVector(Data); - AssertHelpers.AlmostEqual(16.0346843f, (float)vector.Norm(1), 7); + AssertHelpers.AlmostEqual(16.0346843f, vector.Norm(1).Real, 7); } [Test] public void CanComputeSquareNorm() { var vector = CreateVector(Data); - AssertHelpers.AlmostEqual(60f, (float)vector.Norm(2) * (float)vector.Norm(2), 6); + AssertHelpers.AlmostEqual(60f, vector.Norm(2).Real * vector.Norm(2).Real, 6); } [Test] @@ -65,14 +65,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 public void CanComputeNormP(int p, float expected) { var vector = CreateVector(Data); - AssertHelpers.AlmostEqual(expected, (float)vector.Norm(p), 7); + AssertHelpers.AlmostEqual(expected, vector.Norm(p).Real, 7); } [Test] public void CanComputeNormInfinity() { var vector = CreateVector(Data); - AssertHelpers.AlmostEqual(5.0990195, (float)vector.Norm(Double.PositiveInfinity), 7); + AssertHelpers.AlmostEqual(5.0990195, vector.Norm(Double.PositiveInfinity).Real, 7); } [Test] diff --git a/src/UnitTests/LinearAlgebraTests/Double/UserDefinedVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Double/UserDefinedVectorTests.cs index 7b21fe80..37c51450 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/UserDefinedVectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/UserDefinedVectorTests.cs @@ -3,9 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// // Copyright (c) 2009-2010 Math.NET -// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +12,8 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -30,14 +26,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double { - using System; using System.Collections.Generic; - using Distributions; + using LinearAlgebra.Double; using LinearAlgebra.Generic; - using Properties; - using Threading; - internal class UserDefinedVector : Vector + internal class UserDefinedVector : Vector { private readonly double[] _data; @@ -75,125 +68,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double { return new UserDefinedVector(size); } - - public override Vector Negate() - { - var result = new UserDefinedVector(Count); - CommonParallel.For( - 0, - _data.Length, - index => result[index] = -_data[index]); - - return result; - } - - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (UserDefinedVector)CreateVector(length); - for (var index = 0; index < v._data.Length; index++) - { - v._data[index] = randomDistribution.Sample(); - } - - return v; - } - - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (UserDefinedVector)CreateVector(length); - for (var index = 0; index < v._data.Length; index++) - { - v._data[index] = randomDistribution.Sample(); - } - - return v; - } - - public override int MinimumIndex() - { - var index = 0; - var min = _data[0]; - for (var i = 1; i < Count; i++) - { - if (min > _data[i]) - { - index = i; - min = _data[i]; - } - } - - return index; - } - - public override int MaximumIndex() - { - var index = 0; - var max = _data[0]; - for (var i = 1; i < Count; i++) - { - if (max < _data[i]) - { - index = i; - max = _data[i]; - } - } - - return index; - } - - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(1.0 / norm, clone); - - return clone; - } - - protected sealed override double AddT(double val1, double val2) - { - return val1 + val2; - } - - protected sealed override double SubtractT(double val1, double val2) - { - return val1 - val2; - } - - protected sealed override double MultiplyT(double val1, double val2) - { - return val1 * val2; - } - - protected sealed override double DivideT(double val1, double val2) - { - return val1 / val2; - } - - protected sealed override double AbsoluteT(double val1) - { - return Math.Abs(val1); - } } public class UserDefinedVectorTests : VectorTests @@ -214,4 +88,4 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double return vector; } } -} \ No newline at end of file +} diff --git a/src/UnitTests/LinearAlgebraTests/Single/Solvers/StopCriterium/ResidualStopCriteriumTest.cs b/src/UnitTests/LinearAlgebraTests/Single/Solvers/StopCriterium/ResidualStopCriteriumTest.cs index bacdda85..83b8931e 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Solvers/StopCriterium/ResidualStopCriteriumTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Solvers/StopCriterium/ResidualStopCriteriumTest.cs @@ -12,7 +12,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [ExpectedArgumentOutOfRangeException] public void CreateWithNegativeMaximum() { - new ResidualStopCriterium(-0.1); + new ResidualStopCriterium(-0.1f); Assert.Fail(); } @@ -28,10 +28,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [MultipleAsserts] public void Create() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); - Assert.AreEqual(1e-8, criterium.Maximum, "Incorrect maximum"); + Assert.AreEqual(1e-6f, criterium.Maximum, "Incorrect maximum"); Assert.AreEqual(50, criterium.MinimumIterationsBelowMaximum, "Incorrect iteration count"); } @@ -39,7 +39,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [MultipleAsserts] public void ResetMaximum() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.ResetMaximumResidualToDefault(); @@ -50,7 +50,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [MultipleAsserts] public void ResetMinimumIterationsBelowMaximum() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.ResetMinimumIterationsBelowMaximumToDefault(); @@ -61,7 +61,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [ExpectedArgumentOutOfRangeException] public void DetermineStatusWithIllegalIterationNumber() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(-1, @@ -75,7 +75,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [ExpectedArgumentNullException] public void DetermineStatusWithNullSolutionVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -89,7 +89,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [ExpectedArgumentNullException] public void DetermineStatusWithNullSourceVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -103,7 +103,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [ExpectedArgumentNullException] public void DetermineStatusWithNullResidualVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -116,7 +116,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [ExpectedArgumentException] public void DetermineStatusWithNonMatchingSolutionVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -129,7 +129,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [ExpectedArgumentException] public void DetermineStatusWithNonMatchingSourceVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -142,7 +142,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [ExpectedArgumentException] public void DetermineStatusWithNonMatchingResidualVector() { - var criterium = new ResidualStopCriterium(1e-8, 50); + var criterium = new ResidualStopCriterium(1e-6f, 50); Assert.IsNotNull(criterium, "There should be a criterium"); criterium.DetermineStatus(1, @@ -155,7 +155,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [MultipleAsserts] public void DetermineStatusWithSourceNaN() { - var criterium = new ResidualStopCriterium(1e-3, 10); + var criterium = new ResidualStopCriterium(1e-3f, 10); Assert.IsNotNull(criterium, "There should be a criterium"); var solution = new DenseVector(new[] { 1.0f, 1.0f, 2.0f }); @@ -170,7 +170,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [MultipleAsserts] public void DetermineStatusWithResidualNaN() { - var criterium = new ResidualStopCriterium(1e-3, 10); + var criterium = new ResidualStopCriterium(1e-3f, 10); Assert.IsNotNull(criterium, "There should be a criterium"); var solution = new DenseVector(new[] { 1.0f, 1.0f, 2.0f }); @@ -203,7 +203,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [MultipleAsserts] public void DetermineStatus() { - var criterium = new ResidualStopCriterium(1e-3, 10); + var criterium = new ResidualStopCriterium(1e-3f, 10); Assert.IsNotNull(criterium, "There should be a criterium"); // Note that the solution vector isn't actually being used so ... @@ -226,7 +226,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [MultipleAsserts] public void ResetCalculationState() { - var criterium = new ResidualStopCriterium(1e-3, 10); + var criterium = new ResidualStopCriterium(1e-3f, 10); Assert.IsNotNull(criterium, "There should be a criterium"); var solution = new DenseVector(new[] { 0.001f, 0.001f, 0.002f }); @@ -244,7 +244,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.StopCrite [MultipleAsserts] public void Clone() { - var criterium = new ResidualStopCriterium(1e-3, 10); + var criterium = new ResidualStopCriterium(1e-3f, 10); Assert.IsNotNull(criterium, "There should be a criterium"); var clone = criterium.Clone(); diff --git a/src/UnitTests/LinearAlgebraTests/Single/UserDefinedVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Single/UserDefinedVectorTests.cs index 6e8b316e..7b484eef 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/UserDefinedVectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/UserDefinedVectorTests.cs @@ -35,10 +35,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single using System.Linq; using Distributions; using LinearAlgebra.Generic; + using LinearAlgebra.Single; using Properties; using Threading; - internal class UserDefinedVector : Vector + internal class UserDefinedVector : Vector { private readonly float[] _data; @@ -76,162 +77,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single { return new UserDefinedVector(size); } - - public override Vector Negate() - { - var result = new UserDefinedVector(Count); - CommonParallel.For( - 0, - _data.Length, - index => result[index] = -_data[index]); - - return result; - } - - public override Vector Random(int length, IContinuousDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (UserDefinedVector)CreateVector(length); - for (var index = 0; index < v._data.Length; index++) - { - v._data[index] = (float)randomDistribution.Sample(); - } - - return v; - } - - public override Vector Random(int length, IDiscreteDistribution randomDistribution) - { - if (length < 1) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "length"); - } - - var v = (UserDefinedVector)CreateVector(length); - for (var index = 0; index < v._data.Length; index++) - { - v._data[index] = randomDistribution.Sample(); - } - - return v; - } - - public override int MinimumIndex() - { - var index = 0; - var min = _data[0]; - for (var i = 1; i < Count; i++) - { - if (min > _data[i]) - { - index = i; - min = _data[i]; - } - } - - return index; - } - - public override int MaximumIndex() - { - var index = 0; - var max = _data[0]; - for (var i = 1; i < Count; i++) - { - if (max < _data[i]) - { - index = i; - max = _data[i]; - } - } - - return index; - } - - public override double Norm(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - if (1.0 == p) - { - return CommonParallel.Aggregate( - 0, - Count, - index => Math.Abs(this[index])); - } - - if (2.0 == p) - { - return _data.Aggregate(0.0f, SpecialFunctions.Hypotenuse); - } - - if (Double.IsPositiveInfinity(p)) - { - return CommonParallel.Select( - 0, - Count, - (index, localData) => Math.Max(localData, Math.Abs(_data[index])), - Math.Max); - } - - var sum = CommonParallel.Aggregate( - 0, - Count, - index => Math.Pow(Math.Abs(_data[index]), p)); - - return Math.Pow(sum, 1.0 / p); - } - - public override Vector Normalize(double p) - { - if (p < 0.0) - { - throw new ArgumentOutOfRangeException("p"); - } - - var norm = Norm(p); - var clone = Clone(); - if (norm == 0.0) - { - return clone; - } - - clone.Multiply(1.0f / (float)norm, clone); - - return clone; - } - - protected sealed override float AddT(float val1, float val2) - { - return val1 + val2; - } - - protected sealed override float SubtractT(float val1, float val2) - { - return val1 - val2; - } - - protected sealed override float MultiplyT(float val1, float val2) - { - return val1 * val2; - } - - protected sealed override float DivideT(float val1, float val2) - { - return val1 / val2; - } - - protected sealed override double AbsoluteT(float val1) - { - return Math.Abs(val1); - } } public class UserDefinedVectorTests : VectorTests diff --git a/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Norm.cs b/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Norm.cs index 33d4fd15..36719b03 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Norm.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Norm.cs @@ -80,11 +80,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single { var vector = CreateVector(Data); var result = vector.Normalize(2); - AssertHelpers.AlmostEqual(0.134839967f, result[0], 7); - AssertHelpers.AlmostEqual(0.269679934f, result[1], 7); - AssertHelpers.AlmostEqual(0.404519916f, result[2], 7); - AssertHelpers.AlmostEqual(0.539359868f, result[3], 7); - AssertHelpers.AlmostEqual(0.6741998f, result[4], 7); + AssertHelpers.AlmostEqual(0.134839967f, result[0], 6); + AssertHelpers.AlmostEqual(0.269679934f, result[1], 6); + AssertHelpers.AlmostEqual(0.404519916f, result[2], 6); + AssertHelpers.AlmostEqual(0.539359868f, result[3], 6); + AssertHelpers.AlmostEqual(0.6741998f, result[4], 6); } } } \ No newline at end of file