Python中多级线性插值的高效方法

Efficient method for multi-level linear interpolation in Python

我目前正在进行一个估算流量计不确定性的项目。仪表不确定性基于四个不同的值:

  1. 液体流量 (liq)
  2. 流体粘度 (cP)
  3. 水液比 (wlr)
  4. 气体体积分数 (gvf)

第三方为仪表提供 tables 的 liq、cP、wlr 和 gvf 的多个不同值。正如您所猜测的那样,来自仪表的数据永远不会完全落入预定义值之一。例如一分钟的数据可能是:

  1. 液体流量:6532
  2. 流体粘度:22
  3. 水液比:0.412
  4. 气体体积分数: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() 显示了原始问题。

如果您有任何反馈或改进空间,请告诉我