kalmarek / Arblib.jl

Thin, efficient wrapper around Arb library (http://arblib.org/)
MIT License
25 stars 6 forks source link

Plotting and subtyping `AbstractFloat` #92

Open Joel-Dahne opened 3 years ago

Joel-Dahne commented 3 years ago

I want to be able to plot the Number types we define in Arblib.jl directly using Plots.jl. Currently plot([Arb(1)]) fails because float(::Arb) defaults to AbstractFloat(::Arb) which is not defined. There are two ways to solve this, either we define float(x::Arb) = x directly or we make Arb a subtype of AbstractFloat in which case AbstractFloat(::Arb) will work. Similar things happen for Arf.

For Arf I think it makes perfect sense to subtype AbstractFloat, after all it's essentially identical to a BigFloat which does subtype that. For Arb I'm not entirely sure, we could maybe see it as an AbstractFloat with the radius being some extra information. In practice I think this would work quite well, and everything that makes sense on an AbstractFloat I believe makes sense on an Arb as well. The alternative would be to define float(x::Arb) = x, so we would treat it as a float in practice but not subtype it, I'm not sure if that is better in any way. I can note here that float doesn't have to return a subtype of AbstractFloat, for example float(::Complex{Int}) is Complex{Float64}.

For Acb I believe we want to handle it similar to Complex{Arb}, which would mean adding these two plot recipes

RecipesBase.@recipe function f(A::Array{Acb})
    xguide --> "Re(x)"
    yguide --> "Im(x)"
    real.(A), imag.(A)
end

RecipesBase.@recipe function f(x::AbstractArray{T}, y::Array{Acb}) where {T<:Real}
    yguide --> "Re(y)"
    zguide --> "Im(y)"
    x, real.(y), imag.(y)
end

However I realised that many methods, at least in Julia Base, starts by doing float(x) on any input, so it could still be beneficial to define float(x::Acb) = Acb. I can also mention that it doesn't make sense to have Acb <: AbstractFloat since AbstractFloat <: Real.

Finally a comment about plotting balls. In the above plotting an Arb or an Acb would plot the midpoint of the ball. In many cases this makes sense, when working with high precision balls with small radius which is exactly what Arb is optimized for. This covers most of my uses. However in some cases you would like to plot the whole ball, including the radius. It would be nice if we could come up with a convenient method for doing that as well, possibly through some @userplot recipe.

Joel-Dahne commented 3 years ago

I can also mention that at some point I want to add plotting support to the polynomial and series types as well. But that is a future problem.

kalmarek commented 3 years ago

While it makes sense to have Arf <: AbstractFloat (I experimented with this idea while implementing rand) it'd be hard to think how Arbs may fit into AbstractFloat informal interface. From mathematical/cs point of view, probably Arf <: AbstractFloat and float(x::ArbOrRef) = midpoint(x) makes most sense? I have little opinion on this though.

Joel-Dahne commented 3 years ago

The way that Julia treats AbstractFloat seems to be just as something which is good at representing general real numbers, compared to for example Integer or Rational. If we don't subtype AbstractFloat then many of the methods defined in Base will not work out of the box. In some cases this might be good, there is no guarantee that the methods in Base actually give rigorous results. However many of them would work well in that case, for example

rad2deg(z::AbstractFloat) = z * (180 / oftype(z, pi))

from Base would work perfectly fine for Arb. I took a look at a sample of the Julia Base functions that use AbstractFloat at actually it seems like all of them would work fine for Arb.

Joel-Dahne commented 3 years ago

But I agree that it doesn't make perfect sense to have Arb <: AbstractFloat. In my current work I mostly treat it as a high precision floating point with the a little bit of extra information keeping track of the rounding errors, then it makes sense I think. But in other cases you really want to see it as a ball with a larger radius and then it makes less sense.

Defining float(x::Arb) = midpoint(x) would be problematic, for example due to functions like rad2deg above then returning an Arf. However defining float(x::Arb) = x) would also be problematic because we get a stack overflow instead. The only solution to this is to subtype AbstractFloat. But maybe these methods are not the most important to have to work.

I think we can agree that for Arf it makes sense to subtype AbstractFloat at least, I can make a pull request for that.

kalmarek commented 3 years ago

@saschatimme ?

saschatimme commented 3 years ago

I don't have a strong opinion. Since AbstractFloat <: Real we only would choose a more specific supertype for Arb. So this should be fine for me. Maybe more things work then and other things break at other places instead. So just go ahead and make Arf and Arb a suptype of AbstractFloat.

Joel-Dahne commented 3 years ago

Mostly solved by #108