Julia - Healpixを使いたい
JuliaでHealpixを使いたい！！
一つの解はHealpyをbandleする。juliaのみでやりたいなら一つは開発中らしいHealpix.jlを使うこと。ただしたとえばpix2angとか重要な機能がまだ実装されていない。 これは簡単。両方のやり方のメモ
Healpyを使う
Healpyがインストル済みなら
using PyCall @pyimport healpy as hp
これで使えます。まあなんかこれで十分かも、、、
ただしhealpy使う注意点としては、pythonはindexが0から始まるを忘れないことだ！！
Healpix.jl
Healpix.jlをbuildする
https://github.com/mweastwood/LibHealpix.jl
をinstallするのにつまずいた。これはPkg.buildとはbuild.jlというのを実行しているだけなので 上手くいかなかったらそこまでもどろう。
Pkg.add("LibHealpix")
Pkg.build("LibHealpix")
後者Healpix 3.20が
~/.julia/v0.4/LibHealpix/deps/downloads/Healpix_3.20
にダウンロードされるが、私の環境だとOS-10.9.5/gcc-4.2だと g++で
-fno-tree-fre
が引っかかってとまるのだ。 そこで、
~/.julia/v0.4/LibHealpix/deps/downloads/Healpix_3.20> cd src/cxx/config ~/.julia/v0.4/LibHealpix/deps/downloads/Healpix_3.20/src/cxx/config> grep -fno-tree-fre *
として-fno-tree-freのと書いてあるファイルを探し、-fno-tree-fre部分のみ全部コメントアウトして
~/.julia/v0.4/LibHealpix/deps/downloads/Healpix_3.20> make
するとうまくいった。Pkg.build()はbuild.jlなるファイルを実行するだけなので
~/.julia/v0.4/LibHealpix/deps> emacs -nw build.jl
をみて今やった部分をコメントアウト。最後の
# Build the HEALPix wrapper
println("Building the HEALPix wrapper...")
dir = joinpath(depsdir,"src")
run(`make -C $dir`)
run(`make -C $dir install`)
ここだけ残して、
~/.julia/v0.4/LibHealpix/deps> julia build.jl
これで、
julia> Pkg.test("LibHealpix")
INFO: Testing LibHealpix
INFO: LibHealpix tests passed
のようにテストが上手くいくことを確認した。以上。
Julia - fortran 90をwrap (bind) する
FortranのコードをJuliaに変換するのは、ものによっては結構しんどいので、juliaからfortran (90)を呼び出そうという例です。
最も簡単な例
まず、ただ入力値をprintするだけのsubroutine printintを格納したfortran90 module testfを用意します。
module for a simple julia wrapper of fortran90
次にこれをshared libraryにします。
~> gfortran testf.f90 -o testf.so -shared -fPIC
Juliaのwrapperはこんな感じです
simple julia wrapper of fortran 90
ここでポイントはmodule testfのsubroutine printintを呼び出すには
:__testf_MOD_printint
と記述する所です。
~> julia wrap_ftest.jl
7
こんな感じで呼び出せます。
出力も得る
上の例ではinputだけでしたが、fortran90側が
のように、integerを入力して、+3したintegerを返す場合はどうでしょうか？
この場合、julia側では
こんな感じで呼び出せます。
~> julia wrap_ftest.jl Int32[7] Int32[10]
のように返ってきます。
他の型
他の型は、 http://docs.julialang.org/en/release-0.4/manual/calling-c-and-fortran-code/ を参考に対応関係を知るとできます。
- Float32 <-> Real*4
- Float64 <-> Real*8
- Complex64 <-> Complex*8
- Complex128 <-> Complex*16
などなど。Complex*16の場合の、複素共役を返す例を以下に与えた。
fortran 側
julia 側
配列を入出力する
さて次は配列の入出力である。複素数配列の複素共役を返す例を以下に与えた。
fortran 側
julia側
ここでポイントとなるのは配列の要素数を与えるnをccallの中で、&nであたえていることだ。
Generalized Coltrane Changes
- 本記事は、 divrot-hajime.tumblr.com の http://divrot-hajime.tumblr.com/post/125174466518/generalized-coltrane-changes http://divrot-hajime.tumblr.com/post/125341867738/generalized-coltrane-changes-%E3%81%AE%E5%AE%9F%E8%A3%85 初出: Jul 28th, 2015 & Jul 29th, 2015を改訂して再掲したものです。 (gistとの連携をやってみたかったんです)
Coltrane changesは、IIm-V7-I もしくはV7-Iなどが+5度上に転調していく進行のことである。これを数学的に表記してみよう。Giant Stepsを例にとると
B D7 G Bb7 Eb ~
という感じではじまるが、これはトニックだけをとりだすと
B G Eb
なので、+5度上に+5度上に動いているだけである。もう一度+5度上にいくと戻って
B G Eb B
になるわけだから三回の操作で元に戻る循環である。これは+5度とかいうと分かりにくいのだけれど、Cから順に、
C=0, Db=1, D=2, Eb=3, E=4, F=5, Gb=6, G=7, Ab=8, A=9, Bb=10, B=11
という風に数字(index)を割り付けると考えやすい。オクターブを無視すると音階は12回でもとにもどるから、１２以上のindexや負のindexやindexは12で割った余りと等しいとする（つまり12を法とした整数を考える）。
そうすると+5度上にいくという変換F8は、単にindex iに8を足すだけである
F8(i) = i + 8 (mod 12)
具体的に見てみよう。Bはi=11だから一回目の変換は
F8(11) = 11 + 8 (mod 12) = 7 つまりGになる。もう一度変換すると
F8(7) = 7 + 8 (mod 12) = 3 つまりEbとなっている。さらに変換すると
F8(3) = 3 + 8 (mod 12) = 11 でBに戻っている
上では8を足す変換を考えたが、いま一般にnを足す変換
Fn(i) = i + n (mod 12)
は循環をつくる。例えばF5(i)は普通のいわゆるDominant motion (C-F-Bb-Eb-Ab….)を表現できる。
では、もともとのGiant Steps
B D7 G Bb7 Eb ~
はどのような変換で表記できるだろうか？これにはnを偶数回目と奇数回目の変換で異なる物をとればよい。kを変換の回数として
F(i) = i + n_k (mod 12)
n_k = 5 (kが偶数) 、n_k = 3 (k が奇数)
さらにkが偶数のときはt_1=Major、奇数のときは t_2=7thをとるとする。
具体的に見てみよう。Bはi=11だからk=1回目の変換は
F(11) = 11 + 3 (mod 12) = 2 つまりD7になる。k=2回目の変換は F(2) = 2 + 5 (mod 12) = 7 つまりGになる。k=3回目の変換は F(7) = 7 + 3 (mod 12) = 10 つまりBb7になる。k=3回目の変換は F(10) = 10 + 5 (mod 12) = 3 つまりEbになり、Giant Stepsの最初の部分が再現される。ここで3+5=8となっていることが、先ほどの+5度の進行に対応している。
というわけで、このColtrane Changeの法則を一般化してみよう。いまn_kをM回で元に戻る数列{n_1, n_2, … , n_M}
n_k = n_(k mod M)
だとして、コードの部分をL回で元に戻るハーモニー部分{t_1, t_2,…,t_L}とする(ここにはMajorだとか7thだとかdiminishだとかがはいる)
t_k = t_(k mod L)
一般化コルトレーンチェンジ列はこれら{n_1, n_2, … , n_M}と{t_1, t_2,…,t_L}
を用い ルート音: F(i) = i + n_k (mod 12) とハーモニー t_kで構成される。
一番下に簡単なPythonコードを添付しておく。Giant stepsの場合はn_k={5,3}, t_k={Major, 7th}である。計算結果
GCC: for n_k= [5, 3] t_k= [“, ‘7’] B D7 G Bb7 Eb Gb7 B D7 G Bb7 Eb Gb7 B D7 G Bb7 Eb Gb7 B D7 G Bb7 Eb Gb7 B
ちゃんと再現されてますね。
LはMと等しくなくても良い。一般化コルトレーンチェンジで様々なコード進行が作成できるので試されたし、
GCC: for n_k= [5, 3, 1] and t_k= ['mb5’, '7’, ’m’, '7’, ”] Bmb5 D7 Ebm Ab7 B Cmb5 F7 Abm A7 D Fmb5 Gb7 Bm D7 Eb Abmb5 B7 Cm F7 Ab Amb5 D7 Fm Gb7 B Dmb5 Eb7 Abm B7 C Fmb5 Ab7 Am D7 F Gbmb5 B7 Dm Eb7 Ab Bmb5 C7 Fm Ab7 A Dmb5 F7 Gbm B7 D Ebmb5 Ab7 Bm C7 F Abmb5 A7 Dm F7 Gb Bmb5
※ レヴィ・ストロース「親族の基本構造」内にあるアンドレ・ヴェイユの補遺を参考にしました。
追加
pyaudioを用いて、一般化コルトレーンチェンジの音が出るようにしました。下にコードあり。Thanks, matu_mitani !
GCC音声ファイル: http://secondearths.sakura.ne.jp/music/gcc/
ただ、正弦波はイマイチなので波形はそのうち変えたい。
matu_mitaniのご助言に従って倍音成分入れました。コードが聞きやすくなった(30 July 2015)。それからedgeをexpでdampしてプチプチいうのを防ぎました。(31 July 2015)
Hajime_dxdydz (July 29, 2015)
環境
使用法
./gcc.py -key 11 -nk 5 3 -tk Maj 7 -o coltrane.wav
Julia - install, package, and emacs
Juliaを始めることにした(=つまりJuliaが使われることに賭けることにした) -
Macにインストール
Mac版の単体Julia 0.4.2をインストール。
これまでのpythonの仕事環境としては、ターミナル上でemacs開いて書いて、実行するという古いスタイルなので、慣れ親しんだこのスタイル＋Jupyterで仕様する。
bin本体は/Applications/Julia-0.4.2.app/Contents/Resources/julia/bin/julia
PATHに追加する。
~/juliatest> julia test.jl
これで動くようになった。
ソース頭に
#!/Applications/Julia-0.4.2.app/Contents/Resources/julia/bin/julia
とかけば、そのまま実行可能になる。
emacs環境をいれる
Juliaのemacs環境はjulia-mode.elとESSがあるみたいだ。ESSはうまくいかなかったので、前者をいれる。julia-mode.elをgithubからgetして、.emacs.d/site-lisp/にいれて、.emacsに以下を書く
(load “/Users/name/.emacs.d/site-lisp/julia-mode.el” nil t)
(autoload ‘julia-mode “julia-mode” “Emacs mode for Julia” t)
(add-to-list 'auto-mode-alist ’(“\.jl\’” . julia-mode))
これでつかえるようになった。
Package管理
Packageは専用マネージャーみたいのがある
julia> Pkg.init()
こんな感じでaddするようだ
julia> Pkg.add(“Gadfly”)
消す時は
julia> Pkg.rm(“Gadfly”)
何がinstallされているか確認する時は
julia> Pkg.status()
これらは.julia内にどんどんinstallされていく
はまった幾つか
複素行列の転置(transpose)。Juliaでは行列の転置を' (プライム)で省略してかけるが、複素行列の場合、単なる転置ではなく、複素共役転置になる。
Juliaではindexは1から始まる。pythonとか多くの言語では0から始まるので、importするとき面倒くさい。離散フーリエ変換関係でよくでる配列の変換はたとえば以下。
