I’m a little tired of debugging code so any easy solution will be greatly appreciated.
I’m trying to use the AdvancedHMC package (with Zygote for auto differentiation) to draw samples from a known PDF (from my work).
As per my understanding, the derivative is first invoked for the function defined below which then moves backwards to work on “wavefunction” (defined after)
function LogDensityProblems.logdensity(p::LogTargetDensity, Z)
### Z is an array of the form, x1 y1 x2 y2 ...
Zcomp = Z[begin:2:end] + 1.0im*Z[begin+1:2:end] ### Complexifying Z.
return log(abs2(wavefunction(Zcomp)))
end
function wavefunction(Z) ### N not necessary.
#### Given Z -> array of positions,
#### N -> no of electrons, n -> n[i]-> Number of electrons in ith level.
polynomials = Array{Polynomial}(undef, N) ### N
derivatives_collection = Array{ComplexF64}(undef, N)
map(i->polynomials[i]=fromroots(repeat(Z[1:end .!= i], p)),1:N)
slater_det = Matrix{ComplexF64}(undef, N, N)
### Is there a faster way to do this? Vectorization of some sort?? Trust JULIA FOR NOW.
row_index = 0
ll_index = 0
prefactor = 0.0
for i in range(1, length(n))
ll_index = i-1
# println(ll_index)
map(electron_index -> derivatives_collection[electron_index] = derivative(polynomials[electron_index], ll_index)(Z[electron_index]),1:N)
# println(derivatives_collection)
for j in range(1, n[i])
row_index = row_index + 1
#= Just the row in which to put something. m = j-1 -ll index, n[i] -> Landau Level Index + 1
We will calculate derivative for an n-index only once for each particle. That requires a new index.
=#
m = j - 1 - ll_index
prefactor = 1/√(factorial(big(ll_index))*factorial(big(m+ll_index))*2*π*2.0^m)
for column_index in range(1, N)
slater_det[row_index, column_index] = prefactor*(Z[column_index]^m)*derivatives_collection[column_index]
end
end
end
# slater_det = exp(-dot(Z, Z)/4) .* slater_det
return det(slater_det) * exp(-dot(Z, Z)/4)
end
When I try running this, the following is the result from the debugger:
MethodError: no method matching adjoint(::Zygote.Context{false}, ::typeof(sort), ::Vector{ComplexF64}; lt=Polynomials.var"#77#79"())
Closest candidates are:
adjoint(::ZygoteRules.AContext, ::typeof(sort), ::AbstractArray; by) at none:0 got unsupported keyword argument "lt"
adjoint(::ZygoteRules.AContext, !Matched::Base.Fix1, ::Any) at none:0 got unsupported keyword argument "lt"
adjoint(::ZygoteRules.AContext, !Matched::Base.Fix2, ::Any) at none:0 got unsupported keyword argument "lt"
I have seen posts on errors combining Polynomials. evalpoly and Zygote, and how Zygote fails to work with complex-valued functions. I’m unsure whether the error is a result of either of the above or some other implementation error from my end. Beyond this, I have thought of creating a struct for the complexf64 data type in terms of a pair of real numbers and going from there. However, that seems like a lot of work. Also, I could not use the Polynomials package, but that also seems like a little too much effort (I’m essentially trying things out to make things easier for my group and convenience functions, packages would be nice). The above function was the first draft, so feel free to suggest things to boost efficiency.