読者です 読者をやめる 読者になる 読者になる

ローグウェーブソフトウェアのブログ

開発をシンプルに 安全で高品質のコードを 素早くお客様のもとへ

PyNLを使ってIMSL C関数をPythonクラスでラッピングする方法 - CodeBuzzから

IMSL リリース 製品とサービス 記事紹介

IMSLのプロダクトマネージャー Jennifer Lockeが書いたPyNLの記事をご紹介します。連立方程式やコレスキー係数などを題材に使っています。

roguewave.jp

IMSLはローグウェーブ数値計算統計解析ライブラリで、長い間幅広い業界で使われてきました。

blog.roguewave.jp

チュートリアルサポート環境Sphinxで書かれたドキュメンテーションをご覧ください。user forumでも情報を更新していく予定です。

ブログ記事: PyNLを使ってIMSL C関数をPythonクラスでラッピングする方法

blog.klocwork.com

IMSLをお使いの開発者はIMSL Cの関数にPythonからアクセスできるようになりました。IMSL Python数値ライブラリ (PyNL) 1.0がリリースされたからです。IMSL Cの関数の一部にはあらかじめPyNLに用意されているPythonラッパーが使えます。また、追加ラッパーを書いて必要に応じて他の関数へPythonでアクセスすることもできます。このブログ投稿ではPyNLでIMSL C関数のラッパーを作成する際の主要な手順をご紹介します。

PyNLはIMSL Cの関数と統合しやすいようなアーキテクチャで構築されているため、すなおな手続きで新しいラッパーを追加することができます。PyNLにはNumPyとSciPyが同梱されており、Python 2.7.8とPython 3.4.2に対応しています。ローグウェーブのウェブサイトからから依存関係や使用環境についてご確認ください。

既存のPyNLラッパーを活用する

IMSL Cの500を超えるAPIの多くは共通の設計特性を持っているため、IMSL Cの新しいラッパーを作成するには、もしあれば似た関数のラッパーを活用するのがおすすめです。多くの種類のAPIに対してリファレンス実装となるような代表的なアルゴリズムのセットが現在のPyNLには含まれています。例えば lin_sol_posdef() のラッパーを新しく作成したい場合など、たいていの場合は似た関数(この場合は lin_sol_gen() )のラッパーを再利用できます。既存のラッパーはLUという名のPythonクラスです。それでは先に進みましょう。(lin_sol_posdef等の関数の日本語による説明は、[http://roguewave.jp/help-support/documentation#c:title=日本語ドキュメントのIMSL Cライブラリ v8.5 ユーザーズガイド : Math 第一章をご覧ください)

ラッパーの開発

チュートリアルのホワイトペーパー How to access IMSL C functions from Python (PythonからIMSL Cの関数にアクセスする方法) でIMSL C関数をPython関数でラッピングする方法を記述しました。たいていの場合はPython関数で十分ですが、時にはPythonクラスを使ったほうが以下に示す通りPythonオブジェクトの状態を管理する場合などに便利です。参考のためにこのブログの末尾に完全なラッパーコードを付けました。

PyNLをインストールした後は、インストールフォルダ内のsite-package以下、imslパッケージにある全てのラッパーコードが使用可能になります。このソースコードを使って追加のラッパーを書いてPyNLに追加することもできます。

設定

Python site-packages内のPyNLエリアに新しいパッケージを作成すれば、ユーザー作成の複数のPyNLラッパーを共存させることができるためおすすめです。以下の手順にしたがって新しいパッケージを作成してください。ここでPYNL_AREAPythonのインストール先のsite-packages/imslフォルダを参照しています。

  1. PYNL_AREA/user_libという新しいフォルダを作成してください。
  2. PYNL_AREA/user_lib内にinit.pyというファイルを作成し、以下のコードを記述してください。
"""Initialize imsl.user_lib package.""" 
from imsl.user_lib.cholesky import Cholesky 
__all__ = ('Cholesky',) 
  1. パッケージ名user_libをPYNL_AREA/init.pyに追加してください。
  2. ラッピングしたいIMSL C関数をPYNL_AREA/_imsllib.pyに明記してください。例えばIMSL C関数のimsl_d_lin_sol_posdef()を記載するには、すでに記載されている他の関数の後に以下の1行を加えてください。
 self.expose("imsl_d_lin_sol_posdef ", None, _err_check, "math") 

PYNL_AREA/_imsllib.pyに従えばどこにこの行を追加すればよいのかわかります。この指定で、IMSL Cでエラーが検出された時にPythonの例外を投げるようPyNLの例外ハンドリングを設定できました。

注意: もしC/Stat関数をラッピングするのであれば、self.exposeの最後の引数はmathではなくstatになります。

ラッパーのサンプル

先述した通り lin_sol_posdef()のIMSL C APIlin_sol_gen()インターフェイスに類似しているため、PyNL LUクラスの設計とソースコード lin_sol_posdef()のPyNLクラスを作成するのに利用することができます。なお、既存のPyNL LUクラスはlin_sol_gen()の実数版と複素数版の両方をサポートしていますが、今回は便宜上実数版である imsl_d_lin_sol_posdef()のみをラップします。

PyNLのLUクラスには以下の特徴があります。

  • LU(a)クラス
    • Parameters
      • a: An NxN array containing the input matrix
    • Methods
      • cond() Compute the L1 norm condition number of the matrix
      • factor() Compute the pivoted LU factorization of the matrix
      • factor_full() Compute the pivoted LU factorization of the matrix
      • inv() Compute the inverse of the matrix
      • solve(b[, transpose]) Solve a system of linear equations using right-hand side b

こうした特徴はIMSL C関数の imsl_d_lin_sol_posdef()でも同様でなければなりません。PyNL LUクラスを利用して、以下の特徴を持った imsl_d_lin_sol_posdef()用のPythonクラスを作成することができます。

  • class Cholesky(a)
    • Parameters
      • a: An NxN array containing the input matrix
    • Methods
      • cond() Compute the L1 norm condition number of the matrix
      • factor() Compute the Cholesky factorization of the matrix
      • inv() Compute the inverse of the matrix
      • solve(b) Solve a system of linear equations using right-hand side b

このクラスは多くのメソッドや入力配列aの形などがLUクラスと同じであり、このLUクラスのソースコード(site-packages/imsl/linalg フォルダ内)を出発点とすることができます。なお、LU.factor_full() はLUアルゴリズムに特有のものであり、Choleskyクラスには不要です。

新しいクラスの設計を定義し終わったら、次のステップは、

  • 入力データの型をチェックする
  • IMSL C関数を呼び出すためにデータを変換する
  • 引数を設定してIMSL C関数を呼び出す

以下のコードスニペットを使ってこれらがどのように実現されるのか見ていきましょう。

入力データの型をチェックする

IMSL CはCライブラリなので扱うデータは一定の型を持ちます。そのためPythonラッパーは入力データが必要な精度を持つどの型を持つのかを決める必要があります。以下のコードはCholesky.initからの引用です。

   # attempt to promote matrix a to a compatible type.
        common_type = _numpy.promote_types(_numpy.float64, self._a.dtype)
        self._a = _numpy.asarray(self._a, dtype=common_type)

        if (not _numpy.issubdtype(self._a.dtype, _numpy.float64)):
            raise ValueError("array type {} not supported".format(
                             self._a.dtype.name))

IMSL C関数を呼び出すためにデータを変換する

ここでも同様に、IMSL Cは特定のデータ型を必要とするため、IMSL C関数に与えられるデータはNumPyのasarrayメソッドに適した型に変換されなければなりません。たとえばラッパーのCholesky.solveメソッドでは以下のようになります。

  b = _numpy.asarray(b, dtype=self._a.dtype, order='C')

引数を設定してIMSL C関数を呼び出す

Pythonラッパーは引数のリストを作成してIMSL C関数に渡します。この引数リストは求められるアクションに必要な必須とオプション引数(可変数引数)に対応します。なお、IMSL Cは任意引数がゼロで終端します。Cholesky.invの実装は

args = []

row_dim = self._a.shape[0]

args.append(row_dim)
args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
args.append(0)  # b is ignored

self._inverse = _numpy.empty([row_dim, row_dim],
                             dtype=self._a.dtype)

args.append(_constants.IMSL_INVERSE_ONLY)
args.append(_constants.IMSL_INVERSE_USER)
args.append(self._inverse.ctypes.data_as(_ctypes.c_void_p))

args.append(0)

func = _imsllib.imsllib.imsl_d_lin_sol_posdef
func(*args)

サンプルの使い方

IMSL C関数 imsl_d_lin_sol_posdef() のPyhonインターフェイスを見てみましょう。この関数は対称正定値行列の三元線形連立方程式を解くためのものです。またこの係数行列の逆行列、L1ノルムの条件数、コレスキー係数を計算する方法もお見せします。

{ \displaystyle
x_1-3x_2+2x_3=27 \\
-3x_1+10x_2-5x_3=-78 \\
2x_1-5x_2+6x_3=64
}

に対して、

""" Cholesky class example. """
import imsl.user_lib.cholesky as chol
import numpy as np

a = [[1.0, -3.0, 2.0], [-3.0, 10.0, -5.0], [2.0, -5.0, 6.0]]
b = [27.0, -78.0, 64.0]
chol = chol.Cholesky(a)
x = chol.solve(b)
print("Solution:")
print(x)
print("")

inv = chol.inv()
print("Inverse:")
print(inv)
print("")

# Compute a*inv(a) (=I)
prod = np.dot(a, inv)
print("a*inv(a):")
print(prod)
print("")

print("L1 condition number:")
print(chol.cond())
print("")

L = chol.factor()
print("Cholesky factors:")
print(L)
print("")

# Compute a-L*trans(L) (=0)
for i in range(x.size):
    L[i, i+1:] = 0

print("a - L*trans(L):")
print(a - np.dot(L, np.transpose(L)))

出力

Solution:
[ 1. -4.  7.]

Inverse:
[[ 35.   8.  -5.]
 [  8.   2.  -1.]
 [ -5.  -1.   1.]]

a*inv(a):
[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]

L1 condition number:
673.838709677

Cholesky factors:
[[ 1. -3.  2.]
 [-3.  1.  1.]
 [ 2.  1.  1.]]

a - L*trans(L):
[[ 0.  0.  0.]
 [ 0.  0.  0.]
 [ 0.  0.  0.]]

imsl_d_lin_sol_posdef()のラッパーコード完全版

IMSL C関数 imsl_d_lin_sol_posdef()PythonのCholeskyという名のクラスでラップしたコードは以下のとおりです。

"""Cholesky factorization related class, methods."""
import ctypes as _ctypes
import numpy as _numpy

import imsl._constants as _constants
import imsl._imsllib as _imsllib


class Cholesky():

    """Solve a real symmetric positive definite system of linear equations."""

    def __init__(self, a):
        """Instantiate Cholesky class."""
        self._a = a
        self._factor_factor = None
        self._cond = None
        self._inverse = None

        if self._a is None:
            raise TypeError("None not supported")

        self._a = _numpy.array(a, order='C')

        # attempt to promote matrix a to a compatible type.
        common_type = _numpy.promote_types(_numpy.float64, self._a.dtype)
        self._a = _numpy.asarray(self._a, dtype=common_type)

        if (not _numpy.issubdtype(self._a.dtype, _numpy.float64)):
            raise ValueError("array type {} not supported".format(
                             self._a.dtype.name))

        if self._a.ndim != 2:
            raise ValueError("array of dimension {} not"
                             " supported".format(self._a.ndim))
        if self._a.size == 0:
            raise ValueError("empty array not supported")
        if (self._a.shape[0] != self._a.shape[1]):
            raise ValueError("input matrix must be square")

    def solve(self, b):
        r"""Solve a real symmetric positive definite system of linear equations 
            using right-hand side `b`."""
        if b is None:
            raise TypeError("None not supported")

        b = _numpy.asarray(b, dtype=self._a.dtype, order='C')

        if b.ndim != 1:
            raise ValueError("array of dimension {} not"
                             " supported".format(b.ndim))

        if b.size != self._a.shape[0]:
            raise ValueError("dependent variable values length ({}) does not"
                             " match coefficient matrix row count"
                             " ({})".format(b.size, self._a.shape[0]))

        args = []

        row_dim = self._a.shape[0]

        args.append(row_dim)
        args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
        args.append(b.ctypes.data_as(_ctypes.c_void_p))

        solution = _numpy.empty_like(b)
        args.append(_constants.IMSL_RETURN_USER)
        args.append(solution.ctypes.data_as(_ctypes.c_void_p))

        if self._factor_factor is None:
            self._factor_factor = _numpy.empty([row_dim, row_dim],
                                               dtype=self._a.dtype)
        else:
            args.append(_constants.IMSL_SOLVE_ONLY)

        args.append(_constants.IMSL_FACTOR_USER)
        args.append(self._factor_factor.ctypes.data_as(_ctypes.c_void_p))

        if self._cond is None:
            args.append(_constants.IMSL_CONDITION)
            cond = _ctypes.c_double()
            args.append(_ctypes.byref(cond))

        args.append(0)

        func = _imsllib.imsllib.imsl_d_lin_sol_posdef
        func(*args)

        if self._cond is None:
            self._cond = cond.value

        return solution

    def factor(self):
        """Compute the Cholesky factorization of the matrix."""
        if (self._factor_factor is None or self._cond is None):
            args = []

            row_dim = self._a.shape[0]

            args.append(row_dim)
            args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
            args.append(0)  # b is ignored

            self._factor_factor = _numpy.empty([row_dim, row_dim],
                                               dtype=self._a.dtype)

            args.append(_constants.IMSL_FACTOR_ONLY)
            args.append(_constants.IMSL_FACTOR_USER)
            args.append(self._factor_factor.ctypes.data_as(_ctypes.c_void_p))

            args.append(_constants.IMSL_CONDITION)
            cond = _ctypes.c_double()
            args.append(_ctypes.byref(cond))

            args.append(0)

            func = _imsllib.imsllib.imsl_d_lin_sol_posdef
            func(*args)

            self._cond = cond.value

        return (_numpy.copy(self._factor_factor))


    def inv(self):
        """Compute the inverse of the matrix."""
        if self._inverse is None:
            args = []

            row_dim = self._a.shape[0]

            args.append(row_dim)
            args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
            args.append(0)  # b is ignored

            self._inverse = _numpy.empty([row_dim, row_dim],
                                         dtype=self._a.dtype)

            args.append(_constants.IMSL_INVERSE_ONLY)
            args.append(_constants.IMSL_INVERSE_USER)
            args.append(self._inverse.ctypes.data_as(_ctypes.c_void_p))

            args.append(0)

            func = _imsllib.imsllib.imsl_d_lin_sol_posdef
            func(*args)

        return _numpy.copy(self._inverse)

    def cond(self):
        """Compute the L_1 norm condition number of the matrix."""
        if self._cond is None:
            self.factor()

        return self._cond

編集後記

PythonからIMSL Cにアクセスする方法の解説、いかがだったでしょうか。データ解析を行いたいが今更Cのコードを書くのはちょっと、と二の足を踏んでいるデータサイエンティストの皆様は、ぜひこれを機会にIMSLの世界に手軽にアクセスしてみてください。

blog.roguewave.jp

日本オフィスへのお問い合わせはこちらカスタマーサポートフォームはこちらです。お気軽にお問い合わせください。

ローグウェーブ セールスエンジニア 柄澤