Strange type instability in array comprehension

General Usage
1

Hi all,

I am developing StructuralVibration.jl and I am currently try to eliminate all the type instabilities in the code. While some fixes are obvious. I am facing a type instability that I can’t explain easily. The code below aims at computing the stabilization diagram of an experimental modal analysis problem.

import StructuralVibration as SV

struct StabilizationAnalysis{Tc <: Complex, Tf <: Real}
    prob::MdofProblem    
    poles::Vector{Vector{Tc}} 
    modefn::Matrix{Tf}
    mode_stabfn::BitMatrix
    mode_stabdr::BitMatrix
end

function stabilization(prob::MdofProblem, max_order::Int, alg::Union{MdofEMA, MdofOMA} = LSCF(); stabcrit = [0.01, 0.05], progress = false)

    # Extract FRF and frequency from problem
    if prob isa EMAProblem
        frf = prob.frf
    elseif prob isa OMAProblem
        frf = prob.halfspec
    end

    # Initialization
    max_order += 1 
    poles = [similar(frf, i) for i in 1:max_order] # Type instability detected by Cthulhu
    fn = [similar(prob.freq, i) for i in 1:max_order]
    dr = [similar(prob.freq, i) for i in 1:max_order]
    modefn = fill(NaN, max_order, max_order)
    mode_stabfn = falses(max_order, max_order)
    mode_stabdr = falses(max_order, max_order)

    for order in 1:max_order
        poles[order] = SV.poles_extraction(prob, order, alg, stabdiag = true)

        fne, dre = poles2modal(poles[order])
        fn[order] .= fne
        dr[order] .= dre
        modefn[1:order, order] .= fn[order]

        if order > 1
            Nm = length(fn[order-1])
            mode_stabfn[1:Nm, order-1], mode_stabdr[1:Nm, order-1] = SV.check_stability(fn[order], fn[order-1], dr[order], dr[order-1], eltype(prob.freq).(stabcrit))
        end
    end

    # Remove last order (not used for stability check)
    max_order -= 1
    return SV.StabilizationAnalysis(prob, poles[1:max_order], modefn[1:max_order, 1:max_order], mode_stabfn[1:max_order, 1:max_order], mode_stabdr[1:max_order, 1:max_order])
end

Here is an MWE:

import StructuralVibration as SV

# Structure parameters of the beam
L = 1.        # Length
b = 0.03      # Width
h = 0.01      # Thickness
S = b*h       # Cross-section area
Iz = b*h^3/12 # Moment of inertia

# Material parameters
E = 2.1e11  # Young's modulus
ρ = 7850.   # Density
ξ = 0.01    # Damping ratio

# Mesh
xexc = 0:0.05:L
xm = xexc[2]

# Mode calculation - Simply supported boundary conditions
beam = SV.Beam(L, S, Iz, E, ρ)
fmax = 500.

ωn, kn = SV.modefreq(beam, 2fmax)
ms_exc = SV.modeshape(beam, kn, xexc)
ms_m = SV.modeshape(beam, kn, xm)

# FRF calculation
freq = 1.:0.1:fmax
prob = SV.ModalFRFProblem(ωn, ξ, freq, ms_m, ms_exc)
H = SV.solve(prob).u

stab = stabilization(prob_mdof, order, LSCF()) # Type unstable function

According to Cthulhu.jl, the culprit seems to this line poles = [similar(frf, i) for i in 1:max_order]. It is inferred as Vector instead of Vector{Vector{eltype(frf}}. I also have to say that fn and dn are properly inferred by the compiler.

Do you have an idea to solve this type instability?

Thank you

2

After a closer inspection, it seems that the problem comes from these lines

if prob isa EMAProblem
    frf = prob.frf
elseif prob isa OMAProblem
    frf = prob.halfspec
end

Interestingly, the type instability disappears when I use the ternary operator

frf = prob isa EMAProblem ? prob.frf : (prob isa OMAProblem ? prob.halfspec : throw(ArgumentError("Unsupported problem type for stabilization analysis.")))

Do you have some explanation about this?

Thank you

1 Like
3

Comprehensions are based on closures, and the one you highlighted captured the variable frf. The current lowering of that closure to a functor boxes frf because it gets assigned in >1 location, even if it only occurs once at runtime and would otherwise be inferred perfectly at compile-time. You should see a frf::Core.Box in reflection, and that is completely uninferred like a Ref{Any}. If a conditional is too difficult to rewrite into a ?: expression, you can just assign the value to another variable that is never reassigned to be captured instead.

2 Likes
4

I think the following should also work (although I didn’t test whether it’s also stable):

frf = if prob isa EMAProblem
    prob.frf
elseif prob isa OMAProblem
    prob.halfspec
end
3 Likes
5

You are right. Your proposal is type stable. I think it’s in line with @Benny explanations.