Rust转Python卷积数组计算疑问：数组形状与数值溢出问题

我正在将某公式的Rust实现迁移至Python环境，参考公式对应示意图：公式示意图。针对卷积计算流程，有两个核心疑问：

1. 卷积数组的形状如何确定？ #

我在代码中设定convolve[[m,j]]的形状为**(时间步数x可用生产井数)**，但对这个设定的正确性存疑，想明确这类卷积数组的形状应该如何推导确定。

2. 计算结果为何出现极大值？是否与数组形状错误相关？ #

计算结果片段如下：

# Results snippet
[[5.68561108e+001 5.61980823e+001 5.61675861e+001 5.61661613e+001]
 [2.35307528e+002 2.34679642e+002 2.34608585e+002 2.34637574e+002]
 ...
 [4.99748565e+293 4.99665753e+293 4.99716953e+293 4.99634444e+293]]

这些结果将作为损失函数输入，通过scipy.minimize做优化。使用np.mean()计算均方误差时，因数值过大产生NaN，直接导致优化失败。我想知道出现极大值的原因，以及这是否和问题1中的数组形状错误有关。

我的Rust实现代码 #

fn proxy_crm(_py: Python, m: &PyModule) -> PyResult<()> {
    fn q_prim(prod: ArrayView2<'_, f64>, time: ArrayView1<'_, f64>, lambda_prod: ArrayView1<'_, f64>, tau_prim: ArrayView1<'_, f64>) -> Array1<f64> {
        let n_prod: usize = prod.raw_dim()[1];
        let mut result: Array1<f64> = Array1::zeros([n_prod]);

        for j in 0..n_prod {
            let time_decay = (-&time / tau_prim[j]).mapv(f64::exp);
            result[j] = time_decay[j] * prod[[0,j]] * lambda_prod[j]
        }
        result
    }

    fn q_crm(inj: ArrayView2<'_, f64>, time: ArrayView1<'_, f64>, lambda_ip: ArrayView2<'_, f64>, tau: ArrayView1<'_, f64>, mask: ArrayView2<'_, f64> ) -> Array2<f64> {
        let n_t: usize = time.raw_dim()[0];
        let n_inj: usize = lambda_ip.raw_dim()[1];
        let n_prod: usize = mask.raw_dim()[1];
        let mut convolve: Array2<f64> = Array2::zeros([n_t, n_prod]);
        
        let lambda_ip_mod = calc_sh_mask(lambda_ip, mask);

        for i in 0..n_inj {
            for j in 0..n_prod {
                convolve[[0,j]] = (1.0 - ((time[0] - time[1]) / tau[j]).exp()) * lambda_ip_mod[[0,j,i]] * inj[[0,i]];
                for m in 1..n_t {
                    for n in 1..m+1 {
                        let time_decay = (1.0 - ((time[m-1] - time[m]) / tau[j]).exp()) * ((time[n] - time[m]) / tau[j]).exp();

                        convolve[[m,j]] += time_decay * lambda_ip_mod[[m,j,i]] * inj[[m,i]];
                    }
                }
            }
        }
        convolve
    }

    fn q_bhp(time: ArrayView1<'_, f64>, tau: ArrayView1<'_, f64>, press: ArrayView2<'_, f64>, prod_index: ArrayView1<'_, f64>) -> Array2<f64> {
        let n_t = time.raw_dim()[0];
        let n_prod = press.raw_dim()[1];
        let mut convolve: Array2<f64> = Array2::zeros([n_t,n_prod]);
        
        for j in 0..n_prod {
            convolve[[0,j]] = (1.0 - ((time[0] - time[1]) / tau[j]).exp()) * prod_index[0] * tau[j] * (press[[1,0]]-press[[0,0]]) / (time[1]-time[0]);
            for m in 1..n_t {
                for n in 1..m+1 {
                    let time_decay = (1.0 - ((time[m-1] - time[m]) / tau[j]).exp()) * ((time[n] - time[m]) / tau[j]).exp();
                    let delta_bhp = press[[m-1,j]] - press[[m,j]];

                    convolve[[m,j]] += time_decay * prod_index[j] * tau[j] * (delta_bhp / (time[m-1]-time[m]));
                }
            }
        }
        convolve
    }

    fn sh_mask(prod: ArrayView2<'_, f64>) -> Array2<f64> {
        let n_prod: usize = prod.raw_dim()[1];
        let n_t: usize = prod.raw_dim()[0];
        let mut mask = Array2::zeros([n_t, n_prod]);

        for t in 0..n_t {
            for j in 0..n_prod {
                if prod[[t,j]] == 0.0 {
                    mask[[t,j]] = 1.0;
                }
            }
        }
        mask
    }

    fn calc_sh_mask(lambda_ip:ArrayView2<'_, f64>, mask: ArrayView2<'_, f64>) -> Array3<f64> {
        let n_prod: usize = mask.raw_dim()[1];
        let n_inj: usize = lambda_ip.raw_dim()[1];
        let n_t: usize = mask.raw_dim()[0];
        
        let tensor: Array3<f64> = Array3::ones([n_t,n_prod,n_inj]);
        let mut lambda_ip_result: Array3<f64> = Array3::zeros([n_t,n_prod,n_inj]);

        for t in 0..n_t {
            let mut sum_lambda_sh: Array1<f64> = Array1::zeros([n_inj]);
            for j in 0..n_prod {
                for i in 0..n_inj {
                    lambda_ip_result[[t,j,i]] = f64::abs(f64::abs((mask[[t,j]]-1.0) * lambda_ip[[j,i]]) * (tensor[[t,j,i]] + tensor[[t,j,i]] * mask[[t,j]]));
                    if mask[[t,j]] == 1.0 {
                        sum_lambda_sh[i] += lambda_ip[[j,i]];
                        for k in 0..n_prod {
                            if k!=j {
                                lambda_ip_result[[t, k, i]] += f64::abs(f64::abs(lambda_ip[[k,i]]) * (tensor[[t,j,i]] + sum_lambda_sh[i] * tensor[[t,j,i]] * mask[[t,j]]));
                            }
                        }
                    }
                }
            }
        }
        lambda_ip_result
    }
}

上述代码参考了GitHub某仓库，可按需提供链接。

内容的提问来源于stack exchange，提问作者tanosheet