Skip to content

Commit

Permalink
IRID-174: Force QR decomposition to be unique (positive real R diagonal)
Browse files Browse the repository at this point in the history 
IRID-173: Provide permutation matrix on LU decomposition
IRID-172: Complex Number: Unitary - normalize to unit circle
IRID-151: Distribution-specific parameter estimation (PARTIAL ONLY)
Fixing StyleCop rule SA1116, SA1106, SA1107, SA1108
  • Loading branch information
@cdrnet
cdrnet committed Jun 15, 2008
1 parent 6c4ea18 commit 6c86800
Show file tree
Hide file tree
Showing 16 changed files with 355 additions and 118 deletions.
60 changes: 40 additions & 20 deletions src/app/MathNet.Iridium/Library/Collection.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -313,10 +313,14 @@ public object SyncRoot
/// concatenated collection is not built explicitly.</remarks>
public static ICollection Concat(ICollection c1, ICollection c2)
{
if(c1 == null) throw new ArgumentNullException("c1",
string.Format(Resources.ArgumentNull, "c1"));
if(c2 == null) throw new ArgumentNullException("c2",
string.Format(Resources.ArgumentNull, "c2"));
if(c1 == null)
{
throw new ArgumentNullException("c1", string.Format(Resources.ArgumentNull, "c1"));
}
if(c2 == null)
{
throw new ArgumentNullException("c2", string.Format(Resources.ArgumentNull, "c2"));
}

return new ConcatCollection(c1, c2);
}
Expand All @@ -336,10 +340,14 @@ public static ICollection Concat(ICollection c1, ICollection c2)
[Obsolete("Use Set<T>.Intersect() instead.", false)]
public static ICollection Inter(ICollection c1, ICollection c2)
{
if(c1 == null) throw new ArgumentNullException("c1",
string.Format(Resources.ArgumentNull, "c1"));
if(c2 == null) throw new ArgumentNullException("c2",
string.Format(Resources.ArgumentNull, "c2"));
if(c1 == null)
{
throw new ArgumentNullException("c1", string.Format(Resources.ArgumentNull, "c1"));
}
if(c2 == null)
{
throw new ArgumentNullException("c2", string.Format(Resources.ArgumentNull, "c2"));
}

return new InterCollection(c1, c2);
}
Expand All @@ -354,10 +362,14 @@ public static ICollection Inter(ICollection c1, ICollection c2)
[Obsolete("Use Set<T>.Subtract() instead.", false)]
public static ICollection Minus(ICollection c1, ICollection c2)
{
if(c1 == null) throw new ArgumentNullException("c1",
string.Format(Resources.ArgumentNull, "c1"));
if(c2 == null) throw new ArgumentNullException("c2",
string.Format(Resources.ArgumentNull, "c2"));
if(c1 == null)
{
throw new ArgumentNullException("c1", string.Format(Resources.ArgumentNull, "c1"));
}
if(c2 == null)
{
throw new ArgumentNullException("c2", string.Format(Resources.ArgumentNull, "c2"));
}

return new MinusCollection(c1, c2);
}
Expand All @@ -370,10 +382,14 @@ public static ICollection Minus(ICollection c1, ICollection c2)
/// <param name="c2">Should not be null.</param>
public static ICollection Product(ICollection c1, ICollection c2)
{
if(c1 == null) throw new ArgumentNullException("c1",
string.Format(Resources.ArgumentNull, "c1"));
if(c2 == null) throw new ArgumentNullException("c2",
string.Format(Resources.ArgumentNull, "c2"));
if(c1 == null)
{
throw new ArgumentNullException("c1", string.Format(Resources.ArgumentNull, "c1"));
}
if(c2 == null)
{
throw new ArgumentNullException("c2", string.Format(Resources.ArgumentNull, "c2"));
}

return null;
}
Expand All @@ -393,10 +409,14 @@ public static ICollection Product(ICollection c1, ICollection c2)
[Obsolete("Use Set<T>.Union() instead.", false)]
public static ICollection Union(ICollection c1, ICollection c2)
{
if(c1 == null) throw new ArgumentNullException("c1",
string.Format(Resources.ArgumentNull, "c1"));
if(c2 == null) throw new ArgumentNullException("c2",
string.Format(Resources.ArgumentNull, "c2"));
if(c1 == null)
{
throw new ArgumentNullException("c1", string.Format(Resources.ArgumentNull, "c1"));
}
if(c2 == null)
{
throw new ArgumentNullException("c2", string.Format(Resources.ArgumentNull, "c2"));
}

return new UnionCollection(c1, c2);
}
Expand Down
6 changes: 4 additions & 2 deletions src/app/MathNet.Iridium/Library/Combinatorics.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,8 +216,9 @@ int k
)
{
bool[] selection = new bool[n];
if(k * 3 < n) // just pick and try
if(k * 3 < n)
{
// just pick and try
int selectionCount = 0;
while(selectionCount < k)
{
Expand All @@ -230,8 +231,9 @@ int k
}
return selection;
}
else // based on permutation
else
{
// based on permutation
int[] permutation = RandomPermutation(n);
for(int i = 0; i < k; i++)
{
Expand Down
80 changes: 35 additions & 45 deletions src/app/MathNet.Iridium/Library/Complex.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -147,51 +147,6 @@ public static IComparer ArgumentModulusComparer
double real;
double imag;

#region Normalization

// TODO: Never called, remove?
//void
//NormalizeToUnityOrNull()
//{
// if(double.IsPositiveInfinity(real) && double.IsPositiveInfinity(imag))
// {
// real = Constants.Sqrt1_2;
// imag = Constants.Sqrt1_2;
// }
// else if(double.IsPositiveInfinity(real) && double.IsNegativeInfinity(imag))
// {
// real = Constants.Sqrt1_2;
// imag = -Constants.Sqrt1_2;
// }
// else if(double.IsNegativeInfinity(real) && double.IsPositiveInfinity(imag))
// {
// real = -Constants.Sqrt1_2;
// imag = -Constants.Sqrt1_2;
// }
// else if(double.IsNegativeInfinity(real) && double.IsNegativeInfinity(imag))
// {
// real = -Constants.Sqrt1_2;
// imag = Constants.Sqrt1_2;
// }
// else
// {
// //Don't replace this with "Modulus"!
// double mod = Math.Sqrt(real * real + imag * imag);
// if(mod == 0)
// {
// real = 0;
// imag = 0;
// }
// else
// {
// real = real / mod;
// imag = imag / mod;
// }
// }
//}

#endregion

#region Constructors and Constants

/// <summary>
Expand Down Expand Up @@ -585,6 +540,41 @@ public double Argument
}
}

/// <summary>
/// Gets the unity of this complex (same argument, but on the unit circle; exp(I*arg))
/// </summary>
public Complex Unity
{
get
{
if(double.IsPositiveInfinity(real) && double.IsPositiveInfinity(imag))
{
return new Complex(Constants.Sqrt1_2, Constants.Sqrt1_2);
}
if(double.IsPositiveInfinity(real) && double.IsNegativeInfinity(imag))
{
return new Complex(Constants.Sqrt1_2, -Constants.Sqrt1_2);
}
if(double.IsNegativeInfinity(real) && double.IsPositiveInfinity(imag))
{
return new Complex(-Constants.Sqrt1_2, -Constants.Sqrt1_2);
}
if(double.IsNegativeInfinity(real) && double.IsNegativeInfinity(imag))
{
return new Complex(-Constants.Sqrt1_2, Constants.Sqrt1_2);
}

//Don't replace this with "Modulus"!
double mod = Fn.Hypot(real, imag);
if(mod == 0)
{
return Complex.Zero;
}

return new Complex(real / mod, imag / mod);
}
}

#endregion

/// <summary>
Expand Down
27 changes: 26 additions & 1 deletion src/app/MathNet.Iridium/Library/Distributions/Continuous/NormalDistribution.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,13 @@
#endregion

using System;
using System.Collections.Generic;
using MathNet.Numerics.RandomSources;
using MathNet.Numerics.Statistics;

namespace MathNet.Numerics.Distributions
{

/// <summary>
/// Pseudo-random generation of normal distributed deviates.
/// </summary>
Expand All @@ -50,6 +53,7 @@ public sealed class NormalDistribution : ContinuousDistribution
StandardDistribution _standard;

#region Construction

/// <summary>
/// Initializes a new instance, using a <see cref="SystemRandomSource"/>
/// as underlying random number generator.
Expand Down Expand Up @@ -94,6 +98,7 @@ double sigma
_standard = new StandardDistribution(this.RandomSource);
SetDistributionParameters(mu, sigma);
}

#endregion

/// <summary>
Expand All @@ -111,6 +116,7 @@ public override RandomSource RandomSource
}

#region Distribution Parameters

/// <summary>
/// Gets or sets the mu parameter.
/// </summary>
Expand All @@ -121,7 +127,7 @@ public double Mu
}

/// <summary>
/// Gets or sets the sigma parameter.
/// Gets or sets the sigma (standard deviation) parameter.
/// </summary>
public double Sigma
{
Expand Down Expand Up @@ -163,9 +169,25 @@ double sigma
{
return sigma > 0.0;
}

/// <summary>
/// Estimate and set all distribution parameters based on a sample set.
/// </summary>
/// <param name="samples">Samples of this distribution.</param>
public
void
EstimateDistributionParameters(
IEnumerable<double> samples
)
{
Accumulator accumulator = new Accumulator(samples);
SetDistributionParameters(accumulator.Mean, accumulator.Sigma);
}

#endregion

#region Distribution Properties

/// <summary>
/// Gets the minimum possible value of generated random numbers.
/// </summary>
Expand Down Expand Up @@ -251,9 +273,11 @@ double x
{
return _sigma * Constants.Sqrt2 * Fn.ErfInverse(2.0 * x - 1.0) + _mu;
}

#endregion

#region Generator

/// <summary>
/// Returns a normal/gaussian distributed floating point random number.
/// </summary>
Expand All @@ -264,6 +288,7 @@ public override
{
return _mu + _standard.NextDouble() * _sigma;
}

#endregion
}
}
2 changes: 1 addition & 1 deletion src/app/MathNet.Iridium/Library/Distributions/Continuous/ParetoDistribution.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ public static
double shape
)
{
return location > 0.0 && shape > 0.0; ;
return location > 0.0 && shape > 0.0;
}
#endregion

Expand Down
3 changes: 2 additions & 1 deletion src/app/MathNet.Iridium/Library/Distributions/Continuous/StandardDistribution.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,9 +182,10 @@ public override
_extraNormal = null;
return extraNormalCpy;
}
// Generating two new gaussian deviates
else
{
// Generating two new gaussian deviates

double fac, rsq, v1, v2;

// We need a non-zero random point inside the unit circle.
Expand Down
32 changes: 30 additions & 2 deletions src/app/MathNet.Iridium/Library/LinearAlgebra/LUDecomposition.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@ public class LUDecomposition
OnDemandComputation<Matrix> _upperTriangularFactorOnDemand;
OnDemandComputation<int[]> _pivotOnDemand;
OnDemandComputation<Vector> _pivotVectorOnDemand;
OnDemandComputation<Matrix> _permutationMatrixOnDemand;
OnDemandComputation<double> _determinantOnDemand;

/// <summary>
Expand Down Expand Up @@ -140,9 +141,15 @@ Matrix A
{
for(int k = 0; k < _columnCount; k++)
{
double t = LU[p][k]; LU[p][k] = LU[j][k]; LU[j][k] = t;
double t = LU[p][k];
LU[p][k] = LU[j][k];
LU[j][k] = t;
}
int k2 = piv[p]; piv[p] = piv[j]; piv[j] = k2;

int k2 = piv[p];
piv[p] = piv[j];
piv[j] = k2;

pivsign = -pivsign;
}

Expand Down Expand Up @@ -201,6 +208,14 @@ public Vector PivotVector
get { return _pivotVectorOnDemand.Compute(); }
}

/// <summary>
/// Returns the permutation matrix P, such that L*U = P*X.
/// </summary>
public Matrix PermutationMatrix
{
get { return _permutationMatrixOnDemand.Compute(); }
}

/// <summary>
/// Returns pivot permutation vector as a one-dimensional double array.
/// </summary>
Expand Down Expand Up @@ -283,6 +298,7 @@ Matrix B
_upperTriangularFactorOnDemand = new OnDemandComputation<Matrix>(ComputeUpperTriangularFactor);
_pivotOnDemand = new OnDemandComputation<int[]>(ComputePivot);
_pivotVectorOnDemand = new OnDemandComputation<Vector>(ComputePivotVector);
_permutationMatrixOnDemand = new OnDemandComputation<Matrix>(ComputePermutationMatrix);
_determinantOnDemand = new OnDemandComputation<double>(ComputeDeterminant);
}

Expand Down Expand Up @@ -368,6 +384,18 @@ Matrix B
return new Vector(vals);
}

Matrix
ComputePermutationMatrix()
{
int[] pivot = Pivot;
double[][] perm = Matrix.CreateMatrixData(pivot.Length, pivot.Length);
for(int i = 0; i < pivot.Length; i++)
{
perm[pivot[i]][i] = 1.0;
}
return new Matrix(perm);
}

double
ComputeDeterminant()
{
Expand Down
Loading

0 comments on commit 6c86800

Please sign in to comment.