Joining together ordered sample points (xi, yi) when ths xi's are different yields a piecewise linear continuous (P1) function. This snippet defines a handy class to manipulate these functions. It allows
computing the value and the slope of the function at each point, arithmetic operations, absolute value, truncation and linear combinations of such functions, maximum, minimum and integral between bounds. If wxPython is installed on your system, the class also allows plotting the function in a wx plot widget.

``````#!/usr/bin/env python

# Gribouillis at www.daniweb.com

# module P1Function.py
# type "help(P1Function)" in a python shell for class documentation
"""A module implementing piecewise linear continuous functions.
"""
from bisect import bisect_left ,bisect_right
from collections import defaultdict

class P1Function (object ):
"""P1Function(seq) -> new piecewise linear continuous function object
'seq' is an iterable over pairs of numbers (x, y).
Raises ValueError if one of the numbers is not convertible to float
or when 2 points have the same x but not the same y.
"""
_bisect =(bisect_left ,bisect_right )
_eps =1.0e-9
def __init__ (self ,seq ):
pairs =sorted (set (seq ))
pairs =[(float (x ),float (y ))for x ,y in pairs ]
if len (pairs )<2 :
if pairs :
pairs .append ((1.0 +2 *abs (pairs [0 ][0 ]),pairs [0 ][1 ]))
else :
pairs =[(x ,0.0 )for x in self .float_range (0.0 ,1.0 ,2 )]
tmp =[pairs [0 ]]
for i in xrange (1 ,len (pairs )):
x ,y =tmp [-1 ]
if pairs [i ][0 ]==x :
if pairs [i ][1 ]!=y :
raise ValueError ("Sample points aligned vertically")
else :
tmp .append (pairs [i ])
pairs =[tmp [0 ]]
for i in xrange (1 ,len (tmp )-1 ):
a ,c ,b =pairs [-1 ][0 ],tmp [i ][0 ],tmp [i +1 ][0 ]
ya ,yc ,yb =pairs [-1 ][1 ],tmp [i ][1 ],tmp [i +1 ][1 ]
ycc =ya *((c -a )/(b -a ))+yb *((b -c )/(b -a ))
if abs (yc -ycc )>=self ._eps *max (abs (yc ),abs (ycc )):
pairs .append (tmp [i ])
pairs .append (tmp [-1 ])
self ._x ,self ._y =zip (*pairs )
self ._x =[float (t )for t in self ._x ]
def __iter__ (self ):
"""iter(func) -> iterator over pairs of numbers (x, y).
This allows copying by passing a P1Function to the class constructor
"""
return iter (zip (self ._x ,self ._y ))
def __len__ (self ):
"""len(func) -> number of pairs (x, y) contained in the P1Function.
This number may be different from the number of pairs passed to the
constructor.
"""
return len (self ._x )
def pair (self ,i ):
"func.pair(i) -> a sample point (xi, yi) (0 <= i < len(func))"
return self ._x [i ],self ._y [i ]
def x (self ,i ):
"func.x(i) -> the coordinate xi of the sample point at index i"
return self ._x [i ]
def y (self ,i ):
"func.y(i) -> the coordinate yi of the sample point at index i"
return self ._y [i ]
def index (self ,x ,left =False ):
"""func.index(x) -> an index i such that
i = 0    if x < func.x(1)
i = n-2  if func.x(n-2) <= x
func.x(i) <= x < func.x(i+1) otherwise.
Here, x is a float and n = len(func) is the number of pairs in func.
When left is True, <= and < are exchanged in these definitions.
"""

return min (len (self ._x )-2 ,
max (0 ,self ._bisect [not left ](self ._x ,x )-1 ))
def index_bounds (self ,a =None ,b =None ):
"""self->index_bounds(a, b) -> pair of indexes (i, j)
such that a <= self.x(k) <= b for all k in [i, j[.
Returns a pair (i, i) when there is no such k.
a and b default to None, meaning -infinity and +infinity"""
a =a if a is not None else self ._x [0 ]
b =b if b is not None else self ._x [-1 ]
i =self ._bisect [0 ](self ._x ,a )# bisect left
j =self ._bisect [1 ](self ._x ,b )# bisect_right
return (i ,j )if i <j else (i ,i )
def __call__ (self ,x ):
"""func(x) -> the value of the piecewise linear func at point x.
x can be any float number, the function is extended by straight lines
outside it's initial x-bounds"""
i =self .index (x )
return ((self ._x [i +1 ]-x )*self ._y [i ]+
(x -self ._x [i ])*self ._y [i +1 ])/(self ._x [i +1 ]-self ._x [i ])
@staticmethod
def float_range (start ,stop ,n ):
"""P1Function.float_range(start, stop, n) -> list of n float
A range function which returns n equally spaced floating numbers
from start to stop"""
L =[0.0 ]*n
nm1 =n -1
nm1inv =1.0 /nm1
start ,stop =start *nm1inv ,stop *nm1inv
for i in range (n ):
L [i ]=(start *(nm1 -i )+stop *i )
return L
def slope (self ,x ,left =False ):
"""func.slope(x) -> the slope of the function at x
For sample points, this is the slope at the right of the point, unless
left is True.
"""
i =self .index (x ,left )
return (self ._y [i +1 ]-self ._y [i ])/(self ._x [i +1 ]-self ._x [i ])
def max (self ,a =None ,b =None ):
"""func.max(a,b) -> the maximum value of func between a and b
a and b default to None, meaning -infinity and +infinity. In this case,
the maximum can be None, meaning +infinity."""
return self ._extreme (a ,b ,1 )
def min (self ,a =None ,b =None ):
"""func.min(a,b) -> the minimum value of func between a and b
a and b default to None, meaning -infinity and +infinity. In this case,
the minimum can be None, meaning -infinity."""
return self ._extreme (a ,b ,-1 )
def _extreme (self ,a ,b ,eps ):
"""func._extreme(a, b, eps) -> the max or min of func between a and b
depending on eps being 1 or -1"""
opt =(min ,None ,max )[1 +eps ]
if a is None :
if (self ._y [1 ]-self ._y [0 ])*eps <0 :
return None
a =self ._x [0 ]
if b is None :
if (self ._y [-1 ]-self ._y [-2 ])*eps >0 :
return None
b =self ._x [-1 ]
if a >b :
a ,b =b ,a
i ,j =self .index_bounds (a ,b )
m =opt (self (a ),self (b ))
return opt (m ,opt (self ._y [i :j ])if i <j else m )
def integral (self ,a =None ,b =None ):
"""func.integral(a, b) -> the value of the integral of func on [a,b].
a and b default to None, meaning -infinity and +infinity.
May return None when a or b is None, meaning undefined integral."""
if a is None :
if self ._y [0 ]!=0.0 or self ._y [1 ]!=0.0 :
return None
a =self ._x [1 ]
if b is None :
if self ._y [-1 ]!=0.0 or self ._y [-2 ]!=0.0 :
return None
b =self ._x [-2 ]
swap =False
if a >b :
if a ==b :
return 0.0
else :
swap =True
a ,b =b ,a
i ,j =self .index_bounds (a ,b )
j -=1
if i >j :
return ((a -b )if swap else (b -a ))*self (0.5 *(b +a ))
x ,y =self ._x ,self ._y
res =(x [i ]-a )*self (0.5 *(x [i ]+a ))+(b -x [j ])*self (0.5 *(b +x [j ]))
res +=sum (
(x [k +1 ]-x [k ])*(y [k +1 ]+y [k ])*0.5 for k in xrange (i ,j ))
return -res if swap else res
@staticmethod
def combine (pairs ):
"""P1Function.combine(pairs) -> a linear combination of P1Functions
pairs is a sequence (ai, fi) where ai is a number and fi a P1Function.
the sum of the ai * fi is returned as a P1Function"""
d =defaultdict (float )
for c ,f in pairs :
d [f ]+=c
for f ,c in d .items ():
if c ==0.0 :
del d [f ]
xs =sorted (set (t for f in d for t in f ._x ))
data =[(x ,sum (d [f ]*f (x )for f in d ))for x in xs ]
return P1Function (data )
def rescaled (self ,center =(0.0 ,0.0 ),zoom =(1.0 ,1.0 )):
"""func.rescaled(center, zoom) -> a rescaled P1Function.
center is a pair (cx, cy), zoom is a pair (zx, zy).
A P1Function g is returned, such that
g(zx * (x - cx)) = zy * (func(x) - cy)
This allows changing units in x or y (eg celsius to farenheit)"""

cx ,cy =center
zx ,zy =zoom
x ,y =self ._x ,self ._y
return P1Function (
(zx *(x [i ]-cx ),zy *(y [i ]-cy ))for i in xrange (len (x )))
def truncated (self ,a =None ,b =None ):
"""func.truncated(a, b) -> a new P1Function truncated to [a,b].
Truncation means forget the points (xi, yi) when xi is not in [a,b].
a and b default to None, meaning -infinity and +infinity"""
L =[(self ._x [i ],self ._y [i ])for i in self .index_bounds (a ,b )]
L .extend ((x ,self (x ))for x in (a ,b )if x is not None )
return P1Function (L )
def __abs__ (self ):
"""abs(func) -> the absolute value of func (as a P1Function)"""
L =[]
x ,y =self ._x ,self ._y
if (y [0 ]>0.0 and y [1 ]>y [0 ])or (
y [0 ]<0.0 and y [1 ]<y [0 ]):
u =y [0 ]*(x [0 ]-x [1 ])/(y [0 ]-y [1 ])
L .append ((x [0 ]-1.5 *u ,0.5 *abs (y [0 ])))
L .append ((x [0 ]-u ,0.0 ))
for i in range (len (x )-1 ):
L .append ((x [i ],abs (y [i ])))
if (y [i ]<0.0 and y [i +1 ]>0.0 )or (
y [i ]>0.0 and y [i +1 ]<0.0 ):
u =y [i ]*(x [i ]-x [i +1 ])/(y [i ]-y [i +1 ])
L .append ((x [i ]-u ,0.0 ))
L .append ((x [-1 ],abs (y [-1 ])))
if (y [-1 ]>0 and y [-2 ]>y [-1 ])or (
y [-1 ]<0 and y [-2 ]<y [-1 ]):
u =y [-1 ]*(x [-2 ]-x [-1 ])/(y [-1 ]-y [-2 ])
L .append ((x [-1 ]-u ,0.0 ))
L .append ((x [-1 ]-1.5 *u ,0.5 *abs (y [-1 ])))
return P1Function (L )
"""f1 + f2 -> the sum of 2 P1Functions as a P1Function"""
return self .combine ([(1.0 ,self ),(1.0 ,other )])
def __sub__ (self ,other ):
"""f1 - f2 -> the difference of 2 P1Functions as a P1Function"""
return self .combine ([(1.0 ,self ),(-1.0 ,other )])
def __neg__ (self ):
"""-func -> the opposite of func as a P1Function"""
return self .combine ([(-1.0 ,self )])
def __mul__ (self ,factor ):
"""a * func -> a P1Function, func muliplied by the number a"""
return self .combine ([float (factor ),self ])
def plot (self ,**kwd ):
"""func.plot(**options) -> subprocess.Popen object
Starts a new process which plots func using wxPython.
options (with default values) are:
title="P1Function", size=(400,300), line_colour='blue',
line_width=1, marker='plus', caption='',
x_axis='x axis', y_axis='y axis'
Possible marker values:
'circle' 'cross' 'square' 'dot' 'plus' 'triangle'"""
import subprocess
options =dict (self ._plot_defaults )
options .update (dict ((k .upper (),v )for k ,v in kwd .items ()))
options ["POINTS"]=list (self )
program =self ._plot_template %options
child =subprocess .Popen ('python -c "exec(eval(raw_input()))"',
shell =True ,stdin =subprocess .PIPE )
child .stdin .write (repr (program ))
child .stdin .write ("\n")
child .stdin .close ()
return child

_plot_defaults =dict (
TITLE ="P1 Function",
SIZE =(400 ,300 ),
LINE_COLOUR ='blue',
LINE_WIDTH =1 ,
MARKER ='plus',
CAPTION ='',
X_AXIS ='x axis',
Y_AXIS ='y axis'
)

_plot_template ="""
# thanks to vegaseat at www.daniweb.com
import wx
import wx.lib.plot as plot
data = %(POINTS)s

class MyFrame(object):
def __init__(self):
self.frame = wx.Frame( None,
title='%(TITLE)s', id=-1, size=%(SIZE)s)
self.panel = wx.Panel(self.frame)
plotter = plot.PlotCanvas(self.panel)
plotter.SetSize(%(SIZE)s)
plotter.SetEnableZoom(True)
line = plot.PolyLine(data,
colour='%(LINE_COLOUR)s', width=%(LINE_WIDTH)s)
marker = plot.PolyMarker(data, marker='%(MARKER)s')
gc = plot.PlotGraphics([line, marker],
'%(CAPTION)s', '%(X_AXIS)s', '%(Y_AXIS)s')
plotter.Draw(gc)
self.frame.Show(True)
app = wx.PySimpleApp()
f = MyFrame()
app.MainLoop()
"""``````
1
Contributor
1
6
Views
9 Years
Discussion Span
Last Post by Gribouillis

An additional method to plot the function using pylab (matplotlib)

``````def plot_pylab (self ):
"""func.plot_pylab() -> plots func using pylab.
The function returns when the users closes the pylab window."""
import pylab
pylab .plot (self ._x ,self ._y ,"b")
pylab .show ()``````
Have something to contribute to this discussion? Please be thoughtful, detailed and courteous, and be sure to adhere to our posting rules.