How can I fix the type instability

General Usage Performance
1

Hello,

I have the following MWE, where I cannot remove the type instability that appears from the call on the last line.

abstract type FrictionLaw end

struct Cubic{T <: Real} <: FrictionLaw
    γ::T
end

g(x, law::Cubic) = @. law.γ * x^3

function finite_diff_jac(f, x)
    xᵢ₊₁ = similar(x)
    fᵢ = f(x)
    fᵢ₊₁ = similar(fᵢ)
    J = similar(x, length(fᵢ), length(x))

    h = 1e-8
    for j in eachindex(x)
        xᵢ₊₁ .= x
        xᵢ₊₁[j] += h
        fᵢ₊₁ .= f(xᵢ₊₁)
        @. J[:, j] = (fᵢ₊₁ - fᵢ) / h
    end

    return J
end

function perform_aft(E, Eᴴ, X, friction_law)
    return Eᴴ * g(E * X, friction_law)
end

function friction_force_fourier(Xⱼ, xₚ, gₚ, friction_law, E, Eᴴ)
    Gⱼ = perform_aft(E, Eᴴ, [Xⱼ[1] + 2xₚ; Xⱼ[2:end]], friction_law)
    Gⱼ[1] -= 2gₚ
    return Gⱼ
end

H = 10

@code_warntype finite_diff_jac(u -> friction_force_fourier(u,
        0.1,
        0.1,
        Cubic(1.0),
        rand(2H + 1, 2H + 1),
        rand(2H + 1, 2H + 1)),
    ones(2H + 1))

Upon executing the code, I get

MethodInstance for finite_diff_jac(::var"#17#18", ::Vector{Float64})
  from finite_diff_jac(f, x) @ Main ~/ITP_Friction/HBMContinuation.jl/scripts/new/mwe_typeinst.jl:9
Arguments
  #self#::Core.Const(finite_diff_jac)
  f::Core.Const(var"#17#18"())
  x::Vector{Float64}
Locals
  @_4::Union{Nothing, Tuple{Int64, Int64}}
  h::Float64
  J::Union{Vector, Matrix{Float64}}
  fᵢ₊₁::Any
  fᵢ::Any
  xᵢ₊₁::Vector{Float64}
  j::Int64
Body::Union{Vector, Matrix{Float64}}

where the function vectors are of type Any, although I do not understand where this instability is coming from. Help will be appreaciated.

2

Avoid untyped global variables. H is used in the anonymous function, where the type instability then appears. I guess you can fix this by either adding a const when defining H, or adding a let H = H, ... around the definition of the anonymous function.

1 Like
3

While trying this @nsajko already posted the correct answer. Here is a REPL output to showcase the problem. When you just call the anonymous function like you defined here, nothing can be inferred because it captures the untyped global variable H

julia> @code_warntype (u -> friction_force_fourier(u,
               0.1,
               0.1,
               Cubic(1.0),
               rand(2H + 1, 2H + 1),
               rand(2H + 1, 2H + 1)))(
           ones(2H + 1))
MethodInstance for (::var"#6#7")(::Vector{Float64})
  from (::var"#6#7")(u) @ Main REPL[11]:1
Arguments
  #self#::Core.Const(var"#6#7"())
  u::Vector{Float64}
Body::Any
1 ─ %1  = Main.Cubic(1.0)::Core.Const(Cubic{Float64}(1.0))
│   %2  = (2 * Main.H)::Any
│   %3  = (%2 + 1)::Any
│   %4  = (2 * Main.H)::Any
│   %5  = (%4 + 1)::Any
│   %6  = Main.rand(%3, %5)::Any
│   %7  = (2 * Main.H)::Any
│   %8  = (%7 + 1)::Any
│   %9  = (2 * Main.H)::Any
│   %10 = (%9 + 1)::Any
│   %11 = Main.rand(%8, %10)::Any
│   %12 = Main.friction_force_fourier(u, 0.1, 0.1, %1, %6, %11)::Any
└──       return %12

So we can change this by surrounding it with let H=H like so:

fun = let H=H
        u -> friction_force_fourier(u,
               0.1,
               0.1,
               Cubic(1.0),
               rand(2H + 1, 2H + 1),
               rand(2H + 1, 2H + 1))
       end

This effectively converts the global H into a “local” H for the definition. So now everything is type-stable:

julia> @code_warntype fun(ones(2H + 1))
MethodInstance for (::var"#8#9"{Int64})(::Vector{Float64})
  from (::var"#8#9")(u) @ Main REPL[14]:2
Arguments
  #self#::var"#8#9"{Int64}
  u::Vector{Float64}
Body::Vector{Float64}

julia> @code_warntype finite_diff_jac(fun, ones(2H + 1))
MethodInstance for finite_diff_jac(::var"#8#9"{Int64}, ::Vector{Float64})
  from finite_diff_jac(f::F, x) where F @ Main REPL[9]:1
Static Parameters
  F = var"#8#9"{Int64}
Arguments
  #self#::Core.Const(finite_diff_jac)
  f::var"#8#9"{Int64}
  x::Vector{Float64}
Locals
  @_4::Union{Nothing, Tuple{Int64, Int64}}
  h::Float64
  J::Matrix{Float64}
  fᵢ₊₁::Vector{Float64}
  fᵢ::Vector{Float64}
  xᵢ₊₁::Vector{Float64}
  j::Int64
Body::Matrix{Float64}
1 Like
4

Of course. I was actually benchmarking some functions and I thought I had a type instability. It turns out that I didn’t check for global variables… Silly mistake

1 Like