ATLAS Offline Software
Functions | Variables
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 22 of file interpolate.py.

22 def _bisect (a, x, lo=0, hi=None):
23  if hi is None:
24  hi = len(a)
25  while lo < hi:
26  mid = (lo+hi)//2
27  if x <= a[mid]:
28  hi = mid
29  else:
30  lo = mid+1
31  return lo
32 
33 # 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 34 of file interpolate.py.

34 def _bisect2 (a, x, xcol, lo=0, hi=None):
35  if hi is None:
36  hi = len(a)
37  while lo < hi:
38  mid = (lo+hi)//2
39  if x <= a[mid][xcol]:
40  hi = mid
41  else:
42  lo = mid+1
43  return lo
44 
45 
46 #
47 # Polynomial interpolation in a table.
48 #

◆ 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 49 of file interpolate.py.

49 def interpolate (points, x, degree, ycol=1, regions=[], xcol=0):
50  """Polynomial interpolation in a table.
51 
52  points Interpolation table.
53  The x values are in the first column (column 0),
54  and by default, the y values are in the second column (column 1).
55  The x values must be in ascending order, with no duplicates.
56  x The value to interpolate.
57  degree The degree of the interpolating polynomial.
58  ycol The column number of the y values. (0-based.)
59  regions Sometimes, you want to divide the range being interpolated
60  into several distinct regions, where the interpolated
61  function may be discontinuous at the boundaries
62  between regions. To do this, supply this argument,
63  which should contain a list (in ascending order)
64  of the x-coordinates of the region boundaries.
65  When the interpolation runs against a boundary,
66  the algorithm will add a copy of the last point
67  before the boundary, positioned at the boundary.
68  This helps to control runaway extrapolation
69  leading up to the boundary.
70  Returns the interpolated value.
71 
72  The method used is Newtonian interpolation.
73  Based on the cernlib routine divdif.
74  """
75  if len (points) < 2 or degree < 1:
76  raise Exception ('bad args!')
77  degree = min (degree, len (points) - 1)
78 
79 
80  # START. FIND SUBSCRIPT IX OF X IN ARRAY.
81  #ix = _bisect (points, x, lambda x, p, xcol=xcol: cmp (x, p[xcol]))# - 1
82  ix = _bisect2 (points, x, xcol)
83 
84  ir = _bisect (regions, x)
85 
86  # Number of points to try for.
87  # Either degree+1 or degree+2, whichever is even,
88  # to give the same number of points on each side.
89  # If we run up against an edge or a boundary, we'll
90  # fall back to using just degree+1 points (or fewer if we can't
91  # even get that many).
92  # If we end up using degree+2 points, we'll do two interpolations
93  # of degree degree and average them.
94  npts = degree + 2 - (degree%2)
95  l = 0 # noqa: E741
96  t = []
97  d = []
98 
99  # If we run up against the edge of a region boundary,
100  # we'll want to add a psuedopoint right at the boundary
101  # (copying the point closest to the boundary) instead of the
102  # point farthest away from it.
103  extralo = 0
104  extrahi = 0
105 
106  # Starting point index, not considering edges or boundaries.
107  ilo = ix - npts // 2
108 
109  # Make sure this point is within the array range and has not
110  # crossed a region boundary.
111  if ilo < 0:
112  ilo = 0
113  npts = degree+1
114  while ilo < len (points) and ir>0 and points[ilo][xcol] < regions[ir-1]:
115  ilo += 1
116  npts = degree+1
117  extralo = 1
118 
119  # Same deal for the right hand edge.
120  # ihi is one past the last point to use.
121  himoved = 0
122  ihi = ilo + npts
123  if ihi > len (points):
124  ihi = len (points)
125  npts = degree+1
126  himoved = 1
127  while ihi > 0 and ir<len(regions) and points[ihi-1][xcol] >= regions[ir]:
128  ihi -= 1
129  npts = degree+1
130  extrahi = 1
131  himoved = 1
132 
133  lomoved = 0
134  ilo = ihi - npts
135  if ilo < 0:
136  ilo = 0
137  lomoved = 1
138  while ilo < len (points) and ir>0 and points[ilo][xcol] < regions[ir-1]:
139  ilo += 1
140  extralo = 1
141  lomoved = 1
142 
143  npts = ihi - ilo
144  t = []
145  d = []
146 
147  if extralo and points[ilo][xcol] != regions[ir-1]:
148  if not himoved:
149  ihi -= 1
150  else:
151  npts += 1
152  t.append (regions[ir-1])
153  d.append (points[ilo][ycol])
154 
155  if extrahi and points[ihi-1][xcol] != regions[ir]:
156  if not lomoved:
157  ilo += 1
158  else:
159  npts += 1
160  t.append (regions[ir])
161  d.append (points[ihi-1][ycol])
162 
163  t += [points[i][xcol] for i in range (ilo, ihi)]
164  d += [points[i][ycol] for i in range (ilo, ihi)]
165 
166  degree = min (degree, npts-1)
167  extra = npts != degree+1
168 
169  if extra:
170  (t[0], t[npts-2]) = (t[npts-2], t[0])
171  (d[0], d[npts-2]) = (d[npts-2], d[0])
172 
173  # REPLACE D BY THE LEADING DIAGONAL OF A DIVIDED-DIFFERENCE TABLE,
174  # SUPPLEMENTED BY AN EXTRA LINE IF *EXTRA* IS TRUE.
175  for l in range(0, degree):
176  if extra:
177  d[degree+1] = (d[degree+1]-d[degree-1]) / (t[degree+1]-t[degree-1-l])
178  for i in range (degree, l, -1):
179  d[i] = (d[i]-d[i-1]) / (t[i]-t[i-1-l])
180 
181  # EVALUATE THE NEWTON INTERPOLATION FORMULA AT X, AVERAGING TWO VALUES
182  # OF LAST DIFFERENCE IF *EXTRA* IS TRUE.
183  sum = d[degree]
184  if extra:
185  sum=0.5*(sum+d[degree+1])
186  for j in range (degree-1, -1, -1):
187  sum = d[j] + (x - t[j]) * sum
188 
189  return sum
190 
191 

Variable Documentation

◆ __test__

dictionary python.interpolate.__test__ = {}
private

Definition at line 192 of file interpolate.py.

python.interpolate._bisect
def _bisect(a, x, lo=0, hi=None)
Definition: interpolate.py:22
python.interpolate._bisect2
def _bisect2(a, x, xcol, lo=0, hi=None)
Definition: interpolate.py:34
python.interpolate.interpolate
def interpolate(points, x, degree, ycol=1, regions=[], xcol=0)
Definition: interpolate.py:49
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195