nt=10001
tmax=100
t=np.linspace(0,tmax,nt)
ntr=1000
t
array([ 0.00000000e+00, 1.00000000e-02, 2.00000000e-02, ..., 9.99800000e+01, 9.99900000e+01, 1.00000000e+02])
sig=sin(t)
plot(t,sig)
[<matplotlib.lines.Line2D at 0x108bd32d0>]
sigs=array([sig for i in xrange(ntr)])
sigs.shape
(1000, 10001)
noise=np.random.rand(ntr,nt)
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')
plot(t,gather[ntr/2],'k-')
xlabel('Time (s)',fontsize='large')
ylabel('Amplitude',fontsize='large')
savefig('trace.eps')
def stack_py(gather,nt):
stacked=np.zeros(nt)
for trc in gather:
stacked+=trc
return stacked
stack=stack_py(gather,nt)
plot(t,stack)
[<matplotlib.lines.Line2D at 0x10acc94d0>]
%timeit -r 5 -n 3 stack_py(gather,nt)
3 loops, best of 5: 7.02 ms per loop
def stack_np(gather,ntr,nt):
return gather.sum(0)
stack=stack_np(gather,ntr,nt)
plot(t,stack,'k-')
xlabel('Time (s)',fontsize='large')
ylabel('Amplitude',fontsize='large')
savefig('stacked.eps')
%timeit -r 5 -n 3 stack_np(gather,ntr,nt)
3 loops, best of 5: 6.34 ms per loop
from numba import autojit
stack_nb=autojit(stack_py)
%timeit -r 5 -n 3 stack_nb(gather,ntr,nt)
3 loops, best of 5: 19.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 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
%timeit -r 5 -n 3 stack_cy(gather,ntr,nt)
3 loops, best of 5: 8.87 ms per loop
%%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
!f2py -c -m stack_f stack.f90 --f90exec=/opt/local/bin/gfortran-mp-4.9 > log.txt
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)) .
gather_t=gather.T
%timeit -r 5 -n 3 stack_f.stack_f(gather_t)
3 loops, best of 5: 5.72 ms per loop
%timeit -r 5 -n 3 stack_f.stack_fv(gather_t)
3 loops, best of 5: 28.6 ms per loop
%timeit -r 5 -n 3 stack_f.stack_fv2(gather_t)
3 loops, best of 5: 5.75 ms per loop
plot(t,stack_f.stack_f(gather_t))
[<matplotlib.lines.Line2D at 0x11a2e7a10>]