API for choosing power function

[dpsanders]

We need to decide on an API for how to specify whether powers will use the slow but accurate version (via MPFR) or the fast but less accurate version (via power_by_squaring).

Maybe something like

abstract type APIChoice end
abstract type PowerAPIChoice <: APIChoice end

struct AccuratePowers <: PowerAPIChoice end
struct FastPowers <: PowerAPIChoice end

IntervalArithmetic.specify!(AccuratePowers())
IntervalArithmetic.specify!(FastPowers())

Any other suggestions @lbenet @mforets @Kolaru?

[Kolaru]

Could we choose one at pre-compile time based on a specific key in the ENV dict and define the default method at tat point ? (the default would thus be constant unless you rebuild the module)

That's already what we do to determine whether or not to do validity checks. I think it is good unless there is a reason one would like to switch the default during a computation, but I can't see any use case (as long as both versions are available at all time through specific function names).

[mforets]

Since IntervalArithmetic has some setxyz methods, I imagine a function setpower(:MPFR) / setpower(:squaring) (or other names, such as :accurate / :fast).

Could we choose one at pre-compile time based on a specific key in the ENV dict and define the default method at tat point ?

How is that supposed to be used? The following obviously doesn't work, because ENV only checked at precompilation of the package -- is there a REPL command to re-trigger precompilation of the package?

julia> using IntervalArithmetic

julia> haskey(ENV, "IA_VALID")
false

julia> Interval(3, 2)
[3, 2]

julia> ENV["IA_VALID"] = true
true

julia> Interval(3, 2)
[3, 2]

I think it is good unless there is a reason one would like to switch the default during a computation

I agree; mixing different power methods in the same program seems like a rare use case.

[Kolaru]

Best I found is the following (starting with IA set with the validity check):

julia> using IntervalArithmetic

julia> Interval(3, 2)
ERROR: ArgumentError: Interval of form [3.0, 2.0] not allowed. Must have a ≤ b to construct interval(a, b).
Stacktrace:
 [1] Interval at /home/benoit/.julia/dev/IntervalArithmetic/src/intervals/intervals.jl:28 [inlined]
 [2] Interval at /home/benoit/.julia/dev/IntervalArithmetic/src/intervals/intervals.jl:42 [inlined]
 [3] Interval(::Int64, ::Int64) at /home/benoit/.julia/dev/IntervalArithmetic/src/intervals/intervals.jl:49
 [4] top-level scope at REPL[2]:1

julia> haskey(ENV, "IA_VALID")
false

julia> Base.compilecache(Base.PkgId(IntervalArithmetic))
[ Info: Precompiling IntervalArithmetic [d1acc4aa-44c8-5952-acd4-ba5d80a2a253]
"/home/benoit/.julia/compiled/v1.3/IntervalArithmetic/Gjmwo_x70GE.ji"

julia> exit()

Then after restarting julia:

julia> using IntervalArithmetic

julia> Interval(3, 2)
[3, 2]

Apparently the package manager build IntervalArithmetic is not enough, as it only runs a build.jl file (a file we do not explicitely create, so I don't know if it's possible to add stuff to it).

[mforets]

Cool, compilecache is the function that i was looking for, thanks. This should be documented (i opened https://github.com/JuliaIntervals/IntervalArithmetic.jl/issues/356).

[lbenet]

We need to decide on an API for how to specify whether powers will use the slow but accurate version (via MPFR) or the fast but less accurate version (via power_by_squaring).

I fully agree! Surely, not ideal, but we could exploit the already existing IntervalParameters struct, add a new field, and dispatch on that. This means to exploit a global constant (parameters), as we already do, until either we implement something with Cassette.jl, of find another alternative. If I recall correctly, the display options are similarly handled through a global constant.

[dpsanders]

Here's a possible solution / workaround using @dynamo from the IRTools.jl package, as was suggested on Discourse. This allows you to rewrite functions inside a block. The fast_powers function defined below rewrites ^ to pow inside the block:

julia> using IRTools: @dynamo, recurse!, IR

julia> @dynamo function fast_powers(a...)
                ir = IR(a...)
                ir == nothing && return
                recurse!(ir)
                return ir
              end

julia> fast_powers(::typeof(^), a, b) = pow(a, b);

julia> x = 1.1..2.2
Interval(1.0999999999999999, 2.2)

julia> x^10
Interval(2.5937424600999965, 2655.9922791424024)

julia> pow(x, 10)
Interval(2.593742460099994, 2655.992279142405)

julia> fast_powers() do
           x^10
       end
Interval(2.593742460099994, 2655.992279142405)

[dpsanders]

Given this, I would suggest making fast powers the default, and having an accurate_powers function along these lines.

[dpsanders]

Here's the same using Cassette.jl instead:

using Cassette, IntervalArithmetic

Cassette.@context FastPowersCtx
const fast_powers_ctx = Cassette.disablehooks(FastPowersCtx())

Cassette.overdub(::FastPowersCtx, ::typeof(^), a, b) = pow(a,b)

fast_powers(f) = Cassette.overdub(fast_powers_ctx, f)

x = 1.1..2.2

fast_powers() do
    x^10
end

[lbenet]

Very nice and very impressive!

I cannot comment about which approach to follow, but I think the default should use accurate_powers as default, simply because it is the one that complies with the standard.

[Kolaru]

As long as we don't introduce a new (mutable) global variables I'm fine with any API, but I must admit the use of context is very fancy =D.

Also the standard doesn't require anything about naming convention and iirc the function it defines must ensure correctness, but they are not required to be as tight as possible.

[mforets]

Thanks for sharing the different solutions with IRTools and Cassette. This illuminated me about some code transformations we wanted to do in LazySets, which seem to be achievable with this technology :)

Contextual dispatch avoids the more traditional/verbose pow(x, k, algorithm="fast"), but do we always have to write the code inside the fast_powers() do ... end block? I'm a bit confused about the workflow. Can I specify that fast_powers() works for the whole Julia session?

[dpsanders]

Great! Cassette is really super powerful.

Yes it has to be inside a I fast_powers() do block. Or rather, it has to be inside a Cassette context -- which is what fast_powers is using.

I think what you're effectively asking is "can I set a global context to use with Cassette". I think this will end up being equivalent to the kind of global manipulations we are already doing? Hmm, actually maybe that's a neater way of doing what we are trying to do.

This is related to #352.

[Kolaru]

I don't think we need the ability to set a global context. As far as i can see, the worst case scenario would be to have the whole program inside a do block (which can easily be done by putting everything in a single function and calling this one inside a do block). To me this doesn't justifies taking the risk to have a global option possibly causing unexpected side effects.

[OlivierHnt]

Done in PR #593. Let us see if the proposed API holds up to our expectations.