Docs: how to handle branching?

[petvana]

This is more like a question. I haven't found any information in the docs about handling branching while using InrervalArithmetic. For example

julia> using IntervalArithmetic

julia> f(x) = x > 1 ? x-1 : x
f (generic function with 1 method)

julia> f(0.9..1.1)
[0.899999, 1.10001] # expecting [0.0, 1.0]

julia> f.([0.9, 1.1])
2-element Vector{Float64}:
 0.9
 0.10000000000000009

The code produces an incorrect interval because 0.9..1.1 > 1 == false is valid only for part of the interval. Is there any way to handle such situations?

Btw, what about producing a warning?

[lucaferranti]

Hehe this is an oldie but goldie.

Currently, IA uses the definition recommended by the IEEE standard for interval arithmetic, that is

$$ x < y \leftrightarrow \underline{x} < \underline{y} \land \overline{x} < \overline{y} $$

this definition stems from the work of Ulrich Kulisch about axiomatization of interval arithmetic, but in practice is not that useful and indeed it seems to unify developers in disappointment.

The 1.0-dev branch which, as the name suggests, is for the upcoming 1.0 version of IntervalArithmetic.jl, has a configuration mechanism that allows the user to choose how to handle boolean operations. Incidentally, I had a PR related to that yesterday, so you can have a sneak preview of that here.

As a final small comment, in the specific example you posted I wouldn't say the result in incorrect because it still holds $[0.9, 1] \subseteq [0.899999, 1.100001]$, that is the function returned a valid enclosure of the true range, which is the only promise of interval arithmetic. Especially for more complex functions, you should expect some overestimation (this is known as dependecy problem).

To play devil's advocate, here is a function with branching where a valid enclosure is not returned, so your point stands :)

julia> f(x)= x > 0 ? x + 3 : x - 3
f (generic function with 1 method)

julia> f(-1..1)
[-4, -2] # the true range is [-4, 4]

[lucaferranti]

alternatively, you may want to have a look at NumberIntervals.jl which uses a ternary logic for boolean intervals, returning missing for ambiguous cases (this will also be an option in 1.0-dev). Using that, however, your code would simply error for ambiguous cases, so depending on your application that may or may not be useful

[petvana]

Thank you so much for the rapid response. It seems possible to handle conditions manually, but it would lead to quite have machinery to handle it automatically, similarly like in ForwardDiff. Example of a manual solution follows. Would it be worthy to document this behavior?

julia> using IntervalArithmetic

julia> X = 0.9..1.1
[0.899999, 1.10001]

julia> f(x) = x > 1 ? x-1 : x;

julia> f(X)
[0.899999, 1.10001] # incorrect

julia> function f(x::Interval) 
           if x > 1 
               return x-1
           elseif x < 1 
               return x
           else # 1 ∈ x
               union(max(1,x)-1, min(1,x))
           end
       end;

julia> f(X)
[0, 1] # correct, but ...

Further, branching seems to be related to mod (#129) as the most specific output would be [0, 0.100001] ∪ [0.899999, 1].

As a final small comment, in the specific example you posted I wouldn't say the result in incorrect because it still holds

Unfortunately not, f(1.1) = 0.1 ∉ [0.899999, 1.10001].

Finally, thanks for guiding me to NumberIntervals.jl, as they introduced Indeterminate struct for this case in gwater/NumberIntervals.jl#12 .

[lucaferranti]

Would it be worthy to document this behavior?

absolutely!

Unfortunately not, f(1.1) = 0.1 ∉ [0.899999, 1.10001].

🤦‍♂️ I should really learn not to check github before my morning coffee... :)

Further, branching seems to be related to mod (#129) as the most specific output would be [0, 0.100001] ∪ [0.899999, 1].

Yeah that's a beast on itself. Generally, interval functions are defined via interval extension, that is the interval hull of the image. However, that is not always useful, e.g. for the set difference or the mod function as you point out. This paper describes the ideas of using union of intervals and IntervalUnionArithmetic.jl by @AnderGray has an implementation for it.

[AnderGray]

Branching is a tricky one with intervals.

... output would be [0, 0.100001] ∪ [0.899999, 1].

I agree that this should be the answer. If you were to do a Monte Carlo simulation with a distribution in X that union would contain the range of f(X). Automating this is hard, but theoretically possible with metaprogramming. Maybe at the intermediate representation (IR) level which gives you the branch structure:

julia> using IntervalArithmetic, IRTools

julia> f(x) = x > 1 ? x-1 : x;
julia> X = 0.9..1.1

julia> @code_ir f(X)

1: (%1, %2)
  %3 = %2 > 1
  br 3 unless %3
2:
  %4 = %2 - 1
  return %4
3:
  return %2

btw, if you mince, you do start getting an approximation:

julia> hull(f.(mince(X, 1000)))
[0.000199999, 1.0002]

[petvana]

It seems possible to introduce a macro that transforms the basic example. I drafted some code, but I'm not an expert for macros.

y = @interval x > 1 ? x-1 : x

or

y = @interval x > 1 ? f(x) : h(x)

in general, where function f and h must have no side-effect on local variables, i.e., be pure.

https://gist.github.com/petvana/9e583dcc6471368322afec63c0fa81e3

[lbenet]

It seems possible to introduce a macro that transforms the basic example. I drafted some code, but I'm not an expert for macros.

y = @interval x > 1 ? x-1 : x

or

y = @interval x > 1 ? f(x) : h(x)

in general, where function f and h must have no side-effect on local variables, i.e., be pure.

https://gist.github.com/petvana/9e583dcc6471368322afec63c0fa81e3

Could this be (eventually) added to the existing @interval macro? Or do we need a new one?

[petvana]

I was not aware of existing @interval. It is possible to combine, but may be confusing. What about @intervalif?

[lucaferranti]

Actually I think combining with @interval might make sense, The idea of the macro is that I I want to evaluate something like

julia> sin(0.1) + exp(0.3)*cos(0.2+log(0.523))

then to have a rigorous computation you should convert all numbers to intervals before evaluating any function. The @interval macro goes through the tree and replaces every number with the interval version (ok the following is a little messy, but you get the idea

julia> @macroexpand @interval sin(0.1) + exp(0.3)*cos(0.2+log(0.523))
:(IntervalArithmetic.Interval(sin(IntervalArithmetic.atomic(IntervalArithmetic.Interval{(IntervalArithmetic.parameters).precision_type}, 0.1)) + exp(IntervalArithmetic.atomic(IntervalArithmetic.Interval{(IntervalArithmetic.parameters).precision_type}, 0.3)) * cos(IntervalArithmetic.atomic(IntervalArithmetic.Interval{(IntervalArithmetic.parameters).precision_type}, 0.2) + log(IntervalArithmetic.atomic(IntervalArithmetic.Interval{(IntervalArithmetic.parameters).precision_type}, 0.523)))))

as such this could be seen as an enhancement of @interval to work with more complex expressions that have branches.

That being said, if integrating into the existing @interval is too much work, it's perfectly fine to start with a separate macro @intervalif that works only in front of if statements

[petvana]

But if some functions do not fully support intervals then it would be impossible to write @interval mod(a..b, pi) > 0.2 ? sin(a) : sin(b) because pi would be converted to an interval.

[lbenet]

But if some functions do not fully support intervals then it would be impossible to write @interval mod(a..b, pi) > 0.2 ? sin(a) : sin(b) because pi would be converted to an interval.

See this comment about mod...

[dpsanders]

One nice way to handle functions like this is with an implementation of piecewise functions, e.g. this one.

[OlivierHnt]

I like the code you linked.
Maybe we can come up with a small application that uses piecewise functions and include it in the "Examples" section of the docs.

Note, with the new semantics, the code you shared reads:

using IntervalArithmetic
using IntervalArithmetic.Symbols # to use the constructor `..`

piecewise(pieces::Tuple, x::Interval) = reduce(hull, f(intersect_interval(x, region)) for (region, f) ∈ pieces)

f(x) = piecewise((0..3 => x -> x + interval(1), 3..6 => x -> x, 6..Inf => x -> x + interval(2)), x)

f(2..5)

Also, this gives us a good opportunity to highlight the fact that the decoration is trv (i.e. trivial) anytime intersect_interval is called.