リカレントネットワークによる記憶画像の想起
こんばんは．
今日は課題で作ったリカレントネットワークによる記憶画像の想起スクリプトを載せます．画像ファイルはPBM形式のものを扱ってます．課題用なのでかなり課題に特化したものなので，特にファイル操作部分に関しては各自変更が必要でしょう．
import numpy as np import random as rnd import copy def sgn(num): return (1 if(num > 0.0) else (0 if(num == 0.0) else -1)) class AssociativeMemory: def __init__(self, patternNum=1): self.__patternNum = patternNum self.__pattern = [] #ファイル操作 tmpArray = [] lineNum = 0 for line in open('ファイル名','r'): if lineNum >= 2: line.rstrip('\n') for charNum in xrange(len(line)): if ((line[charNum] == "1") or (line[charNum] == "0")): tmpArray.append(1.0 if (int(line[charNum]) == 1) else -1.0) lineNum = lineNum + 1 self.__size = len(tmpArray) self.__lena = np.array(tmpArray) self.__pattern.append(self.__lena) #直交する記憶パターンの生成 for k in xrange(self.__patternNum): tmpNDArray = np.zeros(self.__size) for i in xrange(self.__size): tmpNDArray[i] = (1.0 if(rnd.random() > 0.5) else -1.0) self.__pattern.append(tmpNDArray) #重みの計算 self.__weight = np.zeros((self.__size, self.__size)) for i in xrange(self.__size): for j in xrange(i, self.__size): if i != j : for p in xrange(len(self.__pattern)): self.__weight[i][j] = self.__weight[i][j] + self.__pattern[p][i] * self.__pattern[p][j] self.__weight[i][j] = self.__weight[i][j] / (self.__size) self.__weight[j][i] = self.__weight[i][j] else: self.__weight[i][j] = 0.0 self.__weight[j][i] = 0.0 def recall(self, noiseRate = 0.0, step = 15): state = copy.deepcopy(self.__lena) self.__m = [] #原画像にノイズを加えてそれを初期状態とする for i in xrange(len(state)): if noiseRate > 0.0: state[i] = (state[i] if(rnd.random() > noiseRate)\ else (1.0 if(state[i] == -1)\ else -1.0)) #初期の誤差を計算 prodResult = np.dot(state, self.__lena.T) self.__m.append(prodResult / (self.__size)) for t in xrange(step): nextState = np.zeros(self.__size) #想起 for i in xrange(self.__size): prodResult = np.dot(self.__weight[i], state.T) nextState[i] = sgn(prodResult) #時刻tにおける誤差 prodResult = np.dot(nextState, self.__lena.T) self.__m.append(prodResult / (self.__size)) state = copy.deepcopy(nextState) return self.__m
ドライバはこんな感じ
#!/opt/local/bin/ python # -*- coding: utf-8 -*- import AssociativeMemory test = AssociativeMemory.AssociativeMemory(2)#パターン数を自然数で print test.recall(0.4)#入力画像として原画像にどれだけノイズを混ぜるかを指定（最高1.0）
こんな感じです．できるだけ計算量を削減させようと工夫したけど，結局，かなり遅いです．重みは対称行列なことを利用したりしたんですが，一番時間かかるのは重みの計算です．見通しを立てて，ループ削減が難しそうならC++でやる方がいいよね．
追記があります．次の記事を参照すると幸せになれます．
課題で作ったものいろいろ
こんばんは．明日提出の課題を解いたけどドキュメント化していないアカウントです．
せっかくいろいろ作ったのでなんかの役に立つかもな，もとい，俺こんなの作ったよ褒めて褒めてっという目的の基にスクリプトのソースコードを晒していこうと思います．
まず，ジニ係数算出スクリプトから．
必要なものはNumPy，ScipPyです．
#!/usr/bin/env python import sys,string import numpy as np inputData = [] #データをプログラムで扱うための前処理 #データにあわせて適宜書き換える for line in open('ファイル名','r'): inputData.append(float(line.rstrip('\r\n'))) inputData.sort() N = len(inputData) x = np.array(inputData) mean = np.mean(x) tmp = [] for i in xrange(N): tmp.append(2 * i - N - 1) val = np.array(tmp) gini = np.dot(val, x) / (N**2 * mean) print "Gini Index : " + str(gini)
簡単ですね．
次に平均場近似を用いた画像修復スクリプト．
import sys, copy, string, math, random q = 0.2 h = 0.5 * math.log((1-q)/q) J = 0.4 th = 0.0001 currentError = float(sys.maxint) reconstructRate = 0 errorList = [] baseBitString = [] degBitString = [] newDeg = [] oldDeg = [] baseDeg = [] base = [] #データをプログラムで扱うための前処理 #データにあわせて適宜書き換える lineNum = 1 for line in open('ファイル名','r'): if lineNum >= 3: line.rstrip('\n') for charNum in xrange(len(line)): if ((line[charNum] == "1") or (line[charNum] == "0")): baseBitString.append(1 if (int(line[charNum]) == 1) else -1) if (random.random() < q): degBitString.append((-1 if (int(line[charNum]) == 1) else 1)) else : degBitString.append(1 if (int(line[charNum]) == 1) else -1) lineNum = lineNum + 1 for i in xrange(50, len(baseBitString) - 50, 50): base.append(baseBitString[i:i+50]) baseDeg.append(degBitString[i:i+50]) for i in xrange(len(baseDeg)): tmp = [] for j in xrange(len(baseDeg[i])): newData = math.tanh(J * ((baseDeg[i][j+1] if ((j + 1) < len(baseDeg[i])) else 0)\ + (baseDeg[i][j-1] if ((j - 1) >= 0) else 0)\ + (baseDeg[i+1][j] if ((i + 1) < len(baseDeg)) else 0)\ + (baseDeg[i-1][j] if ((i - 1) >= 0) else 0)) + h * baseDeg[i][j]) tmp.append(newData) oldDeg.append(tmp) print 'Base Image:' for i in xrange(len(base)): print str(base[i]) print 'Degeneration Image :' for i in xrange(len(baseDeg)): print str(baseDeg[i]) print "Error :" while (currentError > th): currentError = 0.0; for i in xrange(len(baseDeg)): tmp = [] for j in xrange(len(baseDeg[i])): newData = math.tanh(J * ((oldDeg[i][j+1] if ((j + 1) < len(oldDeg[i])) else 0)\ + (oldDeg[i][j-1] if ((j - 1) >= 0) else 0)\ + (oldDeg[i+1][j] if ((i + 1) < len(oldDeg)) else 0)\ + (oldDeg[i-1][j] if ((i - 1) >= 0) else 0)) + h*baseDeg[i][j]) currentError = currentError + abs(newData - oldDeg[i][j]) tmp.append(newData) newDeg.append(tmp) oldDeg = copy.deepcopy(newDeg) errorList.append(currentError) print currentError del newDeg[:] for i in xrange(len(oldDeg)): print str(oldDeg[i]) for i in xrange(len(oldDeg)): for j in xrange(len(oldDeg[i])): oldDeg[i][j] = (1 if (oldDeg[i][j] > 0) else (0 if (oldDeg[i][j] == 0) else -1)) reconstructRate = reconstructRate + (1 if (abs(base[i][j] - oldDeg[i][j]) == 0) else 0) print 'Reconstruct Image:' for i in xrange(len(oldDeg)): print str(oldDeg[i]) if (float(len(base) * len(base[0])) > 0.0): print 'Reconstruct Rate : \n \t' + str(reconstructRate) + '/' + str(float(len(base) * len(base[0]))) + \ ' = ' + str(reconstructRate/float(len(base) * len(base[0])))
前処理や出力が多いだけで，操作だけ見るとかなり簡単です．
とにかく，Pythonは
- 簡単に
- 簡潔に
- わかりやすく
プログラムが書ける素晴らしい言語です．ちょっとした課題ならぜひ使うべき．簡潔とわかりやすいって同じじゃね？って思ったあなた，大事なことなので二回言ったんです．
僕なんてちょっとした四則演算のためにインタプリタ起動しちゃうぐらいです（ぇ
次に，Kalmanフィルタをクラスとして実装したものです．
これはC++で書きました．なんでC＋＋かというと，お察しください．
Boostライブラリが必要です．あと，yanoの日記で公開されている逆行列や行列式を計算するコードも利用しました．勝手に拝借してごめんなさい．そしてありがとうございます．
デフォルトでboost::numeric::ublasにはないってことを今回はじめて知ったよ．
Wikipediaのカルマンフィルタの項の設定例が解けるってことは確認しましたが，Kalmanフィルタの理解が充分ではないのでそれ以外には直接使えないクソグラムかもしれません．
#pragma once #include <cstdlib> #include <ctime> #include <cmath> #include <iostream> #include <string> #include <vector> #include <boost/random.hpp> #include <boost/numeric/ublas/io.hpp> #include <boost/numeric/ublas/matrix.hpp> #include "ublasdeterminant.hpp" #include "ublasinvert.hpp" using namespace std; using namespace boost; using namespace boost::numeric::ublas; typedef matrix<double> dmat; class KalmanFilter{ private: int accSigma, mesureSigma; double R; dmat *A,*B,*Q,*H; dmat KalmanGain; public: KalmanFilter(dmat*,dmat*,dmat*, int const); KalmanFilter(dmat*,dmat*,dmat*,dmat*, double const, int const, int const); ~KalmanFilter(); std::vector<dmat> Filtering(dmat*, dmat*, dmat*, dmat*, double const); std::vector<dmat> Filtering(dmat*, dmat*, dmat*, dmat*); std::vector<dmat> Correct(dmat*, dmat*, dmat*, dmat*); std::vector<dmat> Predict(dmat*, dmat*, double const); double calcBelief(dmat*, dmat*, dmat*); }; KalmanFilter::KalmanFilter(dmat* A, dmat* B, dmat* Q, int const accSigma){ this->A = A; this->B = B; this->Q = Q; this->accSigma = accSigma; } KalmanFilter::KalmanFilter(dmat* A, dmat* B, dmat* Q, dmat* H, double const R, int const accSigma, int const mesureSigma){ this->A = A; this->B = B; this->Q = Q; this->H = H; this->R = R; this->accSigma = accSigma; this->mesureSigma = mesureSigma; } KalmanFilter::~KalmanFilter(){ } std::vector<dmat> KalmanFilter::Filtering(dmat* State, dmat* Cov, dmat* Ctrl, dmat* G, double const Mesure){ std::vector<dmat> CorrectResult = KalmanFilter::Correct(State, Cov, Ctrl, G); std::vector<dmat> PredictResult = KalmanFilter::Predict(&(CorrectResult.at(0)), &(CorrectResult.at(1)), Mesure); return PredictResult; } std::vector<dmat> KalmanFilter::Filtering(dmat* State, dmat* Cov, dmat* Ctrl, dmat* G){ mt19937 mesureGen(static_cast<unsigned long>(time(0))); normal_distribution<> mesureDst(0, mesureSigma); variate_generator<mt19937&, normal_distribution<> > GaussNoise(mesureGen, mesureDst); double const Mesure = (prod(*H, *State))(0,0) + GaussNoise(); return KalmanFilter::Filtering(State, Cov, Ctrl, G, Mesure); } std::vector<dmat> KalmanFilter::Correct(dmat* State, dmat* Cov, dmat* Ctrl, dmat* G){ dmat prioriState(State->size1(), State->size2()); dmat prioriCov(Cov->size1(), Cov->size2()); std::vector<dmat> Result; mt19937 accGen(static_cast<unsigned long>(time(0))); normal_distribution<> accDst(0, accSigma); variate_generator<mt19937&, normal_distribution<> > Acceleration(accGen, accDst); prioriState = prod(*A, *State) + prod(*B, *Ctrl) + (*G) * Acceleration(); prioriCov = prod(*A, *Cov); prioriCov = prod(prioriCov, trans(*A)) + *Q; Result.push_back(prioriState); Result.push_back(prioriCov); return Result; } std::vector<dmat> KalmanFilter::Predict(dmat* prioriState, dmat* prioriCov, double const Mesure){ dmat posterioriState(prioriState->size1(), prioriState->size2()); dmat posterioriCov(prioriCov->size1(), prioriCov->size2()); std::vector<dmat> Result; KalmanGain = prod(*H, *prioriCov); KalmanGain = prod(*prioriCov, trans(*H)) / (prod(KalmanGain, trans(*H))(0,0) + R); posterioriState = *prioriState + KalmanGain * (Mesure - (prod(*H, *prioriState))(0,0)); posterioriCov = identity_matrix<double>(posterioriCov.size1()) - prod(KalmanGain, *H); posterioriCov = prod(posterioriCov, *prioriCov); Result.push_back(posterioriState); Result.push_back(posterioriCov); return Result; } double KalmanFilter::calcBelief(dmat* Info, dmat* State, dmat* Cov){ dmat invertCov; ublasinvert::invert(*Cov, invertCov); dmat Gap = prod(trans(*Info - *State), invertCov); return 1 / sqrt(ublasdeterminant::determinant(2.0 * M_PI * (*Cov))) * exp(-0.5 * prod(Gap, (*Info - *State))(0,0)); }
ついでに，ドライバも．
#include <cstdlib> #include <cmath> #include <iostream> #include <string> #include <vector> #include <boost/numeric/ublas/io.hpp> #include <boost/numeric/ublas/matrix.hpp> #include "KalmanFilter.hpp" int main(void); int main (){ int const accSigma = 1; int const mesureSigma = 10; double const R = 10.0; double const deltaT = 1.0; double const Mesure = 5.0; dmat A(2,2), B(2,2), G(2,1), Q(2,2), H(1,2); dmat State(2,1), Cov(2,2), Ctrl(2,1); std::vector<dmat> Result; A(0,0) = 1;A(0,1) = deltaT; A(1,0) = 0;A(1,1) = 1; B(0,0) = 0;B(0,1) = 0; B(1,0) = 0;B(1,1) = 0; G(0,0) = deltaT * deltaT * 0.5; G(1,0) = deltaT; Q = accSigma * prod(G, trans(G)); H(0,0) = 1;H(0,1) = 0; State(0,0) = 0; State(1,0) = 0; Cov(0,0) = 1;Cov(0,1) = 0; Cov(1,0) = 0;Cov(1,1) = 1; Ctrl(0,0) = 0; Ctrl(1,0) = 0; KalmanFilter kf = KalmanFilter(&A, &B, &Q, &H, R, accSigma, mesureSigma); cout << "Initial State : " << State << "\nInitial Cov : " << Cov << endl; for(int i = 0; i < 5; ++i){ if(i < 4){ Result = kf.Filtering(&State, &Cov, &Ctrl, &G); }else{ Result = kf.Filtering(&State, &Cov, &Ctrl, &G, Mesure); } State = Result.at(0); Cov = Result.at(1); cout << "Posteriori State : " << State << "\nPosteriori Cov : " << Cov << endl; } return 0; }
予測だけ使うとかもできます．
#include <cstdlib> #include <cmath> #include <iostream> #include <fstream> #include <string> #include <vector> #include <boost/lexical_cast.hpp> #include <boost/numeric/ublas/io.hpp> #include <boost/numeric/ublas/matrix.hpp> #include "KalmanFilter.hpp" int main(int, char**); int main (int argc, char **argv){ int const accSigma = 1; int const t = (argc == 2)?(lexical_cast<int>(argv[1])):(4); double const deltaT = 1.0; dmat A(2,2), B(2,2), G(2,1), Q(2,2); dmat State(2,1), Cov(2,2), Ctrl(2,1); std::vector<dmat> Result; A(0,0) = 1;A(0,1) = deltaT; A(1,0) = 0;A(1,1) = 1; B(0,0) = 0;B(0,1) = 0; B(1,0) = 0;B(1,1) = 0; G(0,0) = deltaT * deltaT * 0.5; G(1,0) = deltaT; Q = accSigma * prod(G, trans(G)); State(0,0) = 0; State(1,0) = 0; Cov(0,0) = 1;Cov(0,1) = 0; Cov(1,0) = 0;Cov(1,1) = 1; Ctrl(0,0) = 0; Ctrl(1,0) = 0; KalmanFilter kf = KalmanFilter(&A, &B, &Q, accSigma); for(int i = 0; i < t; ++i){ string fileName = "KFCorrect(k=" + lexical_cast<string>(i + 1) + ").txt"; ofstream ofs(fileName.c_str()); Result = kf.Correct(&State, &Cov, &Ctrl, &G); State = Result.at(0); Cov = Result.at(1); for(double location = -20.0; location <= 20.0; location += 0.5){ for(double vel = 0.0; vel <= 0.0; vel += 0.5){ dmat Info(2,1); Info(0,0) = location; Info(1,0) = vel; double prob = kf.calcBelief(&Info, &State, &Cov); ofs << location << "\t" << vel << "\t" << prob << endl; } } } return 0; }
といった感じです．
とにかくC++は
- 冗長
- 冗長
- 冗長
ってことを痛感しましたね（
いいところもありますよ．低級に近い高級言語（shared_pointerとかC++11とかの機能使えばかなり高級かも・・・）なんでとにかく計算が早い．本当に早い．Pythonには早さが足りない！（PyPyあるけどNumPyPyがまだ使い辛い）．
あと，今更だけどboost::numeric::ublasでは
ってことを知りました．余因子とかあの辺りの知識があれば自分で実装できるのかもだけど，素晴らしいコードが既に公開されていたのでそっち使っちゃった///
時間が許せば，せっかくなんでもう2つぐらい課題とにらめっこして作りたいけど・・・まずは今出来ているところまでで提出できる形にしないとね．
Global Game Jam 2012 札幌に参加しました
Global Game Jam 2012（以下，GGJ2012）に札幌会場で参加しました．
そもそも，GGJ2012とは何ぞやというのは公式サイトをご覧ください（英語です）．
http://globalgamejam.org/
簡単に言うと，当日に発表されたテーマを元に，三日間自由にゲーム作って公開しようってお祭りです．
札幌での開催は三回目だかそこらぐらいらしい．運営自体は開催を望む団体自らが行わなければならない．なので，どこでも開催できるけど中々難しいみたい．運営に携わった札幌ゲーム制作者コミュニティkawazの方々，本当にご苦労様でした．
今年のテーマはウロボロスの絵 だ け でした．
僕の属していたグループは，iPadで遊ぶゲームを作ると言うコンセプトのもとに企画を考えて行きました．
ウロボロスの絵から，「ループ」，「間抜け」，「スタイリッシュ」，「ひらめき」，「夢」，「神」などなどイメージを膨らませたり，グラフィッカ担当の方のポートフォリオの世界観を生かそうと模索したり，紆余曲折を経たのちに生まれた僕らのグループの企画は，
「20XX年，かわいいけどなんか殴りたくなるウーパールーパー星人による地球侵略を阻止する新感覚スタイリッシュリズムゲーム Loopy Looper」
というものでした（たしかこんな感じ）．
Objective-Cを用いたiPadアプリケーションの開発ならびにゲームプログラミングの経験のない僕は足を引っ張らないよう全力（ときに脱力）しながら開発を進めました．
後輩であり先輩である@giginetくんがフル稼働してくれたおかげで，無事形にすることができたのは本当によかったです．
とりあえず，
- Objective-Cマジすごい！
- Objective-Cマジ冗長！
- Objective-Cマジ@とか[]とか記号使いすぎ！キモイ！
- Objective-Cマジ信者力必須！
という感想です（ぇ
あと，Cocos2dとKobold2dっていうiOS向けゲーム開発用フレームワークが秀逸過ぎて半端無かったです．マジCocos2dとKobold2d最強．Cocos2dとKobold2d可愛い！（ぇ
最終日の夜には成果発表会がありました．
他のグループのゲームもマジでクオリティ高くて半端なかったス．
Unityが今年は熱いみたいですが，この短期間でアレだけの素晴らしいビジュアルのゲームが作れるとかキモすぎです．
絵師さんも音楽屋さんもこの三日間という短期間で素晴らしい質の物を生産してて，そういう人と一緒に共通の目的持って作業できるGGJパないっす．
発表会後，お片付けしたあとヨシミキッチン 札幌PARCO店で懇親会．
多くの参加者，観覧者の方々とゲームに関する事やそうじゃないこと，いろいろ話せて楽しかったです．
ヨシミキッチン 札幌PARCO店は一昨年前のOSC札幌以来来てなかったのですが，やはりご飯もお酒も美味しかったです．この会費でこれだけ美味しいものを食べれるなんて素晴らしいですね．みなさん，どんどん懇親会とかで利用しましょう．
乱文ではありましたが，僕のGGJはこんな感じでした．
あ！三月に札幌C++勉強会が催されるそうなので，ぜひ興味がみじんでもある人は参加しましょう．
都合が合えば僕も見に行きたい！
Clojure始めました
ClojureはJVM上で動くLISP方言です．関数型言語を触ろうと思い立ち，いろいろ足りない頭で考えた結果，Clojureを触ろうってなった感じです．
いろいろなWebサイトやブログのエントリを見ながら，Clojure入れて，Clojure-contrib入れて，Emacs用にClojure-mode，Swank-Clojure，SLIME入れてとやってました．
しかしSLIMEは上手く動いてくれませんし，REPLも起動前に例外投げて起動しないことも有ります．マジ俺涙目．
想像力を母体に置き去りにしてしまった自分は，自分で問題設定して新しい言語を触ってみるなんて芸当が出来ません．
なので，学校で出たプログラミングを伴う課題をClojureで解いてみました．
;;Simple Perceptron Program ;;2012/01/13 (use 'clojure.contrib.math) (defn staircase-function [var] (if (> var 0) 1 0)) (defn sgn [var] (if (> var 0) 1 (if (= var 0) 0 -1))) (defn learning-theta [lambda sig_t sig_s] (* lambda (* sig_t (* (staircase-function (* (- sig_t) sig_s)))))) (defn main [th_ast l0] (loop [t 0 th (rand)] (let [x (rand) l l0] (print t "\t" th "\t" l "\t" (expt (- th_ast th) 2) "\n") (if (and (< (- th_ast th) 0.001) (> (- th_ast th) -0.0001)) th (recur (inc t) (+ th (learning-theta l (sgn (- th_ast x)) (sgn (- th x))))))))) (def theta 0.5) (def lambda-zero 0.0001) (main theta lambda-zero)
これと
(use 'clojure.contrib.math) (defn staircase-function [var] (if (> var 0) 1 0)) (defn sgn [var] (if (> var 0) 1 (if (= var 0) 0 -1))) (defn attenuation-function [l0 a t] (* l0 (. Math (exp (* (- a) t))))) (defn learning-theta [lambda sig_t sig_s] (* lambda (* sig_t (* (staircase-function (* (- sig_t) sig_s)))))) (defn main [th_ast l0] (loop [t 0 th (rand)] (let [x (rand) l (attenuation-function l0 0.001 t)] (print t "\t" th "\t" l "\t" (expt (- th_ast th) 2) "\n") (if (and (< (- th_ast th) 0.001) (> (- th_ast th) -0.001)) th (recur (inc t) (+ th (learning-theta l (sgn (- th_ast x)) (sgn (- th x))))))))) (def theta 0.5) (def lambda-zero 0.99) (main theta lambda-zero)
これです．
とてもシンプルなパーセプトロンです．前者は学習率係数が固定で後者が減衰型です．
慣れてしまえば結構面白いものだなー．
第2回札幌C++勉強会に参加
タイトルのように，第2回札幌C++勉強会に参加してきました．
3つのセッションと3つのライトニング・トークという内容でした．
間に休憩を挟みつつ，13:15から16:55まで行われました．
最初のセッションでは，必読書や要フォローアカウント（Twitter）の紹介など，初心者C++erのためのセッションという感じでした．
とても魅力的な書物や人物たち（闇の軍団と呼ばれているらしい・・・）を紹介していただけて，初心者C++erとしてはありがたい限りでした．
あ，早速紹介された本を買いました．じっくり読み進めたいと思います．
発表者はlapis_twさんでした．
二つ目のセッションは例外安全入門というものでした．
C++プログラムでの例外に対する対処について，基礎的な事項や技術についてわかりやすく説明していただきました．
C++のコードを基に技術に関する説明は行われたので，実際に自分が使うときにどのように適用すればよいのか想像しやすかったです．
発表者はhotwatermorningさんでした．
最後のセッションは，h_hiro_さんによるC++マクロを使い倒すというお話でした．
時間が押し気味だったため，充分な発表時間を確保できなかったようでした．
ですが，「マクロでここまでできるのか！」という驚きは充分に伝わりました．
LTでは，ignis_fatuusさん，nakayoshixさん，tututenさんが5分間C++に関係あることから無いことまでを喋っていただきました．
それぞれ，Boost.Phoenix，Python札幌とクラウド温泉，コードゴルフ（ショートコーディング？）という内容でした．
ざっくばらんに書いたらこんな感じです．
懇親会・・・行けばよかったかなぁ．
次回勉強会は未定ですが，個人的にPython札幌に出てみたい．
追記：
LT発表者のところ，tututenさんに訂正しました．失礼いたしました．
ProcessingでLotka-Volterraの方程式の実験
学校の課題で作ってみました．クソグラムです．
//################################################ //Lotka-Volterra equation Program 2011/05/31 //Programmer : Takuma KATANOSAKA // Enviroment : 2D 100 x 100 cell (torus) //################################################ import java.util.*; import java.lang.*; //################################################ //Global Variables // predator // true (If the predator exists in (X,Y)) // false (otherwise) // // prey // true (If the prey exists in (X,Y)) // false (otherwise) // // Vector<Predator> predator_vect //################################################ boolean[] predator = new boolean[10000]; boolean[] prey = new boolean[10000]; Vector<Predator> predator_vect = new Vector(); //################################################ //Global Methods // stocasticChoise(float []) // stocasticChoise(float [], Vector<Integer>) // renewCoord(int [], int) //################################################ int stocasticChoise(float [] reward_table){ float sum_gi = 0; Vector oi = new Vector(); for(int i = 0, imax = reward_table.length; i < imax; ++i){ sum_gi += pow(2.0, reward_table[i]); } for(int i = 0, imax = reward_table.length; i < imax; ++i){ oi.add(pow(2.0, reward_table[i]) / sum_gi); } return (int)oi.indexOf(Collections.max(oi)) + 1; } int stocasticChoise(float [] reward_table, Vector<Integer> removeElementsIndex){ float sum_gi = 0; Vector oi = new Vector(); for(int i = 0, imax = reward_table.length; i < imax; ++i){ sum_gi += pow(2.0, reward_table[i]); } for(int i = 0, imax = reward_table.length; i < imax; ++i){ oi.add(pow(2.0, reward_table[i]) / sum_gi); } for(int i : removeElementsIndex){ oi.remove(i); } return (int)oi.indexOf(Collections.max(oi)) + 1; } void renewCoord(int [] coord, int mv){ switch(mv){ case 1: coord[1] -= 1; break; case 2: coord[0] -= 1; coord[1] += 1; break; case 3: coord[0] += 1; coord[1] += 1; break; case 4: coord[0] -= 1; coord[1] -= 1; break; case 5: coord[0] += 1; coord[1] -= 1; break; default: exit(); break; } //Enviroment is torus. if(coord[0] < 0) coord[0] = 99; else if(coord[0] > 99) coord[0] = 0; if(coord[1] < 0) coord[1] = 99; else if(coord[1] > 99) coord[1] = 0; } //################################################ //Predator // -Output : Front, Left, Right, Hardly Left, Hardly Right // -Output Function : Stochastic Choice // -States : Ability Table // -State Transition Function : E = E + delta E (To move neighbor cell need 5 energy) // delta E = 15, Ability(t + 1) = Ability(t) + 1 (when it collides with a prey) // delta E = -5, Ability(t + 1) = Ability(t) - 1 (otherwise) // If predotor's energy is bigger than 1000, its predator segments with mutation. //################################################ class Predator{ private int x, y; private float [] reward_table = {10,10,10,10,10}; private int E; Predator(){ x = (int)random(0, 99); y = (int)random(0, 99); E = 500; } Predator(Predator p){ x = p.getX(); y = p.getY(); E = p.getE() / 2; //Mutation if(p.getE() >= 1000){ float random_value = random(0, 1); //reward_table = Arrays.copyOf(p.getRewardTable(), p.getRewardTable().length); float [] temp = p.getRewardTable(); for(int i = 0, imax = p.getRewardTable().length; i < imax; ++i){ reward_table[i] = temp[i]; } if(0 <= random_value && random_value < 0.2){ reward_table[0] += reward_table[0] * random(-1, 1); }else if(0.2 <= random_value && random_value < 0.4){ reward_table[1] += reward_table[1] * random(-1, 1); }else if(0.4 <= random_value && random_value < 0.6){ reward_table[2] += reward_table[2] * random(-1, 1); }else if(0.6 <= random_value && random_value < 0.8){ reward_table[3] += reward_table[3] * random(-1, 1); }else{ reward_table[4] += reward_table[4] * random(-1, 1); } } } public int getX(){ return x; } public int getY(){ return y; } public int getE(){ return E; } public float [] getRewardTable(){ return reward_table; } public void move(){ int mv = stocasticChoise(reward_table); int [] temp_coord = {x, y}; Vector<Integer> removeElementsIndex = new Vector(); /** Forward : 1 Left : 2 Right : 3 Hardly Left : 4 Hardly Right : 5 */ predator[100*x + y] = false; renewCoord(temp_coord, mv); while(predator[100*temp_coord[0] + temp_coord[1]]){ removeElementsIndex.add(mv - 1); if(removeElementsIndex.size() == 5){ temp_coord[0] = x; temp_coord[1] = y; break; } mv = stocasticChoise(reward_table, removeElementsIndex); renewCoord(temp_coord, mv); } x = temp_coord[0]; y = temp_coord[1]; predator[100*x + y] = true; //Eval Energy energy(mv); } private void energy(int action){ if(prey[x*100 + y]){ reward_table[action - 1] += 0.5; prey[x*100 + y] = false; E += 15; }else{ reward_table[action - 1] -= 3.0; E -= 5; } } public Vector segmentation(){ Vector son_vect = new Vector(); if(E >= 1000){ for(int i = 0; i < 2; ++i){ Predator son = new Predator(this); son_vect.add(son); } } return son_vect; } public boolean isDead(){ if(E <= 0){ return true; }else{ return false; } } } //################################################ // //Main Method Part (Set up & Draw) // //################################################ void setup(){ frameRate(120); size(600, 600); background(0,0,128); smooth(); //PREY fill(255,255,255); for(int i = 0; i < 100; ++i){ for(int j = 0; j < 100; ++j){ if(random(0, 1) < 0.15){ prey[i*100 + j] = true; //predator[i*100 + j] = true; ellipse(6*i+3, 6*j+3 ,6,6); }else{ prey[i*100 + j] = false; } } } //PREDATOR fill(255,0,0); for(int i = 0; i < 5; ++i){ Predator p = new Predator(); if(predator_vect.add(p) == false){ exit(); } predator[p.getX()*100 + p.getY()] = true; prey[p.getX()*100 + p.getY()] = false; ellipse(6*p.getX() + 3, 6*p.getY() + 3, 6, 6); } } //++++++++++++++++++++++++++++++++++++++++ void draw(){ int prey_num = 0; background(0,0,128); //PREDATOR for(int i = 0; i < predator_vect.size(); ++i){ Predator p = predator_vect.get(i); Vector son_vect; p.move(); if(p.isDead()){ predator[p.getX()*100 + p.getY()] = false; predator_vect.remove(i); }else{ //If a predator is at the point of death, change the color of it. if(p.getE() <= 100){ fill(128, 255, 0); }else{ fill(255, 0, 0); } son_vect = p.segmentation(); if(son_vect.isEmpty()){ ellipse(6*p.getX() + 3, 6*p.getY() + 3,6,6); }else{ predator_vect.remove(i); if(!predator_vect.addAll(son_vect)){ exit(); } fill(255, 255, 0); ellipse(6*p.getX() + 3, 6*p.getY() + 3,6,6); } } } //PREY fill(255,255,255); for(int i = 0; i < 100; ++i){ for(int j = 0; j < 100; ++j){ if(prey[i*100 + j] == true){ prey_num++; ellipse(6*i+3, 6*j+3,6,6); if(random(0, 1) < 0.15){ if((j > 0) && prey[i*100 + (j - 1)] == false && predator[i*100 + (j - 1)] == false){ prey[i*100 + (j - 1)] = true; prey_num++; ellipse(6*i+3, 6*(j - 1)+3,6,6); }else if((j < 99) && prey[i*100 + (j + 1)] == false && predator[i*100 + (j + 1)] == false){ prey[i*100 + (j + 1)] = true; prey_num++; ellipse(6*i+3, 6*(j + 1)+3,6,6); }else if((i > 0) && prey[(i - 1)*100 + j] == false && predator[(i - 1)*100 + j] == false){ prey[(i - 1)*100 + j] = true; prey_num++; ellipse(6*(i - 1)+3, 6*(j + 1)+3,6,6); }else if((i < 99) && prey[(i + 1)*100 + j] == false && predator[(i + 1)*100 + j] == false){ prey[(i + 1)*100 + j] = true; prey_num++; ellipse(6*(i + 1)+3, 6*j+3,6,6); } } } } } println(predator_vect.size() + "\t" + prey_num); }
Processing（Proce55ing : P5）はJavaで簡単にアニメーションを作ることのできる面白い代物．
自分の環境（Mac OSX 10.6）では動作が遅く不安定だけど，Windowsで使うとなかなかサクサク動く．
Javaを知ってて，アニメーションを作りたい人は知ってて損はないものだと思います．
本当に簡単に作れるのでプログラミング教育用にも向いていると思います．
・リンク
公式には英語だけど本体とリファレンスが置いてある．Mac（OS X SnowLeopard）でIT++を使う
独立成分分析（Independent Component Analysis : ICA）の勉強をして，いざ自らの手で実装しようと思ったが，「待てよ・・・C++用のICAを関数として含むライブラリがあるんじゃないか？いや，あるだろう．」ということに実装前に気が付いた．実に気付くのが遅かった．
で，調べてみると本当にあった．IT++（アイティープラスプラス）というものだ．
インストール
・・・の前に，自分の環境を提示しておこう．OS : Mac OS X 10.6.4
プロセッサ：Intel Core 2 Duo 2[GHz]
メモリ：DDR-3 2GB
コンパイラ：i686-apple-darwin10-g++-4.2.1
configure実行前
残念な事に，FinkやMac Portsでは提供されていないので，手動でインストールしなければならない．インストールをする前に，必ずRead MeとInstallを読んだ方が良い．
これが今回の俺の得た一番の教訓である（
それらによると，Fortranの処理系（g77，又は最新のもの），ATLAS（Macの場合），FFTW（最新のもの）が必要だそうだ．
自分は今回，Fortranの処理系はg95にした（Mac Portsではそれしか手に入らなそうだから）．
それ以外のものも全てMac Portsで手に入る．マジでMac Ports便利．
configureの実行からmakeまで
それぞれ必要なものをインストールし終えたら，makeをする前にconfigureを回してmakeするための下準備をしなければならない．僕がconfigureを実行したとき，どうしてもFFTWを探し出してくれなかったので，オプションで以下のように指定した．
./configure --prefix="/opt/local" --with-fft=fftw3 --with-fft-include="/opt/local/include"ちなみに，prefixは，インストールするディレクトリを指定するためのオプションである．
こうすると，見事にconfigureは正常に終了して，makeするための準備が整う．
あとは，make installするだけ．
その他インストールに関すること
とりあえず，makeが失敗したらmake cleanとmake distcleanすればいい．atlasとか，その他の外部ライブラリが見つからない場合にも，オプションが用意されているのでRead MeやInstallを読むこと．
何が問題かを把握する為に，吐き出されたメッセージは念のために一からよく読んだ方が良い．
手間を惜しむと時間が無駄に過ぎる．
IT++を実際に使ってみる
無事，インストールも終了したので，早速，自分のプログラムに組み込んでみた．ICAを使う上で必要なヘッダーファイルはitpp/signal/fastica.hです．
fastica.hについてのIT++のページのExampleに従ってプログラムを作ると・・・エラーが出ます．
fastica.hをエディタで開くと分かりますが，クラス名はFastICAではなく，Fast_ICAです．
そんなこんなで作ったプログラムをコンパイルするとき，単純に
g++ -o [実行ファイル名] [ソースファイル名]ではUndefined Symbolとエラーが返されます（俺だけ？）．
なので，-Lオプションでライブラリのディレクトリを指定し，-l（小文字のエル）で使用するライブラリ名を指定してあげなければなりません．
つまり，以下のような感じになると思います．
g++ -L/opt/local/lib -litpp -o [実行ファイル名] [ソースファイル名]これで，コンパイルも通り，実行ファイルもできるはずです．
