[F2PY] 02. Basic Example: Moving Average

in kr-python •  last year 

Screen Shot 2018-04-23 at 1.50.11 PM.png

F2Py란?

  • 수치 계산의 '우사인 볼트'라 할 수 있는 Fortran은 단점으로 다양한 파일의 입출력 같은 부가기능이 약합니다.
  • 다재다능 Python은 단점으로 단순 계산이 조금 느립니다.
  • 그래서 Python 뼈대에 계산 부분만 Fortran을 불러올 수 있게 도와주는 것이 F2PY입니다.
  • F2PYNumpy에 기본으로 내재되어 있습니다.
  • 이 포스팅은 기본적으로 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항목의 이동 평균을 계산하면,

ProgramAlgorithmRunning Time
PythonBy Array1.29 sec
PythonUsing Cumsum0.0145 sec
FortranBy Element0.91 s
FortranUsing Cumsum0.0127 sec

Fortran은 대체로 Python-Numpy보다 빠릅니다. 하지만 그보다 더 중요한 것은 최적화된 알고리즘입니다.


"""
제 개인적 목표는

  1. Object-Oriented Programming과 친해지기
  2. Github와 친해지기 입니다.

이 목표에 닿기 위해 일단 제가 나름 좀 아는 Python, 그 중에서도 NumpyMatplotlib로부터 시작하려 합니다.
"""

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

Authors get paid when people like you upvote their post.
If you enjoyed what you read here, create your account today and start earning FREE STEEM!
Sort Order:  

오, Fortran도 하시는군요!

아뇨, 사실은 포트란을 하는데 요새 파이쏜도 하는 거에요 ㅎㅎ

역시 싸이언티스트! ㅎㅎ

ㅋㅋ 저도 예전에 하드코어 계산할 때 포트란 엄청 돌렸었는데 요새는 다들 파이썬을 주로 쓰더라구요;

요즘은 매틀랩이나 메스메티카 패키지도 많이 개발되어 빠른 결과를 얻는게 아니면 종종 쓴다고는 하는데 아직도 많은 사람들이 포트란을 선호하고 있죠!

병렬 컴퓨팅도 혹시 다루시나요ㅎㅎ
다 추억으로 남기만 한 기억들이라 이젠 기억의 저편으로 ㅋㅋ;

저희 분야 아니면 사실 이제 포트란 사용자 찾기 힘들죠 ㅎㅎ
그러지 않아도 다음편에는 OpenMP 소개하려고 준비하고 있어요~

아.. 머리가 핑글핑글 돕니다.
예전에 잠시 프로그램을 건드릴랑 말랑 했었는데
... 취미로 배워볼만한게 뭐가 있을까여?

취미로 해도 무언가 목적이 있어야 손이 가지 않을까요? ^^
목적에 따라 천차 만별이긴 한데, 일단 가볍고 쉬운건 역시 Python이죠. 한단계 더 어려운걸로는 JAVA가 있고요. Text처리는 Pearl이 좋다는 말도 있구요.
웹페이지나 이쪽은 제가 잘 모르는데, 자바스크립트 이런게 많이 쓰이는듯?

음... 해보고 싶은 것은 있는데 라기보단 필요한 것은 있는데 간단한게 아니라서요.
뭐가 되었든 실제로 해봐야하는데 그게 제일 힘드네요 ㅎㅎ

맞아요. 저도 무언가 새로운 걸 시작하는게 점점 더 힘들어지네요.

어지러워요.ㅎ@@
존경스럽습니다.

충분히 이해합니다.
저는 그저 업무상 필요한 부분을 정리하고자 작성한 글이라 이렇게 여러 분들이 방문해주실지 몰랐어요 ^^