10

I have an application in which I need to define a piecewise function, IE, f(x) = g(x) for [x in some range], f(x)=h(x) for [x in some other range], ... etc.

Is there a nice way to do this in Julia? I'd rather not use if-else because it seems that I'd have to check every range for large values of x. The way that I was thinking was to construct an array of functions and an array of bounds/ranges, then when f(x) is called, do a binary search on the ranges to find the appropriate index and use the corresponding function (IE, h(x), g(x), etc.

It seems as though such a mathematically friendly language might have some functionality for this, but the documentation doesn't mention piecewise in this manner. Hopefully someone else has given this some thought, thanks!

user3587051
  • 473
  • 6
  • 13
  • You might want to look at the implementation of [NumPy's `piecewise` function](http://docs.scipy.org/doc/numpy/reference/generated/numpy.piecewise.html). – jub0bs Dec 27 '14 at 17:58

3 Answers3

5

with a Heaviside function you can do a interval function:

function heaviside(t)
   0.5 * (sign(t) + 1)
end

and

function interval(t, a, b)
   heaviside(t-a) - heaviside(t-b)
end

function piecewise(t)
   sinc(t) .* interval(t,-3,3) + cos(t) .* interval(t, 4,7)
end

and I think it could also implement a subtype Interval, it would be much more elegant

elsuizo
  • 144
  • 1
  • 3
2

I tried to implement a piecewise function for Julia, and this is the result:

function piecewise(x::Symbol,c::Expr,f::Expr)
  n=length(f.args)
  @assert n==length(c.args)
  @assert c.head==:vect
  @assert f.head==:vect
  vf=Vector{Function}(n)
  for i in 1:n
    vf[i]=@eval $x->$(f.args[i])
  end
  return @eval ($x)->($(vf)[findfirst($c)])($x) 
end
pf=piecewise(:x,:([x>0, x==0, x<0]),:([2*x,-1,-x]))
pf(1) # => 2
pf(-2) # => 2
pf(0) # => -1
Reza Afzalan
  • 5,646
  • 3
  • 26
  • 44
1

Why not something like this?

function piecewise(x::Float64, breakpts::Vector{Float64}, f::Vector{Function})
       @assert(issorted(breakpts))
       @assert(length(breakpts) == length(f)+1)
       b = searchsortedfirst(breakpts, x)
       return f[b](x)
end

piecewise(X::Vector{Float64}, bpts, f) = [ piecewise(x,bpts,f) for x in X ]

Here you have a list of (sorted) breakpoints, and you can use the optimized searchsortedfirst to find the first breakpoint b greater than x. The edge case when no breakpoint is greater than x is also handled appropriately since length(breakpts)+1 is returned, so b is the correct index into the vector of functions f.

digbyterrell
  • 3,449
  • 2
  • 24
  • 24