def _infer_interval_breaks(coord, kind=None):
"""
Interpolate the bounds from the data in coord
Parameters
----------
%(CFDecoder.get_plotbounds.parameters.no_ignore_shape)s
Returns
-------
%(CFDecoder.get_plotbounds.returns)s
Notes
-----
this currently only works for rectilinear grids"""
if coord.ndim == 1:
return _infer_interval_breaks(coord)
elif coord.ndim == 2:
from scipy.interpolate import interp2d
kind = kind or rcParams['decoder.interp_kind']
y, x = map(np.arange, coord.shape)
new_x, new_y = map(_infer_interval_breaks, [x, y])
coord = np.asarray(coord)
return interp2d(x, y, coord, kind=kind, copy=False)(new_x, new_y)
python类interp2d()的实例源码
def _interp2_r0(Data, Pow2factor, kind, disp=False):
if disp:
p30 = np.poly1d(np.polyfit([10, 30, 50, 60, 70, 80],
np.log([0.8, 1.3, 6, 12, 28, 56]), 2))
print('Expectd time for NxN data:', np.exp(p30(Data.shape[0])))
x = np.arange(Data.shape[1])
y = np.arange(Data.shape[0])
xv, yv = np.meshgrid(x, y)
f = interpolate.interp2d(xv, yv, Data, kind=kind)
xnew = np.arange(0, Data.shape[1], 1 / (2**Pow2factor))
ynew = np.arange(0, Data.shape[0], 1 / (2**Pow2factor))
Upsampled = f(xnew, ynew)
return Upsampled
def eps_func(self):
interp_real = interpolate.interp2d(self.x, self.y, self.eps.real)
interp_imag = interpolate.interp2d(self.x, self.y, self.eps.imag)
interp = lambda x, y: interp_real(x, y) + 1.j*interp_imag(x, y)
return interp
def n_func(self):
return interpolate.interp2d(self.x, self.y, self.n)
def __init__(self, today, proj, disc, factory, vols, mode):
self.proj = proj
self.disc = disc
self.factory = factory
self.today = today
self.spotdate = factory.roll_date(today, 'spot') # also for marking startdate of hw vol term structure
self.surface = interp2d(vols.columns, vols.index, vols.values, kind='linear')
self.mode = 'capfloor'
class Caplet(BlackLognormalModel if mode == 'lognormal' else BlackNormalModel):
def __init__(self, libor, forward, blackvol, opt_term, annuity):
self.underlier = libor
super(Caplet, self).__init__(forward, blackvol, opt_term, annuity)
self.Caplet = Caplet
def bilin_interp(grid2d, ith, ir, dth, dr, nth, nr):
'''
This function finds an interpolated value in the 2d plume grid, and returns
the value. It's much faster than scipy.interpolate.interp2d
Takes as arguments the 2d grid of plume perturbations, the index in theta,
the index in radius, step size in theta and radius, and number of points
in theta and radius.
'''
value = grid2d[ith,ir]*(1-dth)*(1-dr) + grid2d[ith+1,ir]*dth*(1-dr)
value += grid2d[ith,ir+1]*(1-dth)*dr + grid2d[ith+1,ir+1]*dth*dr
return value
def __init__(self, input_data, output_data, name=""):
# check type and dimensions
assert isinstance(input_data, list)
assert isinstance(output_data, np.ndarray)
# output_data has to contain at least len(input_data) dimensions
assert len(input_data) <= len(output_data.shape)
for dim in range(len(input_data)):
assert len(input_data[dim]) == output_data.shape[dim]
self.input_data = input_data
if output_data.size == 0:
raise ValueError("No initialisation possible with an empty array!")
self.output_data = output_data
self.min = output_data.min()
self.max = output_data.max()
self.name = name
# self._interpolator = si.interp2d(input_data[0], input_data[1], output_data, bounds_error=True)
def __init__(self, eval_data_list, time_point=None, spatial_point=None, ylabel="",
legend_label=None, legend_location=1, figure_size=(10, 6)):
if not ((isinstance(time_point, Number) ^ isinstance(spatial_point, Number)) and \
(isinstance(time_point, type(None)) ^ isinstance(spatial_point, type(None)))):
raise TypeError("Only one kwarg *_point can be passed,"
"which has to be an instance from type numbers.Number")
DataPlot.__init__(self, eval_data_list)
plt.figure(facecolor='white', figsize=figure_size)
plt.ylabel(ylabel)
plt.grid(True)
# TODO: move to ut.EvalData
len_data = len(self._data)
interp_funcs = [si.interp2d(eval_data.input_data[1], eval_data.input_data[0], eval_data.output_data)
for eval_data in eval_data_list]
if time_point is None:
slice_input = [data_set.input_data[0] for data_set in self._data]
slice_data = [interp_funcs[i](spatial_point, slice_input[i]) for i in range(len_data)]
plt.xlabel('$t$')
elif spatial_point is None:
slice_input = [data_set.input_data[1] for data_set in self._data]
slice_data = [interp_funcs[i](slice_input[i], time_point) for i in range(len_data)]
plt.xlabel('$z$')
else:
raise TypeError
if legend_label is None:
show_leg = False
legend_label = [evald.name for evald in eval_data_list]
else:
show_leg = True
for i in range(0, len_data):
plt.plot(slice_input[i], slice_data[i], label=legend_label[i])
if show_leg:
plt.legend(loc=legend_location)
def _interpolate(x1, y1, z1, x2, y2):
"""Helper to interpolate in 1d or 2d
We interpolate to get the same shape than with other methods.
"""
if x1.size > 1 and y1.size > 1:
func = interp2d(x1, y1, z1.T, kind='linear', bounds_error=False)
z2 = func(x2, y2)
elif x1.size == 1 and y1.size > 1:
func = interp1d(y1, z1.ravel(), kind='linear', bounds_error=False)
z2 = func(y2)
elif y1.size == 1 and x1.size > 1:
func = interp1d(x1, z1.ravel(), kind='linear', bounds_error=False)
z2 = func(x2)
else:
raise ValueError("Can't interpolate a scalar.")
# interp2d is not intuitive and return this shape:
z2.shape = (y2.size, x2.size)
return z2
def barycentric_trapezoidial_interpolation(Nx,Ny,p,hexagonalOffset=0.5):
'''
define our function to calculate the position of points from Nx columns and Ny rows.
The vertices are defined by p which is a size(4,2) array.
each row of p are the coordinates or the vertices of our trapezoid
the vertices have to be given in a specific order:
[[1 1]
[1 2]
[2 1]
[2 2]]
an example plot using the barycentric interpolation to regrid data
define number of rows and number of columns and the vertices, then make some plots
Example:
Nx = 20
Ny = 15
coords = [[0,0],[0,1],[1,0],[1,1]] #these are the [x,y] coords of your 4 draggable corners
coords = np.asarray(coords)
f, ax = plt.subplots(2, 2) # sharey=True, sharex=True)
for i,a in enumerate(ax.flatten()):
newCoords = coords[:]
if i > 0:
newCoords = newCoords + np.random.rand(4,2) / 5
xi,yi = barycentric_trapezoidial_interpolation(Nx,Ny,newCoords)
a.plot(xi,yi,'.',markersize=12)
plt.show()
'''
x_basis = np.linspace(0,1,Nx)
y_basis = np.linspace(0,1,Ny)
px = [[p[0,0], p[2,0]],[p[1,0], p[3,0]]] #these are the [2,2] x-coordinates
py = [[p[0,1], p[2,1]],[p[1,1], p[3,1]]] #these are the [2,2] y-coordinates
fx = interpolate.interp2d([0,1], [0,1], px, kind='linear')
xi = fx(x_basis[:],y_basis[:]).flatten()
fy = interpolate.interp2d([0,1], [0,1], py, kind='linear')
yi = fy(x_basis[:],y_basis[:]).flatten()
d1 = (p[2,0] - p[0,0]) / Nx / 2.0
d2 = (p[3,0] - p[1,0]) / Nx / 2.0
offset = (d1 + d2) * hexagonalOffset
#every other row will be shifted in diff(x) * hexagonalOffset
for i in range(0,len(xi)-Nx,Nx*2):
for j in range(Nx):
xi[i+j+Nx] += offset
return xi,yi
def interpolate_curvature(curvatures, times, kind='cubic'):
""" Interpolate the curvatures.
A curvature is a parametrisation `s`, d(angle)/ ds and a parameter between [0,1].
Return a surface f(s,t).
"""
import numpy as np
from scipy import interpolate
from scipy.ndimage import measurements
if len(curvatures) == 3 and not isinstance(curvatures[2], tuple):
curvatures= [curvatures]
curv_abs = [c[0] for c in curvatures]
curves = [c[1] for c in curvatures]
params = [c[2] for c in curvatures]
n = len(params)
if n==1:
curv_abs*=2
curves*=2
curves[0]=np.zeros(len(curves[0]))
params = [0.,1.]
elif False:
if params[0] != 0.:
params.insert(0,0.)
curves.insert(0,curves[0])
curv_abs.insert(0, curv_abs[0])
if params[-1] != 0:
params.append(1.)
curves.append(curves[-1])
curv_abs.append(curv_abs[-1])
# compute a common parametrisation s
# We conserve the last parametrisation because the result is very sensitive
if True:
min_s = min(np.diff(s).min() for s in curv_abs)
s_new = np.unique(np.array(curv_abs).flatten())
ds = np.diff(s_new)
k = np.cumsum(ds >= min_s)
labels = range(k.min(), k.max()+1)
s = np.zeros(len(labels)+1)
s[1:] = measurements.mean(s_new[1:], k, labels)
s[-1]=1.
s = s_new
else:
s = np.array(curv_abs[-1])
# renormalise all the curves
curves = [np.interp(s[1:-1], old_s[1:-1], old_crv) for old_s, old_crv in zip(curv_abs,curves)]
# interpolate correctly the curvatures
x = s[1:-1]
y = np.array(params)
z = np.array(curves)
#f = interpolate.interp2d(y, x, z, kind=kind)
f = interpolate.RectBivariateSpline(x, y, z.T, kx=1, ky=1)
return s, f(x,times)
