Julia: `quadgk` and Unitful numbers

Created on 16 Dec 2016 · 16Comments · Source: JuliaLang/julia

Prompted by https://github.com/ajkeller34/Unitful.jl/issues/44

@stevengj The docstring for quadgk says:

The integrand f(x) can return any numeric scalar, vector, or matrix type,
or in fact any type supporting +, -, multiplication by real values, and a
norm (i.e., any normed vector space).

While Unitful numbers support +, -, multiplication by real values, and vecnorm, quadgk does not work with Unitful. I'd like to suggest making very small modifications to quadgk.jl, or at least clarifying the docstring to say that the result of the norm needs to be a real number.

The field named E in the type Base.QuadGK.Segment is restricted to Reals. This is problematic as unitful numbers are not real numbers, so I would remove the restriction. I presume it was intended to ensure that the type of E implements isless comparisons, but perhaps there's another way to do this without preventing unitful numbers from working?

With this small change, I just need to implement a quadgk method in Unitful to make quadgk compatible with unitful numbers. I've tried it out and it seems to work well. It is however possible that the default abstol will be dimensionally inconsistent with the results of integration. Without knowing the return type of the function being integrated I'm not sure one can do anything cleanly in Base (I guess you could try zero(f(x)*x) where x is some point within an integration interval?). I could always just suggest that Unitful users explicitly provide an abstol.

I'm happy to provide a PR if this change is acceptable.

maths

Source

ajkeller34

Most helpful comment

Virtually none of the numerical code is really "needed" by the core language. By this criterion we could move LinAlg out of Base too.

The whole point of the "default packages" idea is that we can move functionality out of Base without losing the "batteries included" nature of Julia, but this requires that we implement default-packages etc. before moving functionality out.

stevengj on 16 Dec 2016

👍2

All 16 comments

quadgk really doesn't belong in base any more, ref https://github.com/JuliaLang/julia/issues/18795 - would be best to move it to a package asap, then would be easy to make this kind of modification without being tied to Julia's release schedule

tkelman on 16 Dec 2016

@tkelman, you're putting the cart before the horse. Until a working mechanism for default packages is actually implemented, it seems premature to respond to issues by "this should be moved to a package ASAP."

I don't there is any particular need for the Real restriction in this type. It was intended for clarity, I think, but doesn't actually buy us anything. Would be an easy PR to remove this.

Regarding abstol, one possibility would be to define <=(x::Unitful, y::Number) = y == 0 ? x <= zero(x) : error(...), which is correct since 0 is the same for all units.

stevengj on 16 Dec 2016

It's not used by anything in base. If it doesn't need to be here, it shouldn't.

tkelman on 16 Dec 2016

Not having #18795 didn't stop pieces of combinatorics, primes, pkgdev etc from moving. Those all worked fine, so would this.

tkelman on 16 Dec 2016

👍1

I agree that it should move. I am sorry I don't have the bandwidth to help with that, but I do have the bandwidth to make a PR that deletes six characters from quadgk.jl. Am I discouraged from doing that? Should I trust that someone else will take care of this issue?

ajkeller34 on 16 Dec 2016

Thank you @stevengj for the suggestion regarding abstol, I think it is an interesting idea that may solve some problems for me elsewhere.

ajkeller34 on 16 Dec 2016

Making a small PR is fine, but if you want to rely on the result on non-master versions of Julia, moving the code to a package would allow you to do so sooner in this case.

tkelman on 16 Dec 2016

👍1

Virtually none of the numerical code is really "needed" by the core language. By this criterion we could move LinAlg out of Base too.

stevengj on 16 Dec 2016

👍2

The whole point of the "default packages" idea is that we can move functionality out of Base without losing the "batteries included" nature of Julia, but this requires that we implement default-packages etc. before moving functionality out.

Why can't this just be done using a metapackage with Reexport.jl?

ChrisRackauckas on 16 Dec 2016

By this criterion we could move LAPACK out of Base too.

Yes we should. At least make it optional. That's a lot harder though in terms of how much code would need to migrate, than something that's already its own module and completely decoupled.

tkelman on 16 Dec 2016

👍1

To clarify, I think that quadgk should move to a default package, not just anywhere... Otherwise we end up with the situation like linspace, where there is an implementation in base and one in Ranges.jl, and you have to explicitly tell people to use Ranges.jl with unitful numbers (I am of course very grateful for the work that went into that package but it should be a default package).

ajkeller34 on 16 Dec 2016

Otherwise we end up with the situation like linspace, where there is an implementation in base and one in Ranges.jl, and you have to explicitly tell people to use Ranges.jl with unitful numbers (I am of course very grateful for the work that went into that package but it should be a default package).

That's more of a problem due to the fact that the Base linspace was never deprecated. Correct me if I'm wrong, but I believe that's because it would break a lot of code (in Base), and because @timholy 's implementation was not accurate in the last digit.

quadgk is different. In

https://github.com/JuliaLang/julia/issues/18795#issuecomment-258494234

@tkelman mentions that quadgk is already pretty isolated, if a change did occur there may be no reason to have two implementations around since a Quadgk.jl would be the same implementation (as far as I can tell that would be the case?).

ChrisRackauckas on 16 Dec 2016

@ChrisRackauckas, I'm not saying it can't be done, but it hasn't been done yet. Furthermore, I'm concerned that our testing of changes to the Julia will be greatly degraded if we move a lot of functionality out of Base without infrastructure upgrades. And then there's the impact on load time, which again requires improved infrastructure for building "batteries included" system images.

That being said, I agree that the quadrature functionality is among the easiest to move out of base and would be the least missed by casual users. But in the meantime we should still accept patches like this.

stevengj on 16 Dec 2016

@ChrisRackauckas, I'm not saying it can't be done, but it hasn't been done yet. Furthermore, I'm concerned that our testing of changes to the Julia will be greatly degraded if we move a lot of functionality out of Base without infrastructure upgrades. And then there's the impact on load time, which again requires improved infrastructure for building "batteries included" system images.

I agree with you that things like LinAlg need the entirety of what's being discussed in #18795, but I think things more on the periphery like quadgk could already start the move and it wouldn't be a problem.

ChrisRackauckas on 16 Dec 2016

Moving quadgk out of Base allows heap{push,pop}! to be move into a package, e.g. DataStructures. I agree that we should have a default package mechanism but I'm really not entirely convinced that quadgk is a must-have for out-of-the-box Julia, especially considering how often even people needing to do basic quadrature are directed to packages beyond it.