【Common Lisp】py4clでmatplotlib作図

Common Lispで数学的な計算結果をヒートマップで可視化しようとしたとき、vgplotではgnuplotの低水準APIを直接叩く必要があり、コードが複雑になった。
py4clを使うと、Common LispからPythonのサブプロセスを起動してmatplotlibを呼び出せ、imshowなどの高水準APIをそのまま利用できる。
python-execで文を実行し、python-evalで式の結果を受け取り、import-moduleでモジュールをLispパッケージとして扱うのが基本的な使い方になる。

1. なぜpy4clか

Common Lispで数学の問題を考えていたら、vgplotでヒートマップを作図しようとして詰まりました。

(defun compare-powers (x y)
  "x^y と y^x を比較して、x^y<y^x なら -1、等しければ 0、大きければ 1 を返す。"
  (cond ((< (* y (log x)) (* x (log y))) -1)
        ((> (* y (log x)) (* x (log y))) 1)
        (t 0)))

(defun sign-matrix (n)
  "y=1..n, x=1..n の順で -1,0,1 を並べた行列を返す。"
  (loop for y from 1 to n collect
        (loop for x from 1 to n collect
              (compare-powers x y))))Code language: Lisp (lisp)

そこで、Pythonのmatplotlibをpy4clを使って呼び出してみることにしました。

1.1. vgplotでのヒートマップ

Common Lispからgnuplotを扱うライブラリとして vgplot があります¹。
線グラフや散布図なら十分使えますが、「格子の各マスを色で塗りつぶす」ような図示では困りました。
gnuplotの with image や pm3d を vgplot:format-plot 経由で直接叩けば実現できますが、かなりコードが複雑になります²。

;;;; sign(x^y - y^x) を三色グラフとして描く

(ql:quickload :vgplot)

(defun write-sign-matrix-file (n pathname)
  "gnuplot の matrix 用データファイルを書き出す。"
  (with-open-file (out pathname
                       :direction :output
                       :if-exists :supersede
                       :if-does-not-exist :create)
    (dolist (row (sign-matrix n))
      (format out "~{~A~^ ~}~%" row))))

(defun plot-sign (n &optional (file "/tmp/sign-matrix.dat"))
  "1..n の格子点を、画像風の三色プロットで描く。"
  (write-sign-matrix-file n file)
  (vgplot:close-all-plots)

  (vgplot:format-plot t "unset key")
  (vgplot:format-plot t "set size ratio -1")
  (vgplot:format-plot t "set view map")
  (vgplot:format-plot t "set xrange [0.5:~A.5]" n)

  ;; 環境によって上下が逆に見える場合は次の行を使う:
  ;; (vgplot:format-plot t "set yrange [~A.5:0.5]" n)
  (vgplot:format-plot t "set yrange [0.5:~A.5]" n)

  (vgplot:format-plot t "set xtics 1,20,~A" n)
  (vgplot:format-plot t "set ytics 1,20,~A" n)
  (vgplot:format-plot t "set cbrange [-1:1]")

  ;; -1=red, 0=black, 1=blue
  (vgplot:format-plot
   t
   "set palette defined (-1 'red', 0 'black', 1 'blue')")

  ;; 見た目調整
  (vgplot:format-plot t "set tics out")
  (vgplot:format-plot t "set border lw 1")

  (vgplot:format-plot t "plot '~A' matrix with image" file)

  (vgplot:title (format nil "sign(x^y - y^x) for 1 <= x, y <= ~A" n))
  (vgplot:xlabel "x")
  (vgplot:ylabel "y"))

;; 実行例
(plot-sign 120)Code language: Lisp (lisp)

1.2. py4clとmatplotlibでのヒートマップ

一方、Pythonの matplotlib には pcolormesh や imshow といった高水準APIがあり、2次元配列を渡すだけで色付きのグリッドが描けます。

(eval-when (:compile-toplevel :load-toplevel :execute)
  (ql:quickload :py4cl)
  (setf py4cl:*python-command* "python3")
  (py4cl:import-module "matplotlib.pyplot" :as "plt")
  (py4cl:import-module "matplotlib.colors" :as "mcolors"))

(defun tick-list (n &optional (step 20))
  (let ((ticks (loop for k from 1 to n by (if (>= n step) step 1)
                     collect k)))
    (if (= (car (last ticks)) n)
        ticks
        (append ticks (list n)))))

(defun plot-sign/py (n)
  "py4cl 経由で matplotlib を使い、-1=red, 0=black, 1=blue で描画する。"
  (let* ((mat   (sign-matrix n))
         (ticks (coerce (tick-list n) 'vector))
         (cmap  (mcolors:listedcolormap '("red" "black" "blue")))
         (norm  (mcolors:boundarynorm #(-1.5 -0.5 0.5 1.5) 3)))
    (plt:figure :figsize #(8 8))
    (plt:imshow mat
                :origin "lower"
                :interpolation "nearest"
                :cmap cmap
                :norm norm
                :extent (vector 0.5 (+ n 0.5) 0.5 (+ n 0.5)))
    (plt:xticks ticks)
    (plt:yticks ticks)
    (plt:xlim 0.5 (+ n 0.5))
    (plt:ylim 0.5 (+ n 0.5))
    (plt:xlabel "x")
    (plt:ylabel "y")
    (plt:title (format nil "sign(x^y - y^x) for 1 <= x, y <= ~D" n))
    (plt:axis "equal")
    (plt:show)))
Code language: Lisp (lisp)

今回は「数式や計算はCommon Lispで、可視化だけPythonに任せる」と割り切って、py4clを使ってみます。

vgplot で描くなら
SBCL、Quicklisp、vgplot、gnuplot
py4cl で描くなら
SBCL、Quicklisp、py4cl、Python、numpy、matplotlib

なお、SBCL は必須ではなく、Quicklisp は SBCL 以外の複数の Common Lisp 処理系でも使えます。

とはいえ、言うほどpy4clでスッキリ書けるかというと、そうでもない気もする

2. py4clのロードとセットアップ

Quicklispが入っていれば、次でロードできます。

(ql:quickload :py4cl)Code language: Lisp (lisp)

asdf:load-system "py4cl" でも読み込めます³。

(py4cl:python-eval "1 + 2")
; => 3Code language: Lisp (lisp)

ここで END-OF-FILE エラーが出る場合、*python-command* が指す実行ファイルの問題やPython起動時の初期化エラーが考えられます。

環境によってはpythonの起動コマンドを明示しないと、正しいPythonを見つけられないことがあります。
デフォルトで使うコマンドは *python-command* 変数で設定します。

(setf py4cl:*python-command* "python3")Code language: Lisp (lisp)

ターミナルで python3 を単体で起動して正常か確認すると、切り分けが早いです。
フルパスで指定することもできます。

(setf py4cl:*python-command* "/usr/bin/python3")Code language: Lisp (lisp)

py4clは内部でPythonのサブプロセスを起動して通信します⁴。
python-start でサブプロセスを明示的に起動し、簡単な式を評価して確認します⁵。

2.1. numpyやmatplotlibのインストール

ただし、この場合、必要なPythonモジュールをインストールしておく必要があります。

python3 -m pip install -U numpy matplotlib

# 確認
python3 -c "import numpy, matplotlib; print(numpy.__version__, matplotlib.__version__)"Code language: CSS (css)

3. py4clの基本API

py4clの日常的な使い方は、ほぼ次の3つで完結します。

3.1. python-exec

Pythonの「文」を実行します。
statementと呼ばれる種類の構文で、代入や関数定義、import など、値を返さない操作に使います。
戻り値は NIL です。

(py4cl:python-exec "x = 10")Code language: Lisp (lisp)

3.2. python-eval

Pythonの「式」を評価し、結果をLispの値として返します。
expressionと呼ばれる種類の構文です⁶。

(py4cl:python-eval "x * 10")
; => 100Code language: Lisp (lisp)

x = 10 を python-eval に渡すとエラーになります。
代入は文なので python-exec 側の仕事です。
この使い分けだけ最初に意識しておくと詰まりにくくなります。

3.3. import-module

Pythonモジュールを読み込み、Lispパッケージとして扱えるようにします。
:as でパッケージ名のエイリアスを指定できます。

(py4cl:import-module "math" :as "math")
(math:sqrt 2)
; => 1.4142135623730951Code language: Lisp (lisp)

matplotlib.pyplot を plt として読むのが定番の書き方です。

(py4cl:import-module "matplotlib.pyplot" :as "plt")Code language: Lisp (lisp)

4. LispとPythonの型の変換

py4clはLispとPython間でデータを自動変換します⁷。

Pythonのリストは、Lispのベクタになります。

(py4cl:python-eval "[1, 2, 3]")
; => #(1 2 3)  ← Lispのベクタとして返ってくるCode language: Lisp (lisp)

LispのArrayはNumPy配列として渡せます。
これがヒートマップを描くときに効いてきます。

(py4cl:import-module "numpy" :as "np")

(let ((a (make-array '(2 3) :initial-contents '((1 2 3) (4 5 6)))))
  (np:sum a))
; => 21Code language: Lisp (lisp)

NumPyの計算結果は、Lispの数値やベクタとして返ってきます。
ただし、NumPy配列同士の演算をいったんLispの値として取り出してから再び渡すと、単なるリストとして扱われ型情報が失われることがあります。
これを避けるのが次に説明する remote-objects* の役割です⁸。

4.1. matplotlibで線グラフ

最小構成の例として、まず線グラフを描いてみます。

(ql:quickload :py4cl)
(setf py4cl:*python-command* "python3")

(py4cl:import-module "matplotlib.pyplot" :as "plt")

(let ((xs #(1 2 3 4))
      (ys #(1 4 9 16)))
  (plt:plot xs ys)
  (plt:xlabel "x")
  (plt:ylabel "y")
  (plt:title "sample plot")
  (plt:show))Code language: Lisp (lisp)

plt:plot の引数にLispのベクタを直接渡せます。
plt:show でウィンドウが開きます。
ファイルに保存したいなら plt:show の代わりに plt:savefig を使います⁹。

(plt:savefig "/tmp/plot.png")
(plt:savefig "/tmp/plot.svg")Code language: Lisp (lisp)

5. 格子の塗りつぶし

2次元配列を色付きのグリッドとして描きます。

5.1. 行列を作る

今回は x^y と y^x の大小を比べた結果を、-1・0・1で表した行列を作ります。
compare-powers は呼び出し元で定義済みとして話を進めます。

(defun sign-matrix (n)
  "n×n の行列を返す。imshow は左上原点なので、y方向を反転して格納する。"
  (let ((a (make-array (list n n) :element-type 'fixnum)))
    (loop for y from 1 to n do
      (loop for x from 1 to n do
        (setf (aref a (- n y) (1- x))
              (compare-powers x y))))
    a))Code language: Lisp (lisp)

5.2. imshowで描く

imshow は行列の値をカラーマップで色に変換して表示します。
値の範囲が -1・0・1 なので、0・1・2 にずらして3色のカラーマップに対応させます。

カラーマップの定義は python-exec でPython側に書きます。
ListedColormap はインデックスごとに色を対応させる離散的なカラーマップです¹⁰。

NumPy配列の演算は、Lispへの変換を挟まずPython側で完結させるため、remote-objects* の中で行います。

(ql:quickload :py4cl)
(setf py4cl:*python-command* "python3")

(py4cl:import-module "matplotlib.pyplot" :as "plt")
(py4cl:import-module "numpy" :as "np")

(py4cl:python-exec
 "from matplotlib.colors import ListedColormap
cmap = ListedColormap(['red', 'black', 'blue'])")

(defun plot-sign (n)
  (py4cl:remote-objects*
    (let ((z (np:array (sign-matrix n))))
      (plt:figure :figsize #(6 6))
      (plt:imshow (np:add z 1)
                  :origin "lower"
                  :interpolation "nearest"
                  :cmap (py4cl:python-eval "cmap")
                  :vmin 0
                  :vmax 2)
      (plt:xlabel "x")
      (plt:ylabel "y")
      (plt:title (format nil "sign(x^y - y^x), 1 <= x,y <= ~A" n))
      (plt:show))))

(plot-sign 120)Code language: Lisp (lisp)

いくつか補足します。

:origin "lower" を指定すると、y軸の原点が左下になります。
数学的な座標系と合わせたいときに使います¹¹。

:interpolation "nearest" はピクセルを補間せず、各マスをそのままの色で表示します。
指定しないと、隣接するセルの間でグラデーションがかかることがあります¹²。

np:add z 1 は z に1を加えて値の範囲を 0・1・2 にずらす操作です。
ListedColormap はインデックス0・1・2に赤・黒・青を割り当てているので、これで対応が取れます。

5.3. remote-objects* が必要な理由

(np:add z 1) の結果をいったんLispの値として受け取ってから imshow に渡すと、NumPyのndarrayではなくLispのリストとして渡ることがあります。

remote-objects* はPythonオブジェクトをLispに変換せずPython側のまま保持する仕組みで、NumPy演算の結果をそのままNumPy配列として次の関数に渡せます¹³。

6. まとめと注意点

*python-command* の設定は最初に一度やっておきます。
忘れると END-OF-FILE で落ちることがあります。

python-exec は文、python-eval は式です。
代入は exec、計算結果を受け取るのは eval と覚えておくと混乱しません。

NumPy配列どうしの演算を連鎖させるときは、remote-objects* の中でまとめて行います。
vgplotでは with image に降りる必要があった処理が、py4cl と imshow の組み合わせなら高水準APIで書けます¹⁴。

Common Lispの計算資産を活かしながら可視化だけPythonに任せる構成は、両方の強みを引き出せます。

vgplotはVolker Sarodnick氏が開発したCommon Lispライブラリで、gnuplotへのインターフェースとして機能します。APIの設計はOctaveやMATLABのplot系コマンドに似せる方針で作られています。 – volkers/vgplot – GitHub
vgplotのAPIリファレンスを見ると、pcolormesh・imshow・heatmapに相当する高水準関数は公開されていません。gnuplot側の with image や pm3d は format-plot 経由で呼ぶことができます。 – The vgplot Reference Manual
py4clはBen Dudson氏（ヨーク大学）が開発し、GitHubで公開されています。Quicklispからインストールするほか、リポジトリを ~/quicklisp/local-projects/ にクローンしてASDF経由でロードする方法も使えます。 – bendudson/py4cl – GitHub
py4clはストリームを使って別プロセスのPythonと通信する方式を採っています。これはcl4pyが取ったアプローチと同じで、CFFFIを使うburgled-batteriesとは異なる設計です。 – bendudson/py4cl – README
py4clの各関数は、Pythonプロセスが起動していなければ python-start-if-not-alive を自動で呼んで起動します。python-start と対になる python-stop でプロセスを明示的に終了することもできます。 – bendudson/py4cl – README
python-eval はより高水準のラッパー関数（import-module による関数呼び出し等）の内部でも使われています。直接 python-eval を呼ぶのは、式の評価結果をLispで受け取りたい場合の低水準な使い方です。 – bendudson/py4cl – README
py4clのデータ変換はテキスト経由で行われます。Python側では lispify 関数がLispリーダーで読める形式に変換し、Lisp側では pythonize 関数がPythonで評価できる文字列を出力します。なお、Pythonにはないlispの数値型（complex integerなど）は変換できない場合があります。 – bendudson/py4cl – README
py4clのSciPy/NumPy向けドキュメントにも、「NumPy配列の算術演算はremote-objects環境内で行うのがよい。1次元配列をLispに戻してからPythonに送ると単なるリストになる」と説明されています。 – py4cl/docs/scipy.org – GitHub
py4clのmatplotlib向けドキュメントでは、plt:plot、plt:xlabel、plt:savefig、plt:show を呼ぶ形の例が示されています。保存フォーマットはファイル名の拡張子（.png、.pdf等）で自動判定されます。 – py4cl/docs/matplotlib.org – GitHub
ListedColormap は matplotlib.colors モジュールに属するクラスで、色名またはRGB(A)値のリストからカラーマップオブジェクトを生成します。リストの要素数がそのままエントリ数になります。 – matplotlib.colors.ListedColormap — Matplotlib documentation
imshow のデフォルトは :origin "upper" で、行列のインデックス0が左上に配置されます。これは行列表示や画像処理の慣例に合わせたものです。:origin "lower" にすると数学的な座標系（y軸が上向き）と一致させることができます。 – origin and extent in imshow — Matplotlib documentation
matplotlibのデフォルトの補間方式は 'auto' で、アップサンプリング率が3倍以上の場合は 'nearest'、それ以下では 'hanning'（アンチエイリアス）が使われます。格子の各セルを明確に区別して表示したい場合は明示的に 'nearest' を指定します。 – Interpolations for imshow — Matplotlib documentation
remote-objects* は内部でpy4clのPythonオブジェクト管理機能を使い、オブジェクトへの参照をLisp側に返す代わりにPython側のハンドルとして保持します。remote-objects（アスタリスクなし）は最後の式の結果だけをLispに返すバージョンです。 – bendudson/py4cl – README
py4clはSBCL・CCL・ECLでの動作が確認されており、CLISPは uiop:launch-program が未対応のため使えません。また、ASDF3バージョン3.2.0（2017年1月）以降が必要です。 – bendudson/py4cl – README