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
imshow(vel.T)
<matplotlib.image.AxesImage at 0x108bd4190>
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)
vel
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)
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
scaled=scale_py(vel,nx,nz,factor)
scaled
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)
plot_vel(scaled,'BP_meter.eps',h,'m/s',arange(1500,6000,1000))
%timeit -r 3 -n 3 scale_py(vel,nx,nz,factor)
3 loops, best of 3: 2.66 s per loop
def scale_np(vel,factor):
return vel*factor
%timeit -r 3 -n 3 scale_np(vel,factor)
3 loops, best of 3: 6.01 ms per loop
from numba import autojit
scale_nb=autojit(scale_py)
%timeit -r 3 -n 3 scale_nb(vel,nx,nz,factor)
3 loops, best of 3: 12.7 ms per loop
%load_ext cythonmagic
%%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
%timeit -r 3 -n 3 scale_cy(vel,nx,nz,factor)
3 loops, best of 3: 7.06 ms per loop
%%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
!f2py -c -m scale_f scale.f90 --f90exec=/opt/local/bin/gfortran-mp-4.9 > log.txt
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)) .
vel_t=vel.T
%timeit -r 3 -n 3 scale_f.scale_f(vel_t,factor)
3 loops, best of 3: 7.26 ms per loop
%timeit -r 3 -n 3 scale_f.scale_fv(vel_t,factor)
3 loops, best of 3: 7.25 ms per loop