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

Stabilize derivatives in fatiando.gravmag #196

Merged
merged 15 commits into from May 26, 2015

Conversation

Projects
None yet
3 participants
@leouieda
Member

leouieda commented Apr 22, 2015

As pointed out by @mycarta in #194, the FFT-based derivative calculations in fatiando.gravmag are unstable (see the figure below for a comparison with the MATLAB derivatives).

deriv

This PR tries to stabilize the derivatives by padding the input data and using different algorithms.

Fixes #194 and #167

Methods

The input data grids for FFT derivatives are padded using numpy.pad with the values in the borders of the grid. The padding is done to the nearest power of two. This solves the instabilities in the z derivative (specially in the edges). The horizontal derivatives are a bit better but still suffer from the "striping" seen above.

Horizontal derivatives can be stabilized by using a central differences scheme instead of the FFT. This is what matlab uses to get the stable derivatives seen in the figure. The derivx and derivy functions will have options to use the finite differences instead of FFT.

Checklist:

  • Make tests for new code
  • Create/update docstrings
  • Include relevant equations and citations in docstrings
  • Code follows PEP8 style conventions
  • Code and docs have been spellchecked
  • Include new dependencies in docs, requirements.txt, README, and .travis.yml
  • Documentation builds properly
  • All tests pass
  • Check cookbook recipes
  • Can be merged
  • Changelog entry (leave for last)
  • Firt-time contributor? Add yourself to doc/contributors.rst (leave for last)
  • Update I/O functions to read data with shape = (nx, ny)
  • Warn on the mailing list about the changes to gridder.regular
@leouieda

This comment has been minimized.

Member

leouieda commented Apr 24, 2015

I found an error in the way fatiando.gridder.regular works. In Fatiando, we assume that x is North-South and y is East-West. So if data is on a matrix, it would make sense that lines are x and columns are y. In the gridder.regular function we were assuming the opposite. The functions asked for a shape = (ny, nx). The correct form should be shape = (nx, ny).

This might impact some functions but it will be restricted mostly to the FFT derivatives so I'll make the fix in this branch.

leouieda added some commits Apr 17, 2015

Pad the data before doing FFT derivatives
Padding helps reduce edge effects of the derivatives. Use the edge
values to pad the array to the nearest power of two (almost, actually).
Frequencies for x and y were switched
Was using ny, nx = shape when in fact x is north-south so
shape = nx, ny
Started FD derivative. Test fail bc bug in gridder
gridder.regular is generating grids the wrong way. The shape should be
(nx, ny), not (ny, nx). Thus, the derivatives are in the wrong direction
and the tests fail.
Use finite-diff derivatives in Euler tests
The Euler deconvolution tests were using FFT derivatives and failing
because of them. I exchanged for derivatives approximated by finite
differences in the forward modeling. This was the Euler tests don't
depend of the performance of the FFT derivs.
Using nx, ny for shape in all gridder
This is how it should be. It breaks the tests for gravmag.sphere which
were hardcoded numbers as doctests. Now the numbers are out of order so
the test fail. Need a better test for this.
Remove doctests from gravmag.sphere
These tests weren't helping much. They compared the output to previous
output. There are better tests comparing numpy and Cython
implementations in the unit tests.
Include method arg for z deriv. Docstrings
Filled the missing arguments in the docstrings and make explicit using
fft derivs in the tests.

@leouieda leouieda force-pushed the fftderiv branch from 435b464 to 2cf1dc9 Apr 27, 2015

@leouieda leouieda modified the milestones: 0.4, 0.5 Apr 30, 2015

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 30, 2015

@mycarta here are the update results for this branch:

Old results

New results

index

Generated using this notebook.

Notice that the x and y derivatives are exchanged. There is a lot of confusion regarding this because we use x -> North-South in Fatiando and that is not always the case with other software. It's good to check that the results have the proper orientation that you'd expect.

Use finite diff derivatives by default
For the horizontal derivatives, the finite difference solutions are much
more stable than the FFT ones. They underestimate the derivative a bit
where the gradient is high but otherwise behave in a more consistent
manner and are safer to use.
@mycarta

This comment has been minimized.

mycarta commented Apr 30, 2015

Hi Leo

Thanks for the follow-up on derivatives.

Could you please take the coordinates off from the images, and the notebook out altogether or use
mock coordinates in it?

Those coordinates are from my unpublished thesis and not supposed to be out
there, they were for your eyes only.

My apologies I should have said it clearly.

In my tutorial I'd be using mock coordinates.

Thanks

Matteo

On Thu, Apr 30, 2015 at 9:30 AM, Leonardo Uieda notifications@github.com
wrote:

@mycarta https://github.com/mycarta here are the update results for
this branch:
Old results

https://cloud.githubusercontent.com/assets/290082/7272735/9a72e110-e8c3-11e4-907f-5876ded6bd8f.png
New results

[image: index]
https://cloud.githubusercontent.com/assets/290082/7416251/6ed571b8-ef34-11e4-801a-943adc275ea6.png

Generated using this notebook
http://nbviewer.ipython.org/gist/leouieda/ac346557bd28f77861b8.

Notice that the x and y derivatives are exchanged. There is a lot of
confusion regarding this because we use x -> North-South in Fatiando and
that is not always the case with other software. It's good to check that
the results have the proper orientation that you'd expect.


Reply to this email directly or view it on GitHub
#196 (comment).

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 30, 2015

I'm so sorry about that! I should have noticed that your plots have no coordinates. It's fixed now. Also removed from the notebook (which doesn't include any data).

@mycarta

This comment has been minimized.

mycarta commented Apr 30, 2015

If I'd mentioned it to begin with.....
All's well that ends well
:)
Thanks Leo

On Thu, Apr 30, 2015 at 10:54 AM, Leonardo Uieda notifications@github.com
wrote:

I'm so sorry about that! I should have noticed that your plots have no
coordinates. It's fixed now. Also removed from the notebook (which doesn't
include any data).


Reply to this email directly or view it on GitHub
#196 (comment).

leouieda added a commit that referenced this pull request May 26, 2015

Merge pull request #196 from fatiando/fftderiv
Stabilize derivatives in fatiando.gravmag

@leouieda leouieda merged commit 274007e into master May 26, 2015

2 checks passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
continuous-integration/travis-ci/push The Travis CI build passed
Details
@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented May 27, 2015

@leouieda awesome pictures of smootheness

mtb-za added a commit to mtb-za/fatiando that referenced this pull request May 27, 2015

Merge branch 'master' into mstde
Need to pull in some changes from upstream, particularly
fatiando#196 which has
made changes to `fatiando.gridder`, which this needs.

@leouieda leouieda deleted the fftderiv branch May 27, 2015

mtb-za added a commit to mtb-za/fatiando that referenced this pull request May 28, 2015

mtb-za added a commit to mtb-za/fatiando that referenced this pull request May 28, 2015

@leouieda leouieda referenced this pull request Jul 6, 2015

Merged

Fixed bug in PointGrid mesh ordering #209

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