Benchmark and profile the PrismMesh class

In [1]:

    

%pylab inline
from IPython.display import Image
from fatiando.mesher import PrismMesh
from fatiando.utils import SparseList
from fatiando.vis import myv
from line_profiler import LineProfiler


    

    




Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

In [16]:

    

bounds = (0, 1000, 0, 2000, 0, 500)
shape = (20, 20, 20)
mesh = PrismMesh(bounds, shape)
nz, ny, nx = shape
#density = zeros(mesh.size)
#density[::2] = arange(1, mesh.size/2 + 1)
density = SparseList(mesh.size)
density.elements = dict(zip(range(0, mesh.size, 2), arange(1, mesh.size/2)))
mesh.addprop('density', density)
print mesh.size, len(mesh.mask)
backup = copy(density)


    

    




8000 0

In [17]:

    

mesh.mask = arange(mesh.size)[::3]
print mesh.size, len(mesh.mask)


    

    




8000 2667

In [13]:

    

myv.figure()
myv.prisms(mesh, 'density')
myv.axes(myv.outline(bounds))
myv.savefig('mesher_prismmesh.png')
myv.show()
Image(filename='mesher_prismmesh.png')


    

    




/usr/lib/python2.7/dist-packages/mayavi/preferences/preference_manager.py:24: UserWarning: Module docutils was already imported from /usr/lib/python2.7/dist-packages/docutils/__init__.pyc, but /home/leo/.virtualenvs/fatiando/lib/python2.7/site-packages is being added to sys.path
  import pkg_resources
/usr/lib/python2.7/dist-packages/mayavi/preferences/preference_manager.py:24: UserWarning: Module dap was already imported from None, but /usr/lib/python2.7/dist-packages is being added to sys.path
  import pkg_resources






    
Out[13]:

In [23]:

    

%timeit mesh.dump('tmp.msh', 'tmp.den', 'density')


    

    




10 loops, best of 3: 20.5 ms per loop

In [19]:

    

mesh.dump('tmp.msh', 'tmp.den', 'density')
all(mesh.props['density'] == backup)


    

    
Out[19]:





True

In [24]:

    

lp = LineProfiler(PrismMesh.dump)
lp.runcall(mesh.dump, 'tmp.msh', 'tmp.den', 'density')
lp.print_stats()


    

    




Timer unit: 1e-06 s

File: /home/leo/src/fatiando/fatiando/mesher.py
Function: dump at line 1044
Total time: 0.061443 s

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
  1044                                               def dump(self, meshfile, propfile, prop):
  1045                                                   r"""
  1046                                                   Dump the mesh to a file in the format required by UBC-GIF program
  1047                                                   MeshTools3D.
  1048                                           
  1049                                                   Parameters:
  1050                                           
  1051                                                   * meshfile : str or file
  1052                                                       Output file to save the mesh. Can be a file name or an open file.
  1053                                                   * propfile : str or file
  1054                                                       Output file to save the physical properties *prop*. Can be a file
  1055                                                       name or an open file.
  1056                                                   * prop : str
  1057                                                       The name of the physical property in the mesh that will be saved to
  1058                                                       *propfile*.
  1059                                           
  1060                                                   .. note:: Uses -10000000 as the dummy value for plotting topography
  1061                                           
  1062                                                   Examples:
  1063                                           
  1064                                                       >>> from StringIO import StringIO
  1065                                                       >>> meshfile = StringIO()
  1066                                                       >>> densfile = StringIO()
  1067                                                       >>> mesh = PrismMesh((0, 10, 0, 20, 0, 5), (1, 2, 2))
  1068                                                       >>> mesh.addprop('density', [1, 2, 3, 4])
  1069                                                       >>> mesh.dump(meshfile, densfile, 'density')
  1070                                                       >>> print meshfile.getvalue().strip()
  1071                                                       2 2 1
  1072                                                       0 0 0
  1073                                                       2*10
  1074                                                       2*5
  1075                                                       1*5
  1076                                                       >>> print densfile.getvalue().strip()
  1077                                                       1.0000
  1078                                                       3.0000
  1079                                                       2.0000
  1080                                                       4.0000
  1081                                           
  1082                                                   """
  1083         1            6      6.0      0.0          if prop not in self.props:
  1084                                                       raise ValueError("mesh doesn't have a '%s' property." % (prop))
  1085         1            4      4.0      0.0          isstr = False
  1086         1            4      4.0      0.0          if isinstance(meshfile, str):
  1087         1            3      3.0      0.0              isstr = True
  1088         1          143    143.0      0.2              meshfile = open(meshfile, 'w')
  1089         1            5      5.0      0.0          nz, ny, nx = self.shape
  1090         1            3      3.0      0.0          x1, x2, y1, y2, z1, z2 = self.bounds
  1091         1            4      4.0      0.0          dx, dy, dz = self.dims
  1092         1            4      4.0      0.0          meshfile.writelines([
  1093         1           17     17.0      0.0              "%d %d %d\n" % (ny, nx, nz),
  1094         1            9      9.0      0.0              "%g %g %g\n" % (y1, x1, -z1),
  1095         1            9      9.0      0.0              "%d*%g\n" % (ny, dy),
  1096         1            7      7.0      0.0              "%d*%g\n" % (nx, dx),
  1097         1           87     87.0      0.1              "%d*%g" % (nz, dz)])
  1098         1            5      5.0      0.0          if isstr:
  1099         1          156    156.0      0.3              meshfile.close()
  1100         1        29622  29622.0     48.2          values = numpy.fromiter(self.props[prop], dtype='f')
  1101                                                   # Replace the masked cells with a dummy value
  1102         1           51     51.0      0.1          values[self.mask] = -10000000
  1103         1           43     43.0      0.1          reordered = numpy.ravel(numpy.reshape(values, self.shape), order='F')
  1104         1        31261  31261.0     50.9          numpy.savetxt(propfile, reordered, fmt='%.4f')

In [8]:

    




    

In [8]:

    




    

In [8]:

    




    

In [8]: