// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // // Copyright (c) 2009-2013 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. // using System; using System.Collections.Generic; using System.Diagnostics; using System.Globalization; using System.Linq; using IStation.Numerics.Distributions; using IStation.Numerics.LinearAlgebra.Storage; using IStation.Numerics.Providers.LinearAlgebra; using IStation.Numerics.Threading; namespace IStation.Numerics.LinearAlgebra.Double { /// /// A vector using dense storage. /// [Serializable] [DebuggerDisplay("DenseVector {" + nameof(Count) + "}-Double")] public class DenseVector : Vector { /// /// Number of elements /// readonly int _length; /// /// Gets the vector's data. /// readonly double[] _values; /// /// Create a new dense vector straight from an initialized vector storage instance. /// The storage is used directly without copying. /// Intended for advanced scenarios where you're working directly with /// storage for performance or interop reasons. /// public DenseVector(DenseVectorStorage storage) : base(storage) { _length = storage.Length; _values = storage.Data; } /// /// Create a new dense vector with the given length. /// All cells of the vector will be initialized to zero. /// /// If length is less than one. public DenseVector(int length) : this(new DenseVectorStorage(length)) { } /// /// Create a new dense vector directly binding to a raw array. /// The array is used directly without copying. /// Very efficient, but changes to the array and the vector will affect each other. /// public DenseVector(double[] storage) : this(new DenseVectorStorage(storage.Length, storage)) { } /// /// Create a new dense vector as a copy of the given other vector. /// This new vector will be independent from the other vector. /// A new memory block will be allocated for storing the vector. /// public static DenseVector OfVector(Vector vector) { return new DenseVector(DenseVectorStorage.OfVector(vector.Storage)); } /// /// Create a new dense vector as a copy of the given array. /// This new vector will be independent from the array. /// A new memory block will be allocated for storing the vector. /// public static DenseVector OfArray(double[] array) { return new DenseVector(DenseVectorStorage.OfVector(new DenseVectorStorage(array.Length, array))); } /// /// Create a new dense vector as a copy of the given enumerable. /// This new vector will be independent from the enumerable. /// A new memory block will be allocated for storing the vector. /// public static DenseVector OfEnumerable(IEnumerable enumerable) { return new DenseVector(DenseVectorStorage.OfEnumerable(enumerable)); } /// /// Create a new dense vector as a copy of the given indexed enumerable. /// Keys must be provided at most once, zero is assumed if a key is omitted. /// This new vector will be independent from the enumerable. /// A new memory block will be allocated for storing the vector. /// public static DenseVector OfIndexedEnumerable(int length, IEnumerable> enumerable) { return new DenseVector(DenseVectorStorage.OfIndexedEnumerable(length, enumerable)); } /// /// Create a new dense vector and initialize each value using the provided value. /// public static DenseVector Create(int length, double value) { if (value == 0d) return new DenseVector(length); return new DenseVector(DenseVectorStorage.OfValue(length, value)); } /// /// Create a new dense vector and initialize each value using the provided init function. /// public static DenseVector Create(int length, Func init) { return new DenseVector(DenseVectorStorage.OfInit(length, init)); } /// /// Create a new dense vector with values sampled from the provided random distribution. /// public static DenseVector CreateRandom(int length, IContinuousDistribution distribution) { var samples = Generate.Random(length, distribution); return new DenseVector(new DenseVectorStorage(length, samples)); } /// /// Gets the vector's data. /// /// The vector's data. public double[] Values => _values; /// /// Returns a reference to the internal data structure. /// /// The DenseVector whose internal data we are /// returning. /// /// A reference to the internal date of the given vector. /// public static explicit operator double[](DenseVector vector) { if (vector == null) { throw new ArgumentNullException(nameof(vector)); } return vector.Values; } /// /// Returns a vector bound directly to a reference of the provided array. /// /// The array to bind to the DenseVector object. /// /// A DenseVector whose values are bound to the given array. /// public static implicit operator DenseVector(double[] array) { if (array == null) { throw new ArgumentNullException(nameof(array)); } return new DenseVector(array); } /// /// 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) { if (result is DenseVector dense) { CommonParallel.For(0, _values.Length, 4096, (a, b) => { for (int i = a; i < b; i++) { dense._values[i] = _values[i] + scalar; } }); } else { base.DoAdd(scalar, 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 override void DoAdd(Vector other, Vector result) { if (other is DenseVector otherDense && result is DenseVector resultDense) { LinearAlgebraControl.Provider.AddArrays(_values, otherDense._values, resultDense._values); } else { base.DoAdd(other, result); } } /// /// Adds two Vectors together and returns the results. /// /// One of the vectors to add. /// The other vector to add. /// The result of the addition. /// If and are not the same size. /// If or is . public static DenseVector operator +(DenseVector leftSide, DenseVector rightSide) { if (rightSide == null) { throw new ArgumentNullException(nameof(rightSide)); } if (leftSide == null) { throw new ArgumentNullException(nameof(leftSide)); } if (leftSide.Count != rightSide.Count) { throw new ArgumentException("All vectors must have the same dimensionality.", nameof(rightSide)); } return (DenseVector)leftSide.Add(rightSide); } /// /// 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) { if (result is DenseVector dense) { CommonParallel.For(0, _values.Length, 4096, (a, b) => { for (int i = a; i < b; i++) { dense._values[i] = _values[i] - scalar; } }); } else { base.DoSubtract(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) { if (other is DenseVector otherDense && result is DenseVector resultDense) { LinearAlgebraControl.Provider.SubtractArrays(_values, otherDense._values, resultDense._values); } else { base.DoSubtract(other, result); } } /// /// Returns a Vector containing the negated values of . /// /// The vector to get the values from. /// A vector containing the negated values as . /// If is . public static DenseVector operator -(DenseVector rightSide) { if (rightSide == null) { throw new ArgumentNullException(nameof(rightSide)); } return (DenseVector)rightSide.Negate(); } /// /// Subtracts two Vectors and returns the results. /// /// The vector to subtract from. /// The vector to subtract. /// The result of the subtraction. /// If and are not the same size. /// If or is . public static DenseVector operator -(DenseVector leftSide, DenseVector rightSide) { if (leftSide == null) { throw new ArgumentNullException(nameof(leftSide)); } return (DenseVector)leftSide.Subtract(rightSide); } /// /// Negates vector and saves result to /// /// Target vector protected override void DoNegate(Vector result) { if (result is DenseVector denseResult) { LinearAlgebraControl.Provider.ScaleArray(-1.0d, _values, denseResult.Values); } else { base.DoNegate(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 override void DoMultiply(double scalar, Vector result) { if (result is DenseVector denseResult) { LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult.Values); } else { base.DoMultiply(scalar, result); } } /// /// Computes the dot product between this vector and another vector. /// /// The other vector. /// The sum of a[i]*b[i] for all i. protected override double DoDotProduct(Vector other) { return other is DenseVector denseVector ? LinearAlgebraControl.Provider.DotProduct(_values, denseVector.Values) : base.DoDotProduct(other); } /// /// Multiplies a vector with a scalar. /// /// The vector to scale. /// The scalar value. /// The result of the multiplication. /// If is . public static DenseVector operator *(DenseVector leftSide, double rightSide) { if (leftSide == null) { throw new ArgumentNullException(nameof(leftSide)); } return (DenseVector)leftSide.Multiply(rightSide); } /// /// Multiplies a vector with a scalar. /// /// The scalar value. /// The vector to scale. /// The result of the multiplication. /// If is . public static DenseVector operator *(double leftSide, DenseVector rightSide) { if (rightSide == null) { throw new ArgumentNullException(nameof(rightSide)); } return (DenseVector)rightSide.Multiply(leftSide); } /// /// Computes the dot product between two Vectors. /// /// The left row vector. /// The right column vector. /// The dot product between the two vectors. /// If and are not the same size. /// If or is . public static double operator *(DenseVector leftSide, DenseVector rightSide) { if (leftSide == null) { throw new ArgumentNullException(nameof(leftSide)); } return leftSide.DotProduct(rightSide); } /// /// Divides a vector with a scalar. /// /// The vector to divide. /// The scalar value. /// The result of the division. /// If is . public static DenseVector operator /(DenseVector leftSide, double rightSide) { if (leftSide == null) { throw new ArgumentNullException(nameof(leftSide)); } return (DenseVector)leftSide.Divide(rightSide); } /// /// Computes the canonical modulus, where the result has the sign of the divisor, /// for each element of the vector for the given divisor. /// /// The divisor to use. /// A vector to store the results in. protected override void DoModulus(double divisor, Vector result) { if (result is DenseVector dense) { CommonParallel.For(0, _length, 4096, (a, b) => { for (int i = a; i < b; i++) { dense._values[i] = Euclid.Modulus(_values[i], divisor); } }); } else { base.DoModulus(divisor, result); } } /// /// Computes the remainder (% operator), where the result has the sign of the dividend, /// for each element of the vector for the given divisor. /// /// The divisor to use. /// A vector to store the results in. protected override void DoRemainder(double divisor, Vector result) { if (result is DenseVector dense) { CommonParallel.For(0, _length, 4096, (a, b) => { for (int i = a; i < b; i++) { dense._values[i] = _values[i] % divisor; } }); } else { base.DoRemainder(divisor, result); } } /// /// Computes the remainder (% operator), where the result has the sign of the dividend, /// of each element of the vector of the given divisor. /// /// The vector whose elements we want to compute the remainder of. /// The divisor to use, /// If is . public static DenseVector operator %(DenseVector leftSide, double rightSide) { if (leftSide == null) { throw new ArgumentNullException(nameof(leftSide)); } return (DenseVector)leftSide.Remainder(rightSide); } /// /// 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(_values[index]); for (var i = 1; i < _length; i++) { var test = Math.Abs(_values[i]); if (test < min) { index = i; min = test; } } return index; } /// /// 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(_values[index]); for (var i = 1; i < _length; i++) { var test = Math.Abs(_values[i]); if (test > max) { index = i; max = test; } } return index; } /// /// Returns the index of the maximum element. /// /// The index of maximum element. public override int MaximumIndex() { var index = 0; var max = _values[0]; for (var i = 1; i < _length; i++) { if (max < _values[i]) { index = i; max = _values[i]; } } return index; } /// /// Returns the index of the minimum element. /// /// The index of minimum element. public override int MinimumIndex() { var index = 0; var min = _values[0]; for (var i = 1; i < _length; i++) { if (min > _values[i]) { index = i; min = _values[i]; } } return index; } /// /// Computes the sum of the vector's elements. /// /// The sum of the vector's elements. public override double Sum() { var sum = 0.0; for (var index = 0; index < _length; index++) { sum += _values[index]; } return sum; } /// /// Calculates the L1 norm of the vector, also known as Manhattan norm. /// /// The sum of the absolute values. public override double L1Norm() { var sum = 0d; for (var index = 0; index < _length; index++) { sum += Math.Abs(_values[index]); } return sum; } /// /// Calculates the L2 norm of the vector, also known as Euclidean norm. /// /// The square root of the sum of the squared values. public override double L2Norm() { // TODO: native provider return _values.Aggregate(0d, SpecialFunctions.Hypotenuse); } /// /// Calculates the infinity norm of the vector. /// /// The maximum absolute value. public override double InfinityNorm() { return CommonParallel.Aggregate(_values, (i, v) => Math.Abs(v), Math.Max, 0d); } /// /// Computes the p-Norm. /// /// The p value. /// Scalar ret = ( ∑|this[i]|^p )^(1/p) public override double Norm(double p) { if (p < 0d) throw new ArgumentOutOfRangeException(nameof(p)); if (p == 1d) return L1Norm(); if (p == 2d) return L2Norm(); if (double.IsPositiveInfinity(p)) return InfinityNorm(); var sum = 0d; for (var index = 0; index < _length; index++) { sum += Math.Pow(Math.Abs(_values[index]), p); } return Math.Pow(sum, 1.0 / p); } /// /// 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 DoPointwiseMultiply(Vector other, Vector result) { if (other is DenseVector denseOther && result is DenseVector denseResult) { LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values); } else { base.DoPointwiseMultiply(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 override void DoPointwiseDivide(Vector divisor, Vector result) { if (divisor is DenseVector denseOther && result is DenseVector denseResult) { LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseOther._values, denseResult._values); } else { base.DoPointwiseDivide(divisor, result); } } /// /// Pointwise raise this vector to an exponent vector and store the result into the result vector. /// /// The exponent vector to raise this vector values to. /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Vector exponent, Vector result) { if (exponent is DenseVector denseExponent && result is DenseVector denseResult) { LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values); } else { base.DoPointwisePower(exponent, result); } } #region Parse Functions /// /// Creates a double dense vector based on a string. The string can be in the following formats (without the /// quotes): 'n', 'n,n,..', '(n,n,..)', '[n,n,...]', where n is a double. /// /// /// A double dense vector containing the values specified by the given string. /// /// /// the string to parse. /// /// /// An that supplies culture-specific formatting information. /// public static DenseVector Parse(string value, IFormatProvider formatProvider = null) { if (value == null) { throw new ArgumentNullException(nameof(value)); } value = value.Trim(); if (value.Length == 0) { throw new FormatException(); } // strip out parens if (value.StartsWith("(", StringComparison.Ordinal)) { if (!value.EndsWith(")", StringComparison.Ordinal)) { throw new FormatException(); } value = value.Substring(1, value.Length - 2).Trim(); } if (value.StartsWith("[", StringComparison.Ordinal)) { if (!value.EndsWith("]", StringComparison.Ordinal)) { throw new FormatException(); } value = value.Substring(1, value.Length - 2).Trim(); } // parsing var tokens = value.Split(new[] {formatProvider.GetTextInfo().ListSeparator, " ", "\t"}, StringSplitOptions.RemoveEmptyEntries); var data = tokens.Select(t => double.Parse(t, NumberStyles.Any, formatProvider)).ToArray(); if (data.Length == 0) throw new FormatException(); return new DenseVector(data); } /// /// Converts the string representation of a real dense vector to double-precision dense vector equivalent. /// A return value indicates whether the conversion succeeded or failed. /// /// /// A string containing a real vector to convert. /// /// /// The parsed value. /// /// /// If the conversion succeeds, the result will contain a complex number equivalent to value. /// Otherwise the result will be null. /// public static bool TryParse(string value, out DenseVector result) { return TryParse(value, null, out result); } /// /// Converts the string representation of a real dense vector to double-precision dense vector equivalent. /// A return value indicates whether the conversion succeeded or failed. /// /// /// A string containing a real vector to convert. /// /// /// An that supplies culture-specific formatting information about value. /// /// /// The parsed value. /// /// /// If the conversion succeeds, the result will contain a complex number equivalent to value. /// Otherwise the result will be null. /// public static bool TryParse(string value, IFormatProvider formatProvider, out DenseVector result) { bool ret; try { result = Parse(value, formatProvider); ret = true; } catch (ArgumentNullException) { result = null; ret = false; } catch (FormatException) { result = null; ret = false; } return ret; } #endregion } }