Source code for solcore.interpolate

# from scipy import interpolate as I
import numpy as np
from scipy.interpolate import make_interp_spline, splev


[docs]class interp1d(object): """ Interpolate a 1D function. See Also -------- splrep, splev - spline interpolation based on FITPACK UnivariateSpline - a more recent wrapper of the FITPACK routines """ def __init__(self, x, y, kind='linear', axis=-1, copy=True, bounds_error=False, fill_value=np.nan): """ Initialize a 1D linear interpolation class. Description ----------- x and y are arrays of values used to approximate some function f: y = f(x) This class returns a function whose call method uses linear interpolation to find the value of new points. Parameters ---------- x : array A 1D array of monotonically increasing real values. x cannot include duplicate values (otherwise f is overspecified) y : array An N-D array of real values. y's length along the interpolation axis must be equal to the length of x. kind : str or int Specifies the kind of interpolation as a string ('linear', 'nearest', 'zero', 'slinear', 'quadratic, 'cubic') or as an integer specifying the order of the spline interpolator to use. axis : int Specifies the axis of y along which to interpolate. Interpolation defaults to the last axis of y. copy : bool If True, the class makes internal copies of x and y. If False, references to x and y are used. The default is to copy. bounds_error : bool If True, an error is thrown any time interpolation is attempted on a value outside of the range of x (where extrapolation is necessary). If False, out of bounds values are assigned fill_value. By default, an error is raised. fill_value : float If provided, then this value will be used to fill in for requested points outside of the data range. If not provided, then the default is NaN. """ self.copy = copy self.bounds_error = bounds_error self.fill_value = fill_value if kind in ['zero', 'slinear', 'quadratic', 'cubic']: order = {'nearest': 0, 'zero': 0, 'slinear': 1, 'quadratic': 2, 'cubic': 3}[kind] kind = 'spline' elif isinstance(kind, int): order = kind kind = 'spline' elif kind not in ('linear', 'nearest'): raise NotImplementedError("%s is unsupported: Use fitpack " \ "routines for other types." % kind) x = np.array(x, copy=self.copy) y = np.array(y, copy=self.copy) if x.ndim != 1: raise ValueError("the x array must have exactly one dimension.") if y.ndim == 0: raise ValueError("the y array must have at least one dimension.") # Force-cast y to a floating-point type, if it's not yet one if not issubclass(y.dtype.type, np.inexact): y = y.astype(np.float_) # Normalize the axis to ensure that it is positive. self.axis = axis % len(y.shape) self._kind = kind if kind in ('linear', 'nearest'): # Make a "view" of the y array that is rotated to the interpolation # axis. axes = list(range(y.ndim)) del axes[self.axis] axes.append(self.axis) oriented_y = y.transpose(axes) minval = 2 len_y = oriented_y.shape[-1] if kind == 'linear': self._call = self._call_linear elif kind == 'nearest': self.x_bds = (x[1:] + x[:-1]) / 2.0 self._call = self._call_nearest else: axes = list(range(y.ndim)) del axes[self.axis] axes.insert(0, self.axis) oriented_y = y.transpose(axes) minval = order + 1 len_y = oriented_y.shape[0] self._call = self._call_spline self._spline = make_interp_spline(x, oriented_y, k=order) len_x = len(x) if len_x != len_y: raise ValueError("x and y arrays must be equal in length along " "interpolation axis.") if len_x < minval: raise ValueError("x and y arrays must have at " "least %d entries" % minval) self.x = x self.y = oriented_y def _call_linear(self, x_new): # 2. Find where in the orignal data, the values to interpolate # would be inserted. # Note: If x_new[n] == x[m], then m is returned by searchsorted. x_new_indices = np.searchsorted(self.x, x_new) # 3. Clip x_new_indices so that they are within the range of # self.x indices and at least 1. Removes mis-interpolation # of x_new[n] = x[0] x_new_indices = x_new_indices.clip(1, len(self.x) - 1).astype(int) # 4. Calculate the slope of regions that each x_new value falls in. lo = x_new_indices - 1 hi = x_new_indices x_lo = self.x[lo] x_hi = self.x[hi] y_lo = self.y[..., lo] y_hi = self.y[..., hi] # Note that the following two expressions rely on the specifics of the # broadcasting semantics. slope = (y_hi - y_lo) / (x_hi - x_lo) # 5. Calculate the actual value for each entry in x_new. y_new = slope * (x_new - x_lo) + y_lo return y_new def _call_nearest(self, x_new): """ Find nearest neighbour interpolated y_new = f(x_new).""" # 2. Find where in the averaged data the values to interpolate # would be inserted. # Note: use side='left' (right) to searchsorted() to define the # halfway point to be nearest to the left (right) neighbour x_new_indices = np.searchsorted(self.x_bds, x_new, side='left') # 3. Clip x_new_indices so that they are within the range of x indices. x_new_indices = x_new_indices.clip(0, len(self.x) - 1).astype(np.intp) # 4. Calculate the actual value for each entry in x_new. y_new = self.y[..., x_new_indices] return y_new def _call_spline(self, x_new): x_new = np.asarray(x_new) result = splev(x_new.ravel(), self._spline) return result.reshape(x_new.shape + result.shape[1:]) def __call__(self, x_new): """Find interpolated y_new = f(x_new). Parameters ---------- x_new : number or array New independent variable(s). Returns ------- y_new : ndarray Interpolated value(s) corresponding to x_new. """ # 1. Handle values in x_new that are outside of x. Throw error, # or return a list of mask array indicating the outofbounds values. # The behavior is set by the bounds_error variable. x_new = np.asarray(x_new) out_of_bounds = self._check_bounds(x_new) y_new = self._call(x_new) # Rotate the values of y_new back so that they correspond to the # correct x_new values. For N-D x_new, take the last (for linear) # or first (for other splines) N axes # from y_new and insert them where self.axis was in the list of axes. nx = x_new.ndim ny = y_new.ndim # 6. Fill any values that were out of bounds with fill_value. # and # 7. Rotate the values back to their proper place. if nx == 0: # special case: x is a scalar if out_of_bounds: if ny == 0: return np.asarray(self.fill_value) else: y_new[...] = self.fill_value return np.asarray(y_new) elif self._kind in ('linear', 'nearest'): y_new[..., out_of_bounds] = self.fill_value axes = list(range(ny - nx)) axes[self.axis:self.axis] = list(range(ny - nx, ny)) return y_new.transpose(axes) else: y_new[out_of_bounds] = self.fill_value axes = list(range(nx, ny)) axes[self.axis:self.axis] = list(range(nx)) return y_new.transpose(axes) def _check_bounds(self, x_new): """Check the inputs for being in the bounds of the interpolated data. Parameters ---------- x_new : array Returns ------- out_of_bounds : bool array The mask on x_new of values that are out of the bounds. """ # If self.bounds_error is True, we raise an error if any x_new values # fall outside the range of x. Otherwise, we return an array indicating # which values are outside the boundary region. below_bounds = x_new < self.x[0] above_bounds = x_new > self.x[-1] # !! Could provide more information about which values are out of bounds if self.bounds_error and below_bounds.any(): raise ValueError("A value in x_new is below the interpolation " "range.") if self.bounds_error and above_bounds.any(): raise ValueError("A value in x_new is above the interpolation " "range.") # !! Should we emit a warning if some values are out of bounds? # !! matlab does not. out_of_bounds = np.logical_or(below_bounds, above_bounds) return out_of_bounds
[docs]class BilinearInterpolation(object): """docstring for BilinearInterpolation""" def __init__(self, x=None, y=None, z=None, fill_value=0.): super(BilinearInterpolation, self).__init__() self.x = np.array(x) self.y = np.array(y) self.z = np.array(z) # Check the shapes of x,y and z x_s = self.x.shape y_s = self.y.shape z_s = self.z.shape assert x_s[0] >= 3, "Need at least 3 points in x." assert y_s[0] >= 3, "Need at least 3 points in y." assert x_s[0] == z_s[ 0], "The number of rows in z (which is %s), must be the same as number of elements in x (which is %s)" % ( x_s[0], z_s[0]) assert y_s[0] == z_s[ 1], "The number of columns in z (which is %s), must be the same as number of elements in y (which is %s)" % ( y_s[0], z_s[1]) def __call__(self, xvalue, yvalue): """The intepolated value of the surface.""" try: x_i = 0 found = False for i in list(range(0, len(self.x) - 1)): # print self.x[i], "<=", xvalue, "<", self.x[i+1] if self.x[i] <= xvalue < self.x[i + 1]: # print "Found, returning index", i x_i = i found = True break if not found: # print "Not found, returning 0." return 0. x_v = (self.x[x_i], self.x[x_i + 1]) x_i = (x_i, x_i + 1) y_i = 0 found = False for i in list(range(0, len(self.y) - 1)): # print self.y[i], "<=", yvalue, "<", self.y[i+1] if self.y[i] <= yvalue < self.y[i + 1]: # print "Found, returning index", i y_i = i found = True break if not found: # print "Not found, returning 0." return 0. y_v = (self.y[y_i], self.y[y_i + 1]) y_i = (y_i, y_i + 1) D = (x_v[1] - x_v[0]) * (y_v[1] - y_v[0]) ta = self.z[x_i[0], y_i[0]] tb = (x_v[1] - xvalue) tc = (y_v[1] - yvalue) t1 = ta * tb * tc / D # print t1, "=", ta, "*", tb, "*", tc, "/", D ta = self.z[x_i[1], y_i[0]] tb = (xvalue - x_v[0]) tc = (y_v[1] - yvalue) t2 = ta * tb * tc / D # print t1, "=", ta, "*", tb, "*", tc, "/", D ta = self.z[x_i[0], y_i[1]] tb = (x_v[1] - xvalue) tc = (yvalue - y_v[0]) t3 = ta * tb * tc / D # print t1, "=", ta, "*", tb, "*", tc, "/", D ta = self.z[x_i[1], y_i[1]] tb = (xvalue - x_v[0]) tc = (yvalue - y_v[0]) t4 = ta * tb * tc / D # print t1, "=", ta, "*", tb, "*", tc, "/", D # print "Therefore," # print t1, "+", t2, "+", t3, "+", t4 # print self.z return t1 + t2 + t3 + t4 except ValueError: # import pdb; pdb.set_trace() # Is one of the inputs an array? if type(xvalue) == list or type(xvalue) == tuple or type(xvalue) == type(np.array([])): constructed_list = [] for sub_xvalue in xvalue: constructed_list.append(self(sub_xvalue, yvalue)) return np.array(constructed_list)
if __name__ == "__main__": # Some wavelength points x = np.array([400., 600., 601, 1000.]) # Some angle values y = np.array([0., np.pi / 4, np.pi / 2]) # Some reflectivity values # rads @ 400nm, rads @ 600nm, rads @ 601m, rads @ 1000nm z = np.array([[0., 0., 0.], [1e-9, 1e-9, 1e-9], [1 - 1e-9, 1 - 1e-9, 1 - 1e-9], [1., 1., 1.]]) # print z.shape z_i = BilinearInterpolation(x, y, z) print(z_i(610, 0.1))