Functions with branches

[baggepinnen]

Hello! I just tried this package to see if it could be useful for automatic stability proofs of nonlinear feedback systems using the circle criterion. I'm not sure if I'm using it wrong, but I get an infinite enclosure for a simple function (relative saturation)

julia> using RangeEnclosures

julia> domain = RangeEnclosures.Interval(-3, 3)
[-3, 3]

julia> fun = x->clamp(x, -1, 1)/x
#105 (generic function with 1 method)

julia> plot(fun, domain.lo, domain.hi)

julia> enclose(fun, domain, BranchAndBoundEnclosure())
[-∞, ∞]

[lucaferranti]

Hi there,

if you look at the definition of clamp in Julia source code, you will notice that it uses if-else branches, in general that kind of code does not play well with interval arithmetic. For the clamp part, a work-around could be to rewrite it as follows (here I created a new function to highlight the comparison, in real life you probably want to overload Base.clamp)

julia> myclamp(x::Interval, lo, hi) = max(min(x, hi), lo)
myclamp (generic function with 1 method)

julia> clamp(-4..4, 2, 3)
[-4, 4]

julia> myclamp(-4..4, 2, 3)
[2, 3]

That being said, I am not convinced by your result, the function clamp(x, -1, 1)/x is not defined at 0,

julia> fun(x) = clamp(x, -1, 1)/x
fun (generic function with 1 method)

julia> fun(0)
NaN

and if x contains zero, also clamp(x, -1, 1) will contain zero (regardless of what version you use). As a consequence, the ratio will return [-Inf, Inf] and the returned enclosure seems expected to me. Could you clarify why you would expect that result?

The nice plot is because plot(f, -3, 3) doesn't happen to evaluate the function at 0 using plot(-3:0.1:3, fun.(-3:0.1:3)) I get the following picture

[baggepinnen]

Hello and thanks for your answer :)

Is there any reason clamp is not properly defined for intervals? This result seems very easy to improve upon

julia> clamp(-4..4, 2, 3)
[-4, 4]

More generally, this

julia> using IntervalArithmetic

julia> 1..3 > 2
false

seems very dangerous and non-rigorous, if this appears in

if a > 2
    something
else
    something else
end

the code will just pick one branch and return something but the inequality could equally well be true! The code below should most definitely throw an error

julia> 1..3 < 2 ? 1 : 2
2

since the correct and rigorous answer is 1..2, which I see no easy way to reach using operator overloading.

I realized that the example was unfortunate due to the singularity at 0. This particular singularity happens to be removable, but the julia function obviously doesn't know this. Is there any nice strategy to handle removable singularities other than splitting up the domain into two, one for each side of the singularity?

[mforets]

Thanks for trying the package @baggepinnen. I think the short answer is that the library doesn't yet have proper support for functions with branches (I suggest we add a comment about this in the docs). That said, a safer alternative in those cases may be to pass an interval from NumberIntervals (I didn't try if it works out of the box), or wait for a feature in IntervalArithmetic v1.0 that will let you choose the behavior against comparisons iirc.