In [1]:
%pylab inline
from line_profiler import LineProfiler
from fatiando import gravmag, gridder, utils, mesher
from fatiando.vis import mpl
In [2]:
# Make synthetic data
inc, dec = -60, 23
props = {'magnetization':10}
model = [mesher.Prism(-500, 500, -1000, 1000, 1000, 4000, props)]
shape = (50, 50)
x, y, z = gridder.regular([-5000, 5000, -5000, 5000], shape, z=-150)
tf = utils.contaminate(gravmag.prism.tf(x, y, z, model, inc, dec), 5)
# Estimate the magnetization intensity
data = [gravmag.eqlayer.TotalField(x, y, z, tf, inc, dec)]
# Setup the layer
grid = mesher.PointGrid([-5000, 5000, -5000, 5000], 200, (100, 100))
In [3]:
damping, smoothness = 0*10.**-15, 10.**-3
windows, degree = (10, 10), 1
intensity, matrices = gravmag.eqlayer.pel(data, grid, windows, degree, damping, smoothness)
grid.addprop('magnetization', intensity)
predicted = gravmag.sphere.tf(x, y, z, grid, inc, dec)
residuals = tf - predicted
print "Residuals:"
print "mean:", residuals.mean()
print "stddev:", residuals.std()
In [4]:
# Plot the layer and the fit
mpl.figure(figsize=(14,3))
mpl.subplot(1, 3, 1)
mpl.axis('scaled')
mpl.title('Layer (A/m)')
mpl.pcolor(grid.y, grid.x, grid.props['magnetization'], grid.shape)
mpl.colorbar()
mpl.m2km()
mpl.subplot(1, 3, 2)
mpl.axis('scaled')
mpl.title('Fit (nT)')
levels = mpl.contour(y, x, tf, shape, 15, color='r')
mpl.contour(y, x, predicted, shape, levels, color='k')
mpl.m2km()
mpl.subplot(1, 3, 3)
mpl.title('Residuals (nT)')
mpl.hist(residuals, bins=10)
Out[4]:
Run a classic layer to see
In [5]:
classic_damping = 0.02
intensity, predicted = gravmag.eqlayer.classic(data, grid, classic_damping)
grid.addprop('magnetization', intensity)
residuals = tf - predicted[0]
print "Residuals:"
print "mean:", residuals.mean()
print "stddev:", residuals.std()
In [6]:
# Plot the layer and the fit
mpl.figure(figsize=(14,3))
mpl.subplot(1, 3, 1)
mpl.axis('scaled')
mpl.title('Layer (A/m)')
mpl.pcolor(grid.y, grid.x, grid.props['magnetization'], grid.shape)
mpl.colorbar()
mpl.m2km()
mpl.subplot(1, 3, 2)
mpl.axis('scaled')
mpl.title('Fit (nT)')
levels = mpl.contour(y, x, tf, shape, 15, color='r')
mpl.contour(y, x, predicted[0], shape, levels, color='k')
mpl.m2km()
mpl.subplot(1, 3, 3)
mpl.title('Residuals (nT)')
mpl.hist(residuals, bins=10)
Out[6]:
In [7]:
%timeit gravmag.eqlayer.classic(data, grid, classic_damping)
In [8]:
%timeit gravmag.eqlayer.pel(data, grid, windows, degree, damping, smoothness)
In [9]:
%prun -T gravmag_eqlayer.classic.profile gravmag.eqlayer.classic(data, grid, classic_damping)
!cat gravmag_eqlayer.classic.profile
In [10]:
lp = LineProfiler(gravmag.eqlayer.classic)
lp.runcall(gravmag.eqlayer.classic, data, grid, classic_damping)
lp.print_stats()
In [11]:
%prun -T gravmag_eqlayer.pel.profile gravmag.eqlayer.pel(data, grid, windows, degree, damping, smoothness)
!cat gravmag_eqlayer.pel.profile
In [12]:
lp = LineProfiler(gravmag.eqlayer.pel)
lp.runcall(gravmag.eqlayer.pel, data, grid, windows, degree, damping, smoothness)
lp.print_stats()
In [13]:
lp = LineProfiler(gravmag.eqlayer._pel_matrices)
lp.runcall(gravmag.eqlayer.pel, data, grid, windows, degree, damping, smoothness)
lp.print_stats()
In [14]:
lp = LineProfiler(gravmag.eqlayer._gkmatrix)
lp.runcall(gravmag.eqlayer.pel, data, grid, windows, degree, damping, smoothness)
lp.print_stats()
In [15]:
lp = LineProfiler(gravmag.eqlayer.TotalField.sensitivity)
lp.runcall(gravmag.eqlayer.pel, data, grid, windows, degree, damping, smoothness)
lp.print_stats()
Looking for a better way to build the sensitivity than using transpose:
In [16]:
%timeit transpose([arange(100) for i in xrange(100)])
In [17]:
def matrix():
a = zeros((100,100))
for i in xrange(100):
a[:,i] = arange(100)
return a
%timeit matrix()
a = matrix()
print a