Implement piecewise function #689

Kolaru · 2024-11-28T21:41:29Z

Implement piecewise function, as described in #655 (and handling decorations).

For example, the abs function would be defined as

myabs = Piecewise(
    interval(-Inf, 0) => x -> -x,
    interval(0, Inf) => identity ;
    continuous = true
)

It also introduce the Constant constructor, as a workaround for Returns from Base: the later can not convert its output to interval when the input is an interval, which causes problems.

For example:

julia> f = Returns(1.0)
Returns{Float64}(1.0)

julia> f(22.2)
1.0

julia> f(interval(22.2))
1.0

julia> g = Constant(1.0)
Constant{Float64}(1.0)

julia> g(22.2)
1.0

julia> g(interval(22.2))
[1.0, 1.0]_com

It can be used for example to define a step function as

step = Piecewise(
    interval(-Inf, 0) => Constant(0),
    interval(0, Inf) => Constant(1)
)

OlivierHnt · 2024-11-28T23:02:57Z

Oh nice, looks great!

One thing is that there should maybe be a check to prevent overlapping domains. Eg the two following functions do not behave the same, and are essentially ill-defined

julia> foo = Piecewise(
                  interval(-Inf, 0) => x -> -x,
                  interval(-1, 1) => x -> 2x,
                  interval(0, Inf) => x -> x ;
                  continuous = true
              )

julia> goo = Piecewise(
                  interval(-Inf, 0) => x -> -x,
                  interval(0, Inf) => x -> x,
                  interval(-1, 1) => x -> 2x ;
                  continuous = true
              )

Moreover, domain should probably not be a convex hull, maybe a vector of interval (or a union type, related to #424 in some sense)? I also wonder if the resulting interval need to always have the trv decoration. Eg

Julia> julia> f = Piecewise(
                         interval(-1, 0) => identity, 
                         interval(1, 2) => identity);

julia> domain(f)
[-1.0, 2.0]_trv

Julia> julia> g = Piecewise(
                         interval(-1, 0) => identity, 
                         interval(0, 1) => identity);

julia> domain(g) # should it have the decoration `trv`?
[-1.0, 2.0]_trv

One question on the implementation. Currently you can only indicate if a Piecewise function is continuous or not. This may be mitigated by creating nested Piecewise functions... but what would you think about

struct Piecewise2{N,T<:NTuple{N}}
    pieces :: T 
    continuous :: NTuple{N,Bool}
end

Piecewise2(pieces::NTuple{N}, continuous::Bool) where {N} = Piecewise2(pieces, ntuple(_ -> continuous, Val(N)))

This way you do not hit a def decoration when evaluating the piecewise function at continuous points.

dpsanders · 2024-11-28T23:57:09Z

src/piecewise.jl

+
+function (piecewise::Piecewise)(x::Real)
+    for (region, f) in piecewise.pieces
+        in_interval(x, region) && return f(x)


This doesn't seem completely right: you need to return the union of f(x) over all the pieces that x overlaps with, I think.

e.g. the line

return reduce(hull, f(intersect_interval(X, region)) for (region, f) ∈ piecewise.pieces)

in my original implementation in the issue you limked to.

This part is the implementation for non-interval real only, so I think it is fine, as it would be a bit awkward to return an interval for non-interval input.

As @OlivierHnt mentioned, this implementation allows ill-defined functions, but I don't think the answer is the interval-ification of everything.

Kolaru · 2024-11-30T01:32:24Z

One thing is that there should maybe be a check to prevent overlapping domains. Eg the two following functions do not behave the same, and are essentially ill-defined

Yes, that reasonable.

To do that properly, we would need open and closed intervals (otherwise the boundary points are always included, and thus points on the boundary of pieces are always duplicated).

I am wondering if it would be worth to pull intervals from Intervals.jl to do this.

Moreover, domain should probably not be a convex hull, maybe a vector of interval (or a union type, related to #424 in some sense)? I also wonder if the resulting interval need to always have the trv decoration. Eg

Actually yes, here you would have a bug if there is a hole in the middle of the pieces.

For the decoration, I don't know. Those intervals, really have not much in common with the arithmetic intervals...

This way you do not hit a def decoration when evaluating the piecewise function at continuous points.

I thought about this a bit, I am not strongly against implementing it. It is just a bit awkward that the order of the pieces matter in the input.

OlivierHnt · 2024-11-30T10:48:10Z

Ah I did not notice the behaviour when evaluating at the boundary point. So it is also order dependent

julia> foo = Piecewise(
           interval(-Inf, 0) => x -> -x,
           interval(0, Inf) => x -> 3 ;
           continuous = false
       );

julia> foo(0)
0

julia> goo = Piecewise(
           interval(0, Inf) => x -> 3,
           interval(-Inf, 0) => x -> -x;
           continuous = false
       );

julia> goo(0)
3

My first impulse would be not to add a dependency on Intervals.jl. One idea could be to introduce yet an other option to specify the behaviour. Eg Piecewise(...; continuous = (false, boundary_behaviour)), where boundary_behaviour::Symbol could be :left, :right, :average (maybe the latter ought to be the default).

Actually yes, here you would have a bug if there is a hole in the middle of the pieces.

Right, so for f = Piecewise(interval(0, 1) => x -> 1, interval(10, 11) => x -> 2), I would expect domain(f) to return something in the vein of (interval(0, 1, trv), interval(10, 11, trv)).
EDIT: Ah but this is already implemented as subdomains. Now this feels even more confusing to me. Also, the error is different whether the input is an Interval or a Real

julia> Piecewise(interval(0, 1) => x -> 1, interval(10, 11) => x -> 2)(interval(5.))
ERROR: ArgumentError: reducing over an empty collection is not allowed; consider supplying `init` to the reducer

Maybe I am overlooking something, but I think that the behaviour for not being in the domain should be the same regardless whether you are in the convex hull or not. Eg compare the above error with

julia> Piecewise(interval(0, 1) => x -> 1, interval(10, 11) => x -> 2)(interval(-2))
∅_trv

For the decoration, I don't know. Those intervals, really have not much in common with the arithmetic intervals...

Yes you are right, maybe trv is appropriate after all.

I thought about this a bit, I am not strongly against implementing it. It is just a bit awkward that the order of the pieces matter in the input.

Indeed.

OlivierHnt · 2024-11-30T18:52:54Z

Mhm I gave it a bit more thoughts. Here is an attempt (order independent)

struct Piecewise{T<:Tuple}
    pieces :: T
    # The symbol indicates whether or not the boudnary is included in the interval:
    # Symbol(11) means that the interval contains both extremities
    # Symbol(10) means that the interval contains only at the left extremity
    # Symbol(01) means that the interval contains only at the right extremity
    # Symbol(00) means that the interval does not contain either extremity
    # TODO: need to normalize `pieces`
    # TODO: missing `continuous` (to determine the decoration `def`, or `com`)
end

Piecewise(pieces...) = Piecewise(pieces)

_subdomains(piecewise::Piecewise) = first.(piecewise.pieces)

function _intersect_subdomain(x::Interval, region::Tuple{Interval,Symbol})
    y = IntervalArithmetic.setdecoration(intersect_interval(x, region[1]), decoration(x))
    if isempty_interval(y)
        return y
    elseif region[2] === Symbol(11)
        return y
    elseif region[2] === Symbol(01)
        isthin(y, inf(region[1])) && return emptyinterval(y)
        return y
    elseif region[2] === Symbol(10)
        isthin(y, sup(region[1])) && return emptyinterval(y)
        return y
    else # region[2] === Symbol(00)
        isthin(y, inf(region[1])) | isthin(y, sup(region[1])) && return emptyinterval(y)
        return y
    end
end

function (piecewise::Piecewise)(x::Interval)
    subdoms = _subdomains(piecewise)
    vals = _intersect_subdomain.(x, subdoms)
    all(v -> isempty_interval(v[1]), vals) && return throw(DomainError("piecewise function was called with $x which is outside of its domain $subdoms"))
    return _evaluate(piecewise.pieces, vals, nothing)
end

function _evaluate(p, x, r) # maybe not type stable
    isempty_interval(x[1][1]) && return _evaluate(Base.tail(p), Base.tail(x), r)
    fx = last(p[1])(x[1][1])
    return _evaluate(Base.tail(p), Base.tail(x), IntervalArithmetic.setdecoration(hull(fx, r), min(decoration(fx), decoration(r))))
end

function _evaluate(p, x, ::Nothing) # maybe not type stable
    isempty_interval(x[1][1]) && return _evaluate(Base.tail(p), Base.tail(x), nothing)
    fx = last(p[1])(x[1][1])
    return _evaluate(Base.tail(p), Base.tail(x), fx)
end

function _evaluate(p::Tuple{Pair}, x, r) # maybe not type stable
    isempty_interval(x[1][1]) && return r
    fx = last(p[1])(x[1][1])
    return IntervalArithmetic.setdecoration(hull(fx, r), min(decoration(fx), decoration(r)))
end

function _evaluate(p::Tuple{Pair}, x, ::Nothing) # maybe not type stable
    isempty_interval(x[1][1]) && return x[1][1]
    return last(p[1])(x[1][1])
end

then

foo = Piecewise((interval(-Inf, 0), Symbol(11)) => x -> -x,
                (interval( 0, Inf), Symbol(01)) => x -> interval(3));

foo(interval(-1, 0))

foo(interval(-1, 1)) # decoration is wrong, currently not implemented

foo(interval(0, 1))

foo(interval(0))

goo = Piecewise((interval(0, 1), Symbol(00)) => x -> interval(1),
                (interval(2, 3), Symbol(01)) => x -> interval(3));

goo(interval(0)) # now it throws

goo(interval(1, 2)) # now it throws

goo(interval(2, 3))

Kolaru · 2024-12-05T15:57:17Z

I will try to implement the ForwardDiff extension for this, it may give some directions as to how to to handle the boundary points.

Kolaru added 3 commits November 18, 2024 23:29

Prototype

7204462

Add basic tests and Constant

2268e96

Add tests for and correct out-of-domain behavior

35638a8

dpsanders reviewed Nov 28, 2024

View reviewed changes

OlivierHnt linked an issue Dec 12, 2024 that may be closed by this pull request

Piecewise functions #655

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement piecewise function #689

Implement piecewise function #689

Kolaru commented Nov 28, 2024

OlivierHnt commented Nov 28, 2024

dpsanders Nov 28, 2024 •

edited

Loading

dpsanders Nov 29, 2024

Kolaru Nov 30, 2024

Kolaru commented Nov 30, 2024 •

edited

Loading

OlivierHnt commented Nov 30, 2024 •

edited

Loading

OlivierHnt commented Nov 30, 2024

Kolaru commented Dec 5, 2024

Implement piecewise function #689

Are you sure you want to change the base?

Implement piecewise function #689

Conversation

Kolaru commented Nov 28, 2024

OlivierHnt commented Nov 28, 2024

dpsanders Nov 28, 2024 • edited Loading

Choose a reason for hiding this comment

dpsanders Nov 29, 2024

Choose a reason for hiding this comment

Kolaru Nov 30, 2024

Choose a reason for hiding this comment

Kolaru commented Nov 30, 2024 • edited Loading

OlivierHnt commented Nov 30, 2024 • edited Loading

OlivierHnt commented Nov 30, 2024

Kolaru commented Dec 5, 2024

dpsanders Nov 28, 2024 •

edited

Loading

Kolaru commented Nov 30, 2024 •

edited

Loading

OlivierHnt commented Nov 30, 2024 •

edited

Loading