Typstabilität für eine Funktion mit Fallunterscheidungen

Question

Ich schreibe eine Funktion, die die Gewichte für die baryzentrische Interpolationsformel berechnet. Ignoriert Art Stabilität, das ist einfach genug:

function baryweights(x) 
    n = length(x) 
    if n == 1; return [1.0]; end # This is obviously not type stable 

    xmin,xmax = extrema(x) 
    x *= 4/(xmax-xmin) 
    #^Multiply by capacity of interval to avoid overflow 
    return [ 
     1/prod(x[i]-x[j] for j in 1:n if j != i) 
     for i = 1:n 
    ] 
end 

Das Problem für Typ Stabilität ist der Rückgabetyp der n > 1 Fall zu arbeiten, so kann ich ein Array des richtigen Typs im n == 1 Fall zurück. Gibt es einen einfachen Trick, um dies zu erreichen?

Answer 1

Rufen Sie einfach die Funktion rekursiv auf einem Dummy-Argument:

function baryweights(x) 
    n = length(x) 
    if n == 1 
     T = eltype(baryweights(zeros(eltype(x),2))) 
     return [one(T)] 
    end 

    xmin,xmax = extrema(x) 
    let x = 4/(xmax-xmin) * x 
     #^Multiply by capacity of interval to avoid overflow, 
     # and wrap in let to avoid another source of type instability 
     # (https://github.com/JuliaLang/julia/issues/15276) 
     return [ 
      1/prod(x[i]-x[j] for j in 1:n if j != i) 
      for i = 1:n 
     ] 
    end 
end 

Answer 2

Ich bin nicht sicher, ob ich Ihre Pläne zu verstehen. Aber vielleicht könnte so etwas helfen? ->

baryone(t::T) where T<:Real = [1.] 
baryone(t::T) where T<:Complex = [1im] # or whatever you like here 

function baryweights(x::Array{T,1}) where T<:Number 
    n = length(x) 
    n == 1 && return baryone(x[1]) 
    xmin,xmax = extrema(x) # don't forget fix extrema for complex! :) 
    x *= 4/(xmax-xmin) 
    #^Multiply by capacity of interval to avoid overflow 
    return [ 
     1/prod(x[i]-x[j] for j in 1:n if j != i) 
     for i = 1:n 
    ] 
end 

Warnung: Ich bin immer noch Neuling! Wenn ich versuche, @code_warntype baryweights([1]) sehe ich nur viele Warnungen. (Auch wenn ich es vermeide, baryone anzurufen). Zum Beispiel n ist Any !!

Edit: ich asked on discourse und sehe jetzt, dass @code_warn Rückkehr viel besseres Ergebnis, wenn wir eine andere Variable (y) verwenden:

function baryweights(x::Array{T,1}) where T<:Number 
    n = length(x) 
    n == 1 && return baryone(x[1]) 
    xmin,xmax = extrema(x) # don't forget fix extrema for complex! :) 
    let y = x * 4/(xmax-xmin) 
     #^Multiply by capacity of interval to avoid overflow 
     return [ 
      1/prod(y[i]-y[j] for j in 1:n if j != i) 
      for i = 1:n 
     ] 
    end 
end 

Edit2: Ich let hinzugefügt y zu vermeiden sein Core.Box ed