New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Functions to caluclate the radial average power spectrum #303

Merged
merged 13 commits into from Nov 23, 2016

Conversation

Projects
None yet
2 participants
@santis19
Member

santis19 commented Jul 25, 2016

Fixes #224

I've retaken the ideas of the closed #227 and #230.

This time I've written two new functions in gravmag.transform: a simple one that calculates the Power-Density Spectra (PDS) of regular gridded data and the other one radially averages this PDS (or any real 2D-array in the wavenumber domain).

As discussed in #230, this new functions has the following characteristics:

  • The power_density_spectra() function makes use of the _fftfreq() function.
  • The inputs for the radial_average() function are the wavenumbers array (kx and ky) and the pds (or any wavenumber-domain array). The outputs of the power_density_spectra() function are suitable to use as inputs of this one.
  • The are no restrictions of squared x, y grids or equally spaced.
  • All the calculations are made with the kx and ky distances instead of using indexes.
  • The code is more numpy-ish, making the calculation even faster.
  • It allows to optionally change the rings' width (making the averages more or less populated) and changing the radius of the biggest ring (in order to leave highest wavenumbers out of the averaging, or in case that our ky, ky grid isn't square to introduce those points that are left outside by default).

Now the schematic draw looks like this:
rings

I leave an example with its outputs:

import numpy as np
import matplotlib.pyplot as plt
from fatiando.datasets import fetch_bouguer_alps_egm
from fatiando.gridder import load_surfer
from fatiando.gravmag import transform


# Fetch Alps Bouguer Anomaly Data
fetch_bouguer_alps_egm(fname="alps-bouguer.grd")
x, y, bg, shape = load_surfer("alps-bouguer.grd")

plt.contourf(y.reshape(shape), x.reshape(shape), bg.reshape(shape), 150)
plt.axes().set_aspect('equal')
plt.title("Bouguer Anomaly of the Alps Region")
plt.colorbar()
plt.savefig("bouguer.png")
plt.show()

# Calculate the Power Density Spectra
kx, ky, pds = transform.power_density_spectra(x, y, bg, shape)

plt.contourf(ky.reshape(shape), kx.reshape(shape), np.log(pds.reshape(shape)),
             150)
plt.axes().set_aspect('equal')
plt.title("Power Density Spectra of the Bouguer Anomaly")
plt.colorbar()
plt.savefig("pds.png")
plt.show()


# Radially Average the Power Density Spectra
default_width = max(kx[1, 0], ky[0, 1])
default_radius = min(kx.max(), ky.max())

# Default
k_radial, pds_radial = transform.radial_average(kx, ky, pds)
plt.plot(k_radial, np.log(pds_radial), 'o-')
plt.title("Radially Averaged PDS")
plt.savefig("default.png")
plt.show()

# Wider Rings
k_radial, pds_radial = transform.radial_average(kx, ky, pds,
                                                ring_width=4*default_width)

plt.plot(k_radial, np.log(pds_radial), 'o-')
plt.title("Radially Averaged PDS (wider rings)")
plt.savefig("wider_rings.png")
plt.show()

# Bigger Max Radius
k_radial, pds_radial = transform.radial_average(kx, ky, pds,
                                                max_radius=0.5*default_radius)

plt.plot(k_radial, np.log(pds_radial), 'o-')
plt.title("Radially Averaged PDS (max radius reduced to half)")
plt.savefig("max_radius.png")
plt.show()

bouguer
pds
default
wider_rings
max_radius

Checklist:

  • Make tests for new code (at least 80% coverage)
  • Create/update docstrings
  • Include relevant equations and citations in docstrings
  • Docstrings follow the style conventions
  • Code follows PEP8 style conventions
  • Code and docs have been spellchecked
  • Include new dependencies in doc/install.rst, requirements.txt, environment.yml, ci/requirements-conda.txt and ci/requirements-pip.txt.
  • Documentation builds properly (run cd doc; make locally)
  • Changelog entry (leave for last to avoid conflicts)
  • Firt-time contributor? Add yourself to doc/contributors.rst (leave for last to avoid conflicts)

@leouieda leouieda added this to the 0.6 milestone Jul 25, 2016

@leouieda

This comment has been minimized.

Member

leouieda commented Jul 25, 2016

@santis19 thank you for opening this again!

I've taken a quick look at the code and it's very nice and clean! There isn't much to add to the code. I'll leave some comments inline with some suggestions to improve the interface and docstrings. Great job!

A nice thing about the separation is that we can now test the radial_average function on simpler input. I'd suggest testing by creating some data that is made up of rings of integers. Then, if we specify the correct ring_width, we should get the actual integer values as output.

Another test could then be to create a data array that is the distance from the center of the grid. The radial average should be the distance returned to us.

These two test cases should cover most of the code and usage of this function.

I'm having a harder time thinking how we could test the power_density_spectra function. Maybe we could test it on a sin or cos oscillating in one dimension only? Then go and test it on the other dimension? The tests should also check is kx and ky are reasonable. This kind of thing is much harder to test (when you don't really know the output).

Sidetrack

Another point I'm considering is if it would be better to make a sane_fft2d function instead of a power spectra function (or maybe both). The sane_fft2d function would calculate the FFT and return it and the wavenumber arrays. It should also pad the data using the gridder padding functions by @drandykass (optionally), since this is done in most FFT applications anyway. We can offer the power spectra function function as well that takes the space-domain data and uses the sane_fft2d function, much like tga does for the derivatives.

I like this approach better because it gets rid of repeated code in the transform module and would help people develop other functions for FFT processing. But this should be kept for another PR 😄 I'll open an issue for this so that we can keep track.

End sidetrack

For this PR, I think you should focus now on implementing the tests discussed above. They should go in the fatiando/gravmag/tests/test_transform.py file.

Let me know if you need help with these. I recently updated the reduction to the pole tests and the tilt tests are brand new so you should use those as inspiration.

@leouieda

This comment has been minimized.

Member

leouieda commented Jul 25, 2016

@mtb-za and @drandykass would you care to weight in on this as well?

return kx, ky, pds
def radial_average(kx, ky, pds, max_radius=None, ring_width=None):

This comment has been minimized.

@leouieda

leouieda Jul 25, 2016

Member

This function is generic to any 2D array, so it might be better to rename the inputs to x, y, data instead of kx, ky, pds to indicate this. We can offer an example in the docstring to make it clear how to use this to calculate the power spectra.

This comment has been minimized.

@santis19

santis19 Jul 26, 2016

Member

First I thought to make it explicitly generic to any 2D array, but I noticed that the kx and ky arrays (as outputs of _fftfreq) has special properties:

  • they always include the origin (kx, ky) = (0, 0)
  • they are ordered in special way (first the positive and then the negatives)
  • if the original array has even amount of points per axes then the kx and ky arrays are not symmetric to the origin, instead they have one extra point in the positive semiaxe.

So I decided to stick to frequency domain (or wavenumber domain) 2D arrays only. Maybe we can generalize this function. In order to do that I think that we will have to overcome two possible problems:

  • How the distances \Delta k_x and \Delta k_y are calculated.
  • What would happen if the x and y arrays don't include the origin (that would lead to a nan in the radial_pds array)?

This comment has been minimized.

@leouieda

leouieda Jul 27, 2016

Member

OK, I hadn't thought of that. You're right, this has to be specific to something returned by FFT.

I don't think it's worth it trying to make this overly general. It will be a lot of work and there are many subtle ways that bugs can hide in the code. It's better to make it clear that this function only operates on Fourier transforms.

Since this is the case, then it might be better to rename the function to radial_average_spectrum to make it clear that it won't work for general arrays. However, it will work for other quantities in the frequency domain, not only the power spectrum. So removing explicit mentions of it should still be done.

def radial_average(kx, ky, pds, max_radius=None, ring_width=None):
r"""
Calculates the average of the Power Density Spectra points that falls

This comment has been minimized.

@leouieda

leouieda Jul 25, 2016

Member

Same as above.

def radial_average(kx, ky, pds, max_radius=None, ring_width=None):
r"""
Calculates the average of the Power Density Spectra points that falls
inside concentric rings around the zero wavenumber. By default, the

This comment has been minimized.

@leouieda

leouieda Jul 25, 2016

Member

Same as above for "zero wavenumber". Maybe "origin" is a better term?

inside concentric rings around the zero wavenumber. By default, the
biggest ring is the bigger one to have at least one point of both axes,
and the rings' width is set as the maximum of :math:`\Delta k_x` and
:math:`\Delta k_y`, being them the equidistances of the kx and ky arrays.

This comment has been minimized.

@leouieda

leouieda Jul 25, 2016

Member

The defaults should be explained (and are) in the description of each parameter. Here you should explain generally what the function does (what is the method/algorithm used and what you should know when using this function). This should be on a separate paragraph.

This comment has been minimized.

@leouieda

leouieda Jul 25, 2016

Member

The purpose of this text is so that the user understands what you are talking about below when you mention "rings" and "the biggest ring" etc.

This comment has been minimized.

@leouieda

leouieda Jul 25, 2016

Member

The text should make it clear that this only works if the grid is centered on the origin, right? Or did I get this wrong?

respectively.
.. note:: It's recommended to use the outputs of the
power_density_spectra function as input of this one.

This comment has been minimized.

@leouieda

leouieda Jul 25, 2016

Member

Maybe the note could say "To calculate the radially averaged power spectra, use the outputs..."?

* shape : tuple = (nx, ny)
The shape of the grid
* max_radius : float (optional)
Changes the default inner radius of the biggest ring.

This comment has been minimized.

@leouieda

leouieda Jul 25, 2016

Member

It's better to explain the parameters in the order 1) say what this is (Inner radius of the biggest ring) 2) what default is (The largest value of x.max() and y.max()) then 3) what changing this does to the output and recommend a good value.

@santis19

This comment has been minimized.

Member

santis19 commented Jul 26, 2016

I've just finished the radial_average tests: one with integers and another one with distances (as @leouieda proposed).

I have a problem with the last one: the first ring (not the one related to k=0, but the next one) has more than 0.1 relative error. I think that may be related to the fact that this ring is the least populated, so the values are not so randomly spread inside it. But for the rest of the array the test doesn't fail.

I leave the codes here. I don't know if I should start another PR for this test or just push the edited test_transform.py into my branch...please help this github newbie! haha

from __future__ import division
import numpy as np
import numpy.testing as npt
from fatiando.gravmag import transform
from fatiando import gridder


def test_radial_average_distances():
    shape = (201, 201)
    area = (-100, 100, -100, 100)
    x, y = gridder.regular(area, shape)
    x, y = x.reshape(shape), y.reshape(shape)
    x, y = np.fft.ifftshift(x), np.fft.ifftshift(y)
    distances = np.sqrt(x**2 + y**2)
    k, radial_distances = transform.radial_average(x, y, distances)
    npt.assert_allclose(k, radial_distances, rtol=0.1)  # this fails
    npt.assert_allclose(k[2:], radial_distances[2:], rtol=0.1)  # this doesn't


def test_radial_average_integers():
    x = np.arange(-10, 11, 1)
    y = np.arange(-10, 11, 1)
    y, x = np.meshgrid(y, x)
    x, y = np.fft.ifftshift(x), np.fft.ifftshift(y)
    r = np.sqrt(x**2 + y**2)
    z = np.zeros(x.shape)
    integers = np.arange(0, x.max() + 1, 1)
    for i in integers:
        if i == 0:
            inside = r <= 0.5
        else:
            inside = np.logical_and(r > i - 0.5, r <= i + 0.5)
        z[inside] = i
    k, radial_z = transform.radial_average(x, y, z)
    npt.assert_allclose(integers, radial_z, rtol=0.1)

@santis19 santis19 force-pushed the santis19:pds_radial_average branch from d9bf783 to 39b7bb2 Jul 26, 2016

@santis19

This comment has been minimized.

Member

santis19 commented Jul 26, 2016

I've edited the description of the radial_average function following @leouieda suggestions.

@leouieda

This comment has been minimized.

Member

leouieda commented Jul 27, 2016

@santis19 you can push the test here as well. A PR will only be merged when tests are implemented and passing, that's a way to guarantee that all new code works properly. You'll notice that one of the reasons why "All checks have failed" is that the code test coverage decreased (meaning that you added untested code). Just avoid pushing things that do have anything to do with these two functions (typo fixes elsewhere, bugs in other functions, etc). Those should be submitted as separate PRs.

@leouieda

This comment has been minimized.

Member

leouieda commented Jul 27, 2016

Once you push the tests to this branch we can take a look at the TravisCI output and I'll comment on the code.

@leouieda

This comment has been minimized.

Member

leouieda commented Jul 27, 2016

Great work on the tests by the way!

@santis19

This comment has been minimized.

Member

santis19 commented Oct 1, 2016

I've just pushed a few more commits: the tests have been added and I've renamed the radial_average function to radial_average_spectrum as @leouieda proposed. Now all checks have passed!

@leouieda

This comment has been minimized.

Member

leouieda commented Nov 19, 2016

@santis19 perfect! Sorry for the delay, I was away on vacation 🏖

I've sent you a pull request to add the changelog entry. You can merge that and add yourself to doc/contributors.rst and we're done! Let me know as soon as you do this and I'll merge this right away.

Thank you so much for taking the time to implement this! I'm working on ways of making the PR processes faster and less of a burden.

@leouieda leouieda modified the milestones: 0.5, 0.6 Nov 19, 2016

@leouieda

This comment has been minimized.

Member

leouieda commented Nov 23, 2016

@santis19 great work! I'm merging this now. Sorry it took so long. I'm rethinking the contribution process to try to make things faster.

Thanks for putting time into this!

@leouieda leouieda changed the title from Power-Density Spectra and Radial Average to Functions to caluclate the radial average power spectrum Nov 23, 2016

@leouieda leouieda merged commit 81f865a into fatiando:master Nov 23, 2016

3 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
coverage/coveralls Coverage increased (+0.1%) to 73.57%
Details

@santis19 santis19 deleted the santis19:pds_radial_average branch Nov 23, 2017

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment