Gain



In [32]:

    
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)

Python, Numpy



In [33]:

    
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)



In [34]:

    
g=gain_numpy(sm,t,epow)
plot_seismo(g)



In [35]:

    
g=gain_numpy_vec(sm,t,epow)
plot_seismo(g)



In [36]:

    
ncut=100
%timeit -r 3 -n 3 gain_numpy(sm[:,:ncut],t[:ncut],epow)









    



3 loops, best of 3: 50.3 ms per loop



In [37]:

    
ncut=1000
%timeit -r 3 -n 3 gain_numpy(sm[:,:ncut],t[:ncut],epow)









    



3 loops, best of 3: 512 ms per loop



In [38]:

    
%timeit -r 3 -n 3 gain_numpy(sm,t,epow)









    



3 loops, best of 3: 5.06 s per loop



In [38]:



In [39]:

    
ncut=100
%timeit -r 3 -n 3 gain_numpy_vec(sm[:,:ncut],t[:ncut],epow)









    



3 loops, best of 3: 112 µs per loop



In [40]:

    
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



In [41]:

    
%timeit -r 3 -n 3 gain_numpy_vec(sm,t,epow)









    



3 loops, best of 3: 17.4 ms per loop

Numba



In [42]:

    
import numba

gain_numba=numba.autojit(gain_numpy)



In [43]:

    
g=gain_numba(sm,t,epow)
plot_seismo(g)



In [44]:

    
ncut=100
%timeit -r 3 -n 3 gain_numba(sm[:,:ncut],t[:ncut],epow)









    



3 loops, best of 3: 192 µs per loop



In [45]:

    
ncut=1000
%timeit -r 3 -n 3 gain_numba(sm[:,:ncut],t[:ncut],epow)









    



3 loops, best of 3: 2.4 ms per loop



In [46]:

    
%timeit -r 3 -n 3 gain_numba(sm,t,epow)









    



3 loops, best of 3: 23.1 ms per loop

Cython



In [1]:

    
%load_ext cythonmagic



In [2]:

    
%%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



In [49]:

    
g=gain_cy(sm,t,epow)
plot_seismo(g)



In [50]:

    
ncut=100
%timeit -r 3 -n 3 gain_cy(sm[:,:ncut],t[:ncut],epow)









    



3 loops, best of 3: 100 µs per loop



In [51]:

    
ncut=1000
%timeit -r 3 -n 3 gain_cy(sm[:,:ncut],t[:ncut],epow)









    



3 loops, best of 3: 1.81 ms per loop



In [52]:

    
%timeit -r 3 -n 3 gain_cy(sm,t,epow)









    



3 loops, best of 3: 18.1 ms per loop

Fortran



In [53]:

    
%%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



In [54]:

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



In [55]:

    
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))
.



In [56]:

    
smt=sm.T
g=gain_sm.gain(smt,t,epow)
plot_seismo(g.T)



In [57]:

    
ncut=100
%timeit -r 3 -n 3 gain_sm.gain(smt[:ncut,:],t[:ncut],epow)









    



3 loops, best of 3: 115 µs per loop



In [58]:

    
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



In [59]:

    
%timeit -r 3 -n 3 gain_sm.gain(smt,t,epow)









    



3 loops, best of 3: 20.6 ms per loop



In [59]:



In [60]:

    
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



In [61]:

    
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



In [62]:

    
%timeit -r 3 -n 3 gain_sm.gain_v(smt,t,epow)









    



3 loops, best of 3: 20.7 ms per loop