setrounding is slow and unreliable

[JeffBezanson]

setrounding works by redefining a lot of methods. I suppose that can be good for performance, but has the serious costs of (1) setrounding itself is very slow (including file IO as well!), (2) it only works from the top level or if the caller uses invokelatest, (3) it's global and therefore non-composable (another package can change it behind your back). I suggest either removing the function, changing how it works, or at least loudly documenting that it should only be used from the top level e.g. for interactive demos or exploration. Some packages (e.g. ExactPredicates.jl) have already tried to use it within their functions.

[dpsanders]

@JeffBezanson Thanks for the comment -- I totally agree.
I would love to remove / improve it, but I don't have any idea how to do so. Do you have any suggestions?

I fear that Cassette may be the answer, but I'm reluctant to go there.

[JeffBezanson]

Do we know what the overhead would be from checking a global flag on each operation?

Other than that, the complete set of options is probably

• Pass around a rounding mode manually
• Use a macro to insert a rounding mode argument within a lexical block
• Add it to the type of an Interval, so the rounding mode follows the data around, and type specialization takes care of it

[dpsanders]

The thing is it really is a kind of contextual dispatch -- that's even how it's currently implemented, e.g.

julia> +(IntervalRounding{:accurate}(), 0.1, 0.3, RoundDown)
0.39999999999999997

julia> +(IntervalRounding{:none}(), 0.1, 0.3, RoundDown)
0.4

The rounding is a property of the operations, not of the intervals.

[dpsanders]

(So the name IntervalRounding is actually a misnomer.)

[JeffBezanson]

Cassette it is then :)

[Kolaru]

Unless I'm missing something, we could have the rounding integrated to the type (e.g. Interval{RoundingMode} and then define things like

+(a::Interval{R}, b::Interval{R}) where {R <: RoundingMode} = Interval{R}( +(R, a.lo, b.lo, RoundDown), +(R, a.hi, b.hi, RoundUp))

[dpsanders]

We could do that, but it seems conceptually wrong, for the reason I wrote above (the rounding is a property of the implementation of the operation, not of the intervals themselves).

[ericphanson]

Sorry to butt in, but could the R in @Kolaru's suggestion be interpreted just as a default? I.e. an Interval{R} for R <: RoundingMode is an interval on which operations round according to R by default. To me, that seems conceptually a little better; rounding is part of the operation, but the interval itself knows how it generally should be rounded. It seems conceptually contextual dispatch fits well, but maybe it's not practical yet, and type-based dispatch could be a practical solution.

[Kolaru]

As a side note, the IEEE standard advises to consider the flavor of the interval to be a property of the context (section 7.1):

(I think that the trick to define the Interval type used when the package is build is compatible with that btw)

I didn't find anything about the kind of rounding used along this line. So, with regard to the standard, I think we are free here.

[lbenet]

Unless I'm missing something, we could have the rounding integrated to the type (e.g. Interval{RoundingMode}.

I am not sure if I understand the proposal, and for the moment (with my lack of understanding) agree with @dpsanders that it sounds conceptually wrong. For instance, what would it be the rounding mode R for Interval(0.25, 0.5) which needs no rounding?. Or what about Interval(0.1, 0.7)? The latter, in order to contain 1//10 and 7//10, requires that the infimum is rounded down and the supremum rounded up?

I guess we should spend some time learning Cassette.jl, or perhaps use some tricks such as rounding-down a number a is the same as rounding-up -a, which would avoid half the changes of the rounding mode.

[Kolaru]

I confused IntervalRounding and RoundingMode (the fact that IntervalRounding is the rounding mode of the interval doesn't help). So the example becomes:

+(a::Interval{R}, b::Interval{R}) where {R <: IntervalRounding} = Interval{R}( +(R, a.lo, b.lo, RoundDown), +(R, a.hi, b.hi, RoundUp))

where IntervalRounding allow to choose between accurate but slow, and inaccurate but fast (and some others IIRC).

In an example such as @interval(a, b) that could be translated in choosing between @round(a, RoundDown) .. @round(b, RoundUp) (tightest but may be expensive as it checks if a and b needs rounding) and prefloat(a)..nextfloat(b) (fast but not tight if either of the bound didn't need rounding to begin with).

[mewilhel]

I've run into issues 2 & 3 that @JeffBezanson lists in the past and this has been quite painful to recognize and deal with. I'd vote for the parameterizing intervals with a IntervalRounding approach and dispatch based on this parameter value as detailed by @ericphanson.

Isn't the main goal of associating a rounding mode with an execution environment just to avoid mixing rounding mode inappropriately?. If so, we should be able to avoid this just by throwing errors when performing an operation with two interval of distinct rounding modes.

[JeffBezanson]

Here's another possibility. Since it's a property of the operation and not the intervals, just parameterize the operations!

struct Plus{R<:IntervalRounding}
end

(::Plus{RoundUp})(a::Interval, b::Interval) = ...

Then you need e.g. let + = Plus{RoundUp}(), but that could be introduced by a macro. Of course that only handles lexically-scoped use cases, which I guess is not enough.

[dpsanders]

@mewilhel The problem is that the rounding mode is not a property of the intervals, but of the operations between intervals.

@JeffBezanson Interesting idea. But we worked hard to avoid the need to use macros...

[JeffreySarnoff]

Something has to give.

[mewilhel]

@mewilhel The problem is that the rounding mode is not a property of the intervals, but of the operations between intervals.

@dpsanders Conceptually I agree. I just can't think of any use cases were that distinction creates an issue that the parameterize intervals approach can't handle. Right now the functionality of setrounding is to redefine all functions to use a particular rounding mode. Since the IntervalArithmetic package operates by overloading this could be achieved by parameterizing Intervals and dispatching based on this parametric type, no? Plus, the user could alter rounding modes fairly easily in a program if desired.

I am not sure if I understand the proposal, and for the moment (with my lack of understanding) agree with @dpsanders that it sounds conceptually wrong. For instance, what would it be the rounding mode R for Interval(0.25, 0.5) which needs no rounding?. Or what about Interval(0.1, 0.7)? The latter, in order to contain 1//10 and 7//10, requires that the infimum is rounded down and the supremum rounded up?

I guess we should spend some time learning Cassette.jl, or perhaps use some tricks such as rounding-down a number a is the same as rounding-up -a, which would avoid half the changes of the rounding mode.

@lbenet Wouldn't just selecting a default rounding mode be sufficient to address this problem?

Interval(a::T, b::T) = Interval{T,IntervalRounding{:tight}}(a,b)

Is the user wanted to use rounding mode :none a specialized constructor could be added...

[Kolaru]

I think the use of context (typically through Cassette) is unavaoidable. Indeed, we do not want the user to be able to let interact Intervals with different rounding modes (as it would not make sense anyway), which kind of exclude having it as a parameter of the interval.

On the other hand global redefinitions have the undesirable side effect... of having side effects.

That let us with context. Reusing the kind of logic @dpsanders proposed in #355, I think something like the following should do the trick (yet unstested though), making the rounding mode an explicit parameter to the operation:

abstract type RoundingMode end
struct TightRounding <: RoundingMode end
struct AccurateRounding <: RoundingMode end

# Define operation for all rounding modes
Base.:+(::TightRounding, a::Interval, b::Interval) = # Do whatever is needed here for tight rounding
Base.:+(::AccurateRounding, a::Interval, b::Interval) = # Do whatever is needed here for accurate rounding

# Define the default
Base.+(a::Interval, b::Interval) = +(TightRounding(), a, b)

Cassette.@context AccurateRoundingCtx
const accurate_rounding_ctx = Cassette.disablehooks(AccurateRoundingCtx())
Cassette.overdub(::AccurateRoundingCtx, ::typeof(+), R::RoundingMode, a, b) = +(AccurateRounding(), a, b)

accurate_rounding(f) = Cassette.overdub(accurate_rounding_ctx, f)

accurate_rounding() do
    x = (1..2) + (3.4..4.5)
end

There is probably a way to just use type parameter rather than one type per rounding mode, but I need to think a bit more as to how to make it work with Cassette.

Also I suspect it possible to define the overdub function such that any occurrence of a RoundingMode object will be replaced, regardless of the function used, making it very compact.

[mewilhel]

I think the use of context (typically through Cassette) is unavaoidable. Indeed, we do not want the user to be able to let interact Intervals with different rounding modes (as it would not make sense anyway), which kind of exclude having it as a parameter of the interval.

@Kolaru I like this idea. That said, not defining Base.:+(a::Interval{T,R}, b::Interval{T,S}) for R !== S (R<: RoundingMode, S<: RoundingMode) or promotions between Interval{T,R} and Interval{T,R} would result in a method error and prevent the user from mixing different rounding modes. So I think, parameterizing Interval with a RoundingMode is a valid approach if one doesn't want to go the Cassette route.

[dpsanders]

I personally don't like the idea of adding a type parameter to Interval to specify the type of directed rounding. But if somebody wants to make a PR or separate package showing how that would work, we can see if it makes sense.

@Kolaru's Cassette-based solution is certainly a good and immediately-workable solution that we should implement (at the cost of adding a dependency on Cassette.jl to the package).

However, after thinking about this some more, I think maybe the issue has generated more discussion than it maybe warrants: the method used for directed rounding is a trade off between implementation difficulty and execution speed (as usual) that slightly affects some results. But for a user of the package, there is hardly ever (?) a reason to want to modify it once it has been set at the start of a session.

So I think I agree with @JeffBezanson in his suggestion that we should just advertise this more carefully. Note that setrounding was (I believe) never even documented, so I'm pretty surprised that people even ever found out about it!

I also think we might want to change the names setrounding and IntervalRounding, e.g. to set_directed_rounding_mode and DirectedRoundingMode, to avoid confusion and make them sound less like something you want to mess with.

[JeffreySarnoff]

change the names setrounding and IntervalRounding, e.g. to set_directed_rounding_mode and DirectedRoundingMode, to avoid confusion and make them sound less like something you want to mess with.

(on the CB radio) breaker breaker, that's a big 10-4

[mewilhel]

However, after thinking about this some more, I think maybe the issue has generated more discussion than it maybe warrants: the method used for directed rounding is a trade off between implementation difficulty and execution speed (as usual) that slightly affects some results. But for a user of the package, there is hardly ever (?) a reason to want to modify it once it has been set at the start of a session.

So I think I agree with @JeffBezanson in his suggestion that we should just advertise this more carefully. Note that setrounding was (I believe) never even documented, so I'm pretty surprised that people even ever found out about it!

@dpsanders Fair point. I agree that the user hardly ever wants to change the rounding mode in a session. So I think adding a warning (and example usage) to the docstring and documentation would be sufficient. Something like this work?

"""
    setrounding(Interval, rounding_type::Symbol)

THIS FUNCTION SHOULD ONLY BE APPLIED AT THE TOP LEVEL 
AND SHOULD NOT BE CALLED IN ANY DEPENDENT PACKAGES. 

Set the rounding type used for all interval calculations on Julia v0.6 and above.
Valid `rounding_type`s are $rounding_types. 

In order to set rounding modes using in a package Example.jl: 
> using IntervalArithmetic; setrounding(Interval, :accurate)
> using Example
"""

Personally, understanding why an algorithm I was working on suffered a performance hit is what drove me to look into the undocumented set_rounding function. In a few examples,the interval calculations required for constructing relaxations of subproblems in EAGO.jl actually took a 5-6x times longer than an LP solve with CPLEX in specific rounding modes (:tight) and 3x shorter in others (:accurate). This was the limiting step in that nonconvex solver. So this was a pretty big deal for run time comparisons. Once I had this result, it was fairly natural to try to include a setrounding function in that package. Apologies for not posting the exact example, this was a while ago.

[Kolaru]

Interestingly (i was the first surprised), my previous proposition can be easily adapted to avoid Cassette altogether:

abstract type RoundingMode end
struct TightRounding <: RoundingMode end
struct AccurateRounding <: RoundingMode end

# Define operation for all rounding modes
Base.:+(::TightRounding, a::Interval, b::Interval) = # Do whatever is needed here for tight rounding
Base.:+(::AccurateRounding, a::Interval, b::Interval) = # Do whatever is needed here for accurate rounding

# Define the defaults
default_directed_rounding_mode() = TightRounding()
Base.+(a::Interval, b::Interval) = +(default_directed_rounding_mode(), a, b)

# Change the default by redefining the associated function
IntervalArithmetc.default_directed_rounding_mode() = AccurateRounding()

It has the advantage of still explicitly parametrizing the operation, and discouraging its use as it requires explicit import and overwriting a library function. It may comes at the cost of a slight overhead in the default case though.

[dpsanders]

@Kolaru How does that differ from what's currently implemented? I think it's basically the same? (although maybe your version is more readable!)

[Kolaru]

setrounding has been removed in favor of IntervalArithmetic.interval_rounding().