ForwardDiff: Error ForwardDiff.Tag How to transform a Dual Type to a Float

Hallo,

I am struggling to get the simple derivatives works. I looked at many examples of ForwardDiff but still failed to solve my problem.

I am solving a Jacobian of a Triagonal Matrix by using numerical derivatives instead of deriving the derivatives manually which is a pain (worked but lack of flexibility). This is to solve the 1D Hydrology model of the Richards’ equation by computing the derivatives numerically.

The example below fails for Ψ₂[iT,iZ] = ψ . My understanding is that there is a type problem. ForwardDiff uses ψ = Dual Type and Ψ₂[ = Array{Float64}

I Guess that an answer to my problem would be how to simply transform a Dual Type to a Float?

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	#		FUNCTION : ∂R∂Ψ_NUMERICAL
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		function ∂RESIDUAL∂Ψ_NUMERICAL(discret, hydro, iT, iZ, K_Aver, N_iZ, Q, Sorptivity, ΔHpond, ΔPr, ΔSink, ΔT, θ, Ψ₂)

			function FUNCTION(ψ, discret, hydro, iT, iZ, K_Aver, N_iZ, Q, Sorptivity, ΔHpond, ΔPr, ΔSink, ΔT, θ, Ψ₂)
				Ψ₂[iT,iZ] = ψ
				∂R∂Ψ = RESIDUAL_DERIVATIVE(discret, hydro, iT, iZ, K_Aver, N_iZ, Q, Sorptivity, ΔHpond, ΔPr, ΔSink, ΔT, θ,  Ψ₂)
				return ∂R∂Ψ
			end  # function: Function
					
			ψ =Vector{Float64}(undef, 1) 
			∂R∂Ψ_Numerical(ψ::Vector) = FUNCTION(abs(ψ[1]), discret, hydro, iT, iZ, K_Aver, N_iZ, Q, Sorptivity, ΔHpond, ΔPr, ΔSink, ΔT, θ, Ψ₂)

			ψ[1] = Ψ₂[iT,iZ]
			Func_∂R∂Ψ_Numerical = ψ -> ForwardDiff.gradient(∂R∂Ψ_Numerical, ψ)			
			
			∂R∂Ψ2 = Func_∂R∂Ψ_Numerical(ψ)[1]

			return ∂R∂Ψ2 
		end # function: ∂RESIDUAL∂Ψ_NUMERICAL
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	#		FUNCTION : RESIDUAL
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		function RESIDUAL(discret, hydro, iT::Int64, iZ::Int64, Q, ΔSink, ΔT, θ, Ψ)::Float64
			return Residual = discret.ΔZ[iZ] * ((θ[iT,iZ] - θ[iT-1,iZ]) - hydro.So[iZ] * (Ψ[iT,iZ] - Ψ[iT-1,iZ]) * (θ[iT,iZ] / hydro.θs[iZ])) - ΔT[iT] * (Q[iT,iZ] - Q[iT,iZ+1]) + ΔSink[iT,iZ]
		end  # function: RESIDUAL

The error message:

ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{Main.residual.var"#∂R∂Ψ_Numerical#3"{Main.discretization.DISCRET,Main.hydroStruct.KOSUGI,Int64,Int64,Array{Float64,1},Int64,Array{Float64,2},Float64,Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,2},Array{Float64,2},Main.residual.var"function#2"},Float64},Float64,1})
Closest candidates are:
Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
Float64(::T) where T<:Number at boot.jl:716
Float64(::Irrational{:sqrt2}) at irrationals.jl:189

Many thanks for any help you may provide,
Joseph

If Ψ₂ contains dual numbers, they won’t fit into ψ, hence the error. There’s no automatic conversion as this would silently lose the gradient information.

If you must make a 1-element vector, just write [Ψ₂[iT,iZ]] or fill(Ψ₂[iT,iZ], 1), which will infer the correct eltype. But perhaps you want ForwardDiff.derivative?

julia> ForwardDiff.gradient(x -> sin(abs(x[1])), [1])
1-element Vector{Float64}:
 0.5403023058681398

julia> ForwardDiff.derivative(sin∘abs, 1)
0.5403023058681398

You have a lot of named internal functions which you only call once. I suspect the whole thing can be written more clearly as one ForwardDiff.gradient(...) do x ... block, with far fewer opportunities to get the list of arguments slightly wrong.

1 Like

A great thanks for your time and help. To clarify Ψ₂[iT,iZ] is an array which has values of past time step and of values of pressure from adjacent cells which are needed. Therefore fill(Ψ₂[iT,iZ], eltype(ψ)) would unfortunately not be a solution.

What is the difference between ForwardDiff.derivative and ForwardDiff.gradient ?

When I use ForwardDiff.derivative OR ForwardDiff.gradient I still get the sane error.

https://www.juliadiff.org/ForwardDiff.jl/latest/user/api.html

Generally derivatives are for \mathbb{R} \to \mathbb{R} function, while gradients for \mathbb{R}^n \to \mathbb{R}.

Thanks for providing me with this answer.

So how can I transform a Dual Type to a Float?

You can use ForwardDiff.value.

julia> v=ForwardDiff.Dual(1,2)
Dual{Nothing}(1,2)

julia> ForwardDiff.value(v)
1

I

Thanks for providing an improvement. By introducing Ψ₂[iT,iZ] = ForwardDiff.value(ψ) it solves the error problem but the derivative are now zero.

I believe that for it to work one must transform ∂R∂Ψ from Float64 into Dual, but I would not know how to do it?

			function FUNCTION(ψ, discret, hydro, iT, iZ, K_Aver, N_iZ, Q, Sorptivity, ΔHpond, ΔPr, ΔSink, ΔT, θ, Ψ₂)
				Ψ₂[iT,iZ] = ForwardDiff.value(ψ)
				∂R∂Ψ = RESIDUAL_DERIVATIVE(discret, hydro, iT, iZ, K_Aver, N_iZ, Q, Sorptivity, ΔHpond, ΔPr, ΔSink, ΔT, θ,  Ψ₂)
				return ∂R∂Ψ
			end  # function: Function

Yes, this is to be expected. Please reread the above:

Generally, you need to keep the dual information, and preallocate arrays accordingly. It is hard to provide more specific help without an MWE.

Thanks mcabbott, j-fu, and Tamas :smiley:for your suggestions and for confirming that what I was trying to achieve was not possible.

I finally managed to make the derivative work but as you suggested I had to rewrite the code to avoid mixing Vector{Float64} with type dual number, which is somewhat a pain but I guess that it is price to pay to get the fastest derivative package available in JULIA!

I confirmed that in my case using ForwardDiff.derivative to be faster than ForwardDiff.gradient

import ForwardDiff: derivative

function ∂R∂Ψ_FORWARDDIFF(param, Ψ)   
         ψ = Ψ

         ∂R∂Ψ_Func(ψ) = RESIDUAL(param , ψ)[1]

         ∂R∂Ψ_Derivative_1 = ψ -> derivative(∂R∂Ψ_Func, ψ)  

         return ∂R∂Ψ = ∂R∂Ψ_Derivative_1(ψ)
      end # function: ∂RESIDUAL∂Ψ_NUMERICAL ```