Package to work with piecewise functions?

General Usage
1

Hi, @aerappa and I are looking for a package that would allow us to easily represent and work with (mathematical) functions that are defined in a piecewise manner.

To be a bit more specific, assuming we are working with a function such as

f(x) = \begin{cases} e^x & \text{if } x\leq 0 \\ x+1 & \text{otherwise} \end{cases}

we’d like to be able to efficiently represent it as a Julia object and evaluate it or its derivatives at any point using code that could ideally look like (don’t mind the API; I just made it up for the purpose of this post in order to clarify what we need):

f  = PieceWise(exp, 0.0,
               Polynomial([1,1]))
f′ = derivative(f)

f(-1.0), f′(-1.0)  # == exp(-1), exp(-1)
f( 1.0), f′( 1.0)  # == 2, 1

On the one hand, this seems easy enough to implement on our own. On the other hand, it looks general enough that someone might have registered a package to do this (and why reinvent the wheel?). But I haven’t found anything like this (yet!). Does anyone here know of an existing package that would do something this?

1 Like
2 
julia> using ForwardDiff

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

julia> f′(x) = ForwardDiff.derivative(f, x)
f′ (generic function with 1 method)

julia> f(-1.0), f′(-1.0)  # == exp(-1), exp(-1)
(0.36787944117144233, 0.36787944117144233)

julia> f( 1.0), f′( 1.0)  # == 2, 1
(2.0, 1.0)

ApproxFun.jl is another way to work with piecewise functions via piecewise polynomial approximations (see its PiecewiseSegment type).

4 Likes
3

Thanks for your reply! In our use case, we would like to use Symbolics.jl to take derivatives rather than ForwardDiff.jl. In particular, the issue I keep running into is the following issue:

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

julia> using Symbolics

julia> @syms x
(x,)

julia> f(x)
ERROR: MethodError: no method matching isless(::SymbolicUtils.BasicSymbolic{Number}, ::Int64)
Closest candidates are:
  isless(::Union{StatsBase.PValue, StatsBase.TestStat}, ::Real) at ~/.julia/packages/StatsBase/XgjIN/src/statmodels.jl:90
  isless(::ColorTypes.AbstractGray, ::Real) at ~/.julia/packages/ColorTypes/1dGw6/src/operations.jl:38
  isless(::SymbolicUtils.Symbolic{<:Real}, ::Real) at ~/.julia/packages/SymbolicUtils/Vzo2s/src/methods.jl:173
  ...
Stacktrace:
 [1] <(x::SymbolicUtils.BasicSymbolic{Number}, y::Int64)
   @ Base ./operators.jl:356
 [2] <=(x::SymbolicUtils.BasicSymbolic{Number}, y::Int64)
   @ Base ./operators.jl:405
 [3] f(x::SymbolicUtils.BasicSymbolic{Number})
   @ Main ./REPL[9]:1
 [4] top-level scope
   @ REPL[12]:1

julia>

Edit: After having a look at this thread, the following MWE seems pretty good for my needs

julia> f(x) = ifelse(x ≤ 0.0, exp(2*x), x^2+1)
f (generic function with 1 method)

julia> using Symbolics

julia> @syms x::Real
(x,)

julia> f(x)
ifelse(x <= 0.0, exp(2x), 1 + x^2)

julia> Symbolics.derivative(f(x), x)
ifelse(x <= 0.0, 2exp(2x), 2x)

Eventually, I also want to support multiple intervals, so I came up with the (admittedly clunky) solution for the case of two intervals

julia> ifelse2(c1, c2, x, y, z) = ifelse(c1, x, ifelse(!c2, y, z))
ifelse2 (generic function with 1 method)

julia> f2(x) = ifelse2(x ≤ 2, x > 4, x^2, exp(2*x), 4*x)
f2 (generic function with 1 method)

julia> Symbolics.derivative(f2(x), x)
ifelse(x <= 2, 2x, ifelse(!(x > 4), 2exp(2x), 4))