使用口罩时长度不匹配的问题

问题描述 投票:0回答:1

我正在编写代码,并且有一个函数可以用满足条件的值来计算不满足条件的值,但我在管理数组的形状方面遇到了很多麻烦。

我有一个类似的函数,但具有其他逻辑结构来执行此操作(MWE 表示有效的函数)

import numpy as np
def f_est(f, fd, mask):
    """
    This function returns the estimated value given a function f and its derivative fd 
    for the points selected by the mask, estimating forward.
    
    Parameters:
    - f: Array of function values.
    - fd: Array of derivative values of the function.
    - mask: Boolean mask to select points in f and fd.
    """
    h = 0.000001
    
    # Find the last index that satisfies the mask
    last_index = np.max(np.where(mask)[0])
    
    # Create shifted masks (inspired by f_cal), but centered on the last index
    mask_current = mask[:last_index + 1]  # Mask for current positions up to the last index
    mask_prev = mask[:last_index]          # Mask for previous positions (for fd_prev_slice)
    mask_prev2 = mask[:last_index - 1]     # Mask for previous positions (for fd_prev2_slice)
    
    # Apply masks to f and fd (with shifts), centered on the last index
    f_slice = f[:last_index + 1][mask_current]  # Note: adjusted to align with mask_current
    fd_slice = fd[:last_index + 1][mask_current] 
    fd_prev_slice = fd[:last_index][mask_prev]  
    fd_prev2_slice = fd[:last_index - 1][mask_prev2]  
    
    # Perform the calculations with consistent slices, estimating forward
    # Use the last value of f_slice, fd_slice, fd_prev_slice, and fd_prev2_slice for estimation
    last_f = f_slice[-1]
    last_fd = fd_slice[-1]
    last_fd_prev = fd_prev_slice[-1] if len(fd_prev_slice) > 0 else 0
    last_fd_prev2 = fd_prev2_slice[-1] if len(fd_prev2_slice) > 0 else 0
    
    estimated_next_value = (
        last_f 
        + h * last_fd 
        + 1 / 2 * (h * last_fd - h * last_fd_prev) 
        + 5 / 12 * ((h * last_fd - h * last_fd_prev) - (h * last_fd_prev - h * last_fd_prev2))
    )
    
    return estimated_next_value

f = np.array([1, 2, 3, 4, 5, 6, 7])
fd = f 
mask = np.array([True, True, True, False, False, False, False])
print("Original Array:", f)
print("Length of Original Array:", len(f))
print("Masked Array:", f[~mask])
print("Length of Masked Array:", len(f[~mask]))
f[~mask] = f_est(f, fd, mask)

print("Final Array:", f)
print("Length of Final Array:", len(f))

但是有了这个功能(MWE 不起作用):

import numpy as np
def f_cal(var, dvar, mask):
    """
    Calculates next value using trapezoidal method with masks, estimating forward.
    
    Parameters:
        var: array of current values
        dvar: array of derivatives 
        mask: boolean array indicating which positions to calculate
    """
    h = 0.0000001
    # Encontrar el último índice que verifica la máscara
    last_index = np.max(np.where(mask)[0])
    
    # Crear máscaras para posiciones actuales y la siguiente
    mask_current = mask[:last_index]  
    mask_next =  mask[1:last_index+1] # Marcar como True el índice siguiente
    
    # Ajustar los arreglos para alinear con las máscaras
    var_current = var[:last_index+1][mask_current]
    dvar_current = dvar[:last_index+1][mask_current]
    dvar_next = dvar[:last_index+2][mask_next][:1]  # Solo el valor siguiente
    
    # Calculate using trapezoidal method with masks, estimating forward
    result = var_current + h * dvar_next - 1/2*(h*dvar_next-h*dvar_current[-1])
    
    return result

f = np.array([1, 2, 3, 4, 5,6,7])
fd = f 
mask = np.array([True, True, True, False, False, False, False])
print("Original Array:", f)
print("Length of Original Array:", len(f))
print("Masked Array:", f[~mask])
print("Length of Masked Array:", len(f[~mask]))
f[~mask] = f_cal(f, fd, mask)

print("Final Array:", f)
print("Length of Final Array:", len(f))

我在保持数组的长度与不满足条件的元素数量相匹配方面遇到了很多麻烦

arrays numpy numerical-methods numpy-slicing
1个回答
0
投票

我相信您正在尝试使用 泰勒级数 来近似超出最后一个已知值的位置处的函数。函数值存储在

var
数组中,一阶导数存储在
dvar
数组中。看起来您正在使用
mask
来标识最后一个已知的函数值,而
h
代表参数步骤。

从您的代码来看,您似乎正在使用直至二阶导数的项,您将其近似为相邻一阶导数之间的差除以参数步骤。这可以表示为:

f(x+h) = f(x) + h⋅f′(x) + ½⋅h²⋅f″(x) + … ≈
       ≈ f(x) + h⋅f′(x) + ½⋅h²⋅(f′(x+h)-f′(x))/h

至少,这与您作为结果返回的表达式产生共鸣(我稍微调整了):

var_current + h*dvar_next - 1/2*(h*dvar_next-h*dvar_current)

我们可以将上面的逐步改造如下:

 var_current + h*dvar_next - 1/2*h*dvar_next + 1/2*h*dvar_current
== var_current + 1/2*h*dvar_next + 1/2*h*dvar_current
== var_current + h*dvar_current + 1/2*h*dvar_next - 1/2*h*dvar_current
== var_current + h*dvar_current + 1/2*h*(dvar_next - dvar_current)
== var_current + h*dvar_current + 1/2*(h**2)*(dvar_next - dvar_current)/h

最后一行让我想起了泰勒级数,这就是为什么我怀疑你的函数描述中提到的梯形方法可能是一个错误。

如果是这样的话,函数代码可能如下所示:

def f_cal(var, dvar, mask, h=1e-7):
    index = np.asarray(mask).nonzero()[0]
    assert index and index[-1] < len(mask)-1
    pos = index[-1]
    f = var[pos]
    df = dvar[pos]
    df_next = dvar[pos+1]
    return f + h*df + h*(df_next - df)/2

但是,我可能误解了你的意图。例如,如果您故意引用梯形方法,或者您使用不同的约定,则您可能需要澄清您的问题。


附注回答有关使用掩码时长度不匹配的原始问题:您似乎在某些地方错过了处理数组维度。这是一个可供比较的替代实现:

def f_cal(var, dvar, mask):
    h = 0.0000001
    last_index = np.max(np.where(mask)[0])
    mask_current = mask[:last_index+1]  
    mask_next =  mask[:last_index+2]
    var_current = var[:last_index+1][mask_current][-1]
    dvar_current = dvar[:last_index+1][mask_current][-1]
    dvar_next = dvar[:last_index+2][mask_next][-1]
    return var_current + h*dvar_next - 1/2*h*(dvar_next-dvar_current)
© www.soinside.com 2019 - 2024. All rights reserved.