JuliaStats / Loess.jl

Local regression, so smooooth!
Other
103 stars 34 forks source link

type instability in predict #51

Closed miguelbiron closed 1 year ago

miguelbiron commented 2 years ago

Using the example in the readme

using Loess

xs = 10 .* rand(100)
ys = sin.(xs) .+ 0.5 * rand(100)

model = loess(xs, ys, span=0.5)

us = range(extrema(xs)...; step = 0.1)
vs = predict(model, us)

This is the output from @code_warntype for the predict line. The type of many variables cannot be inferred by the compiler.

julia> @code_warntype predict(model, us)
Variables
  #self#::Core.Const(Loess.predict)
  model::Loess.LoessModel{Float64}
  zs::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}
  y2::Float64
  y1::Float64
  u::Any
  v₂::Any
  v₁::Any
  z::Float64
  adjacent_verts::Any
  m::Any

Body::Any
1 ──       Core.NewvarNode(:(y2))
│          Core.NewvarNode(:(y1))
│          Core.NewvarNode(:(u))
│          Core.NewvarNode(:(v₂))
│          Core.NewvarNode(:(v₁))
│          Core.NewvarNode(:(z))
│          Core.NewvarNode(:(adjacent_verts))
│    %8  = Base.getproperty(model, :xs)::AbstractMatrix{Float64}
│          (m = Loess.size(%8, 2))
│    %10 = (m == 1)::Any
└───       goto #4 if not %10
2 ── %12 = Loess.length(zs)::Int64
│    %13 = (%12 > 1)::Bool
└───       goto #4 if not %13
3 ── %15 = Loess.length(zs)::Int64
│    %16 = Core.tuple(%15, 1)::Core.PartialStruct(Tuple{Int64, Int64}, Any[Int64, Core.Const(1)])
│    %17 = Loess.reshape(zs, %16)::Core.PartialStruct(Base.ReshapedArray{Float64, 2, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, Tuple{}}, Any[StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, Core.PartialStruct(Tuple{Int64, Int64}, Any[Int64, Core.Const(1)]), Core.Const(())])
│    %18 = Loess.predict(model, %17)::Vector{Float64}
└───       return %18
4 ┄─ %20 = Loess.length(zs)::Int64
│    %21 = (%20 != m)::Any
└───       goto #6 if not %21
5 ── %23 = m::Any
│    %24 = Loess.length(zs)::Int64
│    %25 = Base.string(%23, "-dimensional model applied to length ", %24, " vector")::String
└───       Loess.error(%25)
6 ┄─ %27 = Base.getproperty(model, :kdtree)::Loess.KDTree{Float64}
│          (adjacent_verts = Loess.traverse(%27, zs))
│    %29 = (m == 1)::Any
└───       goto #15 if not %29
7 ── %31 = Loess.length(adjacent_verts)::Any
│    %32 = (%31 == 2)::Any
└───       goto #9 if not %32
8 ──       goto #10
9 ── %35 = Base.AssertionError("length(adjacent_verts) == 2")::Any
└───       Base.throw(%35)
10 ┄       (z = Base.getindex(zs, 1))
│    %38 = Base.getindex(adjacent_verts, 1)::Any
│    %39 = Base.getindex(%38, 1)::Any
│    %40 = Base.getindex(adjacent_verts, 2)::Any
│    %41 = Base.getindex(%40, 1)::Any
│          (v₁ = %39)
│          (v₂ = %41)
│    %44 = (z == v₁)::Any
└───       goto #12 if not %44
11 ─       goto #13
12 ─ %47 = (z == v₂)::Any
└───       goto #14 if not %47
13 ┄ %49 = Base.getproperty(model, :bs)::Matrix{Float64}
│    %50 = Base.getproperty(model, :verts)::Dict{Vector{Float64}, Int64}
│    %51 = Base.vect(z)::Vector{Float64}
│    %52 = Base.getindex(%50, %51)::Int64
│    %53 = Base.getindex(%49, %52, Loess.:(:))::Vector{Float64}
│    %54 = Loess.evalpoly(zs, %53)::Float64
└───       return %54
14 ─ %56 = (z - v₁)::Any
│    %57 = (v₂ - v₁)::Any
│          (u = %56 / %57)
│    %59 = Base.getproperty(model, :bs)::Matrix{Float64}
│    %60 = Base.getproperty(model, :verts)::Dict{Vector{Float64}, Int64}
│    %61 = Base.vect(v₁)::Vector{_A} where _A
│    %62 = Base.getindex(%60, %61)::Int64
│    %63 = Base.getindex(%59, %62, Loess.:(:))::Vector{Float64}
│          (y1 = Loess.evalpoly(zs, %63))
│    %65 = Base.getproperty(model, :bs)::Matrix{Float64}
│    %66 = Base.getproperty(model, :verts)::Dict{Vector{Float64}, Int64}
│    %67 = Base.vect(v₂)::Vector{_A} where _A
│    %68 = Base.getindex(%66, %67)::Int64
│    %69 = Base.getindex(%65, %68, Loess.:(:))::Vector{Float64}
│          (y2 = Loess.evalpoly(zs, %69))
│    %71 = (1.0 - u)::Any
│    %72 = (%71 * y1)::Any
│    %73 = (u * y2)::Any
│    %74 = (%72 + %73)::Any
└───       return %74
15 ─       Loess.error("Multivariate blending not yet implemented")
└───       Core.Const(:(return %76))
boriskaus commented 2 years ago

In one of my codes, Loess.jl has become a computational bottleneck, most likely due to the type-instability described here and addressed in #53. @andreasnoack: could you perhaps give that PR a look and see if it can be merged?

andreasnoack commented 1 year ago

Closed by #68