Package to work with piecewise functions?

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,
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
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).


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

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
 [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


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

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))