[F2PY] 01. Basic Example: Simple Gaussian 2D Filter
F2Py란?
- 수치 계산의 '우사인 볼트'라 할 수 있는 Fortran은 단점으로 다양한 파일의 입출력 같은 부가기능이 약합니다.
- 다재다능 Python은 단점으로 단순 계산이 조금 느립니다.
- 그래서 Python 뼈대에 계산 부분만 Fortran을 불러올 수 있게 도와주는 것이 F2PY입니다.
- F2PY는 Numpy에 기본으로 내재되어 있습니다.
- 이 포스팅은 기본적으로 Linux 환경에서 Python3를 이용함을 알려드립니다.
1. Intro
오늘 테스트해볼 알고리즘은 "Simple Gaussian 2D Filter" 입니다. 말이 좀 어려운데, 실은 아래 그림과 같이 간단합니다.
이 중 왼쪽 3x3 필터(혹은 스무딩 Smoothing)를 구현해보겠습니다.
SUBROUTINE DFILTER2D(A,B,M,N)
C
DOUBLE PRECISION A(M,N)
DOUBLE PRECISION B(M,N)
INTEGER N, M
CF2PY INTENT(OUT) :: B
CF2PY INTENT(HIDE) :: N
CF2PY INTENT(HIDE) :: M
DO 20 I = 2,M-1
DO 40 J=2,N-1
B(I,J) = A(I,J) +
$ (A(I-1,J)+A(I+1,J) +
$ A(I,J-1)+A(I,J+1) )*0.5D0 +
$ (A(I-1,J-1) + A(I-1,J+1) +
$ A(I+1,J-1) + A(I+1,J+1))*0.25D0
40 CONTINUE
20 CONTINUE
B=B/4.D0
END
출처: https://docs.scipy.org/doc/numpy-1.13.0/user/c-info.python-as-glue.html#a-filtering-example
위는 Fortran77 코드이며, 매우 직관적이며 행렬의 한 원소를 기준으로 사방 팔방에 웨이트 Weight을 곱해 더하면 됩니다. 제일 앞칸의 "C"는 코멘트를 뜻하며, F2PY를 위하여 특별히 3줄이 추가되었습니다.
CF2PY INTENT(OUT) :: B
CF2PY INTENT(HIDE) :: N
CF2PY INTENT(HIDE) :: M
위 3줄은 변수 각각의 의도를 정확히 전달하기 위해 추가되었습니다. 사실 더 정확히 하자면 CF2PY INTENT(IN) :: A도 있어야 하는데, 여기서는 없어도 상관 없어서 생략되었습니다.
위 Fortran 프로그램을 Command Line에서 다음과 같이 컴파일 Compile 해줍니다.
$ f2py3 -c -m filter2d filter2d.f
(Python2를 이용할 때는 f2py라고 하면 됩니다.)
그러면 "filter2d.cpython-34m.so"라는 파일이 생성되며, Python에서 읽을 준비가 되었습니다.
위 함수를 읽을 모체 Python은 다음과 같이 준비했습니다.
import filter2d as f2d
import numpy as np
import timeit
print("*** Module Document ***")
print(f2d.__doc__)
print("*** Function Document ***")
print(f2d.dfilter2d.__doc__)
### Prepare input array for test
nx=2000; ny=1000
array=np.sin(np.arange(nx*ny)).reshape([ny,nx])
print("*** Start Loop ***")
time_record=[]
for i in range(10):
start = timeit.default_timer()
c= f2d.dfilter2d(array)
stop = timeit.default_timer()
time_record.append(stop-start)
print("1:{:.5f}sec 2:{:.5f}sec 3:{:.5f}sec".format(*sorted(time_record)[:3]))
print(c[200,200]) # In order to check the final result
계산 시간을 재기 위해 "timeit" 모듈이 이용되었음을 알려드립니다.
위 프로그램을 실행시키면 다음과 같은 결과가 나옵니다.
*** Module Document ***
This module 'filter2d' is auto-generated with f2py (version:2).
Functions:
b = dfilter2d(a)
.
*** Function Document ***
b = dfilter2d(a)
Wrapper for ``dfilter2d``.
Parameters
----------
a : input rank-2 array('d') with bounds (m,n)
Returns
-------
b : rank-2 array('d') with bounds (m,n)
*** Start Loop ***
1:0.02944sec 2:0.02946sec 3:0.02955sec
-0.227456830239
현재 제 컴퓨터 환경에서 [1000,2000] 행렬을 스무딩하는데 대략 0.03초가 안걸리네요.
2. If using Python only
Python으로만 위 알고리즘을 구현하면 시간이 어느정도 걸릴까요?
그래서 다음과 같은 함수를 만들었습니다.
def dfilter2d(A):
"""
Input: 2D array
Output: Same dimension as input array
To do: Perform simple 2D Gaussian filter
"""
dims=A.shape
if len(dims)!=2:
sys.exit("Error: Input array is not 2D")
n=dims[0]; m=dims[1] # A[n,m]
B=np.copy(A)
for j in range(1,n-1):
for i in range(1,m-1):
B[j,i]+=(A[j-1,i]+A[j+1,i]+A[j,i-1]+A[j,i+1])*0.5
B[j,i]+=(A[j-1,i-1]+A[j-1,i+1]+A[j+1,i-1]+A[j+1,i+1])*0.25
B/=4.
return B
2개의 Loop를 이용하여 위 Fortran 프로그램을 그대로 따라 했습니다. 그 결과는,
1:9.23733sec 2:9.24945sec 3:9.32466sec
-0.227456830239
대략 9.2초입니다.
그런데 위와 같은 구현은 사실 Python에게 불리하죠. 기본적으로 Numpy는 Matlab을 모태로 개발되었고, Matlab의 특징은 개개의 원소단위 계산이 아니라 행렬 통째로 접근하느데 특화되었으니까요. 그래서 다음과 같이 바꿔보았습니다.
def dfilter2d_v2(A):
"""
Input: 2D array
Output: Same dimension as input array
To do: Perform simple 2D Gaussian filter
"""
dims=A.shape
if len(dims)!=2:
sys.exit("Error: Input array is not 2D")
n=dims[0]; m=dims[1] # A[n,m]
B=np.copy(A)
B[1:-1,1:-1]+=(A[0:-2,1:-1]+A[2:,1:-1]+A[1:-1,0:-2]+A[1:-1,2:])*0.5
B[1:-1,1:-1]+=(A[0:-2,0:-2]+A[0:-2,2:]+A[2:,0:-2]+A[2:,2:])*0.25
B/=4.
return B
Loop를 없앤 결과는
1:0.09402sec 2:0.10426sec 3:0.10437sec
-0.227456830239
훨씬 단축되어 대략 0.1초 정도 되겠습니다.
3. Improving for Better Result
위 결과만으로도 Fortran은 약 3배 정도 빠릅니다. 그런데 아직 Fortran의 최고 성능은 나오지 않았습니다. 왜냐하면 위 Fortran 예제에서 Loop의 순서가 올바르지 않기 때문입니다.
Python은 행렬을 접근할 때 C와 같이 "Row-major order"를 따릅니다. 반면 Fortran은 "Column-major order"입니다.
[From Wiki]
이거 사실 외우기 힘들어요. 그래서 전 그냥 다음과 같이 외우고 있습니다.
어떤 행렬 A[m,n]에 대하여
Python: 오른쪽 n 먼저, 그리고 왼쪽 m 나중
Fortran: 왼쪽 m 먼저, 그리고 오른쪽 n 나중
그리고 먼저 나오는 행 또는 렬이 Loop에서 안쪽에 와야 합니다. 그래야 빨라요.
그래서 위 Fortran 예제는
DO 20 I = 2,M-1
DO 40 J=2,N-1
를 다음과 같이 바꿔야 합니다.
DO 20 J = 2,N-1
DO 40 I=2,M-1
그러면 결과는
1:0.02468sec 2:0.02476sec 3:0.02476sec
-0.227456830239
대략 0.025초 정도 되겠습니다.
마지막으로 한가지 더 테스트 해보겠습니다.
만약 Fortran도 Loop를 돌리지 말고, 행렬을 통째로 계산하면 어떨까요?
SUBROUTINE DFILTER2D(A,B,M,N)
C
DOUBLE PRECISION A(M,N)
DOUBLE PRECISION B(M,N)
INTEGER N, M
CF2PY INTENT(OUT) :: B
CF2PY INTENT(HIDE) :: N
CF2PY INTENT(HIDE) :: M
B=A
B(2:M-1,2:N-1)=B(2:M-1,2:N-1)+(A(1:M-2,2:N-1)+A(2:M-1,1:N-2) +
$ A(3:M,2:N-1)+A(2:M-1,3:N) )*0.5D0 +
$ (A(1:M-2,1:N-2) + A(1:M-2,3:N) +
$ A(3:M,1:N-2) + A(3:M,3:N))*0.25D0
B=B/4.D0
END
그러면 결과는
1:0.02951sec 2:0.02954sec 3:0.02955sec
-0.227456830239
오히려 위보다 못한, 약 0.03초 정도 되겠습니다.
4. Summary and Conclusion
[1000,2000] 크기의 행렬을 "Simple 2D Gaussian Filter"를 이용하여 필터링 시키면
| Program | Algorithm | Running Time |
|---|---|---|
| Python | By Element | 9.24 sec |
| Python | By Array | 0.094 sec |
| Fortran | By Element | 0.025 s |
| Fortran | By Element (Wrong Order) | 0.029 sec |
| Fortran | By Array | 0.030 sec |
Python을 쓰려면 행렬 차원에서 접근하시고, 최고의 성능이 필요하다면 Fortran을 이용해서 가장 원초적인 알고리즘으로 프로그램을 작성하시기 바랍니다. Loop 순서만 잘 지켜서요.
"""
제 개인적 목표는
- 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
soosoo님의 [Link & List] “한쿡세계시민” 32차 (66 Steemers)
stylegold님의 『주간어워드』 4차 (부제: 한주의 노력을 보답)
stylegold님의 『오마주』 프로젝트 4차 (부제: 숨겨진 글 발굴하자!!)
오마이갓!! 빙글빙글.ㅋㅋ
『주간어워드』 4차의 수상자가 되어 보팅 지원을 받으셨습니다. 매주 마다 다양한 프로젝트가 진행되니 참여부탁드립니다.^^
고맙습니다.
네 주시하도록 하겠습니다.