Modernize code and tests

Author: jiahao

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

REQUIRE

@@ -1,6 +1,5 @@
-julia 0.3+
-Compat
-Combinatorics
+julia 0.4
+Combinatorics 0.2
 Docile
 Distributions
 ODE

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
 
 

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
 

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

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
 

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)
+

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'))))