Interval arithmetic through functions with branches

[rdeits]

(attempting to summarize our discussion from today).

The particular issue that I was referring to is this:

Suppose we have a function f(x):

function f(x)
  if x > 2
    2x
  else
    -2x
  end
end

We can evaluate this on an interval:

julia> f(3)
6

julia> f(4)
8

julia> f(3..4)
[6, 8]

julia> f(0)
0

julia> f(1)
-2

julia> f(0..1)
[-2, -0]

This gives a correct answer as long as the interval does not include the value 2. When the interval spans 2, we get the wrong answer:

julia> f(1)
-2

julia> f(3)
6

julia> f(1..3)
[-6, -2]

That's because 1..3 > 2 just evaluates to false, so the code path which is taken is the one in which f(x) = -2x across the entire interval.

Given that I just learned about this entire field today, I don't really think I can suggest a "correct" behavior, but this seems like a case where it would be easy for a naive user like myself to get a wrong answer without realizing it.

[lbenet]

Essentially you have to create a method of f for Interval that distinguishes the domains of continuity of the function, and what to do when the discontinuity is inside the interval you consider.

Something like this could be helpful, which uses your definition of f:

function f(x::Interval)
    x.lo > 2 && return Interval(f(x.lo), f(x.hi))
    x.hi ≤ 2 && return Interval(f(x.hi), f(x.lo))
    return Interval( f(2), max(f(x.lo),f(x.hi)))
end

Now, using your examples:

julia> f(0..1)
[-2, -0]

julia> f(3..4)
[6, 8]

julia> f(0..1)
[-2, -0]

julia> f(1..3)
[-4, 6]

julia> f(2..2)
[-4, -4]

For this kind of things, it is worth taking a look on how tan(a::Interval) is coded.

[dpsanders]

It would be nice, but seems hard, to find a way to automatize the generation of the interval version with the correct behaviour.

[dpsanders]

If you write the function as

f(x) = 2x * sign(x - 2)

then it has the correct behaviour, i.e. gives an inclusion of the result:

julia> f(x) = 2x * sign(x - 2)
f (generic function with 1 method)

julia> f(3..4)
[6, 8]

julia> f(-1..1)
[-2, 2]

julia> f(-1)
2

julia> f(1)
-2

julia> f(-1..4)
[-8, 8]

[dpsanders]

You can also use decorated functions. (@rdeits: Decorations are flags that give information about what happened during a calculation.)

julia> displaymode(decorations=true)

julia> y = @decorated(3..4)
[3, 4]_com

julia> sign(y)
[1, 1]_com

julia> f(y)
[6, 8]_com

The com decoration here shows that no discontinuity was encountered during the evaluation.

julia> f(@decorated(3..4))
[6, 8]_com

julia> f(@decorated(-1..1))
[-2, 2]_com

julia> f(@decorated(-1..3))
[-6, 6]_def

The def decoration here shows that a discontinuity was encountered at some point during the evaluation of the history of the calculation by which the interval [-6,6] was calculated.

[dpsanders]

So in principle it may be possible to write a macro to rewrite simple if statements like this to use sign / abs etc.

[dpsanders]

Note that f(@decorated(-1..3)) gives a strict over-estimate of the range. In this particular case, we know that f is monotone, so we could give an exact range using the endpoints. cf #198