渋谷で働くSEのブログ

渋谷で働くSEのお粗末なブログです

中南米音楽の形を見てみたい(3)まずは論文の再現実験から

f:id:kyosuu1:20161202125428p:plainえーと、えーと…
論文の読み間違いで、再現実験に苦戦していました。

実はですね…
このシリーズの(2)で紹介した距離なんですが、
私が紹介した距離は「chord class distance」というもので、
コード間の距離を測るものだそうです。

メロディの距離を測るなら、
足し合わせでいいとのこと。

つまり、f = (f_1,\,f_2,\,\cdots,\,f_n),\,g = (g_1,\,g_2,\,\cdots,\,g_n)について、
 d_{PC}(f,\,g) = \displaystyle\sum _{i=1} ^{n} d_{PC} (f_i,\,\sigma(g_i)),
ということですね。大変失礼いたしました。

それを踏まえて…
今回はいよいよ再現実験です!

まずは「Abbott's Bromley Horn Dance」からです。
こんな曲なんですって。


f:id:kyosuu1:20161201001146p:plain

なんか無印良品っぽいな…と思いながら楽譜を見つめました。

それでですね…
再現実験はRでやろうと思います!
Rで、TDA計算ライブラリのDionysusを用いて計算をします。

Rでは、ユークリッド距離はもちろん、任意の距離行列からもpersistentの計算ができます。
こうやります。

>> library(TDA)
>> rips<-ripsDiag(X, maxdimension = 2, maxscale = 1, dist = "euclidean", library = "GUDHI", location = FALSE, printProgress = FALSE)
>> plot(rips$diagram, barcode = TRUE)
説明です
  • X: 行列です。point cloud もしくは距離行列です。
    • Xがpoint cloudの場合、d次元の点がn個あるなら、n \times d行列を入れます。
    • Xが距離行列の場合は、行列(d(x_i, x_j))_{ij}(つまり、点x_iと点x_jとの距離d(x_i, x_j)を成分として持つ行列。距離の定義より、これは対称行列となる)
  • maxdimension: 文字通り最大次元。persistentを求めたい最大の次元です。
  • maxscale: Rips diagram作成時の、円の最大半径のことです。
  • dist: 距離を指定します。defaultは"euclidean"です。任意の距離空間を使う場合は"arbitrary"と指定します。
  • library: 使うライブラリです。
    • GUDHI, PHAT, Dionysusの3種類があります。
    • 距離行列を使う場合はDionysusのみ使えます。
    • 任意の距離を使う場合、ライブラリが指定しないか、指定が間違っていれば自動的にDionysusを呼び出して計算してくれます。
  • location: TRUEになっていると、birth pointとdeath pointを返してくれます。
  • progress: TRUEにすると、計算がどのぐらい進んでいるかを示してくれます。これをつけておくと、待ち時間もちょい楽しくなります。


我々が使う距離はユークリッド距離ではありません。
メロディの場合、各音の周波数にlog_2をとってその差を求めていました。
それは、ちょうど1オクターブが2倍となっているからだと言っていましたね。

ということで、Rでpersistent homologyを計算するには、距離空間が必要となります。

それで、「楽譜」から距離行列が出力できるソースコード
github.com
を書いておきました。

まだgithub側の説明はチョーテキトーです;;すみません。

使い方は以下のようになります。

python distance_matrix.py --dim 2 --cc 0 input.csv output.csv
コードの説明
  • python: このコードはpython3でのみ動きます。python2.xがdefaultの場合は、必ずpython3にしてください。
  • --dim: メロディのtupleの次元です。連続したメロディを何個まで見たいかという指定です。指定できる数字は1から3です。
  • --cc: chord class distanceを求めたいとき指定します。
    • 0だとchord class distanceを求めてくれます。
    • 次元が2か3の場合のみ有効です。
    • 何も入力しなければ、pitch class distanceで計算してくれます。
  • input.csvの形式: 音階を小文字で書いた、n行1列の文書です。

各音階の名前は、次の楽譜の青文字ように入れます。
f:id:kyosuu1:20161202125428p:plain
githubにサンプルを上げておきましたので、参考にしてください。

それでですね…
いよいよ距離行列も作れたわけです。

その距離行列をRに入れて、persistent図を描くなら次のようにします。

>> library(TDA)
>> dist_mat<-read.csv("path/to/your/input.csv", header=F)
>> rips<-ripsDiag(X=dist_mat, maxdimension = 2, maxscale = 1, dist = "arbitrary", library = "Dionysus", printProgress = TRUE)
>> plot(rips$diagram, barcode = TRUE)

私が読んだ論文では3次元まで計算していますが、Rだと3次元で計算させた途端落ちてしまいます…
というわけで次元は低めに設定していますT_T

でやってみますと…

f:id:kyosuu1:20161203003423p:plain

ほお、論文の

f:id:kyosuu1:20161203000444p:plain

これと似てますね。おおっ、これいけそうじゃん♪

そこで3次元もやってみるわけですね。

f:id:kyosuu1:20161203003459p:plain

これがRの結果で、

f:id:kyosuu1:20161203000448p:plain

こちらはJavaplexの結果です。

DionysusとJavaplexの出力方式が結構違っていて(dionysusは低い次元が下に、javaplexは上にあったり、javaplexは長さごとにちゃんと整理されていたり)違うように見えるかもしれませんが、
よくよく比べてみるとわりと似てます^^

しかし、RでDionysusを動かしていると、次元を3にした途端落ちます。
これじゃ…ちょっと使えないですね。

とりあえずは、評価版のmathematicaにjavaplexを入れて使えたりはしますが、
これも評価版なので期限があります。

やっぱjavaで書いたほうがよさそうですね。

今日は長々と失礼しました^^