From 9ea1328b648f3f25364d02583f841f1f32c42565 Mon Sep 17 00:00:00 2001 From: apkille Date: Thu, 22 May 2025 02:27:52 -0400 Subject: [PATCH 1/4] add /test Project.toml --- Project.toml | 7 ++++++- test/Project.toml | 5 +++++ 2 files changed, 11 insertions(+), 1 deletion(-) create mode 100644 test/Project.toml diff --git a/Project.toml b/Project.toml index f81807d..40feea5 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,6 @@ GSL = "92c85e6c-cbff-5e0c-80f7-495c94daaecd" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Combinatorics = "0.7.0, 1" @@ -19,3 +18,9 @@ FastGaussQuadrature = "0.3, 0.4, 0.5, 1" GSL = "0.4, 1" SpecialFunctions = "0.7, 1, 2" julia = "1.6" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..3cdade3 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,5 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From f3901a34318625a4506154778e1ed0d41598daf6 Mon Sep 17 00:00:00 2001 From: apkille Date: Thu, 22 May 2025 02:28:37 -0400 Subject: [PATCH 2/4] add doctests.jl test file and /docs folder --- docs/.gitignore | 2 ++ docs/Project.toml | 3 +++ docs/make.jl | 28 ++++++++++++++++++++++++++++ docs/src/API.md | 17 +++++++++++++++++ docs/src/index.md | 7 +++++++ test/doctests.jl | 8 ++++++++ test/runtests.jl | 1 + 7 files changed, 66 insertions(+) create mode 100644 docs/.gitignore create mode 100644 docs/Project.toml create mode 100644 docs/make.jl create mode 100644 docs/src/API.md create mode 100644 docs/src/index.md create mode 100644 test/doctests.jl diff --git a/docs/.gitignore b/docs/.gitignore new file mode 100644 index 0000000..da3d337 --- /dev/null +++ b/docs/.gitignore @@ -0,0 +1,2 @@ +build/ +site/ \ No newline at end of file diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..9c892e9 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,3 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +RandomMatrices = "2576dda1-a324-5b11-aa66-c48ed7e3c618" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..f74f42c --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,28 @@ +push!(LOAD_PATH,"../src/") + +using Documenter +using RandomMatrices + +DocMeta.setdocmeta!(RandomMatrices, :DocTestSetup, :(using RandomMatrices, Random); recursive=true) + +function main() + + makedocs( + doctest = false, + clean = true, + sitename = "RandomMatrices.jl", + format = Documenter.HTML( + assets=["assets/init.js"] + ), + modules = [RandomMatrices], + checkdocs = :exports, + warnonly = false, + authors = "Andrew Kille", + pages = [ + "RandomMatrices.jl" => "index.md", + "Full API" => "API.md" + ] + ) +end + +main() \ No newline at end of file diff --git a/docs/src/API.md b/docs/src/API.md new file mode 100644 index 0000000..f254186 --- /dev/null +++ b/docs/src/API.md @@ -0,0 +1,17 @@ +# [Full API](@id Full-API) + +```@raw html + +``` + +## Autogenerated API list + +```@autodocs +Modules = [RandomMatrices] +Private = false +``` \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..1ab30c2 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,7 @@ +# RandomMatrices.jl + +```@meta +DocTestSetup = quote + using RandomMatrices +end +``` \ No newline at end of file diff --git a/test/doctests.jl b/test/doctests.jl new file mode 100644 index 0000000..82c6ec0 --- /dev/null +++ b/test/doctests.jl @@ -0,0 +1,8 @@ +using Documenter +using RandomMatrices +using Test + +@testset "Doctests" begin + DocMeta.setdocmeta!(RandomMatrices, :DocTestSetup, :(using RandomMatrices, Random); recursive=true) + doctest(RandomMatrices) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 923f62a..0212a0f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ using Test include("FormalPowerSeries.jl") include("Haar.jl") include("StochasticProcess.jl") + include("doctests.jl") @testset "densities" begin include("densities/Semicircle.jl") From fa4514e93d5a5e89d472a6917395ad6b35484ab4 Mon Sep 17 00:00:00 2001 From: apkille Date: Thu, 22 May 2025 02:28:59 -0400 Subject: [PATCH 3/4] add doctests for random matrix ensembles --- src/GaussianEnsembles.jl | 101 ++++++++++++++++++++++++++++++++------- src/Ginibre.jl | 37 ++++++++++++-- src/Haar.jl | 25 ++++++++++ src/HaarMeasure.jl | 46 ++++++++++-------- 4 files changed, 166 insertions(+), 43 deletions(-) diff --git a/src/GaussianEnsembles.jl b/src/GaussianEnsembles.jl index ee6ecb0..feb6adb 100644 --- a/src/GaussianEnsembles.jl +++ b/src/GaussianEnsembles.jl @@ -24,20 +24,25 @@ export GaussianHermite, GaussianLaguerre, GaussianJacobi, ##################### """ -GaussianHermite{β} represents a Gaussian Hermite ensemble with parameter β. + GaussianHermite(β::Int) <: ContinuousMatrixDistribution -Wigner{β} is a synonym. +Represents a Gaussian Hermite ensemble with Dyson index `β`. -Example of usage: +`Wigner(β)` is a synonym. - β = 2 #β = 1, 2, 4 generates real, complex and quaternionic matrices respectively. - d = Wigner{β} #same as GaussianHermite{β} +## Examples - n = 20 #Generate square matrices of this size +```jldoctest +julia> d = Wigner(2); n = 3; - S = rand(d, n) #Generate a 20x20 symmetric Wigner matrix - T = tridrand(d, n) #Generate the symmetric tridiagonal form - v = eigvalrand(d, n) #Generate a sample of eigenvalues +julia> Random.seed!(1234); + +julia> rand(d, n) +3×3 LinearAlgebra.Hermitian{ComplexF64, Matrix{ComplexF64}}: + 0.560409+0.0im -0.292145+0.125521im 1.04191+0.513798im + -0.292145-0.125521im -0.346868+0.0im 0.0228835-0.00164725im + 1.04191-0.513798im 0.0228835+0.00164725im -0.118721+0.0im +``` """ struct GaussianHermite{β} <: ContinuousMatrixDistribution end GaussianHermite(β) = GaussianHermite{β}() @@ -47,6 +52,13 @@ Synonym for GaussianHermite{β} """ const Wigner{β} = GaussianHermite{β} +""" + rand(d::Wigner, n::Int) + +Generates an `n × n` matrix randomly sampled from the Gaussian-Hermite ensemble (also known as the Wigner ensemble). + +The Dyson index `β` is restricted to `β = 1,2` or `4`, for real, complex, and quaternionic fields, respectively. +""" rand(d::Type{Wigner{β}}, dims...) where {β} = rand(d(), dims...) function rand(d::Wigner{1}, n::Int) @@ -82,12 +94,15 @@ function rand(d::Wigner{β}, dims::Int...) where {β} end """ -Generates a nxn symmetric tridiagonal Wigner matrix + tridand(d::Wigner, n::Int) + +Generates an `n × n` symmetric tridiagonal matrix from the Gaussian-Hermite ensemble (also known as the Wigner ensemble). -Unlike for `rand(Wigner{β}, n)`, which is restricted to β=1,2 or 4, -`trirand(Wigner{β}, n)` will generate a +Unlike for `rand(Wigner(β), n)`, which is restricted to `β = 1,2` or `4`, +the call `trirand(Wigner(β), n)` will generate a tridiagonal Wigner matrix for any positive +value of `β`, including infinity. -The β=∞ case is defined in Edelman, Persson and Sutton, 2012 +The `β == ∞` case is defined in Edelman, Persson and Sutton, 2012. """ function tridrand(d::Wigner{β}, n::Int) where {β} χ(df::Real) = rand(Distributions.Chi(df)) @@ -146,16 +161,37 @@ end # Laguerre ensemble # ##################### +""" + GaussianLaguerre(β::Real, a::Real)` <: ContinuousMatrixDistribution + +Represents a Gaussian-Laguerre ensemble with Dyson index `β` and `a` parameter +used to control the density of eigenvalues near `λ = 0`. + +`Wishart(β, a)` is a synonym. + +## Fields +- `beta`: Dyson index +- `a`: Parameter used for weighting the joint probability density function of the ensemble + +## References: +- Edelman and Rao, 2005 +""" mutable struct GaussianLaguerre <: ContinuousMatrixDistribution beta::Real a::Real end const Wishart = GaussianLaguerre - -#Generates a NxN Hermitian Wishart matrix #TODO Check - the eigenvalue distribution looks funky #TODO The appropriate matrix size should be calculated from a and one matrix dimension +""" + rand(d::GaussianLaguerre, dims::Tuple) + +Generate a random matrix sampled from the Gaussian Laguerre ensemble (also known as the Wishart ensemble) +with parameters defined in `d` and dimensions given by `dims`. + +The Dyson index `β` is restricted to `β = 1,2` or `4`, for real, complex, and quaternionic fields, respectively. +""" function rand(d::GaussianLaguerre, dims::Dim2) n = 2.0*a/d.beta if d.beta == 1 #real @@ -172,7 +208,11 @@ function rand(d::GaussianLaguerre, dims::Dim2) return (A * A') / dims[1] end -#Generates a NxN bidiagonal Wishart matrix +""" + bidrand(d::GaussianLaguerre, n::Int) + +Generate an `n × n` bidiagonal matrix sampled from the Gaussian Laguerre ensemble (also known as the Wishart ensemble). +""" function bidrand(d::GaussianLaguerre, m::Integer) if d.a <= d.beta*(m-1)/2.0 error("Given your choice of m and beta, a must be at least $(d.beta*(m-1)/2.0) (You said a = $(d.a))") @@ -180,7 +220,11 @@ function bidrand(d::GaussianLaguerre, m::Integer) Bidiagonal([chi(2*d.a-i*d.beta) for i=0:m-1], [chi(d.beta*i) for i=m-1:-1:1], true) end -#Generates a NxN tridiagonal Wishart matrix +""" + tridrand(d::GaussianLaguerre, n::Int) + +Generate an `n × n` tridiagonal matrix sampled from the Gaussian Laguerre ensemble (also known as the Wishart ensemble). +""" function tridrand(d::GaussianLaguerre, m::Integer) B = bidrand(d, m) L = B * B' @@ -218,7 +262,22 @@ end # Jacobi ensemble # ################### -#Generates a NxN self-dual MANOVA matrix +""" + GaussianJacobi(β::Real, a::Real, a::Real)` <: ContinuousMatrixDistribution + +Represents a Gaussian-Jacobi ensemble with Dyson index `β`, while +`a`and `b` are parameters used to weight the joint probability density function of the ensemble. + +`MANOVA(β, a, b)` is a synonym. + +## Fields +- `beta`: Dyson index +- `a`: Parameter used for shaping the joint probability density function near `λ = 0` +- `b`: Parameter used for shaping the joint probability density function near `λ = 1` + +## References: +- Edelman and Rao, 2005 +""" mutable struct GaussianJacobi <: ContinuousMatrixDistribution beta::Real a::Real @@ -226,6 +285,12 @@ mutable struct GaussianJacobi <: ContinuousMatrixDistribution end const MANOVA = GaussianJacobi +""" + rand(d::GaussianJacobi, n::Int) + +Generate an `n × n` random matrix sampled from the Gaussian-Jacobi ensemble (also known as the MANOVA ensemble) +with parameters defined in `d`. +""" function rand(d::GaussianJacobi, m::Integer) w1 = Wishart(m, int(2.0*d.a/d.beta), d.beta) w2 = Wishart(m, int(2.0*d.b/d.beta), d.beta) diff --git a/src/Ginibre.jl b/src/Ginibre.jl index c921396..008c917 100644 --- a/src/Ginibre.jl +++ b/src/Ginibre.jl @@ -1,15 +1,44 @@ export rand, Ginibre import Base.rand -#Samples a matrix from the Ginibre ensemble -#This ensemble lives in GL(N, F), the set of all invertible N x N matrices -#over the field F -#For beta=1,2,4, F=R, C, H respectively +""" + Ginibre(β::Int, N::Int) <: ContinuousMatrixDistribution + +Represents a Ginibre ensemble with Dyson index `β` living in `GL(N, F)`, the set +of all invertible `N × N` matrices over the field `F`. + +## Fields +- `beta`: Dyson index +- `N`: Matrix dimension over the field `F`. + +## Examples + +```jldoctest +julia> Random.seed!(1234); + +julia> rand(Ginibre(2, 3)) +3×3 Matrix{ComplexF64}: + 0.970656-0.763689im -0.0328031-0.0998909im 2.70742+0.942733im + -0.979218-0.534709im -0.600792-0.726142im 1.52445-0.00991272im + 0.901861-0.837116im -1.44518-0.00420647im -0.20563-0.66748im +``` + +## References: +- Edelman and Rao, 2005 +""" struct Ginibre <: ContinuousMatrixDistribution beta::Float64 N::Integer end +""" + rand(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) beta, n = W.beta, W.N if beta==1 diff --git a/src/Haar.jl b/src/Haar.jl index 7c67e8e..3d0761b 100644 --- a/src/Haar.jl +++ b/src/Haar.jl @@ -64,7 +64,32 @@ end data(P::Ptr{gsl_permutation}) = [convert(Int64, x)+1 for x in pointer_to_array(permutation_data(P), (convert(Int64, permutation_size(P)) ,))] +""" + Haar(β::Int) <: ContinuousMatrixDistribution +Represents a Haar measure with Dyson index `β`, in which values of `β = 1,2` or `4` +correspond to 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. + +## Fields +- `beta`: Dyson index + +## Examples + +```jldoctest +julia> Random.seed!(1234); + +julia> rand(Haar(2), 3) +3×3 Matrix{ComplexF64}: + -0.587093-0.106607im -0.109587-0.0619909im 0.666498+0.428819im + 0.120119+0.525468im 0.205226+0.641757im -0.0408499+0.503801im + -0.591627-0.0582071im 0.725559+0.0611728im -0.28036-0.194448im +``` + +## References: +- Edelman and Rao, 2005 +""" mutable struct Haar <: ContinuousMatrixDistribution beta::Real end diff --git a/src/HaarMeasure.jl b/src/HaarMeasure.jl index f4a81d5..417302c 100644 --- a/src/HaarMeasure.jl +++ b/src/HaarMeasure.jl @@ -1,25 +1,29 @@ -# 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. -# 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. -# -# The last parameter specifies whether or not the piecewise correction -# is applied to ensure that it truly of Haar measure -# This addresses an inconsistency in the Householder reflections as -# implemented in most versions of LAPACK -# Method 0: No correction -# Method 1: Multiply rows by uniform random phases -# Method 2: Multiply rows by phases of diag(R) -# References: -# Edelman and Rao, 2005 -# Mezzadri, 2006, math-ph/0609050 #TODO implement O(n^2) method -#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::Int) + +Computes samples of real or complex Haar matrices of size `n`×`n`. + +For `beta = 1,2,4`, generates random orthogonal, unitary and symplectic matrices +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. + + rand(W::Haar, n::Int, doCorrection::Int = 1) + +The additional argument `doCorrection` specifies whether or not the piecewise correction +is applied to ensure that it is truly of Haar measure. +This addresses an inconsistency in the Householder reflections as +implemented in most versions of LAPACK. +- Method 0: No correction +- Method 1: Multiply rows by uniform random phases +- Method 2: Multiply rows by phases of diag(R) + +## References: +- Edelman and Rao, 2005 +- Mezzadri, 2006, math-ph/0609050 +""" function rand(W::Haar, n::Int, doCorrection::Int=1) beta = W.beta M=rand(Ginibre(beta,n)) From d2b3451867c70eff824eac040984c5f4550d3ed5 Mon Sep 17 00:00:00 2001 From: apkille Date: Thu, 22 May 2025 03:41:57 -0400 Subject: [PATCH 4/4] switch jldoctest to @example --- src/GaussianEnsembles.jl | 14 +++++--------- src/Ginibre.jl | 10 ++++------ src/Haar.jl | 10 ++++------ 3 files changed, 13 insertions(+), 21 deletions(-) diff --git a/src/GaussianEnsembles.jl b/src/GaussianEnsembles.jl index feb6adb..ddb39f7 100644 --- a/src/GaussianEnsembles.jl +++ b/src/GaussianEnsembles.jl @@ -32,16 +32,12 @@ Represents a Gaussian Hermite ensemble with Dyson index `β`. ## Examples -```jldoctest -julia> d = Wigner(2); n = 3; - -julia> Random.seed!(1234); - -julia> rand(d, n) +```@example +julia> rand(Wigner(2), 3) 3×3 LinearAlgebra.Hermitian{ComplexF64, Matrix{ComplexF64}}: - 0.560409+0.0im -0.292145+0.125521im 1.04191+0.513798im - -0.292145-0.125521im -0.346868+0.0im 0.0228835-0.00164725im - 1.04191-0.513798im 0.0228835+0.00164725im -0.118721+0.0im + 0.383322+0.0im -0.0452508+0.491032im -0.313208-0.330435im + -0.0452508-0.491032im -0.264521+0.0im -0.131337+0.0904235im + -0.313208+0.330435im -0.131337-0.0904235im -0.481758+0.0im ``` """ struct GaussianHermite{β} <: ContinuousMatrixDistribution end diff --git a/src/Ginibre.jl b/src/Ginibre.jl index 008c917..21f9df2 100644 --- a/src/Ginibre.jl +++ b/src/Ginibre.jl @@ -13,14 +13,12 @@ of all invertible `N × N` matrices over the field `F`. ## Examples -```jldoctest -julia> Random.seed!(1234); - +```@example julia> rand(Ginibre(2, 3)) 3×3 Matrix{ComplexF64}: - 0.970656-0.763689im -0.0328031-0.0998909im 2.70742+0.942733im - -0.979218-0.534709im -0.600792-0.726142im 1.52445-0.00991272im - 0.901861-0.837116im -1.44518-0.00420647im -0.20563-0.66748im + 0.781329+2.00346im 0.0595122+0.488652im -0.323494-0.35966im + 1.11089+0.935174im -0.384457+1.71419im 0.114358-0.360676im + 1.54119+0.362003im -0.693623-2.50141im -1.42383-1.06341im ``` ## References: diff --git a/src/Haar.jl b/src/Haar.jl index 3d0761b..03a8459 100644 --- a/src/Haar.jl +++ b/src/Haar.jl @@ -77,14 +77,12 @@ classical orthogonal, unitary and symplectic groups `O(n)`, `U(n)` and ## Examples -```jldoctest -julia> Random.seed!(1234); - +```@example julia> rand(Haar(2), 3) 3×3 Matrix{ComplexF64}: - -0.587093-0.106607im -0.109587-0.0619909im 0.666498+0.428819im - 0.120119+0.525468im 0.205226+0.641757im -0.0408499+0.503801im - -0.591627-0.0582071im 0.725559+0.0611728im -0.28036-0.194448im + -0.275126-0.112754im -0.217139-0.293544im 0.299633-0.829756im + 0.48675-0.575106im 0.226526-0.445825im -0.406164-0.131472im + -0.245835-0.532433im 0.375591+0.689594im -0.0243468-0.197175im ``` ## References: