python.interpolate Namespace Reference

Functions
def	_bisect (a, x, lo=0, hi=None)
def	_bisect2 (a, x, xcol, lo=0, hi=None)
def	interpolate (points, x, degree, ycol=1, regions=[], xcol=0)

Variables
dictionary	__test__ = {}

Function Documentation

_bisect()

def python.interpolate._bisect (   a,   x,   lo = 0,   hi = None  )

private

Definition at line 21 of file interpolate.py.

def _bisect (a, x, lo=0, hi=None):
    if hi is None:
        hi = len(a)
    while lo < hi:
        mid = (lo+hi)//2
        if x <= a[mid]:
            hi = mid
        else:
            lo = mid+1
    return lo
 
# Helper: do a binary search in a 2D table.

_bisect2()

def python.interpolate._bisect2 (   a,   x,   xcol,   lo = 0,   hi = None  )

private

Definition at line 33 of file interpolate.py.

def _bisect2 (a, x, xcol, lo=0, hi=None):
    if hi is None:
        hi = len(a)
    while lo < hi:
        mid = (lo+hi)//2
        if x <= a[mid][xcol]:
            hi = mid
        else:
            lo = mid+1
    return lo
 
 
#
#  Polynomial interpolation in a table.
#

interpolate()

def python.interpolate.interpolate	(		points,
			x,
			degree,
			ycol = 1,
			regions = [],
			xcol = 0
	)

Polynomial interpolation in a table.

  points   Interpolation table.
       The x values are in the first column (column 0),
       and by default, the y values are in the second column (column 1).
       The x values must be in ascending order, with no duplicates.
  x        The value to interpolate.
  degree   The degree of the interpolating polynomial.
  ycol     The column number of the y values.  (0-based.)
  regions  Sometimes, you want to divide the range being interpolated
       into several distinct regions, where the interpolated
       function may be discontinuous at the boundaries
       between regions.  To do this, supply this argument,
       which should contain a list (in ascending order)
       of the x-coordinates of the region boundaries.
       When the interpolation runs against a boundary,
       the algorithm will add a copy of the last point
       before the boundary, positioned at the boundary.
       This helps to control runaway extrapolation
       leading up to the boundary.
  Returns the interpolated value.
 
  The method used is Newtonian interpolation.
  Based on the cernlib routine divdif.

Definition at line 48 of file interpolate.py.

def interpolate (points, x, degree, ycol=1, regions=[], xcol=0):
    """Polynomial interpolation in a table.
    
  points   Interpolation table.
           The x values are in the first column (column 0),
           and by default, the y values are in the second column (column 1).
           The x values must be in ascending order, with no duplicates.
  x        The value to interpolate.
  degree   The degree of the interpolating polynomial.
  ycol     The column number of the y values.  (0-based.)
  regions  Sometimes, you want to divide the range being interpolated
           into several distinct regions, where the interpolated
           function may be discontinuous at the boundaries
           between regions.  To do this, supply this argument,
           which should contain a list (in ascending order)
           of the x-coordinates of the region boundaries.
           When the interpolation runs against a boundary,
           the algorithm will add a copy of the last point
           before the boundary, positioned at the boundary.
           This helps to control runaway extrapolation
           leading up to the boundary.
  Returns the interpolated value.
 
  The method used is Newtonian interpolation.
  Based on the cernlib routine divdif.
  """
    if len (points) < 2 or degree < 1:
        raise Exception ('bad args!')
    degree = min (degree, len (points) - 1)
    
 
    # START.  FIND SUBSCRIPT IX OF X IN ARRAY.
    #ix = _bisect (points, x, lambda x, p, xcol=xcol: cmp (x, p[xcol]))# - 1
    ix = _bisect2 (points, x, xcol)
 
    ir = _bisect (regions, x)
 
    # Number of points to try for.
    # Either degree+1 or degree+2, whichever is even,
    # to give the same number of points on each side.
    # If we run up against an edge or a boundary, we'll
    # fall back to using just degree+1 points (or fewer if we can't
    # even get that many).
    # If we end up using degree+2 points, we'll do two interpolations
    # of degree degree and average them.
    npts = degree + 2 - (degree%2)
    l = 0  # noqa: E741
    t = []
    d = []
 
    # If we run up against the edge of a region boundary,
    # we'll want to add a psuedopoint right at the boundary
    # (copying the point closest to the boundary) instead of the
    # point farthest away from it.
    extralo = 0
    extrahi = 0
    
    # Starting point index, not considering edges or boundaries.
    ilo = ix - npts // 2
 
    # Make sure this point is within the array range and has not
    # crossed a region boundary.
    if ilo < 0:
        ilo = 0
        npts = degree+1
    while ilo < len (points) and ir>0 and points[ilo][xcol] < regions[ir-1]:
        ilo += 1
        npts = degree+1
        extralo = 1
 
    # Same deal for the right hand edge.
    # ihi is one past the last point to use.
    himoved = 0
    ihi = ilo + npts
    if ihi > len (points):
        ihi = len (points)
        npts = degree+1
        himoved = 1
    while ihi > 0 and ir<len(regions) and points[ihi-1][xcol] >= regions[ir]:
        ihi -= 1
        npts = degree+1
        extrahi = 1
        himoved = 1
 
    lomoved = 0
    ilo = ihi - npts
    if ilo < 0:
        ilo = 0
        lomoved = 1
    while ilo < len (points) and ir>0 and points[ilo][xcol] < regions[ir-1]:
        ilo += 1
        extralo = 1
        lomoved = 1
 
    npts = ihi - ilo
    t = []
    d = []
 
    if extralo and points[ilo][xcol] != regions[ir-1]:
        if not himoved:
            ihi -= 1
        else:
            npts += 1
        t.append (regions[ir-1])
        d.append (points[ilo][ycol])
 
    if extrahi and points[ihi-1][xcol] != regions[ir]:
        if not lomoved:
            ilo += 1
        else:
            npts += 1
        t.append (regions[ir])
        d.append (points[ihi-1][ycol])
 
    t += [points[i][xcol] for i in range (ilo, ihi)]
    d += [points[i][ycol] for i in range (ilo, ihi)]
 
    degree = min (degree, npts-1)
    extra = npts != degree+1
 
    if extra:
        (t[0], t[npts-2]) = (t[npts-2], t[0])
        (d[0], d[npts-2]) = (d[npts-2], d[0])
 
    # REPLACE D BY THE LEADING DIAGONAL OF A DIVIDED-DIFFERENCE TABLE,
    # SUPPLEMENTED BY AN EXTRA LINE IF *EXTRA* IS TRUE.
    for l in range(0, degree):
        if extra:
            d[degree+1] = (d[degree+1]-d[degree-1]) / (t[degree+1]-t[degree-1-l])
        for i in range (degree, l, -1):
            d[i] = (d[i]-d[i-1]) / (t[i]-t[i-1-l])
 
    # EVALUATE THE NEWTON INTERPOLATION FORMULA AT X, AVERAGING TWO VALUES
    # OF LAST DIFFERENCE IF *EXTRA* IS TRUE.
    sum = d[degree]
    if extra:
        sum=0.5*(sum+d[degree+1])
    for j in range (degree-1, -1, -1):
        sum = d[j] + (x - t[j]) * sum
 
    return sum
 
 

Variable Documentation

__test__

dictionary python.interpolate.__test__ = {}

private

Definition at line 191 of file interpolate.py.