matplotlibでアニメーションを作成，保存
一口にアニメーションといっても時間毎に図の更新がみたいだけの場合とその様子を動画に保存したい場合とがある。
閲覧用のアニメーション
plt.pause(interval)を用いることでノンブロッキングで現在の画像を更新していくことができる。
subplotsを用いた場合
ちょっと長くなるのが難点
import numpy as np import matplotlib.pyplot as plt # make data data = [[i, np.sin(i/50)] for i in range(100)] # prepare figure fig,ax = plt.subplots(1,1) ax.set_xlim(min([x[0] for x in data]),max([x[0] for x in data])) ax.set_ylim(min([x[1] for x in data]),max([x[1] for x in data])) xdata = [] ydata = [] line, = ax.plot(xdata,ydata) # iterate for dat in data: print(dat) xdata.append(dat[0]) ydata.append(dat[1]) line.set_data(xdata,ydata) plt.pause(0.1) # 10ms late
plt.figure()でやる方法
こっちのほうが楽。
import numpy as np import matplotlib.pyplot as plt # make data data = [[i, np.sin(i/50)] for i in range(100)] # prepare figure xdata = [] ydata = [] # iterate for dat in data: fig = plt.figure(1) xdata.append(dat[0]) ydata.append(dat[1]) plt.plot(xdata,ydata) plt.pause(0.1) # 10ms late plt.clf()
matplotlib.animationを使って図を保存する方法
matplotlibにはanimationという図を作成，保存するためのモジュールが存在する。
そのうち2つの手法が存在しイメージとしては以下のような感じである。
- ArtistAnimation：図をリスト形式で溜め込んで最後にまとめる。重い…？
- FuncAnimation：逐次的に図を作って追加していく。軽い…？
こちらはリアルタイムに描写というよりは，その後の保存に重点を置いている。一応plt.pause()とも連携できるはずだが，結構重いのでその辺は覚悟したほうが良いかも。
imagemagicのインストール
windows環境でやったのでLinuxの人は別の箇所を探して欲しいが，ffmpegとimagemagicをインストールすることでgifやmp4などの形式で動画を保存することができる。 実はimagemagicのインストール時にffmpegをインストールできるので実質imagemagicをインストールすれば良い。
以下のブログの手順に基づいてすすめていこう。 matplotlib のanimation を保存 - はしくれエンジニアもどきのメモ
１．imagemagicのインストール
こちらのダウンロードページに行き，Windows Binaryの一番上を選択してダウンロードすれば良い。 ダウンロード後はいろいろ聞かれるが基本的にYesで進めていけば良い。
２．matplotlibにimagemagicの場所を教える
お使いの環境のmatplotlibrcというファイルを操作する。ファイルの場所を確かめるには以下のコマンドを実行すると良い。
import matplotlib matplotlib.matplotlib_fname()
例えば以下のような感じの場所に置いてある。
'C:\\<Python path>\\lib\\site-packages\\matplotlib\\mpl-data\\matplotlibrc'
そしたらメモ帳などで開いて，# animation.convert_path:となっているところの＃をはずしてコメントアウトを解除し，
animation.convert_path: C:\Program Files\ImageMagick-7.0.1-Q16\magick.exe #みんなの環境だと多分違うぞ
のようにimagemagicのexeファイルのパスを与える。（エクスプローラから自分で調べること！）
ArtistAnimationで図を作成
ArtistAnimationの図の作り方は非常にわかりやすく，matlabのようにplot情報をリストに入れていき，最後にそれらをつなげて動画にする。 自分でもコードなど書いたが，こちらのブログが大変わかりやすく，自分で書くまでもないと思ったのでコードを借りて以下に表記する。どうもありがとうございます。
import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation fig = plt.figure() #figure objectを取得． ax = fig.add_subplot(111) x = np.arange(0,10,1) ims = [] #Line2D objectを入れるリストを用意 for time in range(x.shape[0]): im = ax.plot(x[0:time] ,x[0:time] ** 2) #Line2D objectを取得 ims.append(im) #imsにappendする #ArtistAnimation機能で，imsの中の画像を繋ぎ合わせる． ani = animation.ArtistAnimation(fig, ims, interval = 100) #imagemagickを使って，gif画像を保存 ani.save("test.gif", writer="imagemagick")
FuncAnimationで動画を作成保存する
こちらは，Animationの生成規則を与えて，逐次これを計算させるものになる。先程と同様の方がわかりやすい解説を上げていてくださり，まずはこちらを試してみると良いだろう。
以下に軽く私の理解を述べておく。URL先のコードで以下のように呼ばれているFuncAnimationでは，
- 最初にキャンバスとなるfigureのハンドル
figをわたし - 逐次処理を実行する
update()関数を渡す - updateに渡す引数を決める
fargsというのがオプションで存在し， - 画像間のinterval時間
- 最大繰り返しフレーム数
framesを渡す
ani = animation.FuncAnimation(fig, update, fargs = (x, y, ax1, ax2), interval = 100, frames = 10)
という風になっている。従ってframeが0から9に至るまで，すなわち，
for i in range(frame): update(i,args)
という感じになっているはず。
この他にもinit()という初期化関数を渡したりできる。
自分Verとしては例えば以下のようなものになる。ただし，ランダムな点をPlotしようとするので挙動はきっとクソ。 （Originalではちゃんとplotすべきデータがあったが簡略化させていただいた）
fig, ax = plt.subplots() line, = ax.plot([], [], 'r',label='trajectory') line2, = ax.plot([], [], 'kx',label='current position') ax.grid() xdata, ydata = [], [] plt.xlabel('x [m]') plt.ylabel('y [m]') plt.grid ='on' #def data_gen(): data = [[x[0],x[1]] for x in np.random.rand(100,2)] def init(): ax.set_ylim(-0.2, 1.2) ax.set_xlim(-1.2, 0.2) line.set_data(xdata, ydata) #line2.set_data(xdata, ydata) return line, def run(i): # update the data x, y = data[i] xdata.append(x) ydata.append(y) line.set_data(xdata, ydata) line2.set_data(x,y) plt.legend(loc='upper left') ani = FuncAnimation(fig, run, blit=False, interval=10, frames = int(len(data))-1, repeat=False, init_func=init) ani.save("plot2.gif",writer='imagemagick') plt.show()
Cannon Power Shot G9X MarkⅡで月を撮ってみた
最近手持ちのノートPCが死んで個人ファイル以外初期化しました。 ドキュメントやピクチャは消えませんがProgramFilesやAppDataが消えたので結構不便してます。
Cannon Power Shot G9X MarkⅡについて
訳あってこのコンデジを家で保管しています。
Canon コンパクトデジタルカメラ PowerShot G9 X Mark II ブラック 1.0型センサー/F2.0レンズ/光学3倍ズーム PSG9XMARKIIBK
- 発売日: 2017/02/23
- メディア: エレクトロニクス
この機に少しコンパクトデジタルカメラについて調べたので軽く紹介すると、
- 「一型センサ」と呼ばれる受光センサがコンデジとしては大きい部類
- 電池を入れて218gと非常に軽量
というのが魅力的な点だと思います。
センササイズの大きさはそのまま画質に影響する1とのことで、夜景などを撮るときに如実に差が出るそうな。 何よりSonyのRX100シリーズより安めなのがよいです。
月を撮った
光学ズームは3倍までと弱いのですが撮れるらしいとのことで、以下のブログで月を撮った人のパラメータを参考にしました。 Canon PowerShot G7Xで、月面の撮影に挑戦：南の島旅
撮影条件
最初は手でもって試しましたが手振れの影響がひどいのでTHETAを買った時の三脚を使いました。そのままだとシャッターボタンを押したときの振動で三脚が振動してぶれたのでタイマー機能を使って振動を抑えてから写真を撮っています。
パラメータは
謎現象でフォーカスがなぜか勝手に近いところによるためマニュアルで固定しました。
結果
撮ってきた画像を拡大してトリミングしました。
確かに月の模様が見えています。
拡大するとちょっと荒いのがばれてしまいますが初めてでもこれくらいのは撮れるようですね。
なお、個人的興味でwaifuで高精細化できるか試してみましたがローパスがかかったような画像になってしまいました。
やはり複数撮って合成するしかないようですね。
EEIC後期実験演習例 Python Control版
前にjupyterで書いてみたは良いが重くてGithubで公開をためらった代物です。 Markdown変換してみました。
はじめに：Environment settings
Using anaconda you just need to activate following command.
pip install python-control
For further information see here. folllowing command is additional.
- Textbook for this experiment is here.
# !python -m pip install --upgrade pip # !pip install slycot # I failed to install these command
Initialization
Import matlab like package from control.matlab
from control.matlab import * # MATLAB-like functions
basic feature
At first you can use tf function like matlab.
- 下記は演習2の(1)の問題を例に取り上げます。
# Q2-1 num1 = [1,2] den1 = [3,1] sys = tf(num1,den1) print(sys) # this is the system output
s + 2
-------
3 s + 1
Load matplotlib for drawing figure
import matplotlib.pyplot as plt %matplotlib nbagg # import scipy from scipy import arange
# step responce # Generate step responses plt.figure(1) (y1a, T1a)=step(sys,T = arange(0, 10, 0.01)) plt.plot(T1a,y1a)
<IPython.core.display.Javascript object>
[<matplotlib.lines.Line2D at 0x15ae2ea4438>]
Error メモ
sysの関数を変更すると？ときどき'tuple' object is not callableというエラーが起きます。
イマイチ自信がないですがどこか変更してはいけない部分を変更しようとしているようです。
一度カーネルをRestartすると治るのでとりあえずはこれで勘弁してください。
# bode plot
bode1 = bode(sys)
<IPython.core.display.Javascript object>
演習4：2次系の応答
数式も打ち込みたいけど時間がないので割愛。 減衰比（ダンピング）を変えたらどうなる？という話です。
共振周波数を変えるVerもありましたが，ダンピングの例だけ示します。
Gs = [] K = 1 w_n = 10 Ys = [] Ts =[] legend = [] # changing zeta plt.figure(3) for z in range(0,8): zeta = z * 0.2 Gs.append(tf([K*w_n*w_n],[1,2*zeta*w_n,w_n*w_n])) (y1, T1)=step(Gs[-1],T = arange(0, 1, 0.001)) Ys.append(y1) Ts.append(T1) # plot plt.plot(Ts[-1],Ys[-1]) legend.append('\zeta is '+str(zeta)) plt.legend(legend) plt.show()
<IPython.core.display.Javascript object>
- どうせならbode線図も見てみましょう
plt.figure(4) for G in Gs: bode(G)
<IPython.core.display.Javascript object>
演習5：ゼロ点の性質
最終値の定理から行き着く値は1/15で0.06くらいになることはわかりきっている。 ではゼロ点はどこに寄与するのか？
よしんば過渡応答であったとしてどのように？
という設問であるが，私が思うに分子がs+1になるような極端な例がないとこの設問の意図は組みにくいであろう。
早急に実験レポートに反映されるべきである。
# sys5 = [tf([1],[1,8,15]), tf([1,4],[4,32,60]) , tf([1,50],[50,400,750]), tf([1,1],[1,8,15])] plt.figure(5) for sys in sys5: (y1, T1)=step(sys,T = arange(0, 1, 0.001)) plt.plot(T1,y1) plt.legend(['G1','G2','G3','Additional']) plt.show()
<IPython.core.display.Javascript object>
レポートの解答で無難に言うならば，「遅いゼロ点はオーバーシュートを引き起こす。早い零点はまぁ無視しても良い。」ということになるが，これはどうしてか。 せっかくなのでBode線図も見てみよう
plt.figure(6) for sys in sys5: print(sys) bode(sys) plt.legend(['G1','G2','G3','Additional'])
<IPython.core.display.Javascript object>
1
--------------
s^2 + 8 s + 15
s + 4
-----------------
4 s^2 + 32 s + 60
s + 50
--------------------
50 s^2 + 400 s + 750
s + 1
--------------
s^2 + 8 s + 15
<matplotlib.legend.Legend at 0x15ae3a986a0>
遅いゼロ点ではゲインが不用意に盛り上がっていることが確認できる。 不要な周波数成分を増幅することにより変な振動が起きてしまう → オーバーシュートが出てくる といった感じである。
直感的に説明すると，「分子が1の場合，安定極であれば分母は単調増加して高周波ほどゲインが下がるのに対し，分子にsの項が加わると分子も増加する要素が加わるためsの寄与が相対的に高い遅い零点ではある周波数周辺でゲインが盛り上がる → オーバーシュート」という感じであろうか。
演習８
結構難しい気がする。問題の意図が汲めない… 一応低域の位相を遅らせることと引き換えに高域の位相を少し早めているって言うことなんだろうが。誰か解説手伝って。
余談：zpk記述
因数分解されたシステムは記述が面倒くさい。 zpk2tfというscipyの関数を用いることでゼロ点，極，DCゲインの３つからtfを生成可能だが，生成結果はtupleであり２変数を受け取るtfに直接食わせられない。
そんなときは*を用いた展開表現を用いることで，tupleを複数argumentとして認識させることが出来る。
Gp =tf(*zpk2tf([],[0,-2,-10],40)) Gc =tf(*zpk2tf([-0.2,-1],[-0.02,-5],5)) plt.figure(7) bode(Gp) bode(Gp*Gc) bode(feedback(Gp*Gc)) plt.legend(['Plant','Plant+C', 'Feedback']) plt.show()
<IPython.core.display.Javascript object>
根軌跡
rlocusやってrlockfindがない！ってOctaveで毎年なるんですがpythonにもなかった…
tau = 0.01 plant = tf([1],[1,3,2,0]) gc = tf([-1.5],[-1.0/tau]) x1,y1 = rlocus(plant*gc)
<IPython.core.display.Javascript object>
Saving
データとか保存用になにかする（予定）
MATLAB＆YALMIPでLMIを解けるようになるまでの手順
Software Installation
ソルバはYALMIP以外にCVXgen等もありますが，どっちがいいんですかね？議論ありますが有識者の意見を聞きたいです。
Yalmipまわりのインストールは以下のスライドが比較的参考になります。2018年9月の時点ですでに情報は古いですがURL等は参考になるでしょう。 http://www.maizuru-ct.ac.jp/control/kawata/iscie/doc/install_yalmip_new.pdf
SeDuMi
単にダウンロード MATLAB直下に解凍。
SDP3
SDPT3 Ver3-4をダウンロード。 MATLAB直下に解凍。
その後，解凍したディレクトリに移動して以下のコマンドを実行します。
Installmex(1)
YALMIP
実はYALMIP単体でもそれなりにソルバはありますが，SDP（半正定値計画問題）を解くには上記のソルバが必要です。
addpath(genpath('C:\Users\<ユーザ名>\Documents\MATLAB\SeDuMi_1_3'))
addpath(genpath('C:\Users\<ユーザ名>\Documents\MATLAB\SDpt3-4.0'))
addpath(genpath('C:\Users\<ユーザ名>\Documents\MATLAB\YALMIP-master'))
yalmiptest
実行結果は以下のとおりです。とりあえずSDPは解けているようなのでOKです。
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | Test| Solution| Solver message| +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | Core functionalities| N/A| Successfully solved (YALMIP)| | LP| Correct| Successfully solved (LINPROG)| | LP| Correct| Successfully solved (LINPROG)| | QP| Correct| Successfully solved (QUADPROG)| | QP| Correct| Successfully solved (QUADPROG)| | SOCP| Correct| Successfully solved (SeDuMi-1.3)| | SOCP| Correct| Successfully solved (SeDuMi-1.3)| | SOCP| Correct| Successfully solved (SeDuMi-1.3)| | SDP| Correct| Successfully solved (SeDuMi-1.3)| | SDP| Correct| Successfully solved (SeDuMi-1.3)| | SDP| Correct| Successfully solved (SeDuMi-1.3)| | SDP| Correct| Successfully solved (SeDuMi-1.3)| | MAXDET| Correct| Successfully solved (SeDuMi-1.3)| | MAXDET| Correct| Successfully solved (SeDuMi-1.3)| | Infeasible LP| N/A| Infeasible problem (LINPROG)| | Infeasible QP| N/A| Infeasible problem (QUADPROG)| | Infeasible SDP| N/A| Infeasible problem (SeDuMi-1.3)| | Moment relaxation| Correct| Successfully solved (SeDuMi-1.3)| | Sum-of-squares| Correct| Successfully solved (SeDuMi-1.3)| | Bilinear SDP| N/A| No suitable solver| +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
LMIの試験プログラム
状態FBでのレギュレーター設計問題を考えます。 LMIでは極配置では出来ない入出力制約等も考えられます。（解があるとは言っていない）
今回はただ単に極が-αより小さいという条件を設定します。
式の解説
ちょっと書くのが面倒なのでしばらくは（省略）。 2次形式のリアプノフ関数の安定化のための条件を解きます。
以下の教科書4章p80~を参照のこと。
LMIによるシステム制御 - ロバスト制御系設計のための体系的アプローチ
- 作者: 蛯原義雄
- 出版社/メーカー: 森北出版
- 発売日: 2012/03/08
- メディア: 単行本（ソフトカバー）
- クリック: 1回
- この商品を含むブログを見る
プログラム
-パラメータ設定（状態方程式）
k = 0.2; f = 0.01; A = [ 0 0 1 0; 0 0 0 1; -k k -f f; k -k f -f]; B = [ 0; 0; 1; 0];
- LMI 立式
αは0.1とします。あんまり大きくすると解けなくなるっぽい。極をある一定以上早くなるよう設定した上でγでフィードバックゲインKをなるべく小さくなるように縛ります。
X = sdpvar(4,4,'symmetric'); Y = sdpvar(1,4,'full'); gamma = sdpvar(1); LMI = [X>eye(4), A*X+B*Y+(A*X+B*Y).'+ 2*alpha*X<0,[gamma,Y;Y.', eye(4)]>0]; sol = solvesdp(LMI,gamma) X = double(X); Y = double(Y); gamma = double(gamma); K = Y/X; % LMI2 result pole = eig(A+B*K); display('Solved LMI pole') display(pole);
配置結果
ゲインを上げずにという制約があるので極の制約のギリギリを攻めるのは予想が付きます。 はいこんな感じ。
pole = -0.1000 + 0.6672i -0.1000 - 0.6672i -0.1000 + 0.0979i -0.1000 - 0.0979i
X>eyeとしているのは安全のためなのでこれを本来のX>0としてしまえばもう少しゆるい条件で解けそうと思っているのですがどうでしょう？。
やってみた所感
とりあえず動作を確認したまでで，LMIの利点を活かした設計などはまだ触れられていません。 ここまでのプログラムで一応ロバスト安定化制御器などは導けますが。
例えばalphaを3とかにするだけでInfeasibleになります。
あと，Successfully Solvedと出ても解けてなかったりして意外といろいろ気を配ります。
そもそもソルバに全部投げているのでCVXgenとかだと解けたりするんでしょうか？ イマイチわかりきっていないところがあります。
Latex 表組みのサイズを微調整
表のサイズの微調整
従来，Latexの表が大きすぎた時 \smallなどのコマンドで表の文字を小さくしてサイズを収めるといったことをしていました。 しかしこの調整は非常に表を醜くします。
正しくは\resizebox{}環境を用いるべきです。
\columnwidthで横幅ぴったりに合わせます。2カラムぶち抜くときは\textwidthを用いると良いです。
% \usepackage{graphics} \resizebox{\columnwidth}{!}{% %table 環境を中に }
GUIに依る表の作成
Latexの表の作成は最もだるい事の一つですが，Tablegeneratorを用いることで簡略化出来ます。 Excelとかで表を作ってしまってGUIで整形，その後にTexに変換して貼り付けるのが最近のやり方です。
お試しあれ。
Google ColaboratryでGPU（TPU）付の機械学習環境が使えるらしい
趣味人としてはそろそろ独自のコードを書いて公開したいところ。（なお本業…） [:contents]
Google Colaboratry事始め
どうも今年の初めにGoogleがCloudで機械学習環境を提供しはじめたようで，それに関する記事が林立しています。 流行りには乗り遅れましたが，気になったのでフォローしてみます。
引用
以下の記事を基にしてテストしました。図なども借りています。
[1] Google Colaboratoryの無料GPU環境を使ってみた – MISO
[2] 【秒速で無料GPUを使う】深層学習実践Tips on Colaboratory
仕様
文献[2]によると重要そうな仕様は以下のとおり。 Ubuntu環境です。
- n1-highmem-2 instance
- Ubuntu 17.10
- 2vCPU @ 2.2GHz
- 13GB RAM
- 40GB Free Space
- GPU NVIDIA Tesla K80
- アイドル状態が90分続くと停止
- 連続使用は最大12時間
- Notebookサイズは最大20MB
Teslaというやつは浮動小数点に強い大規模計算機用のより強いGeforceみたいなやつです。単純にRAMだけでも1080Tiの倍以上の24Gあります。 容量が40Gという点以外は個人が趣味で買えるスペックは大幅に超えていますね。
Jupyterで開発をするという仕様上，コマンドライン操作は!をつけて実行します。
!pip freeze
等でパッケージを確認できますが，TensorflowやKeras，matplotlibなんかは当然のように入っています。環境構築も不要です。
Ubuntuベースにしているのでaptgetも使えるのですが今日は試しませんでした。
準備することとか
以下のURLにGoogleアカウントを持った状態でアクセスすれば環境が手に入ります。 https://colab.research.google.com/notebooks/welcome.ipynb
「ファイル」からPython3のノートを作成して実行していきましょう。 最初に上げたURLの記事からたどればよいでしょう。
CPU → GPU環境への移行
「ランタイム」の欄から「ランタイムのタイプ」を開き，図のように選択します。
今はTPUも選択できるそうな。
作成したファイルの保存
ファイルの保存方法は
- ダウンロード
- Githubへ保存
- GoogleDriveへ保存
がありました。Githubに保存するときはポップアップが出るので事前にポップアップブロックを無効にしておいてください。
KerasでのMNISTのテストと所感
KerasでのMNISTサンプルを試してみました。 ノートブックを上げておきます。 https://colab.research.google.com/gist/YoshiRi/b7a984902222a14d3f4fba31625060ca/kerastest.ipynb
6万Stepを12Epocやって自前のGeforce1050iを載せたPCで計算した時は13s/Epocくらいでしたが，この環境では11s/Epocなのでちょっと早くなっていますね。コードの書き方で変わるんでしょうか。
[http://ossyaritoori.hatenablog.com/entry/2018/03/27/Keras%26_Tensorflow%28GPU%E6%9C%89%29%E3%81%AE%E7%92%B0%E5%A2%83%E6%A7%8B%E7%AF%89_on_Windows_with_Anaconda:embed:cite]
結果
このあたりはサンプルプログラムなのでどうせほぼ同じですが一応。
x_train shape: (60000, 28, 28, 1) 60000 train samples 10000 test samples Test loss: 0.02876817006157935 Test accuracy: 0.9913
所感
GPUマシンを持っていなくても（環境構築しなくても）機械学習の学習が出来るのは初学者にはかなり良いと思います。 特に環境構築や機材は初学者にとって大きな壁なので回線さえあればきちんとした環境で実行できるというのは教育効果が非常に高そうです。
とりあえず，コードを試して満足したい方はここから初めて良いと思いますが，一方で今後ちゃんとやっていくならどこかで環境構築はしておくべきだと思います。 結局自分でデータ集めてモデル組まないと意味ないので僕もこんなサンプルコードばかり実行してないで自分のやりたいプロジェクトを進めていこうとおもいます。
MATLAB背景色変更
Texstudioの背景色変更に加え，matlabも背景色変更します。 ossyaritoori.hatenablog.com
ここで配ってます。 github.com
windows環境では以下のコマンドをmatlabのホームで実行しておしまい。私はdarksteel推しです。
!git clone https://github.com/scottclowe/matlab-schemer.git
cd matlab-schemer\
schemer_import('schemes\darksteel.prf')
