Rotational invariance of spatial scores is not present for BMSE

[sidekock]

During our efforts to convert the pysteps datamodel to xarray, we needed to change tons of tests.

One of the things we know can happen in xarray is that if you take the values out of the dataset, the resulting numpy array can flip compared to when you keep the metadata. We noticed that this flip made the test of test_verification_spacialscores.py fail. When flipping the data over the y axis, we do indeed get the correct score. It seems strange to us that this score does not seem to be the same under rotational invariance (where all fields are rotated similarly, of course). Since we don't know the details of "BMSE" we want to check if our assumption of rotational invariance are valid?

[sidekock]

A long chat with ChatGPT gave me some solutions but I don't know how needed this solution is:

How to make it (almost) invariant

If you want invariance to global flips (and shifts), use a transform that doesn’t decimate or that handles boundaries in a flip-friendly way:

Use the stationary/undecimated wavelet transform: pywt.swt2 (or pywt.dwtn/idwtn with no decimation). SWT is (near) shift/flip-invariant. You’d replace _wavelet_decomp to use swt2 and compute per-level energies from the undecimated detail images directly.

Or force periodic boundaries: call wavedec2(..., mode='periodization') and waverec2(..., mode='periodization'). Periodization reduces boundary asymmetries. It still downsamples, so you may see small differences, but they’re usually much smaller.

Ensure sizes are even at all levels: crop/pad so m and n are multiples of 2**J (max level). That reduces alignment artifacts from leftovers at the borders.

J = pywt.dwt_max_level(min(m, n), pywt.Wavelet(wavelet).dec_len)
m2 = (m // (2**J)) * (2**J)
n2 = (n // (2**J)) * (2**J)
X = X[:m2, :n2]