2013年4月17日水曜日

Numba




Pythonでforをまわすような数値計算を行うとmatlabなどと同じで非常に遅い。
それを克服するためにpypycythonなどがある。しかし基本的にコードを書き換える必要があるので、もっと楽をしたい。

最近numbaというものを知った。ここのサンプルにある例だと、

from numba import autojit

@autojit

のおまじないだけで関数が劇的に早くなる。

少し、自分で例で実験してみた。
全球のSSTに対して、NINO3の相関をとるということをやってみる(これはnumbaの実験のために作った例であり、相関をとるだけならもっとエレガントな方法があるだろう)。

そのノートブックが
http://nbviewer.ipython.org/5396344

説明
(1) 図をプロットするためのおまじない。

(2) SSTアノマリとNINO3を読む関数を定義。
  元データは NOAA OISST v2でmonthly dataのsst.mnmean.ncを使う。
  データはPython Interface to GrADSを使って読む。
  1982から2010年のデータから気候値からのanomaly(Bin Guan's GrADS Script Libraryのdeseasonを使った。古いバージョンを使ったので、今は使い方が違うかも知れない)とNINO3 SSTアノマリを定義し、返す。
 ここらへんは以前の記事も参照。
 /usr/local/lib/python2.7 なんちゃらというメッセージは私の環境でインストール上のゴミなので気にしないで。

(3) 上で定義した関数を呼んでいる。

(4)  時系列xと二次元時系列(水平2次元と時間の3次元)から2次元の相関係数を返す関数。

from numba import autojit
@autojit
の部分がnumba.

(5) ベンチマーク
  この場合、5.7秒かかっている。
   
  @autojitの代わりに@jit('f8[:,:](f8[:],f8[:,:,:])')としても速度はかわらなかった。
  
  もし上で@autojitを付けない場合、7.34秒だった。
 
  つまりnumbaを使っていることで、確かに速くなっているが、劇的にではなかった。

  これは、forでまわしている中身がプリミティブな式ではなくて、numpyの関数corrcoefであるからだろうか。

  ちなみに、欠損値処理は何もやっていないので、warningが出ている。この場合はプロットするときに陸地のデータはマスクするので、これで構わない。

  欠損値処理もできるma.corrcoefを使うとかなり遅くなる。
  numbaを使わない場合で50.2秒、numbaを47.5秒だった。

(6) 実際に相関係数を求める。

(7)  land sea maskを読む関数

(8) land sea maskを読む

(9) プロット  

  エルニーニョに典型的な馬蹄形が見られる。インド洋ダイポールへの相関も見られる。


結論
今回試したケースでは、劇的というほどではないが、ちょっとした追加で、確かに速くはなった。場合によっては劇的に速くなる場合もあるだろう。。今後の進展に期待したい技術である。


  
  

2013年4月12日金曜日

Lanczos filter と Wakari




Lanczos filterによるlow-pass filterを理解するためにIPythonノートブックを作成した。
https://www.wakari.io/nb/tmiyama/Lanczos_test

結論は、少ない点でもadjustすれば低周波で応答関数が1に近づくようにできるが、基本は点数を大きくする必要があるということだ。

前回はIPython Viewerを使ってノートブックを公開したが、今回はWakariを使用してみた。

Wakariは端末がなくとも、クラウドでいろいろ計算できるので、なかなかおもしろい。

参考文献
Duchon, C. E., 1979: Lanczos Filtering in One and Two Dimensions. Journal of Applied Meteorology, 18, 1016–1022.

気象研究ノート 221号 気象学と海洋物理学で用いられるデータ解析法 伊藤 久徳・見延 庄士郎/著 日本気象学会/刊 A4判 253p

GFD電脳Ruby小物置き場 (Library) Lanczosフィルター
http://davis.gfd-dennou.org/rubygadgets/ja/?(Library)+Lanczos%A5%D5%A5%A3%A5%EB%A5%BF%A1%BC

2013年4月5日金曜日


Torrence and Compo (1998)に出てくる カイ2条分布を理解するためのiPython Notebookを作った。
http://nbviewer.ipython.org/5311742

(追記するかも)



参考文献
Torrence, C., and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bulletin of the American Meteorological Society, 79, 61–78. http://journals.ametsoc.org/doi/abs/10.1175/1520-0477%25281998%2529079%253C0061%253AAPGTWA%253E2.0.CO%253B2

2013年4月2日火曜日

IPython Notebook




IPython Notebookで解析して、それをパブリックに表示するデモンストレーション。

作ったnotebookはIPytnon Notebook viewerを使って、ここで見られる。
http://nbviewer.ipython.org/urls/dl.dropbox.com/u/439394/ipyhon/seaice.ipynb

  1. データとしては、北極海海氷面積データのcsvでダウンロード
  2. pandasでデータを読み、9月の平均データを作る。
  3. IPythonのRmagicを使い、線形回帰を行いトレンドラインを引く



参考文献