포트란만큼 빠른 C++ 벡터 라이브러리를 파이썬에서 사용하기

GMES의 기하 관련 모듈은 원래 파이썬으로 구현되어 있었는데, 커다란 3차원 구조를 다루기에는 속도가 너무 느렸다. 속도 개선을 위해서 기하 관련 부분을 C++로 포팅하기로 결심했다. 그리곤, 기하 루틴에 많이 들어 있는 벡터 관련 함수를 어떻게 구현해야 할지 고민에 빠졌다. 파이썬에선 벡터의 자료형으로 numpy의 ndarray를 쓰면 되지만, C++에선 어떻게 구현하는 게 가장 좋을까? 시스템 호환성과 관리의 편리성을 고려해 널리 쓰이는 방법 중에서 선택할 수 있는 것들은 아래와 같다. 사용할 벡터는 실수를 원소로 갖는 3차원 벡터로 한정한다.

1. double *v, int size
2. std::array<double, 3>
3. boost::numeric::ublas::bounded_array<double, 3>
4. blitz::TinyVector<double, 3>

1번은 가장 쉽고 가장 빠른 방법이다. SWIG를 쓴다면, numpy.i를 이용해서 numpy와 연동하기도 쉽다. 하지만 벡터 연산을 일일이 직접 구현해야 하고, 함수의 인자 부분이 못생긴 형태가 되며, 왠지 C++답지 않은 방법이다.

2번은 표준 템플릿 라이브러리에 포함된 방법이며, 훨씬 C++스럽다. 하지만 1번처럼 벡터 관련 연산을 일일이 구현해야 하며 1번보다 느리다는 문제가 있다.

3번은 boost 라이브러리에 들어있는 컨테이너로 사실상 표준이라 할 수 있으며, 벡터 관련 연산도 모두 구현되어 있다. 하지만 속도가 조금 느리다는 단점이 있다. 간단한 벡터 연산뿐만 아니라 다양한 수학연산을 사용해야 한다면 boost 라이브러리를 선택하는 것이 좋다.

4번은 blitz++ 라이브러리에 들어있는 컨테이너로 포트란에 견줄 만한 속도가 장점이다. 과학계산 분야에선 거의 표준이라 할 수 있다.

보는 바와 같이 3번이나 4번이 좋은 선택이다. GMES 구현에는 4번을 선택했는데, 그 이유는 GMES가 SWIG를 이용해서 파이썬과 연동을 하고 있기 때문이다. boost는 자체적인 파이썬 연동 방식을 제공하고 있어서 한 프로그램에서 SWIG와 boost가 동시에 사용하는 것은 뭔가 자연스럽지 못한 느낌이 들었다. 그뿐만 아니라, blitz++가 아직은 더 빠르다.

blitz::TinyVector<double, 3>을 입력이나 출력으로 사용하는 함수를 파이썬에 연결하려면 어떻게 해야 할까? 우선 아래와 같은 간단한 C++ 함수를 생각해보자.

// util.hh

#include <blitz/tinyvec.h>

namespace util {
  typedef blitz::TinyVector<double, 3> Vector3;
  
  double accumulate(const Vector3& v);
  void weight(Vector3& v, double s = 1);
}

util::accumulate는 와 벡터 성분의 합을 구하는 함수이고, util::weight는 벡터에 상수를 곱하는 함수이다. 이들 함수는 아래와 같이 별로 어렵지 않게 구현해볼 수 있다.

// util.cc

#include "util.hh"

using namespace util;

double
util::accumulate(const Vector3& v)
{
  return v[0] + v[1] + v[2];
}

void 
util::weight(Vector3& v, double s)
{
  v *= s;
}

util::accumulate는 상수 인자를 받으니, 파이썬 측의 시퀀스(tuple, list, numpy.ndarray 등)를 C++ 쪽으로 변환해 주기만 하면 된다. 반면, util::weight의 인자는 값이 변경되므로 util::Vector 형을 다시 파이썬의 시퀀스로 변경해 줘야 한다. 일단 파이썬 시퀀스 형태만 만족하면 tuple, list형은 문제가 되지 않는다. numpy.ndarray 형도 numpy가 알아서 처리한다. 우선 blitz::TinyVector<>를 고려하지 않고 기본적인 인터페이스 파일을 만들어 보자.

// util.i

%module util

%{
#include "util.hh"
%}

%include "util.hh"

이렇게만 해도 컴파일이 되고, 파이썬 모듈을 읽어들일 수 있다. 하지만 함수를 실행할 수는 없다.

>>> import util
>>> a=[1,2,3]
>>> util.accumulate(a)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: in method 'accumulate', argument 1 of type 'util::Vector const &'

역시나 함수 인자의 형태가 맞지 않다는 에러를 내보낼 뿐이다. 즉, 파이썬의 시퀀스 형을 util::Vector 형으로 변환하는 방식을 알려줘야 한다. 문제는 SWIG가 blitz::TinyVector<>를 자동으로 처리해 주지 않는다는 점이다. 이건 %typemap을 써서 직접 변환 규칙을 만들어 줘야 한다.

%typemap(in) util::Vector3&
{
  $1 = convert_to_vector3($input);
  if (!$1) return NULL;
}

%typemap(argout) util::Vector3&
{
  for (int i = 0; i < 3; i++) {
    PyObject *o = PyFloat_FromDouble((*$1)[i]);
    PySequence_SetItem($input, i, o);
  }
}

%typemap(freearg) util::Vector3&
{
  if ($1) delete $1;
}

기본적인 역할은 %typemap(in)이 맡는다. 함수 인자의 util::Vector& 부분을 찾아 이 규칙을 적용한다. $1은 C++ 함수의 첫 번째 인자를 뜻하고, $input은 파이썬 측 함수의 인자를 가리킨다. convert_to_vector에 실제 변환 코드가 들어가는데, 이 코드는 아래에서 다룬다. util::weight 함수가 제대로 작동하려면 %typemap(argout)을 써야 한다. %typemap(argout)은 출력 역할을 맡은 인자를 처리할 수 있도록 해준다. 마지막으로 %typemap(freearg)는 형 변환 중에 힙에 생성한 자료 영역을 지울 수 있도록 해준다. 아래의 convert_to_vector를 보면 %typemap(freearg)를 사용한 이유를 알 수 있다.

%{
static util::Vector3*
convert_to_vector3(PyObject* input)
{
  if (!PySequence_Check(input)) {
    PyErr_SetString(PyExc_ValueError, "Expected a sequence");
      return NULL;
  }
  if (PySequence_Length(input) != 3) {
    PyErr_SetString(PyExc_ValueError, "Size mismatch. Expected 3 elements");
    return NULL;
  }

  util::Vector* v = new util::Vector3();
  for (int i = 0; i < 3; i++) {
    PyObject *o = PySequence_GetItem(input, i);
    if (PyNumber_Check(o))
      (*v)[i] = PyFloat_AsDouble(o);
    else {
      PyErr_SetString(PyExc_ValueError, "Sequence elements must be number");
      return NULL;
    }
  }
  return v;
}
%}

파이썬측의 입력이 시퀀스인지, 길이가 3인지를 확인한 후에 util::Vector를 임시로 생성하고 여기에 인자를 하나씩 복사하도록 했다. 복사하기 전에 각 원소가 숫자인지 확인하도록 했다.

아직 조금 더 남아있다. 형식 변환 규칙이 하나 더 필요한데, 바로 %typemap(typecheck)이다. util::weight 함수의 두 번째 인자가 기본값을 가지고 있는데, SWIG은 이럴 경우, util::weight(util::Vector3&)와 util::weight(util::Vector3&, double)을 각각 만들고, 동적 자료형 검사를 통해서 어느 함수를 사용할지 결정하도록 코드를 생성한다. 그러므로 다음과 같이 간단한 자료형 비교 코드도 만들어 줘야 한다.

%typemap(typecheck, precedence=SWIG_TYPECHECK_POINTER) util::Vector3&
{
  $1 = false;
  if (PySequence_Check($input) && PySequence_Length($input) == 3)
    $1 = true;
}

이제 distutils를 이용해서 컴파일하려면, 다음과 같은 setup.py를 만들고, python setup.py build_ext -i라고 입력하면 파이썬에서 불러올 수 있는 모듈을 만들어준다.

from distutils.core import setup, Extension

util_module = Extension(name = '_util',
                        sources = ['util.i', 'util.cc'],
                        depends = ['util.hh'],
                        swig_opts = ['-c++'])

setup(ext_modules = [util_module])

이제 파이썬 해석기를 다시 시작하고 모듈을 읽어오면(reload로는 안 된다.) 잘 작동하는 것을 확인할 수 있다.

>>> import numpy, util
>>> v1 = [1,2,3]
>>> v2 = numpy.array(v1, float)
>>> util.accumulate(v1)
6.0
>>> util.weight(v1, 2)
>>> v1
[2.0, 4.0, 6.0]
>>> util.weight(v2, 3)
>>> v2
array([ 3.,  6.,  9.])
>>> 

C++의 반복자(iterator)를 이용하여 파이썬 생성자(generator) 만들기

최근에 GMES의 C++ 모듈 설계를 변경하면서, 내부 자료형으로 STL의 컨네이너를 사용하는 클래스가 필요하게 되었다. 이 글에서 설명하려는 내용과 관계 없는 부분을 제거하면, 대략적인 형태는 아래와 같다.

// itergen.hh
#include <vector>

class A
{
public:
  A(): a(3)
  {
    for (size_t i = 0; i < a.size(); i++)
      a[i] = i;
  }

  std::vector<int>::const_iterator
  begin() const
  {
    return a.begin();
  }

  std::vector<int>::const_iterator
  end() const
  {
    return a.end();
  }

private:
  std::vector<int> a;
};

// itergen.cc
#include "itergen.hh"

이 코드는 파이썬으로 래핑하여 사용하게 되는데, 아래와 같이 생성자로 이용하고자 한다. 어떻게 해야 할까?

>>> for i in A():
...    print i
1
2
3

파이썬 생성자에 관한 API를 살펴봤는데, 이거 잠깐 보고 구현할 난이도는 아닌 것 같다. 이왕 SWIG를 이용한다면, 뭔가 더 간단한 방법이 없을까? Reference Man의 블로그에 방법이 나와 있다. 이 블로그에 나와 있는 내용을 간단하게 정리해 보자.

일단, C++ 클래스에선 시작과 끝 부분의 반복자를 반환하는 begin, end 메서드만 제공하면 된다. 이는 일반적으로 STL에서 사용하는 방식이므로 C++ 클래스에서 제공하는 것이 자연스럽다. C++ 코드에서 할 일은 여기까지다. SWIG가 제공하는 %extend 키워드를 사용하면 C++ 코드를 더럽히지 않고도 파이썬 래핑에 필요한 코드를 추가할 수 있다. swig의 인터페이스 파일은 아래와 같다.

// itergen.i
%module itergen

%{
#include "itergen.hh"
%}

%extend A {
  std::vector<int>::const_iterator* _begin()
  {
    return new std::vector<int>::const_iterator($self->begin());
  }

  std::vector<int>::const_iterator* _end()
  {
    return new std::vector<int>::const_iterator($self->end());
  }

  int _deref(std::vector<int>::const_iterator* it)
  {
    return **it;
  }

  bool _compref(std::vector<int>::const_iterator* lhs,
  std::vector<int>::const_iterator* rhs)
  {
    return *lhs == *rhs;
  }

  void _incref(std::vector<int>::const_iterator* it)
  {
    ++(*it);
  }

  %pythoncode %{
    def __iter__(self):
      it = self._begin()
      while not self._compref(it, self._end()):
          yield self._deref(it)
          self._incref(it)
  %}
};

%include "itergen.hh"

파이썬은 STL의 반복자를 제대로 다루지 못하므로, 이를 처리해 주는 모든 메서드를 C++ 측에서 만들어줘야 한다. C++ 코드를 수정할 것 없이 인터페이스 파일에서 %extend 키워드를 이용해서 관련 메서드를 추가한다.

파이썬 클래스는 __iter__ 메서드를 선언해주면, 생성자로 이용할 수 있다. 이 메서드는 앞에서 만들어준 C++ 메서드를 이용하여 C++의 반복자를 처리한다. 원리는 간단하다. STL에서 흔히 사용하는, 반복자를 이용한 컨테이너 순회를 이용하는 데, 차이점이라면 값의 반환으로 return 대신에 yield를 사용한다는 점이다. yield를 사용하면, 함수의 값은 반환되는 데, 함수는 종료되지 않고 next 메서드가 호출될 때마다 yield를 호출한 다음 지점부터 연속해서 실행한다. __iter__ 메서드는 %pythoncode 키워드를 이용하여 파이썬 프록시 클래스에 추가한다.

컴파일 과정은 아래와 같다. 인터페이스 파일로부터 프록시 클래스와 파이썬 인터페이스를 생성하고, C++ 코드를 컴파일 한 후, 공유 라이브러리로 만들어 주면 된다.

$ swig -c++ -python itergen.i
$ g++ -Wall -g -I/usr/include/python2.6 -c itergen_wrap.cxx
$ g++ -shared itergen.o itergen_wrap.o -o _itergen.so

이제 파이썬 해석기에서 모듈을 올리고 생성자를 사용해보자.

$ python
Python 2.6.5 (r265:79063, Apr 16 2010, 13:09:56)
[GCC 4.4.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import itergen
>>> a=itergen.A()
>>> iter(a)
<generator object __iter__ at 0xb7720234>
>>> for i in a: print i
...
0
1
2
>>>

완료!