#!/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 )
def __add__ (self ,other ):
"""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()
"""
