New RECFAST doesn't like dual-numbers

[JaimeRZP]

Dear Jamie and Zack

After the latest update I get the following error when I try to ForwardDiff.jl through RECFAST:

Got exception outside of a @test
  MethodError: no method matching Bolt.RECFAST(::Bolt.Background{ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, Interpolations.ScaledInterpolation{ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, 1, Interpolations.BSplineInterpolation{ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, 1, OffsetArrays.OffsetVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Int64, ::Int64, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::Float64, ::Float64, ::Float64)

  Closest candidates are:
    Bolt.RECFAST(::AB, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::Any, ::Any, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T) where {T, AB<:(Bolt.AbstractBackground{T, IT} where IT<:(Interpolations.AbstractInterpolation{T, 1}))}

Code to replicate the error:

function get_Bolt_pk0(Ω_c)
    𝕡 = Bolt.CosmoParams(Ω_c = Ωc)
    bg = Background(𝕡)
    𝕣 = Bolt.RECFAST(bg=bg, Yp=𝕡.Y_p, OmegaB=𝕡.Ω_b,OmegaG=𝕡.Ω_r)
    ih = IonizationHistory(𝕣, 𝕡, bg)
    # bg.H₀ = 0.00023349 [1/Mpc] 
    kmin,kmax,nk = 0.5*0.00023349,30000*0.00023349,70 
    ks_Bolt = log10_k(kmin,kmax,nk)
    pL = [plin(k,𝕡,bg,ih) for k in ks_Bolt]
    return ks_Bolt, pL
end
ForwardDiff.derivative(get_Bolt_pk0, 0.25)

(not exactly what I am running but should the trick)

The speed up is amazing though!!!

Thanks a lot,
Jaime

[jmsull]

Thanks Jaime!

@xzackli I can fix this if I switch the types of the recfast struct arguments to ::Float64 for all the things not involving CosmoPararms. Is there a reason e.g. m_H is of type ::T in your new version?

[jmsull]

cf this version of the file (which works)

[xzackli]

I think the issue is that the new RECFAST wants the background as a non-kwarg,

𝕣 = Bolt.RECFAST(bg; Yp=𝕡.Y_p, OmegaB=𝕡.Ω_b,OmegaG=𝕡.Ω_r)

[JaimeRZP]

Problem solved using the new syntax. Thank you so much!
Anoher related question though, are we surprised that the derivative of Bolt.jl looks so different than the one from Eis&Hu?

[xzackli]

We need to fix this, something here is not right. Let's chat more about this, I'll get Uros to get me onto BCCP slack.

[jmsull]

Just taking derivatives wrt $\Omega_{c}$, I'm seeing no problems from our side?

Here is a plot and gist comparing these things:

https://gist.github.com/jmsull/d00eebca111eccaa6f93349cca9b35d5

For the Julia code I just minimially modified the plot_deriv_pk.jl examples/ file:

using Bolt,ForwardDiff,DelimitedFiles

function pkc(Ω_c::DT,k_grid) where DT
    𝕡 = CosmoParams{DT}(Ω_c=Ω_c)
    bg = Background(𝕡; x_grid=-20.0:0.01:0.0, nq=15)
    𝕣 = Bolt.RECFAST(bg; Yp=𝕡.Y_p, OmegaB=𝕡.Ω_b, OmegaG=𝕡.Ω_r)
    ih = IonizationHistory(𝕣, 𝕡, bg)
    nk=length(k_grid)
    a=zeros(DT,nk)
    for i = 1:nk
              a[i] = plin(k_grid[i],𝕡,bg,ih)[1] #not sure why this is returning a vector and not a T...
          end
   return a
end


Ωc0 = 0.27-0.046
𝕡 = CosmoParams(Ω_c=Ωc0,
                h = 0.7 ,
                Ω_r = 5.042e-5 ,
                Ω_b = 0.046  ,
                A = 2.097e-9 ,
                n = 1.0 ,
                Y_p = 0.24,  
                N_ν = 3.046 ,
                Σm_ν = 0.0 
                )
bg = Background(𝕡; x_grid=-20.0:0.01:0.0, nq=15)
kmin,kmax= 0.1bg.H₀*100,5000bg.H₀
k_grid = log10_k(kmin,kmax,33)
k_grid_hMpc = k_grid/(bg.H₀*3e5/100)
fc(Ω_c) = pkc(Ω_c,k_grid)
pk = fc(Ωc0)
∂pk = ForwardDiff.derivative(fc, Ωc0)  
Δ = 1e-2
finitediff_∂pk = (fc(Ωc0 + Δ) .- fc(Ωc0 - Δ)) ./ 2Δ

writedlm("/Users/jsull/Documents/berkeley/bolt/compare/bolt_camb_class_bolt_dOmegac.dat", 
            [k_grid_hMpc, pk, ∂pk, finitediff_∂pk])

[jmsull]

Gonna close this given my test here - there may be a problem with the Limberjack implementation that causes the discrepancy above