MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8})

[marouanehanhasse]

It is my first time trying this package and i tried to read the doc carefully and when i try to use on my code, i have a target function and a customized proposal distribution called qrs.psetproposal that takes a vector as an argument

using AdvancedHMC, ForwardDiff
alpha = pi/4
delta = pi
n=2
name="1"
chi = fill(sqrt(0.06), n)                  # the parameter chi 
phi = im*tanh.(chi)
omega = -log(prod(cosh.(chi)))
syms, op = qrelay_op(n, phi, alpha, delta)
op_a, op_ab, mat, coef = op_mat(op)

op_q2 = [syms.apH[1], syms.apV[1], syms.bpH[end], syms.bpV[end]]
op_q1 = [syms.apH[2:end]..., syms.apV[2:end]..., syms.bpH[1:end-1]..., syms.bpV[1:end-1]...]
mask_q1 = [op in op_q1 for op in op_a];

mask_q2 = [op in op_q2 for op in op_a];
qq = [x in syms.apH || x in syms.bpV ? 1 : 0 for x in op_a]

pdet0 = pdet_maker(0.04, 1e-5)
qrs = QRelaySampler(mat, coef, omega, pdet0)       #calling a function inside a module
targetcache=Dict{Vector{Int}, Float64}()

D= length(qq)
qq= float.(qq)
initial_θ= rand(qrs.psetproposal(qq))

ℓπ(θ) = log(qrs.prob(qq, θ, mask_q1) * qrs.prob(θ))      #the target function of MCMC

n_samples, n_adapts = 2_000, 1_000

metric = DiagEuclideanMetric(D)

hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)

initial_ϵ = find_good_stepsize(hamiltonian, initial_θ)

integrator = Leapfrog(initial_ϵ)

proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)

adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))

samples, stats = sample(hamiltonian, proposal, initial_θ, n_samples, adaptor, n_adapts; progress=true)

when i run it i get this error:

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8})
Closest candidates are:
  Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  Float64(::T) where T<:Number at boot.jl:732
  Float64(::Irrational{:invsqrt2π}) at irrationals.jl:189
  ...

Stacktrace:
 [1] convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8}) at ./number.jl:7
 [2] set_normalized_rhs(::JuMP.ConstraintRef{JuMP.Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}},JuMP.ScalarShape}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8}) at /Users/hanhasse/.julia/packages/JuMP/e0Uc2/src/constraints.jl:483
 [3] _broadcast_getindex_evalf at ./broadcast.jl:648 [inlined]
 [4] _broadcast_getindex at ./broadcast.jl:621 [inlined]
 [5] getindex at ./broadcast.jl:575 [inlined]
 [6] copy at ./broadcast.jl:876 [inlined]
 [7] materialize at ./broadcast.jl:837 [inlined]
 [8] setc at /Users/hanhasse/Downloads/scan.jl:15 [inlined]
 [9] (::Main.Quantum.var"#prob#8"{Array{Complex,1},Float64,Array{Float64,2},Array{Float64,2},Main.Quantum.var"#setc#5"{Array{JuMP.ConstraintRef{JuMP.Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}},JuMP.ScalarShape},1}},Main.Quantum.var"#scan#6"{JuMP.Model,Int64,Array{Int64,1},Array{Int64,1},Array{Int64,1},Array{JuMP.VariableRef,1}},Main.Quantum.var"#pdet0#2"{Float64,Float64}})(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8},1}) at ./In[2]:60
 [10] ℓπ(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8},1}) at ./In[4]:28
 [11] vector_mode_dual_eval at /Users/hanhasse/.julia/packages/ForwardDiff/QOqCN/src/apiutils.jl:37 [inlined]
 [12] vector_mode_gradient!(::DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}}, ::typeof(ℓπ), ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8},1}}) at /Users/hanhasse/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:113
 [13] gradient! at /Users/hanhasse/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:37 [inlined]
 [14] gradient!(::DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}}, ::typeof(ℓπ), ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8},1}}) at /Users/hanhasse/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:35 (repeats 2 times)
 [15] ∂ℓπ∂θ_forwarddiff at /Users/hanhasse/.julia/packages/AdvancedHMC/MIxdK/src/contrib/forwarddiff.jl:5 [inlined]
 [16] ∂ℓπ∂θ at /Users/hanhasse/.julia/packages/AdvancedHMC/MIxdK/src/contrib/forwarddiff.jl:46 [inlined]
 [17] ∂H∂θ at /Users/hanhasse/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:31 [inlined]
 [18] phasepoint(::Hamiltonian{DiagEuclideanMetric{Float64,Array{Float64,1}},typeof(ℓπ),AdvancedHMC.var"#∂ℓπ∂θ#44"{typeof(ℓπ)}}, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/hanhasse/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:69
 [19] find_good_stepsize(::Random._GLOBAL_RNG, ::Hamiltonian{DiagEuclideanMetric{Float64,Array{Float64,1}},typeof(ℓπ),AdvancedHMC.var"#∂ℓπ∂θ#44"{typeof(ℓπ)}}, ::Array{Float64,1}; max_n_iters::Int64) at /Users/hanhasse/.julia/packages/AdvancedHMC/MIxdK/src/trajectory.jl:813
 [20] #find_good_stepsize#13 at /Users/hanhasse/.julia/packages/AdvancedHMC/MIxdK/src/trajectory.jl:873 [inlined]
 [21] find_good_stepsize(::Hamiltonian{DiagEuclideanMetric{Float64,Array{Float64,1}},typeof(ℓπ),AdvancedHMC.var"#∂ℓπ∂θ#44"{typeof(ℓπ)}}, ::Array{Float64,1}) at /Users/hanhasse/.julia/packages/AdvancedHMC/MIxdK/src/trajectory.jl:873
 [22] top-level scope at In[4]:36
 [23] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091

and i really don't know how to fix it .

[Red-Portal]

seems like you're using JuMP somewhere within the target distribution. Can you show us qrs.prob?

[marouanehanhasse]

Sorry for the late answer, yes i am using JuMP ans here is the qrs.prob:

function prob(na)
            @assert count(!iszero, ui2*na) == 0
            b = T0*na
            setc(-b)
            total = 0.0
            for x in Channel(scan)
                nab = vi2*x + b #the photon numbers for each item in the sum in the note (10)
                total += prod([c.^complex(n)/factorial(n) for (c, n) in zip(coef, nab)])
            end
            return abs2(total*omega)
        end

or you can find it here https://github.com/marouanehanhasse/Quantum/blob/main/src/Quantum.jl as for the JuMP part it is here: https://github.com/marouanehanhasse/Quantum/blob/main/src/scan.jl

[Red-Portal]

Hi @marouanehanhasse , I'm not sure if current AD frameworks in Julia are expected to differentiate through JuMP. Is it normal in your domain to use gradient based MCMC? Also, what is the dimensionality (number of parameters) of your problem?

[marouanehanhasse]

In this problem, I have my target function also a customized distribution, as for the parameters when it comes to MCMC I need the initial vector as the input is an array and not a number

[xukai92]

It looks like AD fails. Can you double check if you can get the gradient of ℓπ outside AdvancedHMC.jl?

[marouanehanhasse]

ℓπ is working every time i don't use AdvancedHMC

[xukai92]

I mean the gradient of ℓπ (e.g. ForwardDiff.jl), not ℓπ itself.

[marouanehanhasse]

I'm really not that familiar of how to do that

[xukai92]

Please look at the README of ForwardDiff.jl (https://github.com/JuliaDiff/ForwardDiff.jl#forwarddiffjl). You should try replacing the f in the example with your ℓπ and see if ForwardDiff.jl can backprop through your function.

[marouanehanhasse]

when i do this

g = x -> ForwardDiff.gradient(ℓπ, x)
g(zeros(8))

i get this error

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8})
Closest candidates are:
  Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  Float64(::T) where T<:Number at boot.jl:732
  Float64(::Irrational{:invsqrt2π}) at irrationals.jl:189
  ...

Stacktrace:
 [1] convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8}) at ./number.jl:7
 [2] set_normalized_rhs(::JuMP.ConstraintRef{JuMP.Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}},JuMP.ScalarShape}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8}) at /Users/hanhasse/.julia/packages/JuMP/e0Uc2/src/constraints.jl:483
 [3] _broadcast_getindex_evalf at ./broadcast.jl:648 [inlined]
 [4] _broadcast_getindex at ./broadcast.jl:621 [inlined]
 [5] getindex at ./broadcast.jl:575 [inlined]
 [6] copy at ./broadcast.jl:876 [inlined]
 [7] materialize at ./broadcast.jl:837 [inlined]
 [8] setc at /Users/hanhasse/Downloads/scan.jl:15 [inlined]
 [9] (::Main.Quantum.var"#prob#8"{Array{Complex,1},Float64,Array{Float64,2},Array{Float64,2},Main.Quantum.var"#setc#5"{Array{JuMP.ConstraintRef{JuMP.Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}},JuMP.ScalarShape},1}},Main.Quantum.var"#scan#6"{JuMP.Model,Int64,Array{Int64,1},Array{Int64,1},Array{Int64,1},Array{JuMP.VariableRef,1}},Main.Quantum.var"#pdet0#2"{Float64,Float64}})(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8},1}) at ./In[2]:60
 [10] ℓπ(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8},1}) at ./In[5]:28
 [11] vector_mode_dual_eval at /Users/hanhasse/.julia/packages/ForwardDiff/QOqCN/src/apiutils.jl:37 [inlined]
 [12] vector_mode_gradient(::typeof(ℓπ), ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8},1}}) at /Users/hanhasse/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:106
 [13] gradient(::Function, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8},1}}, ::Val{true}) at /Users/hanhasse/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:19
 [14] gradient(::Function, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(ℓπ),Float64},Float64,8},1}}) at /Users/hanhasse/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:17 (repeats 2 times)
 [15] (::var"#15#16")(::Array{Float64,1}) at ./In[37]:2
 [16] top-level scope at In[37]:3
 [17] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091

so i cannot get the gradient of my target function

[xukai92]

Yes so it looks like your target is not differentiable, at least by ForwardDiff.jl. It may be due to the fact that @Red-Portal has pointed out: you cannot differentiate through JuMP. You should make sure you can get gradient even outside AdvancedHMC first as it's done by other packages.

In case you can hand-code your gradient function, you can pass it directly as well. Just to keep in mind that it should "return a tuple containing both the log-posterior and its gradient"; see https://github.com/TuringLang/AdvancedHMC.jl#gradients