Skip to content

Commit bc8599e

Browse files
jiahaojiahao
committed
Prototype of symbolic Haar integrator.
Currently has a bug where the indices of the Haar matrices are assigned incorrectly.
1 parent c0456cf commit bc8599e

File tree

1 file changed

+293
-0
lines changed

1 file changed

+293
-0
lines changed

‎src/Haar.jl

Lines changed: 293 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,293 @@
1+
using GSL
2+
3+
function permutations_in_Sn(n::Integer)
4+
P = permutation_calloc(n)
5+
while true
6+
produce(P)
7+
try permutation_next(P) catch break end
8+
end
9+
end
10+
11+
function compose(P::Ptr{gsl_permutation}, Q::Ptr{gsl_permutation})
12+
#Compose the permutations
13+
n=convert(Int64, permutation_size(P))
14+
@assert n==permutation_size(Q)
15+
x=permutation_alloc(n)
16+
Pv = [x+1 for x in pointer_to_array(permutation_data(P), (n,))]
17+
Qv = [x+1 for x in pointer_to_array(permutation_data(Q), (n,))]
18+
Xv = [Qv[i] for i in Pv]
19+
for i=1:n
20+
unsafe_assign(permutation_data(x), Xv[i]-1, i)
21+
end
22+
permutation_valid(x)
23+
x
24+
end
25+
26+
#Compute cycle structure, i.e. lengths of each cycle in the cycle factorization of the permutatation
27+
function cycle_structure(P::Ptr{gsl_permutation})
28+
n=convert(Int64, permutation_size(P))
29+
Pcanon = permutation_linear_to_canonical(P)
30+
PCANON = data(Pcanon)
31+
HaveNewCycleAtPos = [PCANON[i]>PCANON[i+1] for i=1:n-1]
32+
cycleindices = [1]
33+
for i=1:n-1 if HaveNewCycleAtPos[i] cycleindices = [cycleindices; i+1] end end
34+
cycleindices = [cycleindices; n+1]
35+
cyclestructure = [cycleindices[i+1]-cycleindices[i] for i=1:length(cycleindices)-1]
36+
@assert sum(cyclestructure) == n
37+
cyclestructure
38+
end
39+
40+
#Returns a vector of indices (starting from 1 in the Julia convention)
41+
data(P::Ptr{gsl_permutation}) = [convert(Int64, x)+1 for x in pointer_to_array(permutation_data(P), (convert(Int64, permutation_size(P)) ,))]
42+
43+
44+
# In random matrix theory one often encounters expressions of the form
45+
#
46+
#X = trace(Q * A * Q' * B)
47+
#
48+
#where A and B are deterministic matrices with fixed numerical matrix entries and Q is a random matrix that does not have explicitly defined matrix elements. Instead, one takes an expectation over of expressions of this form and this "integrates out" the Qs to produce a numeric result.
49+
#
50+
#expectation(X) #= some number
51+
#
52+
#I'm curious how far one can implement, in Julia, a data type for Q and expectation() method for which expressions of this form can be written natively as code. This probably strays into the larger question of how much symbolic algebra can be written and handled in Julia, but I think this would be have real payoffs for making possible extremely elegant manipulations of random matrices that doesn't exist in any known computer language.
53+
#
54+
#Any thoughts on how feasible this is? Perhaps to begin with, it would be nice to have an expectation() that would work on arbitrary strings of products of matrices. Probably this would involve taking an input expression argument Would like to hear some stuff on t would be nice to start with just how one can write arbitrary products of matrices My initial thoughts are to declare Q::MyRandomMatrix with the type MyRandomMatrix<:AbstractMatrix and to write an expectation(Ex::Expr) function that at first blush traverses the expression Ex, finds a * that takes at least one MyRandomMatrix in its arguments and applies the appropriate substitution.
55+
#
56+
#So... 1) How feasible is this and 2) How do I traverse the expression in the way described in the previous paragraph? The "Metaprogramming" manual is somewhat sparse on this other than saying .args[], but that doesn't quite help
57+
58+
function Expectation(X::Expr)
59+
if X.head != :call
60+
error(string("Unexpected type of expression: ", X.head))
61+
end
62+
63+
n = length(X.args) - 1
64+
if n < 1 return eval(X) end #nothing to do Haar-wise
65+
66+
if X.args[1] != :*
67+
error("Unexpected expression, only products supported")
68+
end
69+
70+
# Parse expression involving products of matrices to extract the
71+
# positions of Haar matrices and their ctransposes
72+
Qidx=[] #Indices for Haar matrices
73+
Qpidx=[] #Indices for ctranspose of Haar matrices
74+
Others=[]
75+
MyQ=None
76+
for i=1:n
77+
thingy=X.args[i+1]
78+
if typeof(thingy)==Symbol
79+
matrixtype = typeof(eval(thingy))
80+
if matrixtype == HaarMatrix
81+
if MyQ==None MyQ=thingy end
82+
if MyQ == thingy
83+
Qidx=[Qidx; i]
84+
else
85+
println(string("only one instance of HaarMatrix supported, skipping the other guy ", thingy))
86+
end
87+
else
88+
Others = [Others; (thingy, i, i+1)]
89+
end
90+
println(i, ' ', thingy, "::", typeof(eval(thingy)))
91+
elseif typeof(thingy)==Expr
92+
println(i, ' ', thingy, "::Expr")
93+
if thingy.head==symbol('\'') && length(thingy.args)>=1 #Maybe this is a Q'
94+
if typeof(thingy.args[1])==Symbol && typeof(eval(thingy.args[1]))==HaarMatrix
95+
println("Here is a Qtranspose")
96+
Qpidx=[Qpidx; i]
97+
end
98+
end
99+
else
100+
error(string("Unexpected token ", thingy ," of type ", typeof(thingy)))
101+
end
102+
end
103+
104+
if length(Qidx) == length(Qpidx) == 0 return eval(X) end #nothing to do Haar-wise
105+
println(MyQ, " is in places ", Qidx)
106+
println(MyQ, "' is in places ", Qpidx)
107+
108+
n = length(Qidx)
109+
#If there are different Qs and Q's, the answer is a big fat 0
110+
if n != length(Qpidx) return zeros(size(eval(X.args[2]),1)...) end
111+
112+
##################################
113+
# Step 2. Enumerate permutations #
114+
##################################
115+
AllTerms = :(+(A*B))
116+
delete!(AllTerms.args, 2)
117+
println("BEGINS WITH ", AllTerms)
118+
for sigma in @task permutations_in_Sn(n)
119+
for tau in @task permutations_in_Sn(n)
120+
sigma_inv=permutation_inverse(sigma)
121+
#Compose the permutations
122+
perm=compose(sigma_inv, tau)
123+
124+
#Compute cycle structure, i.e. lengths of each cycle in the cycle
125+
#factorization of the permutatation
126+
cyclestruct = cycle_structure(perm)
127+
Qr = Int64[n+1 for n in Qidx]
128+
Qpr= Int64[Qpidx[n]+1 for n in data(tau)]
129+
#Consolidate deltas
130+
Deltas=Dict()
131+
for i=1:n
132+
Deltas[Qr[i]] = Qpr[i]
133+
Deltas[Qidx[i]] = Qpidx[data(sigma)[i]]
134+
end
135+
#Print deltas
136+
print("V(", cyclestruct, ") ")
137+
for k in keys(Deltas)
138+
print("δ(",k,",",Deltas[k],") ")
139+
end
140+
for (Symb, col_idx, row_idx) in Others
141+
print(Symb,"(",col_idx,",",row_idx,") ")
142+
end
143+
println()
144+
print("= V(", cyclestruct, ") ")
145+
#Substitute
146+
ReindexedSymbols = []
147+
for (Symb, col_idx, row_idx) in Others
148+
new_col, new_row = col_idx, row_idx
149+
if has(Deltas, col_idx) new_col = Deltas[col_idx] end
150+
if has(Deltas, row_idx) new_row = Deltas[row_idx] end
151+
print(Symb,"(",new_col,",",new_row,") ")
152+
ReindexedSymbols = [ReindexedSymbols; (Symb, new_col, new_row)]
153+
end
154+
println()
155+
println(ReindexedSymbols)
156+
println()
157+
#Reconstruct expression
158+
println("START PARSING")
159+
Symb, left_idx, right_idx = delete!(ReindexedSymbols, length(ReindexedSymbols))
160+
Expression={{{Symb}, left_idx, right_idx}}
161+
println("E =", Expression)
162+
while length(ReindexedSymbols) > 0
163+
pop_idx = expr_idx = do_transpose = is_left = nothing
164+
for expr_iter in enumerate(Expression)
165+
expr_idx, expr_string = expr_iter
166+
_, left_idx, right_idx = expr_string
167+
for iter in enumerate(ReindexedSymbols)
168+
idx, data = iter
169+
Symb, col_idx, row_idx = data
170+
if row_idx == left_idx
171+
pop_idx = idx
172+
do_transpose = false
173+
is_left = true
174+
println("Case A")
175+
break
176+
elseif col_idx == right_idx
177+
pop_idx = idx
178+
do_transpose = false
179+
is_left = false
180+
println("Case B")
181+
break
182+
elseif row_idx == right_idx
183+
pop_idx = idx
184+
do_transpose = true
185+
is_left = false
186+
println("Case C")
187+
break
188+
elseif col_idx == left_idx
189+
pop_idx = idx
190+
do_transpose = true
191+
is_left = true
192+
println("Case D")
193+
break
194+
end
195+
end
196+
if pop_idx != nothing break end
197+
end
198+
println("Terms left = ", ReindexedSymbols)
199+
if pop_idx == nothing #Found nothing, start new expression blob
200+
println("NEW EXPRBLOB: The term is =", ReindexedSymbols[1])
201+
Symb, left_idx, right_idx = delete!(ReindexedSymbols, 1)
202+
insert!(Expression, length(Expression)+1, {{Symb}, left_idx, right_idx})
203+
else #Found something
204+
println("The term is =", ReindexedSymbols[pop_idx])
205+
Symb, col_idx, row_idx = delete!(ReindexedSymbols, pop_idx)
206+
Term = do_transpose ? Expr(symbol("'"), Symb) : Symb
207+
println(Expression[expr_idx])
208+
println(Expression[expr_idx][1])
209+
if is_left
210+
println("insert left")
211+
insert!(Expression[expr_idx][1], 1, Term)
212+
Expression[expr_idx][2] = do_transpose ? row_idx : col_idx
213+
else
214+
println("insert right")
215+
insert!(Expression[expr_idx][1],
216+
length(Expression[expr_idx][1])+1, Term)
217+
Expression[expr_idx][3] = do_transpose ? col_idx : row_idx
218+
end
219+
end
220+
println("E=",Expression)
221+
end
222+
println("DONE PARSING TREE")
223+
224+
# Evaluate closed cycles
225+
NewExpression=:(+(A*B))
226+
delete!(NewExpression.args, 2)
227+
println("HIHIHI", NewExpression.args)
228+
for ExprBlob in Expression
229+
if ExprBlob[2]==ExprBlob[3] #Have a cycle; this is a trace
230+
println(ExprBlob[1])
231+
ex=:(trace(A*B))
232+
println(ex.args[2].args)
233+
ex.args[2].args={:*; ExprBlob[1]...}
234+
println(ex)
235+
println(ex.args)
236+
insert!(NewExpression.args, length(NewExpression.args)+1, ex)
237+
else #Not a cycle, regular chain of multiplications
238+
ex=:(A*B)
239+
ex.args={:*; ExprBlob[1]...}
240+
insert!(NewExpression.args, length(NewExpression.args)+1, ex)
241+
end
242+
end
243+
println("Final expression is ", NewExpression)
244+
insert!(AllTerms.args, length(AllTerms.args)+1, NewExpression)
245+
end
246+
end
247+
println("THE ANSWER IS ", AllTerms)
248+
eval(AllTerms)
249+
end
250+
251+
252+
#Iterate over partitions of n in lexicographic order
253+
function part(n::Integer)
254+
if n==0 produce([]) end
255+
if n<=0 return end
256+
for p in @task part(n-1)
257+
p = [p; 1]
258+
produce(p)
259+
p = p[1:end-1]
260+
if length(p) == 1 || (length(p)>1 && p[end]<p[end-1])
261+
p[end] += 1
262+
produce(p)
263+
p[end] -= 1
264+
end
265+
end
266+
end
267+
268+
269+
270+
###############
271+
# SANDBOX
272+
N=5
273+
A=randn(N,N)
274+
B=randn(N,N)
275+
type HaarMatrix
276+
beta::Real
277+
end
278+
Q = HaarMatrix(2)
279+
280+
281+
println("Case 1")
282+
println("E(A*B) = ", Expectation(:(A*B)))
283+
println("A*B/N = ", A*B)
284+
285+
println("Case 2")
286+
println("tr(A)*tr(B)/N = ", trace(A)*trace(B)/N)
287+
#println("E(A*Q*B*Q') = ", Expectation(:(A*Q*B*Q')))
288+
println("E.tr(A*Q*B*Q') = ", trace(Expectation(:(A*Q*B*Q'))))
289+
290+
#println("Case 3")
291+
#println(Expectation(:(A*B)) == A*B)
292+
#println("E(A*Q*B*Q'*A*Q*B*Q') = ", Expectation(:(A*Q*B*Q'*A*Q*B*Q')))
293+

0 commit comments

Comments
 (0)