カテゴリー: 科学

というなかなか自己矛盾なタイトルでの投稿。　pythonはスクリプトの書く時間を短縮するのが利点であって、pythonに計算をやらせるのはそんなに早くないことが多いので。
だからと言ってふと思い立ったときにCを書くかと言うとちょっと手が渋るのも事実。　この投稿ではそんな大量計算をpythonで頑張ってやってみた時の経験について書きたい。

私的結論から書いておくと道は二つ。　一つはやはり面倒臭がらずにC系統で書く。　pythonを使うことにこだわるならpython的な高速化は求めずに、愚直にスクリプトを書いてからcythonに書き直したりnumba@jitで高速化を目指す。　このどちらかの方が結果として早くなるんじゃないかなと思った。

しかしこの投稿で書くのはそう言う方法ではなく、なるべくpython的に頑張ろうとした時の試行錯誤をいくつか書いていきたい。　もちろんpythonianな人たちはもっとエクセレントな方法を知っているだろう。　気づいたことがあったら教えていただけるととても助かる。

pythonでどのような網羅的計算をしたかったか

ちょっとここの説明が長いので、この章を飛ばしてもらってもなんとなくわかると思う。
４つのパラメータa,b,c,dがある式。　この式は３つの整数変数h,k,lの総組み合わせに対してそれぞれ解を吐き出すような式。　例えば

s = h*a+k*b+c*l+d

といった感じで。　実際はもう少し複雑な計算をさせているが。　いずれにしろ１０個ずつ整数変数があるとしたら一つのa,b,c,dに対して１０００個のs値が計算されることになる。

一方でこの時の実験値が８個。　この実験値と上記で計算された値sとが良く合うパラメータa,b,c,dのセットを探したいというのがやりたかったこと（どのsと実験値を合わせても良い）。
ある程度情報があって最小二乗法などで解けるなら何も問題がない。　そうでない場合はちょっと考えないといけない。

いくつかある方法の中で一番愚直な方法は、とりうるパラメータa,b,c,dを細かい刻みで動かして全ての値のを網羅的に計算する。　その中で一番実験値と合う値を引き算で算出。　その値を叩き出したパラメータa,b,c,dを求めるという方法だ。

プログラムの構成

a,b,c,dのパラメータ変数を適当な間隔で定義する。　各a,b,c,dから式の値を計算して、実験値との差分を計算する関数を回す。　差分が例えば0.001以下のパラメータセットを出力する。　そんな単純プログラム。　しかしa,b,c,dを増やせば増やすほどガンガンと計算量が増えるプログラムだった。

for文を避けてnumpyで頑張る

上記のようなプログラムを書きたい時はCだったらfor文をガスガス回すと言うのは正しい戦略だと思う。　わかりやすいし早いし。
ただpythonはfor文が割と遅いという特徴がある。　for文を如何に避けるかはpythonで重たい計算をするときの一つの方針にはなる。
と言うわけでまず４つのパラメータに対して式を計算するのにfor文を回すのではなく、numpyのarrayに収納してみることをトライ。　どの程度速度が変わるかをまずは単純計算で比較。　４つのパラメータを１０個ずつ取ったとき。　まずは単なる網羅的足し算での比較を。

for文：

import time
import numpy as np

def main():
 start = time.time()
 for a in range(10):
 for b in range(10):
 for c in range(10):
 for d in range(10):
 a+b+c+d
 elapsed_time=time.time()-start
 print ('elapsed_time',elapsed_time)
 return
結果
elapsed_time 0.00111 [s]

array文：

import time
import numpy as np

def main():
 start = time.time()
 a = np.linspace(0,9,10)
 b = np.linspace(0,9,10)
 c = np.linspace(0,9,10)
 d = np.linspace(0,9,10)

 ar1 = np.array(np.meshgrid(a,b,c,d)).T.reshape(-1,4)
 np.sum(ar1,1)

 elapsed_time=time.time()-start
 print ('elapsed_time',elapsed_time)
 return
結果
elapsed_time 0.00107 [s]

あれ？　ちょっとだけ早いけど１００００回の計算だと大して変わらない。　しかもこのnumpyのarrayの組み方は弱点がある。　網目を大きくするとarrayがとっても大きくなる。　その分だけメモリを食うという弱点。
というわけで網目を３０x３０x３０x３０にした場合。

for: elapsed_time 0.0681 [s]
array: elapsed_time 0.1256 [s]

網目を５０x５０x５０x５０にした場合（６２５００００回計算）。

for: elapsed_time 0.51124 [s]
array: elapsed_time 1.07741 [s]

というわけでこれだけだと実はforの方が良かった。

さてこれをちょっとだけ複雑な計算にしてあげるとどうなるか。　各３０パラメータの例で、変更は二乗の和をとるだけ。　変更と結果は下記の通り。
for:
a**2+b**2+c**2+d**2
elapsed_time 0.90076 [s]

array:
ar1 = ar1*ar1
ar2 = np.sum(ar1,1)
elapsed_time 0.13798 [s]

forは二乗を取るだけでかなり遅くなる。　一方でarrayは大して変わらない。　というわけでここでだいぶarrayの方が早くなった。

とはいえこういう単純計算だったらもっと早く書く方法があるのかな。　pythonianな方に教えて頂きたいところ。

ループで使用する関数内での変数宣言を避ける

これは回数の少ないループだとむしろ中に書いておいた方がわかりやすいかもしれない。　しかし１００万回計算するループの中に変数宣言が入っていたりするとpythonはとても遅くなる。
普段はあまり使いたくないgrobal変数だけど、この場合は使用を検討してみると早くなるかもしれない。　grobal変数は、関数の外でgrobalに続けて宣言すればよい。

grobal h
h = np.linspace(0,5,6,dtype=int)

など。

argminよりargpartitionを使う

本題に戻って私がしたかったことは、実験値と計算値を比べること。　これは単純に差分を取ってあげて、そこから最小値を抜き取ればよい。　ここでパッと考えるとargminを使いたくなるところ。　しかしargpartitionという関数があり、こちらの方が早い。

dif2 = np.absolute(slist - dobs)
min =np.argpartition(dif2,0,axis=0)

こんな感じで差分を取ったarrayの最小のインデックスが抜き出せる。

@numba.jit

時折とても強力に高速化してくれる。　けれどむしろ遅くなる場合もある。　全体にかけて効くときもあるし、関数１個だけにかけた方がよく効くときもある。　よくわからないけれど、効いてくれたら嬉しい機能。　試す価値はあり。

数学関数はmathかnumpyか

物によってはmathのが早いし、numpyの方が早かったこともある。　というかあまり変わらない。　最後の最後にどっちのが早いか見てみてもよいかも。　１００万回も計算すれば違いが出るかもしれない。

anacondaのpythonを使う

私がanacondaのpython3に乗り換えたのはこれが原因。　なぜかちょっとだけ計算が速くなる。　もっと早いのがあるなら教えて欲しい。

計算量を減らす努力

小さな努力が大きな結果に跳ね返る場合もある。　計算に対称性などがあって計算が減らせる場合は減らしておく。　計算量が多い場合だとパラメータが減ったり式が単純化されるだけで結構な違いに跳ね返る。　なんだかんだでマニュアル努力が大切。

終わりに

というわけでいくつか私がpythonの高速化を計った時に使った方法を書いてみた。　私は結局目的に使い物になるくらいのところ（１回の試行が帰宅前に回し始めて朝には終わっている）で諦めてしまったのだけど。
最終盤のプログラムでもCで適当に書いた似たようなプログラムに数十倍の速度差をつけられていた。　本職pythonの人がどのくらい早くできるのかは気になるところではある。

cythonへの変換も試した。　しかしいろいろnumpyで書いていたせいか、変数の宣言などcythonを単純に使っても速度を上げることはできなかった。
というわけである程度計算量が必要なプログラムを書く時は、最初からcythonを意識した書き方をしておく。　もしくはcythonフォーマットで書き始めるというのが正しい道なのかな。

関連記事

pythonのまとめ

D

エクセルにpythonが導入される？

最近エクセルにpythonが導入されるかもしれないという記事をいくつか読んだ。　実際に導入されるかはまだ決定はされていない模様。

python as an Excel scripting language

ことの発端は上記のマイクロソフトのユーザーフォーラム。　pythonの導入がだいぶ前にリクエストされていたそうなのだ。　もともとのポストは２０１５年の１１月４日だから２年前。　それに対してマイクロソフトチームが昨年末の２０１７年１１月１４日に突如コメントをつける。　どういう需要があるのかを調査したいのでアンケートに協力してねっ、てことなのでポジティブなコメントと言えるだろう。　しかしアンケートしてるぐらいだしまだまだ初期も初期段階。

マイクロソフトエクセルチームからのアンケート

せっかくなのでアンケートに答えてきた。

1. どこ住んでるの？
2. お仕事の種類はお金目的？
3. 会社の種類はどんな感じ？
4. 研究開発的な専門職についてるの？
5. どんだけの間python使ってる？
6. pythonは何に使ってるの？
7. エクセルとpythonが統合するとしたらどんなことに一番便利かな？
8. ７と同じ質問だけど、こっちは一番じゃないので何個も答えて良い
9. エクセルにpython導入されたら君の仕事のどんな役にたつの？　特に７、８の質問と関連づけてくれるとわかりやすいな。
10. pythonスクリプトは自分で使ってるの？　他の人に提供してるの？
11. いままでに何かエクセルとpythonを繋げるようなソフト使ってた？
12. 他にどんなプログラミングする？
13. 何か他にも意見ある？
14. もうちょっと協力してくれるならemailアドレス教えて

重要な質問は７−９と１１くらいかな。

エクセルとpython統合に関する私見

さてエクセルもpythonもデータ解析にちょろっとずつ使っている私の意見。
自分用にだけならば今のところこの二つを統合する利点は感じていない。　大量計算を簡単にしたいだけならpythonを使う。　もっと大量に早く計算したいとなるとむしろCとかにいく必要がある。
元データがエクセルデータなら簡単な計算はエクセルでそのまま計算。　そうじゃないならばpandasでpythonにインポートして計算させる。

ただし他の人のためのスクリプトを書くとなるとちょっと話しは別かな。　サイエンス界隈の人たちでもエクセルを使ってる人は多い。　裾野を広げたらもっと多くなるだろう。　何かしらスクリプトを書くとしたらデータがエクセルフォーマットでやってくることは多い。　しかも中身が統一されていない様々なフォーマットだったり。　そんな時にいちいちスクリプトを書き直すのは面倒臭い。　エクセルにネイティブにpythonが組み込まれていれば、そこらへんの融通は効きやすくなるんじゃないかなと期待はしている。

といわけでエクセルとpythonの統合で何かが便利になりそうという期待はあまりしていない。　しかし統合されたならば使ってみたい、と思うくらいの興味はあるというのが現状。　上のアンケートなどでとても便利でエクセレントなアイディアを出す人たちがいると良いのだけど。
またエクセルのような大御所がpythonを導入するならば、pythonコミュニティーがまた広がるんじゃないかなと思う。　pythonってわかりやすいし書きやすい良い言語だと思うのだよね。　遅いけど。

関連記事

1. pythonでのテキストデータの入出力方法いろいろ

2. pythonのまとめ

D

pythonを使って積分計算をする機会があった。　備忘にここにまとめておく。

データは上の図の青ラインのようなガビガビの散布データ。　numpyの３６０行のarrayに入っている。　横軸は３６０度の角度で縦軸は強度。　やりたかったことは１度刻みで強度を持つデータに、その角度ごとに三角関数をかけてからその総和を取るということ。
過去の文献をのぞいて見るとシグマで計算している人とインテグラルで計算している人がいる。　もともとの式の定義はインテグラル。　なのでインテグラルを選択するのが順当だろうけど、せっかくなのでpythonを使ってできる計算を色々と書いてみようと思い立った。
まずできるのはデータポイントごとに、三角関数の計算をしてそれからnumpy.sumをとってのシグマ計算。　それからやっぱり積分計算の方が良いかと思い、同じプロセスからnumpy.trapzを使っての積分計算。　とはいってもポイントごとの値のばらつきが大きいデータなのでやっぱりピークフィットしてから、その関数を使って積分計算をかけた方が良いかなとも思い、scipy.integral.quadを使っての積分計算。
というわけでこれらの方法での結果と実行速度を比べてみることにした。

生データでのシグマ計算

まずは横軸を０度から３５９度まで作っておく。　ついでに三角関数用にradianも準備しておく。　例のごとくimportは省略。　numpyとscipy、lmfit、それからTimeくらいかな今回は。

pix = np.linspace(0,359,360)
AzTheta = np.radians(pix)

コンスタントバックグラウンドの準備。　ガビガビのデータなので外れ値を引いてもしょうがない。　numpy.nanminを使って適当な範囲で適当な最小値を選択。

const = np.ones(90)*np.nanmin(Int[140:220])

バックグラウンドを引いたピークの裾からピークトップまでのデータ。　それから各データポイントにそれぞれ０から９０度までのsinとcosをかけた式。

eqsum = (Inttheta[88:179]-const)
eqsin = (Inttheta[88:179]-const)*(np.sin(AzTheta[0:91]))
eqsin = (Inttheta[88:179]-const)*(np.cos(AzTheta[0:91]))

そうしたらsumを使ってarrayの総和をとる。

sum1 = np.sum(eqsum)
sum2 = np.sum(eqsin)
sum3 = np.sum(eqcos)

printで出力。

sum1 66885.0328331
sum2 61156.5549622
sum3 15573.1626916
Calculation time: 0.0002231597900390625[s]

生データからのインテグラル計算

上と同じものをsumではなくnumpy.trapzで計算。　trapzはtrapezoidal ruleでの積分計算。　関数じゃなくても生データから積分計算ができる。

sum1 = np.trapz(eqsum)
sum2 = np.trapz(eqsin)
sum3 = np.trapz(eqcos)

１度刻みのデータなので少しだけ総和よりも小さくなるのは想定通り。　しかしtrapzの方が速度が早かったのは少し意外。

sum1 63719.308712
sum2 58081.6139909
sum3 15482.3795419
Calculation time: 0.00015211105346679688[s]

フィッティングからのインテグラル計算

こちらはscipyのintegrate.quadを用いて。 　フィッティングは以前書いたので省略。　私はいろいろ便利なのでlmfitを使っている。　投稿下にカーブフィッティングへの関連記事リンクを貼っておいた。　ここではconstantバックグラウンドとareaGaussianでフィットした（上の図の赤点線）。
それから積分用の関数を書く。　シンプルなX依存の関数ならそのままintegrate.quadの中に式を書けば良い。
他のパラメータが入る式を使いたい時は先に式を定義しておく。　それからintegrate.quadのargsでパラメータの値を指定することができる。　今回のように最初に関数のフィットをして、その結果の値を積分に突っ込みたいときにはこの方法が便利だと思う。

def eq1(x, amp, cen, hf):
return amp*(2/hf)*np.sqrt((np.log(2)/np.pi))*np.exp((-4*np.log(2))*((cen-x)/hf)**2)

def eq2(x, amp, cen, hf):
return (amp*(2/hf)*np.sqrt((np.log(2)/np.pi))*np.exp((-4*np.log(2))*((cen-x)/hf)**2))*(np.sin(np.radians(x)))

def eq3(x, amp, cen, hf):
return (amp*(2/hf)*np.sqrt((np.log(2)/np.pi))*np.exp((-4*np.log(2))*((cen-x)/hf)**2))*(np.cos(np.radians(x)))

式のパラメータにカーブフィッティングのベストの値をresult.best_values.getを使って代入。　それから０から９０まで積分計算を行う。

s1 = sp.integrate.quad(eq1, 0, 90, args = (result.best_values.get('g1_amp'),90,result.best_values.get('g1_hf')))
s2 = sp.integrate.quad(eq2, 0, 90, args = (result.best_values.get('g1_amp'),90,result.best_values.get('g1_hf')))
s3 = sp.integrate.quad(eq3, 0, 90, args = (result.best_values.get('g1_amp'),90,result.best_values.get('g1_hf')))

生データはピークの裾が持ち上がっているので、フィットしてからの積分値はだいぶ値が異なっている。　それをどう解釈するかはデータ次第。
フィットをもうちょっと真面目にするとしたらVoigtなどを使ってフィットしても良いし。　このデータの場合だと生データのバックグラウンドの取り方を変えても値は近くなるかな。

sum1 48043.715584253885
sum2 47717.13305405471
sum3 4456.932291343271

時間はフィッティングも合わせるとだいぶ遅い。　積分だけでも生データより１０倍遅くなっている。　とはいえよっぽど大量のデータ処理じゃなければ問題はない速度。

total time: 0.011394023895263672[s]
integration time: 0.002457857131958008[s]

というわけで結構値が変わるので、どう計算するかを結構慎重に考える必要がある。　大体の場合先行研究をいくつか見れば答えが書いてあるのだけれど。

関連記事

1. pythonのlmfitを使ったカーブフィッティング

2. pythonのまとめ

D

tkinter、wxpython、pyqtあたりが有名どころのようである。　私はシンプルなことがしたいだけで、見た目も割とどうでも良いということでtkinterを使うことにした。　見た目にこだわったり複雑な処理をしたい場合には、後者二つの方が良いらしい。
使う機会があればこれらについても後ほど書くかもしれないがそんな複雑なスクリプトを書く機会には恵まれないんじゃないかなとも思う。　今の所はtkinterでハッピー。

tkinterをimportするときのエラー

最初にtkinterを使うにあたって一つ最初に躓いたことある。　たぶんインストール設定次第だと思うのだけど、matplotlibと一緒にインポートするとクラッシュする。　matplotlibでTkAggを使うように指定することによって回避できた。

import matplotib
matplotlib.use('TkAgg')
from matplitlib import pyplot as plt
from tkinter import *

といった感じで。

tkinterを使ってユーザー入力補助のスクリプト

ここでは適当なスペクトルデータを表示するプログラムを例に。

tkinterでファイルの読み込み

さてtkinterを使ってまずしたかったのは入力データファイルの読み込み。　自分で使うのにはlsなりdirなりで確認してargvsで読み込めば良い。　もしくはインプットファイル名を決めておいて、勝手に読み込ませるとか。
しかし今回プログラムをphdの学生に渡すのにちょっと手間取った。　特に彼女はあまりpythonというかスクリプトを使わないので。　またwindowsユーザーなのでコマンドライン操作にも慣れがない。
そうなるとプログラムを走らせた時にウインドウを開き、ファイルをクリックして開けるようにした方が間違いがないかなと思ったのだ。
前置きが長くなったが、まずはtkinterの必要なものをインポート。　他のインポートは省略。　brewer2mpl以外はありふれたものしかインポートしてないはず。

from tkinter import Tk, ttk
from tkinter import filedialog as tkfd

それからダミーのルートファイルを開いて隠しておく。　これをしておかないと意味のない小窓が開くので煩わしい。

tk = Tk()
tk.withdraw()

続いてaskopenfilenamesを使ってファイルのパスを取得しfilenamesに格納。　ファイルタイプは適当に。　初期ディレクトリは私はファイルと同じ場所にプログラムを持ってくことが多いので./でその場を指定。　もちろんデスクトップなど好きなところで良い。

filenames = tkfd.askopenfilenames(filetypes= [('txt','*.txt')], initialdir='./')

この操作でこの場合だとテキストファイルを複数開くことができる。　このfilenamesはタプルになるので、必要がある場合はリストに変換する。　取得するのはファイルパス。　なのでファイル名などが欲しい場合はosを使ってコンバージョンする。　他の方法もあるかも。
さてファイルパスが取得できたら一度メインウインドウを閉じておく。　これをしないで他の関数で再びメインウインドウを開くと止まってしまうことがある。

tk.destroy()

filedialogは色々あるのでちょうど良いものを使うと良い。

tkinterでユーザー入力の数値の取得

さてファイルを取得したら、そのファイルをどうしたいかなどをユーザーに決めてもらいたい場合がある。　ここではプロットする範囲を取得するという例で。

num1,num2 = get_range('plot')

get_rangeは以下の自己定義関数。　インデントは省略。

def get_range(word):
#インプット窓の動作を設定するクラス。
class Getr:
#クラス変数の初期化。
input1 =0
input2 =0
#コンストラクタでルートウインドウの修飾をする。
def __init__(self,root):
#インプット窓のラベルの設置。　
#場所の指定はgridが便利だと思うのだけど他はどうなんだろう。
label0 = ttk.Label(root,text = 'input')
label0.grid(row=0,column=1)
#インプット窓のラベルと入力欄の設置。
word1 = 'low ' + word + ' limit'
label = ttk.Label(root,text = word1)
label.grid(row=1,column=0)
self.ent = ttk.Entry(root)
self.ent.grid(row=1,column=1)
self.ent.focus_set()
#同じようにもう一つラベルと入力欄の設置。
word2 = 'high ' + word + ' limit'
label2 = ttk.Label(root,text = word2)
label2.grid(row=2,column=0)
self.ent2 = ttk.Entry(root)
self.ent2.grid(row=2,column=1)
self.ent2.focus_set()
#ボタンの設置と、ボタンを押した時の動作の設定。
#lambdaの理由は忘れたけど、これないと動かなかったはず。
but = ttk.Button(root, text='OK',width=10,command = lambda: self.get_quit())
but.grid(row=3,column=0)

#ボタン動作用のメソッドをいくつか。
def gettext(self):
Getr.input1 = self.ent.get()
Getr.input2 = self.ent2.get()
def quit(self):
root.destroy()
def get_quit(self):
self.gettext()
self.quit()
#クラスはここでおしまい。
#インプットルートウインドウ。　ここら辺はテンプレ。
root = Tk()
#メインウインドウでクラスGetrを。
range=Getr(root)
root.mainloop()
#インプットの取得。
return Getr.input1, Getr.input2

さてこうしてユーザー入力で得られる数値とプログラムで使いたいデータの列番号は異なっていたりするので、列番号に直しておく。　ここでargpartitionを使っている理由は早いからだけど、この場合はあまり関係ない。

def conv_line(plot,num):
num = float(num)
dif = np.absolute(plot-num)
mini = np.argpartition(dif,0,axis=0)
return mini[0][0]

こうして得られたライン１からライン２までのプロットを打つ。　ちなみに想定してるのはデータはx軸とy軸の２列データ。

def plot_raw(filenames,num1,num2):
num = len(filenames)
print (num," files read")
if num >2 and num <11:
bmap = brewer2mpl.get_map('spectral','diverging', num)
bmap2 = brewer2mpl.get_map('OrRd','sequential', num)
colors =bmap.mpl_colors
colors2 =bmap2.mpl_colors
for i in range(num):
plot = np.loadtxt(filenames[i])
lin1 = conv_line(plot,num1)
lin2 = conv_line(plot,num2)
plt.plot(plot[lin1:lin2,0],plot[lin1:lin2,1], '-',color=colors[i],mec=[1,0,0],ms=8,lw=0.5, alpha=1)
return

といった感じのスクリプトでした。　あえてpythonでプロットを打つ必要はないのだけど、打てると便利な場合もある。　私はバックグラウンド補正とか強度のキャリブレーションを読み込んだ複数ファイルにまとめてかけて、それからプロットを出力するのに使っている。　というわけでtkinterの一例でした。

関連記事

pythonのまとめ

D

最近phdの学生のwindowsコンピュータへのpythonのセットアップを手伝うことがあった。　ただWindowsにインストールするだけなら簡単なはずだが、私が最近Windowsを使っていなかったこともありそう簡単でもなかった。
今回の投稿ではマックでのpyenvを使ったpython3のセットアップと、Windowsにpythonをインストールした体験を記録しておきたい。

Macでのpyenvを使ったpythonのインストール

さてMacでのセットアップはとても簡単だ。　というかxcodeとか入れてばどれかのバージョンが入っているだろう。　また適当なソフトウエアパッケージを使えば好きなバージョンのpythonもインストールできるだろう。
しかしpythonはpython２、python３をはじめとしてバージョンが大量にあるのでいくつかダウンロードしてしまうとその管理が大変。　いろんなソフトウエアで一緒にダウンロードしてたりすると、今どのpythonを使っているかがわからなくなってしまったりもする。
そういった問題を避けるために私はpyenvだけ使ってpythonを管理している。　pyenvのインストールはhomebrewからできる。　fink、macportじゃないとダメなソフトもあるけど、最近はhomebrewもかなり頑張っている印象。　homebrew のインストールは公式ホームページのコマンドをコピーしてきて、ターミナルに貼り付けて実行するだけで良い。　続いてhomebrewを使ってpyenvのインストール。

brew install pyenv

それからUsers/”user”/下の.bash_profileに以下を記述しておく。

export PYENV_ROOT="${HOME}/.pyenv"
export PATH="${PYENV_ROOT}/bin:$PATH"
eval "$(pyenv init -)"

続いてpyenvを使って必要なpythonのバージョンをインストールしておく。

pyenv install "version"

続いてメインで使うバージョンをグローバルにセットしておく。　私はanacondaのpython3をセットしているけどまあお好みで。

pyenv global "version"

あとはあるフォルダの中だけで別のpythonを使いたかったりしたら、そのフォルダに移動してから

pyenv local "version2"

といった具合で設定する。　ちなみにpipなどで入れるパッケージはバージョンごとに入れなおす必要がある。
その他たまに使うコマンド。
pyenv install –list　インストール可能なpythonバージョン一覧
pyenv versions　　インストール済みのバージョン一覧

Windowsコンピュータでのセットアップ

普通のWindowsでのインストールならそんなに面倒くさくもないと思う。　公式でも良いし何か適当なpythonをダウンロードしてインストールすればよい。　そうしたらpythonのcommand lineかPython GUIを使って、書いたりプログラムを走らせたりすることができる。
ただpythonのGUIから入るとディレクトリの移動とかがちょっと面倒臭い。　osをimportしてディレクトリ移動とかできるけど、ちょっとコマンドが長いし覚えにくい。
そういうわけでpythonをターミナルから使いたかったら、コマンドプロンプトを開いてC:\python32\pythonのように実行すれば良い。　もしくはパスを通しておくか。　コマンドプロンプトならcdとかdirとか割とみたことがあるようなコマンドで操作できるので便利だ。　それからpipでも使って好きなようにパッケージを入れるなどカスタマイズしていけば良い。

問題は最近大学や研究所などソフトウエアのインストールを管理している場合がある。　そうすると追加のソフトウエアのインストールを毎度ITチームにお願いしないといけない。　そうなると上記の方法だと私のプログラムを走らせるセットアップにするのがちょっと面倒くさい。
なので最小限のソフトウエアインストールで切り抜ける方法を探してみた。　私はanacondaのpython3を使っているので、phdの学生にはとりあえずそれをインストールしておくように頼んだ。
最初はpythonをコマンドプロンプトから使おうかと思い、実行ファイルの場所を探していのだ。　しかしWindowsにanaconda-python3をインストールすると、実行ファイルの位置がいまいちわかりにくかった。　そんな探している最中にanacondaのパッケージの中にanaconda promptなるものを発見。
なんとなく使ってみると、このターミナルで大体問題なく書いておいたプログラムを使えることが判明した。　pipもこのターミナルから直接使えたので、追加のパッケージのインストールも簡単にできた。
というわけでWindowsコンピュータではanaconda-python3だけインストールしとけばだいたいなんとかなりそうな気がしている。

関連記事

pythonのまとめ

D

この投稿ではpythonでテキストデータを処理するときの、データの入出力方法いろいろについてまとめたい。　いやまとめというか私がよく使っている方法を書くだけだが。　ついでにそんなにいろいろもなかったが。

データの読み込みの方法

スクリプトを書くときは、何かしらデータを読み込んで処理を加えて結果なりなんなりを出力させるってプロセスを経ることは多いと思う。　そうするには最初にデータを読み込む必要がある。
データといっても千差万別なので、データによって読み込み方法を変える必要がある。　既にraw画像の読み込みの仕方について投稿を一つ書いてあるので、興味があればそちらもご覧いただきたい。
この投稿ではテキストデータを読み込む場合について書きたいと思う。　テキストデータといっても、数値や文字などいろいろな種類がある。　数値データ一つをとっても、整数や小数などデータタイプの違いがある。　さらにそれらの混合タイプのデータもあるので、読み込むためには色々と設定が必要なこともある。

数値データの読み込み方法

さてまずはシンプルな区切り数値データを読み込む場合。　便利なのはnumpyのloadtxt。　いずれにしろpythonで計算などをする場合はnumpyでプロセスしたい場合が多いので、numpyのarrayに直接入れられる方が便利だろう。

pwd = np.loadtxt(argvs[1])

といった感じで読み込む。　データの形式を指定する必要がある場合は、dtypeでint,float,doubleなど。　delimiterも必要ならばtabやspaceなど指定できるが、あまり必要ないことが多い。
問題点は同じ行列数のデータでないとエラーを吐き出す。　ヘッダーやコメントが付いているくらいなら、skiprowsで最初の行を飛ばしたりcommentsを利用して指定文字で始まる行をスキップしたりはできる。
これを書くのに公式マニュアルを確認していたら、convertersっていうオプションが追加されていた。　これを使うともう少し便利なのかもしれないが、まだ試してはいない。

文字と数値データが混ざっている場合

上記の方法は一行に数値と文字データが混ざっているようなテキストには普通は使えないことが多い。　文字データだけを読み込みたい場合は、目的に合わせてreadやreadline、readlinesを使って読み込めば良いだろう。　文字データと数値データが混ざっている場合も基本的には同じ方針で読み込める。
さて数値データを抜き出したい場合だが、例えばスペース区切りのテキストだったら以下のように読み込むことができる。　この投稿ではインポートは省略しているので、適時必要なパッケージはインポートしていただきたい。

try:
f = open('input.txt')
except IndexError:
print ('cannot open input file')
else:
input = f.readline().split()
f.close()

このようにして、スペース区切りのテキストがinput[i]に収納することができる。　こうして得られるinput[i]は文字列なので、数字として使いたい場合は変換する必要がある。　int(input[1])やfloat(input[2])などといった感じで。　Read系統を使っていく利点は、間違いなく読み込みはできること。　あとはどう読み込んだデータを使って行くかを工夫して行けば良い。
さてもうちょっとマシなデータフォーマットが整っている場合だと直接numpyに突っ込みたい時もある。　そんな時に便利なのがnumpyのgenfromtxtだ。

temp = np.genfromtxt(argvs[2],dtype=['U8','int','U8','U8','int','float','float','float','float','float'],skip_header=4)

例えばこんな感じで色々混ざったデータを読み込んで、その中の数値データだけを計算に使うなどといったことができる。　genfromtxtは多くのオプションがあるので、色々と読み込みの設定ができる。

２０１８年３月３１日追記：

genfromtxtについてもう少し書いた記事を投稿した。

リンク：pythonを使って座標変換計算　〜回転行列を使って座標回転

EXCELデータの読み込みの方法

２０１７年１２月１０日追記：
色々試したけれど、EXCEL系のデータは有名パッケージのpandasがどうにも優秀だった。　まずはインポート。

import pandas as pd

それからエクセルシートをread_excelを使ってダイレクトに読み込む。　filenamesはtkinterで取得。　tkinterの使い方に興味がある方は投稿下部の関連記事から。　それから読み込むsheet名の指定。　それからsheetのデータを読み込む。

excel = pd.read_excel(filenames[0], sheetname='Sheet1')
data = excel.values

これで各行列が読み込まれている。　例えば何かしらプロットデータだったら下のように列を指定してあげればプロットできる。　便利。

plt.plot(data[:,0], data[:,1])

データの書き出しの方法

numpyの数値データを書き出す場合、単純にsavetxtというオプションを使ってあげれば良い。

np.savetxt('name.txt', DATA)

といった具合。
genfromtxtなどで読み込んだ文字列と数値が混ざったテキストを出力したい場合もsavetxtで良い。　ただしフォーマットを指定してあげる必要がある。　これは結晶のpdbフォーマットでの出力の例。

np.savetxt('test.pdb',dft_pdb,fmt='%4s %6d %3s %-4s %2d %11.3f %7.3f %7.3f %4.3f %5.2f')

numpyデータでない場合は、readと同じように書き込みモードでファイルを開てwriteやwritelinesを使って出力する。
あと便利なのは同じように書き込みモードでファイルを開いておいて、プリントでファイルに書き込んでいく方法。　ファイルにスクリプトの違う場所から次々と書き込んでいきたいときなどに便利。　この例では文字列だがもちろんデータも書き込める。

res = open('result.txt','w')
print ("test test test",file=res)

以上。　新しい入出力を試す機会があったら追記する。

関連記事

pythonのまとめ

D

今回の投稿ではpythonのmatplotlibを使ったプロットの描画の仕方をまとめておきたいと思う。

プロットを作るためのソフトウェア

さて数値データを人に説明しようかと思うと、やはり何らかのプロット（絵）に変換して見せるのがわかりやすいかと思う。　数値データからプロットを作るソフトウエアはエクセルをはじめとして各種の優れた有料ソフトウエアが存在している。　普通はエクセルか所属機関が購入しているソフトを使うことが多いのだろうけど、見た目にこだわりがあるなら自分の予算次第で買うのもありだろう。　そんなに高いわけでもない。
無料ソフトで作りたいと思ったら、gnuplotなどはかなり高性能な上に見た目も良いプロットを作ることができる。　プログラミングと組み合わせて使っている人が多いソフトでもある。　使い勝手に対する意見は分かれるところかもしれないが。

pythonを使ったプロットの描画（matplotlib）

さて我らがpythonはどうなんだというと、matplotlibというパッケージが用意されている。　pythonを使い慣れていればそれなりに使い勝手は良いし、絵も努力次第ではかなり綺麗な絵を出力することができる。　pythonは使っている人が多いので、色合いやプロット形式などのフォーマットを拾ってくれば簡単にクオリティーの高い絵を出力することができる。　もちろん自分の好みを追求していくことも可能だ。
とはいえpythonでプロット出力のスクリプトを書くのには、それなりに時間もかかる。　一つのプロットを作るのにスクリプトを書くのは、あまり賢いことではないだろう（私は自分に練習だからと言い聞かせて書いてしまうこともあるが）。　一方で似たような形式のデータから多数のプロットを作製したい時などは、スクリプトを書くことで時間の短縮ができることもある。　さらに何かのプログラムを書いた時に、ついでにプロットを一緒に出力してしまうのも良い方法だ。　例えば私は良く２Dイメージを強度vsピクセルのプロットに変換するのだが、テキストデータを作ると同時にプロットも出力してしまうようにスクリプトを書いている。

matoplotlibを使ったスクリプトの例

matplotlibはanacondaのpythonを使っていれば標準インストールされていると思う。　そうでない場合はpipでインストール。

pip install matplotlib

スクリプトの最初にmatplotlibのpyplotをインポートしておく。　基本的にはこれだけで良いと思うけど、必要に応じて他の関数もインポートする。　私はグレイスケールの画像を使うのにmatplotlib.cmなどもインポートしている。

import matplotlib.pyplot as plt

私の場合他のパッケージとの競合でのバグ避けようにおまじない。

matplotlib.use('TkAgg')

さてプロットの基本設定は紛らわしくないように、メイン関数の中にでも書いておくと良いかと思う。上からフォント・サイズ・各軸の線の太さ・x軸のメモリの太さ・x軸のラベルサイズ・y軸のメモリの太さ・y軸のラベルサイズ。　最後のは絵を出力する際のサイズ比。　調べてないけどrcParamsはもっと色々設定できるんじゃないだろうか。

plt.rcParams['font.family']='Arial'
plt.rcParams['font.size']= 20
plt.rcParams['axes.linewidth']=3.5
plt.rcParams['xtick.major.width']=2.5
plt.rcParams['xtick.labelsize']=20
plt.rcParams['ytick.major.width']=2.5
plt.rcParams['ytick.labelsize']=20
plt.rcParams['figure.figsize']=(8,6)

さてデータは適当に読み込む。　numpyだったらこんな感じ？

plot = np.loadtxt(argvs[1])

そうしたらプロットを作る関数なりクラスなりを書く。　この場合yが縦軸のデータでx軸が横軸になる。　–は線グラフになるが、oなどを使うと散布図で出力できる。

plot_data(plot)
def plot_data(plot):
y = plot[50:2000,1]
x = plot[50:2000,0]
plt.plot(x, y, '--', colour = [1,1,1], linewidth =2, label='test')

といった感じ。
細かいプロットの調整は関数の中でしても良いかも。　特に一つのプログラムの中で複数のプロットを出力するのに、それぞれ設定を変えたい時なんかに。　そんな時は一つプロットが終わって絵を吐き出したらplt.clf()でリセットしておく。
話を戻してX軸y軸のラベルの調整。

plt.tick_params(axis ='y', which ='both', left='off', right='off', labelleft='on', pad=10)
plt.tick_params(axis ='x', which ='both', top='off',pad=8)

縦横軸の設定。　フォントは初期設定から変えたい場合はここでも変えられる。　特殊記号も下記のようにtex形式やらコードやらで記入可能。

plt.xlabel(r"$\theta$"' (\u00B0)', fontsize = 30, fontname='Arial')
plt.ylabel("Int", fontsize = 30, fontname ='Arial')

レジェンドを図中に入れるなら。　locで場所の指定ができる。　bestみたいなお任せで適当に入れてもらう機能もあるが、必ずしもよくない。

plt.legend(fontsize =14, loc='upper right')

図の保存。　サイズ指定で絵を出力するとはみ出ることがあるので、bbox_inches=’tight’を指定しとくときっちり収まる。

plt.savefig(name, bbox_inches='tight')

あとはx軸とy軸の範囲指定はよく使う。

plt.xlim(5, 100)
plt.ylim(0,10000)

などなど。　最後に保存した絵に加えて、pythonにもプロット出力させておきたかったら。

plt.show()

基本的な使い方はこんな感じでしょうか。　あとはプロットする色の組み合わせなど外部パッケージを使うと少し楽できる。　私が使ったことがあるのはbrewer2mpl。　まずは色のグラデーションを作ってあげる。

import brewer2mpl
bmap = brewer2mpl.get_map('Blues','sequential', 7)
bmap2 = brewer2mpl.get_map('OrRd','sequential', 7)
colors =bmap.mpl_colors
colors2 =bmap2.mpl_colors

それからカラーをcolor[]で指定して出力してあげる。

plt.plot(arr1[50:1000,0],17000+arr1[50:1000,2], '-',color=colors[6], lw=2, alpha=0.7, label='WATER')

てな具合に。
今思いついたのはこのぐらい。　何か有用なことでも思い出したら追記するが、あまり思いつきそうなことももないかなとも思っている。

関連記事

pythonのまとめ

D