add doctests for random matrix ensembles

Author: apkille

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,15 +208,23 @@ 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))")
   end
   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,14 +262,35 @@ 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
   b::Real
 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)

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

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

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