Python中多级线性插值的高效方法
Efficient method for multi-level linear interpolation in Python
我目前正在进行一个估算流量计不确定性的项目。仪表不确定性基于四个不同的值:
- 液体流量 (liq)
- 流体粘度 (cP)
- 水液比 (wlr)
- 气体体积分数 (gvf)
第三方为仪表提供 tables 的 liq、cP、wlr 和 gvf 的多个不同值。正如您所猜测的那样,来自仪表的数据永远不会完全落入预定义值之一。例如一分钟的数据可能是:
- 液体流量:6532
- 流体粘度:22
- 水液比:0.412
- 气体体积分数:0.634
使用以上数据对 table 执行四向插值以找出不确定性。
我想出了一个解决方案,但它看起来很笨拙,我想知道是否有人有任何想法。我还是 pandas 游戏的新手,非常感谢看到其他人的解决方案。
最初我对数据进行排序以将 table 降低到高于和低于我正在寻找的实际点的值。
aliq = 6532 # stbpd
avisc = 22 # centipoise
awlr = 0.412 # water liquid ratio
agvf = 0.634 # gas volume fraction
def findclose(num, colm):
arr = colm.unique()
if num in arr:
clslo = num
clshi = num
else:
clslo = arr[arr > num].min() # close low value
clshi = arr[arr < num].max() # close high value
return [clslo, clshi]
df = tbl_vx52[
(tbl_vx52['liq'].isin(findclose(aliq,tbl_vx52['liq']))) &
(tbl_vx52['visc'].isin(findclose(avisc,tbl_vx52['visc']))) &
(tbl_vx52['wlr'].isin(findclose(awlr,tbl_vx52['wlr']))) &
(tbl_vx52['gvf'].isin(findclose(agvf,tbl_vx52['gvf'])))
].reset_index(drop=True)
table 值从 2240 减少到 16。而不是包括所有数据 (tbl_vx52)。我创建了一些要加载的代码,因此您可以看到名为 df 的子数据框的外观,其中仅包含此示例区域上方和下方的值。
df = pd.DataFrame({'liq':[5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 7000, 7000, 7000, 7000, 7000, 7000, 7000, 7000],
'visc':[10, 10, 10, 10, 30, 30, 30, 30, 10, 10, 10, 10, 30, 30, 30, 30],
'wlr':[0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5],
'gvf':[0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75],
'uncert':[0.0707, 0.0992, 0.0906, 0.1278, 0.0705, 0.0994, 0.091, 0.128, 0.0702, 0.0991, 0.0905, 0.1279, 0.0704, 0.0992, 0.0904, 0.1283],
})
完成了一些非常粗略的循环,以开始根据各个输入(liq、visc、wlr 或 gvf)对值进行配对。下面显示的是 gvf.
上的第一个循环
pairs = [
slice(0,1),
slice(2,3),
slice(4,5),
slice(6,7),
slice(8,9),
slice(10,11),
slice(12,13),
slice(14,15)]
for pair in pairs:
df.loc[pair,'uncert'] = np.interp(
agvf,
df.loc[pair,'gvf'],
df.loc[pair,'uncert']
)
df.loc[pair,'gvf'] = agvf
df = df.drop_duplicates().reset_index(drop=True)
删除了重复值,从 16 行减少到 8 行。然后再次为 wlr 重复此操作。
pairs = [
slice(0,1),
slice(2,3),
slice(4,5),
slice(6,7)
]
for pair in pairs:
df.loc[pair,'uncert'] = np.interp(
awlr,
df.loc[pair,'wlr'],
df.loc[pair,'uncert']
)
df.loc[pair,'wlr'] = awlr
df = df.drop_duplicates().reset_index(drop=True)
对 visc(四行)和最后的 liquid(两行)重复上述结构,直到子数组中只剩下一个值。这给出了您操作点的仪表不确定性。
我知道它很笨重。对不同方法的任何意见或想法表示赞赏。
好的,我能够找到并应用基于矩阵的解决方案。它基于可扩展为四线性插值的三线性插值的矩阵方法。维基百科提供了一篇关于 trilinear interpolation 的精彩文章。维基百科文章中的 8x8 矩阵可以扩展为 16x16 以进行四线性插值。下面写了一个函数来使矩阵中的每一行。
def quad_row(x, y, z, k):
"""
Generate a row for the quad interpolation matrix
x, y, z, k are scalar input values
"""
qrow = [1,
x, y, z, k,
x*y, x*z, x*k, y*z, y*k, z*k,
x*y*z, x*y*k, x*z*k, y*z*k,
x*y*z*k]
return qrow
很明显,这只是三线性矩阵内部行的扩展。该函数可以循环 16 次以生成整个矩阵。
旁注:如果你想变得更有趣,你可以使用 itertools 组合来完成 quad_row 功能。优点是您可以输入任意大小的数组,它 returns 是插值矩阵的格式正确的行。该功能更灵活,但最终速度较慢。
from itertools import combinations
def interp_row(values):
values = np.asarray(values)
n = len(values)
intp_row = [1]
for i in range(1, n+1):
intp_row.extend([np.product(x) for x in list(combinations(values, i))])
return intp_row
接受输入 table 的函数,找到接近您的插值的值,构建插值矩阵并执行矩阵数学运算,如下所示。
def quad_interp(values, table):
"""
values - four points to interpolate across, pass as list or numpy array
table - lookup data, four input columns and one output column
"""
table = np.asarray(table)
A, B, C, D, E = np.transpose(table)
a, b, c, d = values
in_vector = quad_row(a, b, c, d)
mask = (
np.isin(A, findclose(a, A)) &
np.isin(B, findclose(b, B)) &
np.isin(C, findclose(c, C)) &
np.isin(D, findclose(d, D)))
quad_matrix = []
c_vector = []
for row in table[mask]:
x, y, z, v, w = row
quad_matrix.append(quad_row(x, y, z, v))
c_vector.append(w)
quad_matrix = np.matrix(quad_matrix)
c_vector = np.asarray(c_vector)
a_vector = np.dot(np.linalg.inv(quad_matrix), c_vector)
return float(np.dot(a_vector, in_vector))
例如,调用函数如下所示。
df = pd.DataFrame({'liq':[5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 7000, 7000, 7000, 7000, 7000, 7000, 7000, 7000],
'visc':[10, 10, 10, 10, 30, 30, 30, 30, 10, 10, 10, 10, 30, 30, 30, 30],
'wlr':[0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5],
'gvf':[0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75],
'uncert':[0.0707, 0.0992, 0.0906, 0.1278, 0.0705, 0.0994, 0.091, 0.128, 0.0702, 0.0991, 0.0905, 0.1279, 0.0704, 0.0992, 0.0904, 0.1283],
})
values = [6532, 22, 0.412, 0.634]
quad_interp(values, df)
如上所示,上述函数不存在错误处理。如果尝试以下操作,它将崩溃:
1。 table 边界外的插值。
2。输入已在 table 中的查找值,导致选择的点少于 16 个。
此外,我承认以下几点:
1。命名约定可能会更好
2。可能存在创建掩码函数的更快方法
函数 findclose() 显示了原始问题。
如果您有任何反馈或改进空间,请告诉我。
我目前正在进行一个估算流量计不确定性的项目。仪表不确定性基于四个不同的值:
- 液体流量 (liq)
- 流体粘度 (cP)
- 水液比 (wlr)
- 气体体积分数 (gvf)
第三方为仪表提供 tables 的 liq、cP、wlr 和 gvf 的多个不同值。正如您所猜测的那样,来自仪表的数据永远不会完全落入预定义值之一。例如一分钟的数据可能是:
- 液体流量:6532
- 流体粘度:22
- 水液比:0.412
- 气体体积分数:0.634
使用以上数据对 table 执行四向插值以找出不确定性。
我想出了一个解决方案,但它看起来很笨拙,我想知道是否有人有任何想法。我还是 pandas 游戏的新手,非常感谢看到其他人的解决方案。
最初我对数据进行排序以将 table 降低到高于和低于我正在寻找的实际点的值。
aliq = 6532 # stbpd
avisc = 22 # centipoise
awlr = 0.412 # water liquid ratio
agvf = 0.634 # gas volume fraction
def findclose(num, colm):
arr = colm.unique()
if num in arr:
clslo = num
clshi = num
else:
clslo = arr[arr > num].min() # close low value
clshi = arr[arr < num].max() # close high value
return [clslo, clshi]
df = tbl_vx52[
(tbl_vx52['liq'].isin(findclose(aliq,tbl_vx52['liq']))) &
(tbl_vx52['visc'].isin(findclose(avisc,tbl_vx52['visc']))) &
(tbl_vx52['wlr'].isin(findclose(awlr,tbl_vx52['wlr']))) &
(tbl_vx52['gvf'].isin(findclose(agvf,tbl_vx52['gvf'])))
].reset_index(drop=True)
table 值从 2240 减少到 16。而不是包括所有数据 (tbl_vx52)。我创建了一些要加载的代码,因此您可以看到名为 df 的子数据框的外观,其中仅包含此示例区域上方和下方的值。
df = pd.DataFrame({'liq':[5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 7000, 7000, 7000, 7000, 7000, 7000, 7000, 7000],
'visc':[10, 10, 10, 10, 30, 30, 30, 30, 10, 10, 10, 10, 30, 30, 30, 30],
'wlr':[0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5],
'gvf':[0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75],
'uncert':[0.0707, 0.0992, 0.0906, 0.1278, 0.0705, 0.0994, 0.091, 0.128, 0.0702, 0.0991, 0.0905, 0.1279, 0.0704, 0.0992, 0.0904, 0.1283],
})
完成了一些非常粗略的循环,以开始根据各个输入(liq、visc、wlr 或 gvf)对值进行配对。下面显示的是 gvf.
上的第一个循环pairs = [
slice(0,1),
slice(2,3),
slice(4,5),
slice(6,7),
slice(8,9),
slice(10,11),
slice(12,13),
slice(14,15)]
for pair in pairs:
df.loc[pair,'uncert'] = np.interp(
agvf,
df.loc[pair,'gvf'],
df.loc[pair,'uncert']
)
df.loc[pair,'gvf'] = agvf
df = df.drop_duplicates().reset_index(drop=True)
删除了重复值,从 16 行减少到 8 行。然后再次为 wlr 重复此操作。
pairs = [
slice(0,1),
slice(2,3),
slice(4,5),
slice(6,7)
]
for pair in pairs:
df.loc[pair,'uncert'] = np.interp(
awlr,
df.loc[pair,'wlr'],
df.loc[pair,'uncert']
)
df.loc[pair,'wlr'] = awlr
df = df.drop_duplicates().reset_index(drop=True)
对 visc(四行)和最后的 liquid(两行)重复上述结构,直到子数组中只剩下一个值。这给出了您操作点的仪表不确定性。
我知道它很笨重。对不同方法的任何意见或想法表示赞赏。
好的,我能够找到并应用基于矩阵的解决方案。它基于可扩展为四线性插值的三线性插值的矩阵方法。维基百科提供了一篇关于 trilinear interpolation 的精彩文章。维基百科文章中的 8x8 矩阵可以扩展为 16x16 以进行四线性插值。下面写了一个函数来使矩阵中的每一行。
def quad_row(x, y, z, k):
"""
Generate a row for the quad interpolation matrix
x, y, z, k are scalar input values
"""
qrow = [1,
x, y, z, k,
x*y, x*z, x*k, y*z, y*k, z*k,
x*y*z, x*y*k, x*z*k, y*z*k,
x*y*z*k]
return qrow
很明显,这只是三线性矩阵内部行的扩展。该函数可以循环 16 次以生成整个矩阵。
旁注:如果你想变得更有趣,你可以使用 itertools 组合来完成 quad_row 功能。优点是您可以输入任意大小的数组,它 returns 是插值矩阵的格式正确的行。该功能更灵活,但最终速度较慢。
from itertools import combinations
def interp_row(values):
values = np.asarray(values)
n = len(values)
intp_row = [1]
for i in range(1, n+1):
intp_row.extend([np.product(x) for x in list(combinations(values, i))])
return intp_row
接受输入 table 的函数,找到接近您的插值的值,构建插值矩阵并执行矩阵数学运算,如下所示。
def quad_interp(values, table):
"""
values - four points to interpolate across, pass as list or numpy array
table - lookup data, four input columns and one output column
"""
table = np.asarray(table)
A, B, C, D, E = np.transpose(table)
a, b, c, d = values
in_vector = quad_row(a, b, c, d)
mask = (
np.isin(A, findclose(a, A)) &
np.isin(B, findclose(b, B)) &
np.isin(C, findclose(c, C)) &
np.isin(D, findclose(d, D)))
quad_matrix = []
c_vector = []
for row in table[mask]:
x, y, z, v, w = row
quad_matrix.append(quad_row(x, y, z, v))
c_vector.append(w)
quad_matrix = np.matrix(quad_matrix)
c_vector = np.asarray(c_vector)
a_vector = np.dot(np.linalg.inv(quad_matrix), c_vector)
return float(np.dot(a_vector, in_vector))
例如,调用函数如下所示。
df = pd.DataFrame({'liq':[5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 7000, 7000, 7000, 7000, 7000, 7000, 7000, 7000],
'visc':[10, 10, 10, 10, 30, 30, 30, 30, 10, 10, 10, 10, 30, 30, 30, 30],
'wlr':[0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5, 0.375, 0.375, 0.5, 0.5],
'gvf':[0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75, 0.625, 0.75],
'uncert':[0.0707, 0.0992, 0.0906, 0.1278, 0.0705, 0.0994, 0.091, 0.128, 0.0702, 0.0991, 0.0905, 0.1279, 0.0704, 0.0992, 0.0904, 0.1283],
})
values = [6532, 22, 0.412, 0.634]
quad_interp(values, df)
如上所示,上述函数不存在错误处理。如果尝试以下操作,它将崩溃:
1。 table 边界外的插值。
2。输入已在 table 中的查找值,导致选择的点少于 16 个。
此外,我承认以下几点:
1。命名约定可能会更好
2。可能存在创建掩码函数的更快方法
函数 findclose() 显示了原始问题。
如果您有任何反馈或改进空间,请告诉我。