Python RuntimeWarning overflow encountered in double_scalars

Question

I'm trying to port a "C" code to Python, this is what I coded:

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 and parray are two flattened arrays with 262144 values (a 512x512 matrix). When I'm perfoming this operation I get:

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

So, I started to "debug" the code.

1) I put a print to max(parray) before the loops and I obtained:

Maximum value for parray  15.2665322926

2) I added an if statement inside the loop to watch parray[k1] behavior if parray[k1] > maxparray: print k1 and I obtained a lot of "k1's":

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

So, the question: If I never changed "parray" why I'm getting different values exceeding max(parray)?

This is the C code:

  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]);
    }
  }

I'm just passing dxwts (one set of weights) so the others are not needed, that's why I'm using one if statement, where SIMIN is:

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

Answer 1

The problem is the way I'm copying the arrays, the correct way is:

zarray[:] = rarray

Answer 2

Most likely is that parray and zarray are the same array.

If you do this:

zarray = DCTF()
parray = zarray

then zarray and parray are the same array and changing an element in zarray will be visible in uses of parray. (This is not significantly different from C.)

If you want a copy, you can do this:

zarray = DCTF()
parray = zarray[:]