Tesseroid forward calculation in Cython

[santisoler]

Due to instabilities with Numba, lack of proper documentation, and hard-to-debug codes, @leouieda and I thought that would be better to migrate the Tesseroid forward calculations to Cython. This new Cython code have been recovered from an old commit (d2ef02d) and adapted to work with the new tesseroid.py, and the new stack algorithm (see Uieda2016).

Although Numba is a lot easier to write, I think is only better than Cython for high performance calculations where variables have simple types. I'm currently working on improvements in Tesseroid forward calculation and I need to manipulate more sophisticated variables, such as passing functions as arguments.

On the other hand, comparing the speed of the Numba code with this new Cython code, the later is faster than the former. The following codes takes around 18s to run using the Numba code, and 10s using this Cython code (in a Intel i3 2.4Ghz processor). This big difference is mainly due to the fact that the Numba code must be compiled on the first run.

import time
from fatiando import gridder
from fatiando.gravmag import tesseroid
from fatiando.mesher import Tesseroid


model = [Tesseroid(-60, -55, -30, -27, 0, -500000, props={'density': 200}),
         Tesseroid(-66, -62, -18, -12, 0, -300000, props={'density': -500})]

# Create the computation grid
area = (-80, -30, -40, 10)
shape = (100, 100)
lons, lats, heights = gridder.regular(area, shape, z=250000)

start = time.time()
fields = [
    tesseroid.potential(lons, lats, heights, model),
    tesseroid.gx(lons, lats, heights, model),
    tesseroid.gy(lons, lats, heights, model),
    tesseroid.gz(lons, lats, heights, model),
    tesseroid.gxx(lons, lats, heights, model),
    tesseroid.gxy(lons, lats, heights, model),
    tesseroid.gxz(lons, lats, heights, model),
    tesseroid.gyy(lons, lats, heights, model),
    tesseroid.gyz(lons, lats, heights, model),
    tesseroid.gzz(lons, lats, heights, model)
    ]
print "Time it took: %s" % (time.time() - start)

Comparing both codes with a mesh of Tesseroid gives more similar times, but Cython keeps winning.
Numba takes 136s and Cython 116s (tested on the same CPU).

import time
import numpy as np
from fatiando import gridder
from fatiando.gravmag import tesseroid
from fatiando.mesher import Tesseroid, TesseroidMesh

         
model = TesseroidMesh((-66, -55, -30, -12, 0, -500000), (10, 10, 2))
model.addprop('density', np.ones(model.size)*-500)

# Create the computation grid
area = (-80, -30, -40, 10)
shape = (100, 100)
lons, lats, heights = gridder.regular(area, shape, z=250000)

start = time.time()
fields = [
    tesseroid.potential(lons, lats, heights, model),
    tesseroid.gx(lons, lats, heights, model),
    tesseroid.gy(lons, lats, heights, model),
    tesseroid.gz(lons, lats, heights, model),
    tesseroid.gxx(lons, lats, heights, model),
    tesseroid.gxy(lons, lats, heights, model),
    tesseroid.gxz(lons, lats, heights, model),
    tesseroid.gyy(lons, lats, heights, model),
    tesseroid.gyz(lons, lats, heights, model),
    tesseroid.gzz(lons, lats, heights, model)
    ]
print "Time it took: %s" % (time.time() - start)

The tests will fail because pytest does not recognize when Cython raises an OverflowError as such, but it throws a Segmentation Fault. The code raises the proper OverflowError when the tesseroid stack tries to be bigger than the STACK_SIZE, so I assume that the problem is in the pytest function. If anyone knows how to overcome this, any advise would be helpful.

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, environment.yml, ci/requirements.txt.
• Documentation builds properly (run cd doc; make locally)
• Changelog entry (leave for last to avoid conflicts)
• First-time contributor? Add yourself to doc/contributors.rst (leave for last to avoid conflicts)