Correctly forming nested dual numbers.

[vpuri3]

Correctly forming nested dual numbers.

[mcabbott]

I think this is perturbation confusion. Making Dual{Nothing} means you are circumventing this package's mechanisms for keeping track of multiple perturbations. The user-facing functions seem to work fine:

julia> ForwardDiff.derivative(x -> ForwardDiff.derivative(f, x), 3.0)
2.0

julia> ForwardDiff.derivative(x -> ForwardDiff.derivative(f, x), 4.0)
2.0

[vpuri3]

I see. I've changed my code to use Dual.(Dual.(x, true), true). Thank you.

[vpuri3]

Update: z = Dual.(Dual.(x, 1v), Dual.(1v, 0v)) gives the right results.

For x -> x .^ 2

4-element Vector{Dual{Nothing, Dual{Nothing, Float64, 1}, 1}}:
  Dual{Nothing}(Dual{Nothing}(1.0,2.0),Dual{Nothing}(2.0,2.0))
  Dual{Nothing}(Dual{Nothing}(4.0,4.0),Dual{Nothing}(4.0,2.0))
  Dual{Nothing}(Dual{Nothing}(9.0,6.0),Dual{Nothing}(6.0,2.0))
 Dual{Nothing}(Dual{Nothing}(16.0,8.0),Dual{Nothing}(8.0,2.0))

For x -> exp.(x)

julia> include("examples/burgers_fourier/visc_burg_param_ic/autodecode.jl")
4-element Vector{Dual{Nothing, Dual{Nothing, Float64, 1}, 1}}:
  Dual{Nothing}(Dual{Nothing}(2.718281828459045,2.718281828459045),Dual{Nothing}(2.718281828459045,2.718281828459045))
  Dual{Nothing}(Dual{Nothing}(7.38905609893065,7.38905609893065),Dual{Nothing}(7.38905609893065,7.38905609893065))
 Dual{Nothing}(Dual{Nothing}(20.085536923187668,20.085536923187668),Dual{Nothing}(20.085536923187668,20.085536923187668))
 Dual{Nothing}(Dual{Nothing}(54.598150033144236,54.598150033144236),Dual{Nothing}(54.598150033144236,54.598150033144236))

[vpuri3]

@mcabbott I am getting similar behavior with different tags.

using ForwardDiff                                                               
using ForwardDiff: Dual, value, partials                                        

f = x -> x .^ 2                                                                 
x = [1.0, 2.0, 3.0, 4.0]                                                        
v = ones(4)                                                                     

# 2st order                                                                     
z = Dual{:FD_D2Tag}.(                                                           
    Dual{:FD_D2TagInt}.(x, 1v),                                                 
    Dual{:FD_D2TagInt}.(1v, 1v)                                                 
)                                                                               
fz = f(z)                                                                       
fx = value.(value.(fz))                                                         
df = value.(partials.(fz, 1))                                                   
d2f = partials.(partials.(fz, 1), 1)                                            

display(fz)                                                                              

4-element Vector{Dual{:FD_D2Tag, Dual{:FD_D2TagInt, Float64, 1}, 1}}:
 Dual{:FD_D2Tag}(Dual{:FD_D2TagInt}(1.0,2.0),Dual{:FD_D2TagInt}(2.0,4.0))
 Dual{:FD_D2Tag}(Dual{:FD_D2TagInt}(4.0,4.0),Dual{:FD_D2TagInt}(4.0,6.0))
 Dual{:FD_D2Tag}(Dual{:FD_D2TagInt}(9.0,6.0),Dual{:FD_D2TagInt}(6.0,8.0))
 Dual{:FD_D2Tag}(Dual{:FD_D2TagInt}(16.0,8.0),Dual{:FD_D2TagInt}(8.0,10.0))

[mcabbott]

Maybe it's not perturbation confusion, just wrong inputs? If you print out what the perturbations created by the user-facing function, you get this:

julia> f(x) = x .^ 2; f(x::Real) = @show(x) ^ 2;

julia> ForwardDiff.derivative(x -> ForwardDiff.derivative(f, x), 3.0)
x = Dual{ForwardDiff.Tag{typeof(f), Dual{ForwardDiff.Tag{var"#41#42", Float64}, Float64, 1}}}(Dual{ForwardDiff.Tag{var"#41#42", Float64}}(3.0,1.0),Dual{ForwardDiff.Tag{var"#41#42", Float64}}(1.0,0.0))
2.0

julia> Dual{:b}.(Dual{:a}.(x, 1), Dual{:a}.(1, 0)) |> f
4-element Vector{Dual{:b, Dual{:a, Float64, 1}, 1}}:
  Dual{:b}(Dual{:a}(1.0,2.0),Dual{:a}(2.0,2.0))
  Dual{:b}(Dual{:a}(4.0,4.0),Dual{:a}(4.0,2.0))
  Dual{:b}(Dual{:a}(9.0,6.0),Dual{:a}(6.0,2.0))
 Dual{:b}(Dual{:a}(16.0,8.0),Dual{:a}(8.0,2.0))

julia> Dual.(Dual.(x, 1), 1) |> f
4-element Vector{Dual{Nothing, Dual{Nothing, Float64, 1}, 1}}:
  Dual{Nothing}(Dual{Nothing}(1.0,2.0),Dual{Nothing}(2.0,2.0))
  Dual{Nothing}(Dual{Nothing}(4.0,4.0),Dual{Nothing}(4.0,2.0))
  Dual{Nothing}(Dual{Nothing}(9.0,6.0),Dual{Nothing}(6.0,2.0))
 Dual{Nothing}(Dual{Nothing}(16.0,8.0),Dual{Nothing}(8.0,2.0))

[vpuri3]

I see. So I was forming z incorrectly. It should be done like Dual.(Dual.(x, true), true) which is equivalent to what I was doing.

julia> Dual.(Dual.(x, true), true) == Dual.(Dual.(x, v), Dual.(v, v))
true

julia> Dual.(Dual.(x, true), true)
4-element Vector{Dual{Nothing, Dual{Nothing, Float64, 1}, 1}}:
 Dual{Nothing}(Dual{Nothing}(1.0,1.0),Dual{Nothing}(1.0,0.0))
 Dual{Nothing}(Dual{Nothing}(2.0,1.0),Dual{Nothing}(1.0,0.0))
 Dual{Nothing}(Dual{Nothing}(3.0,1.0),Dual{Nothing}(1.0,0.0))
 Dual{Nothing}(Dual{Nothing}(4.0,1.0),Dual{Nothing}(1.0,0.0))

Thanks @mcabbott for figuring this out

[vpuri3]

made a PR to add docs describing how to form nested duals.

[mcabbott]

which is equivalent to what I was doing.

No, this is not true. It's just that == on the tagged version of this package ignores duals. On master it does not:

julia> Dual.(Dual.(x, true), true)
4-element Vector{Dual{Nothing, Dual{Nothing, Float64, 1}, 1}}:
 Dual{Nothing}(Dual{Nothing}(1.0,1.0),Dual{Nothing}(1.0,0.0))
 Dual{Nothing}(Dual{Nothing}(2.0,1.0),Dual{Nothing}(1.0,0.0))
 Dual{Nothing}(Dual{Nothing}(3.0,1.0),Dual{Nothing}(1.0,0.0))
 Dual{Nothing}(Dual{Nothing}(4.0,1.0),Dual{Nothing}(1.0,0.0))

julia> Dual.(Dual.(x, v), Dual.(v, v))  # original guess above
4-element Vector{Dual{Nothing, Dual{Nothing, Float64, 1}, 1}}:
 Dual{Nothing}(Dual{Nothing}(1.0,1.0),Dual{Nothing}(1.0,1.0))
 Dual{Nothing}(Dual{Nothing}(2.0,1.0),Dual{Nothing}(1.0,1.0))
 Dual{Nothing}(Dual{Nothing}(3.0,1.0),Dual{Nothing}(1.0,1.0))
 Dual{Nothing}(Dual{Nothing}(4.0,1.0),Dual{Nothing}(1.0,1.0))

julia> Dual.(Dual.(x, true), true) == Dual.(Dual.(x, v), Dual.(v, v))  # on master
false