In [1]:

    

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)
savefig("seismo.eps")


    

In [2]:

    

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 [3]:

    

g=gain_numpy(sm,t,epow)
plot_seismo(g)
savefig("seismo_gain.eps")


    

In [4]:

    

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


    

In [5]:

    

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


    

    




3 loops, best of 3: 5.22 s per loop

In [6]:

    

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


    

    




3 loops, best of 3: 17.7 ms per loop

In [10]:

    

import numba

gain_numba=numba.autojit(gain_numpy)


    

In [11]:

    

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


    

In [12]:

    

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


    

    




3 loops, best of 3: 24.8 ms per loop

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 [16]:

    

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


    

In [17]:

    

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


    

    




3 loops, best of 3: 18.4 ms per loop

In [20]:

    

%%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 [21]:

    

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


    

In [22]:

    

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 [23]:

    

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


    

In [24]:

    

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


    

    




3 loops, best of 3: 21.6 ms per loop

In [25]:

    

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


    

    




3 loops, best of 3: 23.4 ms per loop

In [25]: