ナンクル力学系

学んだ事を書き連ねていこう。

Posts Tagged ‘python

RailGun: C+Pythonでお手軽に数値計算プログラミングを加速させるライブラリ

leave a comment »

シミュレーションのプログラム書く時に,毎回同じようなコード書くのが嫌なので作ってみました.とりあえずドキュメント書いてみたから使ってみてよ!

Written by tkf

October 15, 2010 at 4:08 am

Posted in 研究日誌

Tagged with , ,

Python で RNN (PyRNN) を書いたので公開します

with 2 comments

ソースはBitbucketに置いている. > tkf / PyRNN / overview — bitbucket.org

何が出来るかというと,こんなのとか(インパクトが欲しかったので,学習の様子をアニメーションにしてみた):

左上がエラーの学習曲線,右上がパラメタのRMSの学習曲線,左下が教示信号とネットワーク出力の相空間プロット,右下がコンテキストの相空間プロット.このアニメを作るソースはこれ(が吐いたpngをconvert -delay 5 *.png nn.mpgで変換).

これを作った理由は,Pythonだと簡単に式を書けるからバグ入りにくい,だからCで書いた本番用のテストに使える!と思ったから.だから,かなり計算効率は悪いけど分かりやすい書き方になっている(はず).これを使って本番用のネットの一つバグが落とせたのはかなりうれしかったけど,その本番用のはラボにいる別の人のより性能悪いっぽいので両方共にバグがあるかも(おいw

という訳で,バグレポートお待ちしてます!←ココ!

あ,簡単なネットワークの式はPyRNN v0.0 documentation で説明してます.

追記

  • README.rst にインストール方法書いてるけど,実はインストールしなくても使えます. toy/ElmanNet に色々遊べるスクリプトがあるので,それを一番根元のディレクトリ(setup.py がある場所)にもってきて実行すればおk(なはず.
  • numpy と matplotlib が必要です.
  • アニメの再生速度はlog scale(っぽく)速くしてます.実際は後半待つのが超ダルいです.
  • momentumという黒魔術項を入れて,ネットワークの学習を加速しますけどそのせいで暴れています.でも暴れているのを見るのが楽しいです.

Written by tkf

November 18, 2009 at 5:17 pm

Posted in 研究日誌, PC

Tagged with , ,

matplotlib でタイル状にグラフをしきつめるようなプロットをする

leave a comment »

こういうグラフを書きたい:

https://i0.wp.com/farm3.static.flickr.com/2568/3977125592_ae7c87d90e.jpg
最近気づいたのは pylab は遅いけど matplotlib.pyplot だと早いということ. pylab を使うと簡単にかけそうだけど,たくさん書きたいので matplotlib.pyplot で頑張ることに.

どこの処理が時間かかってるのか計るコードも入れた (tm なんちゃらってのが計測している部分). 最新の matplotlib (0.99) が必要.

コード:

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid import AxesGrid
import numpy

tm = []
def tm_append(tm):
    import time
    tm.append(time.time())
def tm_print_diff(tm):
    for (t0,t1) in zip(tm[:-1],tm[1:]):
        print t1 - t0

list_func = [ numpy.sin,
              lambda x: numpy.exp(-0.1*x),
              lambda x: numpy.sin(x)*numpy.exp(-0.1*x) ]
list_tmax = [5, 10, 15, 20]
nrows_ncols = (len(list_func), len(list_tmax))

tm_append(tm) ## time
F = plt.figure(1)
F.clf()
tm_append(tm) ## time
grid = AxesGrid(F, 111, # similar to subplot(111)
                nrows_ncols = nrows_ncols,
                axes_pad = 0.0,
                add_all=True,
                aspect=False,
                )
tm_append(tm) ## time

for ax in grid:
    ax.cla()
grid.set_label_mode('L')
tm_append(tm) ## time

i = 0
for f in list_func:
    for tmax in list_tmax:
        t = numpy.arange(0,tmax,0.1)
        grid[i].plot(t, f(t))
        i += 1
tm_append(tm) ## time

plt.draw()
plt.savefig('axesgrid.png')
#plt.show()
tm_append(tm) ## time

tm_print_diff(tm)

結果:

1.4825220108
2.09476304054
0.745758056641
0.0568108558655
1.95767402649

やっぱり, plt.figure (1行目) や AxesGrid (2行目) は時間かかるっぽいな. 二回目以降は for ax in grid 以下を実行してやれば良いので,かなりお得かな. ax.cla() した後で grid.set_label_mode(‘L’) して一番下と一番左以外の軸にラベルがつかないように しているのがミソ.

Written by tkf

October 3, 2009 at 10:50 pm

Posted in PC

Tagged with ,

PySmell で標準ライブラリ用の PYSMELLTAGS を作る

leave a comment »

pysmell /usr/lib/python2.5/ -o ~/PYSMELLTAGS とやると重すぎて使えなかったけど, pysmell /usr/lib/python2.5/ -x site-packages test -o ~/PYSMELLTAGS とやるとまともな速度で動いた.

ただ, os. とかの補完はちゃんと効かないっぽい? os.py のソース見てみたら 解析が難しそう(動的に読み込んでいる)感じだったから,それが原因かな?

参考

Written by tkf

September 24, 2009 at 3:18 pm

Posted in PC

Tagged with , ,

python の C 拡張をスレッドに対応させる方法 (with callback)

with 3 comments

python の C 拡張をスレッドに対応させつつ,時々は Python の callback 関数を呼びたい時は,グローバルインタプリタロック(GIL) を取得しないといけないですよ,という例.

基本的には 「8.1 スレッド状態 (thread state) とグローバルインタプリタロック (global interpreter lock)」 に書いてある方法.

こんなことをする:
Py_BEGIN_ALLOW_THREADS
...ブロックが起きるような何らかの I/O 操作...
PyGILState_STATE gstate;
gstate = PyGILState_Ensure();
...コールバック呼び出し...
PyGILState_Release(gstate);
Py_END_ALLOW_THREADS

コンパイルに必要なファイル:

callbackthreadmodule.c
 #include <Python.h>

static PyObject *
set_and_call_callback(PyObject *dummy, PyObject *args)
{
  PyObject *temp;
  PyObject *my_callback = NULL;
  PyObject *callback_arg;
  PyObject *callback_result;
  PyGILState_STATE gstate;
  int i;

  if (! PyArg_ParseTuple(args, "OO:set_callback", &temp, &callback_arg)) {
    return NULL;
  }
  if (!PyCallable_Check(temp)) {
    PyErr_SetString(PyExc_TypeError, "parameter must be callable");
    return NULL;
  }
  Py_XINCREF(temp);         /* 新たなコールバックへの参照を追加 */
  Py_XDECREF(my_callback);  /* 以前のコールバックを捨てる */
  my_callback = temp;       /* 新たなコールバックを記憶 */

  Py_BEGIN_ALLOW_THREADS
  for(i=0; i<1000000000; ++i){
    if( !(i % 200000000) ){
      printf("i=%d\n", i);
      gstate = PyGILState_Ensure();
      callback_result = PyEval_CallObject(my_callback, callback_arg);
      if (callback_result == NULL){
        return NULL; /* エラーを返す */
      }
      PyGILState_Release(gstate);
    }
  }
  Py_END_ALLOW_THREADS

  /* "None" を返す際の定型句 */
  Py_INCREF(Py_None);
  return Py_None;
}

static PyMethodDef module_methods&#91;&#93; = {
  {"set_and_call_callback", (PyCFunction)set_and_call_callback,
   METH_VARARGS, "set callback"},
  {NULL}  /* Sentinel */
};

 #ifndef PyMODINIT_FUNC /* declarations for DLL import/export */
 #define PyMODINIT_FUNC void
 #endif
PyMODINIT_FUNC
initcallbackthread(void)
{
  (void) Py_InitModule("callbackthread", module_methods);
}
&#91;/sourcecode&#93;</div>
</div>
<div id="outline-container-12.4.1.2" class="outline-5">
<h5>setup.py</h5>
<div id="text-12.4.1.2">
from distutils.core import setup, Extension

module1 = Extension(
        'callbackthread',
        sources = ['callbackthreadmodule.c']
)

setup(
        name = 'my',
        ext_modules = [module1],
)
コンパイル方法 (その場にcallbackthread.soを作る):

python setup.py config build build_ext –inplace

まずは, thread モジュールを使ってみた

test_thread.py
import callbackthread
import thread

x = (1,2,3)
y = dict(a=2,c=2)

def print_anything(a):
    print a

thread.start_new_thread(callbackthread.set_and_call_callback, (print_anything, (x,)))
thread.start_new_thread(callbackthread.set_and_call_callback, (print_anything, (y,)))
実行
% python test_thread.py
i=0
zsh: segmentation fault  python test_thread.py
% python test_thread.py
i=0
i=0
% python test_thread.py
i=0
i=0
zsh: segmentation fault  python test_thread.py
ん?

確率的に segmentation fault している... そして,時々しないけどなぜだか一回目しかループしてないし. 謎.よく分からなかったけど,threading モジュールだとうまくいった↓

threading モジュールだとうまく行く

test_threading.py
import callbackthread
import threading

x = (1,2,3)
y = dict(a=2,c=2)

def print_anything(a):
    print a

class TestThread(threading.Thread):
    def run(self):
        callbackthread.set_and_call_callback(print_anything, (x,))

TestThread().start()
TestThread().start()
結果
% python test_threading.py
i=0
(1, 2, 3)
i=0
(1, 2, 3)
i=200000000
(1, 2, 3)
i=200000000
(1, 2, 3)
i=400000000
(1, 2, 3)
i=400000000
(1, 2, 3)
i=600000000
(1, 2, 3)
i=600000000
(1, 2, 3)
i=800000000
(1, 2, 3)
i=800000000
(1, 2, 3)
うまくいった!

ちゃんと CPU 二つ使ってることも確かめられた!

どうやら, 「sh1.2 pyblosxom : pythonのthread/threadingモジュールを見てみた」 と同じ問題くさい. しかし segmentation fault って...

やりたかったことは:

test_functor.py
import callbackthread
import threading

x = (1,2,3)
y = dict(a=2,c=2)

class Func():
    i = 0
    def call_back(self, a):
        print "Func.i =",self.i
        print a
        self.i += 1

func = Func()

class TestThread(threading.Thread):
    def run(self):
        callbackthread.set_and_call_callback(func.call_back, (x,))

TestThread().start()
TestThread().start()
実行
% python test_functor.py
i=0
Func.i = 0
(1, 2, 3)
i=0
Func.i = 1
(1, 2, 3)
i=200000000
Func.i = 2
(1, 2, 3)
i=200000000
Func.i = 3
(1, 2, 3)
i=400000000
Func.i = 4
(1, 2, 3)
i=400000000
Func.i = 5
(1, 2, 3)
i=600000000
Func.i = 6
(1, 2, 3)
i=600000000
Func.i = 7
(1, 2, 3)
i=800000000
Func.i = 8
(1, 2, 3)
i=800000000
Func.i = 9
(1, 2, 3)
Func.i ってのを二つのスレッドが交互に書き換えてるのが分かる

python 側のオブジェクトを操作するタイミングを C で決めらる!

Written by tkf

May 11, 2009 at 9:05 pm

Posted in PC, programming

Tagged with ,

python の C 拡張から python の関数を呼ぶ方法の例

with one comment

まんま「1.6 C から Python 関数を呼び出す」にかかれている方法だけど.

C拡張, setup.py, テスト用コード↓

callbackmodule.c

 #include <Python.h>

static PyObject *my_callback = NULL;

static PyObject *
my_set_callback(PyObject *dummy, PyObject *args)
{
  PyObject *result = NULL;
  PyObject *temp;

  if (PyArg_ParseTuple(args, "O:set_callback", &temp)) {
    if (!PyCallable_Check(temp)) {
      PyErr_SetString(PyExc_TypeError, "parameter must be callable");
      return NULL;
    }
    Py_XINCREF(temp);         /* 新たなコールバックへの参照を追加 */
    Py_XDECREF(my_callback);  /* 以前のコールバックを捨てる */
    my_callback = temp;       /* 新たなコールバックを記憶 */
    /* "None" を返す際の定型句 */
    Py_INCREF(Py_None);
    result = Py_None;
  }
  return result;
}

static PyObject *
my_call_callback(PyObject *dummy, PyObject *args)
{
  PyObject *callback_arg;
  PyObject *result;

  if (! PyArg_ParseTuple(args, "O:call_callback", &callback_arg)) {
    return NULL; /* PyArg_ParseTuple has raised an exception */
  }

  /* ここでコールバックを呼ぶ */
  result = PyEval_CallObject(my_callback, callback_arg);
  if (result == NULL){
    return NULL; /* エラーを返す */
  }
  /* 場合によってはここで結果を使うかもね */
  /* Py_DECREF(result); //これやったら,segmentation fault で落ちた */
  return result;
}

static PyObject *
my_call_callback2(PyObject *dummy, PyObject *args)
{
  PyObject *callback_arg;

  if (! PyArg_ParseTuple(args, "O:call_callback", &callback_arg)) {
    return NULL; /* PyArg_ParseTuple has raised an exception */
  }
  return PyEval_CallObject(my_callback, callback_arg);
}

static PyMethodDef module_methods[] = {
  {"my_set_callback", (PyCFunction)my_set_callback, METH_VARARGS, "set callback"},
  {"my_call_callback",(PyCFunction)my_call_callback, METH_VARARGS, "call callback"},
  {"my_call_callback2",(PyCFunction)my_call_callback2, METH_VARARGS, "call callback"},
  {NULL}  /* Sentinel */
};

 #ifndef PyMODINIT_FUNC /* declarations for DLL import/export */
 #define PyMODINIT_FUNC void
 #endif
PyMODINIT_FUNC
initcallback(void)
{
  (void) Py_InitModule("callback", module_methods);
}

メモ

初めに作ったのが my_call_callback で,かなりコードが絞れることが分かったので my_call_callback2 も書いてみた.やることは同じ.

setup.py

from distutils.core import setup, Extension

module1 = Extension(
        'callback',
        sources = ['callbackmodule.c']
)

setup(
        name = 'my',
        ext_modules = [module1],
)

test.py

import callback

x = (1,2,3)
y = dict(a=2,c=2)

callback.my_set_callback(list)
print callback.my_call_callback((x,))
print callback.my_call_callback((y,))

callback.my_set_callback(type)
print callback.my_call_callback((x,))
print callback.my_call_callback((y,))

callback.my_set_callback(list)
print callback.my_call_callback2((x,))
print callback.my_call_callback2((y,))

callback.my_set_callback(type)
print callback.my_call_callback2((x,))
print callback.my_call_callback2((y,))

コンパイル (その場にcallback.soを作る)

python setup.py config build build_ext –inplace

テスト:

% python test.py
[1, 2, 3]
['a', 'c']
<type 'tuple'>
<type 'dict'>
[1, 2, 3]
['a', 'c']
<type 'tuple'>
<type 'dict'>

Written by tkf

May 11, 2009 at 8:28 pm

Posted in PC, programming

Tagged with ,

Gnuplot.py で頑張って pm3d map を使う方法

leave a comment »

Gnuplot.py に入っているオプションは,matrixを受け付けてくれないので ちょっとした(?)ハックが必要だった.かなり試行錯誤したけど.

ソース:

import Gnuplot
import numpy

gnp = Gnuplot.Gnuplot(debug=1)

def plot_mat(mat0):
    # 端の行と列がなくなってしまうので,ダミーを入れる
    (x,y) = mat0.shape
    mat1 = numpy.zeros((x+1, y+1))
    mat1[:x,:y] = mat0
    # splot に 'matrix' オプションを加えるためのhack
    # /usr/lib/python2.5/site-packages/Gnuplot/PlotItems.py を読んだ
    plot = Gnuplot.PlotItems.Data(mat1, using='($2):($1):($3)')
    if plot._option_sequence[0] != 'matrix':
        plot._option_sequence.insert(0,'matrix')
    plot._options['matrix'] = (None,"matrix")
    # plot
    gnp('set pm3d map corners2color c1')
    gnp('set xrange [*:*] reverse')
    gnp.splot(plot)
    gnp('unset pm3d')
    gnp('set xrange [*:*] noreverse')

plot_mat(
    numpy.array([[0, 0  , 0  , 0  , 1],
                 [0, 0.4, 0.1, 0.0, 0],
                 [0, 1.0, 0.6, 0.1, 0],
                 [0, 0.1, 1.0, 0.8, 0],
                 [0, 0  , 0  , 0  , 0]])
    )

実行結果

Written by tkf

April 18, 2009 at 10:48 pm

Posted in 研究日誌

Tagged with ,