nx=779
nt=10000
dt=0.001
tmax=(nt-1)*dt
t=np.linspace(0.,tmax,nt)
seismo=np.fromfile("seismo.45.drt",dtype=float32)
seismo.shape=(nx,nt)
sm=np.array(seismo,dtype=float64)
def plot_seismo(sm):
seismo_plot=sm.clip(-1.2e-7,1.2e-7)
figure(figsize=[10,8])
dx=0.02
xmax=(nx-1)*dx
imshow(seismo_plot.T,cmap="gray_r",aspect=2,extent=[0,xmax,tmax,0])
ylabel("Time (s)",fontsize="large")
xlabel("Distance (km)",fontsize="large")
epow=0.2
plot_seismo(sm)
def gain_numpy(sm,t,epow):
g=np.empty_like(sm)
table=np.empty_like(t)
for it in xrange(t.shape[0]):
table[it]=exp(epow*t[it])
for ix in xrange(sm.shape[0]):
for it in xrange(t.shape[0]):
g[ix,it]=sm[ix,it]*table[it]
return g
def gain_numpy_vec(sm,t,epow):
return sm*exp(epow*t)
g=gain_numpy(sm,t,epow)
plot_seismo(g)
g=gain_numpy_vec(sm,t,epow)
plot_seismo(g)
ncut=100
%timeit -r 3 -n 3 gain_numpy(sm[:,:ncut],t[:ncut],epow)
3 loops, best of 3: 50.3 ms per loop
ncut=1000
%timeit -r 3 -n 3 gain_numpy(sm[:,:ncut],t[:ncut],epow)
3 loops, best of 3: 512 ms per loop
%timeit -r 3 -n 3 gain_numpy(sm,t,epow)
3 loops, best of 3: 5.06 s per loop
ncut=100
%timeit -r 3 -n 3 gain_numpy_vec(sm[:,:ncut],t[:ncut],epow)
3 loops, best of 3: 112 µs per loop
ncut=1000
%timeit -r 3 -n 3 gain_numpy_vec(sm[:,:ncut],t[:ncut],epow)
3 loops, best of 3: 2.41 ms per loop
%timeit -r 3 -n 3 gain_numpy_vec(sm,t,epow)
3 loops, best of 3: 17.4 ms per loop
import numba
gain_numba=numba.autojit(gain_numpy)
g=gain_numba(sm,t,epow)
plot_seismo(g)
ncut=100
%timeit -r 3 -n 3 gain_numba(sm[:,:ncut],t[:ncut],epow)
3 loops, best of 3: 192 µs per loop
ncut=1000
%timeit -r 3 -n 3 gain_numba(sm[:,:ncut],t[:ncut],epow)
3 loops, best of 3: 2.4 ms per loop
%timeit -r 3 -n 3 gain_numba(sm,t,epow)
3 loops, best of 3: 23.1 ms per loop
%load_ext cythonmagic
%%cython
import cython
import numpy as np
cimport numpy as np
from libc.math cimport exp
ctypedef np.double_t DTYPE_t
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
def gain_cy(np.ndarray[DTYPE_t,ndim=2] sm,np.ndarray[DTYPE_t] t,double epow):
cdef np.ndarray[DTYPE_t,ndim=2] g=np.empty_like(sm)
cdef np.ndarray[DTYPE_t] table=np.empty_like(t)
cdef int nx=sm.shape[0]
cdef int nt=sm.shape[1]
cdef int ix,it
for it in xrange(nt):
table[it]=exp(epow*t[it])
for ix in xrange(nx):
for it in xrange(nt):
g[ix,it]=sm[ix,it]*table[it]
return g
g=gain_cy(sm,t,epow)
plot_seismo(g)
ncut=100
%timeit -r 3 -n 3 gain_cy(sm[:,:ncut],t[:ncut],epow)
3 loops, best of 3: 100 µs per loop
ncut=1000
%timeit -r 3 -n 3 gain_cy(sm[:,:ncut],t[:ncut],epow)
3 loops, best of 3: 1.81 ms per loop
%timeit -r 3 -n 3 gain_cy(sm,t,epow)
3 loops, best of 3: 18.1 ms per loop
%%file gain_seismo.f90
subroutine gain(nt,nx,sm,t,epow,g)
implicit none
integer,intent(in):: nt,nx
real(kind=8),intent(in):: sm(nt,nx),t(nt),epow
real(kind=8),intent(out):: g(nt,nx)
real(kind=8):: table(nt)
integer it,ix
do it=1,nt
table(it)=dexp(epow*t(it))
enddo
do ix=1,nx
do it=1,nt
g(it,ix)=sm(it,ix)*table(it)
enddo
enddo
end subroutine
subroutine gain_v(nt,nx,sm,t,epow,g)
implicit none
integer,intent(in):: nt,nx
real(kind=8),intent(in):: sm(nt,nx),t(nt),epow
real(kind=8),intent(out):: g(nt,nx)
real(kind=8):: table(nt)
integer ix
table(:)=dexp(epow*t(:))
do ix=1,nx
g(:,ix)=sm(:,ix)*table(:)
enddo
end subroutine
Overwriting gain_seismo.f90
!f2py -c -m gain_sm gain_seismo.f90 --f90exec=/opt/local/bin/gfortran-mp-4.9 > log.txt
import gain_sm
print gain_sm.__doc__
This module 'gain_sm' is auto-generated with f2py (version:2). Functions: g = gain(sm,t,epow,nt=shape(sm,0),nx=shape(sm,1)) g = gain_v(sm,t,epow,nt=shape(sm,0),nx=shape(sm,1)) .
smt=sm.T
g=gain_sm.gain(smt,t,epow)
plot_seismo(g.T)
ncut=100
%timeit -r 3 -n 3 gain_sm.gain(smt[:ncut,:],t[:ncut],epow)
3 loops, best of 3: 115 µs per loop
ncut=1000
%timeit -r 3 -n 3 gain_sm.gain(smt[:ncut,:],t[:ncut],epow)
3 loops, best of 3: 4.12 ms per loop
%timeit -r 3 -n 3 gain_sm.gain(smt,t,epow)
3 loops, best of 3: 20.6 ms per loop
ncut=100
%timeit -r 3 -n 3 gain_sm.gain_v(smt[:ncut,:],t[:ncut],epow)
3 loops, best of 3: 142 µs per loop
ncut=1000
%timeit -r 3 -n 3 gain_sm.gain_v(smt[:ncut,:],t[:ncut],epow)
3 loops, best of 3: 4.11 ms per loop
%timeit -r 3 -n 3 gain_sm.gain_v(smt,t,epow)
3 loops, best of 3: 20.7 ms per loop