`quadgk` and Unitful numbers

[ajkeller34]

Prompted by JuliaPhysics/Unitful.jl#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.

[tkelman]

quadgk really doesn't belong in base any more, ref #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

[stevengj]

@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.

[tkelman]

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

[tkelman]

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

[ajkeller34]

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]

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

[tkelman]

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.

[stevengj]

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.

[ChrisRackauckas]

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?

[tkelman]

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.

[ajkeller34]

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).

[ChrisRackauckas]

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

#18795 (comment)

@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?).