02. グラフを描こう

Photo by Chris Liverani on Unsplash

プログラミング言語について語るとき,通常は「文法が…」「この言語の新しい機能は…」とやるところだが,そういう流れだと教わる方はかなり退屈だし,なかなかピンとこないだろう. そこで,まずはデータの可視化手段の一つである「グラフを描く」話をしよう.

実際に計算をしてもその結果はグラフなどで可視化しないと理解できないことがほとんどなわけだから,そういう意味でもこの話は重要だしね.

さて,Julia でグラフを描くツール, ライブラリ(package)はたくさんあるが、現在デフォルト扱いになっているものが "Plots" package だ. インストール済みでなければ, Julia 環境で

1
2
using Pkg
Pkg.add("Plots")

とすればインストールできる. 大きめの Package なのでインストールにはやや時間がかかるが,まあ待とう.

いったんインストールしてしまえば,その後は

1
using Plots

とすればそのセッション中では Plots package が有効になり(後述の) plot 命令などが使える.

実は,上のように using Plots とするだけではなく

1
2
using Plots
plotly()

とすると, 以降の同じコマンドでより便利なプロットが出力される.マウスで拡大縮小とか,3次元グラフだとマウスで回転もできる. ただこの機能は web browser との相性問題がある上,少し遅い(特に最後の gif 作成は大変遅い)ので,適宜使い分けよう.

さて,では、この Package をインストールしていない人はインストールしよう. その後,この「Plots を使うよ宣言」(using Plots)をしておこう.

そして、Julia で基本的なグラフを描くには、だいたいは plot 命令を使えば良い. 初歩的なところだとおおよそ、

  • 関数を plot に渡すと、普通の 2次元グラフを描いてくれる
  • ベクトルを plot に渡すと、普通の 2次元グラフを描いてくれる
  • 行列を contour に渡すと(3次元データの)等高線グラフを、
  • 行列を heatmap に渡すと(3次元データの)色付けによる z軸上方向から見たグラフを、
  • 行列を wireframe に渡すと(3次元データの)立体ワイヤーフレームグラフを、
  • 行列を surface に渡すと(3次元データの)立体のきれいなグラフを、描いてくれる.

あたりを押さえておけば良い感じだろう.

プログラミングについて知っている人向けに言うと,Julia でベクトルといえば 1次元配列そのもののことで,行列は2次元配列そのものである.
上の contour(A)plot(A, st = :contour) の省略形だ.heatmap 等も同様.

さっそく,例えば三角関数のグラフを描いてみよう.

1
plot(sin)

そのセッションで初めて使う plot は実行まで少し時間がかかる.待とう.

とても簡単だね!!

次に,例えば、データぽいベクトルを作ってみる.

1
2
v = [ n^2 for n in -2:0.1:2 ]  
# -2から2まで0.1おきの数字を二乗して並べたベクトルだね.
41-elementArray{Float64,1}:
4.0
3.61
3.24
2.8899999999999997
2.5600000000000005
2.25
1.9599999999999997
1.6900000000000002
1.44
1.2100000000000002
1.0
0.81
0.6400000000000001
⋮
0.81
1.0
1.2100000000000002
1.44
1.6900000000000002
1.9599999999999997
2.25
2.5600000000000005
2.8899999999999997
3.24
3.61
4.0

このベクトル vplot 命令に渡してみよう

1
plot(v)

無事にグラフが描かれることがわかってもらえると思う.
基本はこのノリで、データをベクトル形式に突っ込んでそれを plot に渡せば良い.

グラフの各点にマーカーがあったほうが良いな、という場合は例えば次のようにすれば良い.

1
plot(v, marker = :auto)

ちなみにこうした情報は、Julia 自身のもつ内部マニュアルでひいていくことができる(英語だけどね). 今回は次のような感じでわかる.

1
?plot  # plot 命令について、尋ねてみる
The main plot command. Use plot to create a new plot object, and plot! to add to an existing one:
  plot(args...; kw...)                  # creates a new plot window, and sets it to be the current
  plot!(args...; kw...)                 # adds to the `current`
  plot!(plotobj, args...; kw...)        # adds to the plot `plotobj`

There are lots of ways to pass in data, and lots of keyword arguments... just try it and it will likely work as expected.
When you pass in matrices, it splits by columns. To see the list of available attributes, use the plotattr([attr]) 
function, where attr is the symbol :Series:, :Subplot:, :Plot or :Axis. Pass any attribute to plotattr as a String to 
look up its docstring; e.g. plotattr("seriestype").
1
2
plotattr(:Series)  
# 最初にでてきたキーワード :Series について、上の通りに調べてみる.
Defined Series attributes are:
arrow, bar_edges, bar_position, bar_width, bins, contour_labels, contours, fill_z, fillalpha, fillcolor, 
fillrange, group, hover, label, levels, line_z, linealpha, linecolor, linestyle, linewidth, marker_z, 
markeralpha, markercolor, markershape, markersize, markerstrokealpha, markerstrokecolor, markerstrokestyle, 
markerstrokewidth, match_dimensions, normalize, orientation, primary, quiver, ribbon, series_annotations, 
seriesalpha, seriescolor, seriestype, smooth, stride, subplot, weights, x, xerror, y, yerror, z
1
plotattr("markershape")  # マーカーの形状についてなにかわかりそうだ
markershape {Symbol, Shape, or AbstractVector}
markershapes, shape
  Choose from Symbol[:none, :auto, :circle, :rect, :star5, :diamond, :hexagon, :cross, :xcross, :utriangle, 
  :dtriangle, :rtriangle, :ltriangle, :pentagon, :heptagon, :octagon, :star4, :star6, :star7, :star8, 
  :vline, :hline, :+, :x].
Series attribute,  default: none

これによると、マーカーの形状を自分で指定することも簡単だ.例えば六角形でやってみると…

1
plot(v, marker = :hexagon)

よく見ると確かに六角形がマーカーになっている.

次に、3次元グラフにもチャレンジしてみよう

1
2
A = [ x^2 + y^2 for x in -1:0.1:1, y in -1:0.2:1 ]  
# x^2 + y^2 の値に相当する行列を作る
21×11 Array{Float64,2}:
2.0   1.64  1.36  1.16  1.04  1.0   1.04  1.16  1.36  1.64  2.0
1.81  1.45  1.17  0.97  0.85  0.81  0.85  0.97  1.17  1.45  1.81
1.64  1.28  1.0   0.8   0.68  0.64  0.68  0.8   1.0   1.28  1.64
1.49  1.13  0.85  0.65  0.53  0.49  0.53  0.65  0.85  1.13  1.49
1.36  1.0   0.72  0.52  0.4   0.36  0.4   0.52  0.72  1.0   1.36
1.25  0.89  0.61  0.41  0.29  0.25  0.29  0.41  0.61  0.89  1.25
1.16  0.8   0.52  0.32  0.2   0.16  0.2   0.32  0.52  0.8   1.16
1.09  0.73  0.45  0.25  0.13  0.09  0.13  0.25  0.45  0.73  1.09
1.04  0.68  0.4   0.2   0.08  0.04  0.08  0.2   0.4   0.68  1.04
1.01  0.65  0.37  0.17  0.05  0.01  0.05  0.17  0.37  0.65  1.01
1.0   0.64  0.36  0.16  0.04  0.0   0.04  0.16  0.36  0.64  1.0
1.01  0.65  0.37  0.17  0.05  0.01  0.05  0.17  0.37  0.65  1.01
1.04  0.68  0.4   0.2   0.08  0.04  0.08  0.2   0.4   0.68  1.04
1.09  0.73  0.45  0.25  0.13  0.09  0.13  0.25  0.45  0.73  1.09
1.16  0.8   0.52  0.32  0.2   0.16  0.2   0.32  0.52  0.8   1.16
1.25  0.89  0.61  0.41  0.29  0.25  0.29  0.41  0.61  0.89  1.25
1.36  1.0   0.72  0.52  0.4   0.36  0.4   0.52  0.72  1.0   1.36
1.49  1.13  0.85  0.65  0.53  0.49  0.53  0.65  0.85  1.13  1.49
1.64  1.28  1.0   0.8   0.68  0.64  0.68  0.8   1.0   1.28  1.64
1.81  1.45  1.17  0.97  0.85  0.81  0.85  0.97  1.17  1.45  1.81
2.0   1.64  1.36  1.16  1.04  1.0   1.04  1.16  1.36  1.64  2.0

さあ、プロットしてみよう.最初は contour (等高線表示になる)だ.

1
contour(A)

ちなみに、plot に行列 A だけを渡すと、(この場合は)「縦ベクトル11個のデータを一度に渡した」と解釈され、11個のグラフを同時描画する. これはこれで使えるテクニックなので覚えておこう.

1
plot(A)

次は heatmap だ.

1
heatmap(A)

heatmap

まあどちらかというと、contour に fill オプションをつけたほうが良さそうだ.

1
contour(A, fill = true )

次は wireframe だ.

web browser への負荷は高いが…

1
wireframe(A)

wireframe

最後は sufrace だ.

1
surface(A)

surface

wireframesurface などの 3次元プロットでは、camera = (平面角, 高さ角) というオプションをつけると、その位置にカメラを置いた表示にしてくれる. デフォルトは camera = (30,30) だ. ちなみに、平面角も高さ角も (なぜか) 0 から90までに制限されているようだ. 詳しくは plotattr("camera") とすると簡単なマニュアルが出るので読んでみよう.

ちなみに、下の画像はこの camera オプションを使ってみた例だ. デフォルトのものと、グラフを眺めている位置が異なるのが分かるだろう.

1
wireframe(A, camera=(22,58))

wireframe

plot したグラフを保存しよう

役に立つグラフが得られたならば保存しておきたい. plot したグラフの保存はいくつかの方法があるが,一番簡単な方法はグラフを描いた後,savefig コマンドを使う方法である. それを紹介しておこう.

例えば,適当な4次関数を以下のように描いたとしよう.

1
plot( x -> x^4 - 20*x^2 + 1 )

このグラフを画像ファイルに保存したければ,この後に,

1
2
3
4
5
6
plot!(size = (1200,800)) 
# 画像ファイルのサイズ指定.指定しない場合は 600x400 になる.

savefig("sample.png")
# これで sample.png という名前で png ファイルとして保存される.
# 拡張子で画像ファイルの種類が決まる.

などとすれば良い. jupyter のファイラー画面などから,確かに画像ファイルが保存されていることを確認しておこう.

savefig コマンドで出力できる画像ファイルの種類は描画 backend (実際に描画を行うライブラリ)の種類によって異なる. とくに設定を変えずに Plots を使っているデフォルトだと GR と呼ばれる backend が用いられ,pdf, png, ps, svg が出力可能だ. 詳しくは Plots マニュアルの Supported outpout file formats を見よう.

パラメータを変えて試行錯誤するときに便利なパッケージ Interact を使って、グラフを動かそう

なかなか便利な package Interact をここで紹介しよう. 詳しくはまた解説する機会を設けるが,まずは簡単に,という感じだ.

この package のインストール方法は環境によって異なる. ここでは Julia + jupyter 環境の話をする. 他の環境にインストールしたい人は 開発者のGitHub を見よう.

まず,できれば Julia の REPL 環境(黒い文字だけの環境)で(まあ,たぶんjupyter 環境上でも大丈夫だが),

1
2
3
4
5
6
using Pkg
Pkg.add("Interact")
Pkg.add("WebIO")

using WebIO
WebIO.install_jupyter_nbextension()

として 2つの package をインストールするとともに jupyter 用の拡張機能をインストールする. ちなみに,jupyter 上でこの作業を行うと上の 5行目の後で「WebIO が動かない」旨のエラー表示が出るが気にしないで良い.そのための拡張機能インストールが 6行目でなされるので.

その後,この Julia + jupyter 環境をいったん終了し,再び Julia + jupyter を立ち上げよう. こうしないとインストールした拡張機能がうまく動かない.

これでインストールは終了だ.

ちなみに,jupyter 用の拡張機能は上でインストールした WebIO 対応のもの以外にもいろいろある. 興味がある人は Anaconda3 のプロンプト(jupyter が動く時点で既にインストールされているだろうからメニューから探してみよう)を起動し,そこで

1
conda install -c conda-forge jupyter_contrib_nbextensions

とすると,「典型的な拡張機能のインストール & 設定画面表示がONに設定」がなされる. jupyter のファイラー画面の上方のタブを見てみよう."Nbextensions" というタブが増えていて,そこをクリックすると拡張機能の一覧が出るはずだ. あとは触って解説を読むなどしてみると良い.ただし,python でしか有効でなく Julia では動かない機能もあったりするので,解説はよく読むこと.

さて,話を戻そう. このようにして上記のインストールが済んだ状態で,この package を読み込む.

1
using Interact

これでこのセッション中は Interact package が使えるわけだ. 試しに,グラフと関係ない次の簡単な例でこの Interact package の機能を味わってみよう.

1
2
3
@manipulate for x in 0:10
  x^2
end

さて,この3行を実行すると という感じの表示が出て,バーをマウスで掴んで動かすと変数 $x$ の値を変えられ,出力 $x^2$ の値もそれに応じて変わるはずだ.

このように,プログラムの中の変数/パラメータの数字をマウスで自在に変えられ,その結果をリアルタイムでみることができるのがこの @manipulate 命令のありがたい機能なのだ.

@manipulate 命令がうまく働かない場合:
jupyter の拡張機能に少し不具合があるようで, @manipulate 命令がきちんと動かないことがある1(バーをマウスで掴んで動かしても $x^2$ の表記が変わらないなど). そうした場合,Web browser でページ再読み込みをし,再度同じセルの @manipulate 命令を実行すれば良い. 例えば firefox だと下図の赤丸の部分のボタンを押して再読み込みしてから,同じセルの @manipulate 命令を実行すれば良い.

再読み込みを行うと Interact package などで出力した表記は消えてしまうが,再読み込みを行う前からの Julia の計算は背後で続いたままなので,そのまま計算や作業の続きができる.

さて話をグラフ描画に戻して, 次の内容を実行してみよう(最初の2行はもう見たことがあるだろう)..

1
2
3
4
5
6
using Plots
A = [ x^2 + y^2 for x in -1:0.1:1, y in -1:0.2:1 ]  

@manipulate for θ in 0:2:90, h in 0:2:90
  wireframe(A, camera=(θ,h))
end

すると下のように,バーが 2本とワイヤーフレームのグラフが表示される. バーを掴んで動かしてみよう.「このグラフを見る角度」を変えられるはずだ.便利だね!

bars

wireframe

簡易アニメーション

まず,上の camera オプションに絡めて、簡単なアニメーションをその場で作る方法も書いておこう(本格的なものについてはまた教えよう). これは上の @manipulate にそっくりの使い方で、下記の例をやってみればわかるだろう.

1
2
3
@gif for θ in 0:1:120
    wireframe(A, camera=(θ,30))
end

とすると、

┌ Info: Saved animation to 
│   fn = c:\home\julia-programs\v1.5\tmp.gif
└ @ Plots c:\julia\PKG\v1.5\packages\Plots\D7Ica\src\animation.jl:104

と出力されて、"tmp.gif" という名前でアニメーション gif ファイルが作られたことがわかる. そこで、web browser でさきほどのファイラー画面に行き、このファイルがないか見てみよう. 見つけたら、あとはそのファイルを例えばクリックするだけで、browser がその動画を見せてくれるだろう.

gif-sample

レポート

下記要領でレポートを出してみよう.

  • e-mail にて,
  • 題名を 2020-numerical-analysis-report-02 として,
  • 教官宛(アドレスは web の "TOP" を見よう)に,
  • 自分の学籍番号と名前を必ず書き込んで,
  • 内容はテキストでも良いし,pdf などの電子ファイルを添付しても良いので,

下記の内容を実行して,結果や解析,感想等をレポートとして提出しよう.

  1. 授業用 web に出てくる例と異なる 2次元グラフ,3次元グラフを plot し,そのデータの解説とグラフそのものをレポートにて提出しよう.
  2. 授業用 web に出てくる例と異なる gif グラフを作成し,レポートに添付して提出しよう.

  1. より正確に書くと,「最初の @manipulate 命令が拡張機能を起動し,ブラウザが?それとの連携を試みる.このタイミングが悪いなどの理由で上手く連携できないとそれっきり」という感じだ.だからブラウザでページを再読み込みすると,既に起動している拡張機能と再び連携を試みることになってうまくいく,という様子に見える.まあ,拡張機能の設計が少しまずいのかね. ↩︎