Python で数学的なことを試すときに利用される Matplotlib/SciPy/pandas/NumPy についてサンプルコードを記載します。

• Matplotlib
• SciPy
• pandas
• NumPy

チュートリアル

• Installing packages
• Quickstart tutorial

インストール

pip でインストールする場合

pip のインストール

sudo apt install python-pip

ユーザ毎にインストールする場合

pip install --user numpy scipy matplotlib ipython jupyter pandas sympy nose
ls ~/.local/lib/python2.7/site-packages/

システム全体にインストールする場合

sudo pip install numpy scipy matplotlib ipython jupyter pandas sympy nose

apt/yum でインストールする場合

pip と比較してインストールされるバージョンが古くなります。

sudo apt install python-numpy python-scipy python-matplotlib ipython python-pandas python-sympy python-nose

新規 ndarray 配列の作成

import numpy as np
a = np.array([[1,0],
              [0,1]])

a = np.arange(15).reshape(3, 5)
a.resize(3,5) # 破壊的に自身を変更

型の明記

np.array( [ [1,2], [3,4] ], dtype=np.float64 )

すべて 0/1/未定 (dtype 未指定時は np.float64 )

np.zeros( (2,2), dtype=np.int64 )
np.ones( (2,2), dtype=np.int64 )
np.empty((2,2))

乱数

np.random.random((2,2))

10 以上 30 未満で +5

np.arange(10, 30, 5)
np.arange(0, 2, 0.3)

10 以上 30 以下で 10 要素

np.linspace(10, 30, 10)
np.linspace(0, 2*np.pi, 100)

コピー

b = a.copy()

一次元配列の取得

b = a.ravel()

条件を満たすインデックスを取得

a = np.array([1, 2, 3], dtype=np.uint8)
np.where(a == 2)[0][0]  #=> 1

配列としての演算

要素の合計値

np.sum(a)
a.sum()

要素の最大値

np.max(a)
a.max()

次元数

a.ndim

各次元の要素数

a.shape

全次元の要素数の合計

a.size

要素の型およびバイト数

a.dtype.name
a.dtype.itemsize
a.itemsize

各要素の値を二乗

a**2
a * a

各要素に対して条件判定

a < 2

各次元毎の演算

a.sum(axis=0) # 列j の合計
a.sum(axis=1) # 行i の合計

関数の各要素への適用

np.exp(a)
np.sqrt(a)
np.sin(a)
np.floor(a)

結合/分割

np.vstack((a,a))
np.hstack((a,a))
np.vsplit(a,999)
np.hsplit(a,999)

基本的な線形代数の演算

加減算

a + a
a - a

行列の積

a.dot(a)

転置行列

a.T
a.transpose()

逆行列

np.linalg.inv(a)

単位行列

np.eye(サイズ)

トレース

np.trace(a)

外積、内積、ノルム

from numpy import cross
from numpy import array
from numpy.linalg import norm

a = array([1,2,3])

cross(a, a)
a.dot(a)
norm(a)

軸を指定した、多次元配列 (テンソル) の転置

x = np.random.rand(1, 2, 3, 4)
y = x.transpose(1, 0, 3, 2)
print(x.shape)  #=> (1, 2, 3, 4)
print(y.shape)  #=> (2, 1, 4, 3)

可視化

画像データの表示 imshow

import numpy as np
import matplotlib.pyplot as plt

myimg = np.array([
  [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], 
  [[0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]])

plt.imshow(myimg)
plt.show()

グラフの表示 plot

import numpy as np
import matplotlib.pyplot as plt

a = np.array([9,1,9,1])

plt.plot(a)
plt.show()

移動平均による平滑化

import numpy as np
a = [1,1,1,1,1]

def MovingAverage(a, n=3):
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n - 1:] / n

MovingAverage(a)

同次変換

ある座標系の点 (x,y,z) を別の座標系で表現したときの座標 (X,Y,Z) は、二つの座標系の平行移動と回転移動を表現する変換行列を用いて以下のようになります。

|X|   |      |tx|   |x|
|Y| = |  R   |ty| . |y|
|Z|   |_ _ _ |tz|   |z|
|1|   |0 0 0 | 1|   |1|

例えば、ある座標系が、ワールド座標系から (10, 10, 10) 平行移動されて更にワールド座標系の z 軸まわりに -90 度回転したものであった場合、その座標系における点 (1,0,0) はワールド座標系において以下のように計算できます。

import numpy as np

theta = 90 * (np.pi/180.0)

transform = np.array([
    [np.cos(theta), -np.sin(theta), 0, 10],
    [np.sin(theta),  np.cos(theta), 0, 10],
    [            0,              0, 1, 10],
    [            0,              0, 0,  1]
])

pos = np.array([
    [1],
    [0],
    [0],
    [1]
])

posWorld = transform.dot(pos)

# array([[10.],
#        [11.],
#        [10.],
#        [ 1.]])

回転行列は一方の座標系から見たときの他方の座標系の単位ベクトルを並べたものです。単位ベクトルの長さは 1 であり互いに直行するため、回転行列は正規直交行列です。正規直交行列には必ず逆行列が存在することが保証されており、転置行列が逆行列となります。

9 要素のうち、単位ベクトルの長さが 1 であることの制約が三つ、単位ベクトルが互いに直行することによる制約が三つ、残りの自由度は 3 となります。この自由度 3 で回転が表現されます。例えば三つの軸回りそれぞれで回転すれば任意の回転を表現できます。

したがって変換行列 $TT$ にも逆行列 $T^{-1}T^\{-1\}$ が存在することが保証されており、以下のようになります。

OpenRAVE 等で座標変換を行う場合は、存在が保証されている逆行列を用いて以下のような記述が可能になります。回転行列の転置行列から逆行列を計算してもよいですが、以下では numpy.linalg.inv を用いています。

from numpy.linalg import inv

manip = robot.GetActiveManipulator()
tManipToWorld = manip.GetEndEffectorTransform()
tWorldToManip = inv(tManipToWorld)

target = env.GetKinBody('mug1')
tBodyToWorld = target.GetTransform()

tBodyToManip = tWorldToManip.dot(tBodyToWorld)

把持された対象とマニピュレータとの位置関係 tBodyToManip はロボットの姿勢が変化しても把持されている間は変化しないことも重要です。

その他の用語

• AABB (Axis Aligned Bounding Box)

バイト処理

8 バイトを 16 バイトで確認すると以下のようになりました。リトルエンディアンのマシンということになります

In [26]: np.array([255, 0], dtype=np.uint8).view(dtype=np.uint16)
Out[26]: array([255], dtype=uint16)

In [27]: np.array([0, 1], dtype=np.uint8).view(dtype=np.uint16)
Out[27]: array([256], dtype=uint16)

In [44]: np.array([255], dtype=np.uint16).view(dtype=np.uint8)
Out[44]: array([255,   0], dtype=uint8)

In [45]: np.array([0b0000000011111111], dtype=np.uint16).view(dtype=np.uint8)
Out[45]: array([255,   0], dtype=uint8)

In [46]: np.array([0x00FF], dtype=np.uint16).view(dtype=np.uint8)
Out[46]: array([255,   0], dtype=uint8)

バイト列をビッグエンディアンに変換することができます。

In [56]: np.array([255, 0], dtype=np.uint8).view(dtype=np.uint16).tobytes()
Out[56]: b'\xff\x00'

In [57]: np.array([255, 0], dtype=np.uint8).view(dtype=np.uint16).byteswap().tobytes()
Out[57]: b'\x00\xff'

文字列のバイト処理

バイト列を ndarray に変換 frombuffer()

In [51]: byteArray = b'\x61\x62\x63'

In [52]: data = np.frombuffer(byteArray, dtype=np.uint8)

In [53]: data
Out[53]: array([97, 98, 99], dtype=uint8)

ndarray をバイト列に変換 tostring()

In [8]: encoding = 'utf-8'

In [9]: data.tostring()
Out[9]: 'abc'  <-- python2
Out[9]: b'abc'  <-- pyhton3

In [10]: data.tostring().decode(encoding)
Out[10]: u'abc'  <-- python2
Out[10]: 'abc'  <-- pyhton3

In [11]: data.tostring().decode(encoding).encode(encoding)
Out[11]: 'abc'  <-- python2
Out[11]: b'abc'  <-- pyhton3

python2 の str は bytes です。

In [2]: str == bytes
Out[2]: True

In [3]: str == unicode
Out[3]: False

python3 の str は python2 の unicode です。

In [1]: str == bytes
Out[1]: False

In [2]: str == unicode
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-2c1b95aac9eb> in <module>
----> 1 str == unicode

NameError: name 'unicode' is not defined

デバイスファイルの読み書き

memmap を利用すると ndarray をファイルシステム上のバイナリファイルにマッピングできます。ブロックデバイスを作成することができます。

import numpy as np

filePath = '/dev/shm/mydev'

fp = np.memmap(filename=filePath, dtype=np.uint8, mode='w+', shape=(5))
fp[:] = np.array([1, 2, 3, 4, 5], dtype=np.uint8)

fp2 = np.memmap(filename=filePath, dtype=np.uint8, mode='r+', shape=(5))
In [14]: fp2
Out[14]: memmap([1, 2, 3, 4, 5], dtype=uint8)

$ du -b /dev/shm/mydev
5       /dev/shm/mydev