change isotropize to satisfy Parseval's thereom

[miniufo]

This PR tries to solve #219.

Also, we add 3D detrending codes, which may be useful when calculating 3D (2D space and 1D time) frequency-wavenumber spectrum.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[roxyboy]

Thanks for the PR @miniufo !
I think changing the mean to sum makes sense. Could you run black on your code for formatting and also check whether the tests (here) pass by running pytest locally? It seems that the automatic Github workflow for tests isn't working at the moment...

[miniufo]

Hi @roxyboy, following your suggestions, I've reformatted the modified codes using black. But the assertion of all-close for isotropize at here:

xrft/xrft/tests/test_xrft.py

Lines 953 to 959 in 9ac46f9

	def _test_iso(theta):
	ps = xrft.power_spectrum(theta, spacing_tol=spacing_tol, dim=dims)
	ps = np.sqrt(ps.freq_x**2 + ps.freq_y**2)
	ps_iso = xrft.isotropize(ps, fftdim, nfactor=nfactor, truncate=truncate)
	assert len(ps_iso.dims) == 1
	assert ps_iso.dims[0] == "freq_r"
	npt.assert_allclose(ps_iso, ps_iso.freq_r**2 * 2 * np.pi, atol=0.02)

cannot be passed. The slopes of the two arrays are the same but differ by a constant that may depend on the number of points within each radial band, as sum is used instead of mean. Do you have any suggestion on changing the test codes?

[roxyboy]

Hi @roxyboy, following your suggestions, I've reformatted the modified codes using black. But the assertion of all-close for isotropize at here:

xrft/xrft/tests/test_xrft.py

Lines 953 to 959 in 9ac46f9

	def _test_iso(theta):
	ps = xrft.power_spectrum(theta, spacing_tol=spacing_tol, dim=dims)
	ps = np.sqrt(ps.freq_x**2 + ps.freq_y**2)
	ps_iso = xrft.isotropize(ps, fftdim, nfactor=nfactor, truncate=truncate)
	assert len(ps_iso.dims) == 1
	assert ps_iso.dims[0] == "freq_r"
	npt.assert_allclose(ps_iso, ps_iso.freq_r**2 * 2 * np.pi, atol=0.02)

cannot be passed. The slopes of the two arrays are the same but differ by a constant that may depend on the number of points within each radial band, as sum is used instead of mean. Do you have any suggestion on changing the test codes?

Apologies for my slow reply, I was away for the past week. Thanks for the Black formatting!

I understand that the assertion isn't passing because the radial wavenumber is taken using a mean but the spectrum is isotropized using sum. Since your PR is about making sure the isotropic spectra satisfy the Parseval's theorem, instead of trying to use the wavenumber, we should probably remove line 955 and change the assertion to something like npt.assert_allclose( (ps_iso*ps_iso.freq_r_diff).sum(), ps.sum() ).

[miniufo]

I've changed the assertion to npt.assert_allclose( (ps_iso).sum(), ps.sum() ). Since sum is used azimuthally in isotropize, another radial sum is enough to get the total energy. Adding an attribute ps_iso.freq_r.spacing would be much better, as there are ps.freq_x.spacing and ps.freq_y.spacing. But I find ps_iso.freq_r.diff('freq_r') is not a constant.

I also changed atol to 1 in test_isotropic_ps_slope:

xrft/xrft/tests/test_xrft.py

Line 1012 in 9ac46f9

npt.assert_allclose(a, s, atol=0.1)

Not sure if it is fine. It cannot get passed when atol=0.1.

[roxyboy]

Does this mean that this PR results in the slopes differing by a factor of one...? If so, this is not good.

[roxyboy]

I also changed atol to 1 in test_isotropic_ps_slope:

xrft/xrft/tests/test_xrft.py

Line 1012 in 9ac46f9

npt.assert_allclose(a, s, atol=0.1)

Not sure if it is fine. It cannot get passed when atol=0.1.

Yes, this is what I was personally worried that could happen when switching the isotropization from mean to sum. In a discrete (and not continuous) system as we have for data analyses, the number of points within each radial band differs. A mean isn't sensitive to this but a sum is. The spectral slopes disagreeing by a factor of one is not going to be acceptable unfortunately...

[miniufo]

Hi @roxyboy, I just looked into the test, and found that the frequency range to be fitted is crucial.

xrft/xrft/tests/test_xrft.py

Line 1011 in 9ac46f9

y_fit, a, b = xrft.fit_loglog(iso_ps.freq_r.values[4:], iso_ps.values[4:])

Previously, indices from 4 to end are used, which include the highest frequencies that are steeper than expected. Now I switch the indices from 0 to -35, which overlap with the expected range, and this test passed even with a smaller `atol=0.06.

The results looks like:

The red curve is the fitted line over indices from 0 to -35. You may want to check if this is OK.

[roxyboy]

@miniufo Ok, this is great! Can you run the Black formatting once more and I can merge this :)

[miniufo]

Oops, I forgot to do this after any changes. Can you check it again?

[roxyboy]

@miniufo Sorry, it seems that some files are still missing the Black formatting...

[miniufo]

That's strange. I format all .py files. Possibly missing the .ipynb files. Please try again and see if it is fine now.

[roxyboy]

@miniufo I agree that this is kind of unexpected behavior but the GitHub log is saying that test_xrft-checkpoint.py is missing the formatting...

> Run black --check xrft
  black --check xrft
  shell: /usr/bin/bash -e {0}
  env:
    pythonLocation: /opt/hostedtoolcache/Python/3.8.18/x64
    LD_LIBRARY_PATH: /opt/hostedtoolcache/Python/3.8.18/x64/lib
would reformat xrft/tests/.ipynb_checkpoints/test_xrft-checkpoint.py

[roxyboy]

Thanks for your work @miniufo ! The checks are all passing now so I will merge this.

[miniufo]

That's great~