From 357c48eb80aac2ee63f597fb70e577d84b4665b6 Mon Sep 17 00:00:00 2001 From: Jiahao Chen Date: Sat, 7 Nov 2015 18:07:18 -0500 Subject: [PATCH] Modernize code and tests --- LICENSE.md | 2 +- REQUIRE | 5 ++- src/FormalPowerSeries.jl | 42 +++++++++++------------ src/HaarMeasure.jl | 34 +++++++++++-------- src/InvariantEnsembles.jl | 9 ++--- src/RandomMatrices.jl | 10 ++---- test/FormalPowerSeries.jl | 71 +++++++++++++++++++++------------------ test/Haar.jl | 5 ++- 8 files changed, 93 insertions(+), 85 deletions(-) diff --git a/LICENSE.md b/LICENSE.md index 73e5fd9..19e9275 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,6 +1,6 @@ The RandomMatrices Julia module is licensed under the MIT License: -> Copyright (c) 2013-2014: Jiahao Chen, Alan Edelman, Jameson Nash, Sheehan +> Copyright (c) 2013-2015: Jiahao Chen, Alan Edelman, Jameson Nash, Sheehan > Olver and other contributors > > Permission is hereby granted, free of charge, to any person obtaining diff --git a/REQUIRE b/REQUIRE index 1a614ba..7df70da 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,6 +1,5 @@ -julia 0.3+ -Compat -Combinatorics +julia 0.4 +Combinatorics 0.2 Docile Distributions ODE diff --git a/src/FormalPowerSeries.jl b/src/FormalPowerSeries.jl index b2e1cf4..8fcd490 100644 --- a/src/FormalPowerSeries.jl +++ b/src/FormalPowerSeries.jl @@ -193,12 +193,12 @@ end #TODO implement the FFT-based algorithm of one of the following # doi:10.1109/TAC.1984.1103499 # https://cs.uwaterloo.ca/~glabahn/Papers/tamir-george-1.pdf -function reciprocal{T}(P::FormalPowerSeries{T}, n :: Integer) - n<0 ? error(sprintf("Invalid inverse truncation length %d requested", n)) : true +function reciprocal{T}(P::FormalPowerSeries{T}, n::Integer) + n<0 && error(sprintf("Invalid inverse truncation length %d requested", n)) a = tovector(P, 0:n-1) #Extract original coefficients in vector - a[1]==0 ? (error("Inverse does not exist")) : true - inv_a1 = T<:Number ? 1.0/a[1] : inv(a[1]) + a[1]==0 && error("Inverse does not exist") + inv_a1 = inv(a[1]) b = zeros(n) #Inverse b[1] = inv_a1 @@ -210,13 +210,13 @@ end #Derivative of fps [H. Sec.1.4, p.18] function derivative{T}(P::FormalPowerSeries{T}) - c = Dict{BigInt, T}() - for (k,v) in P.c - if k != 0 && v != 0 - c[k-1] = k*v + c = Dict{BigInt, T}() + for (k,v) in P.c + if k != 0 && v != 0 + c[k-1] = k*v + end end - end - FormalPowerSeries{T}(c) + FormalPowerSeries{T}(c) end #[H, Sec.1.4, p.18] @@ -233,17 +233,17 @@ end # Power of fps [H, Sec.1.5, p.35] function ^{T}(P::FormalPowerSeries{T}, n :: Integer) - if n == 0 - return FormalPowerSeries{T}([convert(T, 1)]) - elseif n > 1 - #The straightforward recursive way - return P^(n-1) * P - #TODO implement the non-recursive formula - elseif n==1 - return P - else - error(@sprintf("I don't know what to do for power n = %d", n)) - end + if n == 0 + return FormalPowerSeries{T}([one(T)]) + elseif n > 1 + #The straightforward recursive way + return P^(n-1) * P + #TODO implement the non-recursive formula + elseif n==1 + return P + else + error("I don't know what to do for power n = $n") + end end diff --git a/src/HaarMeasure.jl b/src/HaarMeasure.jl index 68e8eef..d19b344 100644 --- a/src/HaarMeasure.jl +++ b/src/HaarMeasure.jl @@ -1,10 +1,7 @@ -export NeedPiecewiseCorrection - - -# Computes samplex of real or complex Haar matrices of size nxn +# Computes samples of real or complex Haar matrices of size nxn # # For beta=1,2,4, generates random orthogonal, unitary and symplectic matrices -# of uniform Haar measure. +# of uniform Haar measure. # These matrices are distributed with uniform Haar measure over the # classical orthogonal, unitary and symplectic groups O(N), U(N) and # Sp(N)~USp(2N) respectively. @@ -20,13 +17,16 @@ export NeedPiecewiseCorrection # Edelman and Rao, 2005 # Mezzadri, 2006, math-ph/0609050 #TODO implement O(n^2) method -function rand(W::Haar, n::Integer, doCorrection::Integer) +#By default, always do piecewise correction +#For most applications where you use the HaarMatrix as a similarity transform +#it doesn't matter, but better safe than sorry... let the user choose else +function rand(W::Haar, n::Int, doCorrection::Int=1) beta = W.beta M=rand(Ginibre(beta,n)) q,r=qr(M) if doCorrection==0 q - elseif doCorrection==1 + elseif doCorrection==1 if beta==1 L = sign(rand(n).-0.5) elseif beta==2 @@ -50,11 +50,6 @@ function rand(W::Haar, n::Integer, doCorrection::Integer) end end -#By default, always do piecewise correction -#For most applications where you use the HaarMatrix as a similarity transform -#it doesn't matter, but better safe than sorry... let the user choose else -rand(W::Haar,n::Integer) = rand(W,n, 1) - #A utility method to check if the piecewise correction is needed #This checks the R part of the QR factorization; if correctly done, #the diagonals are all chi variables so are non-negative @@ -94,8 +89,9 @@ function Stewart(::Type{Float64}, n) for i = 0:n-2 larfg!(n - i, pβ + (n + 1)*8i, pH + (n + 1)*8i, 1, pτ + 8i) end - LinAlg.QRPackedQ(H,τ) + Base.LinAlg.QRPackedQ(H,τ) end + function Stewart(::Type{Complex128}, n) τ = Array(Complex128, n) H = complex(randn(n, n), randn(n,n)) @@ -107,7 +103,17 @@ function Stewart(::Type{Complex128}, n) for i = 0:n-2 larfg!(n - i, pβ + (n + 1)*16i, pH + (n + 1)*16i, 1, pτ + 16i) end - LinAlg.QRPackedQ(H,τ) + Base.LinAlg.QRPackedQ(H,τ) end +export randfast +function randfast(W::Haar, n::Int) + if W.beta==1 + Stewart(Float64, n) + elseif W.beta==2 + Stewart(Complex128, n) + else + error("beta = $beta not implemented") + end +end diff --git a/src/InvariantEnsembles.jl b/src/InvariantEnsembles.jl index 3713709..23dab73 100644 --- a/src/InvariantEnsembles.jl +++ b/src/InvariantEnsembles.jl @@ -139,7 +139,7 @@ function adaptiveie(w,μ0,α,β,d) end #Reads in recurrence relationship and constructs OPs -function InvariantEnsemble(str::@compat(AbstractString),V::Function,d,n::Integer) +function InvariantEnsemble(str::AbstractString,V::Function,d,n::Integer) file = Pkg.dir("RandomMatrices/data/CoefficientDatabase/" * str * "/" * string(n)) μ0=readdlm(file * "norm.csv")[1] A=readdlm(file * "rc.csv",',') @@ -151,7 +151,7 @@ end # For constructing InvariantEnsembles that do not scale with n -function InvariantEnsembleUnscaled(str::@compat(AbstractString),V::Function,d,n::Integer) +function InvariantEnsembleUnscaled(str::AbstractString,V::Function,d,n::Integer) file = Pkg.dir("RandomMatrices/data/CoefficientDatabaseUnscaled/" * str * "/") μ0=readdlm(file * "norm.csv")[1] A=readdlm(file * "rc.csv",',') @@ -166,7 +166,7 @@ end #Decides whether to use built in recurrence or read it in # Also contains ensemble data -function InvariantEnsemble(str::@compat(AbstractString),n::Integer) +function InvariantEnsemble(str::AbstractString,n::Integer) if(str == "GUE") InvariantEnsemble(str,x->x.^2,[-3.,3.],n) elseif(str == "Quartic") @@ -204,7 +204,8 @@ function iekernel(q::Array{Float64,2},d,plan::Function) end -samplespectra(str::@compat(AbstractString),n::Integer,m::Integer)=samplespectra(InvariantEnsemble(str,n),m) +<<<<<<< 16c414b28d81a04c589c8ec5e66aa724f4e110f0 +samplespectra(str::AbstractString,n::Integer,m::Integer)=samplespectra(InvariantEnsemble(str,n),m) # Sample eigenvalues of invariant ensemble, m times diff --git a/src/RandomMatrices.jl b/src/RandomMatrices.jl index 453740d..e3a4cb9 100644 --- a/src/RandomMatrices.jl +++ b/src/RandomMatrices.jl @@ -4,10 +4,6 @@ using Combinatorics,Compat import Distributions: ContinuousMatrixDistribution -if VERSION < v"0.4-" - using Docile -end - import Base.rand #If the GNU Scientific Library is present, turn on additional functionality. @@ -24,17 +20,15 @@ export bidrand, #Generate random bidiagonal matrix eigvalrand, #Generate random set of eigenvalues eigvaljpdf #Eigenvalue joint probability density -typealias Dim2 @compat(Tuple{Int,Int}) #Dimensions of a rectangular matrix +typealias Dim2 Tuple{Int, Int} #Dimensions of a rectangular matrix #Classical Gaussian matrix ensembles include("GaussianEnsembles.jl") - # Classical univariate distributions include("densities/Semicircle.jl") include("densities/TracyWidom.jl") - # Ginibre include("Ginibre.jl") @@ -58,7 +52,7 @@ include("StatisticalTests.jl") include("StochasticProcess.jl") #Invariant ensembles -if filemode(Pkg.dir("ApproxFun")) != 0 +if Pkg.installed("ApproxFun")!== nothing include("InvariantEnsembles.jl") end diff --git a/test/FormalPowerSeries.jl b/test/FormalPowerSeries.jl index 2bc9a3b..48603d4 100644 --- a/test/FormalPowerSeries.jl +++ b/test/FormalPowerSeries.jl @@ -1,69 +1,74 @@ +using RandomMatrices +using Base.Test + # Excursions in formal power series (fps) -MaxSeriesSize=10 -MaxRange = 50 -MatrixSize=15 +MaxSeriesSize=8 +MaxRange = 20 +MatrixSize=10 TT=Int64 -(n1, n2, n3) = round(Int,rand(3)*MaxSeriesSize) +(n1, n2, n3) = rand(1:MaxSeriesSize, 3) -X = FormalPowerSeries{TT}(round(Int,rand(n1)*MaxRange)) -Y = FormalPowerSeries{TT}(round(Int,rand(n2)*MaxRange)) -Z = FormalPowerSeries{TT}(round(Int,rand(n3)*MaxRange)) +X = FormalPowerSeries{TT}(rand(1:MaxRange, n1)) +Y = FormalPowerSeries{TT}(rand(1:MaxRange, n2)) +Z = FormalPowerSeries{TT}(rand(1:MaxRange, n3)) -c = round(Int,rand()*MatrixSize) #Size of matrix representation to generate +c = rand(1:MatrixSize) #Size of matrix representation to generate -nzeros = round(Int,rand()*MaxSeriesSize) -@assert X == trim(X) +nzeros = rand(1:MaxSeriesSize) +@test X == trim(X) XX = deepcopy(X) for i=1:nzeros - idx = round(Int,rand()*MaxRange) + idx = rand(1:MaxRange) if !haskey(XX.c, idx) XX.c[idx] = convert(TT, 0) end end -@assert trim(XX) == X +@test trim(XX) == X #Test addition, p.15, (1.3-4) -@assert X+X == 2X -@assert X+Y == Y+X -@assert MatrixForm(X+Y,c) == MatrixForm(X,c)+MatrixForm(Y,c) +@test X+X == 2X +@test X+Y == Y+X +@test MatrixForm(X+Y,c) == MatrixForm(X,c)+MatrixForm(Y,c) #Test subtraction, p.15, (1.3-4) -@assert X-X == 0X -@assert X-Y == -(Y-X) -@assert MatrixForm(X-Y,c) == MatrixForm(X,c)-MatrixForm(Y,c) +@test X-X == 0X +@test X-Y == -(Y-X) +@test MatrixForm(X-Y,c) == MatrixForm(X,c)-MatrixForm(Y,c) #Test multiplication, p.15, (1.3-5) -@assert X*Y == Y*X -@assert MatrixForm(X*X,c) == MatrixForm(X,c)*MatrixForm(X,c) -@assert MatrixForm(X*Y,c) == MatrixForm(X,c)*MatrixForm(Y,c) -@assert MatrixForm(Y*X,c) == MatrixForm(Y,c)*MatrixForm(X,c) -@assert MatrixForm(Y*Y,c) == MatrixForm(Y,c)*MatrixForm(Y,c) +@test X*Y == Y*X +@test MatrixForm(X*X,c) == MatrixForm(X,c)*MatrixForm(X,c) +@test MatrixForm(X*Y,c) == MatrixForm(X,c)*MatrixForm(Y,c) +@test MatrixForm(Y*X,c) == MatrixForm(Y,c)*MatrixForm(X,c) +@test MatrixForm(Y*Y,c) == MatrixForm(Y,c)*MatrixForm(Y,c) -@assert X.*Y == Y.*X +@test X.*Y == Y.*X #The reciprocal series has associated matrix that is the matrix inverse #of the original series #Force reciprocal to exist X.c[0] = 1 -discrepancy = (norm(inv(float(MatrixForm(X,c)))[1, :]'[:, 1] - tovector(reciprocal(X, c),0:c-1))) -if discrepancy > c*sqrt(eps()) - error(@sprintf("Error %f exceeds tolerance %f", discrepancy, c*sqrt(eps()))) +discrepancy = (norm(inv(float(MatrixForm(X,c)))[1, :]'[:, 1] - tovector(reciprocal(X, c), 0:c-1))) +tol = c*√eps() +if discrepancy > tol + error(string("Error ", discrepancy, " exceeds tolerance ", tol)) end -#@assert norm(inv(float(MatrixForm(X,c)))[1, :]'[:, 1] - tovector(reciprocal(X, c),c)) < c*sqrt(eps()) +#@test norm(inv(float(MatrixForm(X,c)))[1, :]'[:, 1] - tovector(reciprocal(X, c),c)) < c*sqrt(eps()) #Test differentiation XX = derivative(X) for (k, v) in XX.c - k==0 ? continue : true - @assert X.c[k+1] == v/(k+1) + k==0 && continue + @test X.c[k+1] == v/(k+1) end #Test product rule [H, Sec.1.4, p.19] -@assert derivative(X*Y) == derivative(X)*Y + X*derivative(Y) +@test derivative(X*Y) == derivative(X)*Y + X*derivative(Y) #Test right distributive law of composition [H, Sec.1.6, p.38] -@assert compose(X,Z)*compose(Y,Z) == compose(X*Y, Z) +@test compose(X,Z)*compose(Y,Z) == compose(X*Y, Z) #Test chain rule [H, Sec.1.6, p.40] -@assert derivative(compose(X,Y)) == compose(derivative(X),Y)*derivative(Y) +@test derivative(compose(X,Y)) == compose(derivative(X),Y)*derivative(Y) + diff --git a/test/Haar.jl b/test/Haar.jl index 39af17f..7f3c22b 100644 --- a/test/Haar.jl +++ b/test/Haar.jl @@ -1,3 +1,6 @@ +using RandomMatrices +using Base.Test + if isdefined(:_HAVE_GSL) N=5 A=randn(N,N) @@ -7,7 +10,7 @@ Q=UniformHaar(2, N) #Test case where there are no symbolic Haar matrices @test_approx_eq eval(expectation(N, :Q, :(A*B))) A*B #Test case where there is one pair of symbolic Haar matrices -@test_approx_eq eval(expectedtrace(N, :Q, :(A*Q*B*Q'))) trace(A)*trace(B)/N +@test_approx_eq eval(expectedtrace(N, :Q, :(A*Q*B*Q'))) trace(A)*trace(B)/N println("Case 3") println("E(A*Q*B*Q'*A*Q*B*Q') = ", eval(expectation(N, :Q, :(A*Q*B*Q'*A*Q*B*Q'))))