カテゴリー: プログラミング

1変数連続関数の根(値が0になる点)を求める、二分法というアルゴリズムとそれをScipyで実装する方法を紹介します。

アルゴリズムの内容

二分法というのは中間値の定理をベースとした求根アルゴリズムです。アイデアは非常に単純で、連続関数$f$に対して、$f(x_1)$と$f(x_2)$の符号が異なるように、$x_1, x_2$を選びます。この時点で、中間値の定理より区間$(x_1, x_2)$に根があることがわかりますのでさらに細かくみていきます。次は、$x_1, x_2$の中点$x_M = \frac{x_1+x_2}{2}$を取り、$f(x_M)$の符号を調べます。$f(x_M)$と$f(x_1)$の符号が同じであれば、$x_1$を$x_M$で置き換え、逆に$f(x_M)$と$f(x_2)$の符号が同じであれば、$x_2$を$x_M$で置き換えると、区間の幅が半分になった$(x_1, x_2)$が得られますが根はこの中にあることがわかります。これを繰り返すと、根が存在する範囲を狭めていくことができ、$f(x_M)$の絶対値が0になるか、もしくは十分0に近づいたらその値を数値的に求めた根とします。

以上が、一般的な二分法のアルゴリズムの説明です。ただし、後に紹介するSciPyではどうやら区間が十分狭くなったかイレーション回数が上限に達したか等の基準でループを打ち切り、$f(x)$の値を確認していないようです。

二分法のメリットとデメリット

方法が単純でわかりやすい、というのが個人的に感じている１番のメリットです。

また、連続関数であれば使えるため、ニュートン法などのアルゴリズムのように元の関数の微分を必要とせず、微分が難しい関数や微分不可能な関数でも使えます。

また、根が存在しうる区間を狭めながら探索するため、最初の区間の幅と繰り返し回数により、結果の精度を保証できることも大きな利点です。

逆にデメリットとしては、ニュートン法等と比較して収束が遅いこととか、初期値として関数が異符号になる2点を探して与える必要があること、もしその区間に複数の根が存在した場合にどれに収束するか不確定なことなどが挙げられます。

ただし僕の経験上では、ある程度根の目処がついていたり、単調な関数に対して使う場面が多くこれらのデメリットを深刻に感じることは少ないです。

SciPyによる実装

SciPyではscipy.optimizeというモジュールで実装されています。bisectという専用メソッドを使うか、root_scalarという汎用的なメソッドで(method=’bisect’)を指定して使うことになります。

参考:
– scipy.optimize.bisect — SciPy v1.11.4 Manual
– root_scalar(method=’bisect’) — SciPy v1.11.4 Manual

試しに、$\sin$関数の 3と4の間にある根を探させて見せましょう。既知の通りそれは円周率$\pi$になるはずです。

from scipy import optimize
import numpy as np


root1 = optimize.bisect(np.sin, 3, 4)
print(root1)
# 3.1415926535901235

root_result = optimize.root_scalar(np.sin, bracket=[3, 4], method='bisect')
# 結果は複数の情報を含むRootResults形式で戻る。
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 41
     iterations: 39
           root: 3.1415926535901235
"""
print(type(root_result))
# class 'scipy.optimize._zeros_py.RootResults'

# 根の値へのアクセス方法
print(root_result.root)
# 3.1415926535901235

いかにも円周率ぽい結果が得られましたね。root_scalarの方では収束したことを示すフラグや、イテレーション回数なども得られています。

失敗事例1. 区間の両端での関数の値が同符号の場合

二分法は初期設定を誤ってると失敗するのでその場合のSciPyの挙動も見ておきましょう。失敗パターンの一つは、最初に指定した区間の両端で符号が一致していた場合です。もちろん関数の形によってはその区間内に根がある可能性もあるのですが、存在は保証されなくなります。

また$\sin$関数で、その間に根が存在しない区間$(1, 2)$と、実は両端で同符号だけど根が存在する区間$(1, 7)$でやってみましょう。bisectとroot_scalarで全く同じエラーメッセージ出るのでbisectの方だけ載せます。

try:
    optimize.bisect(np.sin, 1, 2)
except Exception as e:
    print(e)
# f(a) and f(b) must have different signs

try:
    optimize.bisect(np.sin, 1, 7)
except Exception as e:
    print(e)
# f(a) and f(b) must have different signs

はい、$f(a)$と$f(b)$は違う符号にせよとのことでした。根が存在しない1個目の例はさておき、2個目の例は根は区間内に2個存在するのですが探さずにエラーになりました。

失敗事例2. 連続関数ではなかった場合

もう一つ失敗するのは、関数が連続関数ではないケースです。

例えば$\tan$の、$\frac{\pi}{2}$近辺の挙動で見てみましょう。数学的に厳密な話をすると、$\frac{\pi}{2}$では$\tan$は定義されないので、$\tan$は数学的には連続関数(定義域内のすべての点で連続)なのですが、数値計算的には不連続と考えた方が都合が良いです。

話が脇に逸れたので実例の話に移ります。実は、bisectメソッドは結果を返して来ちゃうんですよね。そしてそれが全然根ではないということも見ておきましょう。

root = optimize.bisect(np.tan, 1, 2)
# pi/2に近い結果が得られている
print(root)
# 1.5707963267941523

# 元の関数に代入した結果は全く0に近くない
np.tan(root)
# 1343445450736.3804

root_scalarの方だったら、結果のフラグ等もあるのでアラート等あげてくれるのかと期待したのですがそういう機能はなさそうです。

root_result = optimize.root_scalar(np.tan, bracket=[1, 2], method='bisect')
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 41
     iterations: 39
           root: 1.5707963267941523
"""
np.tan(root_result.root)
# 1343445450736.3804

要するに、SciPyに渡す関数が本当に連続関数であるかどうかは利用者が責任を持たないといけないということです。また、結果が本当に根なのかどうかは代入して確認した方が良いでしょう。

まとめ

以上が二分法の説明とSciPyで利用する方法、その注意点でした。

はじめに

今週の記事もPythonスクリプトで引数を受け取って使う話です。前回はsys.argvつかって受け取る方法を紹介していましたが、今回は便利な専用モジュールのargparseを紹介します。

これを使うと、引数を変数に自動的に格納したり、オプション引数やフラグを作成したり、ヘルプ機能を自動的に作ってくれたりします。

よくUnix/Linux コマンドでは -o filename みたいな感じで出力先ファイルを指定できたりしますが、これをsys.argvで実装しようとすると、配列を全部見て-oがあるかどうか確にして、その次の値をfilenameとして取得して、みたいな結構面倒な処理を自分で作る必要があります。-oが複数出てきたらどうするかとか、-oの次にファイル名がなかった場合のハンドリングとか色々考えないといけないのでとても面倒です。こういう手間を削減してくれます。

順に使い方書いていきますが、ドキュメントはこちらです。
参考: argparse — コマンドラインオプション、引数、サブコマンドのパーサー — Python 3.12.1 ドキュメント

基本的な使い方

ざっくりいうと、argparseは次の3手順で使います。

• ArgumentParserオブジェクトの作成
• 必要な引数をパーサーオブジェクトに追加する
• 引数を解析して結果を取得する

一回単純なサンプル作ってやってみましょう。上記ドキュメントの例をそのまま使います。
sample.py というファイル名で次のスクリプトを作成し、実行権限を `$ chmod u+x sample.py` でつけておきます。

#!/usr/bin/env python
import argparse


# パーサーオブジェクトの作成
parser = argparse.ArgumentParser(
    prog="ProgramName",
    description="What the program does",
    epilog="Text at the bottom of help"
)

# 必要な引数の追加
parser.add_argument("filename")  # 位置引数
parser.add_argument("-c", "--count")  # 値を取るオプション　
parser.add_argument("-v", "--verbose", action="store_true")  # on/off フラグ

# 引数の解析
args = parser.parse_args()
print(args.filename, args.count, args.verbose)

だいたいイメージできると思うのですが、./sample.py を実行する時、最初の引数が filename に格納されて、 -c か –count で指定した値が count変数に格納され、 -v を選択したかどうかがTrue/False で verbos に入ります。 色々やってみましょう。

$ ./sample.py test.txt
test.txt None False

$ ./sample.py test.txt -c 4 -v
test.txt 4 True

$ ./sample.py -v --count abc test.txt
test.txt abc True

$ ./sample.py test1.txt test2.txt 
usage: ProgramName [-h] [-c COUNT] [-v] filename
ProgramName: error: unrecognized arguments: test2.txt

$ ./sample.py                    
usage: ProgramName [-h] [-c COUNT] [-v] filename
ProgramName: error: the following arguments are required: filename

$ ./sample.py --help
usage: ProgramName [-h] [-c COUNT] [-v] filename

What the program does

positional arguments:
  filename

options:
  -h, --help            show this help message and exit
  -c COUNT, --count COUNT
  -v, --verbose

Text at the bottom of help

はい、最初の3例が正しくコマンドを打ったケースでしたが、だいたいイメージ通りに引数を受け取れていることが確認できると思います。
4つ目は位置引数を過剰に設定、5つ目は逆に指定しませんでしたが、それぞれちゃんとエラー文を出してくれていますね。
6個目の例は–helpをつけていますが、なんと自動的にヘルプメッセージを作成して表示してくれています。

コマンド名が ProgramName になっていますが、これはパーサーを作成したときのprog 引数をプログラム名として使っているからです。progを省略すると、ファイル名が使われます。

これは大事なことなのですが、プログラム名=ファイル名のことが多いと思うので、基本的に省略した方がいいと思います。(さっきの例は公式ドキュメントをただ真似しただけ。)

descriptionでプログラム中身の説明、epilogでヘルプの最後に表示するメッセージを指定できますが、これらもどちらも省略可能です。ただ、descriptionは何か書いていておいた方がいいと思います。

ここから細かく仕様を見ていきます。

引数の種類

引数の種類としては、コマンドの後に何番目に渡されたかどうかで扱いが決まる位置引数と、-(ハイフン)付きの名前で始まるオプション引数があります。

argparseは接頭辞の”-“を特別な文字として扱って、これによって挙動を変えています。

上の例でもわかりますが、次のように複数の名前を指定することもできますし、1種類だけの名前でも良いです。このとき注意しないといけないのは、参照するときの変数名です。

-c みたいな短い名前だけの時は – をとって c として参照しますが、-c, –count という2種類の名前を指定した場合は、一番最初に登場する長い名前、が採用されます。長い名前というのは文字列の長さの話ではなく、 – ではなく、 — で始まる引数ということです。
つまり、次のように3つの名前をつけたら、長い名前の中で最初に登場した countが採用されるということです。

#!/usr/bin/env python
import argparse


# パーサーオブジェクトの作成
parser = argparse.ArgumentParser()

# 必要な引数の追加
parser.add_argument("-b")
parser.add_argument("-c", "--count", "--cnt")

# 引数の解析
args = parser.parse_args()
# 短い名前 -b しかないのでbでアクセス
print(args.b)

# --count と --cnt が長い名前だが、先に登場したcountの方が優先
print(args.count)

add_argument の引数

add_argument には様々な引数を指定でき、各種の設定を行うことができます。
全部紹介するのも大変なので一部抜粋して紹介しますが、公式ドキュメントの該当欄の一読をお勧めします。
参考: add_argument() メソッド

default ･･･ コマンドラインに対応する引数が存在せず、さらに namespace オブジェクトにも存在しない場合に利用されるデフォルト値。
type ･･･ データ型。 int や float、ユーザー定義の型など色々指定できる。省略すると文字列(str)。
choices – 引数として許される値のシーケンス。
help ･･･ 引数の説明。-h や –help を使用した時に使われる。
nargs ･･･ 受け取れるコマンドライン引数の数。後で説明します。
action ･･･ コマンドラインにこの引数があったときの動作。後で説明します。

だいたいはイメージ通りの挙動をしてくれるのですが、nargsとactionについてはこの後説明します。

受け取れるコマンドライン引数の数について

nargs という値を使って、受け取れる引数の数を指定できます。

nargs の指定は正規表現風になっています。 整数Nを指定すればその個数、?なら1個か0個で、0だったらdefalut値が使われます。*とすると任意の数受け取れます。また、+だとこちらも任意の数受け取れますが、0個だった場合にエラーが起きます。

例えば、任意の数のファイルのデータを入力として、1個のファイルに結果を書き出すようなコマンドがあったとしましょう。(というより、tar コマンドでアーカイブ作る時はそういう指定しますよね。 )

次のような形です。

#!/usr/bin/env python
import argparse


# パーサーオブジェクトの作成
parser = argparse.ArgumentParser()

# 必要な引数の追加
parser.add_argument("-o", "--out_file")  # nargsを省略しているので1
parser.add_argument("-i", "--in_file", nargs="+")

# 引数の解析
args = parser.parse_args()
print(args.in_file)
print(args.out_file)

# 以下実行結果
$ ./sample.py -o out.txt  -i in1.txt in2.txt in3.txt
['in1.txt', 'in2.txt', 'in3.txt']
out.txt

— in_file の方は複数の結果を受け取れるようにしたので、Python上は配列で結果が来るようになりましたね。

actionによる動作の指定について

actionを使って、オプション引数が存在したときの挙動を指定できます。

デフォルトは store でこれは要するに変数を値に格納するという挙動です。さっきまで見てるのがこれですね。

ただし、Linux/Unixコマンドではこのような値を受け取る引数ばかりではありません。皆さんがよく使う ls コマンドの -l とか -a は別に何か引数を受け取ったりせず、その存在の有無だけが重要ですよね。

この記事の冒頭のコードの `parser.add_argument(“-v”, “–verbose”, action=”store_true”) # on/off フラグ`
もまさにそうで、 -v の有無だけが問題になります。これを実現しているのが、action=”store_true”の部分です。

要するに -v が見つかったら verboseにTrueを格納するよ、という挙動になります。
逆に見つからなかったらFalseが格納されます。

これと逆にオプションがあったらFalseでなかったらTrueになるのが、”store_false”です。

このほか、キーワード引数の登場回数を数えて格納する”count”とか、複数回登場したら結果を都度配列に追加していく”append”などもあります。

これらも一通り公式ドキュメントの一読をお勧めします。
参考: action

ヘルプの作成について

最後にヘルプ機能についてです。自動的に、-h と –helpがヘルプ機能として実装されます。

これはもう実際に試していただくのが一番早いのですが、description等で指定されたプログラムの説明や、受け取れるコマンドライン引数の情報などが表示でき大変便利です。

気をつけないといけないのは、 -h と –help を上書きしないようにすることですね。もちろんどうしてもこれらの引数名を別用途で使いたいとか、自作のヘルプメッセージを実装したいとか事情があれば話は別ですが、普通はデフォルトのヘルプを使った方が良いと思います。

argparseをコマンドライン引数以外の文字列のパースに使う

最後にちょっとマニアックな使い方を紹介します。
このargparseですが、何も指定せずに、 parser.parse_args() すると コマンドライン引数をパースしにいきますが、ここで配列を渡すとその配列をパースします。

sample_str = “-i filename -c 5” みたいな文字列があった時にsample_str.split()して配列に分解して、 parser.parse_args(sample_str.split())と渡すとそれをコマンドライン引数と見立ててパースしてくれるのです。

そんな技術いつ使うねん、と思われるかもしれませんが、僕はマジックコマンドを作る時などに使ってます。
参考: Snowflakeに手軽にSQLを打てるJupyterマジックコマンドを作ってみた｜ホンディー | ライフイズテック 　(このブログ書いてる人のnote記事です。)

これをやると、parse_args は渡された配列をパースしてるのでコマンドラインから渡した引数は全部無視する点には注意してください

まとめ

長くなりましたが、以上がargparseの説明になります。argparseはシンプルに利用することもできますし、多くの引数を活用して細かいカスタマイズもでき、大変柔軟にツールを作ることができます。 自前ツールを作成する際の大変有益な武器になりますので是非触ってみてください。

久々にコマンドラインツールを作った時に、引数の扱い方をド忘れしてしまっていたのでそのメモです。

ちなみに、argpaseっていう大変便利なライブラリもあるのでそれの使い方も紹介したいのですが、それは別記事に回して今回はもっとシンプルな方について書きます。

記事タイトルにも書いていますが、基本的には、sys.argvというのを使います。これは関数ではなく配列型のデータで、argv[0]にそのスクリプトのファイル名、argv[1]以降に、コマンドライン引数が入ります。

参考: sys — システムパラメータと関数 — Python 3.12.1 ドキュメント

また、Pythonのバージョン3.10 以降では、sys.orig_argv ってのも追加されています。

動かしてみるのが一番確認しやすいので、やってみましょう。
次のようなファイルを作ります。ファイル名は sample.py としました。

import sys


for i in range(len(sys.argv)):
    print(f"sys.argv[{i}]=", sys.argv[i])

for i in range(len(sys.orig_argv)):
    print(f"sys.orig_argv[{i}]=", sys.orig_argv[i])

ではこれを実行してみましょう。一旦実行権限とかつけてないので、明示的にpythonコマンドとして実行します。

引数はただ表示するだけで何も処理に影響はないので、デタラメにつけてみました。

% python sample.py file.txt -l aaa.txt "bbb" --name=xyz
sys.argv[0]= sample.py
sys.argv[1]= file.txt
sys.argv[2]= -l
sys.argv[3]= aaa.txt
sys.argv[4]= bbb
sys.argv[5]= --name=xyz
sys.orig_argv[0]= /Users/{macのユーザー名}/.pyenv/versions/3.11.1/bin/python
sys.orig_argv[1]= sample.py
sys.orig_argv[2]= file.txt
sys.orig_argv[3]= -l
sys.orig_argv[4]= aaa.txt
sys.orig_argv[5]= bbb
sys.orig_argv[6]= --name=xyz

はい、argvの方は、[0]番目にスクリプト名が表示され、[1]以降に引数が格納されていましたね。”bbb”のところに注目ですが、ダブルクオーテーションは外されています。

orig_argv は、ファイル名の前の、pythonパスまで入っています。the original command line arguments というから、コマンドに打った python がそのまま入ってるかと思ったらフルパスに展開されています。

ちなみに、sys.argv[0]= sample.py と sys.orig_argv[1]= sample.py のスクリプトパスですが、これはコマンドラインで打ったファイルパスがそのまま表示されいます。どういうことかというと、現在はカレントディレクトリで実行したからこのように表示されているだけで、違う場所から相対パスや絶対パスで指定して実行したらこの中身は変わります。フルパスで指定したらフルパスが入ります。

続いて、このファイルに実行権限をつけて頭のpythonを外してみましょう。
ファイルの先頭に、 #!/usr/bin/env python のシバンを挿入して、chmod で実行権限つけときます。

これで実行すると次のようになります。(引数少し減らしました)

% ./sample.py file.txt -l aaa.txt
sys.argv[0]= ./sample.py
sys.argv[1]= file.txt
sys.argv[2]= -l
sys.argv[3]= aaa.txt
sys.orig_argv[0]= /Users/{macのユーザー名}/.pyenv/versions/3.11.1/bin/python
sys.orig_argv[1]= ./sample.py
sys.orig_argv[2]= file.txt
sys.orig_argv[3]= -l
sys.orig_argv[4]= aaa.txt

sys.argv の方はさっきと大して変わらないですね。カレントディレクトリにパスが通ってないので、スクリプトの指定が./ファイル名 になったのが反映されています。

orig_argv の方も先の結果とほぼ変わりませんが、これはちょっと意外でした。the original command line arguments というから、コマンドラインで python 省略したらここからも省略されると思ってたら相変わらずpython本体のパスが登場しています。

orig_argv が何を意図して追加されたのかがよくわからないのですが、とりあえず argv の方を使っておけば良さそうです。

[1]以降の引数たちは冒頭でも書いたargpaseってライブラリで取ることが多く、sys.argvで取得するのはよほど単純なスクリプトの場合に限られるかなとも思います。となると、もしかしたらargv[0]のファイル名の方がよく使うかもしれないですね。

おまけ

記事の本題と逸れるのですが、argv[0]でスクリプト名が取れるとしたら、スクリプトのフルパスやカレントディレクトリの取得方法が気になるかもしれないのでその撮り方もメモしておきます。

結論を言うと、ファイル名は __file__ って特別な変数に格納されています。
また、カレントディレクトリは os.getcwd() で取れます。

個人的な好みの話ですが、僕は割合の可視化は円グラフより100%積み上げ棒グラフの方が好みです。前職時代はこの種の可視化はもっぱらTableauで描いていたのですが転職して使えなくなってしまったので、Pythonでサクッと描ける方法をメモしておきます。

最初、matplotlibのbarやbarhでbottomやleftを逐一指定してカテゴリごとに描いていく方法を紹介しようと思っていたのですが、pandas.DataFrameのメソッドを使った方が簡単だったのでそちらを先に紹介します。

データの準備

サンプルデータが必要なので適当にDataFrameを作ります。地域ごとに、どの製品が売れているか、みたいなイメージのデータです。

import pandas as pd


data = {
    '製品A': [20, 30, 50],
    '製品B': [60, 70, 40],
    '製品C': [20, 0, 10]
}
df = pd.DataFrame(data, index=["地域1", "地域2", "地域3"])

print(df)
"""
     製品A  製品B  製品C
地域1   20   60   20
地域2   30   70    0
地域3   50   40   10
"""

こちらを地域別に、売れている製品の内訳を可視化していきます。
これも個人的な好みの話ですが、棒グラフは縦向きではなく横向きにします。(ラベルが日本語だとその方が自然な結果になりやすいからです。)

シンプルにやるには、jupyter環境であれば以下のコードだけで目的のグラフが表示されます。

# データの正規化
df_normalized = df.div(df.sum(axis=1), axis=0)
# 横向きの積み上げ棒グラフの描画
df_normalized.plot(kind='barh', stacked=True)

少し脇道に逸れますが、みなさん、この行ごとに正規化するdf.div(df.sum(axis=1), axis=0) ってやり方ご存知でした？昔、「Pandasのデータを割合に変換する」って記事を書いたときは知らなかった方法でした。このときは転置して列ごとに正規化してもう一回転置するという手順を踏んでましたね。

記事の本題に戻ると、plotメソッドのstackedっていう引数をTrueにしておくと棒グラフを積み上げて表示してくれるため目的のグラフになります。
参考: pandas.DataFrame.plot.barh — pandas 2.1.3 documentation

ただ、これで出力すると、グラフの棒の並び順が上から順に地域3, 地域2, 地域1 となり、ちょっと嫌なのと、あとラベル等のカスタマイズも行えた方が良いと思うのでもう少し丁寧なコードも紹介しておきます。(こだわりなければさっきの2行で十分です。)

import matplotlib.pyplot as plt


# データの正規化
df_normalized = df.div(df.sum(axis=1), axis=0)
# 行の順序を逆にする
df_normalized = df_normalized.iloc[::-1]

# 横向きの積み上げ棒グラフの描画
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
df_normalized.plot(kind="barh", stacked=True, ax=ax)

# タイトルとラベルの設定(任意)
ax.set_title("地域別の販売製品割合")
ax.set_ylabel("地域")
ax.set_xlabel("割合")
plt.show()

出来上がりがこちら。

はい、サクッとできましたね。

ちなみに、pandasのplotメソッドを使わない場合は棒グラフの各色の部分の左端の位置を指定しながら順番に描くため、次のコードになります。

import numpy as np


# 積み上げ棒グラフを作成するための基準位置
left = np.zeros(len(df_normalized))

fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)


# 各カテゴリごとにバーを追加
for column in df_normalized.columns:
    ax.barh(
        df_normalized.index,
        df_normalized[column],
        left=left,
        label=column,
        height=0.5
    )
    left += df_normalized[column]

# タイトルとラベルの設定
ax.set_title("地域別の販売製品割合")
ax.set_ylabel('地域')
ax.set_xlabel('割合')

# 凡例の表示
ax.legend()

plt.show()

明らかに面倒なコードでしたね。pandasがいい感じに調整してくれていた、棒の太さ(横向きなので今回の例では縦幅)なども自分で調整する必要があります。
Pythonでコーディングしてるのであれば大抵の場合は元データはDataFrameになっているでしょうし、もしそうでなくても変換は容易なので、Pandasに頼った方を採用しましょう。

いつも使っているJupyterに関する話です。Jupyterの中では、InteractiveShellというオブジェクトが使われていて、これを利用することでさまざまな高度な操作を行うことができます。

多くのことができるのですが、最近これを使って変数を操作する機会があったのでその辺を紹介していきます。

そもそも、InteractiveShell オブジェクトとは？という話ですが、これはIPython（Jupyter Labのバックエンドとして動作する対話的なPythonシェル）の中心的なクラスです。このオブジェクトを通じて、現在のセッションの状態を取得したり、マジックコマンドを実行したりすることができます。

このInteractiveShellはget_ipython() という関数を使うことでそのnotebook自身で取得し、アクセスすることができます。

マジックコマンドの動的な実行

このInteractiveShellオブジェクトを使うと例えば、マジックコマンドをプログラムで実行したりできます。

例えば、%pwd っていう現在のでディレクトリを取得するマジックコマンドがありますが、これが次のようにして実行できます。

ip = get_ipython()
ip.run_line_magic('pwd', '')

(pwdをメソッドとして実行するメリットはなかなか薄いので例が適切でなくてすみません。)

マジックコマンドをそのままマジックコマンドととして実行するのと比べると、結果を文字列として取得して以降の処理で使えるとか、実行するコマンドをプログラムで動的に変化させて実行させることができるといった利点があります。

もちろん、 run_cell_magic() って関数もあるのでセルマジックも実行できますよ。

変数の取得や設定

今回紹介するもう一つの使い方は変数の取得や追加です。名前空間の操作といった方がわかりやすいかもしれません。

InteractiveShellオブジェクトを使って、upyter Labのセッション内で定義されている変数を取得したり、新しい変数を設定したりすることができます。

変数の取得は user_ns 属性を使います。現在定義されている変数が全部辞書として入っています。デフォルトで入ってるJupyter/IPython固有の変数も多いので、自分が宣言した覚えのないものもあると思います。

例えば次のようにして使います。

x = 10
y = "Hello"

ip = get_ipython()
print(ip.user_ns['x'])
# 10
print(ip.user_ns['y'])
# "Hello"

変数の取得としてはあまり用途ないかもしれないですね。どちらかというと既に宣言されているかどうかのチェックとかで使えるかもしれません。

逆に、次のようにして新しい変数を定義することもできます。

ip = get_ipython()
ip.push({'z': "新しい変数"})

print(z)
# 新しい変数

この例だけ見てると、 z=”新しい変数” って直接書くのと同じやん、と思えますが、この方法の利点は変数名の方もプログラムで動的に指定できることです。

まぁ、 exec(“z=’新しい変数'”) ってすれば変数名を動的に変えられるので固有のメリットというわけではないのですが、execは本当にリスクの大きな手段なのでそれよりは安全に使えるという利点があると思います。

この変数の名の指定は、関数やマジックコマンドの中で使うと、例えばユーザーが指定した変数名に結果を格納するといった応用ができます。

最近、別の場所で書いているnote記事の中で自分が自作したマジックコマンドを紹介したのですが、そこで使っているクエリの結果を引数で指定した変数に格納するのにこのip.pushを使っています。

参考: Snowflakeに手軽にSQLを打てるJupyterマジックコマンドを作ってみた

ip.push を使わずに、
ip.user_ns[‘z’] = “新しい変数”
みたいに、user_nsオブジェクトに直接突っ込んでも同じことができますので合わせて覚えておいてください。

InteractiveShellオブジェクトでできることのごく一部ではありますが、今回の記事としては以上となります。

はじめに

ちょっと重い記事が続いてたので今回は小ネタです。

二つの区間データ、要するに数値の区間であれば最小値と最大値のペア、時刻の区間であれば開始時間と終了時間のデータが二つあったときに、それらの区間に重複があるかどうかをサクッと判定するアルゴリズムを紹介します。昔どこかでみた覚えがあるのですが、久々に実装で必要になったときにちょっとだけ悩んだのでそのメモです。

結論

結果だけ知りたい人向けに先に結論だけ書いときます。

区間[s1, e1]と[s2, e2]に重複があるかどうかは次の2つの方法で判定できます。

方法1: $(s1 \leq e2) \land (s2 \leq e1)$ ならばその2区間には重複があります。

方法2: $\max(s1, s2) \leq \max(e1, e2)$ ならばその2区間には重複があります。

両方の方法のメリットデメリットですが、方法1の方が計算速度が速いです。大小比較2個と論理演算ですからmax関数呼び出すより有利ですね。ただ、方法2の方は、応用として、二つの区間に重複があった場合、$\max(e1, e2) – \max(s1, s2)$で重複区間の長さを算出できます。(値が負になったら重複無し。)

(max関数がそんな極端に遅いってこともなく、わざわざ時間測定しない限り体感することもない速度差ではあるのであくまでも参考にどうぞ。)

方法1の導出

これだけで終えてしまうと記事量としてあまりに少ないのと、結論だけ書いてるとすぐ忘れそうなので、真面目に導出を説明しておきます。まずは方法1の方からです。

まず、二つの区間の相対的な位置関係は次の6パターンが存在します。(説明の単純化のため。s1,s2,e1,e2は全て異なる値とします。)

1. 区間1全体が区間2より小さい。つまり、 $s1 < e1 < s2 < e2$.
2. 区間1全体が区間2より大きい。つまり、 $s2 < e2 < s1 < e1$.
3. 区間1が区間2に内包される。つまり、$s2 < s1 < e1 < e2$.
4. 区間1が区間2を内包する。つまり、$s1 < s2 < e2 < e1$.
5. 区間1の後半と区間2の前半が重なる。つまり、$s1 < s2 < e1 < e2$.
6. 区間1の前半と区間2の後半が重なる。つまり、$s2 < s1 < e2 < e1$.

上記の6パターンのうち、1,2 が区間に重複がなく、3,4,5,6は区間に重複があります。

ここで、「そうか！3,4,5,6のどれかを満たすことを確認すればいいんだ！」と判断して4パターンの論理式を書いてorで繋ぐ、としないことがコツです。

区間に重複がないパターンの方が2種類しかなくて判定が単純なんですよね。そして、区間の定義から$s1<e1$と$s2<e2$はわかってるので、あと確認するのはe1とs2, e2とs1の代償関係だけなんですよ。

ということで、 $(e1 < s2) \lor (e2 < s1)$ であれば二つの区間に重複はないと言えます。

これの否定を考えると、ドモルガンの法則から方法1としてあげた式が導かれます。

方法2の導出

続いて方法2の導出です。

これは、先ほど挙げた3,4,5,6の4パターンの不等式を眺めると見つかる法則性なのですが、左の2辺はs1,かs2で、右の2辺はe1, e2なんですよね。

これを素直に数式に落とすと、$\max(s1, s2) \leq \max(e1, e2)$ となります。等号が成立するのはただ1点だけが重複する場合。

そして、区間が重複する場合は、$\max(s1, s2)$は重複区間の開始点であり、$\max(e1, e2)$は重複区間の終了点なので、この2項の差を取ると重複部分の長さも得られます。

まとめ

そんなに難しい計算ではなく、覚えてさえればサクッと実装できる式ではありますが、重複のパターンって4種類あるよなぁと考え始めてしまうと意外に手間取ります。

結果を丸暗記しなくても、区間が重複しないパターンは2個しかなくてそれを否定したら簡単だってことを頭の片隅にでも置いといてもらえると幸いです。

エントロピーや相互情報量の記事が続いていますが、今回の記事で計算の実装方法を紹介して一旦区切りとします。

エントロピーも相互情報量も数式はそこまで難しくないので、numpyで定義通り計算しても良いのですが、エントロピー関係はSciPyに、相互情報量はscikit-learnに用意されているので今回の記事ではそれを使っていきます。

計算対象のデータは、[“a1”, “a2”, “a1”, “a1”, “a2”] みたいにローデータの一覧で保有している場合もあれば、”a1″が3個で”a2″が2個のようにカウント済みのものがある場合もあると思うのでそれぞれ説明していきます。

エントロピーの計算

まず一番基本的なエントロピーの計算からです。これは、scipy.stats.entropy メソッドを使います。
参考: scipy.stats.entropy — SciPy v1.11.3 Manual

基本的な引数はpkなので、確率の一覧を渡すのが想定されていますが、和が1でないなら1になるように正規化してくれるのでサンプルがある場合は個数を渡しても大丈夫です。また、base引数で対数関数の底を指定でき、デフォルトが$e$なので、情報理論で使う場合は$2$を指定しましょう。

やってみます。

import numpy as np  # データ作りに利用
import pandas as pd  # データ作りに利用
from scipy.stats import entropy


pk = np.array([1/2, 1/3, 1/6])  # 確率の一覧が得られた場合。
print(entropy(pk, base=2))
# 1.459147917027245

count_list = np.array([3, 2, 1])  # データの個数の場合
print(entropy(count_list, base=2))
# 1.4591479170272446

# カウント前のデータの一覧がある場合
data_sr = pd.Series(["a1", "a1", "a1", "a2", "a2", "a3"])
# value_counts()で数えあげたものをentropyに渡す
print(entropy(data_sr.value_counts(), base=2))
# 1.4591479170272446

結合エントロピーの計算

次は結合エントロピーです。エントロピーを単純に2次元に拡張したやつですね。(条件付きエントロピーではないので注意してください、

例えば次のような例を考えましょうか。

	b1	b2
a1	4	1
a2	2	3

結合エントロピーの場合はですね、元のカウントデータを2次元から1次元に並び替えて渡します。

matrix_data = np.array([[4, 1], [2, 3]])
print(matrix_data)
"""
[[4 1]
 [2 3]]
"""

# ravel か flattenで1次元化して計算する
print(entropy(matrix_data.ravel(), base=2))
# 1.8464393446710157

# 標本データがある場合
df = pd.DataFrame({
    "A": ["a1", "a1", "a1", "a1", "a1", "a2", "a2", "a2", "a2", "a2"],
    "B": ["b1", "b1", "b1", "b1", "b2", "b1", "b1", "b2", "b2", "b2"],
})

# カウントしたデータを使う
print(df.groupby(["A", "B"]).size())
"""
A   B
a1  b1    4
    b2    1
a2  b1    2
    b2    3
dtype: int64
"""

print(entropy(df.groupby(['A', 'B']).size(), base=2))
# 1.8464393446710157

条件付きエントロピー

次は条件付きエントロピーです。残念なことなのですが、メジャーなライブラリでは条件付きエントロピー専用の関数は提供されていません。

そこで、$H(A|B) = H(A, B) – H(B)$などのエントロピー間の関係式を使って計算することになります。相互情報量も含めて、$H(A|B) = H(A) – I(A; B)$などで計算してもいいのですが、SciPyで完結できるので最初の式のほうが良いでしょう。

先ほどの表データをサンプルとします。$H(B)$については、表データを縦に足し合わせてBだけのカウントデータを作って計算します。

data_B = matrix_data.sum(axis=0)
print(data_B)
# [6 4]

# H(B)の計算
entropy_B = entropy(data_B, base=2)
print(entropy_B)
# 0.9709505944546688

# H(A, B)の計算
joint_entropy = entropy(matrix_data.ravel(), base=2)
print(joint_entropy)
# 1.8464393446710157

# H(A|B) = H(A, B) - H(B)
conditional_entropy_A_given_B = joint_entropy - entropy_B
print(conditional_entropy_A_given_B)
# 0.8754887502163469

# 標本データがある場合
entropy_B = entropy(df["B"].value_counts(), base=2)
joint_entropy = entropy(df.groupby(["A", "B"]).size(), base=2)

conditional_entropy_A_given_B = joint_entropy - entropy_B
print(conditional_entropy_A_given_B)
# 0.8754887502163469

以上で、3種類のエントロピーが計算できました。

相互情報量

最後に相互情報量の計算方法です。

$I(A; B) =H(A)-H(A|B)$など複数の表現方法があるので、ここまでに計算してきた値から算出することもできます。

entropy_A = entropy(df["A"].value_counts(), base=2)
print(entropy(df["A"].value_counts(), base=2) - conditional_entropy_A_given_B)
# 0.12451124978365313

ただ、scikit-learnに専用のメソッドがあるのでこちらの使い方も見ておきましょう。
参考: sklearn.metrics.mutual_info_score — scikit-learn 0.18.2 documentation

引数は、mutual_info_score(labels_true, labels_pred, contingency=None)
となっており、標本データを受け取るのが標準的な使い方で、その第一,第二引数をNoneにしてcontingency引数にカウントデータを渡すこともできます。(contingencyがNoneでない場合はこれが優先されて、先の二つの引数が無視されます。)

1点注意しないといけないのは、entropyと違って対数の底が指定できず、自然対数に固定されてしまうことです。底を$2$で考えたい場合は、$\ln{x}/\ln{2} = \log_2{x}$を使って変換が必要です。

from sklearn.metrics import mutual_info_score


# np.log(2)で割ることを忘れない
# カウントした表データがある場合
print(mutual_info_score(None, None, contingency=matrix_data)/np.log(2))
# 0.12451124978365345

# 標本データがある場合
print(mutual_info_score(df["A"], df["B"])/np.log(2))
# 0.12451124978365345

計算の都合上超軽微な誤差がありますが、それ以外は想定通りの値が得られていますね。

以上で、相互情報量も計算できるようになりました。

自分はあまり使ってこなかったのですが、構造化したデータを保存するのに利用されるYAMLというデータ形式が存在します。いかにもPythonと相性が良さそうなフォーマットです。
参考: YAML – Wikipedia

最近、そこそこ大きいYAMLファイルの総チェックを行う必要があったので、Pythonでスクリプトを書いてチャチャっとやりました。その時のYAMLの読み込み方法のメモと、ついでに逆にデータをYAMLに書き出す方法を残しておきます。

サンプルとしては、次の内容が保存した sample.yaml というファイルを使います。

name: John Doe
age: 30
address:
  city: Tokyo
  country: Japan

必要なライブラリはいくつかありますがPyYAMLが代表的です。インストール時と利用時で名前が違うので注意してください。
参考: PyYAML · PyPI

# インストールコマンド
$ pip install PyYAML

# Pythonでインポートする時
>>> import yaml

それではやっていきましょう。まずはYAMLファイルの読み込みからです。

YAMLを読み込むメソッドは、load(), load_all(), safe_load(), safe_load_all() と4つも用意されています。普通はsafe_load()を使えば良いです。safeがついてない2メソッドにはセキュリティ面のリスクもあるので使用を避けることを推奨します。

import yaml


with open("sample.yaml", "r") as f:
    data = yaml.safe_load(f)

print(data)
# {'name': 'John Doe', 'age': 30, 'address': {'city': 'Tokyo', 'country': 'Japan'}}

簡単ですね。気をつけないといけないのは、safe_load()にファイルパスを渡しても読み込んではくれず、openでファイルを開いて、それを渡すって点です。

YAMLファイルは、”—“というセパレーターを使うことで単一ファイルの中に複数のYAMLオブジェクトを記述できますが、そのようなファイルを読み込むときはsafe_load_all()の方を使うので注意してください。読み込んだ結果はジェネレーターで返ってきます。

ファイルにするが手間だったのでセパレーターで区切られたYAML文字列を直接作ってやって見ました。

yaml_data = yaml.safe_load_all("""name: John Doe
age: 30
---
name: Jane Smith
age: 25
---
name: Bob Johnson
age: 40""".strip())

for y in yaml_data:
    print(y)
"""
{'name': 'John Doe', 'age': 30}
{'name': 'Jane Smith', 'age': 25}
{'name': 'Bob Johnson', 'age': 40}
"""

以上が、YAMLの読み込みです。

次は書き出しです。最初のコード例の、data変数をYAML形式ファイルに保存してみます。

保存には yaml.safe_dump() というメソッドを使います。これにPythonのオブジェクトを渡すとYAMLのテキストに変換してくれます。

print(yaml.safe_dump(data))
"""
address:
  city: Tokyo
  country: Japan
age: 30
name: John Doe
"""

ただし、普通はファイルに保存することになると思いますので、次のように使うことになると思います。

with open("output.yaml", "w") as f:
    yaml.safe_dump(data, f)

これで先ほどのテキストが、”output.yaml”ファイルに保存されます。

YAMLには一般的なブロックスタイルと、JSONに近い見た目のフロースタイルという形式がありますが、defalut_flow_style という引数をTrueにすると、フロースタイルで出力できます。ただ、これは使うことあまり無いかな。

また、インデントの数はindentという引数で指定できます。デフォルトは2ですね。

以上がPythonを用いたYAMLファイルの読み書き方法でした。

Jupyterには便利なマジックコマンド(%や%%を付けて呼び出すアレです)がたくさんありますが、あれを自作する方法を紹介します。

ドキュメントは IPythonのドキュメントのこちらを参照します。
参考: Defining custom magics — IPython 8.14.0 documentation

簡単な方法は、register_line_magic, register_cell_magic, register_line_cell_magic の3種のデコレーターをマジックコマンドとして使いたい関数につけることです。

register_line_magicはその行の文字列を格納する引数を1個だけ、register_cell_magicとregister_line_cell_magicは、マジックコマンドと同じ行の文字列を格納する引数と、セル内の文字列を格納する引数の2個をもちます。

ざっと、受け取った文字列をprintするだけのコマンドを作ってみましょう。3種類それぞれのサンプルです。

from IPython.core.magic import register_line_magic
from IPython.core.magic import register_cell_magic
from IPython.core.magic import register_line_cell_magic


@register_line_magic
def line_magic(line):
    print(line)


@register_cell_magic
def cell_magic(line, cell):
    print(f"line: {line}")
    print(f"cell:\n{cell}")


@register_line_cell_magic
def line_cell_magic(line, cell=None):
    print(f"line: {line}")
    if cell:
        print(f"cell:\n{cell}")

順番に使ってみます。

%line_magic ラインマジックテスト
print()
# 以下出力
# ラインマジックテスト

%%cell_magic セルマジックと同じ行のテキスト
セルマジック内のテキスト
その２行目

＃ 以下出力。
"""
line: セルマジックと同じ行のテキスト
cell:
セルマジック内のテキスト
その２行目
"""

%line_cell_magic ラインマジックとして動作させた場合

# 以下出力
# line: ラインマジックとして動作させた場合

%%line_cell_magic セルマジックとして動作させた場合。
セルの中身

# 以下出力
"""
line: セルマジックとして動作させた場合。
cell:
セルの中身
"""

めっちゃ簡単ですね。

最初のマジックコマンドを定義したコードをPythonファイルとして保存して、import可能なディレクトリに置いておくと、インポートして使うこともできる様になります。例えば、 my_magic.py というファイル名で保存しておけば次の様に使えます。

import my_magic


%line_magic 読み込んだモジュールのマジックコマンドが使える
# 読み込んだモジュールのマジックコマンドが使える

my_magic.line_magic("普通の関数としても呼び出せる")
# 普通の関数としても呼び出せる

さて、通常マジックコマンドをライブラリ等から読み込んで使う場合、この様にimport するのではなく、%load_ext して使うことが多いと思います。これは、先ほどあげたドキュメントのページでベストプラクティスとされているのがその方式だからです。@register_* のデコレーターで直接登録する上記の方法は推奨されてないんですね。

その代わりにどうするかというと、 load_ipython_extension というメソッドを持つpythonファイルを作り、このメソッドの中で定義した関数たちを register_magic_function でマジックコマンドへ登録していきます。

引数は順に、登録する関数本体、コマンドの種類(省略したら’line’)、マジックコマンドとして呼び出す時の名前(省略したら元の関数名)です。

例えば、 my_ext.py というファイルを作りその中を次の様にします。

def load_ipython_extension(ipython):
    ipython.register_magic_function(
        line,
        magic_kind='line',
        magic_name='line_magic'
    )

    ipython.register_magic_function(
        cell,
        magic_kind='cell',
        magic_name='cell_magic'
    )

    ipython.register_magic_function(
        line_cell,
        magic_kind='line_cell',
        magic_name='line_cell_magic'
    )


def line(line):
    print(line)


def cell(line, cell):
    print(f"line: {line}")
    print(f"cell:\n{cell}")


def line_cell(line, cell=None):
    print(f"line: {line}")
    if cell:
        print(f"cell:\n{cell}")

各メソッドそれぞれにはデコレーターはつきません。

この様なファイルを用意すると、load_ext で読み込んだ時に、load_ipython_extension が実行されて、その中でマジックコマンドの登録が行われます。結果、次の様に使えます。

%load_ext my_ext


%line_magic ロードしたマジックコマンドが使えた
# ロードしたマジックコマンドが使えた

先ほどのimportした場合との挙動の違いとしては、これは明示的にマジックコマンドの読み込みだけを行っているので、各メソッドはインポートはされておらず、個々のメソッドの、line, cell, line_cell は名前空間に登録されてないということです。(マジックコマンドとして登録された、line_magic, cell_magic, line_cell_magic の名前でなら通常のメソッドと同じ様に使うことも可能です)

以上が簡単なマジックコマンドの作り方になります。

もう結構前なのですが、メモ化というテクニックを紹介しました。
参考: pythonの関数をメモ化する

これは@functools.lru_cacheというデコレーターを使って、関数の戻り値を記録しておいて何度も同じ関数を実行するコストを削減するのでしたね。計算コストが削減される代わりに、結果を保存しておく分メモリを消費します。

僕はこれを結構使ってたのですが、最近、これをクラスのメソッドに対して利用しているとメモリリークを引き起こすことがあるという気になる情報を得ました。ブログで紹介しちゃった責任もあるので、今回はその問題について調べました。

この問題は何箇所かで指摘されていて、一例を挙げるとこのissueなどがあります。
参考: functools.lru_cache keeps objects alive forever · Issue #64058 · python/cpython

こっちのYoutubeでも話されていますね。
参考: don’t lru_cache methods! (intermediate) anthony explains #382 – YouTube

具体的に説明していくために超単純なクラスを作って実験していきましょう。
まず、そのまま返すメソッドを持ってるだけのシンプルクラスを作ります。そして、このクラスがメモリから解放されたことを確認できる様に、デストラクターが呼び出されたらメッセージを表示する様にしておきます。これをインスタンス化して関数を1回使って、delで破棄します。

class sample1:
    def __init__(self, name):
        self.name = name

    def __del__(self):
        print(f"インスタンス: {self.name} を破棄しました。")

    def identity(self, x):
        return x


a = sample1("a")
print(a.identity(5))
# 5
del a
# インスタンス: a を破棄しました。

デストラクターがちゃんと呼び出されていますね。

これが、メソッドがメモ化されていたらどうなるのかやってみます。

from functools import lru_cache


class sample2:
    def __init__(self, name):
        self.name = name

    def __del__(self):
        print(f"インスタンス: {self.name} を破棄しました。")

    @lru_cache(maxsize=None)
    def identity(self, x):
        return x


b = sample2("b")
print(b.identity(5))
# 5
del b
# 何も表示されない。

今度はガベージコレクターが動きませんでしたね。これは、変数bを削除したことによって変数bからの参照は消えたのですが、メソッドのidentity の一つ目の引数がそのインスタンス自体をとっていて、これを含めてキャッシュしているので、キャッシュがインスタンス自身への参照を保存しているためガベージコレクションの対象にならなかったのです。そのため、インスタンスbが確保していたメモリは解放されず、占拠されたままになります。

ちなみに、循環参照の状態なので、明示的にガベージコレクターを動かすと消えます。

import gc


gc.collect()
# インスタンス: b を破棄しました。

ちなみに、最初のキャッシュが発生した時に循環参照が生まれているのでメモ化したメソッドを一回も使わなかったら普通に消えます。

c = sample2("c")
del c
# インスタンス: c を破棄しました。

以上の様な問題があるので、クラスメソッドで lru_casheを使う時は気をつけて使うことをお勧めします。

とはいえ、最近のMacBookくらいのメモリ量であれば、インスタンスが何個か過剰に残ったとしてそれでメモリが枯渇する様なことはないんじゃないかなとも思います。仮にメモリがピンチになる様な使い方をしていたとしても、maxsizeを適切に設定してメモリサイズを押さえておくとか、明示的にgc/collect()動かすとかの対応が取れるかと。

僕としては、メモリが解放されないことよりも、デストラクターが動かなくてそこに仕込んだ後始末形の処理が動かないのが気になりましたね。