matplotlib.pyplot.contourを上手に使うことで、等値線の座標(x,y)を取得することができます。
本記事では、解析的な扱いが難しい2変数関数f(x,y)の等値線を表す方程式
$$ f(x,y)=\mathrm{c} $$
の解(x,y)を
- 数行のコードで
- 広く用いられるライブラリベースで
簡単に得る方法について解説します。
matplotlib.pyplot.contourの補間
matplotlib.pyplot.contourはmatplotlibで等値線を描くメソッドです。
デフォルトでは「marching squares法」と呼ばれる補間手法を用いて、等値線を滑らかに描画します。
この補間手法は等値線(とその頂点)の位置や向きをたくさんのパターンで考え、
それぞれの点ごとに最も似たパターンを探す方法です。
自力で実装するのは大変そうですが、matplotlibではcontourが内部で計算してくれます。
実例：「Rosenbrock関数の等値線の座標を取得」
では、最適化の分野では有名な「Rosenbrock関数」を例に、等値線の座標を取得します。
Rosenbrock関数とは
Rosenbrock関数は非線形最適化アルゴリズムの性能評価に用いられるベンチマーク関数です。
湾曲した谷がある点が特徴です。
2変数の場合は以下の式で表されます。
$$ f(x,y) = (1-x)^2 + 100(y-x^2)^2 $$
Pythonの関数として、以下のように定義します。
def rosenbrock(x, y):
return (1 - x)**2 + 100*(y - x**2)**2試しに、plt.contourで等値線を描画
始めに、単純なplt.contour、plt.contourfで等値線プロットを描いてみます。
contourfは対数スケールで、cmap=terrainで描画します(cmapについてはこちら)。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
def rosenbrock(x, y):
return (1-x)**2 + 100*(y-x**2)**2
# 格子
x = np.linspace(-2, 2, 500)
y = np.linspace(-1, 3, 500)
X, Y = np.meshgrid(x, y)
Z = rosenbrock(X, Y)
# 等値線（線形）
levels_line = np.arange(0,1001,50)
# 塗り（対数）
levels_fill = np.logspace(0, 4, 40)
norm = colors.LogNorm(vmin=levels_fill[0], vmax=levels_fill[-1])
plt.figure(figsize=(7,6))
# 塗り：log
cf = plt.contourf(X, Y, Z,
levels=levels_fill,
cmap="terrain",
norm=norm,
extend="both")
# 線：linear
cs = plt.contour(X, Y, Z,
levels=levels_line,
colors="black",
linewidths=0.6)
plt.clabel(cs, fmt="%.0f")
plt.plot(1,1,"ro")
plt.gca().set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.show()
以下の図が出力されます。青色、緑色の領域が”谷”に当たります。
では本題の、特定の値の等値線の座標を取得してみましょう。
任意の値の等値線を取得(本題)
f(x,y)=100となる等値線は以下の赤線です。
const = 100
plt.title(f"f(x,y)={const}")
plt.contourf(X, Y, Z,
levels=levels_fill,
cmap="terrain",
norm=norm,
extend="both") # normの設定などは前節のコードと共通。
plt.contour(X, Y, Z, levels=[const],colors="r",linewidths=3)
plt.show()
まずはplt.contourを取得します。
ここからはグラフ画像は必要ないので、plt.close()を付けます。
cs = plt.contour(X, Y, Z, levels=[const],colors="r",linewidths=3)
plt.close()
paths = cs.get_paths()plt.contourの返り値「cs」のメソッドget_paths()を使うと、
「その値を取る等値線の座標(複数本ある場合は全て)のリスト」を得られます。
print(len(paths))
... 1今回は1本です(「2本あるのでは？」という疑問については、後ほど種明かしがあります)。
では、リストの要素を取り出します。
array = paths[0]
array
... Path(array([[-1.42473828, 3. ],
[-1.42284569, 2.99468769],
[-1.42189176, 2.99198397],
...,
[-1.98452853, 2.98396794],
[-1.98650879, 2.99198397],
[-1.9885053 , 3. ]], shape=(2351, 2)), array([1, 2, 2, ..., 2, 2, 2], shape=(2351,), dtype=uint8))リストの要素はmatplotlibで線分や曲線を表す基本的なクラス「matplotlib.path.Path」です。
shapeは(頂点の数、2)です。
したがって、今回は頂点の数が2351個あることが分かります。
matplotlib.path.Pathクラスには頂点の座標を返すメソッド「vertices」が定義されています。
これを用いて、等値線の頂点(通過点)の座標(x,y)を取得できます。
x = array.vertices[:,0]
y = array.vertices[:,1]x,yはそれぞれ2351次元のnumpy.ndarrayです。
最後に、x,yのグラフを描いて、本当に等値線が取れたか確かめてみましょう。
plt.plot(x,y,c="r")
plt.scatter(x[0],y[0],c="r",s=50)
plt.show()
赤い丸から谷に平行な線を描いています。
ただ、描画領域の端であるy=3.0の部分で、2本の等値線が繋がっています。
これが等値線の本数が2本でなく1本となった理由です。
まとめ
グラフ描画ライブラリであるmatplotlibのcontourを使うと、等値線の座標を取得することができます。
解析的には解けないデータの等値線を簡単に得られる点が
matplotlib.pyplot.contourを使うメリットですね。
以上、読んでいただきありがとうございました。
コメント