Compressed Sensingを理解するために。

最近世の中ではCompressed Sensing(CS)が流行っている。

CSはサンプリング手法を大きく変える。

例えば1ピクセルしかなくても画像を取得できる。

http://dsp.rice.edu/cscamera

この技術のいいところは、検出器が高価でピクセル配列を準備できないものなんかでも、空間情報を取得できる点だ。(X線とか?)

大学で習うサンプリングといえばNyquist-Shannon Sampling Theoremで、これは取得したい周波数の2倍の周波数でサンプリングすれば信号を完全に復元できるというもの。

CSはもっと少ないサンプルでも信号を復元できるよ、というものらしい。

例えばカメラはRAWからJPEGに変換する際、画像は多くの情報を捨てている。DCTによって重要な係数のみを保存する。あまり重要でないものは捨ても問題ない。むしろ、保存や転送には容量が小さいほうがよい。

なんかものすごく将来性のありそうな技術で勉強を始めたけど、数学が難しいorz

自分のためにCSに学んだことをアップしていこうと思う。

以下目次

• Restricted Isometry Property
• Matching Pursuit
• Basis Pursuit
• Soft Thresholding
• など

RISEでインタラクティブなIPython notebookのスライド作成!

前にも書いたけど、ちょっとしたプレゼン資料を作るときに以下のようなことが多々ある。

手軽に書きたい
Texで数式
図を入れたい。
プログラムを入れたい。
かっこ良くしたい。
PDFにしたい。
個人で使う分はフリーで済ませたい。

しかもreveal.jsを使っていかにもギークっぽくしたい。
そんな要求に答えてくれるライブラリがある!
https://github.com/damianavila/RISE

まずはIPython notebookをインストール。

IPython notebookはAnacondaを入れると一発で入る。Pythonを使う人は入れておいて損はないけど、他の言語でIPython notebook使う人にはちょっと大きすぎるかも。

RISEをインスール

https://github.com/damianavila/RISE
をどこかに保存して,そのディレクトリに移動した後に

python setup.py install

するとIPython notebook上に新しいボタンが現れる。

ボタンを押すとスライドになって☓を押すと元に戻る。

そのままではすべて1枚のスライドになってしまうのでCell ToolbarでSub slideとかに設定すればよい。

これから会社ではこれ使って発表しよう。

​

ipython notebookでVim(バインディング)

ipython notebook(Jupyter)でもvim bindingが使えるらしい。
https://github.com/lambdalisue/jupyter-vim-binding

ipython notebookで使いにくかった点が1つ改善された。

環境はMacにAnacondaを導入している。
基本的な導入は上のGithubのリンクに書いてあるけど…一応メモしておこう。

まずはipythonをアップデートして、さらにjupyterを入れる。
(Jupyterを入れないとnotebookモジュールがないと怒られる)

$ conda update ipython
$ conda install jupyter

次に以下を実行

from notebook.nbextensions import install_nbextension
install_nbextension('https://goo.gl/5TK96v', user=True, destination="vim_binding.js")

一時的に使うときは

%%javascript
Jupyter.utils.load_extensions('vim_binding')

ずっと有効にしたいときは

%%javascript
Jupyter.notebook.config.update({
  'load_extensions': { 'vim_binding': true },
});

これでipython notebook上でVimバインディングが使える。いやー、便利になった。

​

Iteratively Reweighted Least Squares についてサクッと。

Iteratively Reweighted Least Squares(IRLS)ってたまに出てくるけど、何か分かってなかった.

Least Squares

LSでは線形モデル、

というモデルを考えたときにデータとモデルの2乗和誤差を最小にするようにを決定する。このとき誤差を1つの定数で書ける(データのばらつきが次元で同じ)ときは擬似逆行列を使って係数を決定することができる。導出はないけど、微分が0になるようにを決定すれば

Weighted Least Squares

先程はデータの次元でばらつきが同じとしたけど、次元間の分散が違うことだってある。そんな場合にWLSをつかう。調べた感じWikipediaが一番分かりやすい。
https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Weighted_linear_least_squares
この場合、

となる。ただしは対角成分に分散の逆数を並べた行列。例えば、あるデータ次元についてはほぼ確定的なら始めから分散は小さくできるので、それを反映させることができる。式を見ると、データを標準偏差で割っているのと同じになっているので、データが事前にある場合は標準化(standardization)しとけば同じことになりそうだ。

Iteratively Reweighted Least Squares(IRLS)

上の2つは一発で最適な係数を求めることができた。これは2乗和誤差がの2次で書けるからだ。IRLSでは、2次でない場合でも2次で近似しながら重みWを更新して、最適な係数を求めてくれる。これまたWikipediaがわかりやすかった。
https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares

Lpノルムを最小化したいときは、

は対角行列。と初期化しておき、対角成分の各要素は以下のように更新していく。

IRLSはスパースな係数を選択するためにL1最適化を使う場合に用いられたりするようだ。
コンピュータビジョン最先端ガイド6巻の3章に少しだけIRLSの記述がある。
Scikit-learnなどではOMPが実装されているみたいだが。

​

Python でPRML7.2章 RVMを実装。

いつ使うかは分からんが、RVMを実装しておこう。
7.2章の数式をそのまま書いた。動作確認は1次元のガウス基底で行った。

def rvm(t, PHI, alphas0, beta0, maxloop=10):
    """
    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.
    """
    beta = beta0
    alphas = alphas0
    
    for i in range(maxloop):
        A = np.diag(alphas)
        V = np.linalg.inv(A+beta*np.dot(PHI.T, PHI))
        m = beta*np.dot(V, np.dot(PHI.T, t))
        gammas = 1 - alphas * np.diag(V)
        
        #new alphas and beta
        alphas = gammas / m**2 
        beta =  (N - np.sum(gammas)) / np.linalg.norm(t - np.dot(PHI, m))
        
    return alphas, beta, m, V

10個のガウス基底でフィッティングをした。
今回書いたコードではループの回数が多すぎるとzero divisionでエラーとなった。でもループは10回くらいで十分フィットしているように見えた。

Python3 で OpenCV3を使う(Windows)。

最近やりたいことはほとんどPython3で出来ているが、OpenCVを使う時だけpython2に切り替えている。OpenCVもPython3に対応したようだが、配布バイナリはPython2のものしか提供されていない。オフィシャルでないけど、サクッとPython3にOpenCV3をインストールする方法を見つけたので紹介する。


まず、ここから.whlファイルをダウンロードする。
(このサイト経由だと余計なものをダウンロードしなくて済むのでよい)
http://www.lfd.uci.edu/~gohlke/pythonlibs/#opencv

whlはpipでインストールできるので、

pip install opencv_python-3.0.0-cp34-none-win_amd64.whl

とする。

前はソースからコンパイルとかしたけど、こんなにラクになるとは。ありがとうございます。

matplotlibのcolormap

scipy2015に興味深い講演があった。
A Better Default Colormap for Matplotlib
https://www.youtube.com/watch?v=xAoljeRJ3lU


Matplotlibの現在のデフォルトの色はjet。
Matplotlib2.0からはデフォルトの色が'viridis'に変わる。
http://betterfigures.org/2015/07/10/a-welcome-development-for-matplotlib/
'viridis'は4つの候補のうちの1つで、投票で選ばれたらしい(緑が入ってるから)。4つともMatplotlib1.5から導入はされるらしい。

なにが良くなったかというと、
色空間を連続的に変化させるだけでなく、人間の視覚を通した時にスムーズに変化するようになっているらしい。こうすることによって、より細かい変化に気付くことができる。たとえば、医師の誤診断とかも減るとか。更に、白黒にしたとしても分かるし、色盲の人でもわかりやすいように設計されているらしい。

ここには動画もある。
https://bids.github.io/colormap/

ちなみに今使えるcolromapサンプルはここにある。
http://matplotlib.org/examples/color/colormaps_reference.html

Pythonで多項式回帰

pythonで最小二乗法で係数を求めるのには関数を定義しなければいけないが、高次の多項式だとクロスタームやらいっぱい出てきて書くのが大変だ。

scikit-learnのPolynomialFeaturesはそこを自動で計算してくれる。
例えば2次の場合degree=2とすれば、

[x1, y1] -> [x1**2, x1*y1, y1**2, x1, y1]

を計算してくれる。

それを使って多項式近似をする。
以下メモ。
Pipelineについては
http://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html

from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline

#original function
def f(x):
    return 2*x[0]**2 + 3*x[0]*x[1] + 4*x[1]**2 + 2*x[0] + 3*x[1] + 4

#observed data
xx, yy = np.meshgrid(np.linspace(-2, 2, 10), np.linspace(-2, 2, 10))
X = np.vstack((xx.reshape(-1), yy.reshape(-1))).T
y = list(map(f, X))

#fit the data
degree = 2
model = make_pipeline(PolynomialFeatures(degree), Ridge())
model.fit(X, y)

#accessing to the coeffs
#(model is a Pipeline.)
print(model.steps[1][1].coef_)

Python でサクッとLeast Square fit

Pythonで最小二乗を使うときはサクッと以下のように行おう。

import numpy as np
from scipy.optimize import leastsq

def func(c, x):
    return c[0]*x**2 + c[1]*x + c[2] 

def errfunc(c, x, t):
    return t - func(c, x)

#observed data and initial guess
x = np.linspace(-5, 5, 100)
t = 4*x**2 + 1 
c0 = [0, 0, 0]

leastsq(errfunc, c0, args=(x, t))

AnacondaでPython2と3のIPython notebookを起動する。

Pythonで計算をするときはAnacondaを入れとけばとりあえず色々できる。

Python2と3の両方を使うためにAnacondaとAnaconda3を2つインストールしていたが、そんなことしなくてもよいことがわかった。Anaconda3でPython2を扱う方法をメモしておく。


まず、AnacondaでPython2の環境を作る。これはいろいろなサイトに載っている通りに行う。

>conda create -n python2 python=2.7 anaconda

IPython notebookをPython2で使うときは

>activate python2
>ipython notebook

SublimetextでPython2を使いたいときは

{
    "shell": true,
    "cmd": ["C:\\Users\\myname\\Anaconda3\\envs\\python2\\python.exe", "-u", "$file"],
    "file_regex": "^[ ]*File \"(...*?)\", line ([0-9]*)"
}

Lispを保存

Land of Lispを読み始めた。 といっても実は2回目の挑戦である。一回目に読んだときはパソコンの電源を落としたら、インタプリタに打ち込んだ変数やら関数がなくなってしまって挫折した。

今回はそうならないようにLispの環境を保存できるようにする。
インストールしたのはGNU CLisp.

保存する場合は

(ext:saveinitmem "lispinit.mem")

読み込んで始める場合は、

clisp -norc -M lispinit.mem 

インタプリタで開発は書きづらいのでエディタで書いてそれを上のコマンド保存、読み出しをするのがやりやすそうだ（個人的には）

Xperia Z SO-02E MVNOのSIMでもテザリングする。

ドコモのXperia Zを使って2年がたった。
乗り換えばかりを優遇する3大キャリアに嫌気がさしたのと、普段月2Gも使って居ないのでMVNOにMNPした。

たまに電車の中とかでテザリングしたいときがあったけど、そのままのXperiaではテザリングできないらしい。いやーひどい。

調べてみると、WifiテザリングはRoot化しないと無理とかいう情報ばかりだったけど、したのリンクの記事にRoot化しなくても行けると書いてあった。
http://androplus.org/Entry/282

ここに書いてある通りにやればadb.exeもダウンロードできてうまくいった。
http://xperia-freaks.org/2015/03/21/docomo-tethering/


テザリングはうまくいったけど、たまにspモードがどうとかいう警告がでるのがうざい。
これはRoot化しないときえないのか…

Juliaで外部のプログラムを動かすときのメモ

Juliaで外部のプログラムを動かすときのメモ。
PythonではPopenとかを使ってたので、Juliaでもそれと同じようなものがあるかなーと思ったけど、なかなか見つからなかった。

公式マニュアルはこれ。
http://julia.readthedocs.org/en/latest/manual/running-external-programs/

外部プログラムを起動して、そこから標準入力を扱いたいときは、こっちが参考になった。http://blog.leahhanson.us/running-shell-commands-from-julia.html

概要は、
プログラムが起動して、
write()で書き込み。書き込んだらその出力がsoにたまっていくので、
readall()でその出力を得る。

(so,si,pr) = readandwrite(`yourprogram`)
write(si, "command1")
write(si, "command2")
close(si)
output = readall(so)

これも一応メモ。
http://julialang.org/blog/2013/04/put-this-in-your-pipe/

Novel Class

アメリカ滞在の後半はNovel Classを取った。文字通り小説の授業だ。

半期に3つの小説を読んだ。1つはベトナム戦争についての話。もうひとつはThe Catcher In The Rye。最後にSidartha。Sidarthaはヘルマン・ヘッセの小説。どれも一字一句辞書で調べながらじっくり読んだんで、ものすごく思い出に残っている。

ベトナム戦争についての話はあんまり覚えてない。ただ、そのなかの参考資料にBreaking Awayという映画があった気がする。本当にこの小説の授業で出会ったものかは定かではないえけど、ものすごく感傷にひたったのを覚えている。Breaking Awayは青春を舞台にしたもので4人の登場人物がそれぞれ精一杯行きている。あーおれはこんなところで独りでなにをやっているんだろう。と心の通った友達とのチームプレイがいかに大事かが身に染みたのを覚えている。

2つ目のThe Catcher In the Ryeは有名な小説だが、村上春樹の日本語訳はまったくおすすめしない。帰国後日本語訳を手にとって読んでみたけど、まったく感情移入することが出来んかった。こっちはあまり前向きな話ではない。最後にDon't tell anybody anything. You will miss somebody.というのがすごく思い出に残っている。他にもI have to admit itとかThat kills me.とかものすごく当時のこころに響くセリフばかりだった。

Sidarthaは難しい本が読めなかったので、優しい英訳の本を読んだ。Sidarthaが悟りをひらくまでの苦悩が思い出に残っているが、詳しい内容は忘れたw


それぞれの本を読んだ後の課題は自由。どんな方法でも読んだ感想を表現すればよい。模型をつくった人もいれば、冊子を作った人もいた。日本のいわゆる読書感想文みたいに決まった書式はない。自分の得意な方法で表現すればよかった。いま日本にこんな教育をするところはあるんかな?

Julia Array操作の基本をメモ。

Pythonでやっていたlistの操作にすっかり慣れてしまったけど、JuliaのArrayではどうやるか確かめる。
Julia0.3で試して動いたものを掲載。

配列の生成

a = [i for i in 1:10]     #1から10を要素に持つ配列生成
a = [1:10]                #同上
b = reshape([1:9], (3,3)) #3x3の2次元配列。
b = [1 2 3; 4 5 6; 7 8 9] #3x3の2次元配列。上とは行、列が逆(転置の関係)

要素数

length(a)

要素アクセス

a[1]
a[1:end]   #要素1から最後まで
a[2:3:8] #要素2から8を3個おきに

追加,挿入

push!(a, 13)         #要素と同じ型追加
append!(a, [1,2,3])  #末尾に[1,2,3]を追加
prepend!(a, [1,2,3]) #先頭に[1,2,3]を追加
insert!(a, 3, 15)    #3番目に15を挿入する。(要素数は増える)
vcat(a, [11,12,13])  #aに末尾に[11,12,13]を追加(aは非破壊になる)
hcat(b, [10,11,12])  #bに列を追加
vcat(b, [10,11,12]') #bに行を追加('は転置)

要素数があらかじめわかっているなら、始めに配列を確保して代入したほうが速いっぽい。

置換

a[1] = 15
a[1:4] = [4,4,4,4]
#a[2:5]を削除して[0,0,0]を挿入。戻り値は削除した元のaの[2:5]の値。
splice!(a, 2:5, [0,0,0]) 

削除

shift!(a)  #最初の要素を削除
pop!(a)    #最後の要素を削除

検索

3 in a
findfirst(a, 4) #はじめに4が出る要素番号、なければ0を返す
findin(a, 4)    #4の存在する要素数をArrayに入れて返す)

Lytroが日本に上陸してアクセスが増えたのでLytroのすごいと思うところをのべる。

Lytroが日本にやってくるらしい。
http://dc.watch.impress.co.jp/docs/news/20141105_674603.html

なんか記事のアクセスがすごい。
これは良い機会とばかりにLytroがすごいと思うところを挙げる。1つは文化的内容。もう一つは技術的内容。

1.　アナログとデジタルの融合が進んでいる。

上の記事にレンズの断面図がのってある。35mm判換算で30-250でFナンバーは2.0らしい。最近のカメラと比べると広角側が少し足りないと思われるかもしれない。注目すべきはレンズのスペックが高いにも関わらず、枚数12枚で非球面レンズを使っていないことだ。
http://dc.watch.impress.co.jp/img/dcw/docs/674/603/html/17.jpg.html

　確かIllumの発表の時に言及しとったと思うけど、レンズの収差もデジタルで補正できるような最適設定にしているらしい。収差のデジタル補正と言えば富士フイルムとかキヤノンにそういう機能(デコンボリューション)がついていたと思うけど、こういう技術はどちらかというと、まだレンズの最適化で取れない部分をデジタルで補正しているイメージ。しかし、Lytroの場合発想が逆で、デジタルで補正できるところまでレンズの収差を出しているようだ。Ngさんの博士論文にもそういう研究はあったし。もちろん日本の企業でも「倍率色収差」とか「歪曲収差」とかは、デジタルで補正できるようにセンサーとレンズの最適化を行っとるやろうけど、Lytroの最適化は1つ先の次元に行っとる気がする。

2.　技術をもって本気で新しい表現を模索している

Lytroは新しい表現を探している。ホームページのムービーでも言っとるけど、Light Fieldカメラで一瞬を切り取るとどういう表現ができるか模索している。最初に出した製品はセンサーも小さく使い勝手も悪かったが、今回のIllumはクリエイターをターゲットにしている。Light Fieldカメラが生き残るかどうかは、新しい表現方法が見つかるかどうかにかかっているといってもいい気がする。ボケ味を追求する日本のメーカーとは次元が違う。(もちろんどっちがいいとかいう話は抜きにして)
残念ながらぼくにはどういう表現がよいのか全くわかりません。

やっぱりSilicon Valleyの企業は日本とは考え方が違う。
ただLytroのIllumは使い勝手が悪いらしい。やっぱり手に馴染む機器の作り込みには時間がかかるようだ。

Juliaで固定配列

JuliaでArrayを使えばベクトルや行列計算できるけど、2とか3次元の点を表すのに可変長の配列を使うのはどうも無駄が多い気がする。自前で

immutable Vec
    x::Float64
    y::Float64
    z::Float64
end

とかやったけど、演算子とかを定義するのが面倒。

そんな問題を解決してくれるのがImuutableArrays。
Julia0.3にはすでに採用されているので、

Using ImuutableArrays
v3 = Vector3{Float64}(1, 2, 3)
v2 = Vector2{Float64}(4, 5)

と書ける。デフォルトでは読み込んだ時に4次元まで作ってくれて、演算子の定義とかもやってくれている。(これがメタプログラミングの力か…)
しかも上で自分で定義したのより速い…

行列は

typealias Mat3d Matrix3x3{Float64}

とすれば使える。

IPython notebook でプレゼンテーション用スライドを作る

プレゼンテーションを作るときに困ることがある。

1. 手軽に書きたい
2. Texで数式
3. 図を入れたい。
4. プログラムを入れたい。
5. かっこ良くしたい。
6. PDFにしたい。
7. 個人で使う分はフリーで済ませたい。

ざっと思いついたままに書くとこうなる。

いろいろ探していたら、最近IPython notebookからスライドが作成できることが判明した!

しかもreveal.jsを使っているのでいかにもギークっぽいかっこいいスライドが作れる。

こんな感じのやつ→http://lab.hakim.se/reveal-js/#/

まえからreveal.jsは使ってみたかったけど、jsとかHTMLとかよく分からんし、Web上で使えるやつもあったけど、PDFにするのに有料会員にならんといけんとかで挫折した。

スライドの作り方

1.　IPython Notebookをインストールする。

IPython NotebookはAnacondaをつかってインストールするのがおすすめ。一式全部はいるし、アップデートも

conda update anaconda

で完了する。

2.　新しいNotebookを開いて設定する。

New Notebookを作成したら以下の図のところをSlideに設定する。

右のタブはSlide、Sub-Slideとか選べる。これはreveal.jsでいう横にスライドか、縦にスライドかを指定できる。いろいろ試してみて。

3.　スライドへ変換する

IPython Notebookの編集が終わったら、NotebookからSlideへの変換をする。スライドの作成はここにあるとおりのコマンドを実行するとできる。--postでHTMLサーバーを指定する必要があるらしい。

ipython nbconvert notebook.ipynb --to slides --post serve

IPythonのバージョンが古いと文字がものすごく小さいので最新版にアップデートするのをおすすめ。
※Chromeで開くと数式がレンダリングされなかったのでmathjaxのプラグインをいれた。
ちなみにPDFにするときは

ipython nbconvert notebook.ipynb --to latex --post pdf

問題は記録用にPDFとかHTMLにしたらちょっとダサくなってしまうことか。
まぁよい。

Opencv3.0 ALPHA Python版をMacにソースからインストールしたときのメモ

Mac OSXでPythonはAnacondaでインストールしたPython3.4を使用

1. OpenCV for Linux/MacをダウンロードしてZipを解凍する
2. ダウンロードフォルダに移動して、作業用フォルダを作っておく。ここではbuildとした。
3. CMakeをダウンロードして実行。
4. ダウンロードしたOpenCVのフォルダとbuildフォルダを指定する。
5. Searchにpythonと入れてConfigを押す
6. Configを押すとPython3の設定がいろいろあるので上の図や参考サイトを参考にして埋める。赤い表示がなければとりあえずOK。
7. Generateを押す。
8. buildフォルダで
   make -j8
   を実行。失敗するとエラーが出る。Python.hがないとか言われたけど、6の設定したPATHを見直すとうまくいった。
9. sudo make install
10. pythonを立ち上げてimport cv2
11. 今回は「 Library not loaded: libpython3.4m.dylib」とか言われたのでシンボリックリンクを貼った。
sudo ln -/Users/name/anaconda/lib/libpython3.4m.dylib /usr/local/lib

これでとりあえずopencvが動いた。

参考サイト

1. http://luigolas.com/blog/2014/09/15/install-opencv3-with-python-3-mac-osx/
2. https://stackoverflow.com/questions/20953273/install-opencv-for-python-3-3

Light fieldの原理をLytroがじきじきに解説

久しぶりにLytroのblogをのぞいたら、Lytroが論文を公表していた。
まだ目を全部に目を通してないけど、見たことない頭があったりした。
http://blog.lytro.com/post/97168596625/modeling-the-light-field-camera

最近はGoogle Cameraでもボケを作れるので静止画でのLytroの優位性はなくなってきている気が。動画だと処理が大変そうだしなー。