Piecewise linear continuous functions.

Gribouillis 1 Tallied Votes 583 Views Share

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 )
  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()
"""
Gribouillis 1,391 Programming Explorer Team Colleague

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 ()
Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.