Add docstrings for random matrix ensembles (#89)

* add /test Project.toml

* add doctests.jl test file and /docs folder

* add doctests for random matrix ensembles

* switch jldoctest to @example

Author: apkille

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"]

docs/.gitignore

@@ -0,0 +1,2 @@
+build/
+site/

docs/Project.toml

@@ -0,0 +1,3 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+RandomMatrices = "2576dda1-a324-5b11-aa66-c48ed7e3c618"

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

docs/src/API.md

@@ -0,0 +1,17 @@
+# [Full API](@id Full-API)
+
+```@raw html
+<style>
+    .content table td {
+        padding-top: 0 !important;
+        padding-bottom: 0 !important;
+    }
+</style>
+```
+
+## Autogenerated API list
+
+```@autodocs
+Modules = [RandomMatrices]
+Private = false
+```

docs/src/index.md

@@ -0,0 +1,7 @@
+# RandomMatrices.jl
+
+```@meta
+DocTestSetup = quote
+    using RandomMatrices
+end
+```

src/GaussianEnsembles.jl

@@ -24,20 +24,21 @@ 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
-
-    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
+```@example
+julia> rand(Wigner(2), 3)
+3×3 LinearAlgebra.Hermitian{ComplexF64, Matrix{ComplexF64}}:
+   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
 GaussianHermite(β) = GaussianHermite{β}()
@@ -47,6 +48,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 +90,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 +157,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 +204,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 +258,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,42 @@
 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
+
+```@example
+julia> rand(Ginibre(2, 3))
+3×3 Matrix{ComplexF64}:
+ 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:
+- 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,30 @@ 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
+
+```@example
+julia> rand(Haar(2), 3)
+3×3 Matrix{ComplexF64}:
+ -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:
+- 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))