ForwardDiff & "complicated" functions?

BLI · November 14, 2019, 5:59pm

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:

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

Tamas_Papp · November 14, 2019, 6:17pm

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

http://www.juliadiff.org/ForwardDiff.jl/latest/user/api.html#ForwardDiff.derivative

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)

ssfrr · November 14, 2019, 6:22pm

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

ssfrr · November 14, 2019, 6:31pm

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.

BLI · November 14, 2019, 9:00pm

Ah. A typo.

BLI · November 15, 2019, 12:59pm

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?

Tamas_Papp · November 15, 2019, 1:15pm

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).

BLI · November 15, 2019, 1:42pm

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?

DrPapa · November 15, 2019, 1:59pm

@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.

BLI · November 15, 2019, 2:25pm

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.

ChrisRackauckas · November 15, 2019, 3:13pm

It’s part of the tutorial now:

http://docs.juliadiffeq.org/latest/tutorials/advanced_ode_example/#Automatic-Derivation-of-Jacobian-Functions-1

BLI · November 15, 2019, 3:42pm

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.

ChrisRackauckas · November 15, 2019, 3:50pm

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.

longemen3000 · November 16, 2019, 5:52am

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

BLI · November 16, 2019, 10:15am

Thanks! I’ll check it out!