当前位置: 首页 > news >正文

CQF预备知识:Python相关库 -- 插值过渡指南 scipy.interpolate

文中内容仅限技术学习与代码实践参考,市场存在不确定性,技术分析需谨慎验证,不构成任何投资建议。

插值过渡指南

1. 如何从使用 interp2d 过渡

interp2d 会在二维规则网格插值和二维散乱数据插值之间静默切换。切换基于 (raveled) xyz 数组的长度。简而言之,对于规则网格,使用 scipy.interpolate.RectBivariateSpline;对于散乱插值,使用 bisprep/bisplev 组合。下面我们将给出逐点过渡的示例代码,这些代码应该能够精确地保留 interp2d 的结果。

1.1 interp2d 在规则网格上的应用

我们从(稍作修改的)文档示例开始。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import interp2d, RectBivariateSpline>>> x = np.arange(-5.01, 5.01, 0.25)
>>> y = np.arange(-5.01, 7.51, 0.25)
>>> xx, yy = np.meshgrid(x, y)
>>> z = np.sin(xx**2 + 2.*yy**2)
>>> f = interp2d(x, y, z, kind='cubic')

这是“规则网格”代码路径,因为

>>> z.size == len(x) * len(y)
True

另外,注意 x.size != y.size

>>> x.size, y.size
(41, 51)

现在,我们构建一个方便的函数来构造插值器并绘制它。

>>> def plot(f, xnew, ynew):
...    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
...    znew = f(xnew, ynew)
...    ax1.plot(x, z[0, :], 'ro-', xnew, znew[0, :], 'b-')
...    im = ax2.imshow(znew)
...    plt.colorbar(im, ax=ax2)
...    plt.show()
...    return znew
...
>>> xnew = np.arange(-5.01, 5.01, 1e-2)
>>> ynew = np.arange(-5.01, 7.51, 1e-2)
>>> znew_i = plot(f, xnew, ynew)

img

替代方法:使用 RectBivariateSpline,结果完全一致

注意转置:首先,在构造函数中,其次,你需要转置评估的结果。这是为了撤销 interp2d 所做的转置。

>>> r = RectBivariateSpline(x, y, z.T)
>>> rt = lambda xnew, ynew: r(xnew, ynew).T
>>> znew_r = plot(rt, xnew, ynew)

img

>>> from numpy.testing import assert_allclose
>>> assert_allclose(znew_i, znew_r, atol=1e-14)
插值阶数:线性、三次等

interp2d 默认为 kind="linear",即在 x-y- 两个方向上均为线性。
RectBivariateSpline 则默认为三次插值,kx=3, ky=3

以下是精确的等价关系:

interp2dRectBivariateSpline
无关键字参数kx = 1, ky = 1
kind='linear'kx = 1, ky = 1
kind='cubic'kx = 3, ky = 3

1.2 interp2d 与点的完整坐标(散乱插值)

这里,我们将上一练习中的网格数组展平,以说明其功能。

>>> xxr = xx.ravel()
>>> yyr = yy.ravel()
>>> zzr = z.ravel()
>>> f = interp2d(xxr, yyr, zzr, kind='cubic')

注意,这是“非规则网格”代码路径,用于散乱数据,满足 len(x) == len(y) == len(z)

>>> len(xxr) == len(yyr) == len(zzr)
True
>>> xnew = np.arange(-5.01, 5.01, 1e-2)
>>> ynew = np.arange(-5.01, 7.51, 1e-2)
>>> znew_i = plot(f, xnew, ynew)

img

替代方法:直接使用 scipy.interpolate.bisplrep / scipy.interpolate.bisplev
>>> from scipy.interpolate import bisplrep, bisplev
>>> tck = bisplrep(xxr, yyr, zzr, kx=3, ky=3, s=0)
# 为了方便起见:从 bisplev 构造一个可调用的函数
>>> ff = lambda xnew, ynew: bisplev(xnew, ynew, tck).T   # 注意转置,以模仿 interp2d 的行为
>>> znew_b = plot(ff, xnew, ynew)

img

>>> assert_allclose(znew_i, znew_b, atol=1e-15)
插值阶数:线性、三次等

interp2d 默认为 kind="linear",即在 x-y- 两个方向上均为线性。
bisplrep 则默认为三次插值,kx=3, ky=3

以下是精确的等价关系:

interp2dbisplrep
无关键字参数kx = 1, ky = 1
kind='linear'kx = 1, ky = 1
kind='cubic'kx = 3, ky = 3

2. interp2d 的替代方法:规则网格

对于新代码,推荐的替代方法是 RegularGridInterpolator。它是一个独立的实现,不基于 FITPACK。支持最近邻、线性插值以及奇数阶张量积样条。

样条节点保证与数据点重合。

注意,这里:

  1. 元组参数为 (x, y)
  2. z 数组需要转置
  3. 关键字名称为 method,而不是 kind
  4. bounds_error 参数默认为 True
>>> from scipy.interpolate import RegularGridInterpolator as RGI
>>> r = RGI((x, y), z.T, method='linear', bounds_error=False)

评估:创建一个二维网格数组。使用 indexing='ij'sparse=True 以节省一些内存:

>>> xxnew, yynew = np.meshgrid(xnew, ynew, indexing='ij', sparse=True)

评估,注意元组参数:

>>> znew_reggrid = r((xxnew, yynew))
>>> fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
# 再次注意转置,以撤销 `interp2d` 的约定
>>> znew_reggrid_t = znew_reggrid.T
>>> ax1.plot(x, z[0, :], 'ro-', xnew, znew_reggrid_t[0, :], 'b-')
>>> im = ax2.imshow(znew_reggrid_t)
>>> plt.colorbar(im, ax=ax2)

img

3. 散乱二维线性插值:优先选择 LinearNDInterpolator 而非 SmoothBivariateSplinebisplrep

对于二维散乱线性插值,SmoothBivariateSplinebiplrep 可能会发出警告,或者无法插值数据,或者产生远离数据点的样条节点。相反,优先选择 LinearNDInterpolator,它基于 QHull 的数据三角剖分。

# TestSmoothBivariateSpline::test_integral
>>> from scipy.interpolate import SmoothBivariateSpline, LinearNDInterpolator
>>> x = np.array([1,1,1,2,2,2,4,4,4])
>>> y = np.array([1,2,3,1,2,3,1,2,3])
>>> z = np.array([0,7,8,3,4,7,1,3,4])

现在,使用基于 Qhull 数据三角剖分的线性插值:

>>> xy = np.c_[x, y]   # 或者直接使用 list(zip(x, y))
>>> lut2 = LinearNDInterpolator(xy, z)
>>> X = np.linspace(min(x), max(x))
>>> Y = np.linspace(min(y), max(y))
>>> X, Y = np.meshgrid(X, Y)

结果易于理解和解释:

>>> fig = plt.figure()
>>> ax = fig.add_subplot(projection='3d')
>>> ax.plot_wireframe(X, Y, lut2(X, Y))
>>> ax.scatter(x, y, z,  'o', color='k', s=48)

img

注意,bisplrep 的行为是不同的!它可能会将样条节点放置在数据之外。

为了说明,考虑上一个示例中的相同数据:

>>> tck = bisplrep(x, y, z, kx=1, ky=1, s=0)
>>> fig = plt.figure()
>>> ax = fig.add_subplot(projection='3d')
>>> xx = np.linspace(min(x), max(x))
>>> yy = np.linspace(min(y), max(y))
>>> X, Y = np.meshgrid(xx, yy)
>>> Z = bisplev(xx, yy, tck)
>>> Z = Z.reshape(*X.shape).T
>>> ax.plot_wireframe(X, Y, Z, rstride=2, cstride=2)
>>> ax.scatter(x, y, z,  'o', color='k', s=48)

同样,SmoothBivariateSpline 无法插值数据。再次使用上一个示例中的相同数据。

>>> lut = SmoothBivariateSpline(x, y, z, kx=1, ky=1, s=0)
>>> fig = plt.figure()
>>> ax = fig.add_subplot(projection='3d')
>>> xx = np.linspace(min(x), max(x))
>>> yy = np.linspace(min(y), max(y))
>>> X, Y = np.meshgrid(xx, yy)
>>> ax.plot_wireframe(X, Y, lut(xx, yy).T, rstride=4, cstride=4)
>>> ax.scatter(x, y, z,  'o', color='k', s=48)

注意,SmoothBivariateSplinebisplrep 的结果都存在伪影,与 LinearNDInterpolator 的结果不同。这里说明的问题是针对线性插值的,但 FITPACK 的节点选择机制并不能保证在更高阶(例如三次)样条曲面上避免这些问题。

风险提示与免责声明
本文内容基于公开信息研究整理,不构成任何形式的投资建议。历史表现不应作为未来收益保证,市场存在不可预见的波动风险。投资者需结合自身财务状况及风险承受能力独立决策,并自行承担交易结果。作者及发布方不对任何依据本文操作导致的损失承担法律责任。市场有风险,投资须谨慎。

相关文章:

  • 接口测试常用工具及测试方法(基础篇)
  • [SKE]CPU 与 GPU 之间数据加密传输的认证与异常处理
  • 触觉智能RK3576核心板工业应用之软硬件全国产化,成功适配开源鸿蒙OpenHarmony5.0
  • aws s3 sdk c++使用指南、适配阿里云oss和aws
  • OCCT 中 BRepBuilderAPI_MakePolygon与BRepBuilderAPI_MakeWire
  • 5种常见的网络保密通信协议
  • 如何从 Ansys SpaceClaim 模型中提取 CAD 数据,该模型是在我计算机上安装的未来版本中创建的?
  • 亚马逊云服务器配置推荐
  • SMB协议在Windows内网中的核心地位
  • 华为:eSight网管平台使用snmp纳管交换机
  • React---Hooks深入
  • Shell 脚本:系统管理与任务自动化的利器
  • 基于Jenkins与Kubernetes的系统化变更管理实践
  • 【C/C++】创建文件夹
  • 单元测试基本步骤
  • vue3 全局过滤器
  • Bash (Bourne Again SHell)
  • Seedance:字节发布视频生成基础模型新SOTA,能力全面提升
  • 文本的数据分析
  • Redis+Kafka实现动态延时任务
  • 网站建设推广公司排名/关键词推广seo怎么优化
  • 网站轮换图片怎么做/百度招聘官网首页
  • 合山网站建设/如何建造自己的网站
  • 香港公司注册代理/seo外包公司专家
  • 中山做百度网站的公司/1688黄页大全进口
  • 门户网站建设制作/自己如何制作一个网站