[F2PY] 02. Basic Example: Moving Average

3년 전

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!
STEEMKR.COM IS SPONSORED BY
ADVERTISEMENT
Sort Order:  trending

오, Fortran도 하시는군요!

·

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

·
·

역시 싸이언티스트! ㅎㅎ

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

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

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

·

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

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

·

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

·
·

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

·
·
·

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

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

·

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