Incorrect root of function involving `max`

[perrutquist]

julia> roots(x -> max(x,-x)+1, -2..2)
1-element Array{Root{Interval{Float64}},1}:
 Root([1, 1], :unique)

This seems like an incorrect result. The function (equivalent to abs(x)+1) does not have any roots.

Providing a derivative explicitly: roots(x -> max(x,-x)+1, sign, -2..2) or changing the contractor: roots(x -> max(x,-x)+1, -2..2, Bisection) makes the problem go away.

Maybe related to #98 ?

[Kolaru]

It has to do with the fact that the function is not differentiable, which is an implicit requirement for Newton (the default) and Krawczyk methods. I don't know what ForwardDiff.jl does in such cases (it is used if the derivative is not explicitly given), but even if it can find the correct derivative, IntervalArithmetic.jl has currently no way of dealing with arbitrary piecewise continuous functions, so this failure is expected.

We should probably document that, with a big warning sign.

[dpsanders]

IntervalArithmetic.jl is I think not the problem; rather, it is the lack of information about derivatives of the functions. This should be solvable with the new ChainRules.jl and ForwardDiff2.jl packages.

[Kolaru]

For a function like f(x) = max(x, -3x), the derivative would be

f'(x < 0) = -3
f'(x >= 0) = 1

As far as I can tell, we don't currently have a mechanism to conclude that e.g. f'(-2 .. 1) = (-3 .. 1), a problem that remains even if we fully know the derivative. The example with sign being the derivative works because the sign function is explicitly define in IntervlArithmetic.jl.

ChainRules.jl and ForwardDiff2.jl looks very promising though, I'll check them in more details.

[dpsanders]

Here I showed how to define a general piecewise function. Copying (with a slight improvement) here for simplicity:

"Evaluate interval function with piecewise region--function pairs"
function piecewise(pieces::Vector{<:Pair}, X)
    return reduce(∪, f(X ∩ region) for (region, f) in pieces)
end

piecewise(pieces) = X -> piecewise(pieces, X)

X = 2..5
f2 = piecewise([0..3 => x->x+1,
                3..6 => x->x,
                6..∞ => x->x+2])

f(x) = sum(f2 ∘ abs, x)

EDIT: We should probably just export this function from IntervalArithmetic.jl.

[dpsanders]

If we introduce a Piecewise type then we should be able to define ForwardDiff.derivative for objects of that type.

[dpsanders]

Hmm, making this efficient could require FunctionWrappers.jl.

[dpsanders]

Something like this, ~but it is currently not working correctly~ FIXED.

Note that ideally we would somehow generate and store only once the derivative; currently it is re-constructed each time:

using IntervalArithmetic

struct PiecewiseFunction{P}
    pieces::P
end

"Evaluate f over the region X if X is non-empty"
function evaluate(f, X)
    isempty(X) && return X
    return f(X)
end

function (f::PiecewiseFunction)(X::Interval)
    return reduce(∪, evaluate(f, X ∩ region) for (region, f) in f.pieces)
end

pieces = (0..3 => x->x+1,
 3..6 => x->x,
 6..∞ => x->x+2)

f = PiecewiseFunction(pieces)

X = 2..5
f(X)

using ForwardDiff
import ForwardDiff: derivative

ForwardDiff.derivative(f) = x -> ForwardDiff.derivative(f, x)

ForwardDiff.derivative(f::PiecewiseFunction) =
    PiecewiseFunction( Tuple(region => ForwardDiff.derivative(g) for (region, g) in f.pieces))

ForwardDiff.derivative(f::PiecewiseFunction, x::Interval) = ForwardDiff.derivative(f)(x)

myabs = PiecewiseFunction( (
                        (-Inf..0) => x -> -x,
                        (0..Inf) => x -> x ) )

ForwardDiff.derivative(myabs, -1..1)
ForwardDiff.derivative(myabs, 2..3)

myabs( (-1..1) ∩ (2..3))

myabsderiv = ForwardDiff.derivative(myabs)

myabsderiv(-1..1)   # -1..1
myabsderiv(2..3)   # 1..1
myabsderiv(-3.. -2)   # -1.. -1