Type-stability in Julia for creating SVector with ::Function

104 Views Asked by At

I'm trying to figure out the proper way to write code in a type-stable manner. I know that some limitations, such as ::Function and ::Val{} arguments are an issue, and so I'm trying to figure out the best way of going around those. I'm still deep into the learning phase, so please feel free to provide as much information and comments as you feel like.

As part of that, I was trying to create an SVector that is type-stable, however I keep getting it wrong. here is a MWE (and pardon the stupid names for functions and variables, I tried to maitain the actual structure of the code while removing meaningful content), where the output is SVector{_A, Float64} where _A. How do I get a proper SVector here? Feel free to modify the code as you please, as long as the function order and structure remains.

Supporting functions:

using StaticArrays
@inline function HyperTanDistribSymm(α::T,Lx::T=one(T)) where T<:AbstractFloat
    return x -> convert(T,0.5)*Lx*( one(T) - tanh( α*(one(T)-convert(T,2.0)*x) ) / tanh(α) )
end

function GenerateVec(L::TF, indxrng::AbstractUnitRange{TI}, ntotal::TI=length(indxrng), Stretchfun::Function=x::TF->x::TF, nghost::TI=1) where {TI<:Integer, TF<:AbstractFloat}
    nterms::TI=2*nghost+length(indxrng)+one(TI)
    return SVector{nterms,TF}( map( Stretchfun,L*convert.(TF, indxrng[1]-one(TI)-nghost:indxrng[end]+nghost)/convert(TF, ntotal)  ) )
end


function madness(varD::TF, indx::AbstractUnitRange{TI}, varN::TI, stretch::Function, nghost::TI) where {TF<:AbstractFloat, TI<:Integer}
    testvar=GenerateVec(varD, indx, varN, stretch, nghost)
    nterms::TI=length(testvar)-one(TI)
    diffvar=SVector{nterms,TF}( diff(testvar) )
end

And running examples:

using Cthulhu
@descend_code_warntype madness(1.0, 1:6, 6, x->x, 1)
@descend_code_warntype madness(2.0, 2:8, 10, HyperTanDistribSymm(2.20,6.0), 2)

as GenerateVec returns SVector{nterms::Int64,TF::Type{Float64}}::Type{SVector{_A, Float64}} where _A. I will comment that I tried other versions, including the use of Val approach on something on the lines of:

@inline CreateSVector(input::Base.AbstractVecOrTuple{T}) where {T<:Number} = CreateSVector(Val(length(input)), input)
@inline CreateSVector(::Val{N}, input::Base.AbstractVecOrTuple{T}) where {N,T<:Number} = SVector{N,T}( input )

this is just one example, I also tried using N just as a regular N<:Integer input and such. All result in the same issue.

Note: I have read around and figured that possibly the use of Tuples pass arguments and set sizes might be a solution - if that is the case, why?

1

There are 1 best solutions below

0
On

You might find a text I wrote useful in your quest for understanding type stability: https://www.juliabloggers.com/writing-type-stable-julia-code/

But the short version is as follows: a) the type of a StaticVector depends on its length. b) type stability means that types (and so StaticVectors') length are defined at compile time c) so the length must be either hard coded, or be given as a type parameter. Examples, where the parameter is called VecLen:

... ::Val{VecLen}}...) where {VecLen} 

or

... v::StaticArray{VecLen,Float64}...) where {VecLen}

You want type stability 'cause you want speed 'cause you're in a loop. If you can compute VecLen outside the loop (see the concept of barrier function) and pass it down to your function (using Val, or other technique) you win.

If you cannot precompute VecLen (e.g. because it changes at each iteration) then we are busted, you will not get type stability. You would then StaticArrays are not for you today.

I hope this helps!

: )

Philippe