Piecewise linear continuous functions.

#!/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()
"""