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

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

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

チュートリアル: Python から IMSL ライブラリへアクセスする方法

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

先日リリースされたIMSL CのPythonラッパー(PyNL)を使うためのチュートリアルが公開されましたので、概要をご紹介します。IMSLの経験的分位点分析機能を例にとって解説していきます。

blog.roguewave.jp

roguewave.jp

技術チュートリアル: Python から IMSL ライブラリへアクセスする方法

Tech tutorial: How to access IMSL C functions from Python

このチュートリアルはIMSL Cの統計ライブラリ(Stat)のimsls_d_empirical_quantiles()にもとづいています。PyNLにはラッパー自体の他にもわかりやすいエラーハンドリングやSphinxで書かれたドキュメントが付属しており(Sphinx日本ユーザ会)、ラッパーを独自開発していくためのよいガイドとなっています。

初期設定

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

  1. PYNL_AREA/user_libという新しいフォルダを作成してください。
  2. PYNL_AREA/user_lib内にinit.pyというファイルを作成してください。既存のファイル PYNL_AREA/cluster/init.py が参考になるでしょう。
  3. user_libというパッケージ名をPYNL_AREA/init.pyに追加してください。
  4. ラッピングしたいIMSL C関数をPYNL_AREA/_imsllib.pyに明記してください。例えばIMSL C関数の imsls_d_empirical_quantiles()を記載するには、すでに記載されている他の関数の後に以下の1行を加えてください。
self.expose("imsls_d_empirical_quantiles", None, _err_check, "stat")

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

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

ラッパーのサンプル

以下はIMSL Cの関数をPyNLにフレームワークとしてラップするための3つの段階です。

  1. データ変換を行う
  2. IMSL Cに渡す引数を設定する
  3. IMSL C関数からの返り値をパッケージ化する

このサンプルではこの3つのステップに従ってIMSL C/Statの関数imsls_d_empirical_quantiles()をラップしています。 empirical_quantiles.py ファイルを PYNL_AREA/user_libに置いてください。

"""Empirical Quantiles related functions."""
import ctypes as _ctypes
import numpy as _numpy
import imsl._constants as _constants
import imsl._imsllib as _imsllib
import collections as _collections

def empirical_quantiles(x, qprop):
 r"""Compute empirical quantiles."""
 ref_type = _numpy.float64
 # Convert Data
 _x = _numpy.asarray(x, order='C', dtype=ref_type)
 if _x.ndim != 1:
 raise ValueError("x must be a 1-D array")
 _qprop = _numpy.asarray(qprop, order='C', dtype=ref_type)
 if _qprop.ndim != 1:
 raise ValueError("qprop must be a 1-D array")
 _n_obs = _x.shape[0]
 _n_qprop = _qprop.shape[0]
 # Set up arguments to be passed to the CNL function.
 args = []
 # Required input argument list
 args.append(_ctypes.c_int(_n_obs))
 args.append(_x.ctypes.data_as(_ctypes.c_void_p))
 args.append(_ctypes.c_int(_n_qprop))
 args.append(_qprop.ctypes.data_as(_ctypes.c_void_p))
 # Create space for results
 _result = _numpy.empty(_n_qprop, dtype=ref_type)
 args.append(_constants.IMSLS_RETURN_USER)
 args.append(_result.ctypes.data_as(_ctypes.c_void_p))
 _n_missing = _ctypes.c_int32()
 args.append(_constants.IMSLS_N_MISSING)
 args.append(_ctypes.byref(_n_missing))
 _xlo = _numpy.empty(_n_qprop, dtype=ref_type)
 args.append(_constants.IMSLS_XLO_USER)
 args.append(_xlo.ctypes.data_as(_ctypes.c_void_p))
 _xhi = _numpy.empty(_n_qprop, dtype=ref_type)
 args.append(_constants.IMSLS_XHI_USER)
 args.append(_xhi.ctypes.data_as(_ctypes.c_void_p))
 args.append(0) # Terminating zero for CNL argument list.
 # Call the CNL function
 _imsllib.imsllib.imsls_d_empirical_quantiles(*args)
 # Package results from the CNL function
 result = _collections.namedtuple("empirical_quantiles",
 ["quantiles", "n_missing",
 "xlo", "xhi"])
 result.quantiles = _result
 result.n_missing = _n_missing
 result.xlo = _xlo
 result.xhi = _xhi
 return result

サンプルの使い方

このサンプルでは5つの経験的分位点(empirical quantiles)をサイズ30のサンプルから取り出しました。0.5の分位点がサンプルの中央値(median)に相当します。使用したデータはHinkley (1977)*1 と Velleman and Hoaglin (1981)*2 に基づいており、ミネアポリス・セントポール都市圏の3月の降水量を30年にわたって計測したものです。

import imsl.user_lib as user_lib
x = [0.77, 1.74, 0.81, 1.20, 1.95,
 1.20, 0.47, 1.43, 3.37, 2.20,
 3.00, 3.09, 1.51, 2.10, 0.52,
 1.62, 1.31, 0.32, 0.59, 0.81,
 2.81, 1.87, 1.18, 1.35, 4.75,
 2.48, 0.96, 1.89, 0.90, 2.05]

qprop = [0.01, 0.5, 0.9, 0.95, 0.99]
result = user_lib.empirical_quantiles(x, qprop)
print(" Smaller Empirical Larger")
print("Quantile Datum Quantile Datum")
for i in range(5):
 print("{:6.2f} {:9.2f} {:9.2f} {:9.2f}".format(
 qprop[i], result.xlo[i], result.quantiles[i],
 result.xhi[i]))

出力結果

          Smaller  Empirical Larger
Quantile  Datum    Quantile  Datum
  0.01     0.32      0.32     0.32
  0.50     1.43      1.47     1.51
  0.90     3.00      3.08     3.09
  0.95     3.37      3.99     4.75
  0.99     4.75      4.75     4.75

編集後記

いかがでしょうか。これは先にご紹介した記事の内容を、別の関数について改めて説明したものです。あらためてこれを機会にIMSLの世界に手軽にアクセスしてみてください。

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

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

*1: Hinkley, David (1977), On quick choice of power transformation, Applied Statistics, 26, 67-69.

*2:Velleman, Paul F. and David C. Hoaglin (1981), Applications, Basics, and Computing of Exploratory Data Analysis, Duxbury Press, Boston