In [1]:

    

nt=10001
tmax=100
t=np.linspace(0,tmax,nt)
ntr=1000
t


    

    
Out[1]:





array([  0.00000000e+00,   1.00000000e-02,   2.00000000e-02, ...,
         9.99800000e+01,   9.99900000e+01,   1.00000000e+02])

In [2]:

    

sig=sin(t)
plot(t,sig)


    

    
Out[2]:





[<matplotlib.lines.Line2D at 0x108bd32d0>]






    

In [3]:

    

sigs=array([sig for i in xrange(ntr)])
sigs.shape


    

    
Out[3]:





(1000, 10001)

In [4]:

    

noise=np.random.rand(ntr,nt)


    

In [5]:

    

gather=sigs+noise
imshow(gather,aspect=0.1,extent=[0,tmax,ntr,0],cmap=cm.gray_r)
xlabel('Time (s)',fontsize='large')
ylabel('Trace number',fontsize='large')
savefig('gather.eps')


    

In [6]:

    

plot(t,gather[ntr/2],'k-')
xlabel('Time (s)',fontsize='large')
ylabel('Amplitude',fontsize='large')
savefig('trace.eps')


    

In [9]:

    

def stack_py(gather,nt):
    stacked=np.zeros(nt)
    for trc in gather:
        stacked+=trc
    return stacked


    

In [8]:

    

stack=stack_py(gather,nt)
plot(t,stack)


    

    
Out[8]:





[<matplotlib.lines.Line2D at 0x10acc94d0>]






    

In [11]:

    

%timeit -r 5 -n 3 stack_py(gather,nt)


    

    




3 loops, best of 5: 7.02 ms per loop

In [10]:

    

def stack_np(gather,ntr,nt):
    return gather.sum(0)


    

In [11]:

    

stack=stack_np(gather,ntr,nt)
plot(t,stack,'k-')
xlabel('Time (s)',fontsize='large')
ylabel('Amplitude',fontsize='large')
savefig('stacked.eps')


    

In [12]:

    

%timeit -r 5 -n 3 stack_np(gather,ntr,nt)


    

    




3 loops, best of 5: 6.34 ms per loop

In [13]:

    

from numba import autojit
stack_nb=autojit(stack_py)


    

In [14]:

    

%timeit -r 5 -n 3 stack_nb(gather,ntr,nt)


    

    




3 loops, best of 5: 19.7 ms per loop

In [1]:

    

%load_ext cythonmagic


    

In [2]:

    

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

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
def stack_cy(np.ndarray[double,ndim=2] gather,int ntr,int nt):
    cdef np.ndarray[double] stacked=np.zeros(nt)
    cdef int itr,it
    for itr in xrange(ntr):
        for it in xrange(nt):
            stacked[it]+=gather[itr,it]
    return stacked


    

In [17]:

    

%timeit -r 5 -n 3 stack_cy(gather,ntr,nt)


    

    




3 loops, best of 5: 8.87 ms per loop

In [18]:

    

%%file stack.f90

subroutine stack_f(stacked,gather,ntr,nt)
implicit none
integer,intent(in):: ntr,nt
real(kind=8),intent(in):: gather(nt,ntr)
real(kind=8),intent(out):: stacked(nt)
integer itr,it
do it=1,nt
    stacked(it)=0.d0
enddo
do itr=1,ntr
do it=1,nt
    stacked(it)=stacked(it)+gather(it,itr)
enddo
enddo
end subroutine

subroutine stack_fv(stacked,gather,ntr,nt)
implicit none
integer,intent(in):: ntr,nt
real(kind=8),intent(in):: gather(nt,ntr)
real(kind=8),intent(out):: stacked(nt)
stacked=sum(gather,2)
end subroutine

subroutine stack_fv2(stacked,gather,ntr,nt)
implicit none
integer,intent(in):: ntr,nt
real(kind=8),intent(in):: gather(nt,ntr)
real(kind=8),intent(out):: stacked(nt)
integer itr
stacked=0.d0
do itr=1,ntr
    stacked=stacked+gather(:,itr)
enddo
end subroutine


    

    




Overwriting stack.f90

In [19]:

    

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


    

In [20]:

    

import stack_f
print stack_f.__doc__


    

    




This module 'stack_f' is auto-generated with f2py (version:2).
Functions:
  stacked = stack_f(gather,ntr=shape(gather,1),nt=shape(gather,0))
  stacked = stack_fv(gather,ntr=shape(gather,1),nt=shape(gather,0))
  stacked = stack_fv2(gather,ntr=shape(gather,1),nt=shape(gather,0))
.

In [21]:

    

gather_t=gather.T


    

In [22]:

    

%timeit -r 5 -n 3 stack_f.stack_f(gather_t)


    

    




3 loops, best of 5: 5.72 ms per loop

In [23]:

    

%timeit -r 5 -n 3 stack_f.stack_fv(gather_t)


    

    




3 loops, best of 5: 28.6 ms per loop

In [24]:

    

%timeit -r 5 -n 3 stack_f.stack_fv2(gather_t)


    

    




3 loops, best of 5: 5.75 ms per loop

In [25]:

    

plot(t,stack_f.stack_f(gather_t))


    

    
Out[25]:





[<matplotlib.lines.Line2D at 0x11a2e7a10>]






    

In [25]: