Fixes exports and abstracts Haar random matrix types

Author: jiahao

src/HaarMeasure.jl

@@ -1,8 +1,11 @@
+export rand, Ginibre, NeedPiecewiseCorrection, jpdfGinibre,
+    CircularOrthogonal, CircularUnitary, CircularSymplectic
+import Base.rand
+
 # Computes samplex of real or complex Haar matrices of size nxn
 #
-# For beta=1,2,4, generates random matrices from the circular
-# orthogonal, unitary and symplectic ensembles (COE, CUE, CSE) respectively
-# These are collectively known as the Dyson circular ensembles.
+# 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.
@@ -12,14 +15,15 @@
 # This addresses an inconsistency in the Householder reflections as
 # implemented in most versions of LAPACK
 # Method 0: No correction
-# Method 1: Edelman correction - multiply rows by uniform random phases
-# Method 2: Mezzadri correction - multiply rows by phases of diag(R)
+# 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
-function HaarMatrix(n::Integer, beta::Integer, doCorrection::Integer)
-    M=GinibreMatrix(n, beta)
+function rand(W::UniformHaar, doCorrection::Integer)
+    n, beta = W.N, W.beta
+    M=rand(Ginibre(n, beta))
     q,r=qr(M)
     if doCorrection==0
         q
@@ -49,16 +53,15 @@ 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
-HaarMatrix(n::Integer, beta::Integer) = HaarMatrix(n, beta, 1)
-
+#it doesn't matter, but better safe than sorry... let the user choose else
+rand(W::UniformHaar) = rand(W, 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
 function NeedPiecewiseCorrection()
     n=20
-    R=qr(randn(n,n))[2]
+    R=qr(randn(Ginibre(2,n)))[2]
     return any([x<0 for x in diag(R)])
 end
 
@@ -67,7 +70,12 @@ end
 #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
-function GinibreMatrix(n::Integer, beta::Integer)
+immutable Ginibre <: ContinuousMatrixDistribution
+   beta::Float64
+   N::Integer
+end
+
+function rand(W::Ginibre)
     if beta==1
         randn(n,n)
     elseif beta==2
@@ -87,6 +95,28 @@ function jpdfGinibre{z<:Complex}(Z::AbstractMatrix{z})
     pi^(size(Z,1)^2)*exp(-trace(Z'*Z))
 end
 
+#Dyson Circular orthogonal, unitary and symplectic ensembles
+immutable CircularOrthogonal
+    N :: Int64
+end
+
+function rand(M::CircularOrthogonal)
+    U = rand(Ginibre(2, M.N))
+    U * U'
+end
+
+immutable CircularUnitary
+    N :: Int64
+end
+
+rand(M::CircularUnitary) = rand(Ginibre(2, M.N))
+
+immutable CircularSymplectic
+    N :: Int64
+end
+
+rand(M::CircularSymplectic) = error("Not implemented")
+
 #TODO maybe, someday
 #Haar measure on U(N) in terms of local coordinates
 #Zyczkowski and Kus, Random unitary matrices, J. Phys. A: Math. Gen. 27,

src/Haar.jl renamed to src/HaarSymbolic.jl

@@ -1,6 +1,28 @@
 using GSL
 using Catalan
 
+export permutations_in_Sn, compose, cycle_structure, data, part #Functions working with partitions and permutations
+    UniformHaar, expectation, expectationtrace, WeingartenUnitary
+
+#Functions working with partitions and permutations
+# TODO Migrate these functions to Catalan
+
+#Iterate over partitions of n in lexicographic order
+function part(n::Integer)
+    if n==0 produce([]) end
+    if n<=0 return end
+    for p in @task part(n-1)
+        p = [p; 1]
+        produce(p)
+        p = p[1:end-1]
+        if length(p) == 1 || (length(p)>1 && p[end]<p[end-1])
+            p[end] += 1
+            produce(p)
+            p[end] -= 1
+        end
+    end
+end
+
 function permutations_in_Sn(n::Integer)
     P = permutation_calloc(n)
     while true 
@@ -42,7 +64,11 @@ end
 data(P::Ptr{gsl_permutation}) = [convert(Int64, x)+1 for x in
     pointer_to_array(permutation_data(P), (convert(Int64, permutation_size(P)) ,))]
 
-
+immutable UniformHaar <: ContinuousMatrixDistribution
+    beta::Float64
+    N::Integer
+end
+    
 # In random matrix theory one often encounters expressions of the form
 #
 #X = Q * A * Q' * B
@@ -58,7 +84,7 @@ data(P::Ptr{gsl_permutation}) = [convert(Int64, x)+1 for x in
 #of matrices over the symmetric group that Q is uniform Haar over.
 #It takes an expression consisting of a product of matrices and replaces it
 #with an evaluated symbolic expression which is the expectation.
-function Expectation(X::Expr)
+function expectation(X::Expr)
     if X.head != :call
         error(string("Unexpected type of expression: ", X.head))
     end
@@ -238,6 +264,8 @@ function Expectation(X::Expr)
     eval(X)
 end
 
+expectationtrace(X::Expr)=trace(expectation(X))
+
 #Computes the Weingarten function for permutations
 function WeingartenUnitary(P::Ptr{gsl_permutation})
     C = cycle_structure(P)
@@ -262,43 +290,3 @@ function WeingartenUnitary(P::partition)
     thesum
 end
 
-#Iterate over partitions of n in lexicographic order
-function part(n::Integer)
-    if n==0 produce([]) end
-    if n<=0 return end
-    for p in @task part(n-1)
-        p = [p; 1]
-        produce(p)
-        p = p[1:end-1]
-        if length(p) == 1 || (length(p)>1 && p[end]<p[end-1])
-            p[end] += 1
-            produce(p)
-            p[end] -= 1
-        end
-    end
-end
-
-###############
-# SANDBOX
-N=5
-A=randn(N,N)
-B=randn(N,N)
-type HaarMatrix
-    beta::Real
-end
-Q = HaarMatrix(2)
-
-println("Case 1")
-println("E(A*B) = ", Expectation(:(A*B)))
-println("A*B/N = ", A*B)
-
-println("Case 2")
-println("tr(A)*tr(B)/N = ", trace(A)*trace(B)/N)
-println(trace(B)*A)
-println("E(A*Q*B*Q') = ", Expectation(:(A*Q*B*Q')))
-println("E.tr(A*Q*B*Q') = ", trace(Expectation(:(A*Q*B*Q'))))
-
-println("Case 3")
-println(Expectation(:(A*B)) == A*B)
-println("E(A*Q*B*Q'*A*Q*B*Q') = ", Expectation(:(A*Q*B*Q'*A*Q*B*Q')))
-

src/RandomMatrices.jl

@@ -1,4 +1,6 @@
 require("Distributions")
+require("Catalan")
+require("GSL")
 
 module RandomMatrices
     importall Distributions
@@ -18,4 +20,8 @@ module RandomMatrices
 
     #Formal power series
     include("FormalPowerSeries.jl")
+    
+    #Symbolic integrator of uniform Haar matrices
+    include("HaarSymbolic.jl")
+    include("HaarMeasure.jl")
 end #module