Matplotlibのデフォルト配色
Pythonのグラフ描画パッケージであるmatplotlibは便利なツールである。 僕は最近、論文やプレゼン用のグラフ作成はもっぱらmatplotlibで行っている。
ただ、matplotlibのデフォルトの配色はちょっと気に食わない。
import matplotlib.pyplot as plt for i in range(15): plt.plot( [0,3] , [-i,-i] , lw=4 , label=str(i) ) plt.xlim(0,10) plt.legend( ncol=2 , fontsize=16 ) plt.show()
でプロットした結果は以下のような具合である。
これについて思うこととして、
- 色が全体的に優しすぎる。学術研究のグラフで最も重要なのは結果を明瞭に伝えることなので、もっと強い色を使いたい。
- 好みの問題だが、赤と青は最初に使いたい。もっとも、原色である必要はない。微妙な中間色の方が今風でかっこいい気がする。
- 色がやや薄い。白色背景のグラフで視認性を高めるため、もっと濃い色を使いたい。
といった不満が挙げられる。
配色を変更する一つの手段として、plt.plot()などによる線や点の描画の都度、 引数colorで色を指定することも可能だが、 デフォルトの配色を上の不満を解消したものにしておけばいちいち指定する手間が省ける。 そこで僕は、以下のような手順で変更したデフォルト配色を使っている。
まず、matplotlibのデフォルト配色は、matplotlibの各種デフォルト値を記録したplt.rcParamsの中の
plt.rcParams['axes.prop_cycle']
に記録されており、デフォルト値は、
>>> import matplotlib.pyplot as plt
>>> plt.rcParams['axes.prop_cycle']
cycler('color', ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf'])
となっている。
この変数を変更すればデフォルト配色が変更できる。 例として、
import matplotlib.pyplot as plt from cycler import cycler plt.rcParams['axes.prop_cycle'] = cycler('color' , ['red','green','blue']) for i in range(15): plt.plot( [0,3] , [-i,-i] , lw=4 , label=str(i) ) plt.xlim(0,10) plt.legend( ncol=2 , fontsize=16 ) plt.show()
とすれば以下のように赤、緑、青が順番に使われるようになる。
これで好きな色が指定できるのだが、このようにして色を何種類も指定するのは骨が折れるし、 色の素人が配色を決めてもろくなことにならないので、何らかの既存の配色を使いたい。 そこで、そのような配色としてmatplotlibの中に入っているカラーマップなるものを拝借することにする。 カラーマップは matplotlib.org あたりに一覧があるので眺めてみると、 僕の場合はSet1というやつが最も好みに近い。 ちなみに、デフォルト配色はこの中のtab10というやつである。
カラーマップやcyclerの使い方はなかなかに難しいのだが、
plt.rcParams['axes.prop_cycle'] = cycler( color=[plt.get_cmap('Set1')(i) for i in range(9)])
とすればSet1の配色全９色をmatplotlibのデフォルトに指定できる。
import matplotlib.pyplot as plt from cycler import cycler plt.rcParams['axes.prop_cycle'] = cycler( color=[plt.get_cmap('Set1')(i) for i in range(9)]) for i in range(15): plt.plot( [0,3] , [-i,-i] , lw=4 , label=str(i) ) plt.xlim(0,10) plt.legend( ncol=2 , fontsize=16 ) plt.show()
の結果が、以下のようなものである。
これを見ると、6色目以降の黄色やピンクは視認性が悪くていまいちだが、 最初の5色の配色は僕のイメージにかなり近い。 しかし残念ながら、これらの5色も微妙な中間色でやや薄いので、この配色で論文用のグラフを作ってみると全体的に視認性が低いことに気がつく。 従って、もう少し濃くてはっきりした色が使えればそれに越したことはない。
その一方でSet1の配色自体は悪くないので、ここではSet1をベースにして、 この配色を全体的に濃くすることにする。 そんなことができるのかというと、以下で説明するようにいとも簡単にできてしまうのがPythonの素晴らしいところである。
まず、「色を濃くする」という操作を行うには、RGBではなくHSV (色相、彩度、明度)で色を取得して、 明度(V)の値を調整すればよい。 一方、カラーマップの値はRGBで与えられているので、やるべきこととしては、
- カラーマップSet1の各要素のRGB値をHSVに変換して、
- 明度（と彩度）を適宜調整したうえで
- その色を再びRGBに戻す、
- という操作を行ったcyclerなリストを作成し、
- このリストをplt.rcParams['axes.prop_cycle']に代入する、
ということになる。
こう書くと一見複雑だが、この中のRGBとHSVの交互の変換という操作が
colorsysというパッケージの関数colorsys.rgb_to_hsv()、colorsys.hsv_to_rgb()を使うことで
それぞれ一行に実現できるあたりがさすがPythonである。
以上の操作を行っているのが以下のコードである。
import colorsys clist = [] for i in range(8): c = colorsys.rgb_to_hsv( *plt.get_cmap('Set1')(i)[:3] ) c = colorsys.hsv_to_rgb( c[0] , 1.-(1.-c[1])*.2 , c[2]*.85 ) clist.append( c ) plt.rcParams['axes.prop_cycle'] = cycler( color=clist )
matplotlibによるプロット操作の前にこのコードを挿入しておけばデフォルト配色が変更される。 こうして作った配色によるプロットが以下である。
なお、上の調整ではHSVの第2引数である彩度を上げて(1に近づけて)、 第3引数の明度を0.85倍に小さくしてあるのだが、このあたりはいろいろ試して決めたものである。
この設定で、僕としては5色目までは全く不満がなくなった。 6色目以降はやや不満だが、一つのグラフで６種類以上の色を使うことはまずないし、 そんな場合にはそれこそ配色を真面目に検討して個別にcolorを指定するべきなので、 デフォルト配色としてはこれで十分である。
なお、カラーサイクルを最初の5色のみで回すという選択肢も考えられ、 その場合は上のコードの3行目の'range(8)'を'range(5)'にすれば良い。
この配色を使ったグラフの作例は、
のFigs. 6, 8, 10-13などである。
僕は、matplotlibの各種設定をまとめて1つのpyファイルに書いておいて グラフ描画の際に読み込むことにしており、上の配色設定もこのファイルに書いてある。
僕の場合、これでmatplotlibの配色に関する不満はほぼ完全に解消された。
Jackknife on Python
はじめに
Pythonの便利さに震撼させられたことの一つが、 ジャックナイフ法の汎用的なコードが驚くほど簡素に実装できたことだった。
ジャックナイフ法とは誤差解析法のひとつで、 特にマルコフ鎖モンテカルロ計算における誤差評価で重宝する。 モンテカルロ計算では必須のツールである。
ところが、ジャックナイフのコードを書くのは初心者にとっては意外と難題で、
周囲の学生などを見ているとここでつまづく人も少なくない。
そこでこの記事では、僕がPythonで書いたジャックナイフ解析のルーチンJackKnife()を紹介したい。
ジャックナイフ法
コードの説明の前に、ジャックナイフ法の利点をまとめておこう。 マルコフ鎖モンテカルロ計算の誤差解析でこの方法を使うご利益は、主には以下の二点である。
1) 自己相関の考慮
マルコフ鎖モンテカルロの各測定で得られる測定値は自己相関と呼ばれる相関を持つ。 自己相関はモンテカルロ更新に伴って減衰するため、 各測定間の更新回数を十分大きく取ることで抑制できるのだが、 現実の計算ではこれが難しいことが多々あり、 その場合は誤差評価の段階で自己相関を考慮する必要がある。
ジャックナイフ法では、ビンサイズの調整によって自己相関の効果を誤差評価に取り込める。 より具体的には、ビンサイズを自己相関長より十分長く取るだけでよい。 また、誤差推定値のビンサイズ依存性を調べることで、自己相関長の推定もできる。
2) 誤差伝搬の評価
モンテカルロ計算では、数値計算によって直接測定される量だけではなく、 それらを引数に持つ関数の誤差評価が必要になることがよくある。 つまり、直接の測定量を x1, x2, ... , xn としたとき、 関数 f( x1 , x2 , ... , xn ) の誤差評価が必要になるケースである。
この場合の誤差評価を誤差伝搬法で行おうとすると、 関数 f( x1, x2, ... , xn ) の微係数が必要となることに加え、 各測定値 x1, x2, ... , xn が相関している場合には それも考慮する必要があり、非常にややこしい。 ところが、ジャックナイフ法だと各ジャックナイフサンプルで評価したfの値の分散を使うことで、 これらの効果を全て取り込んだ誤差評価が適切かつ機械的に行える。
使い方
以上のご利益の恩恵にあずかるため、JackKnife()では、
という三種類の操作が簡易なインターフェースで実現できるようになっている。
JackKnife()でこれらの処理を行うには、
まず、モンテカルロ計算の結果として得られる測定値をnumpyの2次元配列
data[i,j]
として準備する。 ただしここで、最初の引数iが異なる測定を表すラベル、 二つ目の引数jが各測定で得られた測定値x1, ... , xnを並べたラベルである。
このようなデータが準備できたら、あとは
average , error = JackKnife( data )
とするだけでx1, ... , xnの期待値とジャックナイフ誤差が計算され、 1次元配列としてaverageとerrorに記録される。
デフォルトではビンサイズ1のジャックナイフ解析が行われるが、ビンサイズを例えば10にしたければ
average , error = JackKnife( data , binsize=10 )
のようにデフォルト引数binsizeを指定すればよい。 ただし、全測定数がbinsizeで割り切れない場合には最後の余った測定値が無視される仕様になっている。
さらに、直接測定量の関数として与えられる量f(x1,...xn)の誤差をジャックナイフ解析したい場合は、 まずそのような関数
def func( arg ): ...
を準備する。 ただしfuncの引数argは直接測定量x1, ... , xnをあらわす1次元配列であり、 funcの戻り値は1次元配列でもよい。 そのうえで、
average , error = JackKnife( data , func , binsize=10 )
とすれば、前と同じようにaverageとerrorに関数funcの期待値とジャックナイフ誤差が記録される。 複数の関数の誤差解析を同時に行いたい場合にはfuncの戻り値を1次元配列にすればよく、 その場合には上の一行で各戻り値の期待値と誤差が一挙に計算され、averageとerrorに1次元配列として記録される。
参考までに、binsizeは第3引数でデフォルト値は1なので、 指定しなかった場合はこの場合もビンサイズ1となるのと、binsizeは第3引数として指定してもよい。
例として、直接測定量x1とx2の比x1/x2のジャックナイフ誤差を解析したければ
average , error = JackKnife( data , lambda arg : arg[0]/arg[1] , binsize=10 )
でOK。
よりマニアックな例だが、格子場の理論における時間相関関数の「有効質量解析」を行いたくて、 data[i,j]の引数jが相関関数の虚時間に対応している場合、
func = lambda dat : np.log( np.roll( dat,1) / dat )[1:]
を第2引数に指定すれば一発である。 （分野外の人には全く分からない例でごめんなさい。）
実装
以上の機能を提供するJackKnife()の実装は以下のとおりである。
import numpy as np def JackKnife( input_data , func = lambda arg:arg , bin_size = 1 ): Ntot , Nvar = input_data.shape N = bin_size Nbin = int( Ntot / bin_size ) #number of bins print( "total:" , Ntot , " / vars:" , Nvar ) print( "bin size:" , N , " / num. of bins:" , Nbin ) if( Ntot != Nbin * N ): print( "Ntot != Nbin * bin_size. Some data will not be used!" ) Ntot = Nbin * N sum_all = np.sum( input_data[:Ntot] , axis=0 ) bin_ans = [ func( ( sum_all-np.sum( input_data[c*N:(c+1)*N] , axis=0 ) )/( Ntot-N ) ) for c in range(Nbin) ] ave = np.mean( bin_ans , axis=0 ) err = np.sqrt( np.mean( np.square(bin_ans-ave) , axis=0 ) * (Nbin-1.) ) return ave , err
このコードをプログラムに貼り付ければJackKnife()が使える。
上の実装の中で、最初の約10行は変数の整理と画面出力をしているだけで、 実際にジャックナイフ誤差評価を行っているのは、以下の４行のみである。
sum_all = np.sum( input_data[:Ntot] , axis=0 ) bin_ans = [ func( ( sum_all-np.sum( input_data[c*N:(c+1)*N] , axis=0 ) )/( Ntot-N ) ) for c in range(Nbin) ] ave = np.mean( bin_ans , axis=0 ) err = np.sqrt( np.mean( np.square(bin_ans-ave) , axis=0 ) * (Nbin-1.) )
2行目で、各ジャックナイフサンプルで計算された測定値がbin_ansに記録される。 bin_ansの計算に先立ち、1行目で一旦全体の和を計算したうえで、 2行目ではそこから引き算することで各ジャックナイフサンプル上での 直接測定値の期待値を得ているが、これは効率化のためである。 3,4行目で、2行目で得たbin_ansから期待値と誤差をそれぞれ求めている。
あの面倒なジャックナイフがこれだけで書けるというのも、いかにもPythonらしい。
np.arange() in c++
久しぶりにPythonからc++に戻ったら、 np.arange()やnp.linspace()、そしてnp.append()が使えないことに 強烈なもどかしさを覚えた。 つまり、vector<double>（あるいはvalarrayかただの配列）の等差数列を 始値と差分あたりを指定して生成したり、それらを一発で結合する関数がほしい。
例として、関数double func( double )の値を等間隔で表示したいとき、 もしc++にPythonのnp.arange()に対応する関数ARange()があれば以下のように書ける。
for( auto x : ARange( 10. , 500. , 10. ) ) cout << x << " " << func( x ) << " " << endl;
これだけなら普通のfor文でもあまり変わらないが、 関数を非等間隔なグリッド上で出力しようとすると通常のfor文だけでは対応が難しくなってくる。 そのような需要が、実は数値計算ではかなりある。 特に関数の処理に時間がかかるときは、 関数の評価を上のように等間隔で行うのは時間の無駄で、 関数形に合わせて非等間隔に調整したグリッド上でのみ関数の評価を行うことで計算コストが大幅に削減できる。 Pythonではこのようなグリッドをnp.arange()とnp.append()の組み合わせで簡易的に作っていたので、 c++で同じことをしようとして、はたと手が止まってしまったのである。
それで、c++におけるnp.arange()やnp.append()の代用を探してみたのだが、 標準c++ライブラリの中にはそのような役割を果たす関数はどうもなさそうである。
さらに調べていくと、c++上でnumpyライクな書式を実現するライブラリ"NumCpp"
https://dpilger26.github.io/NumCpp/doxygen/html/index.html
なるものが見つかったりして興味津々なのだが、
今回は等差数列の生成と結合を行いたいだけなので、このライブラリは大げさすぎるように思う。
こういうとき、どうするのがスマートなのだろうか。
いろいろ考えた末に、自分で実装することにして、 ひとまず以下のようなものを書いた。
vector<double> ARange( double bgn , double end , double dif = 1. ) { vector<double> tmp(1,bgn); for( double val=bgn+dif ; val<end ; val+=dif ) tmp.push_back(val); return tmp; } vector<double> LinSpace( double bgn , double end , int div ) { vector<double> tmp(div+1); double val = bgn; double dif = ( end-bgn ) / (double)div; for( auto& t : tmp ){ t = val; val += dif; } return tmp; } vector<double> Append( const vector<double>& v1 , const vector<double>& v2 , const vector<double>& v3 = vector<double>() , const vector<double>& v4 = vector<double>() , const vector<double>& v5 = vector<double>() ) { vector<double> tmp = v1; tmp.insert( tmp.end() , v2.begin() , v2.end() ); tmp.insert( tmp.end() , v3.begin() , v3.end() ); tmp.insert( tmp.end() , v4.begin() , v4.end() ); tmp.insert( tmp.end() , v5.begin() , v5.end() ); return tmp; }
この実装だと、ARangeで差分difを負にできないとか、 Append()の仕様がnp.append()と違うし同時に5つまでしか結合できないとか、 もっと効率の良い実装があるかもしれないとかの問題があるのだが、 これで冒頭のコードは実現できるし、差し当たっての目的は達成できるので目をつぶろう。
早速、執筆中の論文の数値計算用に書いた関数RateM()の値を ファイル出力するのに使ってみる。
ofstream file( filename ); auto grid = Append( ARange( 1. , 10. , 1.5 ) , ARange( 10. , 40. , 5. ) , ARange( 40. , 100. , 10. ) , ARange( 100. , 200. , 20. ) , ARange( 200. , 500. , 40. ) ); for( auto M : grid ) file << M << " " << RateM( M , C ) << " " << endl;
うん、やりたかったことができて大変満足。
さて、こんな感じの関数たちはヘッダファイルにまとめて汎用的に利用するに値するだろうか。 またその場合、関数の命名規約やファイル名はどうするべきだろうか。 さらに、やるならやるで、もっと柔軟な操作も含めて徹底的にやりたくなるのだが、それなら上記のNumCppを使えという話になるので、そのあたりのさじ加減はどうするのだろうか。
楽しい妄想は続く。
はじめに
数値計算はc++一辺倒でやってきた僕が、時代の波に呑まれてPythonを使い出したのは2018年頃だっただろうか。ひとたびこの言語を使い始めると、手軽さにはまって一気に研究現場での使用率が上がり、特に2021年に入ってからはc++が必要になるような重い計算に出くわさなかったこともあってもっぱらPythonのコードばかり書いていた。 一時期は、もうc++には戻れないかとも思ったのだが、12月に入った頃に必要に迫られて久々にc++を書いたら思ったよりスムーズに指が動いて安心した。
両言語を併用しながらc++とPythonを比較すると、どちらの言語でもイケてるコードが書けて嬉しい瞬間があるのだが、Pythonだとそのようなとき、何やらチートしてゴールにたどり着いたような感覚が残るのに対し、c++だと自分の力で高尚な何かを成し遂げたような達成感があり（大概全く高尚ではないのだけれど）、僕の場合にはやっぱりこっちの方が心地よい。
そんなことを考えているうちに、こうした研究現場のよしなしごとをどこかに書き残したい衝動が抑えられなくなり、久しぶりにネット上での情報発信を始めてみることにした。
■
足指骨折による約一ヶ月間の休止期間を経てジョギングを再開したところ、再開二週間で体重が1.5キロ増えた。人体の不思議。
