diff --git a/src/Ginibre.jl b/src/Ginibre.jl index 21f9df2..80767b0 100644 --- a/src/Ginibre.jl +++ b/src/Ginibre.jl @@ -1,5 +1,6 @@ export rand, Ginibre import Base.rand +import Random: AbstractRNG, GLOBAL_RNG """ Ginibre(β::Int, N::Int) <: ContinuousMatrixDistribution @@ -30,30 +31,38 @@ struct Ginibre <: ContinuousMatrixDistribution end """ - rand(W::Ginibre) + rand(rng::AbstractRNG, W::Ginibre) Samples a matrix from the Ginibre ensemble. For `β = 1,2,4`, generates matrices randomly sampled from the real, complex, and quaternion Ginibre ensemble, respectively. """ -function rand(W::Ginibre) +function rand(rng::AbstractRNG, W::Ginibre) beta, n = W.beta, W.N if beta==1 - randn(n,n) + randn(rng,n,n) elseif beta==2 - randn(n,n)+im*randn(n,n) + randn(rng,n,n)+im*randn(rng,n,n) elseif beta==4 - Q0=randn(n,n) - Q1=randn(n,n) - Q2=randn(n,n) - Q3=randn(n,n) + Q0=randn(rng,n,n) + Q1=randn(rng,n,n) + Q2=randn(rng,n,n) + Q3=randn(rng,n,n) [Q0+im*Q1 Q2+im*Q3;-Q2+im*Q3 Q0-im*Q1] else error(string("beta = ", beta, " not implemented")) end end + +rand(W::Ginibre)=rand(GLOBAL_RNG, W) + + + + + + function jpdf(Z::AbstractMatrix{z}) where {z<:Complex} pi^(size(Z,1)^2)*exp(-trace(Z'*Z)) end diff --git a/src/HaarMeasure.jl b/src/HaarMeasure.jl index 417302c..644b1f7 100644 --- a/src/HaarMeasure.jl +++ b/src/HaarMeasure.jl @@ -1,6 +1,6 @@ #TODO implement O(n^2) method """ - rand(W::Haar, n::Int) + rand(rng::AbstractRNG, W::Haar, n::Int) Computes samples of real or complex Haar matrices of size `n`×`n`. @@ -24,19 +24,19 @@ implemented in most versions of LAPACK. - Edelman and Rao, 2005 - Mezzadri, 2006, math-ph/0609050 """ -function rand(W::Haar, n::Int, doCorrection::Int=1) +function rand(rng::AbstractRNG, W::Haar, n::Int, doCorrection::Int=1) beta = W.beta - M=rand(Ginibre(beta,n)) + M=rand(rng,Ginibre(beta,n)) q,r=qr(M) if doCorrection==0 q elseif doCorrection==1 if beta==1 - L = sign.(rand(n).-0.5) + L = sign.(rand(rng,n).-0.5) elseif beta==2 - L = exp.(im*rand(n)*2pi) + L = exp.(im*rand(rng,n)*2pi) elseif beta==4 - L = exp.(im*rand(2n)*2pi) + L = exp.(im*rand(rng,2n)*2pi) else error(string("beta = ",beta, " not implemented.")) end @@ -54,6 +54,10 @@ function rand(W::Haar, n::Int, doCorrection::Int=1) end end +rand(W::Haar, n::Int, doCorrection::Int=1)= rand(GLOBAL_RNG,W,n,doCorrection) + + + #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 diff --git a/test/Haar.jl b/test/Haar.jl index 396ceac..ed3c4a4 100644 --- a/test/Haar.jl +++ b/test/Haar.jl @@ -1,6 +1,7 @@ using RandomMatrices using LinearAlgebra: I, tr, Diagonal, QRPackedQ using Test +using Random @testset "Haar" begin @@ -29,4 +30,23 @@ for T in (Float64, ComplexF64) @test C*A' ≈ C*A2' end + +@testset "Local RNG reproducibility" begin + rng1 = MersenneTwister(1234) + rng2 = MersenneTwister(1234) + + A = rand(rng1, Haar(1), N) + B = rand(rng2, Haar(1), N) + + @test A ≈ B + + # Confirm that rand still works without passing an RNG + C = rand(Haar(1), N) + @test size(C) == (N, N) +end + + end # testset + + +