F2Py란?
- 수치 계산의 '우사인 볼트'라 할 수 있는 Fortran은 단점으로 다양한 파일의 입출력 같은 부가기능이 약합니다.
- 다재다능 Python은 단점으로 단순 계산이 조금 느립니다.
- 그래서 Python 뼈대에 계산 부분만 Fortran을 불러올 수 있게 도와주는 것이 F2PY입니다.
- F2PY는 Numpy에 기본으로 내재되어 있습니다.
- 이 포스팅은 기본적으로 Linux 환경에서 Python3를 이용함을 알려드립니다.
1. Intro
이번에는 Moving Average 이동평균을 계산해 보겠습니다.
가장 간단한 구현 방법은 Loop를 2번 돌리면 되겠습니다.
subroutine ma1d(A,B,m,mpd)
!!!
!!! Moving average from -mpd to +mpd
!!!
implicit none
integer :: m,i,j,mpd
real(8), dimension(m),intent(in) :: A
real(8), dimension(m),intent(out) :: B
!F2PY INTENT(HIDE) :: m
!F2PY INTENT(IN) :: mpd
B=0.d0
do i=1+mpd,m-mpd
do j=-mpd,mpd
B(i)=B(i)+A(i+j)
enddo
enddo
B=B/dble(mpd*2+1)
end subroutine ma1d
위는 Fortran90 코드이며, 다음과 같이 컴파일을 해줍니다.
$ f2py3 -c -m moving_average moving_average.f90
그러면 "moving_average.cpython-34m.so"라는 파일이 생성되며, Python에서 읽을 준비가 되었습니다.
위 함수를 불러들이고 실행할 모체 Python 프로그램은 다음과 같습니다.
import moving_average as mvavg
import numpy as np
import timeit
import sys
print("*** Module Document ***")
print(mvavg.__doc__)
print("*** Function Document ***")
print(mvavg.ma1d.__doc__)
### Prepare input array for test
nt=1000000
array=np.sin(np.arange(nt))
mpd=501; mpd_half=int((mpd-1)/2)
#array=np.arange(8); mpd_half=2
print("*** Start Loop ***")
time_record=[]
for i in range(20):
start = timeit.default_timer()
c= mvavg.ma1d(array,mpd_half)
stop = timeit.default_timer()
time_record.append(stop-start)
print("1st: {:.5f}sec 2nd: {:.5f}sec 3rd: {:.5f}sec".format(*sorted(time_record)[:3]))
print(c[1000:1003])
길이 1,000,000의 1D Vector에 501항목의 이동평균을 계산해서 역시 길이 1,000,000의 결과물을 내놓습니다.
위 프로그램을 실행시키면 다음과 같은 결과가 나옵니다.
*** Module Document ***
This module 'moving_average' is auto-generated with f2py (version:2).
Functions:
b = ma1d(a,mpd)
b = ma1dv2(a,mpd)
.
*** Function Document ***
b = ma1d(a,mpd)
Wrapper for ``ma1d``.
Parameters
----------
a : input rank-1 array('d') with bounds (m)
mpd : input int
Returns
-------
b : rank-1 array('d') with bounds (m)
*** Start Loop ***
1st: 0.91144sec 2nd: 0.91146sec 3rd: 0.91152sec
[-0.00253436 -0.00281975 -0.00051267]
대략 0.91초 정도가 걸립니다.
(참고로, 길이 1,000,000인 Vector에 8Byte-float가 저장되었으므로 데이타 크기는 8MB입니다.)
2. If using Python only
Python만으로 같은 조건을 계산하여 비교해보겠습니다.
전편에서도 얘기했지만, Python - Numpy는 행렬 통째로 계산하는게 빠릅니다.
그래서 다음과 같은 함수를 만들었습니다.
import numpy as np
def ma1d0(A,mpd):
"""
Input: 1D array, and half period of moving average window
Output: Same dimension as input array
To do: Moving average with given window size
"""
m=A.shape[0]
B=np.zeros_like(A)
for i in range(-mpd,mpd+1,1):
B[mpd:m-mpd]+=A[mpd+i:m-mpd+i]
B=B/float(mpd*2+1)
return B
위 함수를 모체에서 불러와 같은 조건으로 실행시키면 다음과 같은 결과가 나옵니다.
*** Start Loop ***
1st: 1.29482sec 2nd: 1.29516sec 3rd: 1.29565sec
[-0.00253436 -0.00281975 -0.00051267]
약 1.29초로써 꽤 선전한 모습입니다.
3. Improving for Better Result
이동 평균 Moving Average에 대해 검색해보니 그 유명한 StackOverflow에서 다음과 같은 글을 찾았습니다.
https://stackoverflow.com/questions/14313510/how-to-calculate-moving-average-using-numpy
요약하자면, Numpy에 내재된 함수 중 cumulative sum이라는 함수를 이용하면 계산이 더 빠르다는 내용입니다.
Numpy.cumsum()
예를 들자면 1D array A에 대하여,
Sum(A[i-n:i+n+1])=Sum(A[:i+n+1])-Sum(A[:i-n]) 임을 이용한다는 것입니다.
그래서 Python 함수를 다음과 같이 고쳤습니다.
def ma1d1(A, mpd) :
mpd2=mpd*2+1
B = np.cumsum(A)
tmp=B[mpd*2]
B[mpd+1:-mpd]=B[mpd2:]-B[:-mpd2]
B[mpd]=tmp
return B/float(mpd2)
그러면 다음과 같이 시간이 줄어듭니다.
*** Start Loop ***
1st: 0.01452sec 2nd: 0.01452sec 3rd: 0.01452sec
[-0.00253436 -0.00281975 -0.00051267]
약 0.0145초입니다.
제일 위의 Fortran 결과보다 훨씬 빨라졌죠?
그럼 마지막으로 같은 알고리즘을 Fortran으로도 구현해볼게요.
subroutine ma1dv2(A,B,m,mpd)
!!!
!!! Moving average from -mpd to +mpd
!!!
implicit none
integer :: m,i,j,mpd,mpd2
real(8), dimension(m),intent(in) :: A
real(8), dimension(m),intent(out) :: B
real(8) :: tmp
!F2PY INTENT(HIDE) :: m
!F2PY INTENT(IN) :: mpd
!!! Cumulative Sum
B=0.d0; B(1)=A(1)
do i=2,m
B(i)=B(i-1)+A(i)
enddo
!!! Moving Average
mpd2=mpd*2+1
tmp=B(mpd2)
B(mpd+2:m-mpd)=B(mpd2+1:)-B(m-mpd2)
B(mpd+1)=tmp
B=B/dble(mpd2)
end subroutine ma1dv2
그러면 결과는
*** Start Loop ***
1st: 0.01274sec 2nd: 0.01275sec 3rd: 0.01275sec
[-0.00253436 -0.00281975 -0.00051267]
약 0.0127초 정도 되겠습니다.
4. Summary and Conclusion
길이 1,000,000인 1D-Array (or Vector)를 가지고 501항목의 이동 평균을 계산하면,
| Program | Algorithm | Running Time |
|---|---|---|
| Python | By Array | 1.29 sec |
| Python | Using Cumsum | 0.0145 sec |
| Fortran | By Element | 0.91 s |
| Fortran | Using Cumsum | 0.0127 sec |
Fortran은 대체로 Python-Numpy보다 빠릅니다. 하지만 그보다 더 중요한 것은 최적화된 알고리즘입니다.
"""
제 개인적 목표는
- Object-Oriented Programming과 친해지기
- Github와 친해지기 입니다.
이 목표에 닿기 위해 일단 제가 나름 좀 아는 Python, 그 중에서도 Numpy와 Matplotlib로부터 시작하려 합니다.
"""
Matplotlib List
[Matplotlib] 00. Intro + 01. Page Setup
[Matplotlib] 02. Axes Setup: Subplots
[Matplotlib] 03. Axes Setup: Text, Label, and Annotation
[Matplotlib] 04. Axes Setup: Ticks and Tick Labels
[Matplotlib] 05. Plot Accessories: Grid and Supporting Lines
[Matplotlib] 06. Plot Accessories: Legend
[Matplotlib] 07. Plot Main: Plot
[Matplotlib] 08. Plot Main: Imshow
[Matplotlib] 09. Plot Accessary: Color Map (part1)
[Matplotlib] 10. Plot Accessary: Color Map (part2) + Color Bar
F2PY List
[F2PY] 01. Basic Example: Simple Gaussian 2D Filter
[F2PY] 02. Basic Example: Moving Average