Skip to content

QuadGK is type unstable for ULogFloat64 #21

Open
Open
QuadGK is type unstable for ULogFloat64#21

Description

@vandenman
vandenman
opened on Aug 27, 2024

I noticed a somewhat odd interaction with QuadGK.jl and a type instability for ULogFloat64.

Below is a MWE.

import LogarithmicNumbers, QuadGK

QuadGK.quadgk(x -> x^2, LogarithmicNumbers.LogFloat64(0),  LogarithmicNumbers.LogFloat64(1))
# (+exp(-1.0986122886681098), 0.0)
QuadGK.quadgk(x -> x^2, LogarithmicNumbers.ULogFloat64(0), LogarithmicNumbers.ULogFloat64(1))
# (+exp(-1.0986122886681098), 0.0)

both work, although I noticed that the error term has type Float64 instead of LogFloat64. However,

@code_warntype QuadGK.quadgk(x -> x^2, LogarithmicNumbers.ULogFloat64(0), LogarithmicNumbers.ULogFloat64(1))

shows type instabilities. Specifically,

import LinearAlgebra
f(x) = x^2
segs = (zero(LogarithmicNumbers.ULogFloat64), one(LogarithmicNumbers.ULogFloat64))
atol=nothing; rtol=nothing; maxevals=10^7; order=7; segbuf=nothing; eval_segbuf=nothing
@code_warntype QuadGK.handle_infinities(f, segs) do f, s, _
    QuadGK.do_quadgk(f, s, order, atol, rtol, maxevals, LinearAlgebra.norm, segbuf, eval_segbuf)
end

is type unstable:

MethodInstance for QuadGK.handle_infinities(::var"#19#20", ::typeof(f), ::Tuple{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64})
  from handle_infinities(workfunc, f, s) @ QuadGK ~/.julia/packages/QuadGK/wnfvV/src/adapt.jl:163
Arguments
  #self#::Core.Const(QuadGK.handle_infinities)
  workfunc::Core.Const(var"#19#20"())
  f::Core.Const(f)
  s::Tuple{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
Locals
  #14::QuadGK.var"#14#23"{LogarithmicNumbers.ULogFloat64}
  #13::QuadGK.var"#13#22"
  #12::QuadGK.var"#12#21"{typeof(f), LogarithmicNumbers.ULogFloat64}
  inf2::Bool
  inf1::Bool
  s2::LogarithmicNumbers.ULogFloat64
  s1::LogarithmicNumbers.ULogFloat64
  #20::QuadGK.var"#20#29"{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
  #19::QuadGK.var"#19#28"{LogarithmicNumbers.ULogFloat64}
  #18::QuadGK.var"#18#27"{typeof(f), LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
  #17::QuadGK.var"#17#26"{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
  #16::QuadGK.var"#16#25"{LogarithmicNumbers.ULogFloat64}
  #15::QuadGK.var"#15#24"{typeof(f), LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
  @_18::Int64
  si::LogarithmicNumbers.ULogFloat64
  s0::LogarithmicNumbers.ULogFloat64
  @_21::Tuple{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
Body::Tuple
1 ──        Core.NewvarNode(:(#14))
│           Core.NewvarNode(:(#13))
│           Core.NewvarNode(:(#12))
│           Core.NewvarNode(:(inf2))
│           Core.NewvarNode(:(inf1))
│    %6   = Base.getindex(s, 1)::LogarithmicNumbers.ULogFloat64
│    %7   = Base.lastindex(s)::Core.Const(2)
│    %8   = Base.getindex(s, %7)::LogarithmicNumbers.ULogFloat64
│           (s1 = %6)
│           (s2 = %8)
│    %11  = QuadGK.realone(s1)::Core.Const(true)
└───        goto #16 if not %11
2 ── %13  = QuadGK.realone(s2)::Core.Const(true)
└───        goto #16 if not %13
3 ── %15  = QuadGK.isinf(s1)::Bool
│    %16  = QuadGK.isinf(s2)::Bool
│           (inf1 = %15)
│           (inf2 = %16)
└───        goto #5 if not inf1
4 ──        goto #6
5 ──        goto #16 if not inf2
6 ┄─        goto #9 if not inf1
7 ──        goto #9 if not inf2
8 ── %24  = QuadGK.:(var"#12#21")::Core.Const(QuadGK.var"#12#21")
│    %25  = Core.typeof(f)::Core.Const(typeof(f))
│    %26  = Core.typeof(s1)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %27  = Core.apply_type(%24, %25, %26)::Core.Const(QuadGK.var"#12#21"{typeof(f), LogarithmicNumbers.ULogFloat64})
│           (#12 = %new(%27, f, s1))
│    %29  = #12::QuadGK.var"#12#21"{typeof(f), LogarithmicNumbers.ULogFloat64}
│           (#13 = %new(QuadGK.:(var"#13#22")))
│    %31  = #13::Core.Const(QuadGK.var"#13#22"())
│    %32  = QuadGK.map(%31, s)::Tuple{Union{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.LogFloat64}, Union{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.LogFloat64}}
│    %33  = QuadGK.:(var"#14#23")::Core.Const(QuadGK.var"#14#23")
│    %34  = Core.typeof(s1)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %35  = Core.apply_type(%33, %34)::Core.Const(QuadGK.var"#14#23"{LogarithmicNumbers.ULogFloat64})
│           (#14 = %new(%35, s1))
│    %37  = #14::QuadGK.var"#14#23"{LogarithmicNumbers.ULogFloat64}
│    %38  = (workfunc)(%29, %32, %37)::Tuple
└───        return %38
9 ┄─        goto #11 if not inf1
10 ─        (@_21 = Core.tuple(s2, s1))
└───        goto #12
11 ─        (@_21 = Core.tuple(s1, s2))
12 ┄ %44  = @_21::Tuple{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
│           Core.NewvarNode(:(#20))
│           Core.NewvarNode(:(#19))
│           Core.NewvarNode(:(#18))
│           Core.NewvarNode(:(#17))
│           Core.NewvarNode(:(#16))
│           Core.NewvarNode(:(#15))
│    %51  = Base.indexed_iterate(%44, 1)::Core.PartialStruct(Tuple{LogarithmicNumbers.ULogFloat64, Int64}, Any[LogarithmicNumbers.ULogFloat64, Core.Const(2)])
│           (s0 = Core.getfield(%51, 1))
│           (@_18 = Core.getfield(%51, 2))
│    %54  = Base.indexed_iterate(%44, 2, @_18::Core.Const(2))::Core.PartialStruct(Tuple{LogarithmicNumbers.ULogFloat64, Int64}, Any[LogarithmicNumbers.ULogFloat64, Core.Const(3)])
│           (si = Core.getfield(%54, 1))
│    %56  = si::LogarithmicNumbers.ULogFloat64
│    %57  = QuadGK.zero(si)::Core.Const(exp(-Inf))
│    %58  = (%56 < %57)::Bool
└───        goto #15 if not %58
13 ─ %60  = QuadGK.:(var"#15#24")::Core.Const(QuadGK.var"#15#24")
│    %61  = Core.typeof(f)::Core.Const(typeof(f))
│    %62  = Core.typeof(s1)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %63  = Core.typeof(s0)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %64  = Core.apply_type(%60, %61, %62, %63)::Core.Const(QuadGK.var"#15#24"{typeof(f), LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64})
│    %65  = s1::LogarithmicNumbers.ULogFloat64
│           (#15 = %new(%64, f, %65, s0))
│    %67  = #15::QuadGK.var"#15#24"{typeof(f), LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
│    %68  = QuadGK.:(var"#16#25")::Core.Const(QuadGK.var"#16#25")
│    %69  = Core.typeof(s0)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %70  = Core.apply_type(%68, %69)::Core.Const(QuadGK.var"#16#25"{LogarithmicNumbers.ULogFloat64})
│           (#16 = %new(%70, s0))
│    %72  = #16::QuadGK.var"#16#25"{LogarithmicNumbers.ULogFloat64}
│    %73  = QuadGK.map(%72, s)::Tuple{LogarithmicNumbers.LogFloat64, LogarithmicNumbers.LogFloat64}
│    %74  = QuadGK.reverse(%73)::Tuple{LogarithmicNumbers.LogFloat64, LogarithmicNumbers.LogFloat64}
│    %75  = QuadGK.:(var"#17#26")::Core.Const(QuadGK.var"#17#26")
│    %76  = Core.typeof(s1)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %77  = Core.typeof(s0)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %78  = Core.apply_type(%75, %76, %77)::Core.Const(QuadGK.var"#17#26"{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64})
│    %79  = s1::LogarithmicNumbers.ULogFloat64
│           (#17 = %new(%78, %79, s0))
│    %81  = #17::QuadGK.var"#17#26"{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
│    %82  = (workfunc)(%67, %74, %81)::Tuple
└───        return %82
14 ─        Core.Const(:(goto %108))
15 ┄ %85  = QuadGK.:(var"#18#27")::Core.Const(QuadGK.var"#18#27")
│    %86  = Core.typeof(f)::Core.Const(typeof(f))
│    %87  = Core.typeof(s1)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %88  = Core.typeof(s0)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %89  = Core.apply_type(%85, %86, %87, %88)::Core.Const(QuadGK.var"#18#27"{typeof(f), LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64})
│    %90  = s1::LogarithmicNumbers.ULogFloat64
│           (#18 = %new(%89, f, %90, s0))
│    %92  = #18::QuadGK.var"#18#27"{typeof(f), LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
│    %93  = QuadGK.:(var"#19#28")::Core.Const(QuadGK.var"#19#28")
│    %94  = Core.typeof(s0)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %95  = Core.apply_type(%93, %94)::Core.Const(QuadGK.var"#19#28"{LogarithmicNumbers.ULogFloat64})
│           (#19 = %new(%95, s0))
│    %97  = #19::QuadGK.var"#19#28"{LogarithmicNumbers.ULogFloat64}
│    %98  = QuadGK.map(%97, s)::Tuple{LogarithmicNumbers.LogFloat64, LogarithmicNumbers.LogFloat64}
│    %99  = QuadGK.:(var"#20#29")::Core.Const(QuadGK.var"#20#29")
│    %100 = Core.typeof(s1)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %101 = Core.typeof(s0)::Core.Const(LogarithmicNumbers.ULogFloat64)
│    %102 = Core.apply_type(%99, %100, %101)::Core.Const(QuadGK.var"#20#29"{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64})
│    %103 = s1::LogarithmicNumbers.ULogFloat64
│           (#20 = %new(%102, %103, s0))
│    %105 = #20::QuadGK.var"#20#29"{LogarithmicNumbers.ULogFloat64, LogarithmicNumbers.ULogFloat64}
│    %106 = (workfunc)(%92, %98, %105)::Tuple
└───        return %106
16 ┄ %108 = (workfunc)(f, s, QuadGK.identity)::Tuple
└───        return %108

the highlighting on GitHub doesn't quite show it but all output from workfunc returns an untyped Tuple.

I played around with some ways to resolve this but it's not entirely clear to me what the best way to fix this would be. In addition, I had a few questions.

  1. Why is it that Logarithmic <: Real, instead of Logarithmic <: AbstractFloat?
  2. Why isn't float defined as float(x::Logarithmic) = x?

As a learning project, I ported the R package Brobdingag to Julia (which does essentially the same as this package). There, I did make those decisions and then QuadGK worked without issues, although I'm sure there are reasons for doing otherwise.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions