Python RuntimeWarning 在 double_scalars 中遇到溢出

Question

我正在尝试将“C”代码移植到Python，这就是我编码的内容:

scale = 1.0 / (rows * cols)
RemoveConstantBiasFrom(rarray, scale)
zarray = rarray[:]
zarray = DCT(zarray, rows, cols)

zarray = zarray.flatten()

beta = np.dot(rarray, zarray)

if iloop == 0:
    parray = zarray[:]
else:
    btemp = beta / beta_prev
    parray = zarray + btemp * parray
RemoveConstantBiasFrom(parray, scale)
beta_prev = beta
for j in range(rows):
        for i in range(cols):
            k = j * rows + i

            k1 = k + 1 if i < cols - 1 else k - 1
            k2 = k - 1 if i > 0 else k + 1
            k3 = k + cols if j < rows - 1 else k - cols
            k4 = k - cols if j > 0 else k + cols


            w1 = SIMIN(wts[k], wts[k1])
            w2 = SIMIN(wts[k], wts[k2])
            w3 = SIMIN(wts[k], wts[k3])
            w4 = SIMIN(wts[k], wts[k4])

            A = (w1 + w2 + w3 + w4) * parray[k]
            B = (w1 * parray[k1] + w2 * parray[k2])
            C = (w3 * parray[k3] + w4 * parray[k4])

            zarray[k] = A - (B + C)

wts 和 parray 是两个具有 262144 个值(512x512 矩阵)的展平数组。当我执行此操作时，我得到:

test.py:152: RuntimeWarning: overflow encountered in double_scalars
  B = (w1 * parray[k1] + w2 * parray[k2])
test.py:154: RuntimeWarning: overflow encountered in double_scalars
  C = (w3 * parray[k3] + w4 * parray[k4])
test.py:155: RuntimeWarning: overflow encountered in double_scalars
  zarray[k] = A - B + C

所以，我开始“调试”代码。

1) 我在循环之前将 print 放入 max(parray) 中，得到:

Maximum value for parray  15.2665322926

2) 我在循环内添加了一条 if 语句来观察 parray[k1] 行为 if parray[k1] > maxparray: print k1 并且我获得了很多“k1”:

...
251902
252414
252926
253438
253950
254462
254974
255486
255998
256510
257022
257534
258046
258558
259070
259582
260094
260606
261118
261630
262142

那么，问题是:如果我从未更改过“parray”，为什么我会得到超过 max(parray) 的不同值？

这是 C 代码:

  int     i, j, k;
  int     k1, k2, k3, k4;
  double  sum, w1, w2, w3, w4, btemp, delta, avg, scale;
  float   *wts;

  scale = 1.0/(xsize*ysize);
  /* remove constant bias from rarray */
  for (k=0, avg=0.0; k<xsize*ysize; k++)  avg += rarray[k];
  avg *= scale;
  for (k=0; k<xsize*ysize; k++)  rarray[k] -= avg;
  /* compute cosine transform solution of Laplacian in rarray */
  for (k=0; k<xsize*ysize; k++)  {
    zarray[k] = rarray[k];
  }
  DCT(zarray, xsize, ysize, xcos, ycos);
  /* calculate beta and parray */
  for (k=0, *beta=0.0; k<xsize*ysize; k++) {
    *beta += rarray[k]*zarray[k];
  }
  printf("beta = %lf\n", *beta);
  if (iloop == 0) {
    for (k=0; k<xsize*ysize; k++) {
      parray[k] = zarray[k];
    }
  }
  else {
    btemp = (*beta)/(*beta_prev);
    for (k=0; k<xsize*ysize; k++) {
      parray[k] = zarray[k] + btemp*parray[k];
    }
  }
  /* remove constant bias from parray */
  for (k=0, avg=0.0; k<xsize*ysize; k++)  avg += parray[k];
  avg *= scale;
  for (k=0; k<xsize*ysize; k++)  parray[k] -= avg;
  *beta_prev = *beta;
  /* calculate Qp */
  for (j=0; j<ysize; j++) {
    for (i=0; i<xsize; i++) {
      k = j*xsize + i;
      k1 = (i<xsize-1) ? k + 1 : k - 1;
      k2 = (i>0) ? k - 1 : k + 1;
      k3 = (j<ysize-1) ? k + xsize : k - xsize;
      k4 = (j>0) ? k - xsize : k + xsize;
      if (dxwts==NULL && dywts==NULL) {  /* unweighted */
        w1 = w2 = w3 = w4 = 1.0;
      }
      else if (dxwts==NULL || dywts==NULL) {  /* one set of wts */
        wts = (dxwts) ? dxwts : dywts;
        w1 = SIMIN(wts[k], wts[k1]);
        w2 = SIMIN(wts[k], wts[k2]);
        w3 = SIMIN(wts[k], wts[k3]);
        w4 = SIMIN(wts[k], wts[k4]);
      }
      else {    /* dxwts and dywts are both supplied */
        w1 = dxwts[k];
        w2 = (i>0) ? dxwts[k-1] : dxwts[k];
        w3 = dywts[k];
        w4 = (j>0) ? dywts[k-xsize] : dywts[k];
      }

      zarray[k] = (w1 + w2 + w3 + w4)*parray[k]
                        - (w1*parray[k1] + w2*parray[k2] 
                                + w3*parray[k3] + w4*parray[k4]);
    }
  }

我只是传递 dxwts (一组权重)，因此不需要其他权重，这就是为什么我使用一个 if 语句，其中 SIMIN 是:

#define SIMIN(x,y) (((x)*(x) < (y)*(y)) ? (x)*(x) : (y)*(y))

Answer 1

最有可能的是，parray 和 zarray 是同一个数组。

如果你这样做:

zarray = DCTF()
parray = zarray

那么 zarray 和 parray 是同一个数组，并且更改 zarray 中的元素将在 use 中可见parray。 (这与 C 没有显着差异。)

如果您想要副本，可以这样做:

zarray = DCTF()
parray = zarray[:]