[F2PY] 02. Basic Example: Moving Average

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보다 빠릅니다. 하지만 그보다 더 중요한 것은 최적화된 알고리즘입니다.

---

"""
제 개인적 목표는

1. Object-Oriented Programming과 친해지기
2. 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