** Vectorized piecewise functions
:PROPERTIES:
:categories: math
:date: 2013/02/23 09:00:00
:updated: 2013/02/27 14:51:57
:END:
[[http://matlab.cheme.cmu.edu/2011/11/05/vectorized-piecewise-functions/][Matlab post]]
Occasionally we need to define piecewise functions, e.g.
\begin{eqnarray}
f(x) &=& 0, x < 0 \\
&=& x, 0 <= x < 1\\
&=& 2 - x, 1 < x <= 2\\
&=& 0, x > 2
\end{eqnarray}
Today we examine a few ways to define a function like this. A simple way is to use conditional statements.
#+BEGIN_SRC python :session
def f1(x):
if x < 0:
return 0
elif (x >= 0) & (x < 1):
return x
elif (x >= 1) & (x < 2):
return 2.0 - x
else:
return 0
print f1(-1)
print f1([0, 1, 2, 3]) # does not work!
#+END_SRC
#+RESULTS:
:
: ... ... ... ... ... ... ... ... >>> 0
: 0
This works, but the function is not vectorized, i.e. f([-1 0 2 3]) does not evaluate properly (it should give a list or array). You can get vectorized behavior by using list comprehension, or by writing your own loop. This does not fix all limitations, for example you cannot use the f1 function in the quad function to integrate it.
#+BEGIN_SRC python :session
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-1, 3)
y = [f1(xx) for xx in x]
plt.plot(x, y)
plt.savefig('images/vector-piecewise.png')
#+END_SRC
#+RESULTS:
:
: >>> >>> >>> >>> >>> []
[[./images/vector-piecewise.png]]
Neither of those methods is convenient. It would be nicer if the function was vectorized, which would allow the direct notation f1([0, 1, 2, 3, 4]). A simple way to achieve this is through the use of logical arrays. We create logical arrays from comparison statements.
#+BEGIN_SRC python :session
def f2(x):
'fully vectorized version'
x = np.asarray(x)
y = np.zeros(x.shape)
y += ((x >= 0) & (x < 1)) * x
y += ((x >= 1) & (x < 2)) * (2 - x)
return y
print f2([-1, 0, 1, 2, 3, 4])
x = np.linspace(-1,3);
plt.plot(x,f2(x))
plt.savefig('images/vector-piecewise-2.png')
#+END_SRC
#+RESULTS:
:
: ... ... ... ... ... ... >>> [ 0. 0. 1. 0. 0. 0.]
: >>> []
[[./images/vector-piecewise-2.png]]
A third approach is to use Heaviside functions. The Heaviside function is defined to be zero for x less than some value, and 0.5 for x=0, and 1 for x >= 0. If you can live with y=0.5 for x=0, you can define a vectorized function in terms of Heaviside functions like this.
#+BEGIN_SRC python :session
def heaviside(x):
x = np.array(x)
if x.shape != ():
y = np.zeros(x.shape)
y[x > 0.0] = 1
y[x == 0.0] = 0.5
else: # special case for 0d array (a number)
if x > 0: y = 1
elif x == 0: y = 0.5
else: y = 0
return y
def f3(x):
x = np.array(x)
y1 = (heaviside(x) - heaviside(x - 1)) * x # first interval
y2 = (heaviside(x - 1) - heaviside(x - 2)) * (2 - x) # second interval
return y1 + y2
from scipy.integrate import quad
print quad(f3, -1, 3)
#+END_SRC
#+RESULTS:
:
: ... ... ... ... ... ... ... ... ... ... >>> ... ... ... ... ... >>> >>> (1.0, 1.1102230246251565e-14)
#+BEGIN_SRC python :session
plt.plot(x, f3(x))
plt.savefig('images/vector-piecewise-3.png')
#+END_SRC
#+RESULTS:
: []
[[./images/vector-piecewise-3.png]]
There are many ways to define piecewise functions, and vectorization is not always necessary. The advantages of vectorization are usually notational simplicity and speed; loops in python are usually very slow compared to vectorized functions.