Boost.Python debugging

Boost.Python で std::valarray をエクスポートする際、__getitem__ に operator[] を渡すと、インデックスが境界を超えてもエラーにならないため、無限ループに陥る現象について書きました。その後の追加情報を紹介します。

1. 非メンバ関数をクラス メソッドとしてエクスポート

前回は、回避策として valarray を継承したクラスを作って operator::at を新たに定義しました。が、Python の wiki を見ていたところ、Boost.Python の class_ で定義するコンバーターのメンバー関数にはメンバー関数ポインターだけでなく、通常の関数も指定できることが分かりました。

boost.python/ExportingClasses – Python Wiki
https://wiki.python.org/moin/boost.python/ExportingClasses

これを使うと、新たにクラスを定義しなくても __getitem__ を正しく定義できます。

#include <boost/python.hpp>
#include <valarray>
#include <random>
#include <stdio.h>

std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(0.0, 1.0);

typedef std::valarray<double> darray;

class Field {
private:
    darray _data;

public:
    Field(int n) : _data(n) {
        printf("Field::Field(): %p\n", this);
    }
    ~Field() {
        printf("Field::~Field() %p\n", this);
    }

    void Dump() const {
        printf("[");
        for (auto &d : _data) {
            printf("% 3.3f", d);
        }
        printf(" ]\n");
    }

    void Churn() {
        for (auto &d : _data) {
            d = distribution(generator);
        }
    }

    const darray &Positions() const {
        return _data;
    }
};

double GetAt(const darray &arr, size_t n) {
    if (n >= arr.size()) {
        std::__throw_out_of_range_fmt(__N("out_of_range"));
    }
    return arr[n];
}

BOOST_PYTHON_MODULE(ems) {
    using namespace boost::python;

    class_<Field>("Field", init<int>())
        .def("Positions", &Field::Positions,
             return_internal_reference<>())
        .def("Dump", &Field::Dump)
        .def("Churn", &Field::Churn)
        ;

    class_<darray>("darray")
        .def("__getitem__", GetAt)
        .def("__len__", &darray::size)
        ;
}

ここから、壮大に話が逸れていくのですが、前回こんなことを書きました。

もう少し見てみないと分かりませんが、ems.Field オブジェクトが iterator を使ってシーケンシャルにアクセスする種類になっているのが悪いのではないかと思います。何らかの方法で、ems.Field は iterator ではなくインデックスを使ってランダム アクセスすべきものに種別できればちゃんと動くのではないかと。そんな種類があるかどうか分かりませんが。後で調べよう。

今回も結論から書くとたぶんこれは無理です。Boost.Python をデバッグしつつ調べたので、その過程を書いておきます。

2. Boost.Python をソースからビルドして使う

まずデバッグ環境の構築です。Boost のコードは GitHub にあるので、クローンしてからビルドします。

git clone –recursive https://github.com/boostorg/boost.git modular-boost
git checkout refs/tags/boost-1.60.0
cd modular-boost/
./bootstrap.sh
./b2 -a –with-python debug-symbols=on
sudo ./b2 -a –prefix=/usr/local/boost/boost-1.60.0 –with-python debug-symbols=on install
sudo ln -s /usr/local/boost/boost-1.60.0 /usr/local/boost/current

コマンドを作るうえで参考にしたページは↓。Boost をビルドしている人は少ない気がする・・。

TryModBoost – Boost C++ Libraries
https://svn.boost.org/trac/boost/wiki/TryModBoost

Installing Boost.Python on your System
http://boostorg.github.io/python/doc/html/building/installing_boost_python_on_your_.html

Builtin features
http://www.boost.org/build/doc/html/bbv2/overview/builtins/features.html

ビルドには、make ではなく b2 (build boost?) を使います。最新のリポジトリをクローンしてそのままビルドしたところ、文法エラーでビルドが失敗したので、boost-1.60.0 のタグの状態でビルドしました。もう直っているかもしれません。

デバッグ情報を生成するため、debug-symbols=on というオプションをつけています。これは、コードの最適化はデフォルトの ON にしたまま、デバッグ情報を生成するという意味です。普通は variant=debug というオプションをつけて、最適化も off にしたほうがいいです。インライン展開された魔境に挑みたい人のみ、debug-symbols=on を使いましょう。

上記コマンドで、libboost_python が /usr/local/boost/current/lib にできるので、C++ の Makefile を以下のように変更します。Boost,Python は Python API を使うので、API の中をデバッグするため、libpython2.7 のデバッグ情報もあったほうが便利です。以下の Makefile は、ソースからビルドした Python 2.7 が /usr/local/python/current/ にインストールされている前提です。Boost.Python を使って共有ライブラリを作るので、libpython2.7 も PIC でコンパイルする必要があります。つまり、Python の configure で –enable-shared オプションをつけておいてください。

CC=g++
RM=rm -f

TARGET=ems.so
SRCS=$(wildcard *.cpp)
OBJS=$(SRCS:.cpp=.o)

override CFLAGS+=-Wall -fPIC -std=c++11 -O2 -g
LFLAGS=

INCLUDES=-I/usr/local/python/current/include/python2.7 -I/usr/local/boost/current/include
LIBDIRS=-L/usr/local/python/current/lib -L/usr/local/boost/current/lib
LIBS=-lpython2.7 -lboost_python

all: clean $(TARGET)

clean:
        $(RM) $(OBJS) $(TARGET)

$(TARGET): $(OBJS)
        $(CC) -shared $(LFLAGS) $(LIBDIRS) $^ -o $@ $(LIBS)

$(OBJS): $(SRCS)
        $(CC) $(INCLUDES) $(CFLAGS) -c $^

デバッガーを起動する前に、python を実行する前のライブラリの検索パスにも追加しておきます。gdb を起動するときはこんな感じにします。

export PYTHONPATH=/usr/local/lib/python2.7/dist-packages
export LD_LIBRARY_PATH=/usr/local/boost/current/lib:/usr/local/python/current/lib
gdb /usr/local/python/current/bin/python

ここから gdb によるデバッグを行いますが、結果として何か得られたわけではありません。考え方の参考になれば幸いです。

まず、Python API に PyObject_GetIter という関数があります。

Object Protocol — Python 2.7.11 documentation
https://docs.python.org/2/c-api/object.html

valarray のコンバーターである darray に対して Python から列挙を行なおうとすると、PyObject_GetIter が呼ばれます。

john@ubuntu1510:~/Documents/pyc$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages
john@ubuntu1510:~/Documents/pyc$ export LD_LIBRARY_PATH=/usr/local/boost/current/lib:/usr/local/python/current/lib
john@ubuntu1510:~/Documents/pyc$ gdb /usr/local/python/current/bin/python
GNU gdb (Ubuntu 7.10-1ubuntu2) 7.10
Copyright (C) 2015 Free Software Foundation, Inc.
(略)
Type "apropos word" to search for commands related to "word"…
Reading symbols from /usr/local/python/current/bin/python…done.
(gdb) r
Starting program: /usr/local/python/python-2.7.11/bin/python
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
Python 2.7.11 (default, Dec 22 2015, 22:30:16)
[GCC 5.2.1 20151010] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import ems
>>> f = ems.Field(4)
Field::Field(): 0x6d09d0
>>> x = f.Positions()
>>> x
<ems.darray object at 0x7ffff7ed5600>
>>>
Program received signal SIGINT, Interrupt.
0x00007ffff74e4723 in __select_nocancel () at ../sysdeps/unix/syscall-template.S:81
81      ../sysdeps/unix/syscall-template.S: No such file or directory.
(gdb) b PyObject_GetIter
Breakpoint 1 at 0x7ffff7a261c0: file Objects/abstract.c, line 3085.
(gdb) i b
Num     Type           Disp Enb Address            What
1       breakpoint     keep y   0x00007ffff7a261c0 in PyObject_GetIter at Objects/abstract.c:3085
(gdb) c
Continuing.

>>> [i for i in x]

Breakpoint 1, PyObject_GetIter (o=o@entry=0x7ffff7ed5600) at Objects/abstract.c:3085
3085    {
(gdb) bt 5
#0  PyObject_GetIter (o=o@entry=0x7ffff7ed5600) at Objects/abstract.c:3085
#1  0x00007ffff7ad7db0 in PyEval_EvalFrameEx (f=f@entry=0x7ffff7f7e050,
    throwflag=throwflag@entry=0) at Python/ceval.c:2790
#2  0x00007ffff7ae068c in PyEval_EvalCodeEx (co=co@entry=0x7ffff7ed6930,
    globals=globals@entry=0x7ffff7f6b168, locals=locals@entry=0x7ffff7f6b168, args=args@entry=0x0,
    argcount=argcount@entry=0, kws=kws@entry=0x0, kwcount=0, defs=0x0, defcount=0, closure=0x0)
    at Python/ceval.c:3582
#3  0x00007ffff7ae07a9 in PyEval_EvalCode (co=co@entry=0x7ffff7ed6930,
    globals=globals@entry=0x7ffff7f6b168, locals=locals@entry=0x7ffff7f6b168) at Python/ceval.c:669
#4  0x00007ffff7b05d81 in run_mod (arena=0x66f100, flags=0x7fffffffe370, locals=0x7ffff7f6b168,
    globals=0x7ffff7f6b168, filename=0x7fffffffe370 "", mod=0x6b9ce8) at Python/pythonrun.c:1370
(More stack frames follow…)

関数のコードはこうなっています。引数として与えられた Python オブジェクトに関連付けられた Iterator をオブジェクトとして返す関数です。

Objects/abstract.c
3083 PyObject *
3084 PyObject_GetIter(PyObject *o)
3085 {
3086     PyTypeObject *t = o->ob_type;
3087     getiterfunc f = NULL;
3088     if (PyType_HasFeature(t, Py_TPFLAGS_HAVE_ITER))
3089         f = t->tp_iter;
3090     if (f == NULL) {
3091         if (PySequence_Check(o))
3092             return PySeqIter_New(o);
3093         return type_error("’%.200s’ object is not iterable", o);
3094     }
3095     else {
3096         PyObject *res = (*f)(o);
3097         if (res != NULL && !PyIter_Check(res)) {
3098             PyErr_Format(PyExc_TypeError,
3099                          "iter() returned non-iterator "
3100                          "of type ‘%.100s’",
3101                          res->ob_type->tp_name);

関数の冒頭で、引数の PyTypeObject で Py_TPFLAGS_HAVE_ITER フラグを判定し、有効な場合は関連付けられた iterator である tp_iter を使うようになっています。これをパッと見たときに、valarray からコンバートした darray オブジェクトからこのフラグを削除してみようと考えました。まずは、現時点での TypeObject を見てます。

(gdb) frame
#0  PyObject_GetIter (o=o@entry=0x7ffff7ed5600) at Objects/abstract.c:3085
3085    {
(gdb) p *o->ob_type
$1 = {ob_refcnt = 5, ob_type = 0x7ffff5a7cc40 <boost::python::class_metatype_object>, ob_size = 0,
  tp_name = 0x7ffff7e94954 "darray", tp_basicsize = 48, tp_itemsize = 1,
  tp_dealloc = 0x7ffff7a8b390 <subtype_dealloc>, tp_print = 0x0, tp_getattr = 0x0,
  tp_setattr = 0x0, tp_compare = 0x0, tp_repr = 0x7ffff7a8f390 <object_repr>,
  tp_as_number = 0x6d3c88, tp_as_sequence = 0x6d3dd8, tp_as_mapping = 0x6d3dc0,
  tp_hash = 0x7ffff7a70be0 <_Py_HashPointer>, tp_call = 0x0, tp_str = 0x7ffff7a8acc0 <object_str>,
  tp_getattro = 0x7ffff7a72410 <PyObject_GenericGetAttr>,
  tp_setattro = 0x7ffff7a72670 <PyObject_GenericSetAttr>, tp_as_buffer = 0x6d3e28,
  tp_flags = 153595, tp_doc = 0x0, tp_traverse = 0x7ffff7a8b0a0 <subtype_traverse>,
  tp_clear = 0x7ffff7a8d290 <subtype_clear>, tp_richcompare = 0x0, tp_weaklistoffset = 32,
  tp_iter = 0x0, tp_iternext = 0x7ffff7a72100 <_PyObject_NextNotImplemented>, tp_methods = 0x0,
  tp_members = 0x6d3e68, tp_getset = 0x0,
  tp_base = 0x7ffff5a7cfa0 <boost::python::objects::class_type_object>, tp_dict = 0x7ffff7e98050,
  tp_descr_get = 0x0, tp_descr_set = 0x0, tp_dictoffset = 24,
  tp_init = 0x7ffff7a91f70 <slot_tp_init>, tp_alloc = 0x7ffff7a8ade0 <PyType_GenericAlloc>,
  tp_new = 0x7ffff5856890
     <boost::python::objects::instance_new(PyTypeObject*, PyObject*, PyObject*)>,
  tp_free = 0x7ffff7b1e2b0 <PyObject_GC_Del>, tp_is_gc = 0x0, tp_bases = 0x7ffff7e96410,
  tp_mro = 0x7ffff7e8e780, tp_cache = 0x0, tp_subclasses = 0x0, tp_weaklist = 0x7ffff7e77f70,
  tp_del = 0x0, tp_version_tag = 0}
(gdb) p/x 153595
$2 = 0x257fb

型の名前は "darray" となっており、Py_TPFLAGS_HAVE_ITER (= 0x80) は ON、tp_iter は NULL になっています。気づくまで時間がかかってしまったのですが、この TypeObject の内容と PyObject_GetIter のコードを見る限り、Py_TPFLAGS_HAVE_ITER の有無による影響はなく、いずれの場合でも f == NULL は true になります。そこで PySequence_Check の定義を見ます。

Objects/abstract.c
1843 int
1844 PySequence_Check(PyObject *s)
1845 {
1846     if (s == NULL)
1847         return 0;
1848     if (PyInstance_Check(s))
1849         return PyObject_HasAttrString(s, "__getitem__");
1850     if (PyDict_Check(s))
1851         return 0;
1852     return  s->ob_type->tp_as_sequence &&
1853         s->ob_type->tp_as_sequence->sq_item != NULL;
1854 }

darray オブジェクトは、3 つの if 文には引っかかりませんが、最後の条件で true になります。sq_item の値は slot_sq_item 関数になっています。slot_sq_item 関数は前回のデバッグ時に、iter_iternext から呼び出されており、これらが境界チェックをしていないのが無限ループの原因でした。

(gdb) p *$1.tp_as_sequence
$8 = {sq_length = 0x7ffff7a93f60 <slot_sq_length>, sq_concat = 0x0, sq_repeat = 0x0,
  sq_item = 0x7ffff7a91090 <slot_sq_item>, sq_slice = 0x0, sq_ass_item = 0x0, sq_ass_slice = 0x0,
  sq_contains = 0x0, sq_inplace_concat = 0x0, sq_inplace_repeat = 0x0}

そんなわけで PySequence_Check は true を返すので、PyObject_GetIter から PySeqIter_New が呼ばれて新しい iterator が作られます。その iterator が↓です。gdb に windbg で言うところの ub コマンドが無いのが辛い・・。

(gdb) x/5i 0x7ffff7ad7da1
   0x7ffff7ad7da1 <PyEval_EvalFrameEx+4225>:    mov    %rcx,%rbp
   0x7ffff7ad7da4 <PyEval_EvalFrameEx+4228>:    mov    -0x8(%rbx),%r13
   0x7ffff7ad7da8 <PyEval_EvalFrameEx+4232>:    mov    %r13,%rdi
   0x7ffff7ad7dab <PyEval_EvalFrameEx+4235>:    callq  0x7ffff7a0f550 <PyObject_GetIter@plt>
=> 0x7ffff7ad7db0 <PyEval_EvalFrameEx+4240>:    subq   $0x1,0x0(%r13)
(gdb) p *(seqiterobject*)$rax
$9 = {ob_refcnt = 1, ob_type = 0x7ffff7d9e9c0 <PySeqIter_Type>, it_index = 0,
  it_seq = 0x7ffff7ed5600}

つまり、__getitem__ だけを実装したオブジェクトはシーケンスとして扱われるが、Python はシーケンスもリストも iterator を使ってアクセスし、それは境界をいちいちチェックしない、という仕様になっています。つまりわりとどうしようもないです。冒頭で書いたように、クラス定義をしなくても通常関数をオブジェクト メソッドとして実装できるので、それが正しい解決策なのだと思います。

Boost のデバッグいつ出てくるんだ、という話ですが、Python オブジェクトから Py_TPFLAGS_HAVE_ITER フラグを消す方法を見つけるためにデバッグしました。

3. Boost.Python のコードをデバッグする

上述の o->ob_type が指すオブジェクトがいつ作られているのか、を追いかけます。Boost や Python のコードを見ても、Py_TPFLAGS_HAVE_ITER を個別に設定しているコードは無く、デフォルトの Py_TPFLAGS_DEFAULT で使われているだけです。darray の TypeObject が作られるヒントとしては、TypeObject の ob_type が ob_type = 0x7ffff5a7cc40 <boost::python::class_metatype_object> になっているのが使えます。

class_metatype_object は、Boost.Python のグローバル変数であり、これが使われる場所は一箇所しかありません。それが class_metatype() です。

src/object/class.cpp
315   BOOST_PYTHON_DECL type_handle class_metatype()
316   {
317       if (class_metatype_object.tp_dict == 0)
318       {
319           Py_TYPE(&class_metatype_object) = &PyType_Type;
320           class_metatype_object.tp_base = &PyType_Type;
321           if (PyType_Ready(&class_metatype_object))
322               return type_handle();
323       }
324       return type_handle(borrowed(&class_metatype_object));
325   }

幸い、この関数はインライン展開されないので、ここを起点にします。ポイントは、この関数はモジュールのインポート時に呼ばれることです。そこで、import 後の状態で関数の完全修飾名を調べてから、再度デバッグ ターゲットを実行します。

john@ubuntu1510:~/Documents/pyc$ gdb /usr/local/python/current/bin/python
GNU gdb (Ubuntu 7.10-1ubuntu2) 7.10
(略)
Reading symbols from /usr/local/python/current/bin/python…done.
(gdb) r
Starting program: /usr/local/python/python-2.7.11/bin/python
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
Python 2.7.11 (default, Dec 22 2015, 22:30:16)
[GCC 5.2.1 20151010] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import ems # まずは何もせずインポート
>>>
Program received signal SIGINT, Interrupt.
0x00007ffff74e4723 in __select_nocancel () at ../sysdeps/unix/syscall-template.S:81
81      ../sysdeps/unix/syscall-template.S: No such file or directory.
(gdb) i func class_metatype # 完全修飾名を調べる
All functions matching regular expression "class_metatype":

File libs/python/src/object/class.cpp:
boost::python::type_handle boost::python::objects::class_metatype();

Non-debugging symbols:
0x00007ffff5845c40  boost::python::objects::class_metatype()@plt
(gdb) r # 再実行
The program being debugged has been started already.
Start it from the beginning? (y or n) y
Starting program: /usr/local/python/python-2.7.11/bin/python
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
Python 2.7.11 (default, Dec 22 2015, 22:30:16)
[GCC 5.2.1 20151010] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>>
Program received signal SIGINT, Interrupt.
0x00007ffff74e4723 in __select_nocancel () at ../sysdeps/unix/syscall-template.S:81
81      ../sysdeps/unix/syscall-template.S: No such file or directory.
(gdb) b boost::python::objects::class_metatype # import 前にブレークポイント
Function "boost::python::objects::class_metatype" not defined.
Make breakpoint pending on future shared library load? (y or [n]) y
Breakpoint 1 (boost::python::objects::class_metatype) pending.
(gdb) c
Continuing.

>>> import ems

Breakpoint 1, 0x00007ffff5845c40 in boost::python::objects::class_metatype()@plt ()
   from /usr/local/boost/current/lib/libboost_python.so.1.60.0
(gdb) i b
Num     Type           Disp Enb Address            What
1       breakpoint     keep y   <MULTIPLE>
        breakpoint already hit 1 time
1.1                         y     0x00007ffff5845c40 <boost::python::objects::class_metatype()@plt>
1.2                         y     0x00007ffff5856b00 in boost::python::objects::class_metatype()
                                                   at libs/python/src/object/class.cpp:317
(gdb) disable 1.1 # plt 上のブレークポイントは無効
(gdb) c
Continuing.

Breakpoint 1, boost::python::objects::class_metatype () at libs/python/src/object/class.cpp:317
317           if (class_metatype_object.tp_dict == 0)
(gdb)

もう一点。この方法だと、モジュールがロードされたときにシンボル名をもとにアドレスをに検索するので、対象が共有ライブラリの場合、関数の先頭のアドレス以外に、plt 上のアドレスにもブレークポイントが設定されます。これは邪魔なので無効にしておきます。

PLT とは Procedure Linkage Table の略で、Windows でいうところのインポート テーブルと似たようなものです。

Technovelty – PLT and GOT – the key to code sharing and dynamic libraries
https://www.technovelty.org/linux/plt-and-got-the-key-to-code-sharing-and-dynamic-libraries.html

コールスタックを見ると、shared.cpp に実装された init_module_ems から、class_ クラスのコンストラクタ (#4) 経由で呼ばれています。テンプレートの入れ子で、シンボル名が複雑怪奇です。

(gdb) bt 10
#0  boost::python::objects::class_metatype () at libs/python/src/object/class.cpp:317
#1  0x00007ffff5856bb8 in boost::python::objects::class_type ()
    at libs/python/src/object/class.cpp:473
#2  0x00007ffff5857d8d in boost::python::objects::(anonymous namespace)::new_class (
    doc=<optimized out>, types=0x7fffffffdd90, num_types=1, name=<optimized out>)
    at libs/python/src/object/class.cpp:561
#3  boost::python::objects::class_base::class_base (this=0x7fffffffdbf0,
    name=0x7ffff5a89bd5 "Field", num_types=1, types=0x7fffffffdd90, doc=0x0)
    at libs/python/src/object/class.cpp:591
#4  0x00007ffff5a874be in boost::python::class_<Field, boost::python::detail::not_specified, boost::python::detail::not_specified, boost::python::detail::not_specified>::class_<boost::python::init<int, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_, mpl_::void_> > (i=…,
    name=0x7ffff5a89bd5 "Field", this=0x7fffffffdbf0)
    at /usr/local/boost/current/include/boost/python/class.hpp:206
#5  init_module_ems () at shared.cpp:51
#6  0x00007ffff58639f3 in boost::function0<void>::operator() (this=0x7fffffffde50)
    at ./boost/function/function_template.hpp:771
#7  boost::python::handle_exception_impl (f=…) at libs/python/src/errors.cpp:25
#8  0x00007ffff5864736 in boost::python::handle_exception<void (*)()> (
    f=0x7ffff5a87440 <init_module_ems()>) at ./boost/python/errors.hpp:29
#9  boost::python::detail::(anonymous namespace)::init_module_in_scope (
    init_function=<optimized out>, m=0x7ffff7e937c0) at libs/python/src/module.cpp:24
(More stack frames follow…)

#3 の class_base のコンストラクタは以下のように定義されています。

src/object/class.cpp
589   class_base::class_base(
590       char const* name, std::size_t num_types, type_info const* const types, char const* doc)
591       : object(new_class(name, num_types, types, doc))
592   {
593       // Insert the new class object in the registry
594       converter::registration& converters = const_cast<converter::registration&>(
595           converter::registry::lookup(types[0]));
596
597       // Class object is leaked, for now
598       converters.m_class_object = (PyTypeObject*)incref(this->ptr());
599   }

class_metatype() は初期化リストのなかで使われています。ここでの object とは、class_base の基底クラスの一つです。object クラスを new_class 関数の戻り値で初期化すると、m_ptr というメンバー関数が初期化されます。引数の types を元に lookup してきたコンバーターの m_class_object に、その m_ptr を代入する、というのがこのコンストラクターの動作です。

ブレークポイントを設定した class_metatype() は、new_class() から 2 回呼ばれます。line:561 の class_type() 経由と、line:575 です。戻り値に直接関係あるのは後者です。

src/object/class.cpp
548     inline object
549     new_class(char const* name, std::size_t num_types, type_info const* const types, char const*     doc)
550     {
551       assert(num_types >= 1);
552
553       // Build a tuple of the base Python type objects. If no bases
554       // were declared, we’ll use our class_type() as the single base
555       // class.
556       ssize_t const num_bases = (std::max)(num_types – 1, static_cast<std::size_t>(1));
557       handle<> bases(PyTuple_New(num_bases));
558
559       for (ssize_t i = 1; i <= num_bases; ++i)
560       {
561           type_handle c = (i >= static_cast<ssize_t>(num_types)) ? class_type() : get_class(type    s[i]);
562           // PyTuple_SET_ITEM steals this reference
563           PyTuple_SET_ITEM(bases.get(), static_cast<ssize_t>(i – 1), upcast<PyObject>(c.release(    )));
564       }
565
566       // Call the class metatype to create a new class
567       dict d;
568
569       object m = module_prefix();
570       if (m) d["__module__"] = m;
571
572       if (doc != 0)
573           d["__doc__"] = doc;
574
575       object result = object(class_metatype())(name, bases, d);
576       assert(PyType_IsSubtype(Py_TYPE(result.ptr()), &PyType_Type));
577
578       if (scope().ptr() != Py_None)
579           scope().attr(name) = result;
580
581       // For pickle. Will lead to informative error messages if pickling
582       // is not enabled.
583       result.attr("__reduce__") = object(make_instance_reduce_function());
584
585       return result;
586     }
587   }

1 回目の呼び出しはスルーして、2 回目の呼び出しで止めます。

(gdb) b boost::python::objects::class_metatype
Function "boost::python::objects::class_metatype" not defined.
Make breakpoint pending on future shared library load? (y or [n]) y
Breakpoint 1 (boost::python::objects::class_metatype) pending.
(gdb) r
Starting program: /usr/local/python/python-2.7.11/bin/python
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
Python 2.7.11 (default, Dec 22 2015, 22:30:16)
[GCC 5.2.1 20151010] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import ems

Breakpoint 1, 0x00007ffff5845c40 in boost::python::objects::class_metatype()@plt ()
   from /usr/local/boost/current/lib/libboost_python.so.1.60.0
(gdb) i b
Num     Type           Disp Enb Address            What
1       breakpoint     keep y   <MULTIPLE>
        breakpoint already hit 1 time
1.1                         y     0x00007ffff5845c40 <boost::python::objects::class_metatype()@plt>
1.2                         y     0x00007ffff5856b00 in boost::python::objects::class_metatype()
                                                   at libs/python/src/object/class.cpp:317
(gdb) disable 1.1
(gdb) c
Continuing.

Breakpoint 1, boost::python::objects::class_metatype () at libs/python/src/object/class.cpp:317
317           if (class_metatype_object.tp_dict == 0)
(gdb) bt 4
#0  boost::python::objects::class_metatype () at libs/python/src/object/class.cpp:317
#1  0x00007ffff5856bb8 in boost::python::objects::class_type ()
    at libs/python/src/object/class.cpp:473
#2  0x00007ffff5857d8d in boost::python::objects::(anonymous namespace)::new_class (
    doc=<optimized out>, types=0x7fffffffdd90, num_types=1, name=<optimized out>)
    at libs/python/src/object/class.cpp:561
#3  boost::python::objects::class_base::class_base (this=0x7fffffffdbf0,
    name=0x7ffff5a89bd5 "Field", num_types=1, types=0x7fffffffdd90, doc=0x0)
    at libs/python/src/object/class.cpp:591
(More stack frames follow…)
(gdb) c
Continuing.

Breakpoint 1, boost::python::objects::class_metatype () at libs/python/src/object/class.cpp:317
317           if (class_metatype_object.tp_dict == 0)
(gdb) bt 3
#0  boost::python::objects::class_metatype () at libs/python/src/object/class.cpp:317
#1  0x00007ffff5857f92 in boost::python::objects::(anonymous namespace)::new_class (
    doc=<optimized out>, types=0x7fffffffdd90, num_types=1, name=<optimized out>)
    at libs/python/src/object/class.cpp:575
#2  boost::python::objects::class_base::class_base (this=0x7fffffffdbf0,
    name=0x7ffff5a89bd5 "Field", num_types=1, types=0x7fffffffdd90, doc=0x0)
    at libs/python/src/object/class.cpp:591
(More stack frames follow…)
(gdb)

line:575 の "object(class_metatype())(name, bases, d);" で、object が関数なのかコンストラクタなのか理解に苦しむところです。後半の (name, bases, d) は operator() 呼び出しだと推測は出来ます。この行は最適化が有効だと全部インライン展開されるので、とりあえず class_metatype() から戻ってきたところでアセンブリを見ます。

(gdb) fin
Run till exit from #0  boost::python::objects::class_metatype ()
    at libs/python/src/object/class.cpp:317
boost::python::objects::(anonymous namespace)::new_class (doc=<optimized out>,
    types=0x7fffffffdd90, num_types=1, name=<optimized out>)
    at libs/python/src/object/class.cpp:575
Value returned is $1 = {m_p = 0x7fffffffdaa0}
(gdb) x/20i $rip
=> 0x7ffff5857f92 <…+690>:   mov    0x60(%rsp),%rbx
   0x7ffff5857f97 <…+695>:   test   %rbx,%rbx
   0x7ffff5857f9a <…+698>:   je     0x7ffff5858582 <…+2210>
   0x7ffff5857fa0 <…+704>:   addq   $0x1,(%rbx)
   0x7ffff5857fa4 <…+708>:   test   %rbp,%rbp
   0x7ffff5857fa7 <…+711>:   mov    0x40(%rsp),%r14
   0x7ffff5857fac <…+716>:   mov    %rbp,%r13
   0x7ffff5857faf <…+719>:   je     0x7ffff5858576 <…+2198>
   0x7ffff5857fb5 <…+725>:   mov    0x18(%rsp),%rdi
   0x7ffff5857fba <…+730>:   callq  0x7ffff5844d40 <_ZN5boost6python9converter19do_return_to_pythonEPKc@plt>
   0x7ffff5857fbf <…+735>:   test   %rax,%rax
   0x7ffff5857fc2 <…+738>:   mov    %rax,%r12
   0x7ffff5857fc5 <…+741>:   je     0x7ffff58585bd <…+2269>
   0x7ffff5857fcb <…+747>:   lea    0x1705a(%rip),%rsi        # 0x7ffff586f02c
   0x7ffff5857fd2 <…+754>:   mov    %r14,%r8
   0x7ffff5857fd5 <…+757>:   mov    %r13,%rcx
   0x7ffff5857fd8 <…+760>:   mov    %r12,%rdx
   0x7ffff5857fdb <…+763>:   mov    %rbx,%rdi
   0x7ffff5857fde <…+766>:   xor    %eax,%eax
   0x7ffff5857fe0 <…+768>:   callq  0x7ffff5845a30 <PyEval_CallFunction@plt>
(gdb)

PyEval_CallFunction という Python API の呼び出しが怪しいのでここで止めます。

(gdb) b PyEval_CallFunction
Breakpoint 2 at 0x7ffff7b00fd0: file Python/modsupport.c, line 544.
(gdb) c
Continuing.

Breakpoint 2, PyEval_CallFunction (
    obj=obj@entry=0x7ffff5a7cc40 <boost::python::class_metatype_object>,
    format=format@entry=0x7ffff586f02c "(OOO)") at Python/modsupport.c:544
544     {
(gdb) bt 5
#0  PyEval_CallFunction (obj=obj@entry=0x7ffff5a7cc40 <boost::python::class_metatype_object>,
    format=format@entry=0x7ffff586f02c "(OOO)") at Python/modsupport.c:544

#1  0x00007ffff5857fe5 in boost::python::call<boost::python::api::object, char const*, boost::python::handle<_object>, boost::python::dict> (a2=…, a1=<synthetic pointer>, a0=<synthetic pointer>,
    callable=0x7ffff5a7cc40 <boost::python::class_metatype_object>) at ./boost/python/call.hpp:66
#2  boost::python::api::object_operators<boost::python::api::object>::operator()<char const*, boost::python::handle<_object>, boost::python::dict> (a2=…, a1=<synthetic pointer>,
    a0=<synthetic pointer>, this=<optimized out>) at ./boost/python/object_call.hpp:19
#3  boost::python::objects::(anonymous namespace)::new_class (doc=<optimized out>,
    types=0x7fffffffdd90, num_types=<optimized out>, name=<optimized out>)
    at libs/python/src/object/class.cpp:575
#4  boost::python::objects::class_base::class_base (this=0x7fffffffdbf0,
    name=0x7ffff5a89bd5 "Field", num_types=<optimized out>, types=0x7fffffffdd90, doc=0x0)
    at libs/python/src/object/class.cpp:591
(More stack frames follow…)
(gdb)

PyEval_CallFunction は、C/C++ から Python の関数を呼ぶための API です。第二引数の format が "(OOO)" になっているので、3 つのオブジェクトからなるタプルを引数として Python の関数を呼んでいます。何の関数を呼ぶのかは、第一引数を見れば分かります。

(gdb) p *obj->ob_type
$4 = {ob_refcnt = 41, ob_type = 0x7ffff7daaf60 <PyType_Type>, ob_size = 0,
  tp_name = 0x7ffff7b4008c "type", tp_basicsize = 872, tp_itemsize = 40,
  tp_dealloc = 0x7ffff7a8b1f0 <type_dealloc>, tp_print = 0x0, tp_getattr = 0x0, tp_setattr = 0x0,
  tp_compare = 0x0, tp_repr = 0x7ffff7a8f600 <type_repr>, tp_as_number = 0x0,
  tp_as_sequence = 0x0, tp_as_mapping = 0x0, tp_hash = 0x7ffff7a70be0 <_Py_HashPointer>,
  tp_call = 0x7ffff7a908a0 <type_call>, tp_str = 0x7ffff7a8acc0 <object_str>,
  tp_getattro = 0x7ffff7a9a2c0 <type_getattro>, tp_setattro = 0x7ffff7a91930 <type_setattro>,
  tp_as_buffer = 0x0, tp_flags = 2148423147,
  tp_doc = 0x7ffff7da8f80 <type_doc> "type(object) -> the object’s type\ntype(name, bases, dict) -> a new type", tp_traverse = 0x7ffff7a8c2c0 <type_traverse>, tp_clear = 0x7ffff7a90370 <type_clear>,
  tp_richcompare = 0x7ffff7a8bac0 <type_richcompare>, tp_weaklistoffset = 368, tp_iter = 0x0,
  tp_iternext = 0x0, tp_methods = 0x7ffff7da9200 <type_methods>,
  tp_members = 0x7ffff7da9500 <type_members>, tp_getset = 0x7ffff7da93e0 <type_getsets>,
  tp_base = 0x7ffff7daadc0 <PyBaseObject_Type>, tp_dict = 0x7ffff7f98280, tp_descr_get = 0x0,
  tp_descr_set = 0x0, tp_dictoffset = 264, tp_init = 0x7ffff7a8d3b0 <type_init>,
  tp_alloc = 0x7ffff7a8ade0 <PyType_GenericAlloc>, tp_new = 0x7ffff7a98200 <type_new>,
  tp_free = 0x7ffff7b1e2b0 <PyObject_GC_Del>, tp_is_gc = 0x7ffff7a8aca0 <type_is_gc>,
  tp_bases = 0x7ffff7f9b090, tp_mro = 0x7ffff7f9a878, tp_cache = 0x0,
  tp_subclasses = 0x7ffff7f4ce60, tp_weaklist = 0x7ffff7f9e050, tp_del = 0x0, tp_version_tag = 11}

名前と doc から、これは Python の組み込み関数 type() であることが分かります。どうやら Boost.Python は、Python の組み込み関数を使って型を作っているようです。テンプレートでこんなことができるとか Boost 凄すぎ・・。

Python のコードの中に、新しい型オブジェクトを作る関数があるはずなので、それを適当なキーワード (Py_TPFLAGS_HAVE_GC とか) を使って探すと、type_new というそれっぽい関数が Objects/typeobject.c に見つかるので、ここで止めます。

(gdb) i b
Num     Type           Disp Enb Address            What
1       breakpoint     keep y   <MULTIPLE>
        breakpoint already hit 3 times
1.1                         n     0x00007ffff5845c40 <boost::python::objects::class_metatype()@plt>
1.2                         y     0x00007ffff5856b00 in boost::python::objects::class_metatype()
                                                   at libs/python/src/object/class.cpp:317
2       breakpoint     keep y   0x00007ffff7b00fd0 in PyEval_CallFunction
                                                   at Python/modsupport.c:544
        breakpoint already hit 1 time
(gdb) disable
(gdb) b type_new
Breakpoint 3 at 0x7ffff7a98200: file Objects/typeobject.c, line 2068.
(gdb) c
Continuing.

Breakpoint 3, type_new (metatype=0x7ffff5a7cc40 <boost::python::class_metatype_object>,
    args=0x7ffff7e8e370, kwds=0x0) at Objects/typeobject.c:2068
2068    {
(gdb) bt 8
#0  type_new (metatype=0x7ffff5a7cc40 <boost::python::class_metatype_object>, args=0x7ffff7e8e370,
    kwds=0x0) at Objects/typeobject.c:2068
#1  0x00007ffff7a908c3 in type_call (type=0x7ffff5a7cc40 <boost::python::class_metatype_object>,
    args=0x7ffff7e8e370, kwds=0x0) at Objects/typeobject.c:729
#2  0x00007ffff7a244e3 in PyObject_Call (
    func=func@entry=0x7ffff5a7cc40 <boost::python::class_metatype_object>,
    arg=arg@entry=0x7ffff7e8e370, kw=<optimized out>) at Objects/abstract.c:2546
#3  0x00007ffff7ad6707 in PyEval_CallObjectWithKeywords (
    func=func@entry=0x7ffff5a7cc40 <boost::python::class_metatype_object>,
    arg=arg@entry=0x7ffff7e8e370, kw=kw@entry=0x0) at Python/ceval.c:4219
#4  0x00007ffff7b01087 in PyEval_CallFunction (
    obj=obj@entry=0x7ffff5a7cc40 <boost::python::class_metatype_object>,
    format=format@entry=0x7ffff586f02c "(OOO)") at Python/modsupport.c:557
#5  0x00007ffff5857fe5 in boost::python::call<boost::python::api::object, char const*, boost::python::handle<_object>, boost::python::dict> (a2=…, a1=<synthetic pointer>, a0=<synthetic pointer>,
    callable=0x7ffff5a7cc40 <boost::python::class_metatype_object>) at ./boost/python/call.hpp:66
#6  boost::python::api::object_operators<boost::python::api::object>::operator()<char const*, boost::python::handle<_object>, boost::python::dict> (a2=…, a1=<synthetic pointer>,
    a0=<synthetic pointer>, this=<optimized out>) at ./boost/python/object_call.hpp:19
#7  boost::python::objects::(anonymous namespace)::new_class (doc=<optimized out>,
    types=0x7fffffffdd90, num_types=<optimized out>, name=<optimized out>)
    at libs/python/src/object/class.cpp:575
(More stack frames follow…)

ここまで来れば、答えは見えました。この関数の中で PyTypeObject を作って返すのですが、tp_flags をセットするときに Py_TPFLAGS_DEFAULT を使っています。したがって、新しい型は必ず Py_TPFLAGS_HAVE_ITER フラグを持っています。というかわざわざデバッグしなくても、Python のドキュメントのどこかに書いてありそう。

2066 static PyObject *
2067 type_new(PyTypeObject *metatype, PyObject *args, PyObject *kwds)
2068 {
2069     PyObject *name, *bases, *dict;
2070     static char *kwlist[] = {"name", "bases", "dict", 0};

2326     /* Initialize tp_flags */
2327     type->tp_flags = Py_TPFLAGS_DEFAULT | Py_TPFLAGS_HEAPTYPE |
2328         Py_TPFLAGS_BASETYPE;

2329     if (base->tp_flags & Py_TPFLAGS_HAVE_GC)
2330         type->tp_flags |= Py_TPFLAGS_HAVE_GC;
2331     if (base->tp_flags & Py_TPFLAGS_HAVE_NEWBUFFER)
2332         type->tp_flags |= Py_TPFLAGS_HAVE_NEWBUFFER;

4. PyTypeObject::tp_flags を動的に変更する方法

Boost.Python や Python の実装を見ても、tp_flags を変更する API は用意されていないようです。しかし boost::python::objects::class_base::class_base の実装から分かるように、各型のひな形となるような TypeObject は、 boost::python::converter::registry::lookup で取ってきたコンバーターの m_class_object から取得できそうです。そこで、以下のようなコードを試してみます。

#include <boost/python.hpp>
#include <valarray>
#include <random>
#include <stdio.h>

std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(0.0, 1.0);

void ChangeType();

typedef std::valarray<double> darray;

class Field {
private:
    darray _data;

public:
    Field(int n) : _data(n) {
        printf("Field::Field(): %p\n", this);
        ChangeType();
    }
    ~Field() {
        printf("Field::~Field() %p\n", this);
    }

    void Dump() const {
        printf("[");
        for (auto &d : _data) {
            printf("% 3.3f", d);
        }
        printf(" ]\n");
    }

    void Churn() {
        for (auto &d : _data) {
            d = distribution(generator);
        }
    }

    const darray &Positions() const {
        return _data;
    }
};

using namespace boost::python;

BOOST_PYTHON_MODULE(ems) {
    class_<Field>("Field", init<int>())
        .def("Positions", &Field::Positions,
             return_internal_reference<>())
        .def("Dump", &Field::Dump)
        .def("Churn", &Field::Churn)
        ;

    class_<darray>("darray")
        .def("__getitem__",
             (const double &(darray::*)(size_t) const)&darray::operator[],
             return_value_policy<copy_const_reference>())
        .def("__len__", &darray::size)
        ;
}

void ChangeType() {
    converter::registration const& converters =
        converter::registry::lookup(type_id<darray>());
    long l = converters.m_class_object->tp_flags;
    converters.m_class_object->tp_flags = l & ~Py_TPFLAGS_HAVE_ITER;
    printf("converters.m_class_object = %p\n", converters.m_class_object);
    printf("tp_flags: %08lx -> %08lx\n", l, converters.m_class_object->tp_flags);
}

ChangeType() という関数を作って、Field のコンストラクターから呼ぶコードを追加しました。これを Python から実行して、darray の tp_flags を確認します。

john@ubuntu1510:~/Documents/pyc$ gdb /usr/local/python/current/bin/python
GNU gdb (Ubuntu 7.10-1ubuntu2) 7.10
(略)
Reading symbols from /usr/local/python/current/bin/python…done.
(gdb) r
Starting program: /usr/local/python/python-2.7.11/bin/python
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
Python 2.7.11 (default, Dec 22 2015, 22:30:16)
[GCC 5.2.1 20151010] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import ems
>>> f = ems.Field(4)
Field::Field(): 0x6d09d0
converters.m_class_object = 0x6d3b00
tp_flags: 000257fb -> 0002577b

>>> x = f.Positions()
>>> x
<ems.darray object at 0x7ffff7ed5600>
>>>
Program received signal SIGINT, Interrupt.
0x00007ffff74e4723 in __select_nocancel () at ../sysdeps/unix/syscall-template.S:81
81      ../sysdeps/unix/syscall-template.S: No such file or directory.
(gdb) p *(PyObject*)0x7ffff7ed5600
$1 = {ob_refcnt = 2, ob_type = 0x6d3b00}
(gdb) p *$1.ob_type
$2 = {ob_refcnt = 5, ob_type = 0x7ffff5a7cc40 <boost::python::class_metatype_object>, ob_size = 0,
  tp_name = 0x7ffff7e94954 "darray", tp_basicsize = 48, tp_itemsize = 1,
  tp_dealloc = 0x7ffff7a8b390 <subtype_dealloc>, tp_print = 0x0, tp_getattr = 0x0,
  tp_setattr = 0x0, tp_compare = 0x0, tp_repr = 0x7ffff7a8f390 <object_repr>,
  tp_as_number = 0x6d3c88, tp_as_sequence = 0x6d3dd8, tp_as_mapping = 0x6d3dc0,
  tp_hash = 0x7ffff7a70be0 <_Py_HashPointer>, tp_call = 0x0, tp_str = 0x7ffff7a8acc0 <object_str>,
  tp_getattro = 0x7ffff7a72410 <PyObject_GenericGetAttr>,
  tp_setattro = 0x7ffff7a72670 <PyObject_GenericSetAttr>, tp_as_buffer = 0x6d3e28,
  tp_flags = 153467, tp_doc = 0x0, tp_traverse = 0x7ffff7a8b0a0 <subtype_traverse>,
  tp_clear = 0x7ffff7a8d290 <subtype_clear>, tp_richcompare = 0x0, tp_weaklistoffset = 32,
  tp_iter = 0x0, tp_iternext = 0x7ffff7a72100 <_PyObject_NextNotImplemented>, tp_methods = 0x0,
  tp_members = 0x6d3e68, tp_getset = 0x0,
  tp_base = 0x7ffff5a7cfa0 <boost::python::objects::class_type_object>, tp_dict = 0x7ffff7e98050,
  tp_descr_get = 0x0, tp_descr_set = 0x0, tp_dictoffset = 24,
  tp_init = 0x7ffff7a91f70 <slot_tp_init>, tp_alloc = 0x7ffff7a8ade0 <PyType_GenericAlloc>,
  tp_new = 0x7ffff5856890
     <boost::python::objects::instance_new(PyTypeObject*, PyObject*, PyObject*)>,
  tp_free = 0x7ffff7b1e2b0 <PyObject_GC_Del>, tp_is_gc = 0x0, tp_bases = 0x7ffff7e96410,
  tp_mro = 0x7ffff7e8e780, tp_cache = 0x0, tp_subclasses = 0x0, tp_weaklist = 0x7ffff7e77f70,
  tp_del = 0x0, tp_version_tag = 0}
(gdb) p/x 153467
$3 = 0x2577b

思惑通り、tp_flags から Py_TPFLAGS_HAVE_ITER を消せました。今回は tp_flags を変えても全く意味がなかったのですが、これで TypeObject の内容は思いのままです。Boost.Python では、初期化時に int など基本型の TypeObject も列挙しているので、基本型の中身を変更する禁じ手も可能です。そんなことをしたらいろいろ壊れそうですが。

広告

How to expose a valarray with Boost.Python

Python のパフォーマンスの不利を補うため、C++ で書いた共有ライブラリを Python のモジュールとして利用することができます。前回の速度比較の結果から、ベクトルの四則演算を std::valarray の形のまま実行すべく C++ のコードを書いたところ、思わぬ結果として Python の罠を見つけたので経緯を含めて紹介します。

1. Boost.Python を使う

まず、C++ で Python モジュールを書く方法ですが、伝統的な Python.h をインクルードして C で書く方法と、Boost.Python のテンプレートを使って C++ で書く方法の 2 通りがあります。Boost の方が圧倒的に楽です。以下のサイトの説明がとても分かりやすい。

c/c++をラップしてpythonで使えるように – Python | Welcome to underground
https://www.quark.kj.yamagata-u.ac.jp/~hiroki/python/?id=19

Boost.Python の機能をざっと紹介してみる – muddy brown thang
http://d.hatena.ne.jp/moriyoshi/20091214/1260779899

その他、ctypes モジュールを使って外部モジュールのエクスポート関数を呼び出す方法 (Windows で言うところの LoadLibrary と GetProcAddress を Python から呼べる) もありますが、これは関数が呼べるだけで、Python のモジュールやオブジェクトを作れるわけではないので、少し趣向が異なります。もし ctypes で間に合うなら一番楽な方法だと思います。

今実現したいのは、次のような動作を実装した C++ クラスを Python モジュールとして見せることです。

  1. std::valarray<double> を内部データとして保持
  2. 演算は valarray のまま実行する
  3. プロット描画のため、valarray を何らかの形で matplotlib.collections.PathCollection.set_offsets(offsets) に渡したい

3. についてですが、matplotlib のリファレンスを見ると、set_offsets の引数は "offsets can be a scalar or a sequence." という妙に曖昧な説明しかありません。

collections — Matplotlib 1.5.0 documentation
http://matplotlib.org/api/collections_api.html#matplotlib.collections.PathCollection

GitHub 上のコードを見ると、引数の offsets はすぐに np.asanyarray によって np.ndarray に変換されて 2 列の行列に変形させられるので、ndarray に変換可能なオブジェクトなら何でもよい、ということになります。

https://github.com/matplotlib/matplotlib/blob/master/lib/matplotlib/collections.py

C++ 側で列数が 2 の numpy.ndarray オブジェクトを作ることができれば速度的にはベストですが、面倒そうなのと、C++ のモジュールが Numpy に依存するのもあまり美しくないので、Python のリストを返すことができれば十分です。そこで、上記のサイトを参考にしつつ動作確認のため以下の C++ コードを書きました。

実行環境は前回と同じです。

  • OS: Ubuntu 15.10
  • gcc (Ubuntu 5.2.1-22ubuntu2) 5.2.1 20151010
  • Python 2.7.10

C++ ソースコード: shared.cpp

#include <boost/python.hpp>
#include <valarray>
#include <random>
#include <stdio.h>

std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(0.0, 1.0);

typedef std::valarray<double> darray;

class Field {
private:
    darray _data;

public:
    Field(int n) : _data(n) {
        printf("Field::Field(): %p\n", this);
    }
    ~Field() {
        printf("Field::~Field() %p\n", this);
    }

    void Dump() const {
        printf("[");
        for (auto &d : _data) {
            printf("% 3.3f", d);
        }
        printf(" ]\n");
    }

    void Churn() {
        for (auto &d : _data) {
            d = distribution(generator);
        }
    }

    const darray &Positions() const {
        return _data;
    }
};

BOOST_PYTHON_MODULE(ems) {
    using namespace boost::python;

    class_<Field>("Field", init<int>())
        .def("Positions", &Field::Positions,
             return_value_policy<copy_const_reference>())
        .def("Dump", &Field::Dump)
        .def("Churn", &Field::Churn)
        ;

    class_<darray>("darray")
        .def("__getitem__",
             (const double &(darray::*)(size_t) const)&darray::operator[],
             return_value_policy<copy_const_reference>())
        .def("__len__", &darray::size)
        ;
}

メイク ファイル: Makefile

CC=g++
RM=rm -f

TARGET=ems.so
SRCS=$(wildcard *.cpp)
OBJS=$(SRCS:.cpp=.o)

override CFLAGS+=-Wall -fPIC -std=c++11 -O2 -g
LFLAGS=

INCLUDES=-I/usr/include/python2.7
LIBS=-lpython2.7 -lboost_python

all: clean $(TARGET)

clean:
        $(RM) $(OBJS) $(TARGET)

$(TARGET): $(OBJS)
        $(CC) -shared $(LFLAGS) $(LIBDIRS) $^ -o $@ $(LIBS)

$(OBJS): $(SRCS)
        $(CC) $(INCLUDES) $(CFLAGS) -c $^

Field クラスを Python オブジェクトとして扱えるようにします。このとき、valarray の const 参照を返すメンバ関数の Field::Positions を Python オブジェクトのメソッドとして実装するため、valarray のコンバーターを別途作って __getitem__ と __len__ メソッドを実装します。これは一種のダック タイピングで、Python 側ではシーケンスとして見えるようになります。

Field::Positions は Field::_data の参照を返すようにしておかないと、Positions を呼ぶたびに valarray をコピーして Python に渡すようになるので、速度が落ちます。しかし、実は上記のコードは間違っています。試しにこのまま Python から使ってみると以下の結果が得られます。

>>> import ems
>>> f = ems.Field(4)
Field::Field(): 0x2011a50
>>> x = f.Positions()
>>> f.Dump()
[ 0.000 0.000 0.000 0.000 ]
>>> x[0]
0.0

>>> f.Churn()
>>> f.Dump()
[ 0.132 0.459 0.219 0.679 ]
>>> x[0]
0.0

>>> quit()
Field::~Field() 0x2011a50
john@ubuntu1510:~/Documents/pyc$

f.Churn() で C++ 側の値を変更した後も、Python 側のシーケンスの値が 0 のままで変わっていません。x = f.Position() したときに、valarray がコピーされたようです。Positions を定義するときの return_value_policy で copy_const_reference ポリシーを使っているのがいけないのです。

    class_<Field>("Field", init<int>())
        .def("Positions", &Field::Positions,
             return_value_policy<copy_const_reference>())

以下のページに記載があるように、このポリシーは "returning a reference-to-const type such that the referenced value is copied into a new Python object." です。ポリシーの名前で勘違いしていましたが、これは、参照をコピーして Python に渡すのではなく、"参照された値" をコピーします。Positions の戻り値となる Python オブジェクトを作るときに、参照に対してコピー コンストラクタを呼んで、そのコピーをオブジェクトに保持させるのだと思います。これでは参照を使う意味がありません。

Boost.Python – <boost/python/copy_const_reference.hpp> – 1.54.0
http://www.boost.org/doc/libs/1_54_0/libs/python/doc/v2/copy_const_reference.html

ここで使うべきポリシーは、return_internal_reference です。説明によると "CallPolicies which allow pointers and references to objects held internally by a free or member function argument or from the target of a member function to be returned safely without making a copy of the referent" となっています。

Boost.Python – <boost/python/return_internal_reference.hpp> – 1.54.0
http://www.boost.org/doc/libs/1_54_0/libs/python/doc/v2/return_internal_reference.html

コードをこう変えます。

class_<Field>("Field", init<int>())
    .def("Positions", &Field::Positions,
         return_internal_reference<>())

これで、C++ オブジェクトへの参照を Python が保持できるようになりました。

>>> import ems
>>> f = ems.Field(4)
Field::Field(): 0xde4a50
>>> f.Dump()
[ 0.000 0.000 0.000 0.000 ]
>>> x = f.Positions()
>>> x[0]
0.0
>>> f.Churn()
>>> f.Dump()
[ 0.132 0.459 0.219 0.679 ]
>>> x[0]
0.13153778773876065

>>> x
<ems.darray object at 0x7fded8919130>
>>> y = f.Positions()
>>> y
<ems.darray object at 0x7fded89193d0>
>>> y[0]
0.13153778773876065
>>> quit()
Field::~Field() 0xde4a50

上記の例で、f.Positions() の戻り値である ems.darray オブジェクトの x と y を単純に表示させた値が異なっています。これも最初勘違いしていましたが、Positions() を呼ぶたびに、ems.darray という Python オブジェクトは常に新しく作られます。それぞれの ems.darray が同じ valarray への参照を内包する構造になっています。

2. Python のハングと原因

ここから本題です。先ほど作った ems.Field を matplotlib の set_offsets に渡すとハングすることが分かりました。上で既に触れたとおり、set_offsets は引数を ndarray に変換しているだけなので、以下のコードで簡単に再現できます。

>>> import ems
>>> import numpy as np
>>> np.__version__
‘1.10.1’
>>> f = ems.Field(4)
Field::Field(): 0x16f90b0
>>> x = f.Positinos()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: ‘Field’ object has no attribute ‘Positinos’
>>> x = f.Positions()
>>> x[0]
0.0
>>> len(x)
4
>>> a = np.array(x) # ここでハング

__getitem__ と __len__ だけだとリストに見せかけるのは無理なのかと諦めかけていたのですが、偶然、Numpy 1.8 だと問題なく動くことが分かりました。

>>> import ems
>>> import numpy as np
>>> np.__version__
‘1.8.2’
>>> f = ems.Field(4)
Field::Field(): 0x233acf0
>>> x = f.Positions()
>>> a = np.array(x)
>>> x[0]
0.0
>>> a
array([ 0.,  0.,  0.,  0.])

Numpy 側のコードを比べると、1.8.x から 1.9.x にバージョンが上がった時に、PyArray_DTypeFromObjectHelper() のロジックが少し変更されています。520 行目あたりで青字の部分が追加されており、ハングは PySequence_Fast から制御が返ってこないことで発生しています。

https://github.com/numpy/numpy/blob/maintenance/1.8.x/numpy/core/src/multiarray/common.c
https://github.com/numpy/numpy/blob/maintenance/1.9.x/numpy/core/src/multiarray/common.c

/*
* fails if convertable to list but no len is defined which some libraries
* require to get object arrays
*/
size = PySequence_Size(obj);
if (size < 0) {
    goto fail;
}

/* Recursive case, first check the sequence contains only one type */
seq = PySequence_Fast(obj, "Could not convert object to sequence");
if (seq == NULL) {
    goto fail;
}

objects = PySequence_Fast_ITEMS(seq);
common_type = size > 0 ? Py_TYPE(objects[0]) : NULL;
for (i = 1; i < size; ++i) {
    if (Py_TYPE(objects[i]) != common_type) {
        common_type = NULL;
        break;
    }
}

結論から先に書くと、ems.darray の __getitem__ を実装するときに、valarray::operator[] を渡しているのが諸悪の根源でした。

    class_<darray>("darray")
        .def("__getitem__",
             (const double &(darray::*)(size_t) const)&darray::operator[],
             return_value_policy<copy_const_reference>())

ネット上でよく出てくるのは、valarray ではなく vector を Python に渡す方法ですが、この場合の __getitem__ は、operator[] ではなく vector::at を使うものばかりです。恥ずかしながら知らなかったのですが、vector::at と vector::operator[] には大きな違いがありました。vector::at は配列の長さを超えた場合に out_of_range 例外を投げてくれますが、operator[] は何もチェックせずに配列のインデックス アクセスを淡々とこなします。

vector::at – C++ Reference
http://www.cplusplus.com/reference/vector/vector/at/

The function automatically checks whether n is within the bounds of valid elements in the vector, throwing an out_of_range exception if it is not (i.e., if n is greater or equal than its size). This is in contrast with member operator[], that does not check against bounds.

Python 側に作られた ems.Field は iteratable なオブジェクトとして作られ、Python の iterator オブジェクト (seqiterobject 構造体) によって iterate されます。この動作は Objects/iterobject.c にある iter_iternext という関数で行われますが、iterate の終了条件は、戻り値が NULL を返すかどうかです。オブジェクトのサイズはチェックしません。

45 static PyObject *
46 iter_iternext(PyObject *iterator)
47 {
48     seqiterobject *it;
49     PyObject *seq;
50     PyObject *result;
51
52     assert(PySeqIter_Check(iterator));
53     it = (seqiterobject *)iterator;
54     seq = it->it_seq;
55     if (seq == NULL)
56         return NULL;
57     if (it->it_index == LONG_MAX) {
58         PyErr_SetString(PyExc_OverflowError,
59                         "iter index too large");
60         return NULL;
61     }
62
63     result = PySequence_GetItem(seq, it->it_index); // getitem を呼ぶ箇所
64     if (result != NULL) {
65         it->it_index++;
66         return result;
67     }
68     if (PyErr_ExceptionMatches(PyExc_IndexError) ||
69         PyErr_ExceptionMatches(PyExc_StopIteration))
70     {
71         PyErr_Clear();
72         Py_DECREF(seq);
73         it->it_seq = NULL;
74     }
75     return NULL;
76 }

Numpy 1.9.x 以上でハングが起こる原因は、Objects/listobject.c にある listextend で無限ループが発生するからです。iter_iternext がどんなインデックスに対しても値を返してくるので、このループが終わりません。

870     /* Run iterator to exhaustion. */
871     for (;;) {
872         PyObject *item = iternext(it);
873         if (item == NULL) {
874             if (PyErr_Occurred()) {
875                 if (PyErr_ExceptionMatches(PyExc_StopIteration))
876                     PyErr_Clear();
877                 else
878                     goto error;
879             }
880             break;
881         }
882         if (Py_SIZE(self) < self->allocated) {
883             /* steals ref */
884             PyList_SET_ITEM(self, Py_SIZE(self), item);
885             ++Py_SIZE(self);
886         }
887         else {
888             int status = app1(self, item);
889             Py_DECREF(item);  /* append creates a new ref */
890             if (status < 0)
891                 goto error;
892         }
893     }

他の人の役に立つか不明ですが、一応コールスタックなどの情報を含むデバッグ ログを載せておきます。

(gdb) r
Starting program: /usr/local/python/python-2.7.11/bin/python
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
Python 2.7.11 (default, Dec 21 2015, 12:40:02)
[GCC 5.2.1 20151010] on linux2
Type "help", "copyright", "credits" or "license" for more information.

>>> import numpy as np
>>> import ems
>>> f = ems.Field(1)
Field::Field(): 0xaa0060
>>> x = f.Positions()
>>> f.Churn()
>>>

Program received signal SIGINT, Interrupt.
0x00007ffff71df723 in __select_nocancel () at ../sysdeps/unix/syscall-template.S:81
81      ../sysdeps/unix/syscall-template.S: No such file or directory.
(gdb) b slot_sq_item
Breakpoint 1 at 0x4816d0: file Objects/typeobject.c, line 4987.
(gdb) c
Continuing.

>>> a = np.array(x)

Breakpoint 1, slot_sq_item (self=0x7ffff1a92670, i=0) at Objects/typeobject.c:4987
4987    {
(gdb) p i
$1 = 0
(gdb) c
Continuing.
Breakpoint 1, slot_sq_item (self=0x7ffff1a92670, i=1) at Objects/typeobject.c:4987
4987    {
(gdb) p i
$2 = 1 <<<< この時点で既に out-of-range
(gdb) c
Continuing.
Breakpoint 1, slot_sq_item (self=0x7ffff1a92670, i=2) at Objects/typeobject.c:4987
4987    {
(gdb) p i
$3 = 2

(gdb) bt
#0  slot_sq_item (self=0x7ffff1a92670, i=2) at Objects/typeobject.c:4987
#1  0x0000000000444e5d in iter_iternext (iterator=0x7ffff7e96910) at Objects/iterobject.c:63
#2  0x0000000000448e9a in listextend (self=self@entry=0x7ffff1ab4cb0, b=b@entry=0x7ffff7e96910)
    at Objects/listobject.c:872
#3  0x000000000044aeb5 in _PyList_Extend (self=self@entry=0x7ffff1ab4cb0, b=b@entry=0x7ffff7e96910)
    at Objects/listobject.c:910
#4  0x000000000042243c in PySequence_List (v=0x7ffff7e96910) at Objects/abstract.c:2264
#5  PySequence_Fast (v=v@entry=0x7ffff1a92670, m=m@entry=0x7ffff5c299e8 "Could not convert object to sequence")
    at Objects/abstract.c:2293
#6  0x00007ffff5b5473c in PyArray_DTypeFromObjectHelper (obj=obj@entry=0x7ffff1a92670, maxdims=maxdims@entry=32,
    out_dtype=out_dtype@entry=0x7fffffffdea8, string_type=string_type@entry=0)
    at numpy/core/src/multiarray/common.c:531
#7  0x00007ffff5b549c3 in PyArray_DTypeFromObject (obj=obj@entry=0x7ffff1a92670, maxdims=maxdims@entry=32,
    out_dtype=out_dtype@entry=0x7fffffffdea8) at numpy/core/src/multiarray/common.c:184
#8  0x00007ffff5b5de81 in PyArray_GetArrayParamsFromObject (op=0x7ffff1a92670, requested_dtype=<optimized out>,
    writeable=<optimized out>, out_dtype=0x7fffffffdea8, out_ndim=0x7fffffffde9c, out_dims=0x7fffffffdeb0,
    out_arr=0x7fffffffdea0, context=0x0) at numpy/core/src/multiarray/ctors.c:1542
#9  0x00007ffff5b5e26d in PyArray_FromAny (op=op@entry=0x7ffff1a92670, newtype=0x0, min_depth=0, max_depth=0,
    flags=flags@entry=112, context=<optimized out>) at numpy/core/src/multiarray/ctors.c:1674
#10 0x00007ffff5b5e63f in PyArray_CheckFromAny (op=0x7ffff1a92670, descr=<optimized out>, min_depth=min_depth@entry=0,
    max_depth=max_depth@entry=0, requires=112, context=context@entry=0x0) at numpy/core/src/multiarray/ctors.c:1852
#11 0x00007ffff5bd850e in _array_fromobject (__NPY_UNUSED_TAGGEDignored=<optimized out>, args=<optimized out>, kws=0x0)
    at numpy/core/src/multiarray/multiarraymodule.c:1773
#12 0x00000000004b6f2a in call_function (oparg=<optimized out>, pp_stack=0x7fffffffe110) at Python/ceval.c:4350
#13 PyEval_EvalFrameEx (f=f@entry=0x7ffff7fc0c20, throwflag=throwflag@entry=0) at Python/ceval.c:2987
#14 0x00000000004b938c in PyEval_EvalCodeEx (co=co@entry=0x7ffff7ecfe30, globals=globals@entry=0x7ffff7f6c168,
    locals=locals@entry=0x7ffff7f6c168, args=args@entry=0x0, argcount=argcount@entry=0, kws=kws@entry=0x0, kwcount=0,
    defs=0x0, defcount=0, closure=0x0) at Python/ceval.c:3582
#15 0x00000000004b9499 in PyEval_EvalCode (co=co@entry=0x7ffff7ecfe30, globals=globals@entry=0x7ffff7f6c168,
    locals=locals@entry=0x7ffff7f6c168) at Python/ceval.c:669
#16 0x00000000004e2b4d in run_mod (arena=0x856e10, flags=<optimized out>, locals=0x7ffff7f6c168,
    globals=0x7ffff7f6c168, filename=0x546f03 "<stdin>", mod=0x8a1c48) at Python/pythonrun.c:1370
#17 PyRun_InteractiveOneFlags (fp=0x7ffff74a6980 <_IO_2_1_stdin_>, filename=0x546f03 "<stdin>", flags=<optimized out>)
    at Python/pythonrun.c:857
#18 0x00000000004e2dfe in PyRun_InteractiveLoopFlags (fp=fp@entry=0x7ffff74a6980 <_IO_2_1_stdin_>,
    filename=filename@entry=0x546f03 "<stdin>", flags=flags@entry=0x7fffffffe3a0) at Python/pythonrun.c:777
#19 0x00000000004e3366 in PyRun_AnyFileExFlags (fp=0x7ffff74a6980 <_IO_2_1_stdin_>, filename=<optimized out>,
    closeit=0, flags=0x7fffffffe3a0) at Python/pythonrun.c:746
#20 0x0000000000416160 in Py_Main (argc=<optimized out>, argv=<optimized out>) at Modules/main.c:640
#21 0x00007ffff7102a40 in __libc_start_main (main=0x415290 <main>, argc=1, argv=0x7fffffffe568, init=<optimized out>,
    fini=<optimized out>, rtld_fini=<optimized out>, stack_end=0x7fffffffe558) at libc-start.c:289
#22 0x00000000004152c9 in _start ()

4. 回避策

valarray には at 関数がなく、境界チェックをしない operator[] があるのみです。したがって回避策として、valarray をラップする darray クラスを C++ 側で作って、at を追加実装しました。記憶が曖昧ですが、boost を使えば .NET のようにクラスを拡張できたような気がしますが、今回はこれで。

#include <boost/python.hpp>
#include <valarray>
#include <random>
#include <stdio.h>

std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(0.0, 1.0);

class darray : public std::valarray<double> {
public:
    darray() : std::valarray<double>() {}
    darray(size_t n) : std::valarray<double>(n) {}

    const double& at(size_t n) const {
        if (n >= this->size()) {
            std::__throw_out_of_range_fmt(__N("out_of_range"));
        }
        return (*this)[n];
    }
};

class Field {
private:
    darray _data;

public:
    Field(int n) : _data(n) {
        printf("Field::Field(): %p\n", this);
    }
    ~Field() {
        printf("Field::~Field() %p\n", this);
    }

    void Dump() const {
        printf("[");
        for (auto &d : _data) {
            printf("% 3.3f", d);
        }
        printf(" ]\n");
    }

    void Churn() {
        for (auto &d : _data) {
            d = distribution(generator);
        }
    }

    const darray &Positions() const {
        return _data;
    }
};

BOOST_PYTHON_MODULE(ems) {
    using namespace boost::python;

    class_<Field>("Field", init<int>())
        .def("Positions", &Field::Positions,
             return_internal_reference<>())
        .def("Dump", &Field::Dump)
        .def("Churn", &Field::Churn)
        ;

    class_<darray>("darray")
        .def("__getitem__", &darray::at,
             return_value_policy<copy_const_reference>())
        .def("__len__", &darray::size)
        ;
}

Numpy 1.10.1 でも問題なく動作するようになりました。

>>> import ems
>>> import numpy as np
>>> np.__version__
‘1.10.1’
>>> f = ems.Field(4)
Field::Field(): 0x1029880
>>> x = f.Positions()
>>> f.Churn()
>>> f.Dump()
[ 0.132 0.459 0.219 0.679 ]
>>> np.array(x)
array([ 0.13153779,  0.45865013,  0.21895919,  0.67886472])
>>> quit()
Field::~Field() 0x1029880

この場合の動作ですが、darray::at で投げた例外が boost_python ライブラリの中の例外ハンドラーによって捕捉されるようです。たぶん #1 か #2 が例外ハンドラー。

#0  PyErr_SetString (exception=0x7a2ce0 <_PyExc_IndexError>, string=0x7ea428 "out_of_range") at Python/errors.c:68
#1  0x00007ffff56f3471 in boost::python::handle_exception_impl(boost::function0<void>) ()
   from /usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.58.0
#2  0x00007ffff56e8719 in ?? () from /usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.58.0
#3  0x0000000000422faa in PyObject_Call (func=func@entry=0x8bce90, arg=arg@entry=0x7ffff0b4f710, kw=kw@entry=0x0)
    at Objects/abstract.c:2546
#4  0x0000000000429dfc in instancemethod_call (func=0x8bce90, arg=0x7ffff0b4f710, kw=0x0) at Objects/classobject.c:2602
#5  0x0000000000422faa in PyObject_Call (func=func@entry=0x7ffff0b3ce60, arg=arg@entry=0x7ffff7e965d0, kw=kw@entry=0x0)
    at Objects/abstract.c:2546
#6  0x0000000000481766 in slot_sq_item (self=<optimized out>, i=<optimized out>) at Objects/typeobject.c:5012
#7  0x0000000000444e5d in iter_iternext (iterator=0x7ffff0b98e50) at Objects/iterobject.c:63
#8  0x0000000000448e9a in listextend (self=self@entry=0x7ffff0b4f5f0, b=b@entry=0x7ffff0b98e50)
    at Objects/listobject.c:872
#9  0x000000000044aeb5 in _PyList_Extend (self=self@entry=0x7ffff0b4f5f0, b=b@entry=0x7ffff0b98e50)
    at Objects/listobject.c:910
#10 0x000000000042243c in PySequence_List (v=0x7ffff0b98e50) at Objects/abstract.c:2264
#11 PySequence_Fast (v=v@entry=0x7ffff0b24670, m=m@entry=0x7ffff4aa09e8 "Could not convert object to sequence")
    at Objects/abstract.c:2293
#12 0x00007ffff49cb73c in PyArray_DTypeFromObjectHelper (obj=obj@entry=0x7ffff0b24670, maxdims=maxdims@entry=32,
    out_dtype=out_dtype@entry=0x7fffffffdea8, string_type=string_type@entry=0)
    at numpy/core/src/multiarray/common.c:531

このハングは、numpy.ndarray に限りません。例えばリスト内包表記で ems::darray を列挙しようとしても再現します。そのときのコールスタックはこのようになります。

#0  darray::at (this=0xb60d30, n=1) at shared.cpp:16
#1  0x00007ffff5e933bb in boost::python::detail::invoke<boost::python::to_python_value<double const&>, double const& (darray::*)(unsigned long) const, boost::python::arg_from_python<darray&>, boost::python::arg_from_python<unsigned long> >
    (rc=…, ac0=…, tc=<synthetic pointer>, f=
    @0x8bcc58: (const double &(darray::*)(const darray * const, unsigned long)) 0x7ffff5e92170 <darray::at(unsigned long) const>) at /usr/include/boost/python/detail/invoke.hpp:88
#2  boost::python::detail::caller_arity<2u>::impl<double const& (darray::*)(unsigned long) const, boost::python::return_value_policy<boost::python::copy_const_reference, boost::python::default_call_policies>, boost::mpl::vector3<double const&, darray&, unsigned long> >::operator() (args_=<optimized out>, this=0x8bcc58)
    at /usr/include/boost/python/detail/caller.hpp:223
#3  boost::python::objects::caller_py_function_impl<boost::python::detail::caller<double const& (darray::*)(unsigned long) const, boost::python::return_value_policy<boost::python::copy_const_reference, boost::python::default_call_policies>, boost::mpl::vector3<double const&, darray&, unsigned long> > >::operator() (this=0x8bcc50, args=<optimized out>,
    kw=<optimized out>) at /usr/include/boost/python/object/py_function.hpp:38
#4  0x00007ffff56eb33d in boost::python::objects::function::call(_object*, _object*) const ()
   from /usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.58.0
#5  0x00007ffff56eb528 in ?? () from /usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.58.0
#6  0x00007ffff56f3363 in boost::python::handle_exception_impl(boost::function0<void>) ()
   from /usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.58.0
#7  0x00007ffff56e8719 in ?? () from /usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.58.0
#8  0x0000000000422faa in PyObject_Call (func=func@entry=0x8bce90, arg=arg@entry=0x7ffff0b4f758, kw=kw@entry=0x0)
    at Objects/abstract.c:2546
#9  0x0000000000429dfc in instancemethod_call (func=0x8bce90, arg=0x7ffff0b4f758, kw=0x0) at Objects/classobject.c:2602
#10 0x0000000000422faa in PyObject_Call (func=func@entry=0x7ffff0b3ce60, arg=arg@entry=0x7ffff7e964d0, kw=kw@entry=0x0)
    at Objects/abstract.c:2546
#11 0x0000000000481766 in slot_sq_item (self=<optimized out>, i=<optimized out>) at Objects/typeobject.c:5012
#12 0x0000000000444e5d in iter_iternext (iterator=0x7ffff0b98e50) at Objects/iterobject.c:63
#13 0x00000000004b2373 in PyEval_EvalFrameEx (f=f@entry=0x7ffff1b3faa0, throwflag=throwflag@entry=0)
    at Python/ceval.c:2806

フレーム#12 と #13 から分かるように、今度は PyEval_EvalFrameEx から iter_iternext が呼ばれています。仮にこのハングを Python 側で直すとして、iter_iternext の中でサイズをチェックする、というアプローチが考えられます。Numpy のコードから確認できますが、PySequence_Size は darray::__len__ を呼ぶので、サイズは正しく取得できます。が、コードを見る限り、Python の iterator オブジェクトはサイズをチェックしなくても iterator するだけで要素を列挙できるのが利点なので、iter_iternext はこのままのほうが自然です。

もう少し見てみないと分かりませんが、ems.Field オブジェクトが iterator を使ってシーケンシャルにアクセスする種類になっているのが悪いのではないかと思います。何らかの方法で、ems.Field は iterator ではなくインデックスを使ってランダム アクセスすべきものに種別できればちゃんと動くのではないかと。そんな種類があるかどうか分かりませんが。後で調べよう。

4. おまけ – Python のデバッグ環境

Python をデバッグするため、Python そのものをソースからビルドしたので、そのときのコマンドを紹介します。

$ wget https://www.python.org/ftp/python/2.7.11/Python-2.7.11.tgz
$ tar -xvf Python-2.7.11.tgz
$ cd Python-2.7.11/
$ ./configure –prefix=/usr/local/python/python-2.7.11 –enable-unicode=ucs4
$ make
$ sudo make install
$ sudo ln -s /usr/local/python/python-2.7.11 /usr/local/python/current

configure するときに ucs4 を指定しておかないと、numpy をインポートしたときに以下のエラーが出ます。

>>> import numpy as np
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/numpy/__init__.py", line 153, in <module>
    from . import add_newdocs
  File "/usr/lib/python2.7/dist-packages/numpy/add_newdocs.py", line 13, in <module>
    from numpy.lib import add_newdoc
  File "/usr/lib/python2.7/dist-packages/numpy/lib/__init__.py", line 8, in <module>
    from .type_check import *
  File "/usr/lib/python2.7/dist-packages/numpy/lib/type_check.py", line 11, in <module>
    import numpy.core.numeric as _nx
  File "/usr/lib/python2.7/dist-packages/numpy/core/__init__.py", line 6, in <module>
    from . import multiarray
ImportError: /usr/lib/python2.7/dist-packages/numpy/core/multiarray.so: undefined symbol: PyUnicodeUCS4_AsUnicodeEscapeString

このままだと、ビルドした Python のモジュール検索パスが空っぽなので、デバッグするときは環境変数を使って、元からあった Python の dist-packages を追加しておきました。たぶんもっとうまいやり方があるはず。

john@ubuntu1510:~/Documents/pyc$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages
john@ubuntu1510:~/Documents/pyc$ gdb /usr/local/python/current/bin/python
GNU gdb (Ubuntu 7.10-1ubuntu2) 7.10
Copyright (C) 2015 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
Type "show configuration" for configuration details.
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>.
Find the GDB manual and other documentation resources online at:
<http://www.gnu.org/software/gdb/documentation/>.
For help, type "help".
Type "apropos word" to search for commands related to "word"…
Reading symbols from /usr/local/python/current/bin/python…done.
(gdb) r
Starting program: /usr/local/python/python-2.7.11/bin/python
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
Python 2.7.11 (default, Dec 21 2015, 12:40:02)
[GCC 5.2.1 20151010] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>>

Numpy は、以前 pip でインストールしたものにデバッグ情報が入っていたので、わざわざソースからビルドすることなくそのまま使いました。というか、Numpy のビルドは追加モジュールがいろいろ必要で面倒そうだったので諦めたともいいます。

Vectorized computation: R vs Python vs C++

前回のコードで、どうやら Python が遅いらしいことが分かったので、実際に C++ でモジュールを書き始める前に、どのぐらい遅いのかを調べてみることにしました。まずは軽く情報収集。

以下のページによると、Python は動的型付けなので、型チェックのせいで四則演算が遅くなるらしい。これは納得できる。さらに、NumPy は SIMD に対応していて、Cython を使うと速くなるとか。

Pythonを高速化するCythonを使ってみた – Kesin’s diary
http://kesin.hatenablog.com/entry/20120306/1331043675

その一方、R や Python の方が C より速くなる、というような不思議な結果を報告しているページもありました。この結論は原理的におかしいような気もしますが・・。

機械学習序曲(仮) > 数値計算の実行速度 – Python, C, Rの比較
http://hawaii.naist.jp/~kohei-h/my/python-numeric-comp.html

(2015/12/20 – 議論が雑だったので、データを取り直して全体的に修正。)

今回計測したいのは、2 次元ベクトルの四則演算 (オイラー法で、微小時間におけるベクトルの変化を計算するため)、及び、ばねの力を計算するときに使うノルムの計算です。というわけで、SSE や AVX といったベクトル演算用のユニット (= SIMD) をちゃんと使いたいところです。具体的には、こんな計算を計測することにしました。

V1 = V1 + d * V2

V1 と V2 は 2 次元ベクトルを横に P 個並べた行列、d はスカラーで、(0, 1) の範囲乱数を使います。この簡単な加算と乗算を N 回繰り返して速度を測ります。ノルムの演算は気が向いたら後で書き足します。

実行環境は以下の通り。

  • OS: Ubuntu 15.10
  • CPU: Intel(R) Core(TM) i5-2520M CPU @ 2.50GHz (Sandy Bridge)
  • gcc (Ubuntu 5.2.1-22ubuntu2) 5.2.1 20151010
  • Eigen 3.2.7
  • Python 2.7.10
  • R version 3.2.3 (2015-12-10) — "Wooden Christmas-Tree"

R と Python はそのまま、C++ では 6 つの処理系で比較します。C++ を使った場合の速度の違いで見ておきたいポイントは 2 つあり、一つは SIMD 演算で従来型の SSE を使うか Sandy Bridge から搭載された AVX 命令を使うかどうか、もう一つは、SSE/AVX 演算ユニットのスカラー演算を使うかベクトル演算を使うかという点です。前者は -mavx オプションで制御できます。後者は intrinsic を使って明示的に実装することもできますが、今回は gcc の自動ベクトル化機能を -ftree-vectorize オプションで有効にして利用します。ただし、最適化オプション -O3 を使うと、自動ベクトル化は勝手に行われます。よって、今回はベクトル化演算を行う前提で SSE と AVX でどれだけ違いが出るのかを見るのが趣旨です。というのも、SIMD は本来ベクトル演算向けのユニットなので、スカラー演算で使ってもあまり意味がないからです。

gcc のベクトル化に関しては ↓ が参考になります。

Auto-vectorization in GCC – GNU Project – Free Software Foundation (FSF)
https://gcc.gnu.org/projects/tree-ssa/vectorization.html

まとめるとこんな感じです。

Language Matrix datatype Gcc option SIMD
R num(,) N/A ?
Python numpy.ndarray N/A ?
C++ std::vector -O3 vectorized, SSE
C++ std:valarray -O3 vectorized, SSE
C++ Eigen::MatrixXd -O3 vectorized, SSE
C++ std::vector -O3 -ftree-vectorize -mavx vectorized, AVX
C++ std::valarray -O3 -ftree-vectorize -mavx vectorized, AVX
C++ Eigen::MatrixXd -O3 -ftree-vectorize -mavx vectorized, AVX

で、書いたコードがこれです。いろいろ詰め込んだせいで無駄にファイルが多い。

https://github.com/msmania/ems/tree/perftest

主要なファイルについてはこんな感じ。

  • perftest.R – R の計算ルーチン
  • perftest.py – Python の計算ルーチン
  • perftest.cpp – C++ の計算ルーチン
  • main.cpp – 実行可能ファイル runc 用の main()
  • shared.cpp – Python モジュール化のためのコード
  • common.h / common.cpp – shared.cpp と main.cpp とで共有するコード
  • draw.R – 結果を折れ線グラフで描く R のコード

コードを make すると、Python モジュールとしてインポートできる perftest.so と、スタンドアロンの実行可能ファイルである runc ができます。いずれのファイルにも、AVX を使った場合の関数と、使っていない場合の関数がそれぞれ avx、noavx という名前空間内に作られます。

Makefile の中で、perftest.cpp を別のレシピでコンパイルするようにしました。

override CFLAGS+=-Wall -fPIC -std=c++11 -O3 -g -Wno-deprecated-declarations
VFLAGS=-ftree-vectorize -ftree-vectorizer-verbose=6 -mavx -D_AVX
perftest.o: perftest.cpp
    $(CC) $(INCLUDES) $(CFLAGS) -c $^
perftest_avx.o: perftest.cpp
    $(CC) $(INCLUDES) $(CFLAGS) $(VFLAGS) -c $^ -o $@

-Wno-deprecated-declarations オプションは、Eigen のヘッダーをインクルードしたときに警告が出るので入れています。

では実測タイム。まずは Python から実行。引数の意味は、2 次元ベクトルを 10,000 個並べた行列を 1,000,000 回計算する、という意味です。

john@ubuntu1510:~/Documents/ems$ python perftest.py 10000 1000000
*** Python: numpy.ndarray ***
time = 26.096 (s)
*** R ***
[1] "*** matrix(n, 2) ***"
[1] time = 69.94 (s)
*** C++: no-vectorize no-avx ***
std::vector   time = 11.398 (s)
std::valarray time = 11.371 (s)
Eigen::Matrix time = 10.547 (s)
*** C++: vectorize avx ***
std::vector   time = 10.208 (s)
std::valarray time = 10.221 (s)
Eigen::Matrix time = 10.549 (s)

R 遅すぎ・・・。コンパイルされないからだろうか。Python も C++ と比べると 2 倍以上遅い。C++ については、avx と no-avx を比べると 10% ぐらい avx の方が速そう。もう少し試行回数を増やすため R と Python には退場してもらって、C++ だけ実行します。そんなときに runc を使います。

john@ubuntu1510:~/Documents/ems$ ./runc 10000 1000000 20 > result-O3-1000k
john@ubuntu1510:~/Documents/ems$ R

> source(‘draw.R’)
> drawResult(‘result-O3-1000k’)

runc の stdout への出力結果は R でそのまま読み込めるようにしているので、これをグラフ化します。レジェンドを入れていないので分かりにくいですが、青 = std::vector, 赤 = std::valarray, 紫 = Eigen、また、明るい色が AVX 有、暗い色が AVX 無になっています。

result-O3-1000k

明らかに、AVX 無の vector と valarray は他より遅いです。AVX を使うと約 10% 速度が向上しています。逆に、Eigen は変化なしといってよさそうです。上位 4 つはほとんど差がありませんが、Eigen よりは AVX を使った vector/valarray の方が有利に見えます。

正直、AVX による速度向上が 10% 程度というのは期待外れでした。理想的には、レジスター サイズが xmm 系の 128bit から ymm 系の 256bit になることによって、速度は 2 倍ぐらいになって欲しいのです。本当に AVX が使われているかどうかを確認するため、生成されたアセンブリを見てみます。gdb を使って逆アセンブル結果をファイルに出力します。

john@ubuntu1510:~/Documents/ems$ gdb runc
(gdb) set height 0
(gdb) set logging file perftest.log
(gdb) i func Test
All functions matching regular expression "Test":

File perftest.cpp:
void avx::TestEigen(int, int);
void avx::TestValArray(int, int);
void avx::TestVector(int, int);
void noavx::TestEigen(int, int);
void noavx::TestValArray(int, int);
void noavx::TestVector(int, int);
(gdb) disas avx::TestEigen
(gdb) disas avx::TestValArray
(gdb) disas avx::TestVector
(gdb) disas noavx::TestEigen
(gdb) disas noavx::TestValArray
(gdb) disas noavx::TestVector
(gdb) set logging  off
Done logging to perftest.log.
(gdb) quit

例えばavx::TestVectorをざっと見ると、以下の部分がベクトル演算をしているところです。pd というサフィックスが Packed Double を意味しています。

0x000000000040236b <+1035>:    vmovddup %xmm0,%xmm9
0x000000000040236f <+1039>:    xor    %ecx,%ecx
0x0000000000402371 <+1041>:    xor    %r15d,%r15d
0x0000000000402374 <+1044>:    vinsertf128 $0x1,%xmm9,%ymm9,%ymm9
0x000000000040237a <+1050>:    vmovupd (%r10,%rcx,1),%xmm1
0x0000000000402380 <+1056>:    add    $0x1,%r15d
0x0000000000402384 <+1060>:    vinsertf128 $0x1,0x10(%r10,%rcx,1),%ymm1,%ymm1
0x000000000040238c <+1068>:    vmulpd %ymm9,%ymm1,%ymm1
0x0000000000402391 <+1073>:    vaddpd (%r8,%rcx,1),%ymm1,%ymm1
0x0000000000402397 <+1079>:    vmovapd %ymm1,(%r8,%rcx,1)
0x000000000040239d <+1085>:    add    $0x20,%rcx
0x00000000004023a1 <+1089>:    cmp    %r11d,%r15d
0x00000000004023a4 <+1092>:    jb  0x40237a <avx::TestVector(int, int)+1050>

上記のコードの前半の vmovddup と vinsertf128 で、xmm0 に入っている 128bit の値 (double  2 つ分) を 2 つにコピーして ymm9 にセットします。おそらくこれは乱数の d です。次に、(%r10,%rcx,1) = r10 + (rcx * 1) のアドレスから 32 バイト (16 バイトずつvmovupd とvinsertf128 を使って) コピーして ymm1 にセットします。これは V2 でしょう。ここで AVX 命令を使って、V2 に d をかけて (%r8,%rcx,1) に加算しています。これが V1 になりそうです。このベクトル演算で、4 つの double を同時に処理していることになります。

avx::TestValArray にも同じ処理が見つかります。avx::TestEigen にも packed double を使っているコードは見つかりますが、ymm レジスタを使っていませんでした。これが、AVX の有無で Eigen の速度が変わらない理由として有力です。

0x0000000000402dcf <+511>:    vmovddup %xmm1,%xmm1
0x0000000000402dd3 <+515>:    xor    %eax,%eax
0x0000000000402dd5 <+517>:    nopl   (%rax)
0x0000000000402dd8 <+520>:    vmulpd (%r12,%rax,8),%xmm1,%xmm0
0x0000000000402dde <+526>:    vaddpd (%rbx,%rax,8),%xmm0,%xmm0
0x0000000000402de3 <+531>:    vmovaps %xmm0,(%rbx,%rax,8)
0x0000000000402de8 <+536>:    add    $0x2,%rax
0x0000000000402dec <+540>:    cmp    %rax,%rbp
0x0000000000402def <+543>:    jg     0x402dd8 <avx::TestEigen(int, int)+520>

AVX 命令の利点の一つとして、従来の SSE 命令のオペランドの数が 2 つだったのに対し、AVX では 3 または 4 オペランドを使用できるようになっています。したがって、計算によっては mov の数を減らすことができるのですが、今回の例ではあまりそれは活きないのかもしれません。

ちなみに -O3 を -O2 に変更した場合、つまり SSE のスカラー演算と AVX のベクトル演算を比較すると、以下のように速度差がより顕著になります。この場合の Eigen のグラフを見るだけでも分かりますが、Eigen は最適化オプションが -O2 であってもベクトル化演算のコードを生成しています。

result-O2

手元に Sandy Bridge 以前の CPU があったので、ベクトル化の有無で速度の違いを比べてみました。コンパイル時に gcc によるベクトル化に失敗したので、コンパイル オプションによる違いはなく、明らかに速度が Eigen > valarray > vector の順になっています。

CPU: AMD Athlon(tm) 64 X2 Dual Core Processor 3800+
result-AMD

というわけで結果からの考察。

  1. 特に工夫をしない場合、R と Python は遅い。
  2. SIMD を意識せずに書いた単純な C++ のループでも、ベクトル化 + AVX によって明白に速度向上が見込める。 というわけで何らかの方法 (既存のライブラリ、intrinsic、gcc のオプション) で必ずベクトル化はしておくべき。
  3. 明白な理由がない限り、 vector を使って自分でループを書くよりは、valarray の四則演算を使っておいた方が利点が多い。ベクトル化 + AVX を使った場合に両者の速度差はほとんどない上、valarray の方がコードが見やすい。 ただし、複雑な演算の時に、vector を使った方がループを工夫できる可能性はあるかも。
  4. Eigen は gcc の自動ベクトル化に頼らず自力でベクトル演算を使える。ただし (現バージョンの 3.2.7 では) AVX を有効にしても ymm レジスタを使わない。Sandy Bridge 以前の CPU では Eigen を使えば間違いない。Sandy Bridge 以降だと、valarray で間に合う演算であれば valarray を使う方がよい。

ただし上記は、SIMD を特に意識せずに C++ を書いた場合なので、intrinsic などを使えばもっと効率よくベクトル化の恩恵を受けられるかもしれません。もちろん、計算アルゴリズムやベクトルのサイズなどが大きく関係してくるので、一概にどれがベストと言い難いのですが、少なくとも、配列要素のアラインメントも特に指定しない単純なループでさえ、自動ベクトル化と AVX はそれなりに効果があるということが分かりました。この速度比較は、それだけで別のテーマになりますが今回はここまでにして、オイラー法は valarray でやってみようかな・・。

Skylake の Xeon モデルには 512bit レジスターを使う AVX3 ユニットが搭載されているらしい。これは使ってみたい。

MNIST with MDS in Python

前回の続きで、MNIST データを MDS で分析してみました。下記のブログの PCA の次のセクション "Optimization-Based Dimensionality Reduction" の最初のグラフ "Visualizing MNIST with MDS" の部分だけです。

Visualizing MNIST: An Exploration of Dimensionality Reduction – colah’s blog
http://colah.github.io/posts/2014-10-Visualizing-MNIST/

発想は PCA よりも単純で、今度は個々の MNIST データの 784 次元空間上の位置は考慮せず、2 つのデータのユークリッド距離をなるべく維持できるように全ての点を 2 次元平面上にプロットし直す、というものです。ただし、元の距離をなるべく維持する点の集合を PCA のようにまともに計算するのはおそらく不可能なので、任意の 2 点間を元の距離と同じ長さを自然長とするばねで繋いで、それを二次元空間のランダムな位置に投げ込み、物理法則に任せることで最適解を得る、というのがよくある方法のようです。以下の日本語の Wikipedia でも言及されているアプローチとほぼ同じです。

力学モデル (グラフ描画アルゴリズム) – Wikipedia
https://ja.wikipedia.org/wiki/%E5%8A%9B%E5%AD%A6%E3%83%A2%E3%83%87%E3%83%AB_(%E3%82%B0%E3%83%A9%E3%83%95%E6%8F%8F%E7%94%BB%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0)

というわけで、今回はこれを Python で書いてみることにしました。Python 超初心者な上、エディタとして初めてまともに Emacs を導入してみたので、悪戦苦闘の毎日。しかもこれから書くように結果が微妙なことこの上ない。

R ではなく Python を選んだ理由は、matplotlib というモジュールを使ってインタラクティブなプロットを作るのが簡単そうだったからです。それと、R と SciPy はよく対比されているのでどちらも使ってみたかった、さらに言えば Python とは友達になっておいたほうが良さそうだった、などなど。

今回の環境はこんな感じ。

  • Ubuntu 15.10 "wily"
  • Python 2.7.10
  • matplotlib (1.4.2)
  • scipy (0.14.1)

というわけで書いたコードがこれ。

msmania/ems at purepython · GitHub
https://github.com/msmania/ems/tree/purepython

プロジェクト名の Euler Method というのは、微分方程式の近似解を得るための方法の一つです。実装は簡単だが精度が悪い、という理由から、基本的には使うべきではないらしい。他の方法としては、ルンゲ=クッタ法とかがあります。

オイラー法 – Wikipedia
https://ja.wikipedia.org/wiki/%E3%82%AA%E3%82%A4%E3%83%A9%E3%83%BC%E6%B3%95

以下のページによると、「『この問題を解く場合、4次 のルンゲクッタだな』と一言いって、プログラムを書き始めると、出来るなと 思われること間違いなし」 らしいです。覚えておこう。

2 数値計算法
http://www.ipc.akita-nct.ac.jp/yamamoto/lecture/2003/5E/lecture_5E/diff_eq/node2.html

オイラー法とは、微小時間における速度、位置の変化を一次近似するだけです。上で紹介した Wikipedia の 「力学モデル」 で出てきた擬似コードがまさにそれです。

今回の Python コードで実装したのは、フックの法則から導き出される運動方程式と、ノードの位置を収束させるための摩擦力だけで、クーロンの法則によるノード間の反発力は完全に無視しました。また、フックの法則において弾性限度はないものとしています。動摩擦係数と静摩擦係数も区別していません。

ems.py の中の Field クラスに、ノード全ての速度と位置を行列として保持させて、オイラー法は行列計算として実装しています。5 フレームごとに計全体の運動エネルギーを計算して、特定の値を下回ったらシミュレーションを停止するようにもしています。

matplotlib の animation を使うと、H.264 形式の動画が簡単に作れるので、これでシミュレーション結果をアップロードしてみました。

MNIST with MDS in Python (N=200)
http://msmania.github.io/ems/purepython.htm

Windows だと、Internet Explorer と Firefox ではこの H.264 動画が再生できません。Chrome だと問題ないみたいです。Ubuntu 上の Firefox も問題なしでした。Windows Media Player 付属のデコーダーが対応していないらしく、Firefox は Windows のデコーダーに依存しているっぽいです。動画も暗号化アルゴリズムも独自実装している Chrome は便利。

image.

MNIST のデータは N=10,000 もしくは 60,000 もあるのに、なんで N=200 なんだ、という話ですが、動作が遅すぎて計算できませんでした。上の N=200 の動画を作るだけでも、手元の環境でそれぞれ 6 分ぐらいかかります。ためしに N=1000 にしたところ、数時間経っても終わらなかったので断念。N=10,000 なんか 1 年経っても終わらないような。使っている CPU は Intel(R) Core(TM) i5-2520M CPU @ 2.50GHz で、アーキテクチャは Sandy Bridge です。

john@ubuntu1510:~/Documents/ems$ export N=200
john@ubuntu1510:~/Documents/ems$ python mnist.py 200.mp4
/usr/lib/python2.7/dist-packages/matplotlib/collections.py:571: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str(‘face’):
time = 360.681 sec
Done!

今回のモデルの場合、ばねのポテンシャル エネルギーが全て摩擦熱として失われるだけなので、系の運動エネルギーを毎回計算しなくてもシミュレーションの停止条件は求められそうです。しかし、試しにエネルギー計算をせずに 400 フレーム (= 20 秒) を生成させてみたところ、結局 7 分かかったので、エネルギー計算を省略しても問題は解決しません。

Python は書き方によって速度が大きく変わるらしいので、ひょっとするとコードを工夫すれば多少は改善するかもしれません。あとは Cython などの高速化ライブラリを使うとかだろうか。

Cython: C-Extensions for Python
http://cython.org/

しかし、ここまでやってようやく気づいてしまったのは、Field クラスは C++ で書いてしまったほうがよかったかもしれない、ということ。Sandy Bridge だから行列計算に AVX 命令を使えば高速化できそう。というか Python の SciPy は AVX 使ってくれないのかな。ノードにかかる運動ベクトルをラムダ式のリストにしてリスト内包表記で計算できるなんてエレガント・・という感じだったのに。

ビット演算ではないけど、↓ の話題に近くなってきている。コンピューター将棋でも、各局面を 81 次元の点として考えれば MNIST と大して変わらないのかもしれない。

プログラムを高速化する話 | やねうら王 公式サイト
http://yaneuraou.yaneu.com/2015/03/19/%E3%83%97%E3%83%AD%E3%82%B0%E3%83%A9%E3%83%A0%E3%82%92%E9%AB%98%E9%80%9F%E5%8C%96%E3%81%99%E3%82%8B%E8%A9%B1/

まだ知見を広められそうなので、もう少し MDS の演算で粘ってみようと思います。

Cooking MNIST with PCA in R

TensorFlow の初心者用チュートリアルを最初から読むと、MNIST (= Mixed National Institute of Standards and Technology) データベースと呼ばれる、手書きのアラビア数字を 28×28 のモノクロ ビットマップで表した画像データが出てきます。

MNIST For ML Beginners
https://www.tensorflow.org/tutorials/mnist/beginners/index.html

チュートリアルの冒頭で、この MNIST データをこねくり回してみたという記事へのリンクがあります。

Visualizing MNIST: An Exploration of Dimensionality Reduction – colah’s blog
http://colah.github.io/posts/2014-10-Visualizing-MNIST/

内容以前に、このウェブページ自体が MathJax という JavaScript ライブラリを使っていたり、ウェブでこんなことできるのか、という凄い代物で驚きます。著者の Christopher Olah という研究者の CV が公開されていて読めますが、まさに若き天才って感じ。流行の分野はやはり天才が多いんだろうなぁ。いいなぁ。

話が逸れました。さて、この記事の内容ですが、MNIST を視覚化する複数の数学的手法を紹介しています。より具体的には、28×28 のビットマップデータを、28×28 = 784 次元空間における 1 つの座標とみなして、いかにして人間が認識できる形にするか、すなわち 3 次元データにマップできるか、ということを試みています。

この発想自体が既に目から鱗だったわけですが、今回の記事では、具体的に出てくる次元下げの最初の手法についてトレースした結果を紹介します。残りの方法についても追って試す予定です。TensorFlow どこいったんだよ、という気もしますが、初めに基礎理論を固めておけば、あとの理解が楽になる気がするので端折らずにやります。

MNIST with PCA in R

そして最初の手法。PCA です。Principal component analysis の略。日本語だと主成分分析。どちらかというと Computer Science よりは、統計学で出てくる多変量解析の一手法です。PCA の次に MDS という手法も出てきますが、これも多変量解析の一つであり、統計学の理解が機械学習においては重要な気がします。

今回の内容を調べるにあたって、統計の初歩 (実は恥ずかしながら、これを調べる前は分散、偏差って何だっけレベルだった・・) から主成分分析の数学的背景までわりと真面目に復習したので、それについても今回の記事の最後に少し触れようと思います。

さて PCA ですが、統計と言えば R 言語だろう、という安易な考えで、Christopher Olah 氏がやったことを R で確かめてみます。デバッグをするわけではないのであまり重要ではないですが、環境は↓です。

  • Ubuntu 15.10 "wily"
  • R version 3.2.2 (2015-08-14) — "Fire Safety"

TensorFlow のチュートリアルに書いてありますが、MNIST データは以下のサイトからダウンロードできます。もしくは、TensorFlow に含まれている python スクリプトを使ってもいいです。

MNIST handwritten digit database, Yann LeCun, Corinna Cortes and Chris Burges
http://yann.lecun.com/exdb/mnist/

ファイルは 4 つあります。ファイル名に images とあるのがビットマップ データ。labels とあるのは、それぞれのビットマップがどのアラビア数字 (0-9) を書いたものなのかを表しています。要するにビットマップの答です。t10k のサンプル数が 10,000、train のサンプル数が 60,000 です。二つのサンプルに分かれている理由についてもチュートリアルのページに書かれていますが、学習用のデータと、学習後に使うデータを分けることで、学習が効果的かどうかを判断できる、ということです。

john@ubuntu1510:~/Documents/MNIST_data$ ls -l
total 53672
-rw-rw-r– 1 john john  7840016 Nov 26 16:16 t10k-images-idx3-ubyte
-rw-rw-r– 1 john john    10008 Nov 26 16:16 t10k-labels-idx1-ubyte
-rw-rw-r– 1 john john 47040016 Nov 26 16:16 train-images-idx3-ubyte
-rw-rw-r– 1 john john    60008 Nov 26 16:16 train-labels-idx1-ubyte

ファイルのフォーマットについては、上記 Yann LeCun 氏のページの FILE FORMATS FOR THE MNIST DATABASE のセクションに書かれています。ヘッダーに次いで、配列を次元毎にべたーっと書き出しただけの形式です。

今回書いた R のコードは以下の通り。

readMNIST <- function(findex) {
    imagefiles <- c("t10k-images-idx3-ubyte", "train-images-idx3-ubyte")
    labelfiles <- c("t10k-labels-idx1-ubyte", "train-labels-idx1-ubyte")
    fname <- paste0("~/Documents/MNIST_data/", imagefiles[findex])
    f <- file(fname, "rb")
    seek(f, 4) # skip magic
    n <- readBin(f, integer(), size=4, n=1, endian="big")
    size.x <- readBin(f, integer(), size=4, n=1, endian="big")
    size.y <- readBin(f, integer(), size=4, n=1, endian="big")
    pixels <- readBin(f, integer(), size=1, n=n * size.x * size.y, signed=F)
    close(f)

    images <- data.frame(matrix(pixels, ncol=28*28, byrow=T))

    fname <- paste0("~/Documents/MNIST_data/", labelfiles[findex])
    f <- file(fname, "rb")
    seek(f, 4) # skip magic
    n <- readBin(f, integer(), size=4, n=1, endian="big")
    labels <- readBin(f, integer(), size=1, n, signed=F)
    close(f)

    l <- list(labels, images)
    names(l) <- c("labels", "images")
    l
}

drawMNIST <- function(mnist, idx) {
    title <- mnist$labels[idx]
    pixels <- as.integer(mnist$images[idx,])
    pixels <- matrix(pixels, 28, 28)
    dev.new(width=5, height=5)
    image(pixels[,28:1], col=grey(seq(0, 1, length = 256)), main=title)
}

pcaMNIST <- function(images, scale) {
    images.iszero <- apply(images, 2, function(col) {
        return (var(col)==0 & var(col)==0)
    })
    validData <- images[,names(subset(images.iszero, !images.iszero))]
    prcomp(validData, scale=scale)
}

drawRotation <- function(images, pca, range) {
    images.iszero <- apply(images, 2, function(col) {
        return (var(col)==0 & var(col)==0)
    })
    zcols <- names(subset(images.iszero, images.iszero))
    l <- length(range)
    zrows <- matrix(rep(0, length(zcols) * l), ncol=l)
    rownames(zrows) <- zcols
    pcs <- pca$rotation[,range]
    pcs <- rbind(pcs, zrows)
    rowindex <- as.integer(substring(rownames(pcs), 2))
    pcs <- pcs[order(rowindex),]
    nul <- apply(pcs, 2, function(col) {
        dev.new(width=5, height=5)
        image(matrix(col, 28, 28)[,28:1],
              col=c(hsv(2/3, seq(1, 0, length=10), 1),
                    hsv(1,   seq(0, 1, length=10), 1)))
    })
}

以下 4 つの関数を定義しています。

  • readMNIST – MNIST ファイルを R のデータ フレームとして読み込む関数
  • drawMNIST – 読み込んだビットマップのデータ フレームをグラフィカル表示する関数
  • pcaMNIST – PCA を実行する関数
  • drawRotation – PCA 結果から、主成分をグラフィカル表示する関数

まずはファイルを読み込んでビットマップを表示してみます。

> source("MNISTCook.R")
> mnist10k <- readMNIST(1)
> mnist60k <- readMNIST(2)
> drawMNIST(mnist60k, 1)
> drawMNIST(mnist60k, 2)
> drawMNIST(mnist60k, 3)
> drawMNIST(mnist60k, 4)

上記のコードで、t10k, train の両方のデータを読み込んだ後、train データに含まれる 60,000 枚のビットマップのうちの最初の 4 枚を表示します。こんな感じ↓。チュートリアルの最初のページにある画像と一致していますね。

01-images-60k

ちなみに、Raw データは、od コマンドで簡単に確かめることも出来ます。例えばこれが一枚目。

john@ubuntu1510:~/Documents$ od -tx1 -v -w28 -j16 -N784 MNIST_data/train-images-idx3-ubyte
0000020 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0000054 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0000110 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0000144 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0000200 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0000234 00 00 00 00 00 00 00 00 00 00 00 00 03 12 12 12 7e 88 af 1a a6 ff f7 7f 00 00 00 00
0000270 00 00 00 00 00 00 00 00 1e 24 5e 9a aa fd fd fd fd fd e1 ac fd f2 c3 40 00 00 00 00
0000324 00 00 00 00 00 00 00 31 ee fd fd fd fd fd fd fd fd fb 5d 52 52 38 27 00 00 00 00 00
0000360 00 00 00 00 00 00 00 12 db fd fd fd fd fd c6 b6 f7 f1 00 00 00 00 00 00 00 00 00 00
0000414 00 00 00 00 00 00 00 00 50 9c 6b fd fd cd 0b 00 2b 9a 00 00 00 00 00 00 00 00 00 00
0000450 00 00 00 00 00 00 00 00 00 0e 01 9a fd 5a 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0000504 00 00 00 00 00 00 00 00 00 00 00 8b fd be 02 00 00 00 00 00 00 00 00 00 00 00 00 00
0000540 00 00 00 00 00 00 00 00 00 00 00 0b be fd 46 00 00 00 00 00 00 00 00 00 00 00 00 00
0000574 00 00 00 00 00 00 00 00 00 00 00 00 23 f1 e1 a0 6c 01 00 00 00 00 00 00 00 00 00 00
0000630 00 00 00 00 00 00 00 00 00 00 00 00 00 51 f0 fd fd 77 19 00 00 00 00 00 00 00 00 00
0000664 00 00 00 00 00 00 00 00 00 00 00 00 00 00 2d ba fd fd 96 1b 00 00 00 00 00 00 00 00
0000720 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 10 5d fc fd bb 00 00 00 00 00 00 00 00
0000754 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 f9 fd f9 40 00 00 00 00 00 00 00
0001010 00 00 00 00 00 00 00 00 00 00 00 00 00 00 2e 82 b7 fd fd cf 02 00 00 00 00 00 00 00
0001044 00 00 00 00 00 00 00 00 00 00 00 00 27 94 e5 fd fd fd fa b6 00 00 00 00 00 00 00 00
0001100 00 00 00 00 00 00 00 00 00 00 18 72 dd fd fd fd fd c9 4e 00 00 00 00 00 00 00 00 00
0001134 00 00 00 00 00 00 00 00 17 42 d5 fd fd fd fd c6 51 02 00 00 00 00 00 00 00 00 00 00
0001170 00 00 00 00 00 00 12 ab db fd fd fd fd c3 50 09 00 00 00 00 00 00 00 00 00 00 00 00
0001224 00 00 00 00 37 ac e2 fd fd fd fd f4 85 0b 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0001260 00 00 00 00 88 fd fd fd d4 87 84 10 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0001314 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0001350 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0001404 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0001440

次に、PCA を実行します。R で PCA を行う関数は複数ありますが、ここでは prcomp 関数を使っています。ググって見つけたサイトが prcomp を使っていたからです。理由なんてそんなもんです。t10k と train のそれぞれのデータを PCA にかけます。train の方は 1 分ぐらいかかりました。

> pca10k <- pcaMNIST(mnist10k$images, F)
> pca60k <- pcaMNIST(mnist60k$images, F)

PCA を実行した後、算出された成分のうち第一主成分から第四主成分を 28×28 のビットマップにプロットしたものを見てみます。というのも、Christopher Olah 氏が第二主成分までのプロットを記事に載せているからです。

> drawRotation(mnist60k$images, pca60k, 1:4)
> drawRotation(mnist10k$images, pca10k, 1:4)

上記一行目のコードで、train データについてのプロットが表示されます。上 2 つの主成分は、Christopher Olah 氏の載せている結果と同じと考えてよさそうです。

02-pca-60k

次に t10k の方から導出した主成分です。なぜかベクトルの向きが反転 (青と赤が逆) したようですが、向きを逆にすれば trainデータから導出された主成分とほぼ同じと見て問題なさそうです。

03-pca-10k

この後、Christopher Olah 氏の記事では、MNIST データを第一主成分と第二主成分の 2 次元座標にプロットした結果を示した上で、あまり nicely じゃないと結論付けて MDS を使った分析に移ります。ちょっと気になったのは、彼が累積寄与率に関する議論をしていないことです。まあ彼からすれば PCA なんて当て馬で大して議論する気もなかったのでしょうが、せっかく R のデータが手元にあるので、寄与率を見てみます。

> summary(pca10k)$importance[,1:4]
                             PC1       PC2       PC3       PC4
Standard deviation     587.63781 509.20428 459.38799 431.82659
Proportion of Variance   0.10048   0.07544   0.06141   0.05426
Cumulative Proportion    0.10048   0.17592   0.23733   0.29158

> summary(pca60k)$importance[,1:4]
                             PC1       PC2       PC3       PC4
Standard deviation     576.82291 493.23822 459.89930 429.85624
Proportion of Variance   0.09705   0.07096   0.06169   0.05389
Cumulative Proportion    0.09705   0.16801   0.22970   0.28359

第二主成分までを採用したとして、寄与率は 17-18% であり、仮に第四主成分まで使っても 70% 以上の情報が失われることが分かります。巷の情報では、寄与率が 70 もしくは 80% を超えるように主成分を使うみたいなので、確かに 2 次元プロットが微妙なのも頷けます。と言っても、元データが 784 次元もあり、変換前の次元 2 つでは 2/784 = 0.255% しか寄与率がないことを考えると PCA はやっぱり強力、ということにはなります。

ここからは余談ですが、最初 pcaMNIST 関数を書いたとき、何をとち狂ったか、MNIST データを標準化した上で PCA にかけていました。すると以下のような全然違う主成分が得られます。この主成分がどんな意味を持つのかについては特に何も考察していませんが、元データがビットマップであり、各次元の単位は共通であることを考慮すると、標準化ダメ絶対、でしょうかね。

04-pca-60k-scale

おまけ – PCA とは何か。どうやって計算するのか。

ここからは TensorFlow はもとより、MNIST とも関係なくなります。

PCA、すなわち主成分分析は、ちょっとググって見たところ一般に広く知られている手法で、需要予測などのマーケティングや、野球選手の戦力分析など、幅広い分野で使われているようです。一般的には次元下げが目的ですが、数学的には単なる座標変換であり、変換後の軸に対する振れ幅、すなわち分散が最大になるように決めた変換です。イメージとしては、以下のサイトの冒頭の図が分かりやすいです。

高校数学の基本問題 – Excelを用いた主成分分析
http://www.geisya.or.jp/~mwm48961/statistics/syuseibun1.htm

一次変換によって N 次元のデータを別の N 次元に変換するわけですが、各次元がデータに寄与する割合を変えて、なるべく少ない次元で多くの情報量を保持するように変換を決めて、その次元を主成分、と呼んでいるわけです。

面白いことに、主成分は、元データの分散共分散行列の固有ベクトルとして求まります。これを数学的に検証している日本語のサイトはさすがに少なかったのですが、以下の 2 つのサイトがとても参考になりました。

統計学復習メモ10: なぜ共分散行列の固有ベクトルが単位主成分なのか (Weblog on mebius.tokaichiba.jp)
http://ynomura.dip.jp/archives/2009/03/10.html

統計学入門-第16章
http://www.snap-tck.com/room04/c01/stat/stat16/stat1601.html

実際にこの計算を自分でも書き下してみたところ、PCA というのは数学的に美しいのが実感できました。ただし、計算の一番のポイントであるラグランジュの未定乗数法の証明はスキップしました。だって以下のサイトに、証明は単純ではないと書いてあるし。

EMANの物理学・解析力学・ラグランジュの未定乗数法
http://homepage2.nifty.com/eman/analytic/lag_method.html

ラグランジュの未定乗数法を除いても、大学で線形代数は苦手だったのでかなり理解するのに苦労しました。あと行列の偏微分はたぶん初見。それ以前に手で数式を解く感覚が完全に失われていたので取り戻すのに時間がかかりましたが。たまには数学をやって頭を柔らかくしないと駄目だなぁ・・あれ、また話が逸れている。

数学的な証明をこのブログに載せるのは面倒なのでパスするとして、では本当に分散共分散行列から主成分が求まるのかどうかを R で確認してみます。というのも、最初 prcomp 関数の戻り値の見方が分かりづらくて苦労したからです。

MNIST は次元が多すぎて見づらいので、R のサンプル データである iris を使うことします。このデータは、アヤメの 3 品種 (setosa, versicolor, virginica]) について、がく片 (Sepal) と花びら (Petal) の幅と長さをそれぞれ計測したものです。つまりデータの次元は 4 で、サンプル数は 150 です。

> str(iris)
‘data.frame’:   150 obs. of  5 variables:
$ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 …
$ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 …
$ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 …
$ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 …
$ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 …
> iris[1:4,]
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa

prcomp 関数を実行します。また、PCA の結果は biplot で、主成分と主成分得点を同じ図にプロットすることが多いみたいなので、そのまま実行します。

> pca <- prcomp(iris[,1:4], scale=T)
> biplot(pca)

05-iris-biplot

これだけだとよく分かりませんが、実は品種ごとにデータが綺麗に分かれています。こちらのサイトで同じプロットを色分けしたものが掲載されています。

主成分分析が簡単にできるサイトを作った – ほくそ笑む
http://d.hatena.ne.jp/hoxo_m/20120106/p1

ここで紹介したいのは、prcomp が返した pca というのは何か、という点です。str 関数で見てみます。

> str(pca)
List of 5
$ sdev    : num [1:4] 1.708 0.956 0.383 0.144
$ rotation: num [1:4, 1:4] 0.521 -0.269 0.58 0.565 -0.377 …
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
  .. ..$ : chr [1:4] "PC1" "PC2" "PC3" "PC4"
$ center  : Named num [1:4] 5.84 3.06 3.76 1.2
  ..- attr(*, "names")= chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
$ scale   : Named num [1:4] 0.828 0.436 1.765 0.762
  ..- attr(*, "names")= chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
$ x       : num [1:150, 1:4] -2.26 -2.07 -2.36 -2.29 -2.38 …
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:4] "PC1" "PC2" "PC3" "PC4"
– attr(*, "class")= chr "prcomp"

どうやら、5 つの要素を持つリストになっているようです。center, scale は iris データを標準化したときの値、rotation が主成分の行列、sdev が成分ごとの標準偏差、x が主成分得点になっています。pca には累積寄与率が含まれていませんが、summary 関数を実行することで、pca の sdev 要素を元に累積寄与率が計算されます。2 つの主成分で寄与率が 96% 近くあることが分かります。

> summary(pca)
Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

試しに summary を使わずに寄与率を求めるとすれば、こうでしょうか。iris は標準化しておいたので、各次元の分散は 1 であり、4 で割れば全体を 1 としたときの累積寄与率が導出されます。

> pca$sdev * pca$sdev / 4
[1] 0.729624454 0.228507618 0.036689219 0.005178709

pca で重要なのは主成分です。主成分は N 次元データを一次変換する行列なので、必ず NxN の正方行列になります。

> pca$rotation
                    PC1         PC2        PC3        PC4
Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863
Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

この主成分が、分散共分散行列の固有ベクトルなのでした。ということで実際に計算します。

> sdata <- scale(iris[,1:4])
> S <- cov(sdata)
> pca2 <- eigen(S)
> pca2
$values
[1] 2.91849782 0.91403047 0.14675688 0.02071484

$vectors
           [,1]        [,2]       [,3]       [,4]
[1,]  0.5210659 -0.37741762  0.7195664  0.2612863
[2,] -0.2693474 -0.92329566 -0.2443818 -0.1235096
[3,]  0.5804131 -0.02449161 -0.1421264 -0.8014492
[4,]  0.5648565 -0.06694199 -0.6342727  0.5235971

> sqrt(pca2$values)
[1] 1.7083611 0.9560494 0.3830886 0.1439265

pca$rotation と同じ行列が得られました。また、固有値の平方根を取ったものが pca$sdev と等しくなっています。つまり分散共分散行列の固有値が主成分の分散になっていることも確かめられました。最後に、主成分を使って主成分得点を求めてみます。150 サンプル全部を計算してもしょうがないので、最初の 5 サンプルを取って・・・

> scale(iris[,1:4])[1:5,] %*% pca2$v
pca2$values   pca2$vectors
> scale(iris[,1:4])[1:5,] %*% pca2$vectors
          [,1]       [,2]        [,3]        [,4]
[1,] -2.257141 -0.4784238  0.12727962  0.02408751
[2,] -2.074013  0.6718827  0.23382552  0.10266284
[3,] -2.356335  0.3407664 -0.04405390  0.02828231
[4,] -2.291707  0.5953999 -0.09098530 -0.06573534
[5,] -2.381863 -0.6446757 -0.01568565 -0.03580287
> pca$x[1:5,]
           PC1        PC2         PC3         PC4
[1,] -2.257141 -0.4784238  0.12727962  0.02408751
[2,] -2.074013  0.6718827  0.23382552  0.10266284
[3,] -2.356335  0.3407664 -0.04405390  0.02828231
[4,] -2.291707  0.5953999 -0.09098530 -0.06573534
[5,] -2.381863 -0.6446757 -0.01568565 -0.03580287

行列積の結果、pca$x と等しい値が得られることが確認できました。

ちなみに、PCA について書かれた英語版の wikipedia がやたら専門的で詳しいのですが、そのページの Derivation of PCA using the covariance method > Iterative computation セクションのところで、実際の高次元データで分散共分散行列を使って主成分を求めるのは効率が悪い、と書かれています。うまくアルゴリズムを書くと、共分散行列を全部計算しなくても主成分が求まるとかなんとか。

Principal component analysis – Wikipedia, the free encyclopedia
https://en.wikipedia.org/wiki/Principal_component_analysis

Export Excel to XML using VBA

CSV などのテーブル データから XML を作りたいことが多々あります。Excel にはもともと XML へのエクスポートを行なう機能がありますが、望みどおりに動いてくれないため、VBA でマクロを書いてみましたのでコードを公開します。しかしこのご時勢では XML より JSON 派のほうが多いんですかね・・。

ちなみに開発環境は、Windows 7 SP1 上の Excel 2010 です。Office 2013 はどうも好きになれない。Windows 10 は入れたいのだが、諸事情でこの PC には入れられない・・。

さて、無味乾燥なサンプルだと面白くないので、楽天 API 経由で取ってきた商品データを使うことにします。といっても簡単で、PowerShell で以下のコマンドレットを実行するだけです。クエリ URL の <APPID> のところは、各自のアプリ ID を入れてください。持っていない人は、https://webservice.rakuten.co.jp/ でアカウントを登録して作ることができます。

PS> $Client = New-Object System.Net.WebClient
PS> $QueryUrl = ‘
https://app.rakuten.co.jp/services/api/BooksTotal/Search/20130522?format=xml&keyword=%E5%AF%BF%E5%8F%B8&booksGenreId=000&hits=10&applicationId=<APPID>’
PS> $Client.DownloadFile($QueryUrl, ‘E:\MSWORK\Generator\Sushi.xml’)

上記は、楽天ブックスから "寿司" というキーワードで 10 件の検索結果を取ってきて、Sushi.xml という XML ファイルで保存するものです。これを Excel で開くと、以下のダイアログが出てくるので、「XML テーブルとして開く」 を選んで OK をクリックします。

image

スキーマ情報がないので以下のダイアログが表示されます。OK をクリックします。

image

テーブルができます。

image

ちなみにこのまま、XML ソースの画面にある [エクスポートする対応付けの確認..] をクリックすると、"例外的なデータ" (= denormalized data) が存在するとかでエクスポートできません。意味不明・・。

image

と、いうわけで、このようなテーブルを XML にエクスポートするためのマクロを書きます。

元の Sushi.xml の構造は、大まかには以下のようになっています。検索結果の各アイテムの情報が root  > items > item というノードに保存されています。簡単のため、items 以外の count や page といった検索のメタ情報は Excel のシートから除外することにします。

<root>
  <count>4926</count>
  <page>1</page>
  <first>1</first>
  <last>10</last>
  <hits>10</hits>
  <carrier>0</carrier>
  <pageCount>100</pageCount>
  <Items>
    <Item>
      <title>検索結果1</title>
      …
    </Item>
    <Item>
      <title>検索結果2</title>
      …
    </Item>
  </Items>
  <GenreInformation />
</root>

列を削除して XMLGenerator.xlsm という名前で保存します。マクロを使うので、拡張子は xlsm という形式にします。

image

リボンに [開発] タブが表示されていない人は、以下の手順で有効にしてください。

[開発] タブを表示する – Office のサポート
https://support.office.com/ja-jp/article/-%E9%96%8B%E7%99%BA-%E3%82%BF%E3%83%96%E3%82%92%E8%A1%A8%E7%A4%BA%E3%81%99%E3%82%8B-e1192344-5e56-4d45-931b-e5fd9bea2d45

コードを書く前にもう一点確認。テーブル内のセルのひとつにカーソルを置いた状態で [デザイン] タブを開いてテーブル名を確認します。以下の画面キャプチャだと、左上の方に "テーブル1" と入力されている部分がそれです。

image

では、[開発] タブにある [Visual Basic[ をクリックして VBA の画面を起動します。ほとんど VB6 と同じですね。

XML の読み書きには MSXML を使いたいので、まずは参照設定を追加します。メニューから [ツール] > [参照設定] を選びます。

image

ライブラリのリストから Microsoft XML というのを探して、チェックをつけてから OK をクリックします。v3.0 と v6.0 の両方がありましたが、新しい方の v6.0 を選んでおきます。

image

コードを追加します。Sheet1 または ThisWorkbook のどちらかのモジュールに紐付ける、もしくは新規にモジュールを作る、という選択肢がありますが、今回は ThisWorkbook にコードを追加します。ThisWorksheet を右クリックして、コンテキスト メニューの [コードの表示] を選びます。

image

コードをごにょごにょ書きます。先ほど確認した "テーブル1" というテーブル名をそのまま使っています。

Option Explicit
Private Function AddElementToNode(xmlObj As MSXML2.DOMDocument60, rootNode As MSXML2.IXMLDOMNode, elementName As String, elementValue As String) As MSXML2.IXMLDOMNode
    Dim newElement As MSXML2.IXMLDOMElement
    Set newElement = xmlObj.createElement(elementName)
    newElement.Text = elementValue
    Set AddElementToNode = rootNode.appendChild(newElement)
End Function

Private Sub AddAttributeToElement(xmlObj As MSXML2.DOMDocument60, rootElement As MSXML2.IXMLDOMElement, attributeName As String, attributeValue As String)
    Dim newAttribute As MSXML2.IXMLDOMAttribute
    Set newAttribute = xmlObj.createAttribute(attributeName)
    newAttribute.Text = attributeValue
    rootElement.setAttributeNode newAttribute
End Sub

Sub ExportTableToXml()
    Dim sheet1 As Worksheet
    Set sheet1 = ThisWorkbook.Worksheets(1)
   
    Dim outputFile As String
    outputFile = ThisWorkbook.Path & "\001.xml"
   
    Dim msg As String
    Dim result As String
    msg = "If [" & outputFile & "] exists, overwrite it?"
    result = MsgBox(msg, vbYesNo, "Just to make sure…")
   
    If result = vbNo Then Exit Sub
   
    Dim table1 As range
    Set table1 = range("テーブル1")
              
    Dim xmlObj As MSXML2.DOMDocument60
    Set xmlObj = New MSXML2.DOMDocument60

    xmlObj.async = False
    xmlObj.setProperty "SelectionLanguage", "XPath"
   
    xmlObj.appendChild xmlObj.createProcessingInstruction("xml", "version=’1.0′ encoding=’utf-8’")

    Dim rootElement As MSXML2.IXMLDOMElement
    Set rootElement = xmlObj.createElement("root")
   
    Dim itemsContainer As MSXML2.IXMLDOMElement
    Set itemsContainer = xmlObj.createElement("items")
    rootElement.appendChild itemsContainer
   
    Dim row As Integer
    For row = 1 To table1.Rows.Count
        Dim itemElem As MSXML2.IXMLDOMElement
        Set itemElem = xmlObj.createElement("item")
       
        Dim col As Integer
        For col = 1 To table1.Columns.Count
            Dim colName As String
            colName = table1.Cells(0, col)
           
            If LCase(colName) = "title" Then
                AddAttributeToElement xmlObj, itemElem, colName, table1.Cells(row, col)
            Else
                AddElementToNode xmlObj, itemElem, colName, table1.Cells(row, col)
            End If
        Next

        itemsContainer.appendChild itemElem
    Next
   
    xmlObj.appendChild rootElement
    xmlObj.Save outputFile
   
    msg = "[" & outputFile & "] has been created/updated."
    MsgBox msg, vbOKOnly, "Done."
End Sub

プロシージャ呼び出しにはカッコをつけない、オブジェクトを New するときは Set を使う、などの意味不明な文法に苦労しましたが、まあこれで動きます。

ユーザビリティを考えるのであれば、ワークシートのほうにボタンを配置して、そこからマクロを呼ぶようにするのが自然でしょうか。というわけで、ワークシートに戻って [開発] タブから [デザイン モード] に入ってボタンを挿入します。[マクロの登録] ダイアログが出てくるので、先ほど作ったプロシージャ ExportTableToXml() を選択して OK をクリックします。

image

こんな感じになりました。

image

ボタンをクリックすると出力ファイルのパスとともに確認ダイアログが出るようにしています。出力ディレクトリは Excel ファイルと同じところ、ファイル名は 001.xml で固定です。

image

出力がうまくいけば、以下のようなダイアログが出ます。

image

出力された XML をブラウザー (ここでは Chrome) で開いて確認します。もとの Sushi.xml とほぼ同じです。変えた部分は、root 直下の items 以外のノードを削除したことと、本のタイトルを <title> ノードの値からく、<item> の name 属性に移動したことです。

image

16bit programming in Windows 95

前回 MS-DOS で QuickAssembler を動かしてみたのはいいものの、32bit のアセンブリ言語をコンパイルできないことに気づいてしまったので、結局 Windows を入れることに。

しかし Hyper-V で Windows 95 のインストーラーを実行すると、仮想マシンがハングしてしまう現象に遭遇。Windows 95 のインストーラーはマウスに対応しているため、このハングの原因はおそらく MS-DOS の MOUSE コマンドを実行するとハングしてしまう現象と同じ気がします。Hyper-V 使えないじゃん、ということで ESXi 5.1 で試してみたら、こちらは問題なし。やっぱり開発用途には VMware に軍配が上がります。

OS のインストール後、開発環境として最初は何も考えずに Visual C++ 6.0 を入れてみたのですが、Visual C++ 6.0 は 16bit アプリケーションの開発に対応していなかった、ということに初めて気づき、以下の 2 つを入れ直しました。何とも世話の焼ける・・

  • VIsual C++ 1.52
  • MASM 6.11

セットアップ中の画面の一部を紹介します。まずは Windows 95。実に懐かしい。

image image

image

十数年ぶりに触ってみると、今では当然のように使っている機能がまだ実装されていなくて驚きます。例えば・・

  • コマンド プロンプトが cmd.exe じゃなくて command.com
  • コマンド プロンプトで矢印キー、 Tab キーが使えない
  • コマンド プロンプトで、右クリックを使ってクリップボードからの貼り付けができない
  • メモ帳で Ctrl+A、Ctrl+S が使えない
  • エクスプローラーにアドレスバーがない

以下は Visual C++ と MASM のインストール画面を抜粋したもの。

image image

MASM 6.11 のインストーラーは CUI。

image

DOS, Windows, WIndows NT の全部に対応しているらしい。

image

MASM はいろいろと注意事項が多い。

image image

image

インストールが完了したら、環境変数の設定を行なうためのバッチ ファイルを作ります。サンプルのスクリプトがインストールされているので、それを組み合わせるだけです。

@echo off
set PATH=C:\MSVC\BIN;C:\MASM611\BIN;C:\MASM611\BINR;%PATH%
set INCLUDE=C:\MSVC\INCLUDE;C:\MSVC\MFC\INCLUDE;C:\MASM611\INCLUDE;%INCLUDE%
set LIB=C:\MSVC\LIB;C:\MSVC\MFC\LIB;C:\MASM611\LIB;%LIB%
set INIT=C:\MSVC;C:\MASM611\INIT;%INIT%

MASM、VC++ ともに NMAKE を持っていますが、VC++ 1.52 に入っている NMAKE のバージョンが新しかったので、VC++ を先頭にしました。

環境ができたところで、今回は 「はじめて読む 486」 の例題を試してみます。本文に掲載されている例題は、Borland C++ 3.1/Turbo Assembler 3.2 用の文法になっているため、細かいところは VC++/MASM 用に書き換えないといけません。というわけでこんな感じ。

まずは C ソースファイル main.c。

#include <stdlib.h>
#include <stdio.h>

extern short GetVer();
extern void RealToProto();
extern void ProtoToReal();

void main(int argc, char **argv) {
    printf("Hello, MS-DOS%d!\n", GetVer());
    RealToProto();
    ProtoToReal();
    printf("Successfully returned from Protected mode.\n");
    exit(0);
}

次がアセンブリ utils.asm。

.386
.MODEL small
.code

;* GetVer – Gets DOS version.
;*
;* Shows:   DOS Function – 30h (Get MS-DOS Version Number)
;*
;* Params:  None
;*
;* Return:  Short integer of form (M*100)+m, where M is major
;*          version number and m is minor version, or integer
;*          is 0 if DOS version earlier than 2.0

_GetVer  PROC

        mov     ah, 30h                 ; DOS Function 30h
        int     21h                     ; Get MS-DOS version number
        .IF     al == 0                 ; If version, version 1
        sub     ax, ax                  ; Set AX to 0
        .ELSE                           ; Version 2.0 or higher
        sub     ch, ch                  ; Zero CH and move minor
        mov     cl, ah                  ;   version number into CX
        mov     bl, 100
        mul     bl                      ; Multiply major by 10
        add     ax, cx                  ; Add minor to major*10
        .ENDIF
        ret                             ; Return result in AX

_GetVer  ENDP

public _RealToProto
_RealToProto    proc    near
                push bp
                mov bp, sp
                ;
                mov eax, cr0
                or eax, 1
                mov cr0, eax
                ;
                jmp flush_q1
flush_q1:
                pop bp
                ret
_RealToProto    endp

public _ProtoToReal
_ProtoToReal    proc    near
                push bp
                mov bp, sp
                ;
                mov eax, cr0
                and eax, 0fffffffeh
                mov cr0, eax
                ;
                jmp flush_q2
flush_q2:
                pop bp
                ret
_ProtoToReal    endp

        END

最後に Makefile。QuickC のときよりは現在の書式に近づきましたが、相変わらずリンカへの入力の渡し方がおかしい。

PROJ = TEST
USEMFC = 0
CC = cl
ML = ml
CFLAGS =/nologo /W3 /O /G3
LFLAGS =/NOLOGO /ONERROR:NOEXE
AFLAGS =
LIBS =
MAPFILE =nul
DEFFILE =nul

all: $(PROJ).EXE

clean:
    @del *.obj
    @del *.exe
    @del *.bnd
    @del *.pdb

UTILS.OBJ: UTILS.ASM
    $(ML) $(AFLAGS) /c UTILS.ASM $@

MAIN.OBJ: MAIN.C
    $(CC) $(CFLAGS) /c MAIN.C $@

$(PROJ).EXE:: MAIN.OBJ UTILS.OBJ
    echo >NUL @<<$(PROJ).CRF
MAIN.OBJ +
UTILS.OBJ
$(PROJ).EXE
$(MAPFILE)
$(LIBS)
$(DEFFILE)
;
<<
    link $(LFLAGS) @$(PROJ).CRF
    @copy $(PROJ).CRF $(PROJ).BND

プログラムはとても単純で、前回と同様に int 21 で DOS のバージョンを表示してから、コントロール レジスタ 0 の PE ビットを変更して特権レベルを変更し、また戻ってくる、というものです。ただし上記のコードでは、C、アセンブリのコンパイルは通るものの、以下のリンクエラーが出てうまくいきません。

UTILS.OBJ(UTILS.ASM) : fatal error L1123: _TEXT : segment defined both 16- and 32-bit

image

ここで 2 時間ぐらいハマりました。リンクエラーの原因は単純で、VC++ 1.52 のコンパイラは 16bit コードを生成しているのに対し、MASM は 32bit コードを生成しているためです。エラー メッセージの通り、16bit と 32bit の 2 つのコード セグメントを含む実行可能ファイルを生成できないためのエラーです。

このプログラムで確かめたいのは、16bit リアル モードから PE ビットを変更してプロテクト モードに移行し、再び 16bit リアル モードに戻ってくる動作です。したがって生成されるべき EXE は 16bit アプリケーションであり、VC++ の生成する main.obj は問題なく、MASM が 16bit コードを生成するように指定したいところです。

MASM が 32bit コードを生成する理由は、utils.asm の先頭の .386 で CPU のアーキテクチャを指定しているためです。これがないと、32bit レジスタの eax などが使えません。CPU のアーキテクチャを指定したことで、その後の .code セグメント指定が自動的に 32bit コード セグメントになってしまっています。

32bit 命令を有効にしつつ、セグメントは 16bit で指定する方法が見つかればよいわけです。NASM や GNU Assembler だと簡単に指定できるようですが・・いろいろと探して、以下のフォーラムを見つけました。10 年前に同じ悩みを抱えている人がいた!

Link Fatal Error 1123
http://www.masmforum.com/board/index.php?PHPSESSID=786dd40408172108b65a5a36b09c88c0&topic=1382.0

ファイルの先頭で .MODEL と CPU 指定を入れ替えればいいらしい。こんなん知らんて。

.MODEL small
.386

.code

;* GetVer – Gets DOS version.
;*
;* Shows:   DOS Function – 30h (Get MS-DOS Version Number)
;*
;* Params:  None
;*
;* Return:  Short integer of form (M*100)+m, where M is major
;*          version number and m is minor version, or integer
;*          is 0 if DOS version earlier than 2.0

_GetVer  PROC

以下、ずっと同じなので省略

これで無事ビルドが成功し、実行できました。文字が出力されているだけなので、本当にプロテクト モードに変わったかどうか疑わしくはありますが、まあ大丈夫でしょう。16bit の道は険しい。

image