grbl1.1+Arduino CNCシールドV3.5+bCNCを使用中。
BluetoothモジュールおよびbCNCのPendant機能でスマホからもワイヤレス操作可能。
その他、電子工作・プログラミング、機械学習などもやっています。
MacとUbuntuを使用。

CNCマシン全般について:
国内レーザー加工機と中国製レーザー加工機の比較
中国製レーザーダイオードについて
CNCミリングマシンとCNCルーターマシンいろいろ
その他:
利用例や付加機能など:
CNCルーター関係:



*CNCマシンの制作記録は2016/04/10〜の投稿に書いてあります。


2018年5月2日水曜日

Bayesian Optimization / Gaussian Process

前回、matplotlibのアニメーションについては大体使えるようになったので、「ベイズ最適化・ガウス過程」を実験してみました。いろいろと専用のライブラリがあるのですが(GPy、GPFlow、Edward、Pyroなど)、今回はscikit-learn Gaussian Processesを使っています。

ベイズ最適化はハイパーパラメータチューニングに使われるようで、ブラックボックス関数における最大値(あるいは最小値)を全領域をスキャンしなくても探し出すことができる方法のようです。

今回の場合は、グラフ内の赤破線の最大値がどこにあるか探索しています。赤破線を見ると、左と右にピークが二箇所ありますが、微妙に右側のほうが高いという設定にしてあります。それを間違えないように探し出せるかという実験です。
右側の赤丸が正解(最大値)の地点です。この地点は、x座標-2〜2までの範囲を0.1刻みで200回スキャンし、その結果から最大値を求めています。
しかし、ベイズ最適化の場合は、200回スキャンしなくても約10回程度で、どこが最大値なのかを効率よく探してくれるようです。

今回の手順として:
・カーネル:Matern、獲得関数:UCBを選択
・xの範囲(-2〜2)からランダムにまず2点(青丸:start印)選んで初期観測値とする
・その観測値から、全域のμ(青破線)とσ(青塗)を求める
・μとσを獲得関数(UCB)に代入
・獲得関数の結果(緑実線)から最大値の箇所を見つける
・獲得関数の最大値から、それに相当するxの位置を見つける
・xの位置から、それに相当するyの位置を求める
・ここで求めたxとyの座標を次回の探索箇所にする
・次回探索箇所を観測値リストに加える
・あとは繰り返し

グラフ下部(緑実線)にあるのが獲得関数UCBの値です。緑実線の一番高いところが、次の探索地のx座標となります。
今回使用したUCBの中身は単純で、

UCB最大値 = mu + kappa * sigma

となっています。係数kappa=1.0の場合は、μ値(青破線)にσ値(青塗)を足した値がそのままUCB最大値となるので、グラフ内の青塗を含めた一番高い箇所が次の探索箇所となります。μ値が高いところほど選ばれがちですが、kappaを大きくするとσ値も大きくなるので、μ値が低い箇所でも探索される可能性があがります。


%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern
import matplotlib.animation as animation
from IPython.display import HTML

np.random.seed(111)

def f(x):
    return -np.sin(x*4.5)*x*0.3+0.02*x                  # 真の関数(赤破線表示)

x = np.linspace(-2,2,200)                               # x値範囲
y = f(x)                                                # 真のy値

fig = plt.figure(figsize=(12,6))                        # 画面サイズ設定
ax = plt.subplot()                                      # サブプロットを用意
ax.axis([x.min(),x.max(),y.min()*2,y.max()*2])          # グラフ範囲設定(x,yの値に依存)
ax.set_aspect(1.0)                                      # アスペクト比設定

ax.plot(x,y,'r--',alpha=0.8)                            # 真の関数曲線をプロット(赤破線)
ax.plot(x[np.argmax(y)],y.max(), 'ro',alpha=0.8)        # y正解値;x全域スキャンによるy最大値(赤丸)

init_pt = 2                                             # 初期観測値個数
xi = np.random.uniform(x.min(),x.max(),init_pt)         # 初期観測値:xの範囲から2箇所ランダム選択
yi = f(xi)                                              # 初期観測値:yの値
ax.plot(xi,yi,'bo')                                     # 初期観測値プロット(青丸)
for i in range(len(xi)):
    ax.text(xi[i],yi[i],'  start'+str(i),color='b')     # 初期観測値:'start'文字表示

def UCB(mu_x, sigma_x, kappa=1.0):                      # 今回使用する獲得関数:upper_confidence_bound
    return mu_x + kappa * sigma_x                       # kappa増大-->σ増大-->探索領域拡大

kernel = Matern(length_scale_bounds="fixed")            # 今回使用するカーネル:Matern
gp = GaussianProcessRegressor(kernel=kernel,alpha=0.01) # scikit-learnガウス過程回帰オブジェクト

ims = []                                                # アニメ用リスト
itr = 15                                                # 探索ループ回数

for i in range(itr):
    gp.fit(xi[:,None],f(xi))                                 # 観測値フィッティング(xiをリシェイプして代入)
    mu_x, sigma_x = gp.predict(x[:, None], return_std=True)  # μ(x)とσ(x)を獲得(事後分布)

    a_x = UCB(mu_x, sigma_x, kappa=2.0)                      # UCBでa(x)最大値を得る
    index = np.argmax(a_x)                                   # a(x)最大値のindexを得る
    x_next = x[index]                                        # a(x)最大値のindexから次のx探索値を決定
    y_next = f(x_next)                                       # y探索値決定
    
    P1 = ax.plot(x, mu_x, 'b--')                                               # μ値(青破線)
    P2 = ax.fill_between(x, mu_x-sigma_x, mu_x+sigma_x, color='b', alpha=0.2)  # σの範囲(青塗)
    P3 = ax.plot(xi,yi,'bx',markersize=8)                                      # 探索済み箇所(青クロス)
    P4 = ax.plot(x_next, y_next, 'mo')                                         # 現在探索箇所(マジェンタ丸)
    P5 = ax.plot([x_next,x_next],[y.min()*2, y_next] ,'m-', alpha=0.6)         # 現在探索箇所(マジェンタ縦線)
    P6 = ax.plot(x_next,(a_x[index]-a_x.min())/2+y.min()*2, 'go',alpha=0.8)    # a(x)最大値:σ(青丸)
    P7 = ax.plot(x, (a_x-a_x.min())/2+y.min()*2, 'g-')                         # 獲得関数値:a(x)
    P8 = ax.plot([xi,xi],[[y.min()*2]*len(yi),yi], 'b-', alpha=0.2)            # 探索済み箇所(青縦線)
       
    xi = np.r_[xi,x_next]                                  # 現在探索箇所をリストに加える(r_: concatenateの代用)
    yi = np.r_[yi,y_next]
    
    ims.append(P1 + [P2] + P3 + P4 + P5 + P6 + P7 + P8)    # アニメ用リストへ各プロットを入れる
    
print('x :',x[y.argmax()], '   y :',y.max())               # 正解値出力(y値)
print('xi:',xi[yi.argmax()], '   yi:',yi.max())            # 最終最大値出力(yi値)

ani = animation.ArtistAnimation(fig, ims, interval=1000)   # アニメ関数
#ani.save("gp_sim.gif", writer = "imagemagick")            # GIFアニメ保存
plt.close()
HTML(ani.to_html5_video())                                 # HTML5 Video 出力(mp4ダウンロード可)
#HTML(ani.to_jshtml())                                     # JavascriptHTML 出力

GaussianProcessRegressor()にはカーネルを指定する箇所があり、今回は、scikit-learnのMaternを使っています。ちなみにカーネルを指定せずデフォルトのままだと、以下のような感じになります。
σ(青塗)の領域が結構変わります。
scikit-learnのこのページにカーネルの違いについて書いてあります。

ハイパーパラメータチューニングと言っても、カーネルの選び方でもかなり違いがでてきます。
また、獲得関数もUCB以外にいくつかあるので、それも選ばなければいけないという感じです。


参考にしたサイト:
Scikit-learn Gaussian Processes
llSourcell/hyperparameter_optimization_strategies
ベイズ最適化入門
機械学習のハイパーパラメータ探索 : ベイズ最適化の活用

関連:
GPyOpt: Digital Recognizer(MNIST) CNN Keras ハイパーパラメータ最適化

2018年4月23日月曜日

Jupyter Notebook: matplotlib /fill() + fill_between() Animation

前回のmatplotlibアニメーションの続き。
今回は領域を塗るfill()とfill_between()関数を使ってみました。静止画なら簡単なのですが、アニメーションになると思った通りにいかないことが多々あったのでメモ。

まず、アニメーションでは:

ArtistAnimation()
FuncAnimation()

がありますが、今回は比較的簡単なArtistAnimationを使用。
追記:
FuncAnimation()の手抜きのやり方も追加しておきました(投稿中程)。

そしてグラフ領域の塗り分けについては、

fill(x, y)
fill_between(x, y1, y2)

があるのですが、
fill()に渡すxとyは、塗る領域を指定するリストにしなければならない。これが少し面倒。
fill_between()の場合は、xはそのままxの数式、y1に下限の数式、y2上限の数式を入れればいいだけなので簡単です。

今回つくりたいグラフは以下のような感じです(今回つくったGIF動画です)。fill_between()を使えば簡単にできそうですが結構ハマりました。
このような領域を塗るサンプルを見てみると、アニメーションのときはfill()を使っていることが多かったのですが、塗る領域を座標指定するのが少々面倒。数式代入だけで済むfill_between()を使うと、

TypeError: unsupported operand type(s) for +: 'PolyCollection' and 'list'

こんな感じのエラーが出て、ポリコレクションとリストを一緒に入れるなみたいに怒られます。便利な分、複雑なデータ形式でそのままでは使えないのかと、ポリコレクションを調べたり、その内部の値にアクセスする方法はないのか探したりしていました。どうやらfill_between()はアニメーションでは使いにくいのか、参考になるサンプルもあまりありませんでした。しかし、非常に単純な理由でfill_between()も使えるということがわかりました。なんだそんなことかという感じです。


結論:
fill_between()が簡単。
ArtistAnimationを使うなら、単にプロットした内容を空のリストにappend()していけばいいのですが、そのままfill_between()のプロットを入れると、上記のようなエラーが出て受け付けません。そこで、[ ]で括って無理やりリスト形式にして入れると動きました。たったそれだけの違いでした。

これに気づくまでかなり時間かかりました。結局は少し面倒だけど、fill()を使うしかないのかなと思っていたところ、ようやく上手くいったという感じです。


Fill_between() + ArtistAnimation()の場合:
一応、Jupyter Notebookを使う前提です。
今回は、HTML()で出力するのでバックエンドは「inline」。前回の投稿に追記しましたが、ipymplがバージョンアップしたので、それを使うこともできます(その場合は「nbagg」になる)。ただし、「inline」のほうが安定しているので、今回はそのまま。


%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation
from IPython.display import HTML

x = np.arange(0,10,0.2)     # x値範囲:(xmin,xmax,pitch)
y = x*np.sin(x)             # メイン関数
sigma = np.cos(x*2) + 2     # シグマ計算式
y1 = y - sigma              # シグマ下限
y2 = y + sigma              # シグマ上限

fig = plt.figure()
plt.axis([0,10,-10,10])     # グラフ範囲:[xmin,xmax,ymin,ymax]

ims = []                    # アニメ用プロット格納リスト

for i in range(len(x)):
    P0 = plt.plot(x,y,'b--',alpha=0.3)                                    # 固定描画(青破線)
    P1 = plt.fill_between(x[:i+1],y1[:i+1],y2[:i+1],color='r',alpha=0.2)  # シグマ値描画 (塗り幅:赤帯)
    P2 = plt.plot(x[:i+1],y[:i+1], 'r-')                                  # 軌跡描画(赤実線)
    P3 = plt.plot(x[i],y[i], 'ro')                                        # 移動描画(赤丸)
    ims.append(P0 + [P1] + P2 + P3)     # 各プロット格納:fill_beteenのプロットP1だけ[]で括って入れる
    
ani = animation.ArtistAnimation(fig=fig, artists=ims, interval=100)       # アニメ関数
plt.close()

HTML(ani.to_html5_video())                                                # HTML5 Video 出力
#HTML(ani.to_jshtml())                                                    # JavascriptHTML 出力

先ほど書いたように、ims=[]に複数のプロットを格納して行く際、fill_between()でつくったプロットP1だけ[]で括って入れます。

ims.append(P0 + [P1] + P2 + P3)

こうすれば、ポリコレクションであっても受け付けてくれるようです。
それから今回気づきましたが、変化しない固定グラフはアニメーションリストに加えなくてもいいと思っていましたが、そうするとアルファ値などが反映されなくなってしまう(どんどん上塗りされて濃くなっていく)ので、リストに加えることにしました。


追加:Fill_between() + FuncAnimation()の場合:
追加しました。
この組み合わせだけ上手くいかず、いろいろ実験していましたが、またもや非常に簡単な方法で表示可能となりました。
今回はフレーム数を少なめにして、そのかわり各座標にポリゴン描画する順番で塗面の座標値に番号を記入してみました。
fill_between()なので、fill()のようにポリゴン座標値は必要ないのですが、位置などずれていないか確かめるために挿入してみました。以下に続く、fill()を使う場合にはこの座標値の順番が参考になるかと。

FuncAnimationの場合、plot関数に予め空の[]を入れておき、update(i)関数内でその値を更新していくというパターンが多いのですが、今回は表示させたいものを全てupdate(i)関数内に入れてしまい無理やりアニメーション化しました。しかしそうすると画面がリフレッシュされないために、最初にax.cla()を入れておきました。特に渡す変数などなくて非常に簡単です。


import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
plt.style.use('ggplot')

fig = plt.figure(figsize=(5,4))
ax = plt.subplot()

x = np.arange(0,10,0.4)
y = x*np.sin(x)
s = np.cos(x) + 2
y1 = y - s
y2 = y + s

def update(i):
    ax.cla()                                                          # 事前に画面クリアで内容更新
    ax.axis([0,10,-8,8])                                              # 座標範囲設定
    ax.plot(x,y,'r--',alpha=0.2)                                      # メイン関数固定プロット(赤破線)
    ax.plot(x[:i+1],y[:i+1],'r-',lw=3)                                # 軌跡描画(赤実線)
    ax.fill_between(x[:i+1],y1[:i+1],y2[:i+1],color='r',alpha=0.1)    # 塗面描画(透過赤:10%)
    
    xs = np.concatenate([x[:i+1],x[i::-1]])                           # 塗面ポリゴン座標番号リスト
    ys = np.concatenate([y1[:i+1],y2[i::-1]])
#     ax.plot(x[:i+1],y1[:i+1],'b+')                                  # 塗面下端描画(青+:現在非表示)
#     ax.plot(x[:i+1],y2[:i+1],'g+')                                  # 塗面上端描画(緑+:現在非表示)
    for j in range(len(x[:i+1])*2):    
        ax.text(xs[j],ys[j],str(j),fontsize=10,color='m')             # 塗面ポリゴン座標番号表示
        ax.plot(xs[j],ys[j],'r+')                                     # 各座標位置を+で表示
    
    return fig,                                                       # とりあえずfigをリターン

ani = animation.FuncAnimation(fig, update,frames=len(x),interval=200) # Funcアニメーション使用
# ani.save("fb_func.gif", writer = "imagemagick")                     # GIFアニメ保存
plt.close()
HTML(ani.to_html5_video())                                            # HTML5 Video 出力
#HTML(ani.to_jshtml())        


特に難しいところはないのですが、update(i)関数からのリターンがないとエラーがでるので、特に影響のでないfigをリターンさせています。ax.cla()で毎フレームクリアしているためか、表示範囲が崩れることがあります。そのため、すぐにax.axis()で表示範囲を設定する必要があります。
基本的に、plot()やfill_between()、そしてtext()など、表示させたいものをupdate(i)関数内に入れてあるだけです。本来のFuncAnimation()のやり方のようにきちんとset関数などで値を更新してあげればax.cla()でリフレッシュする必要もないのでしょうが、とりあえず非常に簡単に表示できるのでこの方法も便利かと。
fill_between()におけるx、y1、y2のset関数が見当たらず、Pathを使って自前で塗面描画するはめになったりと、かえって手間がかかってしまったので、今回のようなやり方が発見できてよかったです。


x[i]とx[:i+1]について:
ここもちょっとハマりました。
ArtistAnimationの場合は毎回ループでi番目をプロットして見せるだけなので、そのまま見せたいプロット関数を並べておけばいいだけです。
しかし今回のように、過去の軌跡も表示するグラフ、変化する現在位置だけを表示するグラフが共存すると表示されるindex番号が一コマずれてしまうことが発生。アニメーションのピッチがこまかければ気づきにくい点ですが、よくみるとx[i]とx[:i]の現在地がずれていました。
そのため、現在地だけを表示するx[i]を基準にして、軌跡も表すx[:i]についてはx[i+1]にして表示する部分を揃えることにしました。
逆に、x[:i]を基準にして、x[i]をx[i-1]にしてしまうと最後のコマが表示されなくなるので、x[i]を基準にしたほうがいいということになったわけです。

このindex番号については、次に書くfill()において結構悩みました。


Fill() + ArtistAnimation()の場合:
fill_between()と比較するために、fill()のコードも書いてみました。やはり少し面倒でした。数式をそのまま代入すればいいfill_between()とは違って、塗面の範囲を座標を使って囲んでいかなければいけません。


%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation
from IPython.display import HTML

x = np.arange(0, 10, 0.2)                                    # x値範囲:(xmin,xmax,pitch)
y = x*np.sin(x)                                              # 描画用基本関数(任意)
sigma = np.cos(x*2) + 2                                      # シグマ値(線幅変化用)

y_high = y + sigma                                           # シグマ上限値
y_low = y - sigma                                            # シグマ下限値

fig = plt.figure()
plt.axis([0,10,-10,10])                                      #グラフ範囲指定:[xmin,xmax,ymin,ymax]

ims = []                                                     # アニメプロット格納リスト

for i in range(len(x)): 
    P0 = plt.plot(x,y,'g--',alpha=0.3)                       # 固定描画(緑破線)
    
    xf = np.concatenate([x[:i+1], (x[:i+1])[::-1]])          # fill用x: [x順列,x逆列]
    yf = np.concatenate([y_low[:i+1], (y_high[:i+1])[::-1]]) # fill用y: [y下限順列,y上限逆列]
    P1 = plt.fill(xf,yf,facecolor='b', alpha=0.2)            # シグマ値描画 (塗り幅:青帯)
    
    P2 = plt.plot(x[:i+1],y[:i+1],'b-')                      # 移動点軌跡描画(青実線)    
    P3 = plt.plot(x[i],y[i],'bo')                            # 移動点描画(青丸)
    
    ims.append(P0 + P1 + P2 + P3)                            # 各プロットをアニメ用リストに格納

ani = animation.ArtistAnimation(fig,ims,interval=100)        # アニメ関数
# ani.save("anim.gif", writer = "imagemagick")               # GIFアニメ保存(Ubuntu用にImagemagick使用)
plt.close()
HTML(ani.to_html5_video())                                   # HTML5 Video 出力(mp4ダウンロード可)
# HTML(ani.to_jshtml())                                      # JavascriptHTML 出力

内容はほぼ同じですが、fill()を使う場合はnp.concatenate()を使って、塗る領域の座標値リストをつくらなければいけないようです。そのまま式をいれただけだときちんと表示されませんでした。
また、ims.append()にはfill()を格納した変数をそのまま入れることができます。fill_between()のときのように[]でくくって入れる必要はないので、その分使いやすかったのですが、今となってはfill_between()の方が簡単。
今回の場合、塗る領域の下端の線と上端の線を連結したリストをつくり、それをfill()に渡してあげればいいようです。
以下がそれの図解です。
赤破線がメイン関数で、それに沿って帯状に塗るわけですが、緑の下端線0〜6と上端線7〜13までをつないだリストをつくるということです。0と13が閉じていませんが、内部的な処理では最後の座標で図形を閉じるので自動的に0と13はつながるようです。
面倒なのは、下端線が左から右、上端線が右から左という順番で並んでいるので、その順番にリスト内の要素も並べないといけないという部分。
ポリゴン図形を描くように、領域をぐるりと回るような座標取りになっています。
x座標は単なる往復でいいのですが、yは下端(順方向)と上端(逆方向)を組み合わせなければならないというわけです。

しかも今回は表示される区間が変化するため、もう少し面倒になります。
ということから、i番目まで表示するなら、

xf = x[:i+1] + x[i::-1]

となります。

x[:i+1] = i番目までの値全てのリスト
x[i::-1] = i番目以下すべてを反転したリスト

ただし、np.concatenate()を使って連結します。

xf = np.concatenate([ x[:i+1], x[i::-1] ])

二重の[]で括られている部分に要注意。
この手続きを踏んでfill()に渡してあげればいいというわけです。

しかし、ここでまたindexをiにするかi+1にするかで一コマずれの表示が発生したり、連結した座標値リストのサイズがxとyとで異なりエラーがでたりと、少々面倒です。最終的に調整した結果、ソースコードのようになっています。


リストのindex番号について:
x[:i]とかx[i::-1]など分かりにくいので以下にまとめてみました。
尚、リバース関数x.reverse()を使うと、元のxが反転してしまうようなので、今回は一時的に反転したりしなかったりするので使わないことに。そのかわり、リスト内全ての値の並びを反転するにはx[::-1]を使っています。
例えば、以下のようなリストがあったとします、


x = np.arange(5)                             # [0 1 2 3 4]
i= 4                                         # 0 1 2 3 4 
print(x[i])                                  # 4            i番目の値
print(x[:i])                                 # [0 1 2 3]    i未満全て(i含まず)
print(x[i:])                                 # [4]          i以上全て

#以下、反転パターン
print(x[::-1])                               # [4 3 2 1 0] 全ての反転
print(x[i::-1])                              # [4 3 2 1 0] i以下全ての反転
print((x[:i])[::-1])                         # [3 2 1 0]    i未満の全ての反転

print(np.concatenate([x[:i],x[i::-1]]))      # [0 1 2 3 4 3 2 1 0]
print(np.concatenate([x[:i],(x[:i])[::-1]])) # [0 1 2 3 3 2 1 0]


fill()を使う際どの組み合わせがいいのかいろいろ試してみました。要素数が一つ足りなくてエラーがでたときもありましたが、表示上一コマずれていたりと微妙に違います。
また、先ほど書いたように、x[i]を基準とすると(i番目の座標を表示しようとすると)、過去の軌跡表示に使用するx[:i]の場合は、i番目未満までの表示となってしまうため、x[:i]をx[:i+1]に変更して、現在位置とそれ以前の軌跡を表示しなければいけなくなり、結構面倒です。ちなみに、現在位置と軌跡を含めた表現に使うx[:i+1]を基準に比較してみると、、


x = np.arange(5)                                 # [0 1 2 3 4]
i= 4                                             # 0 1 2 3 4 
print(x[i])                                      # 4          
print(x[:i+1])                                   # [0 1 2 3 4]

print(x[i+1::-1])                                # [4 3 2 1 0]
print((x[:i+1])[i::-1])                          # [4 3 2 1 0]

print(np.concatenate([x[:i+1],x[i+1::-1]]))      # [0 1 2 3 4 4 3 2 1 0]
print(np.concatenate([x[:i+1],(x[:i+1])[::-1]])) # [0 1 2 3 4 4 3 2 1 0]

このようになり、こちらのほうが今回の場合は都合がよさそうです。特に最後の行はfill()に代入するために連結させた結果、ひとつ前の例だと、

[0,1,2,3,4,3,2,1,0]

のように4が真ん中にひとつしかなかったり、

[0,1,2,3,3,2,1,0]

のように4が抜け落ちていたりすることがあります。
今回必要だったのは、

[0,1,2,3,4,4,3,2,1,0]

だったので、ソースコードに書いてある方法を採用しました。
特に、

x[:i] = [0,1,2,3]
x[:i+1] = [0,1,2,3,4]

の違いだったり、

x[i::-1] = [4,3,2,1,0]
x[:i+1::-1] = [4,3,2,1,0]
(x[:i+1])[::-1] = [4,3,2,1,0]

が同じだったりすることで、連結した結果が微妙に異なりました。最終的には、やや見づらいけれども、

(x[:i+1])[i::-1] = [4,3,2,1,0]

を使うことにしました。
この場合、x[:i+1]を一旦括弧()でくくって、それを[::-1]を使って反転させています。というのは、iに最大値4を入れたときはいいのですが、最小値0を入れるとまた微妙に結果が変わってくるからです。


x = np.arange(5)                                 # [0 1 2 3 4]
i= 0                                             # 0 1 2 3 4 
print(x[i])                                      # 0
print(x[:i+1])                                   # [0]

print(x[i+1::-1])                                # [1 0]
print((x[:i+1])[i::-1])                          # [0]

print(np.concatenate([x[:i+1],x[i+1::-1]]))      # [0 1 0]
print(np.concatenate([x[:i+1],(x[:i+1])[::-1]])) # [0 0]

先ほどは下2行の結果は同じでしたが、今回は違います。求めている結果は[0,1,0]ではなく[0,0]なので最後のやり方を使いました。
この辺の調整は理屈というよりも、やりながら挙動を確かめて修正という感じです。


まとめ:
ということから、fill()を使うといろいろと勉強にはなったけど面倒です。やはり、fill_between()のほうがはるかに簡単でした。

関連:

2018年4月19日木曜日

Jupyter Notebook/JupyterLab: matplotlib GIFアニメーション 作り方+読込み+表示方法

前回に引き続き、Jupyter NotebookとJupyterLabでのmatplotlibをつかったアニメーションについてのメモ。matplotlibは便利なのだけれども、個人的には使いづらく(覚えにくい)、特にアニメーションではかなりハマったので書いておきます。

環境:
Ubuntu 16.10
Python3.5
Jupyter Notebook 5.0.0
JupyterLab 0.32.0
matplotlib 2.2.2


まず、Jupyter NotebookやJupyterLab上に表示する際、
・%matplotlib inline (主に静止画)
・%matplotlib nbagg (動画、インタラクティブ)
この二種類のbackendsを用いるのですが、基本的にJupyterLabだと「nbagg」のほうはダメみたい。

次に、グラフをアニメーション化するには二つの方法があるようで、
・ArtistAnimation (比較的簡単)
・FuncAnimation (やや複雑)
パラメータが少なければどちらもそう難しくはないのですが、場合によっては表示されない場合もあります。表示されなくてもGIFアニメーションに保存はできることもあります。

これはつくったGIFアニメ。比較的簡単につくることができるけれども、JupyterLabに表示させるのが難しい。

表示方法(3種類):
そして、Jupyter上で表示する方法としては3つあり、

・nbaggを用いたインタラクティブウィンドウ(JupyterLabではエラー)
・inlineのままJavascript HTMLで出力する方法(JupyterLabではエラー)
・inlineのままHTML5 Videoとして出力する方法(エラーなし)

バックエンドをinlineのままJavascript HTMLやHTML5 Videoを使うためには、

from IPython.display import HTML


が必要となります。

以下がそれぞれの見た目です(キャプチャ画像なのでボタンを押しても動きません)。

バックエンド:%matplotlib nbagg
表示:animation.ArtistAnimation
JupyterLabでは表示不可


バックエンド:%matplotlib inline
変換表示:JavascriptHTML
JupyterLabでは表示不可
この中では一番多機能、逆再生・往復再生・スピード調整など可能

バックエンド:%matplotlib inline
変換表示:HTML5 Video
Jupyter Notebook、JupyterLabともに表示可能
動画ダウンロード・フルスクリーン表示可能


Jupyter Notebookではどれも大丈夫でしたが、JupterLabではJavascriptのエラーなのか最初の二つはダメそうです。HTML5 Video出力だと大丈夫だったので、これがおすすめかもしれません。
バグなのか、設定方法が悪いのか、インストールが不十分なのかわかりにくく結構手間取ってしまいました。

追記:
尚、ボタン付きのアニメーションとしてWebページに埋め込むには、以下を参考にしてください。
Matplotlib Animation embed on web page:アニメーションのWebページ上への埋め込み


今回試してみたこと:
・グラフをアニメーション化して再生、GIFアニメーション保存
・つくったGIFアニメーションを読み取って再生する方法
・Jupyter上にマークダウンとしてGIFを表示する方法

下準備(インストール:このサイト参照):
Nodejs
npm
ipympl
jupyter-widgets/jupyterlab-manager
(It requires matplotlib 2.0 or and ipywidgets 7.0 more recent.)



グラフのアニメーション化・再生・GIFアニメーション保存:
「ArtistAnimationを使った場合」

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation
from IPython.display import HTML

fig = plt.figure()
x = np.linspace(-5,5,100)
y = np.sin(x)

ims = [] # 動画要素を格納するリスト

for i in range(len(x)):
    ln = plt.plot(x,y, 'b--')     # 動かない線
    pt = plt.plot(x[i],y[i],'ro') # 移動する点
    ims.append(ln + pt)

ani = animation.ArtistAnimation(fig,ims,interval=50) # アニメ関数

ani.save("test.gif", writer = "imagemagick") # GIF保存

plt.close() #余計なウィンドウがでた場合消す

HTML(ani.to_html5_video()) # HTML5 Videoで表示
# HTML(ani.to_jshtml())    # Javascript HTMLで表示


アニメーション化するには、animation.ArtistAnimation()に最低でもfigとimsの二つの引数を渡し、imsに動画になる要素としてのplt.plot()を必要なだけ+を使って入れていくだけ。
その下のani.save()でGIFとして保存。Ubuntuの場合はwriterとしてimagemagickを指定したほうがいいらしい。
plt.close()については、バックエンドで%matplot inlineがはたらいているためか、HTML5 Videoのウィンドウとは別にもうひとつplt.figure()の空ウィンドウが出てくることがあって、それを消しています。
Jupyter Notebookなら、HTML(ani.to_jshtml)に切り替えても表示できました。nbaggを使わずinlineで済むのでお手軽です。ただし、HTMLをインポートしておく必要があります。
%matplotlib inlineは切り替えなければデフォルトになっているようなので、途中でnbaggなどに変えなければ記入する必要はありません。


「FuncAnimationを使った場合」

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation
from IPython.display import HTML

x = np.linspace(-5,5,100)
y = np.sin(x)

fig = plt.figure()
plt.axis([-5,5,-2,2])      # 座標範囲:[xmin,xmax,ymin,ymax]
b, = plt.plot([],[],'b--') # 固定描写(青破線)
t, = plt.plot([],[],'r-')  # 移動描写(赤丸)

# plt.axis('off') # 枠線なしの場合

def update(i):
    b.set_data(x,y)           # 固定描写(青破線)
    t.set_data(x[i], y[i])    # 現在の移動座標のみ表示(点描の場合)
    #t.set_data(x[:i], y[:i]) # 過去の軌跡座標も表示(線描の場合)

ani = animation.FuncAnimation(fig, update, interval=50) # アニメ関数

ani.save("test.gif", writer = "imagemagick") # GIF保存

plt.close()

HTML(ani.to_html5_video())   # HTML5 Videoで表示
# HTML(ani.to_jshtml())      # Javascript HTMLで表示


こちらの方法だとupdate(i)関数を作って、フレーム数となる引数iを入れてFuncAnimation()に渡すことになります。他の方法でも可能なのかもしれませんが、サンプルなどではline2Dオブジェクト用の変数(t, : 変数の後ろにコンマ付き)を用意し、set_data()で更新していくようです。やや複雑なことが可能なのかもしれません。subplotを追加するサンプルが多いですが、plt.plotだけを使って書いてみました。


バックエンド(inline/nbagg)切替方法・確認方法:
以下のコマンドで現在使用されているバックエンドを切替・確認できます。


%matplotlib inline
# %matplotlib nbagg

import matplotlib
matplotlib.get_backend()


この状態だと、inlineが使用されているので以下が出力されます。


'module://ipykernel.pylab.backend_inline'


通常の静止画グラフやHTML5 Video、JavascriptHTMLを使うならばinlineに設定。
また、コメントアウトでnbaggに切り替えれば、

'nbAgg'


が出力されて、nbaggがバックエンドとして使用されていることがわかります。ただし、JupyterLabではJavascriptエラーがでるかもしれません。


GIFアニメの読込と再生:
先ほどの方法でGIFアニメをつくることができますが、それを再生する方法についてです。
最初はそのままGIF画像を読み込めば表示できると思っていましたが、なかなかうまくいきませんでした。静止画であればすぐに表示できますが、動画として連続的に表示してくれずかなりハマりました。

まず、PILを使って以下のようなGIFファイルの内容を調べます。


from PIL import Image

img = Image.open("test.gif")

print("info: ",img.info)
print("animated: ",img.is_animated)
print("n_frames: ",img.n_frames)


そうすると、以下のような情報が出力されます。


info:  {'version': b'GIF89a', 'duration': 50, 'loop': 0, 'background': 0, 'extension': (b'NETSCAPE2.0', 803)}
animated:  True
n_frames:  100

このGIFファイルの場合はanimatedされているというのがわかります。またフレーム数が100あることもわかります。基本的にはこの二つの情報がわかれば大体大丈夫だと思います。
そのままこの二つの変数を使えば、異なるGIFファイルも読み込むことができると思います。


import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation
from IPython.display import HTML
from PIL import Image

img = Image.open(filepath+"-test.gif") # PIL Imageで読込み

# GIFの情報を調べる(今回はフレーム数だけ)
# print("info: ",img.info)
# print("animated: ",img.is_animated)
print("n_frames: ",img.n_frames)       # フレーム数を取得

fig = plt.figure()

# plt.axis('off') # 枠線なしの場合

ims = []    # 各フレームを格納するリスト

for frame in range(0,img.n_frames):  # フレーム数だけループ処理
    img.seek(frame)                  # 順番に各フレーム画像を探し出す
    im = plt.imshow(img)
    ims.append([im])

ani = animation.ArtistAnimation(fig,ims,interval=50)

plt.close()

HTML(ani.to_html5_video())
# HTML(ani.to_jshtml())


まず、内包されているフレーム数を調べるために、PIL Imageをインポートしておきます。img=Image.Open()したあとに、img.n_framesでフレーム数がわかります。
あとはそのフレーム数だけループ処理をし、ひとつずつリストに格納していきますが、img.seek()という関数を使うことで内包されている何番目の画像を取り出すか指示できます。
img.seek(0)なら最初の画像、img.seek(1)ならその次の画像というような感じで次のplt.imshow()に入れていきます。GIFの場合はこのimg.seek()を使うといいようです。
もし連番のPNG画像を複数扱うならば、img.seek()は使わずに、ファイル名の番号を手がかりに読込みと表示をさせればいいことになります。

順番にひとつずつ取り出してはimg.show()やplt.imshow()で表示させてもいいと思ったのですが、そうするとJupyter上では次々と複数の静止画が出てきて動画として見ることができませんでした。出しては消すという方法ももあるのかもしれませんが、今回はArtistAnimation()に取り込んで動画として表示させることにしました。


Jupyter上にマークダウンでGIFアニメを表示する:
先ほどのやり方は、GIFアニメを一枚ずつ取り出して、何らかの加工をするときには使えるテクニックかもしれませんが、ただ出来栄えを確認したり表示するだけなら、マークダウンでJupyter上のセルに表示させればいいと思います。


### Animated GIF
<div class="pull-left">
![Animated GIF](test.gif)
<div>

Jupyter上のセルをマークダウンに切り替えてこんな感じで記入すれば、普通にGIFアニメとして表示されるはずです。この場合は普通の静止画と同じ読込み方です。


追記:ipympl0.1.0から0.1.1へアップデート:
その後、Anacondaでモジュールのインストールをし直したりしていたら、Jupyter用のモジュールである「ipympl」が0.1.0から0.1.1にバージョンアップしていました。
conda install ipympl
でインストールできます。
さっそくアップデートしてみると、0.1.0ではエラーがでていた「%matplotlib ipympl」でしたが、動くようになりました。
「ipympl」にすると、上画像のように「nbagg」がバックエンドとして動くようです。「nbagg」なので従来のインタラクティブウィンドウが使えます。


「%matplotlib ipympl」によるバックエンド「nbagg」の場合:

Jupyter Notebook:
通常のインタラクティブウィンドウ:
表示されるが閉じても止まらない(一部バグ?)。セルのアウトプットクリアで止めるしかない。以前の「%matplotlib nbagg」のほうがいいかも。

JavascriptHTML出力:
動く。エラー無し。

HTML5 Video:
動く。エラー無し。


JupyterLab:
通常のインタラクティブウィンドウ:
表示されるが閉じることができない。表示できるようになったところまではいい。あともう少し。

JavascriptHTML出力:
絵が出ない。まだ使えない。

HTML5 Video:
動く。エラー無し。結局これが一番使える。しかもバックエンドが「inline」でなくても映った。


まとめ:
ということから、少し進展ありましたが、やはりHTML5 Videoで出力するのが一番安定しているようです。またアップデートくれば直るかも。
Jupyter Notebookのほうが安定しているので、いまはNotebookを使ったほうがいいのかもしれません。それほど不便というわけでもないので。

関連:(少し違う作り方)

2018年4月17日火曜日

JupyterLabを使ってみた(便利そうだけど、いまいち)

最近は、Jupyter Notebookを使っていますが、それと似たようなものでJupyterLabというのがAnaconda Navigatorに前からあったので試しにつかってみました。

Jupyter Notebookと同様にブラウザ上で使うことができるので便利。
さらにツリービュー、ターミナル、エディタもあり、画面分割可能でだいたい一通りのことができます。

Anaconda Navigatorからすぐにインストール(不要ならアンインストール)できるので試してみるにはよさそう。
かなりJupyter Notebookからは進化しているのですが、まだベータ版のようです。

右クリックでコピペできない:
ということでしばらく使っていると、やはりベータ版なのか多少不便な点が。
まず、マウス右クリックでコピペがメニューに出てこない。右クリックすると、以下のようなメニュー。

セル単位でのコピペは可能。しかし単語や行だけのコピペはできないために、仕方なくキーボードでCtrl+C、Ctrl+Vを入力するしかない。なんとメニューバーにすらない。

もうひとつとして、せっかくターミナルも同一画面上で使用できるのにもかかわらず、同様に右クリックだけではコピペできない。そのかわり、シフトキーを押しながらの右クリックならメニューが出てくる(以下)。

つい癖で右クリックしてみたものの、コピペできないということが少しストレスを感じてしまいます。改良されないのでしょうか?


%matplotlib nbaggが使えない:
Jupyter notebookではplt.plot()で表示するするとき、%matplotlib inlineを使いますが、アニメーションなどを表示するには、%matplotlib nbaggを使っています。しかしnotebookのほうでは使えるけれども、JupyterLabのほうだと、Javascriptエラーがでてしまいます。

JavaScript output is disabled in JupyterLab

しかし、このページにいくつか対処策があるようです。

conda install notebook jupyterlab # for updating jupyter notebook and lab
pip install ipympl
jupyter labextension install @jupyter-widgets/jupyterlab-manager jupyter-matplotlib

事前にnodejsとnpmをインストールしておき、notebookとjupyterlabをアップデート、さらにipympl、jupyter-matplotlibをインストールとのこと。
しかし、それでもなぜかエラーがでてしまいます。このページでも話されていることなのですが、まだダメっぽい。

また、%matplotlib nbaggを使っても、

こんな感じのエラーがでてしまいます。%matplotlib ipymplを使うともっとひどいエラーがでてしまい、jupyter-matplotlibがダメみたい。何回かインストールし直してみましたが改善されないので、今回は保留。

ということから、結構使い勝手はよくなっているJupyterLabなのですが惜しい。ベータ版なので仕方ないのかもしれませんが、しばらく待ってみることにします。
そもそも、Jupyter Notebookなら特に問題ないので、今まで通りJupyter Notebookでいいかと。

しかし、その後、jupyter-matplotlibのサイトをみると、

It requires matplotlib 2.0 or and ipywidgets 7.0 more recent.

と書いてあることからチェックしてみると、ipywidgetsのバージョンが低いことが判明。さっそくアップデートしてみました。
すると、いままで%matplotlib nbaggではだめだったgifアニメーション制作が可能になりました。
そして、上のようにバックエンドを確認してみたところ、いままではnbAggと出力されていたのが、'module://ipykernel.pylab.backend_inline'になっていました。少しだけ進展したのですが、あいかわらずアニメーション自体は表示されず、もう少しというところ。しばらくすれば、また進展あるかも。
要は、このようなアニメーションをJupyterLab上で表示させたいのですが、一コマごとに出力したpngファイルを順番に表示することができずにいるということです。このアニメーションはpngを束ねてgifアニメにしたものですが、それも表示できません。バージョンアップによって制作はできるようになったので、よかったのですが惜しい。もちろんJupyter Notebook上では表示も制作も問題ないというわけです。

その後もう少しやってみると、どうやらJupyterLabではmatplot.animationを使う場合、アニメーション出力はHTML(ani.to_html5_video())を使うと大丈夫でした。もうひとつHTML(ani.to_jshtml())というのもあるのですが、こちらはダメみたい。Jupyter Notebookではどちらとも大丈夫ですが、JupyterLabの場合はHTML(ani.to_html5_video())を使ったほうがよさそうです。

追記:
その後、GIFアニメーションの作成+読込み+再生について実験してみました


データ分析ツールJupyter入門
Posted with Amakuri at 2018.12.21
掌田津耶乃
秀和システム
販売価格 ¥3,024

2018年4月9日月曜日

Tensorflow 1.7へバージョンアップ

Tensorflowは少し前に1.6へバージョンアップしていたにもかかわらず、pytorchを使っていたため1.5のままでした。しかし1.6のあとすぐに1.7へバージョンアップしており、eager executionが改良された模様。



TensorflowはさすがGoogleだけあって開発スピードがはやいのか、内容に追いついていくのが大変。
Tensorflowといえば、

tf.placeholder()
tf.Variable()
tf.Session()

でコーディングしていくイメージでしたが、より使いやすくするために

tflearn : easy to learn, high level API
slim : a lightweight library for defining, training and evaluating models
estimator : a high-level TensorFlow API that simplifies ML programming
keras : a higher API wrapper
eager : define by run
lite : mobile, RPi

などが次々と導入されたため、どれがいいのかわかりにくくなってきたという感じ。
kerasが使えるようになってからは、シンプルなコーディングが可能になったけれども、pytorchのようなdefine by runタイプが便利になってきたこともあって、Tensorflowもeagerを導入してdefine by runも普通に使えるようになった感じ。

なんでもありとなってきたTensorflowという感じで、やはり機械学習フレームワークとしてはTensorflowがなんだかんだいって便利なのかもしれません。
ということで、使い勝手が改良されたTensorflowも再度使おうと、1.5から一気に1.7へバージョンアップしてみました。

GPUを使用しているので、

pip install --ignore-installed --upgrade https://storage.googleapis.com/tensorflow/linux/gpu/tensorflow_gpu-1.7.0-cp35-cp35m-linux_x86_64.whl

のpipコマンドでAnaconda仮想環境内にインストールしてあったtensorflow-gpu-1.5.0をアップグレード。
eager executionが便利そうなので、少し慣れておこうかと。

最近は、ベイズ推定などの確率的プログラミングに興味があって、

GPy
Pyro
Edward

などを勉強していまずがTensorflowにも、

Tensorflow Probability

というライブラリができたようで、Tensorflowだけであらゆることが可能となってきた感じ。



人気の投稿