# ForwardDiff & "complicated" functions?

I’d need to compute the derivative (and later: the Jacobian) of the following function:

f = x -> sqrt(max(-x,0))


which looks like the following:

Obviously, the derivative doesn’t exist for x=0. Still, I try to use ForwardDiff:

fd = x -> ForwardDiff.derivative(f,x)


I then try to plot this function, leading to:

Upon checking, it turns out that fd(x) gives value NaN for x\in[0,\infty).
I then try to define the function as follows:

g = x -> begin
if x < 0
sqrt(abs(x))
else
0
end
end


and define the derivative as

gd = x -> ForwardDiff(g,x)


OK – this derivative doesn’t evaluate…

julia> gd(-2.)
MethodError: objects of type Module are not callable

Stacktrace:
[1] (::getfield(Main, Symbol("##42#43")))(::Float64) at .\In[113]:1
[2] top-level scope at In[130]:1


Questions:

1. Why doesn’t evaluation of gd(-2.) work?
2. Why does fd(x) evaluate to NaN for x\in [0,\infty)?
3. Is there any other package that can handle this problem?
4. …or is it better to insert a smooth transition from, say f(-10^{-6}), to f(0) or f(10^6)?
1 Like

You are not calling the package functions right. See the docs:

Also, you may want to make g type stable. Eg

using ForwardDiff
g(x) = x < 0 ? sqrt(abs(x)) : zero(x)
ForwardDiff.derivative(g, -1)
ForwardDiff.derivative(g, 1)

3 Likes

This should be ForwardDiff.derivative(g,x), (probably a typo?) which is causing your “Module not callable” error.

You can use the @which macro to see where the function max is defined which might explain the NaN.

As a side note - the style f = x-> ... is not very idiomatic in Julia. Usually that would be written as f(x) = .... For multi-line functions folks generally use the style:

function f(x)
...
end

2 Likes

actually the @edit macro might be more useful. Running @edit max(1.0, 0) brings you to the promote code that turns 0 into a float, so @edit max(1.0, 1.0) gives the code from math.jl:

max(x::T, y::T) where {T<:AbstractFloat} = ifelse((y > x) | (signbit(y) < signbit(x)),
ifelse(isnan(x), x, y), ifelse(isnan(y), y, x))


Something in there (maybe the isnan?) is making ForwardDiff unhappy.

Ah. A typo.

I’m trying to linearize an ODE model…using ForwardDiff, but don’t get the expected result. A simple example:

# Model with input functions
function watertank(dT,T,p,t,u)
# Parameters
g,h,d,per,A,rho,chp,Ux,Po = p
# Inputs
uP,uv,Ti,Ta,VdL = u
#
md = rho*uv*VdL
#
As = 2*A + per*h
V = A*h
m = rho*V
Cp = m*chp
Qd_el = Po*uP
Qd_a = Ux*As*(Ta-T[1])
dT[1] = (md*chp*(Ti-T[1]) + Qd_el + Qd_a)/Cp
#
return dT
end


Data for the model:

# -- constants
g = 9.81
h = 1.5
d = 0.5
per = pi*d
A = pi*d^2/4
rho = 1.0e3
chp = 4.19e3
Ux = 0.45
Po = 15e3
#
p = [g,h,d,per,A,rho,chp,Ux,Po];
#
T_op = [25.]
u_op = [0.5,0.5,28.,20.,8e-3/60];


Here, T_op and u_op constitute the operating point for the linearization. Next, I define functions which only depend on T and u:

wt_T(T) = watertank(copy(T),T,p,0,u_op)
wt_u(u) = watertank(copy(T_op),T_op,p,0,u)


These work as expected:

julia> wt_T(T_op)
1-element Array{Float64,1}:
0.006751564884010684


and similarly for wt_u(u_op).
Next, I introduce Jacobians:

J_T(T) = ForwardDiff.jacobian(wt_T,T)
J_u(u) = ForwardDiff.jacobian(wt_u,u)


The Jacobian wrt. T works as expected:

julia> J_T(T_op)
1×1 Array{Float64,2}:
-0.00022735608347665157


HOWEVER, the Jacobian wrt. u gives an error mesage:

julia> J_u(u_op)
MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5})
Closest candidates are:
Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:194
Float64(::T<:Number) where T<:Number at boot.jl:718
Float64(!Matched::Int8) at float.jl:60
...

Stacktrace:
[1] convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5}) at .\number.jl:7
[2] setindex!(::Array{Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5}, ::Int64) at .\array.jl:766
[3] watertank(::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5},1}) at .\In[2]:16
[4] wt_u(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5},1}) at .\In[4]:3
[5] vector_mode_jacobian(::typeof(wt_u), ::Array{Float64,1}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5},1}}) at C:\Users\Bernt_Lie\.julia\packages\ForwardDiff\N0wMF\src\apiutils.jl:37
[6] jacobian(::Function, ::Array{Float64,1}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5},1}}, ::Val{true}) at C:\Users\Bernt_Lie\.julia\packages\ForwardDiff\N0wMF\src\jacobian.jl:17
[7] jacobian(::Function, ::Array{Float64,1}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(wt_u),Float64},Float64,5},1}}) at C:\Users\Bernt_Lie\.julia\packages\ForwardDiff\N0wMF\src\jacobian.jl:15 (repeats 2 times)
[8] J_u(::Array{Float64,1}) at .\In[4]:4
[9] top-level scope at In[8]:1


What am I doing wrong?

I think that for wt_u, dT is an array of Float64s, so it cannot accommodate the ForwardDiff.Duals.

You could make your code non-modifying (ie reallocate a new dT of the right type), or just calculate dT[1] as a scalar variable and return that (as far as I understand, that’s the only thing that is changed).

Hm. “Jacobian” normally involves a vector valued function in several variables. Doesn’t dT being an array simply mean that the function is vector valued?

The documentation says that the function should return an AbstractArray… [This method assumes that isa(f(x), AbstractArray)]. Is Float64 a subset of AbstractArray?

@BLI this doesn’t directly address your issue with ForwardDiff, but you could consider using ModelingToolkit.jl to specify the ODE and get the Jacobian that way.

2 Likes

I could probably use ModelingToolkit.jl. I’m really interested in ModelingToolkit.jl, but It is not really documented for “the average user” at the moment.

I can, of course, find the Jacobian numerically:

B = zeros(length(T_op),length(u_op))
Iu = Diagonal(ones(length(u_op)))
for j in 1:length(u_op)
B[:,j] = (wt_u(u_op+Iu[:,j]*1e-6) - wt_u(u_op))/1e-6
end


but it would be much more elegant to use ForwardDiff or some other tool.

It’s part of the tutorial now:

4 Likes

I’ll look into it. But I need something more general… I’m considering systems with inputs, \frac{dx}{dt} = f(x,u;t) and I also need the Jacobian wrt. the input, i.e., I need both \frac{\partial f}{\partial x} and \frac {\partial f}{\partial u}. The DiffEq solvers assume that I have made the mapping g(x,t) = f(x,u(x,t),t), so I assume that ModelingToolkit will find the Jacobian of g(x,t) wrt. x, only.

Oh yes, if your parameters are like that you can use ModelingToolkit to derive the Jacobian, but the automatic version needs some work to handle your case.

this works without modification of watertank

wt_T!(dT,T) = watertank(dT,T,p,0.0,u_op)
wt_u!(dT,u) = watertank(dT,T_op,p,0.0,u)
J_T(dT,T) = ForwardDiff.jacobian(wt_T!,dT,T)
J_u(dT,u) = ForwardDiff.jacobian(wt_u!,dT,u)
J_u(u) = J_u(copy(T_op),u)
J_T(T) = J_u(copy(T),T)


the main problem was treating a function ´f!(result,x)´ as ´f(x)´, but ForwardDiff can work with mutable functions if you specify that:

julia> J_u(u_op)
1×5 Array{Float64,2}:
0.012155  0.00135812  0.000226354  1.00239e-6  5.09296


now, i want to see that cp as a function of T

2 Likes

Thanks! I’ll check it out!