2016年9月18日日曜日

パラレルリンクロボットの逆運動学のお話

期末試験や大会などでゴタゴタしていたためか,4ヶ月ぶりの更新です.

現在自分はパラレルリンクロボットというロボットを勉強がてら作っています.以下のような機体です.


よく見かけるロボットアーム(多関節ロボット)とはだいぶ様相が異なっています.

最下部の移動する部分が常に地面に対して水平になるような構造となっています.

ちまちま開発していたのですが,疲れてしまったので気分転換がてらブログを更新することにしました.

実機製作の話やプログラムの話は今後に回すとして,今回はモーターの角度と移動する部分の位置の関係式を導出する逆運動学の話を書くことにします.

1.イントロダクション

実際に機体を作ったとしても,プログラムが書けなければ動かすことが出来ません.
その際にはこのロボットの動きを数式で表す必要があります,

ただ,機構学をまだ履修しておらず,運動学についてはさっぱりわからないので,こちらのPDFをトレースする形になりました.
The Delta Parallel Robot: Kinematics Solutions
Robert L. Williams II, Ph.D
http://www.ohio.edu/people/williar4/html/pdf/DeltaKin.pdf
パラレルリンクロボットの運動学について,数式や図を用いて解説しているPDFです.
多分僕の記事はどこか必ず間違いがあるのでこのPDFを見れば間違いないかと思います.


加えて,今回GitHubの方にリポジトリを作成しました.計算の殆どはsympyに計算は投げてまして,その過程も載せてありますので参考になるかと思います.
パラレルリンクロボットの逆運動学問題の導出過程
PLR_Controller/InversePositionKinematics.ipynb at master · TonyMooori/PLR_Controller

どういう機体であるのかという詳細な説明を書くべきですが,面倒なので省略します.上記PDFの最初で扱われているデルタロボットがそれに当たると思うのですがちゃんとした本を見つけられなかったのでなんとも…….

2.各座標の定義

実際に計算を始める前に座標系や各部の名称を決めていきます.まず各部の名称です.
上のモーターや電装部品がある部分を「胴体」,モーターの取り付け部分を「肩」,その先の関節を「肘」,更にその先を「手首」,一番下の部分を「移動体」と名前をつけます.
また,肩から肘の部分を上腕,肘から手首の部分を下腕と名付けます.
※一般的な名前ではなく僕が勝手につけた名前です

また各腕は,反時計回りに0,1,2と番号を振ります.

座標系は次のように取ります.胴体の中央を原点としていますが,胴体は移動しないためこれがワールド座標系となります.
胴体中央部を原点(0,0,0)としてとり,z方向が上方向が正である点に注意してください.

次に,各部分の原点からの位置ベクトルの名称を定義します.iは添字(0,1,2)です.
\( \vec{A_i} \) …… 肩の位置ベクトル
\( \vec{B_i} \) …… 肘の位置ベクトル
\( \vec{C_i} \) …… 手首の位置ベクトル
\( \vec{D} \) …… 移動体の中央部分の位置ベクトル
\( \vec{e_i} \) …… \( \vec{A_i} \)方向の単位ベクトル.互いに120度向きが異なります.
\( \phi_i \) …… \( \vec{e_i} \) の角度.\( \frac{2 \pi}{3}i \).
\( \vec{e_z} \) …… z方向の単位ベクトル

各部の長さは次のように定義します.
\( A \) …… 原点から肩への距離( \( \|\vec{A_i}\| \) ) 
\( B \) …… 肩から肘への距離(上腕の長さ, \( \| \vec{A_i} - \vec{B_i} \| \) ) 
\( C \) …… 肘から手首への距離(下腕の長さ, \( \| \vec{B_i} - \vec{C_i} \| \) ) 
\( D \) …… 手首から移動体の中央部分への距離( \( \| \vec{C_i} - \vec{D_i} \| \) ) 
サーボモータの角度は次のように定義します.
\( \theta_i \) …… i番目の腕につけられたサーボモーターの角度.±90°とする.正なら上腕が胴体よりも下に行く.


3.逆運動学

ではまず逆運動学から考えていきます.補足ですが,
・サーボモーターを動かしたら移動体の位置がどうなるかを求めるのが順運動学
・移動体の位置を決めたときのサーボモーターの角度を求めるのが逆運動学
と考えて問題ないと思います.

ネットの記事や1で挙げたPDFなどにも書かれているように,パラレルリンクロボットでは逆運動学の方が順運動学よりも求めやすいそうです.
逆に,シリアルリンクロボット(多関節ロボット)では逆運動学のほうが順運動学よりも求めやすいとのこと.

なのでとりあえず逆運動学問題について考えていきます.
ただ,繰り返すようですが運動学を学んでいないので高校レベルのベクトルで考えていきますのでこの道のプロがやる方法とは異なるかもしれないです.

まず\( \vec{e_i} \)について考えると,肩は3つ均等に並べられていることから,
\[ e_i = \left[\begin{matrix}\cos{\left (\phi_{i} \right )}\\\sin{\left (\phi_{i} \right )}\\0\end{matrix}\right] \]
となります.ただし \( \phi_{i} = \frac{2 \pi}{3}i \) , \( i \) は0,1,2のいずれかです.さらに,\( \vec{A_i} \)は \( \vec{e_i} \) と同じ方向で長さが異なるだけなので,
\[ \vec{A_i} = A \vec{e_i} \]
です.スカラーのAとベクトルのAが紛らわしいので注意してください.

次に \( \vec{B_i} \) について考えます.\( \vec{B_i} \) は \( \vec{A_i} \) から伸びており,\( \vec{e_i} \) と \( \vec{e_z} \) 成分しかありません.
したがって,\( \vec{B_i} \) は次のように表されます.\[ \vec{B_i} = \vec{A_i} + B( \vec{e_i}cos(\theta) - \vec{e_z}sin(\theta) ) \]
次に,\( \vec{C_i} \) を考える前\( \vec{D} \)を定義します.逆運動学問題では移動体の位置を決めた上で角度\( \theta_i \) を求めるので,移動体の位置は,\[ \vec{D} = \left[\begin{matrix}x\\y\\z\end{matrix}\right] \]という値が予めわかった上で考えていきます.

これを使うと\( \vec{C_i} \)は移動体がねじれずに(z軸に対して回転することなく)移動するため, \[ \vec{C_i} = \vec{D} + D\vec{e_i} \]となります.

ここで,2で書いた定義,\( C = \| \vec{B_i} - \vec{C_i} \|  \)を用いることで\( \theta_i \)が求まります.この式は求めたい変数\( \theta_i \)ひとつだけなので非常に簡単というか,個人的には驚きです.
この式を両辺を二乗すると次のようになります.
\[ C^{2} = \left(B \sin{\left (\theta_{0} \right )} + z\right)^{2} + \left(- A \sin{\left (\phi_{0} \right )} - B \sin{\left (\phi_{0} \right )} \cos{\left (\theta_{0} \right )} + D \sin{\left (\phi_{0} \right )} + y\right)^{2} + \left(- A \cos{\left (\phi_{0} \right )} - B \cos{\left (\phi_{0} \right )} \cos{\left (\theta_{0} \right )} + D \cos{\left (\phi_{0} \right )} + x\right)^{2} \] 更に展開して, \[ C^{2} = A^{2} + 2 A B \cos{\left (\theta_{0} \right )} - 2 A D - 2 A x \cos{\left (\phi_{0} \right )} - 2 A y \sin{\left (\phi_{0} \right )} + B^{2} - 2 B D \cos{\left (\theta_{0} \right )} - 2 B x \cos{\left (\phi_{0} \right )} \cos{\left (\theta_{0} \right )} - 2 B y \sin{\left (\phi_{0} \right )} \cos{\left (\theta_{0} \right )} + 2 B z \sin{\left (\theta_{0} \right )} + D^{2} + 2 D x \cos{\left (\phi_{0} \right )} + 2 D y \sin{\left (\phi_{0} \right )} + x^{2} + y^{2} + z^{2} \] となります.この辺はsympyという数式処理ソフトに投げてるので手計算ではないです.
この式と2,3分ほどにらめっこすると,こういう式であることに気づきます.
\[ P + Q \sin{\left (\theta_{i} \right )} + R \cos{\left (\theta_{i} \right )} = 0 \] ただし,P,Q,Rは次のような値です. \[ \begin{eqnarray*} P &=& - A^{2} + 2 A D + 2 A x \cos{\left (\phi_{i} \right )} + 2 A y \sin{\left (\phi_{i} \right )} - B^{2} + C^{2} - D^{2} \\ &-& 2 D x \cos{\left (\phi_{i} \right )} - 2 D y \sin{\left (\phi_{i} \right )} - x^{2} - y^{2} - z^{2}\\ Q &=& - 2 B z \\ R &=& - 2 A B + 2 B D + 2 B x \cos{\left (\phi_{i} \right )} + 2 B y \sin{\left (\phi_{i} \right )} \end{eqnarray*} \]

この式から,\( \theta_i \) を求めます.\( t = tan(\frac{\theta_i}{2}) \)を置くことで,この式を2次方程式にすることができ,そこから\( \theta_i \)を求めることが出来ます.

こちらはsympyに投げてるので途中式を省略します.結果は次のようになります. \[ \begin{eqnarray*} \theta_i &=& - 2 \operatorname{atan}{\left (\frac{1}{P - R} \left(Q - \sqrt{- P^{2} + Q^{2} + R^{2}}\right) \right )} \\ \theta_i ' &=& - 2 \operatorname{atan}{\left (\frac{1}{P - R} \left(Q + \sqrt{- P^{2} + Q^{2} + R^{2}}\right) \right )} \end{eqnarray*} \]
2次方程式,と言った時点で気づいた方もいるかもしれませんが,解が2つ出てきます.これは肘が内側に来てる状態も含まれているためです.

因みに,適切な方を選ばなければこのような状態になります.
それぞれの腕に対して2パターンあるので,この条件が成立する状態がこのロボットでは合計8パターン存在します.

加えて,どちらが適切であるかどうかは角度の絶対値が小さい方,すなわち0付近である方を選択することでたいていうまくいきます(理論的な根拠があるわけではないです).
そもそも不適切な方の角度はサーボモーターの移動範囲外である事が多いため,このような状態にはならないはずですが,間違えるとサーボモーターを傷めるので注意してください.

また根号の中が負となる場合がありますが,その場合はこのロボットの移動範囲外であるということを表していると思って問題ないと思います.

4.シミュレーション

シミュレーションというか,可視化というか,ちゃんと数式があっているかの確認を行います.
プログラムは個々に貼り付けると長くなるのでGitHubのリポジトリの方においておきます.

PLR_Controller/IPK_Visualize.ipynb at masterTonyMooori/PLR_Controller

結果はこのようになります.

必要ならば各ベクトルの距離を計算して調べてみるとよいかと思います.


5.実機による実験

Arduino Microを用いて3つのサーボモータを制御します.
Amazonで購入した大きめのトルクの物を動かしています.

回路図はこんな感じです.サーボとマイコンを繋げただけです.
プログラムはこれもGitHubにあげてあります.
PLR_Controller/Arduino/trace_circle at masterTonyMooori/PLR_Controller
実際の動画です.とりあえずサンプルとして円を描くような軌跡で動かしています.
今回は以上ですが,もうしばらく続くかと思います.



2016年5月15日日曜日

GUI自動化ライブラリの使い方のまとめと利用例[Python]

この所スランプ気味でして,気分転換に久しぶりの更新です.

1.GUI自動化について

普段PCを使っていて,単純作業やゲームなどを自動化したいと思うことは結構あるかと思います(チートはダメだよ).

ExcelなどではVBAでマクロを書くことが出来るのですが,それ以外の場合については大抵GUI自動化ライブラリを使って自分でプログラムを書くのが一般的かと思います.

有用に使う例は多々あるのですが,ただ遊ぶだけでもそこそこ楽しいです.頑張ればこんなこともできちゃいます.


自分は普段PyAutoGUIというPythonのライブラリを使っているのですが,日本語で使い方を書いたサイトが少なく,毎回英語の公式ドキュメントをGoogle翻訳片手に読むのはつらいのでブログにまとめることにしました.

2.PyAutoGUIの関数使い方のまとめ

自分が普段使うときに使う関数のまとめです.これである程度使いこなすことは出来るかと思います.

このプログラム自体に意味は無いので実行しないようにしたほうが良いかと思います.

ちなみに間違えて作ったプログラムで暴走した時は左上にカーソルを移動すると例外が発生して止まるようになっています.

import pyautogui as pa

"""
参考:PyAutoGUI公式ドキュメント
[1] Welcome to PyAutoGUI’s documentation! — PyAutoGUI 1.0.0 documentation
    https://pyautogui.readthedocs.io/en/latest/
[2] PyAutoGUI Documentation
    https://media.readthedocs.org/pdf/pyautogui/latest/pyautogui.pdf
"""

if __name__=="__main__":
    # 画面の幅と高さを取得する関数
    w,h = pa.size()
    
    # マウスカーソルの位置を取得する関数
    mouseX,mouseY = pa.position()
    
    # クリック
    pa.click()
    
    # click関数詳細
    pa.click(
        x = 200,y = 200,    # 座標
        clicks = 2,            # クリックの回数
        interval = 0.25,    # クリック間のインターバル
        button = "left",    # マウスのボタン(left,middle,rightの3つ)
        duration = 0.05        # 長押しするかどうか    
    )
    # マウスカーソルの移動(左上を(0,0)とする)
    pa.moveTo(x=100,y=200,duration=0.0)
    
    # マウスカーソルの移動(現在の位置を(0,0)とする)
    pa.moveRel(xOffset=20,yOffset=20,duration=0.0)
    
    # ドラッグ(左上を(0,0)とする)
    pa.dragTo(x=100,y=300,duration=0.0,button="left")
    
    # ドラッグ(現在の位置を(0,0)とする)
    pa.dragRel(xOffset=10,yOffset=30,duration=0.0,button="left")
    
    # マウスの押下(押したまんまになる)
    pa.mouseDown(x=200,y=300,button="left")
    
    # マウスの押下の取り消し
    pa.mouseUp()
    
    # スクロール(上に移動するのが正)
    pa.scroll(clicks=100,x=100,y=100)
    
    # 文字の打ち込み
    pa.typewrite("Hello, World!",interval=0.0)
    
    # 指定したキーを押す
    pa.press(keys="enter",presses=2,interval=0.0)
    
    # 指定したキーを押したまんまにする
    pa.keyDown("shift")
    
    # 指定したキーを離す
    pa.keyUp("shift")
    
    # 複数のキーを同時に押す
    pa.hotkey("ctrl","f")
    
    # スクリーンショット(ファイルに保存,範囲指定)
    # imageFilename,regionは指定しなくても良い
    img = pa.screenshot(
        imageFilename="screenshot.png",    # 保存先ファイル名
        region=(20,20,100,100)    # 撮影範囲(x,y,width,height)
    )
    
    # 指定した画像がスクリーン上のどこにあるかを探す関数
    # 一致箇所をすべて探してタプル(x,y,width,height)の配列で返す
    pos_list = pa.locateAllOnScreen(img)
    
    # ファイル名で与えることも出来る
    pos_list = pa.locateAllOnScreen("screenshot.png")
    
    # 一箇所見つかればよいのならばlocateOnScreen
    pos = pa.locateOnScreen("screenshot.png")


"""
keyDownで指定できるキー一覧
['\t', '\n', '\r', ' ', '!', '"', '#', '$', '%', '&', "'", '(',
')', '*', '+', ',', '-', '.', '/', '0', '1', '2', '3', '4', '5', '6', '7',
'8', '9', ':', ';', '<', '=', '>', '?', '@', '[', '\\', ']', '^', '_', '`',
'a', 'b', 'c', 'd', 'e','f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o',
'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z', '{', '|', '}', '~',
'accept', 'add', 'alt', 'altleft', 'altright', 'apps', 'backspace',
'browserback', 'browserfavorites', 'browserforward', 'browserhome',
'browserrefresh', 'browsersearch', 'browserstop', 'capslock', 'clear',
'convert', 'ctrl', 'ctrlleft', 'ctrlright', 'decimal', 'del', 'delete',
'divide', 'down', 'end', 'enter', 'esc', 'escape', 'execute', 'f1', 'f10',
'f11', 'f12', 'f13', 'f14', 'f15', 'f16', 'f17', 'f18', 'f19', 'f2', 'f20',
'f21', 'f22', 'f23', 'f24', 'f3', 'f4', 'f5', 'f6', 'f7', 'f8', 'f9',
'final', 'fn', 'hanguel', 'hangul', 'hanja', 'help', 'home', 'insert', 'junja',
'kana', 'kanji', 'launchapp1', 'launchapp2', 'launchmail',
'launchmediaselect', 'left', 'modechange', 'multiply', 'nexttrack',
'nonconvert', 'num0', 'num1', 'num2', 'num3', 'num4', 'num5', 'num6',
'num7', 'num8', 'num9', 'numlock', 'pagedown', 'pageup', 'pause', 'pgdn',
'pgup', 'playpause', 'prevtrack', 'print', 'printscreen', 'prntscrn',
'prtsc', 'prtscr', 'return', 'right', 'scrolllock', 'select', 'separator',
'shift', 'shiftleft', 'shiftright', 'sleep', 'stop', 'subtract', 'tab',
'up', 'volumedown', 'volumemute', 'volumeup', 'win','winleft', 'winright', 'yen',
'command', 'option', 'optionleft', 'optionright']
"""

幾つかマニアックな機能などは書きませんでしたが,この程度あれば十分いろんなことが出来るかと思います.


3.実際に使ってみる

関数の使い方だけではつまらないのでちょっと上の関数を使って遊んでみた例を書いてみます.

(あ,一応上で述べたミクのやつはここにコードが貼ってありますがPyAutoGUI以外の知識がいるのであまりオススメではないです)


とりあえずペイントで四角を書くプログラムです.本当はpa.locateOnScreen("画像ファイル")で座標を取り出してボタンをクリックするみたいな処理を入れたいのですが,画像がないと実行できないのでこれで妥協です.
import pyautogui as pa
import time

if __name__=="__main__":
    # 起動してすぐに実行されると困るので5秒待つ
    time.sleep(5)
    
    # 現在の位置を取得する
    x0,y0 = pa.position()
    
    # PyAutoGUIの各動作の間隔を0.1秒に指定
    pa.PAUSE = 0.1
    
    n = 256
    # マウスを押して
    pa.mouseDown(button="left")
    for i in range(n)[::16]:
        # 左上,右上,右下,左下,左上と移動してみる
        pa.moveTo(x0 , y0)
        pa.moveTo(x0 + i , y0)
        pa.moveTo(x0 + i , y0 + i)
        pa.moveTo(x0 , y0 + i )
        pa.moveTo(x0 , y0)
    # マウスを戻す
    pa.mouseUp()
またなんか作ったら記事書きます.


4.参考

公式ドキュメント.割とわかりやすくまとまっています.
Welcome to PyAutoGUI’s documentation! — PyAutoGUI 1.0.0 documentation
https://pyautogui.readthedocs.io/en/latest/

公式ドキュメントのPDF.一度に見れるので良いです.
PyAutoGUI Documentation
https://media.readthedocs.org/pdf/pyautogui/latest/pyautogui.pdf

2016年3月10日木曜日

ルンゲクッタ法の説明と実装例[C言語]

3日かけてルンゲクッタ法を理解する所までにきました.

とは言うものの,この記事では1階常微分方程式についての4次ルンゲクッタについて説明しています.

前回,前々回の記事である

オイラー法と修正オイラー法の説明と実装例
ニュートン・コーツの公式とシンプソンの公式の導出と実装

を読むと今回の記事で説明するルンゲクッタ法の導出部分が理解しやすいかと思います.


1.(4次)ルンゲクッタ法の簡単な説明


詳しい説明はしたところであまり興味がない人が多い気がするので,詳しい説明は最後にすっ飛ばして,実装するにおいて最低限の部分だけ示しておきます.

次の初期値がわかっている1階常微分方程式を解く事を考えます.
\[ \begin{eqnarray*} f(x_0) &=& y_0 \\ \frac{df(x)}{dx} &=& g(x,f(x)) \end{eqnarray*} \] この時,4次ルンゲクッタ法では次のようにして計算します. \[ \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta x \] ただしk1,...,k4は \[ \begin{eqnarray*} k_1 &=& g(x_k,f(x_k)) \\ k_2 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_1 \frac{\Delta x}{2}) \\ k_3 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_2 \frac{\Delta x}{2}) \\ k_4 &=& g(x_{k+1},f(x_k)+k_3 \Delta x) \end{eqnarray*} \]として計算しています.
この辺の導出は一番最後でしているので気になる人は見てってくださいな.


2.実装例


何を実装しようか悩んで適当にWikipediaを漁っていたらロジスティック方程式というモデルがあったので,これのシミュレーションをしてみたいと思います.

Wikipediaによれば,生物の個体数がどのように変化していくかを説明するモデルの一つであるとのこと.

求める常微分方程式はこれ. \[ \frac{dN(t)}{dt} = r \frac{K-N(t)}{K} N(t) \] rは一匹の個体の増加率,N(t)は時間tにおける生物の個体数,Kは個体数のしきい値で,最終的に安定するNの値.

例えば水槽なら最大でも20匹ぐらいしか入らない大きさならK=20みたいな値になると思う.

ちなみに厳密解は普通の変数分離法で多少計算が面倒だけど求まる.初期値をN_0とすると次のようになる. \[ N(t) = \frac{N_0 K e^{rt}}{K - N_0 + N_0 e^{rt}} \]厳密解はさておき,この常微分方程式にルンゲクッタ法を適用する.

まず見通しを良くするために,1で示した式のg(x,f(x))に当たる部分g(t,N(t))を定義しておく. \[ \begin{eqnarray*} \frac{dN(t)}{dt} &=& r \frac{K-N(t)}{K} N(t) \\ &=& g(t,N(t)) \end{eqnarray*} \] そうすると\( \Delta t \)を微小な値として,\( t_k = k \Delta t \)とすると,漸化式はこのようになる. \[ \int_{t_k}^{t_{k+1}} g(t,N(t)) dt \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta t \] ただしk1,...,k4は次の式で計算する. \[ \begin{eqnarray*} k_1 &=& g(t_k,N(t_k)) \\ k_2 &=& g(t_{k+\frac{1}{2}},N(t_k) + k_1 \frac{\Delta t}{2}) \\ k_3 &=& g(t_{k+\frac{1}{2}},N(t_k) + k_2 \frac{\Delta t}{2}) \\ k_4 &=& g(t_{k+1},N(t_k)+k_3 \Delta t) \end{eqnarray*} \] それではC言語による実装.関数g,を定義することで,数式とほぼ変わらないプログラムが書けます.実装の係数は適当なのでアレですが.
シミュレーション結果.
このシミュレーションだとほとんど違いがわかりませんね.誤差をプロットするとこうなります.
オーダーが1e-5といったところでしょうか.誤差が振動しているのがわかります.シミュレーションの精度としては十分なのでしょう.

4次のルンゲクッタ法では4次の項まで級数展開が一致しているとかで,誤差は\( O(\Delta t^5) \)に比例して小さくなるとのことです.

流石にルンゲクッタ法でシミュレーションすれば精度どうなの?とか言われても大丈夫そう……?まだその辺はよく知らないけれど.

3.ルンゲクッタ法の説明


手元の本には書かれていないのではっきりしたことはわからないのですが軽く調べたところによると,修正オイラー法は2次ルンゲクッタ法とも呼ぶらしいです.

なので2次ルンゲクッタについては前回の記事を参照してください.ここでは4次ルンゲクッタについて説明します.

今回はとりあえず1階微分方程式についての4次ルンゲクッタについて説明するので,次の常微分方程式について考えます.
\[ f(x_0) = y_0 \\ \frac{df(x)}{dx} = g(x,f(x)) \] ここで \( x_k = x_0 + k \Delta x \) と置くと \( f(x_{k+1}) \) と \( f(x_{k}) \) の漸化式を得ます. \[ \begin{eqnarray*} \int_{x_k}^{x_{k+1}} \frac{df(x)}{dx} dx &=& f(x_{k+1})-f(x_k) \\ &=& \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \\ \\ ∴f(x_{k+1}) &=& f(x_k) + \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \end{eqnarray*} \]

この \( \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \) の部分をどう計算するかによってオイラー法,修正オイラー法,ルンゲクッタ法が分類できます.

オイラー法では矩形積分,修正オイラー法では台形積分,ルンゲクッタ法ではシンプソンの公式を使います.

シンプソンの公式は次のようになります. \[ \int_{x_k}^{x_{k+1}}h(x) dx \approx \frac{1}{6}(h(x_k)+4h(x_k+\frac{\Delta x}{2}) + h(x_{k+1}) )\Delta x \] これを使って積分値を計算すれば良いわけです.\( x_{k+\frac{1}{2}} = x_k+\frac{1}{2}\Delta x \\ \) として \[ \begin{eqnarray*} g_k &=& g(x_k,f(x_k)) \\ g_{k+\frac{1}{2}} &=& g(x_{k+\frac{1}{2}},f(x_{k+\frac{1}{2}})) \\ g_{k+1} &=& g(x_{k+1},f(x_{k+1})) \\ \end{eqnarray*} \] という値を考えると, \( \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \) は次のように近似できます. \[ \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(g_k + 4g_{k+\frac{1}{2}} + g_{k+1} )\Delta x \]
しかしこの式には未知の値\(f(x_{k+\frac{1}{2}})\)や\(f(x_{k+1})\)を使っています.これらの値を矩形積分で求めて計算することもできますが,良い精度が得られません.

ルンゲクッタ法ではこれらの値を二段階の計算で計算します.
\[ \begin{eqnarray*} k_1 &=& g(x_k,f(x_k)) \\ k_2 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_1 \frac{\Delta x}{2}) \\ k_3 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_2 \frac{\Delta x}{2}) \\ k_4 &=& g(x_{k+1},f(x_k)+k_3 \Delta x) \end{eqnarray*} \] \( g_{k+\frac{1}{2}} \)に当たる部分をk2,k3という形で計算しています.これらを先ほどのシンプソンの公式に代入することで以下のような近似を行うものを4次のルンゲクッタ法と呼ぶそうです. \[ \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta x \]

今日はここまで.

ニュートン・コーツの公式とシンプソンの公式の導出と実装[C言語]

ルンゲクッタ法はどうやって導かれたのか謎だったのですが,『情報処理入門コース 数値計算』によれば,シンプソンの公式から来ているとのこと.
なのでシンプソンの公式がどういう経緯でできたのかメモしていきます.

※数式がうまく表示されない時はPCのブラウザから見ることを推奨します


1.ニュートン・コーツの積分公式


まずラグランジュ補間の公式が与えられているものとして使います.f(x)が近似したい式,p(x)が近似式.x_0,...,x_nが関数f(x)上のわかっている点. \[ p(x)=\sum_{i=0}^{n}f(x_i)N(i,x) \\ N(i,x) = \prod_{j=0,j \neq i}^{n} \frac{x-x_j}{x_i-x_j} \] ラグランジュ補間の出処はよくわかんないけど関数N(x)がクロネッカーのデルタ関数のようになっていることから意味をお察しください.とりあえずラグランジュ補間については前書いたのでそっちも見てくださいな.

次に,このラグランジュ補間で求めた近似式を用いて関数f(x)を区間[a,b]で積分する事を考えます.\[ a = x_0,x_1,x_2,...,x_n = b \]という様に区間をn分割し,上記のラグランジュ補間を適用すると次のようになります. \[ \begin{eqnarray*} S &=& \int_{a}^{b}f(x)dx \\ & \approx & \int_{a}^{b}p(x)dx \\ &=& \int_{a}^{b} \sum_{i=0}^{n}f(x_i)N(i,x) dx \\ \end{eqnarray*} \] ここで少しわかりづらいのですが,x_iは定数なのでf(x_i)はxの関数では無いことから, \[ \begin{eqnarray*} S & \approx & \int_{a}^{b} \sum_{i=0}^{n}f(x_i)N(i,x) dx \\ &=& \sum_{i=0}^{n}f(x_i) \int_{a}^{b}N(i,x) dx \end{eqnarray*} \] という式が得られます.ここで関数N(i,x)はx_0,...,x_nが[a,b]をn等分した値,すなわち \[x_k = a + k\Delta x \\ \Delta x = \frac{b-a}{n} \] と表せることを用いて,次の変数変換をすることを考えます. \[ t = \frac{x-a}{\Delta x}\\ dt = \Delta x \, dx \] ここで,N(i,x)は上記のtを代入して次のように変形したものをM(i,t)とする. \[ \begin{eqnarray*} N(i,x) &=& \prod_{j=0,j \neq i}^{n} \frac{x-x_j}{x_i-x_j} \\ &=& \prod_{j=0,j \neq i}^{n} \frac{(a+ t \Delta x) - (a + j \Delta x)}{(a+i \Delta x) - (a + j \Delta x ) } \\ &=& \prod_{j=0,j \neq i}^{n} \frac{t-j}{i-j} \\ &=& M(i,t) \end{eqnarray*} \]このN(i,x)を変形した関数M(i,t)を用いると先ほどの積分はこうなる. \[ \begin{eqnarray*} S & \approx & \sum_{i=0}^{n}f(x_i) \int_{a}^{b}N(i,x) dx \\ &=& \sum_{i=0}^{n}f(x_i) \Delta x \int_{0}^{n}M(i,t) dt \end{eqnarray*} \] M(i,t)は変数xや,区間[a,b]などに依存しないから,\[ C_i = \int_{0}^{n}M(i,t) dt \]と置いて,先に手計算で計算しておけば使い回しが効いて,最終的に次のようになります. \[ S \approx \sum_{i=0}^{n}f(x_i) C_i \Delta x \] この公式をニュートン・コーツの積分公式と呼ぶとのこと.


2.シンプソンの公式


うえで求めたニュートン・コーツの積分公式にn=2を代入したものが世間一般にシンプソンの公式と呼ばれているものとのこと.C_iの値は上式に突っ込んで \[ \begin{eqnarray*} C_0 &=& \int_0^2 \frac{t-1}{-1} \frac{t-2}{-2} dt &=& \frac{1}{3} \\ C_1 &=& \int_0^2 \frac{t}{1} \frac{t-2}{-1} dt &=& \frac{4}{3} \\ C_2 &=& \int_0^2 \frac{t}{2} \frac{t-1}{1} dt &=& \frac{1}{3} \\ \end{eqnarray*} \] となる.この公式のすごいところはこの係数をあらかじめ求めておけば面倒な計算を端折ることができるところだと思います.
この値をニュートン・コーツの積分公式に代入すると, \[ \begin{eqnarray*} \int_{a}^{b} f(x) dx &=& \sum_{i=0}^{2}f(x_i) C_i \Delta x\\ &=& (\frac{f(x_0)}{3} +\frac{4f(x_1)}{3} +\frac{f(x_2)}{3}) \Delta x\\ &=& (\frac{f(x_0)}{3} +\frac{4f(x_1)}{3} +\frac{f(x_2)}{3}) {\frac{b-a}{2}}\\ &=& \frac{1}{6}(b-a)(f(x_0)+4f(x_1) +f(x_2)) \end{eqnarray*} \] となります.

シンプソンの公式を図で表すとこんな感じ.青色の線が二次式で近似した曲線.その曲線の積分を求めてる.

ちなみにn=3の場合などはニュートン・コーツの公式のWikipediaに載ってたりするのでそちらを参照すると良いかと思います.


3.プログラム


抽象的な話はこの辺にして,実際にプログラムを書いてみます.試しに\( f(x) = \frac{4}{x^2+1} \)を[0,1]区間で積分してみます.

ちなみに\( \int_{0}^{1}f(x)dx \)の厳密解はπです.

まず \( \Delta x = \frac{1}{n} \) とし,\( x_k = k \Delta x \) と置くと, \[ \begin{eqnarray*} \int_{0}^{1}f(x)dx = \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}}f(x) dx \end{eqnarray*} \] という微小区間 \( [x_k,x_{k+1}] \) の積分となります.この微小な区間に対してシンプソンの公式を適用すると次のようになります. \[ \int_{x_k}^{x_{k+1}}f(x) dx \approx \frac{1}{6}(f(x_k)+4f(x_k+\frac{\Delta x}{2}) + f(x_{k+1}) )\Delta x \] 実際はf(x_k)の足し算がかぶっているのでもう少しスマートにできるのですが,今回は愚直にそのままこれを実装します.
実行結果
結果が3.141593なので円周率が概ねちゃんと求まっているので大丈夫そうです.

オイラー法と修正オイラー法の説明と実装例[C言語]

学科的に流体シミュレーションとか多分そのうちやると思うのですが,ルンゲクッタ法を理解していなかったのでとりあえず常微分方程式の数値計算の基本であるオイラー法と修正オイラー法について記事を書いてみます.


1.オイラー法


オイラー法は精度は悪いですが,非常に簡単で計算が高速です.研究などでは恐らく使われることはありませんが,ゲームプログラミング等をしているとシューティングの弾の運動などは大抵オイラー法が使われています.
今回は説明のための例として,次のような常微分方程式があったとします.今回は初期値f(x_0)がわかっているものとします. \[ f(x_0) = y_0 \\ \frac{df(x)}{dx} = g(x,f(x)) \] ここで \( x_k = x_0 + k \Delta x \) と置くと \( f(x_{k}) \) の値は次のようにして求められることがわかります. \[ \begin{eqnarray*} f(x_{k}) &=& f(x_0) + \int_{x_0}^{x_k} g(x,f(x)) dx \\ &=& f(x_0) + \sum_{i=0}^{k-1} \int_{x_{i}}^{x_{i+1}} g(x,f(x)) dx \end{eqnarray*} \] ここで \( f(x_{k+1}) \) と \( f(x_{k}) \) の漸化式を考えます. \[ \begin{eqnarray*} \int_{x_k}^{x_{k+1}} \frac{df(x)}{dx} dx &=& f(x_{k+1})-f(x_k) \\ &=& \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \\ \\ ∴f(x_{k+1}) &=& f(x_k) + \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \end{eqnarray*} \] この \( \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \) の部分がキモです.この部分をどうするかで精度が大きく変わってきます.

オイラー法ではg(x,f(x))が区間\([x_k,x_{k+1}]\)でほとんど変化しないものとして,次のように近似しています. \[ \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx g(x_k,f(x_k)) \Delta x \] これを漸化式に代入したものがオイラー法です. \[ f(x_{k+1}) \approx f(x_k) + g(x_k,f(x_k)) \Delta x \]

2.修正オイラー法


オイラー法は非常に簡単なのですが精度がよろしくないです.なのでもう少し精度を上げる必要があります.

修正オイラー法はその改良版の一つです.もっとも,これもあまり制度が良いというわけではないですが.

先ほど, \( \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \) の部分がキモと言いました.オイラー法では関数g(x,f(x))があまり変化しないものとして計算していました.イメージとしてはこんな感じです.
矩形近似とか長方形近似とか呼ばれているタイプの積分方法です.見るからに精度がよろしくなさそうなのがわかります.

そこで台形公式というものを使って修正オイラー法というものが提案されました.イメージとしてはこんな感じです.
名前の通り台形で近似するというものです.数式で表すと次のようになります.台形の面積は((高さ)×(上辺)+(下辺))÷2という小学生の時に習った公式を用いて,
\[ \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx (g(x_k,f(x_k)) + g(x_{k+1},f(x_{k+1}))) \frac{\Delta x}{2} \] という近似式で計算します.しかしこれを漸化式に代入すると \[ f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},f(x_{k+1}))) \frac{\Delta x}{2} \] となり,\( f(x_{k+1}) \) を求めるために \( f(x_{k+1}) \) を使っているのでこれでは使えません.なので右辺の \( f(x_{k+1}) \) だけは,オイラー法で使った矩形近似を用いて次のように求めておきます. \[ \widehat{f(x_{k+1})} = f(x_k) + g(x_k,f(x_k)) \Delta x \] この \( \widehat{f(x_{k+1})} \) を用いて漸化式に代入したものが修正オイラー法です. \[ f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},\widehat{f(x_{k+1})})) \frac{\Delta x}{2} \]

3.利用例

ここまで抽象的すぎたのでここで一度,次の1階微分方程式を数値的に解いてみたいと思います. \[ \frac{df(x)}{dx} - 2f(x) = x \] まぁ適当に作った式なので意味は無いです.ちなみに厳密解はCを定数として \[ f(x) = -\frac{1}{4}(2x+1) + Ce^{2x} \] と表されます.ちなみに定数Cは関数f(x)の初期値f(0) = aとした時,\( C = a + \frac{1}{4} \)になります.
今回は初期値 f(0) = a = 1 として数値的に求めてみます.

上記の微分方程式を変形してみます. \[ \begin{eqnarray*} \frac{df(x)}{dx} &=& x + 2f(x) \\ &=& g(x,f(x)) \end{eqnarray*} \] これでg(x,f(x))が求まりましたので,先ほど示した漸化式に代入してみれば,オイラー法では \[ f(x_{k+1}) \approx f(x_k) + (x_k + 2f(x_k)) \Delta x \] となります.

修正オイラー法では,まず \( \widehat{f(x_{k+1})} \) が \[ \begin{eqnarray*} \widehat{f(x_{k+1})} &=& f(x_k) + g(x_k,f(x_k)) \Delta x \\ &=& f(x_k) + (x_k + 2f(x_k)) \Delta x \\ \end{eqnarray*} \] となることから,次の漸化式に代入してやれば良いです. \[ f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},\widehat{f(x_{k+1})})) \frac{\Delta x}{2} \\ \] 全部数式を代入すれば,x_kとf(x_k)とΔxだけの式にできますが,数式が長くなるので省略.

もっとも,関数をプログラム内で定義してやればそんなに気にならないと思います.


4.精度の評価


ソースコードは長くなりますので先にオイラー法,修正オイラー法の精度を示しときます.

もっとも,評価と言ってもこの例での真値との差を見るだけですが.以下数値的に求めたものと,真値との差です.
刻み幅Δx=0.001,反復回数1000回です.euler's method 2の方が修正オイラー法です.適切な英語が見つからなかったので…….
1000回繰り返して差が10^-7オーダーですね.とりあえずそれらしい値は取れているのがわかります.
これを見る限り修正オイラー法はオイラー法と比べたら優秀そうではあります.
ただ傾きが指数関数的に増えていきそうなので後日ルンゲクッタ法を実装してそれと比較したいです.


5.ソースコード


C言語による実装.オイラー法はこんな感じです.
修正オイラー法はこんな感じです. 今日はここまで.

2016年3月7日月曜日

自己相関関数で音声の繰り返し部分を抽出する[Python]

フーリエ変換でミクを描いたものの,フーリエ変換をちゃんと勉強したことがなかったので『メカトロニクス入門シリーズ 信号処理入門』を読んでます.

その中で「自己相関関数」というものが紹介されており,非常に面白そうだったので実装してみることにします.


1.自己相関関数とは

「あー」というような定常な音の場合,同じような音の波形が繰り返されています.
すなわち,f(t)を時間tの時の音の大きさだとすると,f(t)には次のような性質を持つことがわかります. \[ f(t) = f(t + nT) \] もっとも,実際の音声では完全にイコールにはならないと思いますが,概ねこのような性質を持っています.こういった関数を周期Tの周期関数と呼びます.例えばサイン関数は周期2πの周期関数です.

この周期Tを求める手法として上記の本では「自己相関関数」を用いています.関数fに対する自己相関関数R_ffを,今回次のように定義します.
\[ R_{ff}(\tau) = \frac{(f(t),f(t+\tau))}{||f(t)||\,\,||f(t+\tau)||}\\ \] ただし,実装するときにどうするかをはっきりしたいので,関数の内積とノルムは次のように定義します. \[ (f(t),f(t+\tau)) = \int_{t_0}^{t_1} f(t)f(t+\tau)dt \\ ||f(t)|| = \sqrt{\int_{t_0}^{t_1}f(t)^2 dt }\\ ||f(t+\tau)|| =\sqrt{ \int_{t_0}^{t_1}f(t+\tau)^2 dt} \\ \]
関数f(t)とf(t+τ)の相関係数みたいなものなので,関数のノルムで割るかどうかや,関数の平均値を減じるかどうかは人それぞれだと思います.

さてこの自己相関関数R_ffを実際に計算していくとこのようになります.
τを変化していき,τ=nT ( n = 0,1,2,... )となる点で似たパターンとなるため,R_ffが極大値を取ることがわかります.

2.抽出アルゴリズム

ソースコードを貼る前に,実際にプログラムを組んで「あー」という音を読み込ませて,その音声に対する自己相関関数を求めるとこんな感じになります.
とまぁ現実の音をやってみるとそんなうまくいかないようです.結構ギザギザしてて嫌な感じです.

恐らくフラクタル図形のように「あー」という音の中にも似た部分があるようです.それはそれで面白いのですが…….

まぁ自分の目で見ればどれが目的の周期なのか一目瞭然なのですが,できれば自動的に判定したいところです.

とりあえずFFTで高周波成分を除去するハイパスフィルタを実装するとこんな感じになります.緑色の部分がハイパスフィルタで高周波成分を除去したもの.移動平均でも良いと思います.
そんでもってscipy.signalの極大値の要素番号をsignal.argrelmaxで求めて,しきい値0.2よりも大きい相関以上かつ極大となっている要素番号を取得すると,{5,  646, 1291, 1938, 2582, 3228, 3872, 4520} みたいな等差数列になっているので差を求めてやればそれっぽい周期が求められる.

その周期で音声の波形を10個ほど取り出すとこんな感じになった.
繰り返されている部分がいい感じに取り出せました.とは言うものの周期Tをint型に型キャストしているのでズレていってしまいますがまぁこんなもんでしょう.

この繰り返されている波形データだけで音声認識とかできないか考え中なところです.

3.プログラム

コメントは多めにしたつもりではありますがうまくいかない場合はコメントを参考にしつつ頑張って下さい.

2016年2月24日水曜日

水玉コラ作成ソフトを公開した話

Bloggerで過去の記事に追記をすると最新の記事の扱いになる現象が勃発している今日この頃です.先日「水玉コラ作成ソフト」をVectorさんで公開しました


一応こんな感じに動作する.
これを使うことで慣れれば2,3分で水玉コラを作ることができるようになる.やったぜ.

一度フリーソフトを公開してみたかったので夢というか目標が達成できてよかった.Vectorさんもフリーソフトの審査を人の目が見てくれているようでとても良い対応だったのでよかった.

ただ現状ダウンロードして実行すると次のような画面が出る……(´・ω・`)
しばらくの辛抱らしい.

Windows SmartScreenは認識されないアプリの起動を停止しました……と.一応弁明させていただきますと僕がこの画面を作ってパスワードを盗み取ろうとしているわけではないです.

この画面,沢山ダウンロードされると消えるらしいのだけど俺だったらこんなソフト使わないな(笑).こういう詐欺もあるだろうし,気をつけねば.

フーリエ級数の可視化プログラムの実装[Processing]

もはや2週間前のネタを真似ただけで新規性に乏しいのですが……(´・ω・`)
一応実装にあたって元ネタを書いておかないと申し訳ないのでツイートを貼らせていただきます……



という経緯です.

今回の実装はProcessingによるものなので上記ツイートのような数式で描くことができませんが,プログラミング言語というだけあって自由度は高いものとなっています.

円の数とか任意に決められますし.関数も自由に決められますし.

追記:鰺坂もっちょ先生がトゥギャってくれました→世はまさに大フーリエ時代 - Togetterまとめ

・ソースコード


例として矩形波を近似しています.関数funcの中身を[0,1)区間で変えてやればその形になりますが,無限大に発散させたりするとダメです.
あと数字が小さすぎる(1.0とか)だと小さすぎて見えませんので気をつけてください.


2016年2月17日水曜日

画像を一筆書きで描くプログラムを実装した[Python,Processing]

1.概要

この動画の最後のやつの,ほぼ完璧なトレースができます.
 

とりあえず順序としては
1.画像を読み込む
2.エッジ部分の座標を取り出す
3.巡回セールスマン問題として一筆書き経路を探索
4.フーリエ級数展開
5.画像を描画
という手順です.

今回サンプルとして画像はこいつを使います.勝手に使ってください.
Processingしかできない人や,Pythonの実行環境ないひとは一筆書きの経路を探索した結果の座標ファイルを置いておきます.こいつをProcessingのソースコードと同一フォルダに入れておけば一応はできます.ちなみにこんな感じ.


プログラムが長過ぎるので注意してください.もっとスマートにかければ良いのだけれど,無駄にクラスを作成して再利用しやすい形にした……(´・3・`)ドウセツカワナイノニ

2.プログラム

まず上記の手順の1-3を行うソースコードです.
蟻コロニー最適化のプログラムなのでこの前書いたやつの書き直しみたいな感じです.
実際はopt2を使って更に最適化してますが,それのコードはネットのやつを丸パクリしたやつなので流石に上げられません(参考URL先に借りたソースコードがあります).


次にProcessingで上記の手順の4,5の描画を行うコード

初音ミクを三角関数で描いてみた

1.動画紹介

フーリエ級数の視覚化を見て動画を作った.



まぁ二番煎じなんですけど.その辺の詳細はこちらで.

2.手法とか

まず序盤のやつはこいつ.Processingで作ったフーリエ級数の視覚化プログラム.こいつでいろんな関数を近似している.ここの先で先駆者たちのツイートが見られます.
フーリエ級数の可視化プログラムの実装[Processing]


次に出てくる二次元平面での三角関数的な奴ってのはこれ.実はフーリエ変換ではなく離散コサイン変換の二次元の基底関数なので嘘といえば嘘だけど,まぁそんなに違いはないし説明が面倒なので省略した.
離散コサイン変換(DCT)をPythonで実装した[Python]


その次の巡回セールスマン問題についてはこいつで解いた.基本的に蟻コロニー最適化を用いているけど,2000個の点で蟻コロニー最適化をすると遅すぎるので,パラメータを変えて貪欲法みたいな方法で経路探索して,その下のサイトで紹介されているopt2法のプログラムを用いてる,

・巡回セールスマン問題を蟻コロニー最適化と遺伝的アルゴリズムで解いてみる[Python]
Algorithms with Python / 巡回セールスマン問題


そんで最後のやつはフーリエ級数の可視化プログラムを使って,巡回セールスマン問題でもとめた一筆書きのルートを描いてる.
・画像を一筆書きで描くプログラムを実装した[Python,Processing]

結構大掛かりなプログラムだった割に動画を見てみると正直あまり面白くない……(´・ω・`)
アイデアおもついて実装したら面白くないのよくあるから仕方ないか…….

2016年2月14日日曜日

巡回セールスマン問題を蟻コロニー最適化と遺伝的アルゴリズムで解いてみる[Python]

巡回セールスマン問題(Traveling Salesman Problem,TSP)は有名だけれどなにげに解いたことがなかったのでちょいと実装してみました.
手法としてはいろいろあると思うのですが,今回は遺伝的アルゴリズム(GA)と蟻コロニー最適化(ACO)の2つを実装・比較しました.

1.巡回セールスマン問題

よく出される例としてはゴミ収集車の経路でしょうか.ゴミ収集所から出発して,100軒の家を最短で全て回ってゴミ収集所にまた戻ってくるまでの経路を求める話です.

まぁ僕がここで云々説明せずともネット上には沢山関連するサイトが有りますからそちらを参照してください.とりあえずWikipediaと検索したらトップに出てきたサイトのURLを貼っときます.

真ん中のサイトは手法について非常に詳しく解説されている上,Pythonで実装されたソースコードがありますのでそちらもごらんになると良いかと思います.

一番下のサイトは巡回セールスマンについての説明が詳しく,フリーソフトの紹介もされています.


この問題はまだ多項式時間で最適解を解くことができるアルゴリズムが無く,うちの大学でも研究してる研究室がたしかあった気がします.

2.蟻コロニー最適化による実装

蟻コロニー最適化に関する数式は今回全てWikipediaに載っていたものを利用しました(蟻コロニー最適化 - Wikipedia).
一応『心はプログラムできるか』という本を読んだ時に出てきてまして非常にわかりやすい解説がなされていましたので読んでみると良いかと思います.

以下ソースコード
とりあえず評価のちゃんとした結果は遺伝的アルゴリズムの実装の後に書きます.

3.遺伝的アルゴリズムによる実装

ぶっちゃけ蟻コロニー最適化よりも精度低いです.30個点でさえも蟻コロニー最適化に負けています.1000個ぐらいの点になるともう最適解求めるなんて夢のまた夢でした.
原因としては恐らく突然変異がランダムすぎて,ほとんど改悪になるばかりという結果でしょうか.遺伝子は巡回経路の順番を表した変数(5つの点なら[2,4,3,1,0]という配列)です.突然変異で中身を入れ替えてます.
※追記:先日図書館にてほかの人の実装例を見ましたがこれとは異なる手法を用いていました.それを踏まえるとこのプログラムは一般的な実装例とは異なると思われますので注意してください.

ソースコード

4.評価

一応同じデータを用いてやった結果です.30個の[0,1)のランダムな二次元の点について行いました.
パラメータはデフォルトです.まずは遺伝的アルゴリズム.
直線が交差してますね.コストは6.056765でした.続いて蟻コロニー最適化.
直線が交差すること無く,コストも5.701502です.しかも計算時間が遺伝的アルゴリズムよりも遥かに実行時間が早いです.

とりあえず今回の実装では蟻コロニー最適化のほうが早かったですが,遺伝的アルゴリズムも突然変異をもう少し工夫してやればもっと良い結果が出るような気がします.いつか焼きなまし法とかやってみたいです.

2016年2月10日水曜日

離散コサイン変換(DCT)をPythonで実装した[Python]

ダブルミーニングでテストが終わりました.春休みに入ったので今年の目標の「数学を勉強する」を早速実行してみます.将来勉強するはめになりそうなものから学ぼうと思うので,まずは主軸変換から勉強します.

フーリエ変換はnumpyやscipyで実装されているので,されていないであろう離散コサイン変換(DCT)を今日は実装してみました……と思っていたらscipy.fftpackにdctの実装ありました(公式サイト).こちらはWikipediaで紹介されている基底関数DCT-I,DCT-II,DCT-IIIが実装されているようです.

1.離散コサイン変換(DCT)の概要

離散コサイン変換(Discrete Cosine Transform)はWikipediaさんを見る限りではいろんな基底関数があるようですが,今手元にある本『画像処理とパターン認識入門』では離散コサイン変換の基底関数を次のように定義しています.Wikipediaで言うところのDCT-IIですね.
\[ \phi_k[i] = \begin{cases} \frac{1}{\sqrt{N}} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(k=0) \\ \sqrt{\frac{2}{N}} \cos {\frac{(2i+1)k\pi}{2N}} (k=1,2,...,N-1) \end{cases} \]

この基底関数ですがこんな感じのグラフになります.N=64,横軸がiに相当し上からk=0,1,2,...,9の時の基底関数のグラフです.この定義だとほぼ普通の三角関数ですね.
離散化された信号をf_iとすると,離散コサイン変換,離散コサイン逆変換の式は以下のとおりです.
\[ C[k] = \sum_{i=0}^{N-1}f_{i}\phi_k[i] \,\,(k=0,1,...,N-1)\\ f_{i}=\sum_{k=0}^{N-1}C[k]\phi_{k}[i] \,\,(i=0,1,...,N-1) \]
とりあえずこいつを実装します.

2.離散コサイン変換(1次元)の実装例

まぁこんな感じです.関数phiは基底関数です.本で紹介されていた数式以外にWikipediaにあったDCT-IVもコメントアウトしてありますが一応実装してあります.
ベンチマークとして乱数で作った配列データを試しに変換・逆変換してみます.

実行結果

まぁ問題なくデータが復元されてますね.乱数なので高周波数成分があるので圧縮は無理ですが.

3.離散コサイン変換(2次元)の概要

さて手元の本のタイトルには「画像処理」の文字が書かれています.すなわち二次元データに対する解説も書かれています(というかむしろそれが本題).なのでpythonで二次元データに対しても実装してみます.

本によれば二次元のDCT係数F_k,lは次のように表されるとのこと.また画像の濃度値f_i,jはDCT係数を用いてこのように逆変換できるとのこと.
\[ F_{k,l}=\sum_{j=0}^{N-1}\sum_{i=0}^{N-1}f_{i,j}\phi_k[i]\phi_l[j] \\ f_{i,j}=\sum_{l=0}^{N-1}\sum_{k=0}^{N-1}F_{k,l}\phi_k[i]\phi_{l}[j] \]
ここでこの基底ベクトル(基底画像)がどんな形をしているか見てみます.
この画像,離散コサイン変換(DCT)で検索すると出てくるやつですね.これらの画像は10x10ピクセルの画像を100個並べたものですが,二次元離散コサイン変換ではこの100個の画像をうまいこと足し合わせることで様々な画像を表現するというのがテンションが上ります.(つまりこの画像の中にあんな画像,こんな画像まで,ありとあらゆる画像が含まれていると言っても過言ではない……?)

4.離散コサイン変換(2次元)の実装例

では実装です.scipyにも2次元の離散コサイン変換は実装されてないのでちょっとテンションが上ります.まぁ2次元のフーリエ変換は実装されてますが.
サンプルとして適当な画像を変換・逆変換してみます.
実行結果です.
画像がちゃんと元の状態に戻っているのがわかります.ちなみにDCT係数(周波数成分)はこんな感じでした.
概ね右下に行くほど色が薄くなっています.ここから低周波成分(左上)が多くを占めていることがわかります.
高周波成分を消していったものがこちらです.
DCT係数cについてc[i:,i:] = 0とした時に復元された画像です.i=3でだいたい認識できそうですね.結構データを圧縮できてるのではないでしょうか.まぁまぁうまく行って良かったです.

2016年2月5日金曜日

相関係数を用いてmnistのラベルを認識してみる[Python]

線形代数はわかるけど,ディープラーニングとかわかる気がしねぇ(゜∀。)ワヒャヒャヒャヒャヒャヒャ
という人向けです.
各ラベルごとの平均画像との相関関係からmnistのラベルを予想します.

1.使用する画像

MNIST は 0 から 9 の手書き数字が書かれた,70,000 点の白黒画像です.こんな感じ.
画像は28×28の784ピクセル.今回は学習データを評価用のデータに含めちゃダメなので,この画像データ70,000点のうち,60,000点を学習用のデータとし,10,000点を評価用のデータとします。
とは言うものの学習という学習はしないのであまり関係がないけれど.

2.相関係数の求め方

まぁ今回は一応相関係数を求める部分も実装したのでTeXの練習がてら書いときます.普通はnumpy.corrcoef関数とか使えばいいと思います.
次のようなベクトルx,yがあったとします.
このとき,x,yの相関係数は次のようになります.

ただし,xバー,yバーは各要素の平均です.
同じことですが,これはx,yの平均ベクトルを次のように定義してやるともう少し見通しが良くなります.
とすると,相関係数rは
と表すことができ,x,yベクトルから要素の平均値を引いたベクトルの方向余弦(cosθ)になるわけです.(割とこれが言いたいだけの記事だったりする.線形代数の先生がドヤ顔して話してた)

今回作ったプログラムではxは調べたい画像のベクトル,yは0,1,2,3...9の平均画像に当たります.最も似ている(相関がある)ものを探していきます.

3.ソースコード


4.結果

こんな感じ.
[[ 878    0    3    6    4   40   36    0   14    2]
 [   0 1095   16    4    1    6    2    3   33    1]
 [  26   37  732   44   23    1   43   24   45    4]
 [   8   14   35  809    2   45    6   17   73   25]
 [   4   14    5    0  818    1   24    4   17  117]
 [  29   41   13  127   21  554   24    9   36   23]
 [  21   23   22    1    9   25  862    0    4    0]
 [  17   27   15    0   17    1    0  823   18   68]
 [   7   44   14   85    8   31   10    8  728   38]
 [  19   19    9   12   66    4    1   53   33  820]]
             precision    recall  f1-score   support

        0.0       0.87      0.89      0.88       983
        1.0       0.83      0.94      0.88      1161
        2.0       0.85      0.75      0.79       979
        3.0       0.74      0.78      0.76      1034
        4.0       0.84      0.81      0.83      1004
        5.0       0.78      0.63      0.70       877
        6.0       0.86      0.89      0.87       967
        7.0       0.87      0.83      0.85       986
        8.0       0.73      0.75      0.74       973
        9.0       0.75      0.79      0.77      1036

avg / total       0.81      0.81      0.81     10000
統計はよく知らない民なのでprecision,recallでググるとwikipediaの画像が非常にわかりやすいので見ると良いかと思います.まぁ普段生活する中で言う「的中率」というのはrecallだと思う….

幸いprecision,recall,f1-scoreも81%なので,とりあえず結果は81%ってことですね.mnistの公式サイトを見ると圧倒的最下位ですね.現状1位はConvolutional neural networkの99.77%ッスか.もうこの域まで来ると人よりも精度がいいらしいのでもう関係ないような気がしますがね.

5.感想とか

こうしてみるとやっぱりニューラルネットすげーってなるわけですが,もう少しニューラルネットを使わずに精度を上げられないかと思われるわけですが,ミスした画像を見ればどうすればいいかはある程度わかります.
一番左は平均画像.周りが判定ミスをした画像たちです.数字が書かれている位置が違っていたり,ある程度傾いていたりすると平均画像との相関係数だけだと限界があるわけですね.
もし更に精度を上げるには傾きや位置によらない特徴量を抽出してやる必要がありそうです.その特徴量ベクトルに対して相関係数を用いるのもよし,ユークリッド距離を用いるのもよし,別の手法でやってみるのも良しです.

ちなみにKaggleにはDigit Recognizerというのがあって,Forumとかを見れば結構わかりやすいコードだとかがありますし,自分で実装してmake a submissionで結果を送信するのも楽しいです.
でも大体みんなForumにある頭のいい人が書いたソースコードを実行して送信してるだけみたいな人が多いのでオリジナルのコードで上位に食い込むにはやっぱりディープラーニングとかを実装するしかない気がします.

6.参考

[1] MNIST handwritten digit database, Yann LeCun, CorinnaCortes and Chris Burges
[2] 酒井幸市 著,『画像処理とパターン認識入門』,森北出版株式会社
[3] 多層パーセプトロンで手書き数字認識 - 人工知能に関する断創録

2016年2月2日火曜日

振り子の動きを画像認識で検出してみる[Python]

前やった特定の色の物体の座標を検出の続きみたいなもんです.
振り子を作ってWebカメラを使って振り子の位置を検出して,グラフにしてみます.

まぁ今日のテストがイミフで,水曜日の第二外国語の勉強したくないだけですけど.

1.振り子

こんなの

文句はなしで行きましょう(笑).
暇つぶしクオリティなので重りはキャップを2つくっつけただけです.
糸も細くて変形しにくい針金などではなく,裁縫道具に入ってたやつです.

2.検出プログラム

こんな感じ.output.csvに時間,x座標,y座標が保存されます.

検出アルゴリズムとしては,HSV色空間にWebカメラの画像を変換し,緑色の成分のみをOpenCVのinRange関数を使って取り出しています.
なので以下のプログラムは「緑色の物体」の座標を認識しているので緑色のものがないところでやるようにしてください.

座標は緑色が一番たくさんあった行と列の要素番号をx,yの座標としてます.
まぁそのへんはプログラムを参照してください.

3.プロット

データがoutput.csvに入ったら次はmatplotlibで可視化してみます.可視化プログラムはこんな感じです.
グラフはこんなんが出てくる.これ横からではなく,下から撮影したのでこうなってます.

x,y座標の振り幅の極大極小はそれなりに取れてますね.振り幅が小さくなるとノイズが目立ちますが.(まぁあとカメラの傾きとか重力の方向とかどうなってんのみたいな話は置いといて……)


今日は以上です.
振り子の長さは大体37cmだったので,T=2*pi*sqrt(l/g)で,周波数はおおよそ0.8[Hz]になっているはずで,本当はここでフーリエ変換で周波数を求めようと思ったのですが……これ一定間隔でサンプリングしてないのです.

numpyやscipyのフーリエ変換の関数はたぶんサンプリング間隔が一定という条件でやってる気がする……時間のデータを与えられないし.
複素数をPythonで扱ったりとかはできず,現状実装できる能力がないので近いうちにフーリエ変換を実装しなければ…….

2016年1月30日土曜日

ラグランジュ補間でルンゲ現象を再現してみた[Python]

テスト期間中だけどラグランジュ補間のプログラムを作成した.で,試しにルンゲ現象を再現してみた.

ベッ……別に単位落とした時の言い訳なんかじゃないんだからねっ(。>﹏<。)

1.前知識

ラグランジュ補間は関数の近似を行うアルゴリズム手法の一つ.例えば値がわかっているデータが10個あったとする(例えばこんな感じ.x_data = { x0,x1,x2,...,x9}の時のy_data={y0,y1,y2,…,y9})

この時,f(x)を近似したい式として,近似式p(x)は次のように計算される.

この時,f(x)を近似したい式として,近似式p(x)は次のように計算される.
N(i,x)はxがi番目のデータに近い数字なら1を返す.もしxがi番目のサンプルデータに近ければ,i番目の時だけN(i,x)が1になるように出来てる.

こいつの利点としては雨の雨量みたいな離散的なデータを積分したいときに,そのサンプリングした奴の間のデータを出せたりする.

2.プログラム

とりあえずプログラムはこんな感じ.例としてsin関数の近似をやってみた.
で,結果はこんな感じ.

6個のデータしか与えてないのにだいぶいい感じに近似出来てる.すごい.まぁサンプルデータが偏ったりしているとうまくいかない.あとルンゲ現象なんてものもある.ちょっとやってみますか.

3.ルンゲ現象

Wikipediaに載ってる式を試しにやってみます.この式です.

プログラムはこいつ.
結果はこんな感じ.

荒ぶってますねぇ~.補間する点を多くすると過学習みたいになってまぁ変な風になってしまうというやつですね.まぁ気をつけましょうってことっすかね.

多変数関数の補間とかできるのか気になる.

2016年1月20日水曜日

MNISTの自己組織化マップ(SOM)を作った[Python]

さっき昔の記事の修正してたら操作を誤ってしまい,古い投稿がさっき投稿されたようになっててちょっとつらいです…….なのでもう少し温めてから出そうと思ってた内容をさっさと書いちゃいます.

自己組織化マップというのは,教師なしクラスタリングの一種で,似た者同士がまとまっていくものです.詳細はプログラムのコメントに書いたリンク先を参照してください.

とりあえずこういった感じの画像が得られます.(実際に使うときはこういった画像を探してくるのではなくて写像した2次元座標を使います.scikit-learnにもクラスが無く,ひょっとしたら可視化以外に使う人はもういないのかもしれません)

MNISTの手書き数字を似た者同士でまとめてる感じです.上のような画像を作るには結構演算が必要なので以下のプログラムのパラメータを幾つか変更してやる必要がありますが,基本的にはこのプログラムで出来ます(n_side=50,n_learn=5000ぐらい).これを見て,画像認識のプログラムを作る際に,どういった特徴量を作るかを考えるのに役立てることができます.

まぁ画像以外にも元素のデータ突っ込んで似た元素でまとめたりとか,動物の特徴から似た者同士をまとめたりとかという使い方もあります.よくネットに上がっている例だと色の自己組織化マップがあります.そんなに時間がかからないので出来る様子を動画にしてみました.

最初は結構変化しますが,時間が立つに連れて落ち着いてくるようすがわかります.

以下ソースコード

PythonでWebカメラの画像を保存する[Python]

 Pythonで画像処理の勉強をちょっとやってみたいと思い,とりあえずWebカメラから画像を取得してみた。

1.OpenCVのインストール


 OpenCVの公式サイトのダウンロードページからOSに合ったバージョンをダウンロードしてインストールする。opencvフォルダ内のbuild\python\2.7\x86\cv2.pydをPython27フォルダ内のLib\site-packagesにコピーすれば大丈夫なはず。

2.プログラム

 俺が書いたみたいに見せてるけれどOpenCVの公式ドキュメントをちょこっと変えただけでなんか申し訳ない。sを押すとphoto.jpgという名前で保存される。qを押すと終了する。

# -*- coding: utf-8 -*-
import numpy as np
import cv2

cap = cv2.VideoCapture(0)

while(True):
    # フレームをキャプチャする
    ret, frame = cap.read()

    # 画面に表示する
    cv2.imshow('frame',frame)

    # キーボード入力待ち
    key = cv2.waitKey(1) & 0xFF

    # qが押された場合は終了する
    if key == ord('q'):
        break
    # sが押された場合は保存する
    if key == ord('s'):
        path = "photo.jpg"
        cv2.imwrite(path,frame)

# キャプチャの後始末と,ウィンドウをすべて消す
cap.release()
cv2.destroyAllWindows()


 左はArduinoの箱。特に意味は無いです。右は画面を撮影したらこうなった。写し鏡みたいな。こちらも特に意味は無いです。

 ちなみにcv2.Cannyという関数を使うとエッジ検出ができます。ラズパイでは難しいかもしれないですが,普通のパソコンならリアルタイムで処理できるはずです。

# -*- coding: utf-8 -*-
import numpy as np
import cv2

cap = cv2.VideoCapture(0)

while(True):
 # フレームをキャプチャする
 ret, frame = cap.read()

 # エッジ検出
 frame = cv2.Canny(frame,100,200)

 # 画面に表示する
 cv2.imshow('frame',frame)

 # キーボード入力待ち
 key = cv2.waitKey(1) & 0xFF

 # qが押された場合は終了する
 if key == ord('q'):
  break
 # sが押された場合は保存する
 if key == ord('s'):
  path = "photo.jpg"
  cv2.imwrite(path,frame)

# When everything done, release the capture
cap.release()
cv2.destroyAllWindows()


 こちらエッジがされた画像。緑茶の二文字が見て取れる。画像処理と言っても既存の関数を使うだけで中身を知らなくてもできるという現実が嬉しいような悲しいような……。