# Why is this expression type unstable?

**URL:** <https://discourse.julialang.org/t/why-is-this-expression-type-unstable/80648>\
**Category:** General Usage\
**Tags:** question, type-stability, symbolics\
**Created:** [May 7, 2022, 4:02am UTC](https://discourse.julialang.org/t/why-is-this-expression-type-unstable/80648 "2022-05-07T04:02:52Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [May 7, 2022, 4:02am UTC](https://discourse.julialang.org/t/why-is-this-expression-type-unstable/80648/1 "2022-05-07T04:02:52Z")

</div>

I have the following function, and I’ve written two return statements which are identical in an input-output sense, but they have different syntax and mechanics. I like the format of the former statement because it’s a generator, and there’s less syntax there so it’s clearer. But the former return statement is not type stable, until I run the second return statement! The latter return statement is always type stable.

```Julia
function lagrangepolynomial(x::AbstractVector, y::AbstractVector, z::Number = Symbolics.variable(:z))

    n = length(x)

    δⱼ(z,x,j,k) = (z - x[k]) / (x[j] - x[k])
    Pⱼ(z,x,j,n) = prod(δⱼ(z,x,j,k) for k in 1:n if k != j)

    # NOT type stable...
    return sum((y[j] * Pⱼ(z,x,j,n)) for j in 1:n) 

    # Type stable
    return sum(map(j -> y[j] * Pⱼ(z,x,j,n), 1:n))

end

```

## Minimum Working Example

```Julia
using Symbolics
@variables z

x = -1:0.1:1
y = cos.(x)
n = length(x)

# Both of these are type stable!
δⱼ(z,x,j,k) = (z - x[k]) / (x[j] - x[k])
Pⱼ(z,x,j,n) = prod(δⱼ(z,x,j,k) for k in 1:n if k != j)

L1(z,x,y,n) = sum(y[j] * Pⱼ(z,x,j,n) for j in 1:n) 
L2(z,x,y,n) = sum(map(j -> y[j] * Pⱼ(z,x,j,n), 1:n))

```

## The first function is type unstable…

```Julia
julia> @code_warntype L1(z,x,y,n)
MethodInstance for L1(::Num, ::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, ::Vector{Float64}, ::Int64)
  from L1(z, x, y, n) in Main at REPL[19]:1
Arguments
  #self#::Core.Const(L1)
  z::Num
  x::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}
  y::Vector{Float64}
  n::Int64
Locals
  #25::var"#25#26"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64}
Body::Any
1 ─ %1 = Main.:(var"#25#26")::Core.Const(var"#25#26")
│ %2 = Core.typeof(z)::Core.Const(Num)
│ %3 = Core.typeof(x)::Core.Const(StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64})
│ %4 = Core.typeof(y)::Core.Const(Vector{Float64})
│ %5 = Core.typeof(n)::Core.Const(Int64)
│ %6 = Core.apply_type(%1, %2, %3, %4, %5)::Core.Const(var"#25#26"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64})
│ (#25 = %new(%6, z, x, y, n))
│ %8 = #25::var"#25#26"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64}
│ %9 = (1:n)::Core.PartialStruct(UnitRange{Int64}, Any[Core.Const(1), Int64])
│ %10 = Base.Generator(%8, %9)::Core.PartialStruct(Base.Generator{UnitRange{Int64}, var"#25#26"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64}}, Any[var"#25#26"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64}, Core.PartialStruct(UnitRange{Int64}, Any[Core.Const(1), Int64])])
│ %11 = Main.sum(%10)::Any
└── return %11

```

## The second function IS type stable…

```Julia
julia> @code_warntype L2(z,x,y,n)
MethodInstance for L2(::Num, ::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, ::Vector{Float64}, ::Int64)
  from L2(z, x, y, n) in Main at REPL[20]:1
Arguments
  #self#::Core.Const(L2)
  z::Num
  x::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}
  y::Vector{Float64}
  n::Int64
Locals
  #27::var"#27#28"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64}
Body::Num
1 ─ %1 = Main.:(var"#27#28")::Core.Const(var"#27#28")
│ %2 = Core.typeof(z)::Core.Const(Num)
│ %3 = Core.typeof(x)::Core.Const(StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64})
│ %4 = Core.typeof(y)::Core.Const(Vector{Float64})
│ %5 = Core.typeof(n)::Core.Const(Int64)
│ %6 = Core.apply_type(%1, %2, %3, %4, %5)::Core.Const(var"#27#28"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64})
│ (#27 = %new(%6, z, x, y, n))
│ %8 = #27::var"#27#28"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64}
│ %9 = (1:n)::Core.PartialStruct(UnitRange{Int64}, Any[Core.Const(1), Int64])
│ %10 = Main.map(%8, %9)::Vector{Num}
│ %11 = Main.sum(%10)::Num
└── return %11

```

## And now, the first function is type stable? Why?

```Julia
julia> @code_warntype L1(z,x,y,n)
MethodInstance for L1(::Num, ::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, ::Vector{Float64}, ::Int64)
  from L1(z, x, y, n) in Main at REPL[8]:1
Arguments
  #self#::Core.Const(L1)
  z::Num
  x::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}
  y::Vector{Float64}
  n::Int64
Locals
  #5::var"#5#6"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64}
Body::Num
1 ─ %1 = Main.:(var"#5#6")::Core.Const(var"#5#6")
│ %2 = Core.typeof(z)::Core.Const(Num)
│ %3 = Core.typeof(x)::Core.Const(StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64})
│ %4 = Core.typeof(y)::Core.Const(Vector{Float64})
│ %5 = Core.typeof(n)::Core.Const(Int64)
│ %6 = Core.apply_type(%1, %2, %3, %4, %5)::Core.Const(var"#5#6"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64})
│ (#5 = %new(%6, z, x, y, n))
│ %8 = #5::var"#5#6"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64}
│ %9 = (1:n)::Core.PartialStruct(UnitRange{Int64}, Any[Core.Const(1), Int64])
│ %10 = Base.Generator(%8, %9)::Core.PartialStruct(Base.Generator{UnitRange{Int64}, var"#5#6"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64}}, Any[var"#5#6"{Num, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, Int64}, Core.PartialStruct(UnitRange{Int64}, Any[Core.Const(1), Int64])])
│ %11 = Main.sum(%10)::Num
└── return %11

```

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [May 9, 2022, 3:02am UTC](https://discourse.julialang.org/t/why-is-this-expression-type-unstable/80648/2 "2022-05-09T03:02:37Z")

</div>

1. On v1.7.2, this phenomenon also happens when `Symbolics` is omitted and the global `z` is changed to `z=0`, so the problem doesn’t seem to be in `Symbolics`. The rest of my comment is about this `z=0` version.

2. Changing `L2` to `sum([...])` doesn’t change anything, so this might be a matter of summing over a generator versus over an array (which `map` makes).

3. The same phenomenon also happens if you run either `Pⱼ(z,x,1,n)`, `L1(z, x, y, n)`, or `L2(z, x, y, n)` instead of `@code_warntype L2(z, x, y, n)` between the 2 runs of `@code_warntype L1(z,x,y,n)`. Running `δⱼ(z,x,1,3)` was not enough.

4. If you run `Pⱼ(z,x,1,n)`, `L1(z, x, y, n)`, `L2(z, x, y, n)`, or `@code_warntype L2(z, x, y, n)` in a fresh session first, `@code_warntype L1(z,x,y,n)` is type-stable on its first try afterward.

5. These Github issues ([18695](https://github.com/JuliaLang/julia/issues/18695), [43330](https://github.com/JuliaLang/julia/issues/43330)) seem related, but I couldn’t find anything about `code_warntype` “changing its mind” like this in those issues or other linked issues. I would really highlight that detail, you could even edit your post’s title to mention it. You could also open a new Github issue, see where it goes.

Due to work done in [another thread](https://discourse.julialang.org/t/need-help-with-incomprehensible-typeinference/81056/6) about inconsistent inference, I also identified:

1. If `@code_warntype L1(z,x,y,n)` reports type-stability, then you redefine either `δⱼ` or `Pⱼ`, it will cause the next `@code_warntype L1(z,x,y,n)` to report type-instability. This re-unstabilizing effect does NOT occur if you redefine `L1` or `L2`.

2. `Test.@inferred L1(z, x, y, n)` can throw a `ERROR: return type Float64 does not match inferred return type Any` even when `@code_warntype` is inferring a `Float64` return type. `@inferred` is successful only if `Pⱼ` is compiled before `L1` is; switch the order (such as by only running `L1`) and `@inferred` gets stuck at that error; if it gets stuck at that error, redefine `L1` to “undo” `L1`’s compilation and it infers successfully again.

3. All of the above is reproducible on v1.4.1, so this isn’t a new regression.

My initial guess was that `Pⱼ` having been compiled seems to be a critical factor, but I could not claim that `@code_warntype L2(z,x,y,n)` compiles `Pⱼ`, let alone claim that’s how it triggers `@code_warntype L1(z,x,y,n)` to report type-stability.