In [1]:

    

h=0.0125
nx=5395
nz=956
vel=np.fromfile("bp0.0125km.drt",dtype="float32")
vel.shape=(nx,nz)
factor=np.float32(1000.0)
print vel.max(), vel.min()


    

    




4.79 1.429

In [2]:

    

imshow(vel.T)


    

    
Out[2]:





<matplotlib.image.AxesImage at 0x108bd4190>






    

In [3]:

    

def plot_vel(vel,img,h,unit='km/s',tick=arange(1.5,6,1)):
    xmax=(vel.shape[0]-1)*h
    zmax=(vel.shape[1]-1)*h
    fs=14
    figure(figsize=(18,4))
    ax=imshow(vel.transpose(),extent=(0,xmax,zmax,0),cmap=cm.gray_r)
    tick_params(labelsize=fs)
    xlabel('Distance (km)',fontsize=fs)
    ylabel('Depth (km)',fontsize=fs)

    ax.axes.set_yticks(arange(0,zmax+1,3)) 
    cb=colorbar(shrink=0.68,pad=0.01,aspect=10,ticks=tick)
    clim([vel.min(),vel.max()])
    cb.set_label(unit,fontsize=fs)
    ct=getp(cb.ax,'ymajorticklabels')
    setp(ct,fontsize=fs)
    
    savefig(img)

plot_vel(vel,'BP.eps',h)


    

In [4]:

    

vel


    

    
Out[4]:





array([[ 1.4860003 ,  1.4860003 ,  1.4860003 , ...,  4.57103157,
         4.57178926,  4.57254696],
       [ 1.4860003 ,  1.4860003 ,  1.4860003 , ...,  4.57103157,
         4.57178926,  4.57254696],
       [ 1.4860003 ,  1.4860003 ,  1.4860003 , ...,  4.57103157,
         4.57178926,  4.57254696],
       ..., 
       [ 1.4860003 ,  1.4860003 ,  1.4860003 , ...,  4.57102346,
         4.57178164,  4.57253933],
       [ 1.4860003 ,  1.4860003 ,  1.4860003 , ...,  4.57098055,
         4.57174253,  4.57250023],
       [ 1.4860003 ,  1.4860003 ,  1.4860003 , ...,  4.57094145,
         4.57170343,  4.57246113]], dtype=float32)

In [5]:

    

def scale_py(vel,nx,nz,factor):
    scaled=empty_like(vel)
    for ix in xrange(nx):
        for iz in xrange(nz):
            scaled[ix,iz]=vel[ix,iz]*factor
    return scaled


    

In [6]:

    

scaled=scale_py(vel,nx,nz,factor)


    

In [7]:

    

scaled


    

    
Out[7]:





array([[ 1486.00024414,  1486.00024414,  1486.00024414, ...,
         4571.03173828,  4571.7890625 ,  4572.546875  ],
       [ 1486.00024414,  1486.00024414,  1486.00024414, ...,
         4571.03173828,  4571.7890625 ,  4572.546875  ],
       [ 1486.00024414,  1486.00024414,  1486.00024414, ...,
         4571.03173828,  4571.7890625 ,  4572.546875  ],
       ..., 
       [ 1486.00024414,  1486.00024414,  1486.00024414, ...,
         4571.0234375 ,  4571.78173828,  4572.53955078],
       [ 1486.00024414,  1486.00024414,  1486.00024414, ...,
         4570.98046875,  4571.74267578,  4572.5       ],
       [ 1486.00024414,  1486.00024414,  1486.00024414, ...,
         4570.94140625,  4571.70361328,  4572.4609375 ]], dtype=float32)

In [8]:

    

plot_vel(scaled,'BP_meter.eps',h,'m/s',arange(1500,6000,1000))


    

In [9]:

    

%timeit -r 3 -n 3 scale_py(vel,nx,nz,factor)


    

    




3 loops, best of 3: 2.66 s per loop

In [10]:

    

def scale_np(vel,factor):
    return vel*factor


    

In [11]:

    

%timeit -r 3 -n 3 scale_np(vel,factor)


    

    




3 loops, best of 3: 6.01 ms per loop

In [14]:

    

from numba import autojit
scale_nb=autojit(scale_py)


    

In [15]:

    

%timeit -r 3 -n 3 scale_nb(vel,nx,nz,factor)


    

    




3 loops, best of 3: 12.7 ms per loop

In [2]:

    

%load_ext cythonmagic


    

In [3]:

    

%%cython
import cython
import numpy as np
cimport numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
def scale_cy(np.ndarray[float,ndim=2] vel,int nx,int nz,float factor):
    cdef np.ndarray[float,ndim=2] scaled=np.empty_like(vel)
    cdef int ix,iz
    for ix in xrange(nx):
        for iz in xrange(nz):
            scaled[ix,iz]=vel[ix,iz]*factor
    return scaled


    

In [20]:

    

%timeit -r 3 -n 3 scale_cy(vel,nx,nz,factor)


    

    




3 loops, best of 3: 7.06 ms per loop

In [21]:

    

%%file scale.f90

subroutine scale_f(scaled,vel,nx,nz,factor)
implicit none
integer,intent(in):: nx,nz
real,intent(in):: vel(nz,nx),factor
real,intent(out):: scaled(nz,nx)
integer iz,ix
do ix=1,nx
do iz=1,nz
    scaled(iz,ix)=vel(iz,ix)*factor
enddo
enddo
end subroutine

subroutine scale_fv(scaled,vel,nx,nz,factor)
implicit none
integer,intent(in):: nx,nz
real,intent(in):: vel(nz,nx),factor
real,intent(out):: scaled(nz,nx)
scaled=vel*factor
end subroutine


    

    




Overwriting scale.f90

In [22]:

    

!f2py -c -m scale_f scale.f90 --f90exec=/opt/local/bin/gfortran-mp-4.9 > log.txt


    

In [23]:

    

import scale_f
print scale_f.__doc__


    

    




This module 'scale_f' is auto-generated with f2py (version:2).
Functions:
  scaled = scale_f(vel,factor,nx=shape(vel,1),nz=shape(vel,0))
  scaled = scale_fv(vel,factor,nx=shape(vel,1),nz=shape(vel,0))
.

In [24]:

    

vel_t=vel.T


    

In [25]:

    

%timeit -r 3 -n 3 scale_f.scale_f(vel_t,factor)


    

    




3 loops, best of 3: 7.26 ms per loop

In [26]:

    

%timeit -r 3 -n 3 scale_f.scale_fv(vel_t,factor)


    

    




3 loops, best of 3: 7.25 ms per loop

In [26]: