tag:blogger.com,1999:blog-47664147338902038662024-03-14T18:02:22.717+09:00光学設計者の学習メモ最近勉強したことをまとめて行きたい。
画像やカメラよりの勉強が多いかも。anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.comBlogger144125tag:blogger.com,1999:blog-4766414733890203866.post-76611386847122706532016-04-15T21:27:00.000+09:002016-04-15T21:27:06.580+09:00Compressed Sensingを理解するために。<div dir="ltr" style="text-align: left;" trbidi="on">
<br />
<div>
<br /></div>
<div>
最近世の中ではCompressed Sensing(CS)が流行っている。</div>
<div>
<br /></div>
<div>
CSはサンプリング手法を大きく変える。</div>
<div>
例えば1ピクセルしかなくても画像を取得できる。</div>
<div>
<a href="http://dsp.rice.edu/cscamera">http://dsp.rice.edu/cscamera</a></div>
<div>
この技術のいいところは、検出器が高価でピクセル配列を準備できないものなんかでも、空間情報を取得できる点だ。(X線とか?)</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
大学で習うサンプリングといえば<a href="https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem" target="_blank">Nyquist-Shannon Sampling Theorem</a>で、これは取得したい周波数の2倍の周波数でサンプリングすれば信号を完全に復元できるというもの。</div>
<div>
CSはもっと少ないサンプルでも信号を復元できるよ、というものらしい。</div>
<div>
<br /></div>
<div>
例えばカメラはRAWからJPEGに変換する際、画像は多くの情報を捨てている。DCTによって重要な係数のみを保存する。あまり重要でないものは捨ても問題ない。むしろ、保存や転送には容量が小さいほうがよい。</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
なんかものすごく将来性のありそうな技術で勉強を始めたけど、数学が難しいorz</div>
<div>
自分のためにCSに学んだことをアップしていこうと思う。</div>
<div>
<br /></div>
<div>
以下目次</div>
<div>
<br /></div>
<div>
<ul style="text-align: left;">
<li>Restricted Isometry Property</li>
<li>Matching Pursuit</li>
<li>Basis Pursuit</li>
<li>Soft Thresholding</li>
<li>など</li>
</ul>
</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-83668255086360140462015-11-03T01:43:00.001+09:002015-11-03T02:00:22.708+09:00RISEでインタラクティブなIPython notebookのスライド作成!<div dir="ltr" style="text-align: left;" trbidi="on">
<div class="markdown-here-wrapper" data-md-url="https://www.blogger.com/blogger.g?blogID=4766414733890203866#editor/target=post;postID=8366825508636014046;onPublishedMenu=posts;onClosedMenu=posts;postNum=0;src=postname">
<div style="margin: 0px 0px 1.2em !important;">
前にも書いたけど、ちょっとしたプレゼン資料を作るときに以下のようなことが多々ある。</div>
<div style="margin: 0px 0px 1.2em !important;">
手軽に書きたい<br />
Texで数式<br />
図を入れたい。<br />
プログラムを入れたい。<br />
かっこ良くしたい。<br />
PDFにしたい。<br />
個人で使う分はフリーで済ませたい。 </div>
<div style="margin: 0px 0px 1.2em !important;">
しかもreveal.jsを使っていかにもギークっぽくしたい。<br />
そんな要求に答えてくれるライブラリがある!<br />
<a href="https://github.com/damianavila/RISE">https://github.com/damianavila/RISE</a></div>
<h2 id="-ipython-notebook-" style="border-bottom-color: rgb(238, 238, 238); border-bottom-style: solid; border-bottom-width: 1px; font-size: 1.4em; font-weight: bold; margin: 1.3em 0px 1em; padding: 0px;">
まずはIPython notebookをインストール。</h2>
<div style="margin: 0px 0px 1.2em !important;">
IPython notebookは<a href="https://www.continuum.io/downloads">Anaconda</a>を入れると一発で入る。Pythonを使う人は入れておいて損はないけど、他の言語でIPython notebook使う人にはちょっと大きすぎるかも。</div>
<h2 id="rise-" style="border-bottom-color: rgb(238, 238, 238); border-bottom-style: solid; border-bottom-width: 1px; font-size: 1.4em; font-weight: bold; margin: 1.3em 0px 1em; padding: 0px;">
RISEをインスール</h2>
<div style="margin: 0px 0px 1.2em !important;">
<a href="https://github.com/damianavila/RISE">https://github.com/damianavila/RISE</a><br />
をどこかに保存して,そのディレクトリに移動した後に</div>
<pre style="font-family: Consolas, Inconsolata, Courier, monospace; font-size: 0.85em; font-size: 1em; line-height: 1.2em; margin: 1.2em 0px;"><code style="background-color: #f8f8f8; border-radius: 3px; border-radius: 3px; border: 1px solid rgb(204, 204, 204); border: 1px solid rgb(234, 234, 234); display: block !important; display: inline; font-family: Consolas, Inconsolata, Courier, monospace; font-size: 0.85em; margin: 0px 0.15em; overflow: auto; padding: 0.5em 0.7em; padding: 0px 0.3em; white-space: pre-wrap; white-space: pre;">python setup.py install
</code></pre>
<div style="margin: 0px 0px 1.2em !important;">
するとIPython notebook上に新しいボタンが現れる。<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgbn2IDhMN6Q-ep0QOF4UmZb6p7By1wYcDTnj-Nc3CDRB06KKw_tBsfG5nCycjbKUIrFFyb2cuKAqksC6_Bhiv0kpAGY027Qs7yvugFrd_AnXS_d5A1n1lEULT01zWpR-RIX4nc8KFxta8/s1600/Untitled1+2015-11-03+01-38-48.png"><img border="0" height="83" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgbn2IDhMN6Q-ep0QOF4UmZb6p7By1wYcDTnj-Nc3CDRB06KKw_tBsfG5nCycjbKUIrFFyb2cuKAqksC6_Bhiv0kpAGY027Qs7yvugFrd_AnXS_d5A1n1lEULT01zWpR-RIX4nc8KFxta8/s320/Untitled1+2015-11-03+01-38-48.png" width="320" /></a></div>
<div style="margin: 0px 0px 1.2em !important;">
ボタンを押すとスライドになって☓を押すと元に戻る。</div>
<div style="margin: 0px 0px 1.2em !important;">
そのままではすべて1枚のスライドになってしまうのでCell ToolbarでSub slideとかに設定すればよい。</div>
<div style="margin: 0px 0px 1.2em !important;">
これから会社ではこれ使って発表しよう。</div>
<div style="font-size: 0em; height: 0; margin: 0; max-height: 0; max-width: 0; overflow: hidden; padding: 0; width: 0;" title="MDH:PHA+PGJyPjwvcD48cD7liY3jgavjgoLmm7jjgYTjgZ/jgZHjganjgIHjgaHjgofjgaPjgajjgZfj
gZ/jg5fjg6zjgrzjg7Pos4fmlpnjgpLkvZzjgovjgajjgY3jgavku6XkuIvjga7jgojjgYbjgarj
gZPjgajjgYzlpJrjgIXjgYLjgovjgII8L3A+PHA+PGJyPjwvcD48b2wgZGF0YS1ibG9nZ2VyLWVz
Y2FwZWQtc3R5bGU9InRleHQtYWxpZ246IGxlZnQ7IiBzdHlsZT0idGV4dC1hbGlnbjogbGVmdDsi
Pgo8bGk+5omL6Lu944Gr5pu444GN44Gf44GEPC9saT4KPGxpPlRleOOBp+aVsOW8jzwvbGk+Cjxs
aT7lm7PjgpLlhaXjgozjgZ/jgYTjgII8L2xpPgo8bGk+44OX44Ot44Kw44Op44Og44KS5YWl44KM
44Gf44GE44CCPC9saT4KPGxpPuOBi+OBo+OBk+iJr+OBj+OBl+OBn+OBhOOAgjwvbGk+CjxsaT5Q
REbjgavjgZfjgZ/jgYTjgII8L2xpPgo8bGk+5YCL5Lq644Gn5L2/44GG5YiG44Gv44OV44Oq44O8
44Gn5riI44G+44Gb44Gf44GE44CCPC9saT4KPC9vbD48ZGl2PjxwPjxicj48L3A+PC9kaXY+PGRp
dj48cD7jgZfjgYvjgoJyZXZlYWwuanPjgpLkvb/jgaPjgabjgYTjgYvjgavjgoLjgq7jg7zjgq/j
gaPjgb3jgY/jgZfjgZ/jgYTjgII8L3A+PHA+44Gd44KT44Gq6KaB5rGC44Gr562U44GI44Gm44GP
44KM44KL44Op44Kk44OW44Op44Oq44GM44GC44KLITwvcD48cD5odHRwczovL2dpdGh1Yi5jb20v
ZGFtaWFuYXZpbGEvUklTRTwvcD48cD48YnI+PC9wPjxwPjxicj48L3A+PGRpdiBkYXRhLWJsb2dn
ZXItZXNjYXBlZC1zdHlsZT0idGV4dC1hbGlnbjogbGVmdDsiIHN0eWxlPSJ0ZXh0LWFsaWduOiBs
ZWZ0OyI+PHA+IyMg44G+44Ga44GvSVB5dGhvbiBub3RlYm9va+OCkuOCpOODs+OCueODiOODvOOD
q+OAgjwvcD48L2Rpdj48cD5JUHl0aG9uIG5vdGVib29r44GvPGEgZGF0YS1ibG9nZ2VyLWVzY2Fw
ZWQtdGFyZ2V0PSJfYmxhbmsiIGhyZWY9Imh0dHBzOi8vd3d3LmNvbnRpbnV1bS5pby9kb3dubG9h
ZHMiPkFuYWNvbmRhPC9hPuOCkuWFpeOCjOOCi+OBqOS4gOeZuuOBp+WFpeOCi+OAglB5dGhvbuOC
kuS9v+OBhuS6uuOBr+WFpeOCjOOBpuOBiuOBhOOBpuaQjeOBr+OBquOBhOOBkeOBqeOAgeS7luOB
ruiogOiqnuOBp0lQeXRob24gbm90ZWJvb2vkvb/jgYbkurrjgavjga/jgaHjgofjgaPjgajlpKfj
gY3jgZnjgY7jgovjgYvjgoLjgII8L3A+PHA+PGJyPjwvcD48ZGl2IGRhdGEtYmxvZ2dlci1lc2Nh
cGVkLXN0eWxlPSJ0ZXh0LWFsaWduOiBsZWZ0OyIgc3R5bGU9InRleHQtYWxpZ246IGxlZnQ7Ij48
cD48Zm9udCBkYXRhLWJsb2dnZXItZXNjYXBlZC1zdHlsZT0iZm9udC13ZWlnaHQ6IG5vcm1hbDsi
IHN0eWxlPSJmb250LXdlaWdodDogbm9ybWFsOyI+IyMgUklTReOCkuOCpOODs+OCueODvOODqzwv
Zm9udD48L3A+PC9kaXY+PHA+aHR0cHM6Ly9naXRodWIuY29tL2RhbWlhbmF2aWxhL1JJU0U8L3A+
PHA+44KS44Gp44GT44GL44Gr5L+d5a2Y44GX44GmLOOBneOBruODh+OCo+ODrOOCr+ODiOODquOB
q+enu+WLleOBl+OBn+W+jOOBqzwvcD48cD5gYGA8L3A+PHA+cHl0aG9uIHNldHVwLnB5IGluc3Rh
bGw8L3A+PHA+YGBgPC9wPjxwPjxicj48L3A+PHA+44GZ44KL44GoSVB5dGhvbiBub3RlYm9va+S4
iuOBq+aWsOOBl+OBhOODnOOCv+ODs+OBjOePvuOCjOOCi+OAgjwvcD48cCBjbGFzcz0ic2VwYXJh
dG9yIiBkYXRhLWJsb2dnZXItZXNjYXBlZC1zdHlsZT0iY2xlYXI6IGJvdGg7IHRleHQtYWxpZ246
IGNlbnRlcjsiIHN0eWxlPSJjbGVhcjogYm90aDsgdGV4dC1hbGlnbjogY2VudGVyOyI+PGEgZGF0
YS1ibG9nZ2VyLWVzY2FwZWQtc3R5bGU9Im1hcmdpbi1sZWZ0OiAxZW07IG1hcmdpbi1yaWdodDog
MWVtOyIgaHJlZj0iaHR0cDovLzMuYnAuYmxvZ3Nwb3QuY29tLy1WaGlVdWV0ai1fVS9WamVSN0dU
T0hHSS9BQUFBQUFBQTBBWS9FNjRFZTd6Nk9yRS9zMTYwMC9VbnRpdGxlZDElMkIyMDE1LTExLTAz
JTJCMDEtMzgtNDgucG5nIiBpbWFnZWFuY2hvcj0iMSIgc3R5bGU9Im1hcmdpbi1sZWZ0OiAxZW07
IG1hcmdpbi1yaWdodDogMWVtOyI+PGltZyBib3JkZXI9IjAiIGhlaWdodD0iODMiIHNyYz0iaHR0
cHM6Ly8zLmJwLmJsb2dzcG90LmNvbS8tVmhpVXVldGotX1UvVmplUjdHVE9IR0kvQUFBQUFBQUEw
QVkvRTY0RWU3ejZPckUvczMyMC9VbnRpdGxlZDElMkIyMDE1LTExLTAzJTJCMDEtMzgtNDgucG5n
IiB3aWR0aD0iMzIwIj48L2E+PC9wPjxwPjxicj48L3A+PHA+44Oc44K/44Oz44KS5oq844GZ44Go
44K544Op44Kk44OJ44Gr44Gq44Gj44Gm4piT44KS5oq844GZ44Go5YWD44Gr5oi744KL44CCPC9w
PjxwPjxicj48L3A+PHA+44Gd44Gu44G+44G+44Gn44Gv44GZ44G544GmMeaemuOBruOCueODqeOC
pOODieOBq+OBquOBo+OBpuOBl+OBvuOBhuOBruOBp0NlbGwgVG9vbGJhcuOBp1N1YiBzbGlkZeOB
qOOBi+OBq+ioreWumuOBmeOCjOOBsOOCiOOBhOOAgjwvcD48cD48YnI+PC9wPjxwPuOBk+OCjOOB
i+OCieS8muekvuOBp+OBr+OBk+OCjOS9v+OBo+OBpueZuuihqOOBl+OCiOOBhuOAgjwvcD48cD48
YnI+PC9wPjxwPjxicj48L3A+PHA+PGJyPjwvcD48cD48YnI+PC9wPjxwPjxicj48L3A+PC9kaXY+
PGRpdj48L2Rpdj4=">
</div>
<div style="text-align: center;">
<iframe frameborder="0" marginheight="0" marginwidth="0" scrolling="no" src="http://rcm-fe.amazon-adsystem.com/e/cm?lt1=_blank&bc1=000000&IS2=1&bg1=FFFFFF&fc1=000000&lc1=0000FF&t=kenske812-22&o=9&p=8&l=as4&m=amazon&f=ifr&ref=ss_til&asins=4873116554" style="height: 240px; width: 120px;"></iframe>
</div>
</div>
</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-63174406999717259442015-09-23T00:11:00.001+09:002015-09-23T00:24:01.324+09:00ipython notebookでVim(バインディング)<div dir="ltr" style="text-align: left;" trbidi="on">
<div class="markdown-here-wrapper" data-md-url="https://www.blogger.com/blogger.g?blogID=4766414733890203866#editor/target=post;postID=6317440699971725944;onPublishedMenu=allposts;onClosedMenu=allposts;postNum=0;src=link">
<div style="margin: 0px 0px 1.2em !important;">
ipython notebook(Jupyter)でもvim bindingが使えるらしい。<br />
<a href="https://github.com/lambdalisue/jupyter-vim-binding">https://github.com/lambdalisue/jupyter-vim-binding</a></div>
<div style="margin: 0px 0px 1.2em !important;">
ipython notebookで使いにくかった点が1つ改善された。</div>
<div style="margin: 0px 0px 1.2em !important;">
環境はMacにAnacondaを導入している。<br />
基本的な導入は上のGithubのリンクに書いてあるけど…一応メモしておこう。</div>
<div style="margin: 0px 0px 1.2em !important;">
まずはipythonをアップデートして、さらにjupyterを入れる。<br />
(Jupyterを入れないとnotebookモジュールがないと怒られる)</div>
<pre style="font-family: Consolas, Inconsolata, Courier, monospace; font-size: 0.85em; font-size: 1em; line-height: 1.2em; margin: 1.2em 0px;"><code style="background-color: #f8f8f8; border-radius: 3px; border-radius: 3px; border: 1px solid rgb(204, 204, 204); border: 1px solid rgb(234, 234, 234); display: block !important; display: inline; font-family: Consolas, Inconsolata, Courier, monospace; font-size: 0.85em; margin: 0px 0.15em; overflow: auto; padding: 0.5em 0.7em; padding: 0px 0.3em; white-space: pre-wrap; white-space: pre;">$ conda update ipython
$ conda install jupyter
</code></pre>
<div style="margin: 0px 0px 1.2em !important;">
次に以下を実行</div>
<pre style="font-family: Consolas, Inconsolata, Courier, monospace; font-size: 0.85em; font-size: 1em; line-height: 1.2em; margin: 1.2em 0px;"><code style="background-color: #f8f8f8; border-radius: 3px; border-radius: 3px; border: 1px solid rgb(204, 204, 204); border: 1px solid rgb(234, 234, 234); display: block !important; display: inline; font-family: Consolas, Inconsolata, Courier, monospace; font-size: 0.85em; margin: 0px 0.15em; overflow: auto; padding: 0.5em 0.7em; padding: 0px 0.3em; white-space: pre-wrap; white-space: pre;">from notebook.nbextensions import install_nbextension
install_nbextension('https://goo.gl/5TK96v', user=True, destination="vim_binding.js")
</code></pre>
<div style="margin: 0px 0px 1.2em !important;">
一時的に使うときは</div>
<pre style="font-family: Consolas, Inconsolata, Courier, monospace; font-size: 0.85em; font-size: 1em; line-height: 1.2em; margin: 1.2em 0px;"><code style="background-color: #f8f8f8; border-radius: 3px; border-radius: 3px; border: 1px solid rgb(204, 204, 204); border: 1px solid rgb(234, 234, 234); display: block !important; display: inline; font-family: Consolas, Inconsolata, Courier, monospace; font-size: 0.85em; margin: 0px 0.15em; overflow: auto; padding: 0.5em 0.7em; padding: 0px 0.3em; white-space: pre-wrap; white-space: pre;">%%javascript
Jupyter.utils.load_extensions('vim_binding')
</code></pre>
<div style="margin: 0px 0px 1.2em !important;">
ずっと有効にしたいときは</div>
<pre style="font-family: Consolas, Inconsolata, Courier, monospace; font-size: 0.85em; font-size: 1em; line-height: 1.2em; margin: 1.2em 0px;"><code style="background-color: #f8f8f8; border-radius: 3px; border-radius: 3px; border: 1px solid rgb(204, 204, 204); border: 1px solid rgb(234, 234, 234); display: block !important; display: inline; font-family: Consolas, Inconsolata, Courier, monospace; font-size: 0.85em; margin: 0px 0.15em; overflow: auto; padding: 0.5em 0.7em; padding: 0px 0.3em; white-space: pre-wrap; white-space: pre;">%%javascript
Jupyter.notebook.config.update({
'load_extensions': { 'vim_binding': true },
});
</code></pre>
<div style="margin: 0px 0px 1.2em !important;">
これでipython notebook上でVimバインディングが使える。いやー、便利になった。<br />
<div style="text-align: center;">
<br /></div>
<div style="text-align: center;">
<br /></div>
</div>
<div style="font-size: 0em; height: 0; margin: 0; max-height: 0; max-width: 0; overflow: hidden; padding: 0; width: 0;" title="MDH:PHA+aXB5dGhvbiBub3RlYm9vayhKdXB5dGVyKeOBp+OCgnZpbSBiaW5kaW5n44GM5L2/44GI44KL
44KJ44GX44GE44CCPC9wPjxwPmh0dHBzOi8vZ2l0aHViLmNvbS9sYW1iZGFsaXN1ZS9qdXB5dGVy
LXZpbS1iaW5kaW5nPC9wPjxwPjxicj48L3A+PHA+aXB5dGhvbiBub3RlYm9va+OBp+S9v+OBhOOB
q+OBj+OBi+OBo+OBn+eCueOBjDHjgaTmlLnlloTjgZXjgozjgZ/jgII8L3A+PHA+PGJyPjwvcD48
cD7nkrDlooPjga9NYWPjgatBbmFjb25kYeOCkuWwjuWFpeOBl+OBpuOBhOOCi+OAgjwvcD48cD7l
n7rmnKznmoTjgarlsI7lhaXjga/kuIrjga5HaXRodWLjga7jg6rjg7Pjgq/jgavmm7jjgYTjgabj
gYLjgovjgZHjganigKbkuIDlv5zjg6Hjg6LjgZfjgabjgYrjgZPjgYbjgII8L3A+PHA+PGJyPjwv
cD48cD7jgb7jgZrjga9pcHl0aG9u44KS44Ki44OD44OX44OH44O844OI44GX44Gm44CB44GV44KJ
44GranVweXRlcuOCkuWFpeOCjOOCi+OAgjwvcD48cD4oSnVweXRlcuOCkuWFpeOCjOOBquOBhOOB
qG5vdGVib29r44Oi44K444Ol44O844Or44GM44Gq44GE44Go5oCS44KJ44KM44KLKTwvcD48cD5g
YGA8L3A+PHA+JCBjb25kYSB1cGRhdGUgaXB5dGhvbjwvcD48cD4kIGNvbmRhIGluc3RhbGwganVw
eXRlcjwvcD48cD5gYGA8L3A+PHA+PGJyPjwvcD48cD7mrKHjgavku6XkuIvjgpLlrp/ooYw8L3A+
PHA+YGBgPC9wPjxwPmZyb20gbm90ZWJvb2submJleHRlbnNpb25zIGltcG9ydCBpbnN0YWxsX25i
ZXh0ZW5zaW9uPC9wPjxwPmluc3RhbGxfbmJleHRlbnNpb24oJ2h0dHBzOi8vZ29vLmdsLzVUSzk2
dicsIHVzZXI9VHJ1ZSwgZGVzdGluYXRpb249InZpbV9iaW5kaW5nLmpzIik8L3A+PHA+YGBgPC9w
PjxwPjxicj48L3A+PHA+5LiA5pmC55qE44Gr5L2/44GG44Go44GN44GvPC9wPjxwPmBgYDwvcD48
cD4lJWphdmFzY3JpcHQ8L3A+PHA+SnVweXRlci51dGlscy5sb2FkX2V4dGVuc2lvbnMoJ3ZpbV9i
aW5kaW5nJyk8L3A+PHA+YGBgPC9wPjxwPjxicj48L3A+PHA+44Ga44Gj44Go5pyJ5Yq544Gr44GX
44Gf44GE44Go44GN44GvPC9wPjxwPmBgYDwvcD48cD4lJWphdmFzY3JpcHQ8L3A+PHA+SnVweXRl
ci5ub3RlYm9vay5jb25maWcudXBkYXRlKHs8L3A+PHA+Jm5ic3A7ICdsb2FkX2V4dGVuc2lvbnMn
OiB7ICd2aW1fYmluZGluZyc6IHRydWUgfSw8L3A+PHA+fSk7PC9wPjxwPmBgYDwvcD48cD48YnI+
PC9wPjxwPuOBk+OCjOOBp2lweXRob24gbm90ZWJvb2vkuIrjgadWaW3jg5DjgqTjg7Pjg4fjgqPj
g7PjgrDjgYzkvb/jgYjjgovjgILjgYTjgoTjg7zjgIHkvr/liKnjgavjgarjgaPjgZ/jgII8L3A+
PHA+PGJyPjwvcD4=">
<div style="text-align: center;">
</div>
</div>
</div>
<div style="text-align: center;">
<iframe frameborder="0" marginheight="0" marginwidth="0" scrolling="no" src="http://rcm-fe.amazon-adsystem.com/e/cm?lt1=_blank&bc1=000000&IS2=1&bg1=FFFFFF&fc1=000000&lc1=0000FF&t=kenske812-22&o=9&p=8&l=as4&m=amazon&f=ifr&ref=ss_til&asins=4774176311" style="height: 240px; width: 120px;"></iframe>
</div>
</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-8132002961228017992015-09-13T01:57:00.001+09:002015-09-13T23:35:49.642+09:00Iteratively Reweighted Least Squares についてサクッと。<div dir="ltr" style="text-align: left;" trbidi="on">
<div class="markdown-here-wrapper" data-md-url="https://www.blogger.com/blogger.g?blogID=4766414733890203866#editor/target=post;postID=813200296122801799;onPublishedMenu=posts;onClosedMenu=posts;postNum=0;src=link">
<div style="margin: 0px 0px 1.2em !important;">
Iteratively Reweighted Least Squares(IRLS)ってたまに出てくるけど、何か分かってなかった.</div>
<h2 id="least-squares" style="border-bottom-color: rgb(238, 238, 238); border-bottom-style: solid; border-bottom-width: 1px; font-size: 1.4em; font-weight: bold; margin: 1.3em 0px 1em; padding: 0px;">
Least Squares</h2>
<div style="margin: 0px 0px 1.2em !important;">
LSでは線形モデル、</div>
<div style="margin: 0px 0px 1.2em !important;">
<img alt="\bf{y} = \bf{\beta}^T \bf{x} + \bf{\epsilon}" src="https://chart.googleapis.com/chart?cht=tx&chl=%5Cbf%7By%7D%20%3D%20%5Cbf%7B%5Cbeta%7D%5ET%20%5Cbf%7Bx%7D%20%2B%20%5Cbf%7B%5Cepsilon%7D" /></div>
<div style="margin: 0px 0px 1.2em !important;">
というモデルを考えたときにデータとモデルの2乗和誤差を最小にするように<img alt="\bf{\beta}" src="https://chart.googleapis.com/chart?cht=tx&chl=%5Cbf%7B%5Cbeta%7D" />を決定する。このとき誤差<img alt="\bf{\epsilon}" src="https://chart.googleapis.com/chart?cht=tx&chl=%5Cbf%7B%5Cepsilon%7D" />を1つの定数で書ける(データのばらつきが次元で同じ)ときは擬似逆行列を使って係数<img alt="\bf{w}" src="https://chart.googleapis.com/chart?cht=tx&chl=%5Cbf%7Bw%7D" />を決定することができる。導出はないけど、微分が0になるように<img alt="bf{\beta}" src="https://chart.googleapis.com/chart?cht=tx&chl=bf%7B%5Cbeta%7D" />を決定すれば</div>
<div style="margin: 0px 0px 1.2em !important;">
<img alt="\bf{\beta}^{*}= (X^TX)^{-1}X^T \bf{y}" src="https://chart.googleapis.com/chart?cht=tx&chl=%5Cbf%7B%5Cbeta%7D%5E%7B*%7D%3D%20(X%5ETX)%5E%7B-1%7DX%5ET%20%5Cbf%7By%7D" /></div>
<h2 id="weighted-least-squares" style="border-bottom-color: rgb(238, 238, 238); border-bottom-style: solid; border-bottom-width: 1px; font-size: 1.4em; font-weight: bold; margin: 1.3em 0px 1em; padding: 0px;">
Weighted Least Squares</h2>
<div style="margin: 0px 0px 1.2em !important;">
先程はデータの次元でばらつきが同じとしたけど、次元間の分散が違うことだってある。そんな場合にWLSをつかう。調べた感じWikipediaが一番分かりやすい。<br />
<a href="https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Weighted_linear_least_squares">https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Weighted_linear_least_squares</a><br />
この場合、</div>
<div style="margin: 0px 0px 1.2em !important;">
<img alt="\bf{\beta}^* = (X^TWX)^{-1} X^TWy" src="https://chart.googleapis.com/chart?cht=tx&chl=%5Cbf%7B%5Cbeta%7D%5E*%20%3D%20(X%5ETWX)%5E%7B-1%7D%20X%5ETWy" /></div>
<div style="margin: 0px 0px 1.2em !important;">
となる。ただし<img alt="W" src="https://chart.googleapis.com/chart?cht=tx&chl=W" />は対角成分に分散の逆数を並べた行列。例えば、あるデータ次元についてはほぼ確定的なら始めから分散は小さくできるので、それを反映させることができる。式を見ると、データを標準偏差で割っているのと同じになっているので、データが事前にある場合は標準化(standardization)しとけば同じことになりそうだ。</div>
<h2 id="iteratively-reweighted-least-squares-irls-" style="border-bottom-color: rgb(238, 238, 238); border-bottom-style: solid; border-bottom-width: 1px; font-size: 1.4em; font-weight: bold; margin: 1.3em 0px 1em; padding: 0px;">
Iteratively Reweighted Least Squares(IRLS)</h2>
<div style="margin: 0px 0px 1.2em !important;">
上の2つは一発で最適な係数<img alt="\bf{\beta}^*" src="https://chart.googleapis.com/chart?cht=tx&chl=%5Cbf%7B%5Cbeta%7D%5E*" />を求めることができた。これは2乗和誤差が<img alt="\bf{\beta}" src="https://chart.googleapis.com/chart?cht=tx&chl=%5Cbf%7B%5Cbeta%7D" />の2次で書けるからだ。IRLSでは、2次でない場合でも2次で近似しながら重みWを更新して、最適な係数<img alt="\bf{\beta}" src="https://chart.googleapis.com/chart?cht=tx&chl=%5Cbf%7B%5Cbeta%7D" />を求めてくれる。これまたWikipediaがわかりやすかった。<br />
<a href="https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares">https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares</a></div>
<div style="margin: 0px 0px 1.2em !important;">
Lpノルムを最小化したいときは、<br />
<img alt="\bf{\beta}^{(t+1)} = (X^TW^{(t)}X)^{-1} X^TW^{(t)}y" src="https://chart.googleapis.com/chart?cht=tx&chl=%5Cbf%7B%5Cbeta%7D%5E%7B(t%2B1)%7D%20%3D%20(X%5ETW%5E%7B(t)%7DX)%5E%7B-1%7D%20X%5ETW%5E%7B(t)%7Dy" /></div>
<div style="margin: 0px 0px 1.2em !important;">
<img alt="W^{(t)}" src="https://chart.googleapis.com/chart?cht=tx&chl=W%5E%7B(t)%7D" />は対角行列。<img alt="W^{(0)}=\bf{1}" src="https://chart.googleapis.com/chart?cht=tx&chl=W%5E%7B(0)%7D%3D%5Cbf%7B1%7D" />と初期化しておき、対角成分の各要素は以下のように更新していく。<br />
<img alt="w_i^{(t)} = | \bf{y}_i - X_i \beta^{(t)}|^{p-2}" src="https://chart.googleapis.com/chart?cht=tx&chl=w_i%5E%7B(t)%7D%20%3D%20%7C%20%5Cbf%7By%7D_i%20-%20X_i%20%5Cbeta%5E%7B(t)%7D%7C%5E%7Bp-2%7D" /></div>
<div style="margin: 0px 0px 1.2em !important;">
IRLSはスパースな係数を選択するためにL1最適化を使う場合に用いられたりするようだ。<br />
コンピュータビジョン最先端ガイド6巻の3章に少しだけIRLSの記述がある。<br />
Scikit-learnなどではOMPが実装されているみたいだが。<br />
<br /></div>
<iframe frameborder="0" marginheight="0" marginwidth="0" scrolling="no" src="http://rcm-fe.amazon-adsystem.com/e/cm?lt1=_blank&bc1=000000&IS2=1&bg1=FFFFFF&fc1=000000&lc1=0000FF&t=kenske812-22&o=9&p=8&l=as4&m=amazon&f=ifr&ref=ss_til&asins=4915851540" style="height: 240px; width: 120px;"></iframe>
<div style="font-size: 0em; height: 0; margin: 0; max-height: 0; max-width: 0; overflow: hidden; padding: 0; width: 0;" title="MDH:PGRpdiBjbGFzcz0ibWFya2Rvd24taGVyZS13cmFwcGVyIiBkYXRhLWJsb2dnZXItZXNjYXBlZC1k
YXRhLW1kLXVybD0iaHR0cHM6Ly93d3cuYmxvZ2dlci5jb20vYmxvZ2dlci5nP2Jsb2dJRD00NzY2
NDE0NzMzODkwMjAzODY2I2VkaXRvci90YXJnZXQ9cG9zdDtwb3N0SUQ9ODEzMjAwMjk2MTIyODAx
Nzk5IiBkYXRhLWJsb2dnZXItZXNjYXBlZC1tYXJrZG93bi1oZXJlLXdyYXBwZXItY29udGVudC1t
b2RpZmllZD0idHJ1ZSI+PGRpdiBkYXRhLWJsb2dnZXItZXNjYXBlZC1zdHlsZT0ibWFyZ2luOiAw
cHggMHB4IDEuMmVtICFpbXBvcnRhbnQ7Ij48cD5JdGVyYXRpdmVseSBSZXdlaWdodGVkIExlYXN0
IFNxdWFyZXMoSVJMUynjgaPjgabjgZ/jgb7jgavlh7rjgabjgY/jgovjgZHjganjgIHkvZXjgYvl
iIbjgYvjgaPjgabjgarjgYvjgaPjgZ8uPC9wPjxwPjxicj48L3A+PHA+IyMgTGVhc3QgU3F1YXJl
czwvcD48cD5MU+OBp+OBr+e3muW9ouODouODh+ODq+OAgTwvcD48cD48YnI+PC9wPjxwPiRcYmZ7
eX0gPSBcYmZ7XGJldGF9XlQgXGJme3h9Jm5ic3A7KyBcYmZ7XGVwc2lsb259JDwvcD48cD48YnI+
PC9wPjxwPuOBqOOBhOOBhuODouODh+ODq+OCkuiAg+OBiOOBn+OBqOOBjeOBq+ODh+ODvOOCv+OB
qOODouODh+ODq+OBrjLkuZflkozoqqTlt67jgpLmnIDlsI/jgavjgZnjgovjgojjgYbjgaskXGJm
e1xiZXRhfSTjgpLmsbrlrprjgZnjgovjgILjgZPjga7jgajjgY3oqqTlt64kXGJme1xlcHNpbG9u
fSTjgpIx44Gk44Gu5a6a5pWw44Gn5pu444GR44KLKOODh+ODvOOCv+OBruOBsOOCieOBpOOBjeOB
jOasoeWFg+OBp+WQjOOBmCnjgajjgY3jga/mk6zkvLzpgIbooYzliJfjgpLkvb/jgaPjgabkv4Lm
lbAkXGJme3d9JOOCkuaxuuWumuOBmeOCi+OBk+OBqOOBjOOBp+OBjeOCi+OAguWwjuWHuuOBr+OB
quOBhOOBkeOBqeOAgeW+ruWIhuOBjDDjgavjgarjgovjgojjgYbjgaskYmZ7XGJldGF9JOOCkuax
uuWumuOBmeOCjOOBsDwvcD48cD48YnI+PC9wPjxwPiRcYmZ7XGJldGF9XnsqfT0gKFheVFgpXnst
MX1YXlQgXGJme3l9JDwvcD48cD48YnI+PC9wPjxwPiMjIFdlaWdodGVkIExlYXN0IFNxdWFyZXM8
L3A+PHA+5YWI56iL44Gv44OH44O844K/44Gu5qyh5YWD44Gn44Gw44KJ44Gk44GN44GM5ZCM44GY
44Go44GX44Gf44GR44Gp44CB5qyh5YWD6ZaT44Gu5YiG5pWj44GM6YGV44GG44GT44Go44Gg44Gj
44Gm44GC44KL44CC44Gd44KT44Gq5aC05ZCI44GrV0xT44KS44Gk44GL44GG44CC6Kq/44G544Gf
5oSf44GYV2lraXBlZGlh44GM5LiA55Wq5YiG44GL44KK44KE44GZ44GE44CCPC9wPjxwPmh0dHBz
Oi8vZW4ud2lraXBlZGlhLm9yZy93aWtpL0xpbmVhcl9sZWFzdF9zcXVhcmVzXyhtYXRoZW1hdGlj
cykjV2VpZ2h0ZWRfbGluZWFyX2xlYXN0X3NxdWFyZXM8L3A+PHA+44GT44Gu5aC05ZCI44CBPC9w
PjxwPjxicj48L3A+PHA+JFxiZntcYmV0YX1eKiA9IChYXlRXWCleey0xfSBYXlRXeSQ8L3A+PHA+
PGJyPjwvcD48cD7jgajjgarjgovjgILjgZ/jgaDjgZckVyTjga/lr77op5LmiJDliIbjgavliIbm
laPjga7pgIbmlbDjgpLkuKbjgbnjgZ/ooYzliJfjgILkvovjgYjjgbDjgIHjgYLjgovjg4fjg7zj
gr/mrKHlhYPjgavjgaTjgYTjgabjga/jgbvjgbznorrlrprnmoTjgarjgonlp4vjgoHjgYvjgonl
iIbmlaPjga/lsI/jgZXjgY/jgafjgY3jgovjga7jgafjgIHjgZ3jgozjgpLlj43mmKDjgZXjgZvj
govjgZPjgajjgYzjgafjgY3jgovjgILlvI/jgpLopovjgovjgajjgIHjg4fjg7zjgr/jgpLmqJnm
upblgY/lt67jgaflibLjgaPjgabjgYTjgovjga7jgajlkIzjgZjjgavjgarjgaPjgabjgYTjgovj
ga7jgafjgIHjg4fjg7zjgr/jgYzkuovliY3jgavjgYLjgovloLTlkIjjga/mqJnmupbljJYoc3Rh
bmRhcmRpemF0aW9uKeOBl+OBqOOBkeOBsOWQjOOBmOOBk+OBqOOBq+OBquOCiuOBneOBhuOBoOOA
gjwvcD48cD48YnI+PC9wPjxwPjxicj48L3A+PHA+PGJyPjwvcD48cD4jIyBJdGVyYXRpdmVseSBS
ZXdlaWdodGVkIExlYXN0IFNxdWFyZXMoSVJMUyk8L3A+PHA+5LiK44GuMuOBpOOBr+S4gOeZuuOB
p+acgOmBqeOBquS/guaVsCRcYmZ7XGJldGF9Xiok44KS5rGC44KB44KL44GT44Go44GM44Gn44GN
44Gf44CC44GT44KM44GvMuS5l+WSjOiqpOW3ruOBjCRcYmZ7XGJldGF9JOOBrjLmrKHjgafmm7jj
gZHjgovjgYvjgonjgaDjgIJJUkxT44Gn44Gv44CBMuasoeOBp+OBquOBhOWgtOWQiOOBp+OCgjLm
rKHjgafov5HkvLzjgZfjgarjgYzjgonph43jgb9X44KS5pu05paw44GX44Gm44CB5pyA6YGp44Gq
5L+C5pWwJFxiZntcYmV0YX0k44KS5rGC44KB44Gm44GP44KM44KL44CC44GT44KM44G+44GfV2lr
aXBlZGlh44GM44KP44GL44KK44KE44GZ44GL44Gj44Gf44CCPC9wPjxwPmh0dHBzOi8vZW4ud2lr
aXBlZGlhLm9yZy93aWtpL0l0ZXJhdGl2ZWx5X3Jld2VpZ2h0ZWRfbGVhc3Rfc3F1YXJlczwvcD48
cD48YnI+PC9wPjxwPkxw44OO44Or44Og44KS5pyA5bCP5YyW44GX44Gf44GE44Go44GN44Gv44CB
PC9wPjxwPiRcYmZ7XGJldGF9XnsodCsxKX0gPSAoWF5UV157KHQpfVgpXnstMX0gWF5UV157KHQp
fXkkPC9wPjxwPjxicj48L3A+PHA+JFdeeyh0KX0k44Gv5a++6KeS6KGM5YiX44CCJFdeeygwKX09
XGJmezF9JOOBqOWIneacn+WMluOBl+OBpuOBiuOBjeOAgeWvvuinkuaIkOWIhuOBruWQhOimgee0
oOOBr+S7peS4i+OBruOCiOOBhuOBq+abtOaWsOOBl+OBpuOBhOOBj+OAgjwvcD48cD4kd19pXnso
dCl9ID0gfCBcYmZ7eX1faSAtIFhfaSBcYmV0YV57KHQpfXxee3AtMn0kPC9wPjxwPjxicj48L3A+
PHA+PGJyPjwvcD48cD5JUkxT44Gv44K544OR44O844K544Gq5L+C5pWw44KS6YG45oqe44GZ44KL
44Gf44KB44GrTDHmnIDpganljJbjgpLkvb/jgYbloLTlkIjjgavnlKjjgYTjgonjgozjgZ/jgorj
gZnjgovjgojjgYbjgaDjgII8L3A+PHA+PGJyPjwvcD48cD48YnI+PC9wPjxwPjxicj48L3A+PC9k
aXY+PC9kaXY+">
</div>
</div>
</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-66834746981270759772015-09-01T20:59:00.001+09:002015-09-01T20:59:36.803+09:00Python でPRML7.2章 RVMを実装。<div dir="ltr" style="text-align: left;" trbidi="on">
いつ使うかは分からんが、RVMを実装しておこう。<br />
7.2章の数式をそのまま書いた。動作確認は1次元のガウス基底で行った。<br />
<br />
<pre style="background: #fdf6e3; color: #586e75;"><span style="color: #748b00;">def</span> <span style="color: #268bd2;">rvm</span>(t, PHI, alphas0, beta0, maxloop<span style="color: #859900;">=</span><span style="color: #d33682;">10</span>):
<span style="color: #269186;"><span style="color: #c60000;">"""</span>
t: observed values.
X: designe matrix
alphas0: initial weights for each parameter.
beta0: an initial precision value.
returns
alpha: prior weights
beta: prior precision
m: prior means
V: prior covariance matrix
See PRML 7.2 for details.
** to many loop causes 0 devision for some reason.
<span style="color: #c60000;">"""</span></span>
beta <span style="color: #859900;">=</span> beta0
alphas <span style="color: #859900;">=</span> alphas0
<span style="color: #859900;">for</span> i <span style="color: #859900;">in</span> <span style="color: #268bd2;">range</span>(maxloop):
A <span style="color: #859900;">=</span> np.diag(alphas)
V <span style="color: #859900;">=</span> np.linalg.inv(A<span style="color: #859900;">+</span>beta<span style="color: #859900;">*</span>np.dot(PHI.T, PHI))
m <span style="color: #859900;">=</span> beta<span style="color: #859900;">*</span>np.dot(V, np.dot(PHI.T, t))
gammas <span style="color: #859900;">=</span> <span style="color: #d33682;">1</span> <span style="color: #859900;">-</span> alphas <span style="color: #859900;">*</span> np.diag(V)
<span style="color: #93a1a1;">#new alphas and beta</span>
alphas <span style="color: #859900;">=</span> gammas <span style="color: #859900;">/</span> m<span style="color: #859900;">*</span><span style="color: #859900;">*</span><span style="color: #d33682;">2</span>
beta <span style="color: #859900;">=</span> (N <span style="color: #859900;">-</span> np.sum(gammas)) <span style="color: #859900;">/</span> np.linalg.norm(t <span style="color: #859900;">-</span> np.dot(PHI, m))
<span style="color: #859900;">return</span> alphas, beta, m, V
</pre>
<br />
<br />
10個のガウス基底でフィッティングをした。<br />
今回書いたコードではループの回数が多すぎるとzero divisionでエラーとなった。でもループは10回くらいで十分フィットしているように見えた。<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgX057Cz_VYFyvlj2GEvlmf3spEC16YO1_B1EOgi4R5DxS_wDAGWP_bxadNtCVk02o-0wbWLXrAEmzCJiNC_AwlG1eGO0lYcJ4mJEp3WOZlBYGS4221M1qJvIVMMcc1ZoyzWYuYaXMoA0k/s1600/download.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" height="214" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgX057Cz_VYFyvlj2GEvlmf3spEC16YO1_B1EOgi4R5DxS_wDAGWP_bxadNtCVk02o-0wbWLXrAEmzCJiNC_AwlG1eGO0lYcJ4mJEp3WOZlBYGS4221M1qJvIVMMcc1ZoyzWYuYaXMoA0k/s320/download.png" width="320" /></a></div>
<br />
<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-59974116215974505512015-08-24T00:31:00.001+09:002015-08-24T00:31:40.882+09:00Python3 で OpenCV3を使う(Windows)。<div dir="ltr" style="text-align: left;" trbidi="on">
最近やりたいことはほとんどPython3で出来ているが、OpenCVを使う時だけpython2に切り替えている。OpenCVもPython3に対応したようだが、配布バイナリはPython2のものしか提供されていない。オフィシャルでないけど、サクッとPython3にOpenCV3をインストールする方法を見つけたので紹介する。<br />
<br />
<br />
まず、ここから.whlファイルをダウンロードする。<br />
(このサイト経由だと余計なものをダウンロードしなくて済むのでよい)<br />
<a href="http://www.lfd.uci.edu/~gohlke/pythonlibs/#opencv">http://www.lfd.uci.edu/~gohlke/pythonlibs/#opencv</a><br />
<br />
whlはpipでインストールできるので、<br />
<br />
<pre style="background: #fdf6e3; color: #586e75;">pip install opencv_python-3.0.0-cp34-none-win_amd64.whl
</pre>
<br />
とする。<br />
<br />
前はソースからコンパイルとかしたけど、こんなにラクになるとは。ありがとうございます。</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-72803274507435212972015-07-19T22:30:00.000+09:002015-07-19T22:30:03.416+09:00matplotlibのcolormap<div dir="ltr" style="text-align: left;" trbidi="on">
scipy2015に興味深い講演があった。<br />
A Better Default Colormap for Matplotlib<br />
<a href="https://www.youtube.com/watch?v=xAoljeRJ3lU">https://www.youtube.com/watch?v=xAoljeRJ3lU</a><br />
<br />
<br />
Matplotlibの現在のデフォルトの色はjet。<br />
Matplotlib2.0からはデフォルトの色が'viridis'に変わる。<br />
<a href="http://betterfigures.org/2015/07/10/a-welcome-development-for-matplotlib/">http://betterfigures.org/2015/07/10/a-welcome-development-for-matplotlib/</a><br />
'viridis'は4つの候補のうちの1つで、投票で選ばれたらしい(緑が入ってるから)。4つともMatplotlib1.5から導入はされるらしい。<br />
<br />
なにが良くなったかというと、<br />
色空間を連続的に変化させるだけでなく、人間の視覚を通した時にスムーズに変化するようになっているらしい。こうすることによって、より細かい変化に気付くことができる。たとえば、医師の誤診断とかも減るとか。更に、白黒にしたとしても分かるし、色盲の人でもわかりやすいように設計されているらしい。<br />
<br />
ここには動画もある。<br />
<a href="https://bids.github.io/colormap/">https://bids.github.io/colormap/</a><br />
<br />
ちなみに今使えるcolromapサンプルはここにある。<br />
<a href="http://matplotlib.org/examples/color/colormaps_reference.html">http://matplotlib.org/examples/color/colormaps_reference.html</a><br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-77342428824782044962015-06-28T01:53:00.001+09:002015-06-28T01:56:23.841+09:00Pythonで多項式回帰<div dir="ltr" style="text-align: left;" trbidi="on">
pythonで最小二乗法で係数を求めるのには関数を定義しなければいけないが、高次の多項式だとクロスタームやらいっぱい出てきて書くのが大変だ。<br />
<br />
scikit-learnのPolynomialFeaturesはそこを自動で計算してくれる。<br />
例えば2次の場合degree=2とすれば、<br />
<br />
[x1, y1] -> [x1**2, x1*y1, y1**2, x1, y1]<br />
<br />
を計算してくれる。<br />
<br />
それを使って多項式近似をする。<br />
以下メモ。<br />
Pipelineについては<br />
<a href="http://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html">http://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html</a><br />
<br />
<pre style="background: #fdf6e3; color: #586e75;"><span style="color: #bd3800;">from</span> sklearn.linear_model <span style="color: #bd3800;">import</span> Ridge
<span style="color: #bd3800;">from</span> sklearn.preprocessing <span style="color: #bd3800;">import</span> PolynomialFeatures
<span style="color: #bd3800;">from</span> sklearn.pipeline <span style="color: #bd3800;">import</span> make_pipeline
<span style="color: #93a1a1;">#original function</span>
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">f</span>(x):
<span style="color: #859900;">return</span> <span style="color: #d33682;">2</span><span style="color: #859900;">*</span>x[<span style="color: #d33682;">0</span>]<span style="color: #859900;">*</span><span style="color: #859900;">*</span><span style="color: #d33682;">2</span> <span style="color: #859900;">+</span> <span style="color: #d33682;">3</span><span style="color: #859900;">*</span>x[<span style="color: #d33682;">0</span>]<span style="color: #859900;">*</span>x[<span style="color: #d33682;">1</span>] <span style="color: #859900;">+</span> <span style="color: #d33682;">4</span><span style="color: #859900;">*</span>x[<span style="color: #d33682;">1</span>]<span style="color: #859900;">*</span><span style="color: #859900;">*</span><span style="color: #d33682;">2</span> <span style="color: #859900;">+</span> <span style="color: #d33682;">2</span><span style="color: #859900;">*</span>x[<span style="color: #d33682;">0</span>] <span style="color: #859900;">+</span> <span style="color: #d33682;">3</span><span style="color: #859900;">*</span>x[<span style="color: #d33682;">1</span>] <span style="color: #859900;">+</span> <span style="color: #d33682;">4</span>
<span style="color: #93a1a1;">#observed data</span>
xx, yy <span style="color: #859900;">=</span> np.meshgrid(np.linspace(<span style="color: #859900;">-</span><span style="color: #d33682;">2</span>, <span style="color: #d33682;">2</span>, <span style="color: #d33682;">10</span>), np.linspace(<span style="color: #859900;">-</span><span style="color: #d33682;">2</span>, <span style="color: #d33682;">2</span>, <span style="color: #d33682;">10</span>))
X <span style="color: #859900;">=</span> np.vstack((xx.reshape(<span style="color: #859900;">-</span><span style="color: #d33682;">1</span>), yy.reshape(<span style="color: #859900;">-</span><span style="color: #d33682;">1</span>))).T
y <span style="color: #859900;">=</span> <span style="color: #859900;">list</span>(<span style="color: #268bd2;">map</span>(f, X))
<span style="color: #93a1a1;">#fit the data</span>
degree <span style="color: #859900;">=</span> <span style="color: #d33682;">2</span>
model <span style="color: #859900;">=</span> make_pipeline(PolynomialFeatures(degree), Ridge())
model.fit(X, y)
<span style="color: #93a1a1;">#accessing to the coeffs</span>
<span style="color: #93a1a1;">#(model is a Pipeline.)</span>
<span style="color: #859900;">print</span>(model.steps[<span style="color: #d33682;">1</span>][<span style="color: #d33682;">1</span>].coef_)
</pre>
<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-42189858053348632642015-06-27T23:41:00.001+09:002015-06-27T23:41:30.013+09:00Python でサクッとLeast Square fit<div dir="ltr" style="text-align: left;" trbidi="on">
Pythonで最小二乗を使うときはサクッと以下のように行おう。<br />
<br />
<pre style="background: #fdf6e3; color: #586e75;"><span style="color: #bd3800;">import</span> numpy <span style="color: #859900;">as</span> np
<span style="color: #bd3800;">from</span> scipy.optimize <span style="color: #bd3800;">import</span> leastsq
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">func</span>(c, x):
<span style="color: #859900;">return</span> c[<span style="color: #d33682;">0</span>]<span style="color: #859900;">*</span>x<span style="color: #859900;">*</span><span style="color: #859900;">*</span><span style="color: #d33682;">2</span> <span style="color: #859900;">+</span> c[<span style="color: #d33682;">1</span>]<span style="color: #859900;">*</span>x <span style="color: #859900;">+</span> c[<span style="color: #d33682;">2</span>]
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">errfunc</span>(c, x, t):
<span style="color: #859900;">return</span> t <span style="color: #859900;">-</span> func(c, x)
<span style="color: #93a1a1;">#observed data and initial guess</span>
x <span style="color: #859900;">=</span> np.linspace(<span style="color: #859900;">-</span><span style="color: #d33682;">5</span>, <span style="color: #d33682;">5</span>, <span style="color: #d33682;">100</span>)
t <span style="color: #859900;">=</span> <span style="color: #d33682;">4</span><span style="color: #859900;">*</span>x<span style="color: #859900;">*</span><span style="color: #859900;">*</span><span style="color: #d33682;">2</span> <span style="color: #859900;">+</span> <span style="color: #d33682;">1</span>
c0 <span style="color: #859900;">=</span> [<span style="color: #d33682;">0</span>, <span style="color: #d33682;">0</span>, <span style="color: #d33682;">0</span>]
leastsq(errfunc, c0, args<span style="color: #859900;">=</span>(x, t))
</pre>
<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-34519715513013005332015-06-02T22:00:00.003+09:002015-06-02T22:01:25.876+09:00AnacondaでPython2と3のIPython notebookを起動する。<div dir="ltr" style="text-align: left;" trbidi="on">
Pythonで計算をするときはAnacondaを入れとけばとりあえず色々できる。<br />
<br />
Python2と3の両方を使うためにAnacondaとAnaconda3を2つインストールしていたが、そんなことしなくてもよいことがわかった。Anaconda3でPython2を扱う方法をメモしておく。<br />
<br />
<br />
まず、AnacondaでPython2の環境を作る。これはいろいろなサイトに載っている通りに行う。<br />
<div>
<br /></div>
<div>
<pre style="background: #fdf6e3; color: #586e75;"></pre>
<pre style="background: #fdf6e3; color: #586e75;"></pre>
<pre style="background: #fdf6e3; color: #586e75;"></pre>
<pre style="background: #fdf6e3; color: #586e75;">>conda create -n python2 python=2.7 anaconda
</pre>
<pre style="background: #fdf6e3; color: #586e75;"></pre>
<pre style="background: #fdf6e3; color: #586e75;"></pre>
<pre style="background: #fdf6e3; color: #586e75;"></pre>
<br />
IPython notebookをPython2で使うときは<br />
<br /></div>
<pre style="background: #fdf6e3; color: #586e75;"></pre>
<pre style="background: #fdf6e3; color: #586e75;"></pre>
<pre style="background: #fdf6e3; color: #586e75;">>activate python2
>ipython notebook
</pre>
<pre style="background: #fdf6e3; color: #586e75;"></pre>
<pre style="background: #fdf6e3; color: #586e75;"></pre>
<div>
<br />
<br />
SublimetextでPython2を使いたいときは<br />
<br /></div>
<pre style="background: #fdf6e3; color: #586e75;">{
<span style="color: #269186;"><span style="color: #c60000;">"</span>shell<span style="color: #c60000;">"</span></span>: <span style="color: #b58900;">true</span>,
<span style="color: #269186;"><span style="color: #c60000;">"</span>cmd<span style="color: #c60000;">"</span></span>: <span style="color: #d31e1e;">[</span><span style="color: #269186;"><span style="color: #c60000;">"</span>C:<span style="color: #cb4b16;">\\</span>Users<span style="color: #cb4b16;">\\</span>myname<span style="color: #cb4b16;">\\</span>Anaconda3<span style="color: #cb4b16;">\\</span>envs<span style="color: #cb4b16;">\\</span>python2<span style="color: #cb4b16;">\\</span>python.exe<span style="color: #c60000;">"</span></span>, <span style="color: #269186;"><span style="color: #c60000;">"</span>-u<span style="color: #c60000;">"</span></span>, <span style="color: #269186;"><span style="color: #c60000;">"</span>$file<span style="color: #c60000;">"</span></span><span style="color: #d31e1e;">]</span>,
<span style="color: #269186;"><span style="color: #c60000;">"</span>file_regex<span style="color: #c60000;">"</span></span>: <span style="color: #269186;"><span style="color: #c60000;">"</span>^[ ]*File <span style="color: #cb4b16;">\"</span>(...*?)<span style="color: #cb4b16;">\"</span>, line ([0-9]*)<span style="color: #c60000;">"</span></span>
}
</pre>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-4474442454916114992015-05-24T22:57:00.003+09:002015-05-24T22:57:52.897+09:00Lispを保存<div dir="ltr" style="text-align: left;" trbidi="on">
Land of Lispを読み始めた。 といっても実は2回目の挑戦である。一回目に読んだときはパソコンの電源を落としたら、インタプリタに打ち込んだ変数やら関数がなくなってしまって挫折した。<br />
<br />
今回はそうならないようにLispの環境を保存できるようにする。<br />
インストールしたのは<a href="http://www.clisp.org/" target="_blank">GNU CLisp</a>.<br />
<br />
保存する場合は<br />
<pre style="background: rgb(253, 246, 227);"><span style="color: #bb3700;">
</span></pre>
<pre style="background: rgb(253, 246, 227);"><span style="color: #bb3700;">(ext:saveinitmem "lispinit.mem")</span></pre>
<pre style="background: rgb(253, 246, 227);"><span style="color: #bb3700;">
</span></pre>
<div>
<br /></div>
読み込んで始める場合は、<br />
<pre style="background: rgb(253, 246, 227);"></pre>
<pre style="background: rgb(253, 246, 227);">clisp -norc -M lispinit.mem </pre>
<pre style="background: rgb(253, 246, 227);"></pre>
<br />
<br />
インタプリタで開発は書きづらいのでエディタで書いてそれを上のコマンド保存、読み出しをするのがやりやすそうだ(個人的には)<br />
<br />
<br />
<br />
<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-71393152283818803262015-03-29T19:55:00.001+09:002015-03-29T19:55:30.473+09:00Xperia Z SO-02E MVNOのSIMでもテザリングする。<div dir="ltr" style="text-align: left;" trbidi="on">
ドコモのXperia Zを使って2年がたった。<br />
乗り換えばかりを優遇する3大キャリアに嫌気がさしたのと、普段月2Gも使って居ないのでMVNOにMNPした。<br />
<br />
たまに電車の中とかでテザリングしたいときがあったけど、そのままのXperiaではテザリングできないらしい。いやーひどい。<br />
<br />
調べてみると、WifiテザリングはRoot化しないと無理とかいう情報ばかりだったけど、したのリンクの記事にRoot化しなくても行けると書いてあった。<br />
<a href="http://androplus.org/Entry/282">http://androplus.org/Entry/282</a><br />
<br />
ここに書いてある通りにやればadb.exeもダウンロードできてうまくいった。<br />
<a href="http://xperia-freaks.org/2015/03/21/docomo-tethering/">http://xperia-freaks.org/2015/03/21/docomo-tethering/</a><br />
<br />
<br />
テザリングはうまくいったけど、たまにspモードがどうとかいう警告がでるのがうざい。<br />
これはRoot化しないときえないのか…<br />
<br />
<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-67796272359576433652014-11-18T20:10:00.000+09:002014-11-18T20:10:02.204+09:00Juliaで外部のプログラムを動かすときのメモ<div dir="ltr" style="text-align: left;" trbidi="on">
Juliaで外部のプログラムを動かすときのメモ。<br />
PythonではPopenとかを使ってたので、Juliaでもそれと同じようなものがあるかなーと思ったけど、なかなか見つからなかった。<br />
<br />
公式マニュアルはこれ。<br />
<a href="http://julia.readthedocs.org/en/latest/manual/running-external-programs/">http://julia.readthedocs.org/en/latest/manual/running-external-programs/</a><br />
<br />
外部プログラムを起動して、そこから標準入力を扱いたいときは、こっちが参考になった。<a href="http://blog.leahhanson.us/running-shell-commands-from-julia.html">http://blog.leahhanson.us/running-shell-commands-from-julia.html</a><br />
<br />
概要は、<br />
プログラムが起動して、<br />
write()で書き込み。書き込んだらその出力がsoにたまっていくので、<br />
readall()でその出力を得る。<br />
<div style="background: #f7f7f7; border-radius: 10px; border: 2px dotted #CCCCCC; padding: 10px;">
(so,si,pr) = readandwrite(`yourprogram`)<br />
write(si, "command1")<br />
write(si, "command2")<br />
close(si)<br />
output = readall(so)</div>
<br />
<br />
これも一応メモ。<br />
<a href="http://julialang.org/blog/2013/04/put-this-in-your-pipe/">http://julialang.org/blog/2013/04/put-this-in-your-pipe/</a><br />
<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-67850689531281032642014-11-11T02:10:00.004+09:002014-11-11T02:10:51.463+09:00Novel Class<div dir="ltr" style="text-align: left;" trbidi="on">
アメリカ滞在の後半はNovel Classを取った。文字通り小説の授業だ。<br />
<br />
半期に3つの小説を読んだ。1つはベトナム戦争についての話。もうひとつはThe Catcher In The Rye。最後にSidartha。Sidarthaはヘルマン・ヘッセの小説。どれも一字一句辞書で調べながらじっくり読んだんで、ものすごく思い出に残っている。<br />
<br />
ベトナム戦争についての話はあんまり覚えてない。ただ、そのなかの参考資料にBreaking Awayという映画があった気がする。本当にこの小説の授業で出会ったものかは定かではないえけど、ものすごく感傷にひたったのを覚えている。Breaking Awayは青春を舞台にしたもので4人の登場人物がそれぞれ精一杯行きている。あーおれはこんなところで独りでなにをやっているんだろう。と心の通った友達とのチームプレイがいかに大事かが身に染みたのを覚えている。<br />
<br />
2つ目のThe Catcher In the Ryeは有名な小説だが、村上春樹の日本語訳はまったくおすすめしない。帰国後日本語訳を手にとって読んでみたけど、まったく感情移入することが出来んかった。こっちはあまり前向きな話ではない。最後にDon't tell anybody anything. You will miss somebody.というのがすごく思い出に残っている。他にもI have to admit itとかThat kills me.とかものすごく当時のこころに響くセリフばかりだった。<br />
<br />
Sidarthaは難しい本が読めなかったので、優しい英訳の本を読んだ。Sidarthaが悟りをひらくまでの苦悩が思い出に残っているが、詳しい内容は忘れたw<br />
<br />
<br />
それぞれの本を読んだ後の課題は自由。どんな方法でも読んだ感想を表現すればよい。模型をつくった人もいれば、冊子を作った人もいた。日本のいわゆる読書感想文みたいに決まった書式はない。自分の得意な方法で表現すればよかった。いま日本にこんな教育をするところはあるんかな?<br />
<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-79681595370490051722014-11-06T22:42:00.000+09:002014-11-06T22:42:19.190+09:00Julia Array操作の基本をメモ。<div dir="ltr" style="text-align: left;" trbidi="on">
Pythonでやっていたlistの操作にすっかり慣れてしまったけど、JuliaのArrayではどうやるか確かめる。<br />
Julia0.3で試して動いたものを掲載。<br />
<h4 id="配列の生成">
配列の生成</h4>
<pre class="prettyprint"><code class=" hljs livecodeserver"><span class="hljs-operator">a</span> = [i <span class="hljs-keyword">for</span> i <span class="hljs-operator">in</span> <span class="hljs-number">1</span>:<span class="hljs-number">10</span>] <span class="hljs-comment">#1から10を要素に持つ配列生成</span>
<span class="hljs-operator">a</span> = [<span class="hljs-number">1</span>:<span class="hljs-number">10</span>] <span class="hljs-comment">#同上</span>
b = reshape([<span class="hljs-number">1</span>:<span class="hljs-number">9</span>], (<span class="hljs-number">3</span>,<span class="hljs-number">3</span>)) <span class="hljs-comment">#3x3の2次元配列。</span>
b = [<span class="hljs-number">1</span> <span class="hljs-number">2</span> <span class="hljs-number">3</span>; <span class="hljs-number">4</span> <span class="hljs-number">5</span> <span class="hljs-number">6</span>; <span class="hljs-number">7</span> <span class="hljs-number">8</span> <span class="hljs-number">9</span>] <span class="hljs-comment">#3x3の2次元配列。上とは行、列が逆(転置の関係)</span></code></pre>
<h4 id="要素数">
要素数</h4>
<pre class="prettyprint"><code class=" hljs livecodeserver"><span class="hljs-built_in">length</span>(<span class="hljs-operator">a</span>)</code></pre>
<h4 id="要素アクセス">
要素アクセス</h4>
<pre class="prettyprint"><code class=" hljs livecodeserver"><span class="hljs-operator">a</span>[<span class="hljs-number">1</span>]
<span class="hljs-operator">a</span>[<span class="hljs-number">1</span>:<span class="hljs-function"><span class="hljs-keyword">end</span>] #要素<span class="hljs-title">1</span>から最後まで</span>
<span class="hljs-operator">a</span>[<span class="hljs-number">2</span>:<span class="hljs-number">3</span>:<span class="hljs-number">8</span>] <span class="hljs-comment">#要素2から8を3個おきに</span></code></pre>
<h4 id="追加挿入">
追加,挿入</h4>
<pre class="prettyprint"><code class=" hljs livecodeserver">push!(<span class="hljs-operator">a</span>, <span class="hljs-number">13</span>) <span class="hljs-comment">#要素と同じ型追加</span>
append!(<span class="hljs-operator">a</span>, [<span class="hljs-number">1</span>,<span class="hljs-number">2</span>,<span class="hljs-number">3</span>]) <span class="hljs-comment">#末尾に[1,2,3]を追加</span>
prepend!(<span class="hljs-operator">a</span>, [<span class="hljs-number">1</span>,<span class="hljs-number">2</span>,<span class="hljs-number">3</span>]) <span class="hljs-comment">#先頭に[1,2,3]を追加</span>
insert!(<span class="hljs-operator">a</span>, <span class="hljs-number">3</span>, <span class="hljs-number">15</span>) <span class="hljs-comment">#3番目に15を挿入する。(要素数は増える)</span>
vcat(<span class="hljs-operator">a</span>, [<span class="hljs-number">11</span>,<span class="hljs-number">12</span>,<span class="hljs-number">13</span>]) <span class="hljs-comment">#aに末尾に[11,12,13]を追加(aは非破壊になる)</span>
hcat(b, [<span class="hljs-number">10</span>,<span class="hljs-number">11</span>,<span class="hljs-number">12</span>]) <span class="hljs-comment">#bに列を追加</span>
vcat(b, [<span class="hljs-number">10</span>,<span class="hljs-number">11</span>,<span class="hljs-number">12</span>]<span class="hljs-string">') #bに行を追加('</span>は転置)</code></pre>
要素数があらかじめわかっているなら、始めに配列を確保して代入したほうが速いっぽい。<br />
<h4 id="置換">
置換</h4>
<pre class="prettyprint"><code class=" hljs livecodeserver"><span class="hljs-operator">a</span>[<span class="hljs-number">1</span>] = <span class="hljs-number">15</span>
<span class="hljs-operator">a</span>[<span class="hljs-number">1</span>:<span class="hljs-number">4</span>] = [<span class="hljs-number">4</span>,<span class="hljs-number">4</span>,<span class="hljs-number">4</span>,<span class="hljs-number">4</span>]
<span class="hljs-comment">#a[2:5]を削除して[0,0,0]を挿入。戻り値は削除した元のaの[2:5]の値。</span>
splice!(<span class="hljs-operator">a</span>, <span class="hljs-number">2</span>:<span class="hljs-number">5</span>, [<span class="hljs-number">0</span>,<span class="hljs-number">0</span>,<span class="hljs-number">0</span>]) </code></pre>
<h4 id="削除">
削除</h4>
<pre class="prettyprint"><code class=" hljs livecodeserver">shift!(<span class="hljs-operator">a</span>) <span class="hljs-comment">#最初の要素を削除</span>
pop!(<span class="hljs-operator">a</span>) <span class="hljs-comment">#最後の要素を削除</span></code></pre>
<h4 id="検索">
検索</h4>
<pre class="prettyprint"><code class=" hljs livecodeserver"><span class="hljs-number">3</span> <span class="hljs-operator">in</span> <span class="hljs-operator">a</span>
findfirst(<span class="hljs-operator">a</span>, <span class="hljs-number">4</span>) <span class="hljs-comment">#はじめに4が出る要素番号、なければ0を返す</span>
findin(<span class="hljs-operator">a</span>, <span class="hljs-number">4</span>) <span class="hljs-comment">#4の存在する要素数をArrayに入れて返す)</span></code></pre>
</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-11855346697035195642014-11-06T00:20:00.000+09:002014-11-06T00:35:14.131+09:00Lytroが日本に上陸してアクセスが増えたのでLytroのすごいと思うところをのべる。<div dir="ltr" style="text-align: left;" trbidi="on">
Lytroが日本にやってくるらしい。<br />
<a href="http://dc.watch.impress.co.jp/docs/news/20141105_674603.html">http://dc.watch.impress.co.jp/docs/news/20141105_674603.html</a><br />
<br />
なんか記事のアクセスがすごい。<br />
これは良い機会とばかりにLytroがすごいと思うところを挙げる。1つは文化的内容。もう一つは技術的内容。<br />
<br />
<h3 style="text-align: left;">
1. アナログとデジタルの融合が進んでいる。</h3>
上の記事にレンズの断面図がのってある。35mm判換算で30-250でFナンバーは2.0らしい。最近のカメラと比べると広角側が少し足りないと思われるかもしれない。注目すべきはレンズのスペックが高いにも関わらず、枚数12枚で非球面レンズを使っていないことだ。<br />
<a href="http://dc.watch.impress.co.jp/img/dcw/docs/674/603/html/17.jpg.html">http://dc.watch.impress.co.jp/img/dcw/docs/674/603/html/17.jpg.html</a><br />
<br />
確かIllumの発表の時に言及しとったと思うけど、レンズの収差もデジタルで補正できるような最適設定にしているらしい。収差のデジタル補正と言えば富士フイルムとかキヤノンにそういう機能(デコンボリューション)がついていたと思うけど、こういう技術はどちらかというと、まだレンズの最適化で取れない部分をデジタルで補正しているイメージ。しかし、Lytroの場合発想が逆で、デジタルで補正できるところまでレンズの収差を出しているようだ。Ngさんの博士論文にもそういう研究はあったし。もちろん日本の企業でも「倍率色収差」とか「歪曲収差」とかは、デジタルで補正できるようにセンサーとレンズの最適化を行っとるやろうけど、Lytroの最適化は1つ先の次元に行っとる気がする。<br />
<br />
<br />
<h3 style="text-align: left;">
2. 技術をもって本気で新しい表現を模索している</h3>
Lytroは新しい表現を探している。ホームページのムービーでも言っとるけど、Light Fieldカメラで一瞬を切り取るとどういう表現ができるか模索している。最初に出した製品はセンサーも小さく使い勝手も悪かったが、今回のIllumはクリエイターをターゲットにしている。<span style="color: red;">Light Fieldカメラが生き残るかどうかは、新しい表現方法が見つかるかどうかにかかっているといってもいい気がする。ボケ味を追求する日本のメーカーとは次元が違う。</span>(もちろんどっちがいいとかいう話は抜きにして)<br />
残念ながらぼくにはどういう表現がよいのか全くわかりません。<br />
<br />
やっぱりSilicon Valleyの企業は日本とは考え方が違う。<br />
ただLytroのIllumは使い勝手が悪いらしい。やっぱり手に馴染む機器の作り込みには時間がかかるようだ。<br />
<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com1tag:blogger.com,1999:blog-4766414733890203866.post-59188968190086699292014-11-03T22:50:00.001+09:002015-01-19T23:08:14.586+09:00Juliaで固定配列<div dir="ltr" style="text-align: left;" trbidi="on">
JuliaでArrayを使えばベクトルや行列計算できるけど、2とか3次元の点を表すのに可変長の配列を使うのはどうも無駄が多い気がする。自前で<br />
<div style="background: #f7f7f7; border-radius: 10px; border: 2px dotted #CCCCCC; padding: 10px;">
immutable Vec<br />
x::Float64<br />
y::Float64<br />
z::Float64<br />
end</div>
とかやったけど、演算子とかを定義するのが面倒。<br />
<br />
そんな問題を解決してくれるのが<span style="color: red;">ImuutableArrays</span>。<br />
Julia0.3にはすでに採用されているので、<br />
<div style="background: #f7f7f7; border-radius: 10px; border: 2px dotted #CCCCCC; padding: 10px;">
Using ImuutableArrays<br />
v3 = Vector3{Float64}(1, 2, 3)<br />
v2 = Vector2{Float64}(4, 5)</div>
<br />
と書ける。デフォルトでは読み込んだ時に4次元まで作ってくれて、演算子の定義とかもやってくれている。(これがメタプログラミングの力か…)<br />
しかも上で自分で定義したのより速い…<br />
<br />
行列は<br />
<div style="background: #f7f7f7; border-radius: 10px; border: 2px dotted #CCCCCC; padding: 10px;">
typealias Mat3d Matrix3x3{Float64}</div>
とすれば使える。<br />
<br />
<br />
<br />
<br />
<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-35190739740736879272014-11-02T23:51:00.005+09:002014-11-02T23:54:50.115+09:00IPython notebook でプレゼンテーション用スライドを作る<div dir="ltr" style="text-align: left;" trbidi="on">
プレゼンテーションを作るときに困ることがある。<br />
<br />
<ol style="text-align: left;">
<li>手軽に書きたい</li>
<li>Texで数式</li>
<li>図を入れたい。</li>
<li>プログラムを入れたい。</li>
<li>かっこ良くしたい。</li>
<li>PDFにしたい。</li>
<li>個人で使う分はフリーで済ませたい。</li>
</ol>
<div>
ざっと思いついたままに書くとこうなる。</div>
<div>
いろいろ探していたら、最近IPython notebookからスライドが作成できることが判明した!</div>
<div>
しかもreveal.jsを使っているのでいかにもギークっぽいかっこいいスライドが作れる。</div>
<div>
こんな感じのやつ→<a href="http://lab.hakim.se/reveal-js/#/">http://lab.hakim.se/reveal-js/#/</a></div>
<div>
<br /></div>
<div>
まえからreveal.jsは使ってみたかったけど、jsとかHTMLとかよく分からんし、<a href="https://slides.com/" target="_blank">Web上で使えるやつ</a>もあったけど、PDFにするのに有料会員にならんといけんとかで挫折した。</div>
<div>
<br />
<br />
<h2 style="text-align: left;">
スライドの作り方</h2>
<h3 style="text-align: left;">
1. IPython Notebookをインストールする。</h3>
IPython NotebookはAnacondaをつかってインストールするのがおすすめ。一式全部はいるし、アップデートも<br />
<div style="background: #f7f7f7; border-radius: 10px; border: 2px dotted #CCCCCC; padding: 10px;">
conda update anaconda</div>
で完了する。<br />
<br />
<h3 style="text-align: left;">
2. 新しいNotebookを開いて設定する。</h3>
New Notebookを作成したら以下の図のところをSlideに設定する。<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjngnebzVBVromYE7GsB7Bqrkd61BHdxa3QXpnlWZIaReTePVODdoq6P9wFDTQMItLFwS-2YKCCthiRq88KpwLS9ML5lB1gCmeS6LWYeZ08exqCzfVAvjiihHlaYHysP_ObaJypA67mDU/s1600/Untitled1+2014-11-02+23-28-26.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjngnebzVBVromYE7GsB7Bqrkd61BHdxa3QXpnlWZIaReTePVODdoq6P9wFDTQMItLFwS-2YKCCthiRq88KpwLS9ML5lB1gCmeS6LWYeZ08exqCzfVAvjiihHlaYHysP_ObaJypA67mDU/s1600/Untitled1+2014-11-02+23-28-26.png" height="118" width="400" /></a></div>
<br />
右のタブはSlide、Sub-Slideとか選べる。これはreveal.jsでいう横にスライドか、縦にスライドかを指定できる。いろいろ試してみて。<br />
<br />
<h3 style="text-align: left;">
3. スライドへ変換する</h3>
</div>
<div>
IPython Notebookの編集が終わったら、NotebookからSlideへの変換をする。スライドの作成は<a href="http://www.damian.oquanta.info/posts/make-your-slides-with-ipython.html" target="_blank">ここ</a>にあるとおりのコマンドを実行するとできる。--postでHTMLサーバーを指定する必要があるらしい。<br />
<div style="text-align: left;">
<br /></div>
<div style="background: #f7f7f7; border-radius: 10px; border: 2px dotted #CCCCCC; padding: 10px;">
ipython nbconvert notebook.ipynb --to slides --post serve</div>
<br />
IPythonのバージョンが古いと文字がものすごく小さいので最新版にアップデートするのをおすすめ。<br />
<b>※Chromeで開くと数式がレンダリングされなかったのでmathjaxのプラグインをいれた。</b><br />
ちなみにPDFにするときは<br />
<div style="background: #f7f7f7; border-radius: 10px; border: 2px dotted #CCCCCC; padding: 10px;">
ipython nbconvert notebook.ipynb --to latex --post pdf</div>
<br />
<br />
問題は記録用にPDFとかHTMLにしたらちょっとダサくなってしまうことか。<br />
まぁよい。<br />
<br />
<br /></div>
<div>
<div>
</div>
</div>
</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com4tag:blogger.com,1999:blog-4766414733890203866.post-38324205750259432212014-09-25T01:59:00.000+09:002014-09-25T02:10:19.470+09:00Opencv3.0 ALPHA Python版をMacにソースからインストールしたときのメモ<div dir="ltr" style="text-align: left;" trbidi="on">
Mac OSXでPythonはAnacondaでインストールしたPython3.4を使用<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiWW_gyfUI_r6dN2zL0REt38rqRSYkv7qJ4oocTtawBL7_CGNVpmD8-4ySzP9Dj9kjskGHiui6SUHVp-iZfooz9xcrpQHJnkKAhiioM6jyjMA8yeJyfgXkkzdXQVrlDFhj-_t0BEmO3nYo/s1600/temp.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiWW_gyfUI_r6dN2zL0REt38rqRSYkv7qJ4oocTtawBL7_CGNVpmD8-4ySzP9Dj9kjskGHiui6SUHVp-iZfooz9xcrpQHJnkKAhiioM6jyjMA8yeJyfgXkkzdXQVrlDFhj-_t0BEmO3nYo/s1600/temp.png" height="199" width="320" /></a></div>
<br />
<ol style="text-align: left;">
<li><a href="http://opencv.org/downloads.html" target="_blank">OpenCV for Linux/Mac</a>をダウンロードしてZipを解凍する</li>
<li>ダウンロードフォルダに移動して、作業用フォルダを作っておく。ここではbuildとした。</li>
<li>CMakeをダウンロードして実行。</li>
<li>ダウンロードしたOpenCVのフォルダとbuildフォルダを指定する。</li>
<li>Searchにpythonと入れてConfigを押す</li>
<li>Configを押すとPython3の設定がいろいろあるので上の図や参考サイトを参考にして埋める。赤い表示がなければとりあえずOK。</li>
<li>Generateを押す。</li>
<li>buildフォルダで
<div style="background: #f7f7f7; border-radius: 10px; border: 2px dotted #CCCCCC; padding: 10px;">
make -j8</div>
を実行。失敗するとエラーが出る。Python.hがないとか言われたけど、6の設定したPATHを見直すとうまくいった。</li>
<li>sudo make install</li>
<li>pythonを立ち上げてimport cv2</li>
<li>今回は「 Library not loaded: libpython3.4m.dylib」とか言われたのでシンボリックリンクを貼った。</li>
<div style="background: #f7f7f7; border-radius: 10px; border: 2px dotted #CCCCCC; padding: 10px;">
sudo ln -/Users/name/anaconda/lib/libpython3.4m.dylib /usr/local/lib</div>
</ol>
<div>
<br /></div>
<div>
これでとりあえずopencvが動いた。</div>
<div>
<br /></div>
<div>
参考サイト</div>
<ol>
<li><a href="http://luigolas.com/blog/2014/09/15/install-opencv3-with-python-3-mac-osx/">http://luigolas.com/blog/2014/09/15/install-opencv3-with-python-3-mac-osx/</a></li>
<li><a href="https://stackoverflow.com/questions/20953273/install-opencv-for-python-3-3">https://stackoverflow.com/questions/20953273/install-opencv-for-python-3-3</a></li>
</ol>
</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-27138110098026150942014-09-21T19:20:00.003+09:002014-09-21T19:20:54.112+09:00Light fieldの原理をLytroがじきじきに解説<div dir="ltr" style="text-align: left;" trbidi="on">
久しぶりにLytroのblogをのぞいたら、Lytroが論文を公表していた。<br />
まだ目を全部に目を通してないけど、見たことない頭があったりした。<br />
<a href="http://blog.lytro.com/post/97168596625/modeling-the-light-field-camera">http://blog.lytro.com/post/97168596625/modeling-the-light-field-camera</a><br />
<br />
最近はGoogle Cameraでもボケを作れるので静止画でのLytroの優位性はなくなってきている気が。動画だと処理が大変そうだしなー。</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-19805521322133580332014-09-19T19:58:00.001+09:002014-09-19T19:58:24.293+09:00Microsoft OCR Library がリリースされた。<div dir="ltr" style="text-align: left;" trbidi="on">
<a href="https://www.nuget.org/packages/Microsoft.Windows.Ocr/">https://www.nuget.org/packages/Microsoft.Windows.Ocr/</a><br />
<br />
WindowsRuntimeってWinRTのことで、Windows Metroのアプリをつくるためのものか。<br />
いつか使ってみたいけど、言語がどれもわからん。<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-8570797664394241132014-09-08T23:42:00.000+09:002014-09-12T00:34:56.630+09:00Juliaで簡単な回帰<div dir="ltr" style="text-align: left;" trbidi="on">
Juliaの練習がてら簡単な線形回帰を行った。<br />
<br />
以下のことが分かる練習にはちょうど良いシンプルなコードだ。<br />
<br />
<ul style="text-align: left;">
<li>乱数をどう発生させるか</li>
<li>計画行列をどう作るか</li>
<li>Gadflyを使って2つのプロットを重ねるにはどうするか</li>
<li>Gadflyを使って軸の設定をどうするか</li>
</ul>
<div>
これらがわかれば自作する当面の計算にはこまらない気がする…</div>
<div>
<a href="http://julia.readthedocs.org/en/latest/stdlib/linalg/#Base.linreg" target="_blank">linreg()</a>はバイアス項を勝手に計算してくれるみたいなので、入れなくてよいみたいだ。</div>
<div>
<br /></div>
<div>
<pre style="background: #fdf6e3; color: #586e75;"><span style="color: #93a1a1;"># #Regression Example</span>
using Gadfly
using Distributions
N <span style="color: #859900;">=</span> <span style="color: #d33682;">15</span>
xdata <span style="color: #859900;">=</span> (rand(N) <span style="color: #859900;">-</span> <span style="color: #d33682;">0.5</span>) <span style="color: #859900;">*</span> <span style="color: #d33682;">2</span> <span style="color: #859900;">*</span> pi
σ² <span style="color: #859900;">=</span> <span style="color: #d33682;">0.</span><span style="color: #d33682;">5</span>
function noisy(x, r=rand(Normal(<span style="color: #d33682;">0</span>, σ²))
return sin(x) + r
end
<span style="color: #93a1a1;">#X results in Nx9 array</span>
t <span style="color: #859900;">=</span> noisy(xdata);
X <span style="color: #859900;">=</span> Float64[x<span style="color: #859900;">^</span>i <span style="color: #859900;">for</span> x<span style="color: #859900;">=</span>xdata, i<span style="color: #859900;">=</span>[<span style="color: #d33682;">1</span>:<span style="color: #d33682;">9</span>]] <span style="color: #93a1a1;">#bias is not needed for linreg(). it will interpolate automatically</span>
<span style="color: #93a1a1;">#regression by linreg</span>
coeffs <span style="color: #859900;">=</span> linreg(X, t)
<span style="color: #93a1a1;">#plotting</span>
f <span style="color: #859900;">=</span> x <span style="color: #859900;">-</span><span style="color: #859900;">></span> coeffs[<span style="color: #d33682;">1</span>] <span style="color: #859900;">+</span> <span style="color: #268bd2;">sum</span>([coeffs[i]<span style="color: #859900;">*</span>x<span style="color: #859900;">^</span>(i<span style="color: #859900;">-</span><span style="color: #d33682;">1</span>) <span style="color: #859900;">for</span> i <span style="color: #859900;">=</span> [<span style="color: #d33682;">2</span>:<span style="color: #d33682;">10</span>]])
x_reg <span style="color: #859900;">=</span> linspace(<span style="color: #859900;">-</span>pi, pi, <span style="color: #d33682;">100</span>)
y_reg <span style="color: #859900;">=</span> [f(xi) <span style="color: #859900;">for</span> xi<span style="color: #859900;">=</span>x_reg]
plot(layer(x<span style="color: #859900;">=</span>x_reg, y<span style="color: #859900;">=</span>y_reg, Geom.line), layer(x<span style="color: #859900;">=</span>xdata, y<span style="color: #859900;">=</span>t, Geom.point),
Coord.Cartesian(xmin<span style="color: #859900;">=</span><span style="color: #859900;">-</span>pi, xmax<span style="color: #859900;">=</span>pi, ymin<span style="color: #859900;">=</span><span style="color: #859900;">-</span><span style="color: #d33682;">1.5</span>, ymax<span style="color: #859900;">=</span><span style="color: #d33682;">1.5</span>))
</pre>
<br /></div>
</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-135163383518188612014-08-22T00:09:00.002+09:002014-08-22T00:09:29.232+09:00グラフ読み取り<div dir="ltr" style="text-align: left;" trbidi="on">
PyQtでグラフを読み取れるアプリを作った。<br />
グラフの画像をCtrl+Vでペーストし、その座標をクリックすると画素値からグラフの値を教えてくれる。<br />
<br />
最初に右クリックで左下と右上の座標を教える必要がる。<br />
なんでこんなんがいるかって、世の中紙でデータを保存する人がまだいるということである。<br />
<br />
(OpenCVは使っていない)<br />
<br />
<pre style="background: #fdf6e3; color: #586e75;"><span style="color: #bd3800;">from</span> PyQt4 <span style="color: #bd3800;">import</span> QtGui, QtCore, QtCore
<span style="color: #bd3800;">import</span> sys
<span style="color: #bd3800;">import</span> cv2
<span style="color: #bd3800;">import</span> numpy <span style="color: #859900;">as</span> np
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">convert_array_to_qimg</span>(cv_img):
<span style="color: #269186;"><span style="color: #c60000;">"""</span>
this is based on https://stackoverflow.com/questions/18676888/how-to-configure-color-when-convert-numpy-array-to-qimage
<span style="color: #c60000;">"""</span></span>
height, width, bytesPerComponent <span style="color: #859900;">=</span> cv_img.shape
bytesPerLine <span style="color: #859900;">=</span> bytesPerComponent <span style="color: #859900;">*</span> width;
<span style="color: #93a1a1;">#cv2.cvtColor(cv_img, cv2.CV_BGR2RGB)</span>
<span style="color: #859900;">return</span> QtGui.QImage(cv_img.data, width, height, bytesPerLine, QtGui.QImage.Format_RGB888)
<span style="color: #748b00;">class</span> <span style="color: #268bd2;">GraphCoordinate</span>:
<span style="color: #748b00;">def</span> <span style="color: #268bd2;"><span style="color: #268bd2;">__init__</span></span>(self):
<span style="color: #859900;">pass</span>
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">setBottomPixelCoordinate</span>(self, x1, y1):
<span style="color: #269186;"><span style="color: #c60000;">"""</span>
x1 and y1 are pixel
<span style="color: #c60000;">"""</span></span>
<span style="color: #268bd2;">self</span>.x1 <span style="color: #859900;">=</span> x1
<span style="color: #268bd2;">self</span>.y1 <span style="color: #859900;">=</span> y1
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">setBottomGraphCoordinate</span>(self, X1, Y1):
<span style="color: #268bd2;">self</span>.X1 <span style="color: #859900;">=</span> X1
<span style="color: #268bd2;">self</span>.Y1 <span style="color: #859900;">=</span> Y1
<span style="color: #859900;">print</span> <span style="color: #269186;"><span style="color: #c60000;">"</span>(X1, Y1) = (<span style="color: #cb4b16;">%.2f</span>, <span style="color: #cb4b16;">%.2f</span>)<span style="color: #c60000;">"</span></span> <span style="color: #859900;">%</span> (<span style="color: #268bd2;">self</span>.X1, <span style="color: #268bd2;">self</span>.Y1)
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">setTopPixelCoordinate</span>(self, x2, y2):
<span style="color: #268bd2;">self</span>.x2 <span style="color: #859900;">=</span> x2
<span style="color: #268bd2;">self</span>.y2 <span style="color: #859900;">=</span> y2
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">setTopGraphCoordinate</span>(self, X2, Y2):
<span style="color: #268bd2;">self</span>.X2 <span style="color: #859900;">=</span> X2
<span style="color: #268bd2;">self</span>.Y2 <span style="color: #859900;">=</span> Y2
<span style="color: #859900;">print</span> <span style="color: #269186;"><span style="color: #c60000;">"</span>(X2, Y2) = (<span style="color: #cb4b16;">%.2f</span>, <span style="color: #cb4b16;">%.2f</span>)<span style="color: #c60000;">"</span></span> <span style="color: #859900;">%</span> (<span style="color: #268bd2;">self</span>.X2, <span style="color: #268bd2;">self</span>.Y2)
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">getGraphCoordinate</span>(self, x, y):
w <span style="color: #859900;">=</span> x <span style="color: #859900;">-</span> <span style="color: #268bd2;">self</span>.x1
h <span style="color: #859900;">=</span> <span style="color: #268bd2;">self</span>.y1 <span style="color: #859900;">-</span> y
W <span style="color: #859900;">=</span> <span style="color: #268bd2;">self</span>.x2 <span style="color: #859900;">-</span> <span style="color: #268bd2;">self</span>.x1
H <span style="color: #859900;">=</span> <span style="color: #268bd2;">self</span>.y1 <span style="color: #859900;">-</span> <span style="color: #268bd2;">self</span>.y2
X <span style="color: #859900;">=</span> <span style="color: #268bd2;">self</span>.X1 <span style="color: #859900;">+</span> w <span style="color: #859900;">*</span> (<span style="color: #268bd2;">self</span>.X2 <span style="color: #859900;">-</span> <span style="color: #268bd2;">self</span>.X1) <span style="color: #859900;">/</span> <span style="color: #859900;">float</span>(W)
Y <span style="color: #859900;">=</span> <span style="color: #268bd2;">self</span>.Y1 <span style="color: #859900;">+</span> h <span style="color: #859900;">*</span> (<span style="color: #268bd2;">self</span>.Y2 <span style="color: #859900;">-</span> <span style="color: #268bd2;">self</span>.Y1) <span style="color: #859900;">/</span> <span style="color: #859900;">float</span>(H)
<span style="color: #859900;">return</span> X, Y
<span style="color: #748b00;">class</span> <span style="color: #268bd2;">MainWindow</span>(QtGui.QMainWindow):
<span style="color: #748b00;">def</span> <span style="color: #268bd2;"><span style="color: #268bd2;">__init__</span></span>(self):
<span style="color: #859900;">super</span>(MainWindow, <span style="color: #268bd2;">self</span>).__init__()
<span style="color: #268bd2;">self</span>.image_frame <span style="color: #859900;">=</span> QtGui.QLabel()
fname <span style="color: #859900;">=</span> <span style="color: #269186;"><span style="color: #073642; font-weight: 700;">r</span><span style="color: #c60000;">"</span>sample08<span style="color: #cb4b16;">.</span>png<span style="color: #c60000;">"</span></span>
img <span style="color: #859900;">=</span> cv2.imread(fname)
<span style="color: #93a1a1;">#img = QtGui.QImage(fname)</span>
qimg <span style="color: #859900;">=</span> convert_array_to_qimg(img)
<span style="color: #268bd2;">self</span>.image_frame.setPixmap(QtGui.QPixmap.fromImage(qimg))
<span style="color: #268bd2;">self</span>.clip <span style="color: #859900;">=</span> QtGui.QApplication.clipboard()
<span style="color: #268bd2;">self</span>.setCentralWidget(<span style="color: #268bd2;">self</span>.image_frame)
<span style="color: #268bd2;">self</span>.move(<span style="color: #d33682;">30</span>,<span style="color: #d33682;">30</span>)
<span style="color: #268bd2;">self</span>.counter <span style="color: #859900;">=</span> <span style="color: #d33682;">0</span>
<span style="color: #268bd2;">self</span>.coord <span style="color: #859900;">=</span> GraphCoordinate()
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">keyPressEvent</span>(self, e):
<span style="color: #859900;">if</span> (e.modifiers() <span style="color: #859900;">&</span> QtCore.Qt.ControlModifier):
<span style="color: #93a1a1;">#selected = self.table.selectedRanges()</span>
<span style="color: #859900;">if</span> e.key() <span style="color: #859900;">==</span> QtCore.Qt.Key_V:<span style="color: #93a1a1;">#past</span>
qimg <span style="color: #859900;">=</span> <span style="color: #268bd2;">self</span>.clip.image()
<span style="color: #268bd2;">self</span>.displayImage(qimg)
<span style="color: #859900;">elif</span> e.key() <span style="color: #859900;">==</span> QtCore.Qt.Key_C: <span style="color: #93a1a1;">#copy</span>
<span style="color: #859900;">pass</span>
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">mousePressEvent</span>(self, e):
<span style="color: #93a1a1;">#print e.x(), e.y()</span>
<span style="color: #859900;">if</span> e.button() <span style="color: #859900;">==</span> QtCore.Qt.RightButton:
<span style="color: #859900;">if</span> <span style="color: #268bd2;">self</span>.counter <span style="color: #859900;">%</span> <span style="color: #d33682;">2</span> <span style="color: #859900;">==</span> <span style="color: #d33682;">0</span>:
text, ok <span style="color: #859900;">=</span> QtGui.QInputDialog.getText(<span style="color: #268bd2;">self</span>, <span style="color: #269186;"><span style="color: #c60000;">'</span>Input Dialog<span style="color: #c60000;">'</span></span>,
<span style="color: #269186;"><span style="color: #c60000;">'</span>Enter X1, Y1 (comma separated):<span style="color: #c60000;">'</span></span>)
X1, Y1 <span style="color: #859900;">=</span> [<span style="color: #859900;">float</span>(v) <span style="color: #859900;">for</span> v <span style="color: #859900;">in</span> text.split(<span style="color: #269186;"><span style="color: #c60000;">'</span>,<span style="color: #c60000;">'</span></span>)]
<span style="color: #859900;">if</span> ok:
<span style="color: #268bd2;">self</span>.counter <span style="color: #859900;">+=</span> <span style="color: #d33682;">1</span>
<span style="color: #268bd2;">self</span>.coord.setBottomPixelCoordinate(e.x(), e.y())
<span style="color: #268bd2;">self</span>.coord.setBottomGraphCoordinate(X1, Y1)
<span style="color: #859900;">else</span>:
text, ok <span style="color: #859900;">=</span> QtGui.QInputDialog.getText(<span style="color: #268bd2;">self</span>, <span style="color: #269186;"><span style="color: #c60000;">'</span>Input Dialog<span style="color: #c60000;">'</span></span>,
<span style="color: #269186;"><span style="color: #c60000;">'</span>Enter X2, Y2 (comma separated):<span style="color: #c60000;">'</span></span>)
X2, Y2 <span style="color: #859900;">=</span> [<span style="color: #859900;">float</span>(v) <span style="color: #859900;">for</span> v <span style="color: #859900;">in</span> text.split(<span style="color: #269186;"><span style="color: #c60000;">'</span>,<span style="color: #c60000;">'</span></span>)]
<span style="color: #859900;">if</span> ok:
<span style="color: #268bd2;">self</span>.coord.setTopPixelCoordinate(e.x(), e.y())
<span style="color: #268bd2;">self</span>.coord.setTopGraphCoordinate(X2, Y2)
<span style="color: #268bd2;">self</span>.counter <span style="color: #859900;">=</span> <span style="color: #d33682;">0</span>
<span style="color: #859900;">elif</span> e.button() <span style="color: #859900;">==</span> QtCore.Qt.LeftButton:
<span style="color: #859900;">try</span>:
x, y <span style="color: #859900;">=</span> <span style="color: #268bd2;">self</span>.coord.getGraphCoordinate(e.x(), e.y())
<span style="color: #859900;">print</span> <span style="color: #269186;"><span style="color: #c60000;">"</span><span style="color: #cb4b16;">%.5f</span> <span style="color: #cb4b16;">%.5f</span><span style="color: #c60000;">"</span></span> <span style="color: #859900;">%</span> (x, y)
<span style="color: #859900;">except</span> <span style="color: #a57800;">AttributeError</span> <span style="color: #859900;">as</span> e:
<span style="color: #859900;">print</span> <span style="color: #269186;"><span style="color: #c60000;">"</span>please set the graph coordinate by right clicks.<span style="color: #c60000;">"</span></span>
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">displayImage</span>(self, qimg):
<span style="color: #268bd2;">self</span>.image_frame.setPixmap(QtGui.QPixmap.fromImage(qimg))
<span style="color: #748b00;">def</span> <span style="color: #268bd2;">main</span>(args):
app <span style="color: #859900;">=</span> QtGui.QApplication(sys.argv)
form <span style="color: #859900;">=</span> MainWindow()
form.show()
sys.exit(app.exec_())
<span style="color: #859900;">if</span> __name__ <span style="color: #859900;">==</span> <span style="color: #269186;"><span style="color: #c60000;">"</span>__main__<span style="color: #c60000;">"</span></span>:
main(sys.argv)
</pre>
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-8833318938829913832014-08-10T17:58:00.001+09:002014-11-16T23:40:34.308+09:00今さら聞けないRichardson-Lucy入門<div dir="ltr" style="text-align: left;" trbidi="on">
デコンボリューションがとても魅力的だ。これを使えばレンズの収差やボケ、手ブレによる画像劣化が除去できるという夢の技術だ。<br />
<div>
<br /></div>
<div>
いろんな論文を読むとかならずといっていいほど以下の2つの手法が挙げられている。</div>
<div>
<b>1.Wiener Filter</b></div>
<div>
Wikipediaに導出方法が載っている。</div>
<div>
<br /></div>
<div>
<b>2.Richardson-Lucy Deconvolution</b></div>
<div>
こちらもWikipediaに載っている。Wikipediaを読んでもわからなかったので、勉強した内容を少し解説してみる。<br />
<br />
わかりやすいように1次元で話を進める。<br />
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
</div>
<div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEirQgi5xIHsBf0Cq5pw8SGcJBjEfBGD1am675ZONZryaqP_yV6dt3fAWi-VrBZYcQfaxipitqoCnnXiDErQQPHxYz5hHQrIz10cje9B9qfrk06R_DCnUo7RXPnAfIn4RXne7WH_R1jMi84/s1600/CodeCogsEqn+(2).png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEirQgi5xIHsBf0Cq5pw8SGcJBjEfBGD1am675ZONZryaqP_yV6dt3fAWi-VrBZYcQfaxipitqoCnnXiDErQQPHxYz5hHQrIz10cje9B9qfrk06R_DCnUo7RXPnAfIn4RXne7WH_R1jMi84/s1600/CodeCogsEqn+(2).png" /></a></div>
をピクセルiの真の値(全くボケていないときの値)<br />
<br /></div>
<div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEglKvffjdFXfrtmFnOiL22RTrUfhSiJOu1V6CSRDMoFwd4yy2I48UGZw_vCOrP9mYKwBR7fBQsZ-vUqCl8sVBMjrZkRX-kNiyopwXIkaKibPeTsdh3DKEpHI4xdeR1sSAeMZLJGQpgbLQ8/s1600/CodeCogsEqn+(3).png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEglKvffjdFXfrtmFnOiL22RTrUfhSiJOu1V6CSRDMoFwd4yy2I48UGZw_vCOrP9mYKwBR7fBQsZ-vUqCl8sVBMjrZkRX-kNiyopwXIkaKibPeTsdh3DKEpHI4xdeR1sSAeMZLJGQpgbLQ8/s1600/CodeCogsEqn+(3).png" /></a></div>
をピクセルiで観測された値<br />
<br />
とする。全くボケていないときはこの2つが一致する。でも今回はPSF(Point Spread Function)によってボケる。ボケがあるとピクセルjへ行く光は広がりをもってしまうので、ピクセルiへ少しおすそ分けしてしまう。<br />
本来jへいくはずの光のボケでのih量を<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgr63kwYZavIS4ZAAtYwSSDu97HPi0XJVfaGzJrXh4isJTJybk2mCMzxBSWtoOItwxm2P_kaB9F5bCGp1-fWsKAhNQN9GIjiscGlG5-BNbfghyQITJ-Ck35F9-CND2jZ23yZB5sgesuSi8/s1600/CodeCogsEqn+(6).png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgr63kwYZavIS4ZAAtYwSSDu97HPi0XJVfaGzJrXh4isJTJybk2mCMzxBSWtoOItwxm2P_kaB9F5bCGp1-fWsKAhNQN9GIjiscGlG5-BNbfghyQITJ-Ck35F9-CND2jZ23yZB5sgesuSi8/s1600/CodeCogsEqn+(6).png" /></a></div>
<br />
<div class="" style="clear: both; text-align: left;">
とする。Richardson-Lucyではこの<span style="color: red;">p(i,j)は既知</span>だとするイメージは下の図。</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjP85aDdw7FKKbSHmnK8tpL0wUeUD91ReaHYpoiTnFba5HgJpYEHe1BNN6T-ig44HZJblvvoO2L0nwzHZesXF-OxiYqEuQ2DiwVspzmoQGFQ1JfxelY1n0-W9sZgbrMbLUeATJ3S87MAMc/s1600/2014-11-161.jpeg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjP85aDdw7FKKbSHmnK8tpL0wUeUD91ReaHYpoiTnFba5HgJpYEHe1BNN6T-ig44HZJblvvoO2L0nwzHZesXF-OxiYqEuQ2DiwVspzmoQGFQ1JfxelY1n0-W9sZgbrMbLUeATJ3S87MAMc/s1600/2014-11-161.jpeg" height="240" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
</div>
<div>
<br />
ここでZ(i,j)なるものを導入する。Z(i,j)は期待値がP(i,j)λjのポアッソン分布に従いう確率変数。<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj7iCoktEIVDICcs7ql7ojpyoLBKYc6AvDSpCyoKKK2wpLFBIOYW-Iu35ERr_fQ2POjkkWWewMYoHtB0rW6FGNg2_yNjEj_ykK0J2p2FaGaBNTewpFV4LbEQ8L458pQjGltBYgf953YOi0/s1600/png.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj7iCoktEIVDICcs7ql7ojpyoLBKYc6AvDSpCyoKKK2wpLFBIOYW-Iu35ERr_fQ2POjkkWWewMYoHtB0rW6FGNg2_yNjEj_ykK0J2p2FaGaBNTewpFV4LbEQ8L458pQjGltBYgf953YOi0/s1600/png.png" height="28" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
これはピクセルiに届く光が確率的に変化し、期待値はp(i, j)λjということを意味する。</div>
Wikipekiaにこの量は出てこないけど、元論文などを読むと出てくる。<br />
もともとデコンボリューションはCTとかで研究されていて、CTは原子の崩壊(ガンマ線?)を捉えるので、観測される値が確率的になる、という経緯がある。<br />
(一様な確率で起こる事象がある期間で観測される回数はポアッソン分布になる。詳しくはWikipediaをみてください。)<br />
<br />
さて、点が1点だけなら復元も簡単かもしれんけど、周りの点の広がりが重なって観測されるので問題が難しくなる。例えば下の図のように重なるのでピクセルiで値は, 周辺のピクセルからの寄与を足しあわせたのものになる。<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjGZjy7DAoO18oiCV7rbuj2q3WQK6ASzUivb6pUb_HkorxeiZRvkWCdemMuBSUjCCQpQXT9BVqT_CKEhoajLbHK6V4u8Xtp34Cray9cq8-x0RtJIQhNs0On1pj4qK6bF0xyDW2F773meSU/s1600/png+(3).png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjGZjy7DAoO18oiCV7rbuj2q3WQK6ASzUivb6pUb_HkorxeiZRvkWCdemMuBSUjCCQpQXT9BVqT_CKEhoajLbHK6V4u8Xtp34Cray9cq8-x0RtJIQhNs0On1pj4qK6bF0xyDW2F773meSU/s1600/png+(3).png" height="70" width="200" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhmwzMZAf-vDOXZlN1jRrdmfr1r8j-ZCFH8Cr2vZvJ6z_lDjgmtSTruGlicb-wnXR0Mzv5toEfME0p7KqL_P7EZng0z5rgcepIQvDY0SkY6w6M5-2sc8ky5aj0cR9XdGoBZrKHMc493Eow/s1600/2014-11-161+(1).jpeg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhmwzMZAf-vDOXZlN1jRrdmfr1r8j-ZCFH8Cr2vZvJ6z_lDjgmtSTruGlicb-wnXR0Mzv5toEfME0p7KqL_P7EZng0z5rgcepIQvDY0SkY6w6M5-2sc8ky5aj0cR9XdGoBZrKHMc493Eow/s1600/2014-11-161+(1).jpeg" height="240" width="320" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
ここでzは確率的な値なので、観測結果と期待値から予測することで話を前に進めよう。<br />
ピクセルiで観測された値をつかうと、その内ピクセルjからきた値は、割合で考えると下のような式になる。右辺の分母がiへくる光の総和で、分子がその内jからきたものを表す。<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg4M3zwoVqnUS3zJRINstAZTp-HpkvinM_kVtAMX_-vHO5Kzuj8Ojg5i28sV9CWK-qMKwsFP8VysD_kB6Gt3cAEfC8H4ldOmfLBFQw94DfchUz6VkgIrtcmwrsiIYOfGxBWOXV-YjF5ay0/s1600/png+(4).png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg4M3zwoVqnUS3zJRINstAZTp-HpkvinM_kVtAMX_-vHO5Kzuj8Ojg5i28sV9CWK-qMKwsFP8VysD_kB6Gt3cAEfC8H4ldOmfLBFQw94DfchUz6VkgIrtcmwrsiIYOfGxBWOXV-YjF5ay0/s1600/png+(4).png" height="78" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiaLHWxmB9PwGN1TpcIsA1x9oivO4LJ1pHVKiQCYV7DM3vJJ_rc_gaX79Y9Py9vL-V_jsL6AcqmWNn_cI3TdrxdSgJZWHaUZzeFoYA3qkevbJxMYck-TNWymEGtR4ye951Y4l6o9qqUnPs/s1600/2014-11-161+(2).jpeg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiaLHWxmB9PwGN1TpcIsA1x9oivO4LJ1pHVKiQCYV7DM3vJJ_rc_gaX79Y9Py9vL-V_jsL6AcqmWNn_cI3TdrxdSgJZWHaUZzeFoYA3qkevbJxMYck-TNWymEGtR4ye951Y4l6o9qqUnPs/s1600/2014-11-161+(2).jpeg" height="320" width="240" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
zはもともとjにあったものが周辺のピクセルiに散らばったものとかんがえることができるので、集めてやると真の値λjが分かる。真のzはわからないので推定値を使うとλを推定できてる。<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiYuKUXrd0p4_ubPHfcV4BaiswxvlyI426bR1A6YBqZWCRwpdSgfcyzXa5TFnyp_-DM8kh4k5KErIha18I9xdq7Rlxjub493zoHMrVg0a9JvoCAZZTAlCDRkcsK2l8wUZXrnzvy_WiZ6os/s1600/2014-11-161+(3).jpeg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiYuKUXrd0p4_ubPHfcV4BaiswxvlyI426bR1A6YBqZWCRwpdSgfcyzXa5TFnyp_-DM8kh4k5KErIha18I9xdq7Rlxjub493zoHMrVg0a9JvoCAZZTAlCDRkcsK2l8wUZXrnzvy_WiZ6os/s1600/2014-11-161+(3).jpeg" height="240" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
この式にさっきのzの推定値を代入すると<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgO5YLiUspFWLz6guGfbgnOLHVrwbYEcfJcOEXGECdvxONeSLAKrbGzTF9yLnfMqgNG4LwqfMWI6Mwul-_DZ-i4OgiLseSwNy6psqJOHkNtt1lYU2ZuQFB1UaYp_e_F6ELUI5xQxfLfpn4/s1600/png+(6).png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgO5YLiUspFWLz6guGfbgnOLHVrwbYEcfJcOEXGECdvxONeSLAKrbGzTF9yLnfMqgNG4LwqfMWI6Mwul-_DZ-i4OgiLseSwNy6psqJOHkNtt1lYU2ZuQFB1UaYp_e_F6ELUI5xQxfLfpn4/s1600/png+(6).png" height="55" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
このλjは真の値なので、知り得ない。そこで前回推定した値を使って逐次的に求めることにすれば</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhhjT87RKVduMvh9F8dux5MLWq5mtLLRR3ktAiZKEPuy0Ifmd-WVxlgzaVI8QQS7wjC72TpxaaXf-1h8rB8g1IKfZ1UmaUev3Oxi8nHlBnpEEnGMzTtcwbq_X5T7J6gl2dU7u1mWHilyhI/s1600/png+(7).png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhhjT87RKVduMvh9F8dux5MLWq5mtLLRR3ktAiZKEPuy0Ifmd-WVxlgzaVI8QQS7wjC72TpxaaXf-1h8rB8g1IKfZ1UmaUev3Oxi8nHlBnpEEnGMzTtcwbq_X5T7J6gl2dU7u1mWHilyhI/s1600/png+(7).png" height="66" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
というふうになる。<br />
<br />
実際に実装するときはWikipediaのようにコンボリューションの関数を使って書いたほうが楽だった。<br />
<br />
<br />
(2014/11/15 式間違ってましたので修正しました。前に書いてたやつではPSFが対称じゃないとうまくいかなかったorz)<br />
<br />
参考<br />
<a href="https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution">https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution</a><br />
<a href="http://people.csail.mit.edu/dgreensp/may/deconvlucy.html">http://people.csail.mit.edu/dgreensp/may/deconvlucy.html</a><br />
<br />
<br />
<br /></div>
</div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0tag:blogger.com,1999:blog-4766414733890203866.post-62942158842084940992014-08-10T15:59:00.001+09:002014-08-10T15:59:05.219+09:00Juliaで自作VectorとArray{Float64,1}の速度比較<div dir="ltr" style="text-align: left;" trbidi="on">
Juliaの勉強を始めた。<br />
Juliaは数値計算の関数を大量に有しているけど、可変長のArrayをもとにしていて、3x3の行列や3x1ベクトルの単純なものだと自分で定義したほうが高速なのだろうか?という疑問が湧いたので実験してみた。<br />
もしかしたらJuliaは裏でものすごいことをやっていてるかもしれんし。<br />
<br />
3x1のベクトルVecを以下のように定義した。<br />
<pre style="background: #fdf6e3; color: #586e75;"><span style="color: #cb4b16;">immutable</span> <span style="color: #268bd2;">Vec</span>
<span style="color: #cb4b16;">x</span>::<span style="color: #268bd2;">Float64</span>
<span style="color: #cb4b16;">y</span>::<span style="color: #268bd2;">Float64</span>
<span style="color: #cb4b16;">z</span>::<span style="color: #268bd2;">Float64</span>
<span style="color: #859900;">end</span>
<span style="color: #268bd2;">dot</span>(<span style="color: #cb4b16;">a</span>::<span style="color: #268bd2;">Vec</span>, <span style="color: #cb4b16;">b</span>::<span style="color: #268bd2;">Vec</span>) <span style="color: #859900;">=</span> (<span style="color: #cb4b16;">a</span>.<span style="color: #cb4b16;">x</span><span style="color: #859900;">*</span><span style="color: #cb4b16;">b</span>.<span style="color: #cb4b16;">x</span> <span style="color: #859900;">+</span> <span style="color: #cb4b16;">a</span>.<span style="color: #cb4b16;">y</span><span style="color: #859900;">*</span><span style="color: #cb4b16;">b</span>.<span style="color: #cb4b16;">y</span> <span style="color: #859900;">+</span> <span style="color: #cb4b16;">a</span>.<span style="color: #cb4b16;">z</span><span style="color: #859900;">*</span><span style="color: #cb4b16;">b</span>.<span style="color: #cb4b16;">z</span>)
<span style="color: #268bd2;">cross</span>(<span style="color: #cb4b16;">a</span>::<span style="color: #268bd2;">Vec</span>, <span style="color: #cb4b16;">b</span>::<span style="color: #268bd2;">Vec</span>) <span style="color: #859900;">=</span> <span style="color: #268bd2;">Vec</span>(<span style="color: #cb4b16;">a</span>.<span style="color: #cb4b16;">y</span><span style="color: #859900;">*</span><span style="color: #cb4b16;">b</span>.<span style="color: #cb4b16;">z</span> <span style="color: #859900;">-</span> <span style="color: #cb4b16;">a</span>.<span style="color: #cb4b16;">z</span><span style="color: #859900;">*</span><span style="color: #cb4b16;">b</span>.<span style="color: #cb4b16;">y</span>, <span style="color: #cb4b16;">a</span>.<span style="color: #cb4b16;">z</span><span style="color: #859900;">*</span><span style="color: #cb4b16;">b</span>.<span style="color: #cb4b16;">x</span><span style="color: #859900;">-</span><span style="color: #cb4b16;">a</span>.<span style="color: #cb4b16;">x</span><span style="color: #859900;">*</span><span style="color: #cb4b16;">b</span>.<span style="color: #cb4b16;">z</span>, <span style="color: #cb4b16;">a</span>.<span style="color: #cb4b16;">x</span><span style="color: #859900;">*</span><span style="color: #cb4b16;">b</span>.<span style="color: #cb4b16;">y</span><span style="color: #859900;">-</span><span style="color: #cb4b16;">a</span>.<span style="color: #cb4b16;">y</span><span style="color: #859900;">*</span><span style="color: #cb4b16;">b</span>.<span style="color: #cb4b16;">x</span>)
</pre>
<br />
ただ外積やらなんやらを計算する関数を作った。<br />
<pre style="background: #fdf6e3; color: #586e75;"><span style="color: #cb4b16;">function</span> <span style="color: #268bd2;">cross_loop</span>(<span style="color: #268bd2;">N</span>::<span style="color: #268bd2;">Int</span>, <span style="color: #cb4b16;">a</span>::<span style="color: #268bd2;">Vec</span>, <span style="color: #cb4b16;">b</span>::<span style="color: #268bd2;">Vec</span>)
<span style="color: #cb4b16;">tmp</span> <span style="color: #859900;">=</span> <span style="color: #d33682;">0.0</span>
<span style="color: #cb4b16;">for</span> <span style="color: #cb4b16;">i</span> <span style="color: #cb4b16;">in</span> <span style="color: #d33682;">1</span>:<span style="color: #268bd2;">N</span>
<span style="color: #cb4b16;">v</span> <span style="color: #859900;">=</span> <span style="color: #268bd2;">cross</span>(<span style="color: #cb4b16;">a</span>,<span style="color: #cb4b16;">b</span>)
<span style="color: #cb4b16;">tmp</span> <span style="color: #859900;">+</span><span style="color: #859900;">=</span> <span style="color: #268bd2;">dot</span>(<span style="color: #cb4b16;">v</span>,<span style="color: #cb4b16;">v</span>)
<span style="color: #859900;">end</span>
<span style="color: #cb4b16;">tmp</span><span style="color: #859900;">/</span><span style="color: #268bd2;">N</span>
<span style="color: #859900;">end</span>
<span style="color: #cb4b16;">function</span> <span style="color: #268bd2;">cross_loop</span>(<span style="color: #268bd2;">N</span>::<span style="color: #268bd2;">Int</span>, <span style="color: #cb4b16;">a</span>::<span style="color: #268bd2;">Array</span>{<span style="color: #268bd2;">Float64</span>,<span style="color: #d33682;">1</span>}, <span style="color: #cb4b16;">b</span>::<span style="color: #268bd2;">Array</span>{<span style="color: #268bd2;">Float64</span>,<span style="color: #d33682;">1</span>})
<span style="color: #cb4b16;">tmp</span> <span style="color: #859900;">=</span> <span style="color: #d33682;">0</span>.
<span style="color: #cb4b16;">for</span> <span style="color: #cb4b16;">i</span> <span style="color: #cb4b16;">in</span> <span style="color: #d33682;">1</span>:<span style="color: #268bd2;">N</span>
<span style="color: #cb4b16;">v</span> <span style="color: #859900;">=</span> <span style="color: #268bd2;">Base</span>.<span style="color: #268bd2;">cross</span>(<span style="color: #cb4b16;">a</span>,<span style="color: #cb4b16;">b</span>)
<span style="color: #cb4b16;">tmp</span> <span style="color: #859900;">+</span><span style="color: #859900;">=</span> <span style="color: #268bd2;">Base</span>.<span style="color: #268bd2;">dot</span>(<span style="color: #cb4b16;">v</span>,<span style="color: #cb4b16;">v</span>)
<span style="color: #859900;">end</span>
<span style="color: #cb4b16;">tmp</span> <span style="color: #859900;">/</span> <span style="color: #268bd2;">N</span>
<span style="color: #859900;">end</span>
</pre>
Arrayの方はデフォルトで入っているBase.dotやcrossをつかった。<br />
<br />
N=10^7くらいにして計算<br />
速度は@timeで計測した。1回目はコンパイルに時間を取られるので2回目の時間を使用。<br />
<pre style="background: #fdf6e3; color: #586e75;"><span style="color: #268bd2;">N</span> <span style="color: #859900;">=</span> <span style="color: #d33682;">10</span>^<span style="color: #d33682;">7</span>
<span style="color: #cb4b16;">a</span> <span style="color: #859900;">=</span> <span style="color: #268bd2;">Vec</span>(<span style="color: #d33682;">1</span>,<span style="color: #d33682;">2</span>,<span style="color: #d33682;">3</span>)
<span style="color: #cb4b16;">b</span> <span style="color: #859900;">=</span> <span style="color: #268bd2;">Vec</span>(<span style="color: #d33682;">4</span>,<span style="color: #d33682;">5</span>,<span style="color: #d33682;">6</span>)
<span style="color: #cb4b16;">v</span> <span style="color: #859900;">=</span> [<span style="color: #d33682;">1</span>.,<span style="color: #d33682;">2</span>.,<span style="color: #d33682;">3</span>.]
<span style="color: #cb4b16;">u</span> <span style="color: #859900;">=</span> [<span style="color: #d33682;">4</span>.,<span style="color: #d33682;">5</span>.,<span style="color: #d33682;">6</span>.]
@<span style="color: #cb4b16;">time</span> <span style="color: #268bd2;">cross_loop</span>(<span style="color: #268bd2;">N</span>, <span style="color: #cb4b16;">a</span>, <span style="color: #cb4b16;">b</span>)
@<span style="color: #cb4b16;">time</span> <span style="color: #268bd2;">cross_loop</span>(<span style="color: #268bd2;">N</span>, <span style="color: #cb4b16;">u</span>, <span style="color: #cb4b16;">v</span>)
</pre>
<br />
結果はダントツで自作のほうが速い。<br />
<pre style="background: #fdf6e3; color: #586e75;"><span style="color: #cb4b16;">elapsed</span> <span style="color: #cb4b16;">time</span>: <span style="color: #d33682;">0.013432361</span> <span style="color: #268bd2;">seconds</span> (<span style="color: #d33682;">96</span> <span style="color: #cb4b16;">bytes</span> <span style="color: #cb4b16;">allocated</span>)
<span style="color: #cb4b16;">elapsed</span> <span style="color: #cb4b16;">time</span>: <span style="color: #d33682;">2.079223219</span> <span style="color: #268bd2;">seconds</span> (<span style="color: #d33682;">1280000096</span> <span style="color: #cb4b16;">bytes</span> <span style="color: #cb4b16;">allocated</span>, <span style="color: #d33682;">43.49</span><span style="color: #93a1a1;">% gc time)</span>
</pre>
<br />
<br />
ちなみにVecをimmutableでなくtypeで定義すると、それだけでかなり遅くなる。GCが働いているようだ。<br />
<pre style="background: #fdf6e3; color: #586e75;"><span style="color: #cb4b16;">elapsed</span> <span style="color: #cb4b16;">time</span>: <span style="color: #d33682;">0.419364716</span> <span style="color: #268bd2;">seconds</span> (<span style="color: #d33682;">320000096</span> <span style="color: #cb4b16;">bytes</span> <span style="color: #cb4b16;">allocated</span>, <span style="color: #d33682;">53.07</span><span style="color: #93a1a1;">% gc time)</span>
</pre>
<br />
<br />
結論:3次元に特化したベクトル計算では自分で定義してやったほうが速い。<br />
<br />
<br />
<br />
<br /></div>
anonymous optical engineerhttp://www.blogger.com/profile/03170825942566594155noreply@blogger.com0