【NumPy】格子点の作成

numpy.meshgrid()

numpy.meshgrid()

数値計算においては、平面上あるいは空間内の各点に定義された量を表すために 格子点（grid point）を扱う機会が多くなります。離散データだけでなく、本来ならば連続実数として扱うべき数値も格子点を使って近似します。たとえば実変数 $x, y$ によって定義された関数 $z = f (x, y)$ を可視化するためには、 $x y$ 平面上から有限個数の点を選んで 2 次元配列に格納し、各要素を関数に入れて値を計算する必要があります。このような目的のために、NumPy には格子点を生成する numpy.meshgrid() という関数が用意されています。

numpy.meshgrid() は一次元配列を受け取って格子点を生成する関数です。以下のコードを実行すると 3 × 3 の格子点を生成します。

# NUMPY_MESHGRID

# In[1]

import numpy as np

# x=[0 1 2]
x = np.arange(3)

# y=[0 1 2]
y = np.arange(3)

# 3×3の格子点を作成
X, Y = np.meshgrid(x, y)

print("配列X\n{}\n".format(X))
print("配列Y\n{}".format(Y))

# 配列X
# [[0 1 2]
#  [0 1 2]
#  [0 1 2]]

# 配列Y
# [[0 0 0]
#  [1 1 1]
#  [2 2 2]]

このように、2 次元格子点は 2 個の 2 次元配列を組合わせて表現されます。下図にあるように、それぞれの配列の同じインデックスの数値を参照することで特定の座標を参照できます。

次は numpy.linspace() を使って浮動小数点型要素をもつ 1 次元配列を生成し、numpy.meshgrid() に渡してみます。

# In[2]

# 配列の要素数
n = 9

x = np.linspace(0, 5, n)
y = np.linspace(0, 5, n)

# 5×5の格子点を作成
X, Y = np.meshgrid(x, y)

print("配列X\n{}\n".format(X))
print("配列Y\n{}".format(Y))

'''
配列X
[[0.    0.625 1.25  1.875 2.5   3.125 3.75  4.375 5.   ]
 [0.    0.625 1.25  1.875 2.5   3.125 3.75  4.375 5.   ]
 [0.    0.625 1.25  1.875 2.5   3.125 3.75  4.375 5.   ]
 [0.    0.625 1.25  1.875 2.5   3.125 3.75  4.375 5.   ]
 [0.    0.625 1.25  1.875 2.5   3.125 3.75  4.375 5.   ]
 [0.    0.625 1.25  1.875 2.5   3.125 3.75  4.375 5.   ]
 [0.    0.625 1.25  1.875 2.5   3.125 3.75  4.375 5.   ]
 [0.    0.625 1.25  1.875 2.5   3.125 3.75  4.375 5.   ]
 [0.    0.625 1.25  1.875 2.5   3.125 3.75  4.375 5.   ]]
　
配列Y
[[0.    0.    0.    0.    0.    0.    0.    0.    0.   ]
 [0.625 0.625 0.625 0.625 0.625 0.625 0.625 0.625 0.625]
 [1.25  1.25  1.25  1.25  1.25  1.25  1.25  1.25  1.25 ]
 [1.875 1.875 1.875 1.875 1.875 1.875 1.875 1.875 1.875]
 [2.5   2.5   2.5   2.5   2.5   2.5   2.5   2.5   2.5  ]
 [3.125 3.125 3.125 3.125 3.125 3.125 3.125 3.125 3.125]
 [3.75  3.75  3.75  3.75  3.75  3.75  3.75  3.75  3.75 ]
 [4.375 4.375 4.375 4.375 4.375 4.375 4.375 4.375 4.375]
 [5.    5.    5.    5.    5.    5.    5.    5.    5.   ]]
'''

細かい話になりますが、np.linspace() の第 3 引数に $2^{n} + 1$ 型の数値 ( $129$ や $259$ など) を渡すと、（整数値で下限と上限を指定された）範囲を綺麗な形で分割できます。2 次元格子点データを 2 変数で定義された関数 $z = f (x, y)$ に渡すと、各格子点における $z$ の値を計算して 2 次元配列に格納します。

# In[3]

# 2変数関数を定義
def f(x, y):
    return x**2 - y**2

# z=f(x,y)を計算
Z = f(X, Y)
Z_round = np.round(Z, 2)

print(Z_round)

'''
[[  0.     0.39   1.56   3.52   6.25   9.77  14.06  19.14  25.  ]
 [ -0.39   0.     1.17   3.12   5.86   9.38  13.67  18.75  24.61]
 [ -1.56  -1.17   0.     1.95   4.69   8.2   12.5   17.58  23.44]
 [ -3.52  -3.12  -1.95   0.     2.73   6.25  10.55  15.62  21.48]
 [ -6.25  -5.86  -4.69  -2.73   0.     3.52   7.81  12.89  18.75]
 [ -9.77  -9.38  -8.2   -6.25  -3.52   0.     4.3    9.38  15.23]
 [-14.06 -13.67 -12.5  -10.55  -7.81  -4.3    0.     5.08  10.94]
 [-19.14 -18.75 -17.58 -15.62 -12.89  -9.38  -5.08   0.     5.86]
 [-25.   -24.61 -23.44 -21.48 -18.75 -15.23 -10.94  -5.86   0.  ]]
'''

Matplotlib を使うと、X, Y, Z を様々な形式で表示することができます（詳細についてはこちらの記事を参照してください）。ここでは、mpl_toolkits.mplot3dモジュールから axes3d をインポートして、曲面 (surface) を描画してみます。

# In[4]

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# Figureと3DAxesを追加
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection="3d")

# 3DAxesの設定
ax.set_xlabel("x", size = 15)
ax.set_ylabel("y", size = 15)
ax.set_zlabel("z", size = 15)
ax.view_init(45, 45)

# 曲面を描画
ax.plot_surface(X, Y, Z, cmap="autumn")

plt.show()

2次元格子点

データを可視化するたびに格子点を作成していると手間がかかります。関数やクラスを定義することによって作業の効率化を図れるかもしれません。その一例として、Mesh2Dクラス を定義してみます。

# NUMPY_MESHGRID_2D

# In[1]

import numpy as np

# 2次元格子点クラス
class Mesh2D:
    def __init__(self, xrange, yrange, n):
        self.n = n
        self.x_min = xrange[0]
        self.x_max = xrange[1]
        self.y_min = yrange[0]
        self.y_max = yrange[1]
        self.x = np.linspace(xrange[0], xrange[1], n)
        self.y = np.linspace(yrange[0], yrange[1], n)
        self.X, self.Y = np.meshgrid(self.x, self.y)

Mesh2Dクラスのインスタンスは、x の範囲と y の範囲、縦横の要素数を指定して生成します。格子点はデータ属性 Mesh2D.X と Mesh2D.Y で参照することができます。Mesh2Dクラスを使って、 $z = \cos (x + y)$ を 3 次元グラフに描いてみましょう。

# In[2]

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# 2変数関数を定義
def f(x, y):
    return np.cos(x + y)

# Mesh2Dクラスのインスタンスを生成
m = Mesh2D([-4,4], [-4, 4], 65)

# 高度式を計算
Z = f(m.X, m.Y)

# Figureと3DAxesを追加
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection="3d")

# 3DAxesの設定
ax.set_xlabel("x", size=15)
ax.set_ylabel("y", size=15)
ax.set_zlabel("z", size=15)
ax.view_init(60, 30)

# 曲面を描画
ax.plot_surface(m.X, m.Y, Z, cmap="winter")

plt.show()

3次元格子点

2 次元格子を 2 個の 2 次元配列で表したように、3 次元格子を生成するためには 3 個の 3 次元格子が必要になります。3 次元格子を使うと、3 次元空間の各点で定義された関数 (3 変数関数) を計算することができます。このような関数はもはやグラフで可視化することはできませんが、機械学習などの高度なプログラムでは、多変数関数は頻繁に登場するので、多次元格子が生成される仕組みを理解しておく必要があります。

# NUMPY_MESHGRID_3D

# In[1]

import numpy as np

x = np.arange(1, 4)
y = np.arange(1, 4)
z = np.arange(1, 4)

# 3×3×3の三次元格子を生成
X, Y, Z = np.meshgrid(x, y, z)

# 3変数関数を定義
def func(x, y, z):
    return x**2 + y**2 + x*y*z

# 三次元格子の各点の値を関数に入れて計算
W = func(X, Y, Z)

print(W)

''
[[[ 3 4 5]
  [7 9 11]
  [13 16 19]]
 
 [[ 7  9 11]
  [12 16 20]
  [19 25 31]]
 
 [[13 16 19]
  [19 25 31]
  [27 36 45]]]
'''

同様に、N 次元格子を生成するには N 個の N 次元配列を組合わせます。

sparse grid (スパースグリッド)

numpy.meshgrid() には sparse というオプション引数にブール値を指定することができます。デフォルトでは False になっていますが、これを True にすると、sparse grid (スパースグリッド) とよばれる簡略化された格子点が生成されます。

# NUMPY_MESHGRID_SPARSE

# In[1]

import numpy as np

x = np.arange(1, 4)
y = np.arange(1, 4)

# sparse grid を作成
X, Y = np.meshgrid(x, y, sparse=True)
print(X)
print(Y)

# [[1 2 3]]

# [[1]
#  [2]
#  [3]]

sparse grid を使うとメモリを節約することができます。sparse grid の値を計算に用いる場合はブロードキャスト機能によって処理されます。

# In[2]

Z = X + Y
print(Z)

# [[2 3 4]
#  [3 4 5]
#  [4 5 6]]