Jongrae

여기의 글은 더 이상 갱신되지 않습니다. 기존의 있던 (특히 python에 대한 설명)글을 새로 옮긴 곳에서 다 볼 수 있습니다. 새로 옮긴 블로그 주소는 

https://sites.google.com/site/tmptestcmd/

Quaternion을 이용한 자세 추정 칼만 필터입니다.

ipython notebook이 이번에 갱신되면서 jupyter notebook으로 바뀌었습니다. 물론 당분간 ipython notebook명령도 작동할것으로 생각됩니다만.. 이제부터는

>> jupyter notebook

이라는 명령으로 notebook을 실행하세요. 자세한 내용은..

http://blog.jupyter.org/2015/08/12/first-release-of-jupyter/

오랜만의 글입니다. 이번 글에서는 ipython notebook이라는 것에 대해 알아보겠습니다.

일단 실행을 하려면 (저같은 경우 Linux에서는) 명령창을 열고 command prompt에 다음과 같이 치면

$ ipython notebook

인터넷 browser에 (저같은 경우 google chrome에) 다음과 같은 창이 하나 생깁니다.

여기서 오른쪽 위에 있는 “new” 버튼을 누르고 “python3″를 선택하면 다음과 같은 창이 생깁니다.

여기에 다음과 같이 a = 1; b =2; c = a+b; print(c) 와 같은 python 프로그램을 쓰고

실행버튼 (|>)을 누르면 아래와 같이 실행된 장면을 볼 수 있습니다.

결론적으로 ipython에서 할 수 있는 모든것을 여기서 할 수 있습니다. 그렇다면 그냥 ipython을 쓰지 굳이 이렇게 notebook을 열어서 쓸 필요가 있을까요?

“필요없다”가 제가 한동안 가지고 있던 생각이었는데, 매우 잘못된 생각이었다는 것을 며칠전에 알았습니다. 

notebook은 이렇게 python 명령을 실행시킬 뿐만 아니라 매우 유용한 한가지 기능이 있습니다. 

각종 버튼이 있는 줄에 “code”라고 되어있는 버튼이 보이는데요. 이걸 누르면 “markdown”이라는 것을 선택한 후 여기에 각종 설명이나 수식을 (심지어 LaTeX 수식을 쓸 수도 있습니다) 쓸 수 있습니다. 아래 그림이 한 예입니다.

위의 그림 맨 밑줄에 있는 다항식은 “markdown” 상자에 다음과 같이 쓰여진것입니다. “$$” 로 “x(t) = -0.206t^4 + 0.9167 t^3 - ...” 둘러싸고 실행 버튼 (|>)을 누누누르면 위와 같이 예쁜 수식으로 바꾸어 줍니다.

notebook의 매우 유용한 기능은 바로 이것입니다. 매우 자세히 여러가지 생각을 프로그램과 함께 적어둘 수 있습니다. 보통 어떤 문제를 풀때 프로그램을 짜기 시작하면 초기 실행 하는 작은 프로그램을 만들고, 그것이 점점 커지고 그러다가 각종 변경을 하는데... 이걸 잘 관리하지 않으면 나중에 도대체 왜 내가 이렇게 했던가? 이 변경은 왜 했지? 이 프로그램하고 저 프로그램하고 뭐가 다른거지? 등등을 완전히 잊어버리게 됩니다.

notebook은 그러한 모든것을 한곳에 관리하게 해 줍니다. 아마도 절대 앞으로는 “이 부분은 왜 이렇게 바꿨지?” 하는 상황은 피할 수 있지 않을까요? 

그뿐만 아니라 자신이 짠 프로그램을 다른 사람에게 설명과 함께 보내줄때도 매우 편리합니다.

이번에는 FORTRAN으로 짜여져 있는 함수를 python에서 부를수 있게 하는 방법을 알아보겠습니다.

21세기에 무슨 FORTRAN이냐 하는 생각을 가질 수도 있지만, 과거에 많은 수치해석 프로그램이 FORTRAN으로 개발되었고 이 프로그램들이 여전히 여기저기에서 사용되고 있습니다. 

하지만 이미 개발이 완성된 함수를 다시 python으로 해석하는 작업은 매우 불필요하고 거추장스러운 일입니다. 이럴때를 위해서 그냥 FORTRAN 함수를 compile해서 일종의 library처럼 만들어 둔 다음 python에서 부르는 방법이 있습니다.

주어진 날짜와 위치에 따른 지구자기장값을 계산해주는 IGRF12 함수를 예제로 사용해보겠습니다. 원래 FORTRAN 프로그램은 여기서 내려받을 수 있습니다.

http://www.ngdc.noaa.gov/IAGA/vmod/igrf12.f

그리고 프로그램에 대한 설명은 여기에서 볼 수 있습니다.

http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html

igrf12.f 내부를 들여다 보면 지구자기장을 계산하는 함수가 

subroutine igrf12syn (isv,date,itype,alt,colat,elong,x,y,z,f)

라는 이름으로 되어있고 그 앞에 있는 내용은 이 함수를 실제 불러서 실행하는 부분입니다. 우리는 이 함수를 python에서 부를것이므로 앞부분 내용은 필요없고 이 함수부분만 필요합니다. 그래서, 앞을 다 잘라내고 위의 함수만 포함한 파일을 다음과 같이 만듭니다.

https://gist.githubusercontent.com/myjr52/62ca6ffc3e9c78ea0411/raw/df1ca256826579fe9810b1058ca3c2b57b928611/igrf12.f

한가지 위의 파일에서 내부에 추가된 부분이 있는데... 다음과 같은 내용입니다:

Cf2py intent(out) x Cf2py intent(out) y Cf2py intent(out) z Cf2py intent(out) f

이 부분이 FORTRAN 프로그램을 python에서 부를 수 있는 함수로 바꿔주는 f2py가 작동할때 subroutine igrf12syn에 있는 매개변수 중에 x, y, z, f는 출력변수라는 것을 알려주기 위함입니다.

operating system prompt에서 이제 f2py를 이용해 변형된 igrf12.f를 다음과 같이 compile합니다. (참고로 python을 저처럼 anaconda를 통해 설치했다면 f2py3라는 프로그램이 이미 같이 깔려있을겁니다)

> f2py -c -m igrf12_fortran igrf12.f

-c 는 compile 하라는 뜻이고, -m은 library 만들때 이름을 igrf12_fortran (특별히 다른 이름을 원한다면 여기서 이름을 바꾸면 됩니다)이라고 하라는 뜻입니다. 위를 실행하고나면 어떤 운영시스템을 쓰느냐에 따라 다른데...

linux에서는 igrf12_fortran.cpython-34m.so라는 library가 생겼고

mac os-x에서는 igrf12_fortran.so라는 library가 생겼습니다.

이를 부르는 방법은 

>> import igrf12_fortran

한 다음에

>> igrf12_fortran.igrf12syn(.....) 

하면 됩니다.

좀 더 프로그램을 간략하게 쓰기 위해서 다음과 같이 python 프로그램을 만들었습니다.

https://gist.githubusercontent.com/myjr52/62ca6ffc3e9c78ea0411/raw/6d8cbc3693f384a3e83f122b7a48155f71fb4b11/igrf12.py

이제는 python에서 다음 순서로 부르면 됩니다.

>> import igrf12          

>>  igrf12.igrf12(0,2016,2,4000,23,32)    

(-1216.8147739252843, 24122.237297648033, 118869.7813460386)

지난번에 만들었던 짝수를 찾는 함수입니다.

def is_even(number):   if ((number % 2 == 0) and (number != 0)):       return True   return False

이것을 이용해서 1부터 100까지 자연수중에 짝수의 루트값을 구해서 다 더하는 일을 해볼까요? 만약 이를 위해서 다음과 같이 프로그램을 짰다면

import math total = 0 for x in range(1,101):    if is_even(x):        total = total + math.sqrt(x)

이건 매우 C나 Fortran언어스러운 프로그램입니다. 이것보다는 다음과 같이 짝수의 루트값을 가지는 generator를 만들고

all_sqrt_even = (math.sqrt(x) for x in range(1,101) if is_even(x))

그 다음 합을 구하면

total = sum(all_sqrt_even)

됩니다.

이제까지는  generator를  “( )”를 이용해서 만들었는데... yield를 함수내에 사용해서 generator를 만들어 보겠습니다.

def get_even(number):    while number > 0:        if is_even(number):            yield number        number -= 1

이 함수는 get_even(10)하면 10보다 작거나 같은 자연수 중에 짝수를 포함하는 generator를 만듭니다. 다음과 같이 함수를 부르면

aa = get_even(10)

aa는 10보다 작거나 같은 짝수인 자연수를 포함하는 generator가 됩니다. 앞에서 이미 말했지만, 주의할 것은 list는 [10, 8, 6, 4, 2] 이렇게 5개의 숫자를 다 만들어서 5군데 저장장소에 보관해 두지만, generator는 필요할때 한개씩 값들을 만들어내고 다 쓰게 되면 텅빈 generator가 됩니다. 한개씩 값을 불러오는 방법은

next(aa)

10

이렇게 next() 함수를 쓰면 됩니다.  연속해서 부르면 8, 6, 4, 2를 출력하고 결국 StopIteration을 출력하면서 모든 요소가 다 소진되었다고 알려줍니다.

어떻게 이런일이 가능할까요? 이건 바로 매우 특별한 yield의 기능때문입니다. 보통함수는 (만약 yield가 아니고 return이라면) 값을 바깥으로 돌려주고 나면 내부의 상태가 초기화 되어 버립니다. 하지만  yield는 값을 돌려주고 또한 내부의 상태를 기억하고서 다음 next( )가 불리우면 전에 상태에서 다음 yield를 만날때까지 실행합니다. 

무슨 말인지 좀 더 자세히 알기위해 위에서 만든 get_even() 함수를 다음과 같이 수정했습니다.

def get_even(number):   while number > 0:               if is_even(number):           print(’before yield’)           yield number           print(’after yield’)       number -= 1       print(’after reducing the number’)

이렇게 수정한 후 

aa = get_even(10)

하고 나서 

next(aa)

하면

before yield Out[57]: 10

이 출력됩니다. 보다시피 위의 함수에서 “yield number”까지 실행되고 현재 number값이 10이므로 이 값을 돌려주고서 함수 실행이 멈추었습니다. 한번 더 

next(aa)

를 부르면

after yield (그전에 yield에서 멈췄으므로 바로 그 다음줄을 실행하고) after decreasing number (number를 하나 줄여서 9로 만들고) after decreasing number (9는 짝수가 아니니까 if안으로 안들어가고 다시 number를 하나 줄여서 8로 만들고) before yield (8은 짝수이니까 if안으로 들어가서 yield number를 만났으니 8을 돌려주고 거기서 멈추니까) Out[58]: 8 (8이 출력됩니다)

오늘은 generator에 대해 알아보겠습니다.

오늘 내용은 아래 링크에 있는 내용을 참조해서 좀 더 쉽게 만들었습니다. 

http://www.jeffknupp.com/blog/2013/04/07/improve-your-python-yield-and-generators-explained/

다음과 같은 경우를 상상해보죠: 주어진 자연수보다 작은 모든 짝수 자연수를 찾으려고 합니다.

그러면 짝수인지 아닌지 구별하는 함수를 일단 다음과 같이 만들 수 있습니다.

def is_even(number):    if ((number % 2 == 0) and (number != 0)):        return True    return False

이 함수를 다음과 같이 실행하면

>> is_even(10)

True

>> is_even(3)

False

주어진 값이 짝수이면 True, 홀수이면 False 값을 돌려줍니다.

그러면 예를 들면 10보다 작은 모든 짝수를 찾으려면 다음과 같이 list를 만들 수 있겠죠.

>> cc=[x for x in range(1,11) if is_even(x)]

그러면 cc = [2, 4, 6, 8, 10]이 됩니다.  그런데 이렇게 list가 짧지 않고 매우 길다면 어떻게 할까요? 예를 들면 range(1,10000000000000000001) 정도 되면 컴퓨터에 있는 메모리를 다 쓰고도 모자랄지도 모릅니다. 이렇게 list가 크면 메모리를 확보하는데 시간이 매우 많이 걸립니다. 대신에 다음과 같이 하면 

>> cc=(x for x in range(1,100000000000000001) if is_even(x))

놀랍게도 이를 실행하는데 거의 시간이 걸리지 않습니다. 어떻게 이런 일이 가능한 걸까요? 

비밀은 “[ ]”와 “( )”의 차이입니다. 이미 알고 있듯이 “[ ]”을 사용하면 list가 만들어집니다. list는 그야말로 그 안에 있는 모든 내용물을 다 만들어서 개별 메모리에 저장해둡니다.

하지만 “( )”를 사용하면 list가 아닌 generator가 만들어집니다. 

>> cc

하면 그전에 list였을때는 갖고있는 내용물을 다 보여주지만, 이제는 generator이기 때문에 단순히 다음과 같은 내용이 출력됩니다.

<generator object <genexpr> at 0x7f3110764fc0>

실제 안에 내용을 하나씩 보고 싶으면 next(...)라는 함수를 사용하면 됩니다. 즉,

>> next(cc)

2

>> next(cc)

4

>> next(cc)

6

와 같이 next(..)를 쓰면 generator cc가 가지고 있는 값에서 그전에 바로 불리워져서 쓰여진 값 바로 그 다음값을 보여줍니다.

즉, 실제로 1부터 100000000000000001까지의 모든 짝수가 한꺼번에 다 필요한 연산은 거의 없습니다. 대부분 각각의 값 하나 하나를 가지고 어떤 연산을 하고 최종 결과값만 알면 되는 경우가 많습니다. 예를 들면 그 숫자 범위내에 있는 모든 짝수의 값을 더한 결과값을 구하려하는 경우 하나하나 짝수값을 2부터 시작해서 만들어가면서 계속 더해나가면 되지, 모든 짝수값을 메모리에 다 불러둘 필요는 없습니다.

좀 더 작은 generator를 만들어보죠.

>> cc=(x for x in range(1,5) if is_even(x))

이제 next(..)로 값을 부르면

>> next(cc)

2

>> next(cc)

4

>> next(cc)

Traceback (most recent call last):

 File "<ipython-input-181-dea75668052a>", line 1, in <module>    next(cc)

StopIteration

여기서 보여주듯이 generator내에 값이 다 소진되면 StopIteraction이라는 오류가 발생하며 더 이상 값이 없다는 것을 알려줍니다.

아주 오랜만에 글을 쓰는데요. 이번 글부터 당분간은 python 문법자체에 집중해서 몇가지를 설명해 보도록 하겠습니다.

특히 저처럼 오래전에 프로그램 짜는법을 배우고 MATLAB같은 응용 프로그램에  대부분의 시간을 보낸 사람에게 약간은 생소한 프로그램 문법에 대해 알아보도록 하겠습니다.

그 첫번째가 list입니다.

참고로 이 내용은 아래 링크에 있는 내용을 참고 했습니다:

http://www.python-course.eu/python3_list_comprehension.php

다음과 같이 “[ ]”안에 내용물을 넣으면 list가 만들어 집니다.

>>  a = [1, 2, “abc”, [3,4,5]]

list의 각 내용물은 방금 예에서 보듯이 반드시 같은 자료 형태일 필요는 없습니다. a[0] = 1이고 a[2]=“abc”이고 a[3] = [3,4,5]입니다. 즉,

>> type(a)

list

>> type(a[0])

int

>> type(a[2])

str

>> type(a[-1])

list

입니다.

이미 만들어진 list를 이용해서 새로운 list를 만들수도 있습니다. a의 각 항을 제곱한 값을 갖는 list를 다음과 같이 만들면

>> b = [ x**2 for x in a]

되는데, 여기 “for x in a” 에서 x가 a에 있는 각각의 내용물로 순차적으로 대치됩니다. 즉, x=a[0] -> x = a[1] -> x = a[2] ...

그런데 위의 내용을 실행하면 오류가 발생합니다. 왜냐하면 a에는 정수만 있는 것이 아니고 문자열도 있고, list도 있는데 이것은 제곱이 안되기 때문입니다. 그래서 다음과 같이 하면

>> b = [x**2 for x in a[0:2]]

b  = [1, 4] 되는데 a[0:2]가 a[0]과 a[1]로  x값에 정수값 1과 2를 넘겨주므로 프로그램이 오류없이 실행됩니다. 이것을 좀 더 자동으로 인식해서 정수일때만 넘겨주게 하려면 다음과 같이 하면 됩니다.

>> b = [x**2 for x in a if type(x)==int]

이와 같이 하면 x의 자료형이 정수일때만 x값이 x**2 연산으로 넘겨지게 됩니다.

세개의 for-loop를 조합해서 쓸 수도 있습니다:

>> [(x,y,z) for x in range(1,31) for y in range(x,31) for z in range(y,31) if x**2 + y**2 == z**2]

이렇게 하면 1부터 30까지 있는 모든 정수의 조합에서 피타고라스의 정리를 만족하는 세개의 정수쌍이 생성됩니다.

2015_1_21(15h 2m 0s){0:04d}.tif'.format(idx)  선충의 이미지를 해석하는 아래 프로그램을 설명하도록 하겠습니다.

https://gist.github.com/myjr52/64b8b80da1e2b4efedde

우리가 그 동안 자주 사용했던 numpy와 matplotlib말고, 다음 두가지 함수 모음집을 이 프로그램에서 사용합니다.

from skimage.segmentation import boundaries from scipy import ndimage

특히 skimage는 여러가지 매우 유용한 그림 처리 함수가 구현되어있습니다. 자세한 설명은 아래에 가면 볼 수 있습니다.

http://scikit-image.org/

일단 네개의 그림 파일 이름이

2015_1_21(15h 2m 0s)0001.tif 

2015_1_21(15h 2m 0s)0002.tif   

2015_1_21(15h 2m 0s)0003.tif  

...

형식으로 되어있습니다. 이 그림을 한장씩 for-loop에서 읽기 위해 다음과 같이 파일 이름을 생성합니다.

file_name = '2015_1_21(15h 2m 0s){0:04d}.tif'.format(idx)

여기서 맨 오른쪽끝에 있는 idx가 1부터 증가하는데 {0:04d}부분이 이 숫자로 대체됩니다. 한개의 입력밖에 없으므로 {0: ...}에서 0이 첫번째 입력 idx를 의미하고 {0:04d}에서 04d는 index 정수 취급을 하는데(d가 정수 취급하라는 뜻입니다), 4자리에 출력하고 빈자리는 0으로 채우라는 (04부분) 뜻입니다.

matplotlib에 imread (image read의 줄임말?)라는 함수에 위의 함수 이름을 대입해서 해당 그림을 읽습니다.

Worm_Image = plt.imread(file_name)

이렇게 읽은 Worm_Image는 3차원 배열인데, [:,:,0], [:,:,1], [:,:,2]가 각각 빨간색, 파란색, 녹색의 정도에 대한 정보를 담고 있는 배열입니다. 즉, 읽은 원래의 그림은 색깔이 보통 있는데, 여기서 이 배열을 흑백그림으로 바꾸기 위해 이 세개의 배열을 다음과 같이 가중치를 두어서 2차원 배열로 만듭니다.

Worm_Image = (0.2989*Worm_Image[:,:,0]/max_val_img                + 0.5870*Worm_Image[:,:,1]/max_val_img                + 0.1140*Worm_Image[:,:,2]/max_val_img)

여기서 mag_val_imag는 원래 그림 배열이 갖고 있는 최대값입니다. 이 최대값으로 원래 그림배열을 나눠서 0과 1사이값으로 만들고, 그 앞에 가중치를 각각 0.2989, 0.5870, 0.1140으로 두어서 더했습니다. 이 가중치는 언듯보아서는 특별한 의미가 없어보이는데요. 색깔있는 그림을 흑백으로 바꿀때 사용되는 표준적인 가중치입니다. 어떻게 이 값이 정해졌는지는 아마도 그림처리 이론에서 알듯합니다.

여기서 이 그림처리의 큰 목적은 선충 몸의 경계선을 찾아내는것인데, 다음과 같이 간단히 find_boundaries라는 함수를 쓰면 됩니다.

Worm_Image = boundaries.find_boundaries(Worm_Image)

그런 다음 큰 폐곡선 안에 작은 폐곡선이 있을 수 있는데 이런 작은 폐곡선을 scipy에 있는 다음 함수로 모두 없애버립니다.

Worm_Image = ndimage.binary_fill_holes(Worm_Image)

공학이나 과학에서 시각적으로 보여주는 것과  시각자료를 제공하는 측정장비의 중요성이 나날이 커지고 있습니다.  그림 자료 처리는 다른 종류의 측정값에 비해 매우 까다롭습니다. 가속도계나 속도계와 같은 전통적인 측정장치는 한개의 특정한 물리량을 측정하지만, 카메라에 의해 주어진 그림 측정치는 무엇을 알고 싶으냐에 따라서 다양한 물리량을 뽑아낼 수 있습니다. 이렇게 간접적으로 그림에서 알고 싶은 물리량을 추출해 내는것은 때때로 매우 성가스러운 일입니다. 

중요한 것은 지난 20여년간 컴퓨터 계산 속도는 매우 빨라져서, 그전에 시간이 너무 많이 걸려서 할 수 없었던 그림 자료 해석이 가능해 졌다는 것입니다. 대부분의 기본적으로 필요한 그림 처리 방법이 python 함수 모음집, scipy와 skimage에 구현되어있습니다.

그림자료 처리를 통해서 얻어진 측정치의 좋고 나쁨은 거의 절대적으로 처음 주어진 그림의 성질과 얻으려고 하는 측정치의 명확한 정의의 존재 여부에 달려있습니다. 아예 처음부터 그림이 좋지 않거나, 그림에서 추출하려는 정보가 매우 미세한 성질이라면 문제가 매우 복잡해 집니다. 예를들면 수천장의 사진에서 흥미로운 사진을 추출한다라고 하면, 흥미롭다는 것 자체가 논리적으로 정확하게 정의하기가 쉽지 않아서 프로그램으로 구현하기가 쉽지 않습니다. 

그림의 좋고 나쁨에 따른 문제는 미리 처리를 좀 해주므로써 해결할 수 있습니다. 프로그램을 새로 짜거나, 이미 구현되어있는 함수모음집에 있는 기능을 쓸때, 할 수 있으면 미리 처리를 해서 그림을 좀 더 명확하게 만들 필요가 있습니다.

이렇게 미리 처리하는 것조차 자동화시킬 수 있지만, 아직까지 풀려고 하는 문제의 전문가보다 더 똑똑한 프로그램은 그렇게 많지 않습니다. 

(이 부분은 앞에서 말한 두번째 문제, 즉 측정치 정의의 명확성 문제와 비슷한데... 나중에 시간이 되면 기계학습이라는 분야의 문제를 다룰때 한번 경험할 수 있을것 같습니다)

여러가지 프로그램이 있겠지만 이렇게 그림을 미리 처리하는데 매우 유용한 프로그램이 ImageJ입니다.

http://imagej.net/

이 프로그램으로 동영상 파일을 읽은 후 명암, 밝기등을 조절해서 알고리듬에서 처리하기 쉽도록 영상을 만든 후 여러개의 tif 파일로 저장할 수 있습니다. 이렇게 저장된 tif 파일을 python에서 읽을 수 있습니다.  

즉, 다음과 같은 동영상 파일이 avi 양식으로 있을때

(이 그림은 C.elegans라고 불리우는 선충의 움직임을 현미경으로 본 동영상의 한 장면입니다. 영상은 대구경북과학기술원의 김규형 교수님 실험실 (http://home.dgist.ac.kr/khkim/) 에서 제공한것으로 사용허가를 해주신것에 대해 매우 감사드립니다.)

위의 영상을 ImageJ에서 읽은 후 밝기와 명암을 조절해서 다음과 같이 만든 후 tif 파일로 저장합니다.

이렇게 하면 배경과 선충의 몸체가 확실히 구별되어서 몸의 경계선을 찾아내는 것이 쉬워집니다.

아래 링크에서 프로그램과 그림 파일을 받아서 실행해보세요.

https://gist.github.com/myjr52/64b8b80da1e2b4efedde

이번에는 참고로 볼 수 있는 몇가지 자료를 소개하는 짧은 글입니다. 인터넷에 python에 관한 자료는 매우 많아서 어떤게 적당한지 판단하는것이 아주 쉬운일은 아닙니다. 그러므로, 아래에 제가 추천하는 자료는 저의 매우 개인적인 선택입니다.

미안스럽게도 추천 자료가 모두 영어로 되어있습니다. 영어 문서에 익숙하지 않다면 매우 불편할 수 있습니다. 사실 제가 이 글을 쓰기 시작한 이유도 영어문서에 익숙치 않은 한국어 사용자를 위한 자료가 부족하다는 것이었습니다.

좀 더 많은 사람들이 영어로 되어있는 자료를 한국어로 옮기는 혹은 자신의 지식을 한국어로 옮기는 작업을 하기를 개인적으로 바랍니다.

1) 공학분야에서 python을 사용한다면 이 문서가 매우 훌륭합니다. 아래 링크에 있는 pdf파일이나 html 링크를 누르세요

http://kitchingroup.cheme.cmu.edu/pycse/

2) 공학 문제나 과학 문제에 사용하는 것이 아니고 python 언어의 일반적인 기능에 대해 공부하고 싶으시다면 dive in python. 이 문서는 구구절절이 for-loop가 어쩌구 while이 어쩌구 하는 설명없이 바로 python 프로그램을 제시하고 한줄 한줄 설명해나갑니다. 제가 그동안 해왔던 방식의 설명도 이 자료에서 배운 것입니다.

http://www.diveintopython3.net/

3) 위의 두 자료는 기본적으로 프로그램의 기본은 되 있다는 가정에서 출발합니다. 그래서, 프로그램을 전혀 해본적이 없는 사람은 아마 따라가기가 힘들 수도 있습니다. python으로 처음 프로그램 짜는걸 배우려고 한다면 다음 자료를 시도해보는 것이 좋을겁니다.

http://www.greenteapress.com/thinkpython/

4) matplotlib는 그림그리는데 사용되는 함수 모음집인데 그 내용이 매우 방대합니다. 아래 자료에서 특별히 뒤쪽에 있는 예제가 매우 유용합니다. 링크에 있는 pdf 파일을 받으면 거기에서 3장, 4장을 읽으면 기본 사용에 문제가 없을겁니다. 특별히 어떤 기능을 하는 함수를 찾을때는 뒤쪽 III절에 있는 예제가 매우 유용합니다.

http://matplotlib.org/contents.html

(matlab 사용자를 위하여)

마지막으로 matlab을 이십여년 넘게 사용한 사람으로서 지난 12월부터 두달 채 못되는 기간동안 python을 numpy, scipy 그리고 matplotlib를 조합하여 사용해보았는데.. 수치 처리와 그림 그리는 기능에 있어서 matlab에 결코 뒤지지 않습니다. 확실히(!!!) matlab m-script는 python으로 완벽하게 대체 가능합니다.

이번에는 지난번에 이어서 아래 프로그램에 있는 그림 그리는 부분을 살펴보겠습니다.

https://gist.github.com/myjr52/141fbdee383c86a47999

plt.ion() 라는 함수를 설명하기 전에 먼저plt.ioff()를 먼저 설명하겠습니다. 이 함수에서 ioff는 interactive-off를 따서 붙인 이름입니다. 즉, 이 함수를 실행시키고 나면 다음과 같이

plt.plot([1,2,3])

해도 그림을 그린 창이 나타나지 않습니다. 즉각즉각 그림을 그리는 기능을 껐기 때문입니다. 한번 그림을 그리는데 시간이 많이 걸리는 데이타를 가지고 있을때, 그림을 그리고 거기에 수평축에 이름 붙이고 (plt.xlabel()), 수직축에 이름붙이고 (plt.ylabel()) 등등을 실행하는것 보다는; 그림에 관련된 모든 명령을 일단 실행시키고 나중에 한번에 다 그리는것이 시간이 절약됩니다. 즉, 순서는..

plt.ioff()

그림에 관련된 모든 명령을 여기서 다 실행한 후

....

plt.show() 

하면 최종적으로 그림을 한번만 그리게 됩니다.

이와 반대로 그때 그때 그림을 갱신하고 싶으면 plt.ion()하면 됩니다.  아마도 그때 그때 갱신하는 기능이 켜져있는것이 기본 설정일겁니다.

우리는 여기서 3차원 상자를 돌리는 효과를 보려고 하는것이므로 그때 그때 갱신 기능을 확실히 켜두기 위해 plt.ion()을 실행했습니다.

그림 창에는 번호가 붙는데; 다음과 같이 1번 그림창을 열었습니다.

fig  = plt.figure(1)

그리고 3차원 그림을 그리기 위해서

ax = fig.add_subplot(111, projection='3d')

그림창 fig에 '3d' 설정을 위와 같이 했습니다. 

(불행하게도 matplotlib에 3차원 그림 그리는 기능은 그렇게 훌륭하지 않습니다. 기본적인 그림은 가능하지만 매우 복잡한 그림을 그리기에는 아직 미흡한 점이 있습니다. 3차원 그림 그리는데 매우 훌륭한 함수 모음집으로 mayavi라는 것이 있는데요.

http://docs.enthought.com/mayavi/mayavi/mlab.html

안타깝게도 이 함수모음집은 아직 python 3로 구현되지 않았습니다. 이 함수모음집을 쓰고 싶으면 python 2.7을 사용해야합니다.)

for-loop 안에서 다음과 같이 상자를 그린 후

ax.plot(XYZ1[0],XYZ1[1],XYZ1[2],'b.') ax.plot(XYZ2[0],XYZ2[1],XYZ2[2],'r.') ax.plot(XYZ3[0],XYZ3[1],XYZ3[2],'g.')

3차원 그림을 보는 각도와 각 축별로 최소값과 최대값을 지정한 후

ax.view_init(30, 30) ax.set_xlim3d(-0.5,0.5) ax.set_ylim3d(-0.5,0.5) ax.set_zlim3d(-0.5,0.5)

다음과 같이 하면 위의 설정이 적용됩니다.

plt.draw()

그림을 그린 다음에 

plt.pause(0.01)

하면 0.01초 동안 프로그램이 멈추게 되는데, 이렇게 하는 이유는 그림이 화면에서 갱신되는 시간을 주기 위해서 입니다. 어떤 이유인지 모르겠지만 Mac OS X에서는 pause 하지 않아도 그림이 잘 갱신되는데 Linux CentOS에서는 이 명령을 넣지 않으면 그림이 전혀 바뀌지 않았습니다.

아뭏든 이렇게 컴퓨터 OS가 그림을 바꿀 시간을 준 후 

ax.clear()

하면, 다음에 그림을 그리는 명령이 실행될때, 그 전 그림을 지웁니다.

앞에서도 말했지만 matplotlib가 3차원 그림 그리는데는 함수 기능이 좀 불만족스럽고, 구현되어있는 함수도 불안정합니다. 이 프로그램도 제가 MacBook에 깔려있는 Anaconda에서는 잘 실행이 되는데, Linux CentOS에서는 그림이 갱신이 되지 않아서 움직이는 효과를 볼 수가 없었습니다.

지난번에 실행해 보았던 3차원 공간에서 상자를 그리고, 그걸 회전시키는 아래 프로그램을 해석해 보도록 하겠습니다.

https://gist.github.com/myjr52/141fbdee383c86a47999

3차원 공간에서 회전하는 물체의 각속도 벡터가 ${\boldsymbol\omega}$라는 $3\times 1$ 벡터로 주어졌을때, 방향 코사인 행렬의 미분 방정식은 다음과 같습니다.

$$\dfrac{dC}{dt} = - [{\boldsymbol\omega}\times] C$$

여기서

$$[{\boldsymbol\omega}\times] = \begin{bmatrix} 0 & -\omega_3 & \omega_2\\\ \omega_3 & 0 & -\omega_1\\\ -\omega_2 & \omega_1 & 0 \end{bmatrix}$$

이고, $\omega_1$, $\omega_2$, $\omega_3$는 각각 $x$, $y$, $z$ 축 방향 각속도입니다. 위의 방정식은 행렬 미분 방정식입니다. $C$는 방향코사인 행렬을 표현하는 $3\times 3$ 행렬입니다.

그런데, scipy를 포함해서 대부분의 수치적분기는 벡터 미분 방정식을 적분하도록 되어있어서, 행렬 미분 방정식은 직접 그 형태로 미분하지 못합니다. 그래서, 아래와 같이 미분 방정식을 정의할때

def DCM_ode(state, t)

state가 $C$의 값을 9차원 배열로 가지고 있다고 가정합니다. 함수 정의 내부에서 이 9차원 배열을 가지고 위의 미분 방정식을 구현할수도 있지만, 미분 방정식 자체가 행렬로 구현되어있으므로,

다음과 같이, 9차원 배열을 다시 행렬로 구성하여

DCM = state.reshape(3,3)

여기서 state는 9차원 배열이고, 이를 reshape 함수를 써서 $3\times 3$행렬로 만들었습니다.

그 다음 행렬 미분 방정식을 구성하기 위해, 다음과 같이 $[{\boldsymbol\omega}\times]$를 구성하고 난 후

wx = np.array([[ 0, -w3, w2],  [ w3, 0, -w1],  [-w2, w1, 0]])

아래와 같이 $dC/dt$를 계산합니다.

dDCM_dt = -np.dot(wx,DCM)

함수는 다시 벡터를 돌려줘야 하므로, 다음과 같이 dDCM_dt를 9차원 배열로 바꾸었습니다.

f = dDCM_dt.reshape(9)

(참고로 matlab 사용자를 위해서) matlab과 python에서 2차원 행렬의 1차원 index를 붙이는 순서가 다릅니다. 예를 들면,

A = np.array([[1,2,3],[4,5,6]])

A.reshape(1,6)

하면 

array([[1, 2, 3, 4, 5, 6]])

입니다. 즉 1행부터 시작해서 2행을 그 뒤에 붙입니다. 그런데 같은 행렬을 matlab에서 만들고 reshape(A,1,6)하면 [1 4 2 5 3 6]과 같이 1열부터 시작해서 2열을 붙이고, 3열을 붙이고 그리고 4열을 붙입니다.

참고로 python에서는 A[0] 하면 matlab에서 A(1)에 해당한다고 생각하여 1 이라고 생각하기 쉬운데, A[0]는 첫번째 행, [1,2,3] 입니다. A[1]은 두번째 행 [4,5,6]이고요. 

그러면 두번째 열을 빼내고 싶으면 어찌해야 할까요. 그건 matlab에서와 비슷하게 A[:,1]하면 됩니다.

그런데 사실, 행렬을 1차원 배열로 바꾸고 그걸 다시 행렬로 바꿔서 연산을 할때 어느 형태로 변환되는지 굳이 알 필요는 없습니다. python에서 다 알아서 일관되게 바꿔줄것입니다.

단지, 구체적으로 몇행 몇열에 있는 값을 배열에서 찾을때는 그게 어떻게 전환되는지 알고 있어야 겠죠. 하지만, 위의 함수에서 처럼 원래 행렬이 배열로 바뀐 변수가 들어오고, 그걸 함수안에서 다시 행렬로 바꾸어 연산할때는 그냥 크기만 맞춰서 reshape 해주면 됩니다. 굳이 어느 값이 어디로 가는지 알 필요는 없습니다. 

"# construct box" 부분은 직육면체 각 모서리를 100개의 점으로 표현하는 부분입니다.

이번에는 그 동안 해보지 않았더 3차원 공간에 그림을 그리는 방법과 이 그림을 지우고 다시 그리는것을 반복해서 마치 움직이는 것 같은 효과를 내는 프로그램을 짜 보겠습니다.

프로그램은 아래에서 받을 수 있습니다.

https://gist.github.com/myjr52/141fbdee383c86a47999

이 프로그램은 직육면체 모서리를 정의하는 점을 먼저 생성한 후, 이 점을 이용해 직육면체를 3차원 공간에 아래와 같이 그립니다.

임의의 각속도로 이 물체가 회전한다고 가정하고 이를 방향 코사인 행렬 미분 방정식으로 표현했습니다.

while 루프에서 주어진 시간을 짧은 구간으로 나눈 후, 각 시간구간마다 적분을 한 다음, 마지막에 얻어진 방향 코사인 행렬로 위의 물체를 표현하는 점들의 회전후 좌표를 구한 후, 위의 그림을 지우고 다시 그리는것을 반복합니다. 

이번에는 두가지 자질구레한 기술에 대해 알아보겠습니다.

1) 2진수, 8진수, 10진수, 16진수 변환

10진수 -> 2진수: bin(234) 하면 '0b11101010' 라는 234의 2진수에 해당하는 문자열(!)이 출력됩니다.

10진수->8진수: oct(234)하면 '0o352' 즉, 8진수에 해당하는 352가 출력됩니다.

10진수->16진수:  hex(234)하면 '0xea' 즉, 16진수에 해당하는 ea가 출력됩니다.

2진수->10진수: int('11101010',2) 하면 해당 10진수 234가 출력됩니다.

마찬가지로 x진수인 yyyy를 10진수로 바꾸려면 int('yyyy',x)하면 됩니다. 주의할점은 int의 첫번째 항은 항상 문자열이라는 것입니다.

ipython prompt에서는 다음과 같이 입력하면 십진수가 바로 출력됩니다.

16진수 11을 십진수로

In [159]: 0x11

Out[159]: 17

2진수 11을 십진수로

In [160]: 0b11

Out[160]: 3

2) 행렬형(matrix) 자료형태

python에서 행렬을 나타내는 자료형태는 배열 (array)입니다. 배열로 표현된 두 행렬 A, B를 A*B하면 항별곱이 됩니다. 행렬곱을 하려면 

import numpy as np

np.dot(A,B) 혹은 A.dot(B) 해야합니다.

그런데 이렇게 행렬곱이 매우 많다면 dot함수를 여러번 괄호안에 써야 하는데요. 예를 들면 A, B, C 세개의 행렬을 행렬곱하려면

np.dot(A,np.dot(B,C)) 아니면 A.dot(B.dot(C))

이거 매우 거추장스럽습니다. 이럴때 쓸 수 있는게 array가 아닌 matrix 입니다. 

다음과 같이 하면

A = np.random.rand(3,3)

3행 3열 무작위 값으로 구성된 배열이 생기는데요.

whos

해보면

A ndarray 3x3: 9 elems, type `float64`, 72 bytes

다음과 같이 A는 ndarray 형이라것을 보여줍니다.

이것을 matrix 자료형으로 바꾸려면

A = np.asmatrix(A)

하면 됩니다. 여기서 whos하면

A matrix [[ 7.49592988e-01 6.00<...>34e-04 9.16307122e-01]]

A가 matrix라는것을 보여줍니다.

마찬가지로 

B = np.asmatrix(np.random.rand(3,3))

한 후 A*B 하면 이 곱은 행렬곱입니다. 즉 자료형이 배열(array)이 아닌 행렬(matrix)이면 모든 연산이 행렬연산이 됩니다.

배열을 거치지 않고 처음부터 행렬형태로 만들고 싶으면

A = np.matrix('1 2 3;4 5 6; 7 8 9')

하고 print(A)하면

[[1 2 3] [4 5 6] [7 8 9]]

가 출력됩니다.

여기서 한가지 의문점을 가질 수 있습니다. 그럼 처음부터 불편하지 않게 행렬형만 쓰면되지 않을까? 이에 대한 답은 자신이 다루려고 하는 문제에 따라 그럴 수도 있고, 아닐수도 있다 입니다.

행렬형으로만 쓰면 행렬연산할때 표현이 간단해지기는 하지만, 많은 python 함수가 배열형태만 변수로 받아들이는 경우가 많습니다. 그러면 행렬형을 다시 배열형으로 바꾸어야 하는데 그 함수는, 쉽게 짐작할 수 있듯이, asarray입니다. 즉, 위에서 만든 행렬형태인 A를 배열형태로 만들고 싶으면, np.asarray(A)하면 됩니다.

참고로 python에서는 행렬형보다 배열형 자료형태를 쓰기를 권장합니다. 하지만, 자신의 문제에 행렬연산이 많다면 행렬형이 더 좋을 수도 있습니다. 그러므로, 둘을 편리한대로 섞어쓰는것이 좋다라고 말할 수도 있지만,  섞어쓸때 한가지 조심할점은 두 배열의 곱 A*B와 두 행렬의 곱 A*B는 매우 다르므로 주의해야합니다.

제 개인적인 의견은 "피할 수 없는 경우가 아니면 배열을 항상 사용하고, 반드시 행렬을 쓰고 싶은 부분이 있으면 매우 제한적으로 몇줄이내에서만 사용하고 그 부분의 프로그램에는 여기는 행렬을 썼다라고 설명문을 눈에 뜨이게 달아두는것이 좋다"입니다.

더 자세한 배열과 행렬에 대한 내용은 아래를 참고하세요.