Skip to content

Commit ea750f9

Browse files
jiahaojiahao
committed
Fixes Haar symbolic integrator to work out of scope
1 parent 44010cf commit ea750f9

File tree

3 files changed

+351
-35
lines changed

3 files changed

+351
-35
lines changed

‎src/Haar.jl

Lines changed: 290 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,290 @@
1+
using GSL
2+
using Catalan
3+
4+
export permutations_in_Sn, compose, cycle_structure, data, part #Functions working with partitions and permutations
5+
UniformHaar, expectation, WeingartenUnitary
6+
7+
#Functions working with partitions and permutations
8+
# TODO Migrate these functions to Catalan
9+
10+
#Iterate over partitions of n in lexicographic order
11+
function part(n::Integer)
12+
if n==0 produce([]) end
13+
if n<=0 return end
14+
for p in @task part(n-1)
15+
p = [p; 1]
16+
produce(p)
17+
p = p[1:end-1]
18+
if length(p) == 1 || (length(p)>1 && p[end]<p[end-1])
19+
p[end] += 1
20+
produce(p)
21+
p[end] -= 1
22+
end
23+
end
24+
end
25+
26+
function permutations_in_Sn(n::Integer)
27+
P = permutation_calloc(n)
28+
while true
29+
produce(P)
30+
try permutation_next(P) catch break end
31+
end
32+
end
33+
34+
function compose(P::Ptr{gsl_permutation}, Q::Ptr{gsl_permutation})
35+
#Compose the permutations
36+
n=convert(Int64, permutation_size(P))
37+
@assert n==permutation_size(Q)
38+
x=permutation_alloc(n)
39+
Pv = [x+1 for x in pointer_to_array(permutation_data(P), (n,))]
40+
Qv = [x+1 for x in pointer_to_array(permutation_data(Q), (n,))]
41+
Xv = [Qv[i] for i in Pv]
42+
for i=1:n
43+
unsafe_assign(permutation_data(x), Xv[i]-1, i)
44+
end
45+
permutation_valid(x)
46+
x
47+
end
48+
49+
#Compute cycle structure, i.e. lengths of each cycle in the cycle factorization of the permutatation
50+
function cycle_structure(P::Ptr{gsl_permutation})
51+
n=convert(Int64, permutation_size(P))
52+
Pcanon = permutation_linear_to_canonical(P)
53+
PCANON = data(Pcanon)
54+
HaveNewCycleAtPos = [PCANON[i]>PCANON[i+1] for i=1:n-1]
55+
cycleindices = Int[1]
56+
for i=1:n-1 if HaveNewCycleAtPos[i] push!(cycleindices, i+1) end end
57+
push!(cycleindices, n+1)
58+
cyclestructure = Int[cycleindices[i+1]-cycleindices[i] for i=1:length(cycleindices)-1]
59+
@assert sum(cyclestructure) == n
60+
cyclestructure
61+
end
62+
63+
#Returns a vector of indices (starting from 1 in the Julia convention)
64+
data(P::Ptr{gsl_permutation}) = [convert(Int64, x)+1 for x in
65+
pointer_to_array(permutation_data(P), (convert(Int64, permutation_size(P)) ,))]
66+
67+
68+
type UniformHaar <: AbstractMatrix
69+
beta::Real
70+
end
71+
72+
# In random matrix theory one often encounters expressions of the form
73+
#
74+
#X = Q * A * Q' * B
75+
#
76+
#where A and B are deterministic matrices with fixed numerical matrix entries
77+
#and Q is a random matrix that does not have explicitly defined matrix
78+
#elements. Instead, one takes an expectation over of expressions of this form
79+
#and this "integrates out" the Qs to produce a numeric result.
80+
#
81+
#expectation(X) #= some number
82+
#
83+
#Here is a function that symbolically calculates the expectation of a product
84+
#of matrices over the symmetric group that Q is uniform Haar over.
85+
#It takes an expression consisting of a product of matrices and replaces it
86+
#with an evaluated symbolic expression which is the expectation.
87+
function expectation(X::Expr)
88+
if X.head != :call
89+
error(string("Unexpected type of expression: ", X.head))
90+
end
91+
92+
n = length(X.args) - 1
93+
if n < 1 return eval(X) end #nothing to do Haar-wise
94+
95+
if X.args[1] != :*
96+
error("Unexpected expression, only products supported")
97+
end
98+
99+
# Parse expression involving products of matrices to extract the
100+
# positions of Haar matrices and their ctransposes
101+
Qidx=[] #Indices for Haar matrices
102+
Qpidx=[] #Indices for ctranspose of Haar matrices
103+
Others=[]
104+
MyQ=None
105+
for i=1:n
106+
thingy=X.args[i+1]
107+
if isa(thingy, Symbol)
108+
if isa(eval(thingy), UniformHaar)
109+
if MyQ==None MyQ=thingy end
110+
if MyQ == thingy
111+
Qidx=[Qidx; i]
112+
else
113+
warning("only one instance of UniformHaar supported, skipping the other guy ", thingy) end
114+
else
115+
Others = [Others; (thingy, i, i+1)]
116+
end
117+
println(i, ' ', thingy, "::", typeof(eval(thingy)))
118+
elseif isa(thingy, Expr)
119+
println(i, ' ', thingy, "::Expr")
120+
if thingy.head==symbol('\'') && length(thingy.args)>=1 #Maybe this is a Q'
121+
if isa(thingy.args[1], Symbol) && isa(eval(thingy.args[1]), UniformHaar)
122+
println("Here is a Qtranspose")
123+
Qpidx=[Qpidx; i]
124+
end
125+
end
126+
else
127+
error(string("Unexpected token ", thingy ," of type ", typeof(thingy)))
128+
end
129+
end
130+
if length(Qidx) == length(Qpidx) == 0 return eval(X) end #nothing to do Haar-wise
131+
println(MyQ, " is in places ", Qidx)
132+
println(MyQ, "' is in places ", Qpidx)
133+
134+
n = length(Qidx)
135+
#If there are different Qs and Q's, the answer is a big fat 0
136+
if n != length(Qpidx) return zeros(size(eval(X.args[2]),1)...) end
137+
138+
##################################
139+
# Step 2. Enumerate permutations #
140+
##################################
141+
AllTerms = {}
142+
for sigma in @task permutations_in_Sn(n)
143+
for tau in @task permutations_in_Sn(n)
144+
sigma_inv=permutation_inverse(sigma)
145+
#Compose the permutations
146+
perm=compose(sigma_inv, tau)
147+
148+
#Compute cycle structure, i.e. lengths of each cycle in the cycle
149+
#factorization of the permutatation
150+
cyclestruct = cycle_structure(perm)
151+
Qr = Int64[n+1 for n in Qidx]
152+
Qpr= Int64[Qpidx[n]+1 for n in data(tau)]
153+
#Consolidate deltas
154+
Deltas=Dict()
155+
for i=1:n
156+
Deltas[Qr[i]] = Qpidx[data(sigma)[i]]
157+
Deltas[Qidx[i]] = Qpr[data(tau)[i]]
158+
end
159+
#Print deltas
160+
print("V(", cyclestruct, ") ")
161+
for k in keys(Deltas)
162+
print("δ(",k,",",Deltas[k],") ")
163+
end
164+
for (Symb, col_idx, row_idx) in Others
165+
print(Symb,"(",col_idx,",",row_idx,") ")
166+
end
167+
println()
168+
print("= V(", cyclestruct, ") ")
169+
#Evaluate deltas over the indices of the other matrices
170+
ReindexedSymbols = []
171+
for (Symb, col_idx, row_idx) in Others
172+
new_col, new_row = col_idx, row_idx
173+
if has(Deltas, col_idx) new_col = Deltas[col_idx] end
174+
if has(Deltas, row_idx) new_row = Deltas[row_idx] end
175+
print(Symb,"(",new_col,",",new_row,") ")
176+
ReindexedSymbols = [ReindexedSymbols; (Symb, new_col, new_row)]
177+
end
178+
println()
179+
#Parse coefficient
180+
Coefficient= WeingartenUnitary(perm)
181+
#Reconstruct expression
182+
println("START PARSING")
183+
println("The term is =", ReindexedSymbols[end])
184+
Symb, left_idx, right_idx = pop!(ReindexedSymbols)
185+
Expression={{{Symb}, left_idx, right_idx}}
186+
while length(ReindexedSymbols) > 0
187+
pop_idx = expr_idx = do_transpose = is_left = nothing
188+
for expr_iter in enumerate(Expression)
189+
expr_idx, expr_string = expr_iter
190+
_, left_idx, right_idx = expr_string
191+
for iter in enumerate(ReindexedSymbols)
192+
idx, data = iter
193+
Symb, col_idx, row_idx = data
194+
if row_idx == left_idx
195+
pop_idx = idx
196+
do_transpose = false
197+
is_left = true
198+
println("Case A")
199+
break
200+
elseif col_idx == right_idx
201+
pop_idx = idx
202+
do_transpose = false
203+
is_left = false
204+
println("Case B")
205+
break
206+
elseif row_idx == right_idx
207+
pop_idx = idx
208+
do_transpose = true
209+
is_left = false
210+
println("Case C")
211+
break
212+
elseif col_idx == left_idx
213+
pop_idx = idx
214+
do_transpose = true
215+
is_left = true
216+
println("Case D")
217+
break
218+
end
219+
end
220+
if pop_idx != nothing break end
221+
end
222+
println("Terms left = ", ReindexedSymbols)
223+
if pop_idx == nothing #Found nothing, start new expression blob
224+
println("New word: The term is ", ReindexedSymbols[1])
225+
Symb, left_idx, right_idx = delete!(ReindexedSymbols, 1)
226+
push!(Expression, {{Symb}, left_idx, right_idx})
227+
else #Found something
228+
println("The term is =", ReindexedSymbols[pop_idx])
229+
Symb, col_idx, row_idx = delete!(ReindexedSymbols, pop_idx)
230+
Term = do_transpose ? Expr(symbol("'"), Symb) : Symb
231+
if is_left
232+
insert!(Expression[expr_idx][1], 1, Term)
233+
Expression[expr_idx][2] = do_transpose ? row_idx : col_idx
234+
else
235+
push!(Expression[expr_idx][1], Term)
236+
Expression[expr_idx][3] = do_transpose ? col_idx : row_idx
237+
end
238+
end
239+
end
240+
println("DONE PARSING TREE")
241+
242+
# Evaluate closed cycles
243+
NewExpression={}
244+
for ExprChunk in Expression
245+
ExprBlob, start_idx, end_idx = ExprChunk
246+
print(ExprChunk, " => ")
247+
if start_idx == end_idx #Have a cycle; this is a trace
248+
ex= length(ExprBlob)==1 ? :(trace($(ExprBlob[1]))) :
249+
Expr(:call, :trace, Expr(:call, :*, ExprBlob...))
250+
else #Not a cycle, regular chain of multiplications
251+
ex= length(ExprBlob)==1 ? ExprBlob[1] :
252+
Expr(:call, :*, ExprBlob...)
253+
end
254+
push!(NewExpression, ex)
255+
println(ex)
256+
end
257+
ex = Expr(:call, :*, Coefficient, NewExpression...)
258+
println("Final expression is: ", ex)
259+
push!(AllTerms, ex)
260+
end
261+
end
262+
X = length(AllTerms)==1 ? AllTerms[1] : Expr(:call, :+, AllTerms...)
263+
println("THE ANSWER IS ", X)
264+
eval(X)
265+
end
266+
267+
#Computes the Weingarten function for permutations
268+
function WeingartenUnitary(P::Ptr{gsl_permutation})
269+
C = cycle_structure(P)
270+
WeingartenUnitary(C)
271+
end
272+
#Computes the Weingarten function for partitions
273+
function WeingartenUnitary(P::partition)
274+
n = sum(P)
275+
m = length(P)
276+
thesum = 0.0
277+
for irrep in @task part(n)
278+
#Character of the partition
279+
S = character(irrep, P)
280+
#Character of the identity divided by n!
281+
T = character_identity(irrep)
282+
#Denominator f_r(N) of (2.10)
283+
f = prod([factorial(BigInt(N + P[i] - i)) / factorial(BigInt(N - i))
284+
for i=1:m])
285+
thesum += S*T/(factorial(n)*f)
286+
println(irrep, " ", thesum)
287+
end
288+
thesum
289+
end
290+

0 commit comments

Comments
 (0)