製造業出身のデータサイエンティストがお送りする記事 cProfile でどこの処理で時間がかかっているかを調べられることを知ったので、メモとして残しておきます。 はじめに cProfile は、プログラムのどこで実行にどれだけ時間がかかったのかを提示してくれます。 普段、あまり使うことは無いのですが、処理を早くしたい時などに使える機能だなと思いました。 使い方 基本的な使い方は、-m cProfile を付けて.pyファイルを実行すると使えます。 sample.py def main(N): for i in range(N): pass a() b() def a(): pass def b(): for i in range(N): pass if __name__ == '__main__': N = 10**6 main(N) 実行結果は下記です。 ~/Desktop ❯ python -m cProfile sample.py 7 function calls in 0.074 seconds Ordered by: standard name ncalls tottime percall cumtime percall filename:lineno(function) 1 0.000 0.000 0.074 0.074 sample.py:1(<module>) 1 0.022 0.022 0.074 0.074 sample.py:1(main) 1 0.020 0.020 0.029 0.029 sample.py:13(b) 1 0.000 0.000 0.000 0.000 sample.py:9(a) 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects} 2 0.032 0.016 0.032 0.016 {range} 簡単に各項目の説明を下記に記載します。 項目 概要 ncalls 呼び出し回数 tottime 関数で消費した合計時間 percall tottime をncalls で割った値 cumtime 関数全てのsubfunction に消費された累積時間 percall cumtime をプリミティブな呼び出し回数で割った値 filename:lineno(function) ファイル名、行番号、関数名 さいごに 最後まで読んで頂き、ありがとうございました。 本日は短いですが、cProfile でプロファイリングできることを知りましたのでメモとして残しておきました。

「PuLPで解く線形計画法と連立方程式とナップサック問題」と「PuLPで解く最適割当問題と最適輸送問題」の続編として、今度は最大流問題と最小費用流問題を解いてみます。 インストール !pip install pulp Collecting pulp [?25l Downloading https://files.pythonhosted.org/packages/14/c4/0eec14a0123209c261de6ff154ef3be5cad3fd557c084f468356662e0585/PuLP-2.4-py3-none-any.whl (40.6MB) [K |████████████████████████████████| 40.6MB 76kB/s [?25hCollecting amply>=0.1.2 Downloading https://files.pythonhosted.org/packages/f3/c5/dfa09dd2595a2ab2ab4e6fa7bebef9565812722e1980d04b0edce5032066/amply-0.1.4-py3-none-any.whl Requirement already satisfied: pyparsing in /usr/local/lib/python3.7/dist-packages (from amply>=0.1.2->pulp) (2.4.7) Requirement already satisfied: docutils>=0.3 in /usr/local/lib/python3.7/dist-packages (from amply>=0.1.2->pulp) (0.17.1) Installing collected packages: amply, pulp Successfully installed amply-0.1.4 pulp-2.4 必要なライブラリのインポート import pulp import numpy as np import matplotlib.pyplot as plt 問題設定 座標 (0, 0) から (7, 7) の 8 x 8 のグリッド状の頂点（ノード）を考えます。これらの頂点間に、次のように縦・横・斜めの辺（エッジ）が存在するものとします。 grid_size = 8 plt.figure(figsize=(12,12)) checkpoints = [] for x1 in range(grid_size): for y1 in range(grid_size): if x1 == y1: pass checkpoints.append((x1, y1)) plt.scatter([x1], [y1]) edges = {} for p1 in checkpoints: for p2 in checkpoints: if p1[0] == p2[0] and abs(p1[1] - p2[1]) == 1: plt.plot([p1[0], p2[0]], [p1[1], p2[1]]) elif p1[1] == p2[1] and abs(p1[0] - p2[0]) == 1: plt.plot([p1[0], p2[0]], [p1[1], p2[1]]) elif p1[1] - p2[1] == 1 and p1[0] - p2[0] == 1: plt.plot([p1[0], p2[0]], [p1[1], p2[1]]) else: continue if p1 not in edges.keys(): edges[p1] = {} if p2 not in edges.keys(): edges[p2] = {} weight = int(np.random.rand() * 50) if p2 not in edges[p1]: edges[p1][p2] = weight if p1 not in edges[p2]: edges[p2][p1] = weight plt.text((p1[0] + p2[0]) / 2, (p1[1] + p2[1]) / 2, edges[p1][p2]) 上の図で、辺に書かれている数字は、その辺の容量（その辺に沿って流すことができる「品物」の最大量）とします。 最大流問題 最大流問題（Maximum flow problem）とは、供給点（ソース) から需要点（シンク）まで流すことができる「品物」の最大量を求める問題です。 ここでは、原点 (0, 0) を供給点（ソース) とし、点 (7, 7) を需要点（シンク）とします。 まず PuLP を使って、問題設定するためのインスタンスを作成しましょう。 import pulp prob = pulp.LpProblem('max_flow_problem', sense = pulp.LpMaximize) 辺ごとの重み（品物の輸送量）の最小値がゼロ、最大値がその辺の容量となるようにします。 val = {} for k1, v1 in edges.items(): val[k1] = {} for k2, v2 in v1.items(): val[k1][k2] = pulp.LpVariable("({},{})-({},{})".format( k1[0], k1[1], k2[0], k2[1] ), lowBound = 0, upBound = v2) 最大化すべき値は、供給点（ソース) から需要点（シンク）まで流すことができる「品物」の最大量ですが、これは「供給点 (0, 0) から出てくる品物の総和」を最大化することと同義です。 prob += pulp.lpSum([v for v in val[(0, 0)].values()]) 「供給点（ソース) と需要点（シンク）以外の頂点（ノード）において、入ってくる品物の量の総和と出ていく品物の量の総和は同じである」という制約条件を加えます。 ここで、ノード (i, j) は (i - 1, j), (i, j - 1), (i - 1, j - 1) から品物が入ってきて、(i + 1, j), (i, j + 1), (i + 1, j + 1) に品物が出ていくものとし、逆流はしないものとします。 for x, y in checkpoints: if x == grid_size - 1 and y == grid_size - 1: continue elif x == 0 and y == 0: continue sahen = 0 if (x, y - 1) in val[(x, y)]: sahen += val[(x, y)][(x, y - 1)] if (x - 1, y) in val[(x, y)]: sahen += val[(x, y)][(x - 1, y)] if (x - 1, y - 1) in val[(x, y)]: sahen += val[(x, y)][(x - 1, y - 1)] uhen = 0 if (x, y + 1) in val[(x, y)]: uhen += val[(x, y)][(x, y + 1)] if (x + 1, y) in val[(x, y)]: uhen += val[(x, y)][(x + 1, y)] if (x + 1, y + 1) in val[(x, y)]: uhen += val[(x, y)][(x + 1, y + 1)] prob += sahen == uhen 今回の問題設定では、実装を簡単にするため、頂点 k1 から k2 に向かう辺と、 k2 から k1 に向かう辺の流量を同じであるとしておきます。 for k1, v1 in val.items(): for k2, v2 in v1.items(): prob += val[k1][k2] == val[k2][k1] これで目的関数の入力と制約条件の入力が完了したので、ソルバーを実行します。 result = prob.solve() pulp.LpStatus[result] 'Optimal' 結果は次のようになりました。 plt.figure(figsize=(12,12)) plt.title("Optimal value = {}".format(prob.objective.value())) for k1, v1 in val.items(): for k2, v2 in v1.items(): if v2.value() > 0: plt.plot([k1[0], k2[0]], [k1[1], k2[1]], linewidth=int(v2.value()), alpha=0.5) #print(k1, k2, v2.value()) plt.text( (k1[0] + k2[0]) / 2, (k1[1] + k2[1]) / 2, "{} / {}".format(int(v2.value()), edges[k1][k2]) ) for x1 in range(grid_size): for y1 in range(grid_size): plt.scatter([x1], [y1]) 品物の流れがゼロになった辺は省略し、たとえば容量が 33 の辺に品物が 2 だけ流れた時は 2 / 33 のように表示しています。流れの多い辺は太く表示しています。 最小費用流問題 最小費用流問題 (Minimum cost flow problem) は、最大流問題が最大の流量を求めるのに対して、流量を固定したときにそのコストの総量を最小化する問題です。 ここでは、上の最大流問題で求まった最大流量で固定することにしましょう。 total_flow = prob.objective.value() 今回新たに加えるのは、各辺において 1 の量の品物を流すためのコストです。 costs = {} for k1, v1 in edges.items(): for k2, v2 in v1.items(): cost = int(np.random.rand() * 100) if k1 not in costs.keys(): costs[k1] = {} if k2 not in costs.keys(): costs[k2] = {} costs[k1][k2] = cost costs[k2][k1] = cost 次のように、最大流問題のときとほぼ同じコードで解けます。違う点は、目的関数が「辺のコスト」と「流量」の積和であるという点です。 import pulp prob = pulp.LpProblem('minimum_flow_problem', sense = pulp.LpMinimize) sum_cost = 0 for k1, v1 in val.items(): for k2, v2 in v1.items(): sum_cost = costs[k1][k2] * v2 prob += sum_cost prob += pulp.lpSum([v for v in val[(0, 0)].values()]) == total_flow for x, y in checkpoints: if x == grid_size - 1 and y == grid_size - 1: continue elif x == 0 and y == 0: continue sahen = 0 if (x, y - 1) in val[(x, y)]: sahen += val[(x, y)][(x, y - 1)] if (x - 1, y) in val[(x, y)]: sahen += val[(x, y)][(x - 1, y)] if (x - 1, y - 1) in val[(x, y)]: sahen += val[(x, y)][(x - 1, y - 1)] uhen = 0 if (x, y + 1) in val[(x, y)]: uhen += val[(x, y)][(x, y + 1)] if (x + 1, y) in val[(x, y)]: uhen += val[(x, y)][(x + 1, y)] if (x + 1, y + 1) in val[(x, y)]: uhen += val[(x, y)][(x + 1, y + 1)] prob += sahen == uhen for k1, v1 in val.items(): for k2, v2 in v1.items(): prob += val[k1][k2] == val[k2][k1] result = prob.solve() pulp.LpStatus[result] 'Optimal' 結果は次のようになりました。 plt.figure(figsize=(12,12)) plt.title("Optimal value = {}".format(prob.objective.value())) for k1, v1 in val.items(): for k2, v2 in v1.items(): if v2.value() > 0: plt.plot([k1[0], k2[0]], [k1[1], k2[1]], linewidth=int(v2.value()), alpha=0.5) #print(k1, k2, v2.value()) plt.text((k1[0] + k2[0]) / 2, (k1[1] + k2[1]) / 2, "{} / {}".format(int(v2.value()), edges[k1][k2])) for x1 in range(grid_size): for y1 in range(grid_size): plt.scatter([x1], [y1]) 最大流問題で得られた結果と、最小費用流問題で得られた結果は、流量は同じですが微妙に異なりますね。どちらも、解が複数ある場合はその中の１つしか出力してくれないと思います。 最小費用流問題の総流量を別の数字に変更すると当然違う結果が得られます。試してみてください。

はじめ Command Prompt, PowerShell, Terminalなどで色付き文字を出力する方法(ライブラリ)は多数あります。 有名どころで言えば Colorama from colorama import Fore print(Fore.GREEN + "This text is green!") termcolor from termcolor import colored print(colored('This text is green!', 'green')) とかですかね。 素晴らしいライブラリ これらは実に便利なライブラリです。簡単に色づき文字を作成したりすることができます。 ただ、実は一つ問題があります。 ネストできない これはつまり、たとえば (赤字の文章)(緑字の文章)(赤字の文章) とか (赤字の文章)(赤字で下線のある文章)(赤字の文章) などを出力しようとすると、結構めんどくさいということです。 以下の画像における各行のようなものです。 たとえばこれをtermcolorで実装しようとすると、以下のようになるでしょう。 from termcolor import colored print(colored('[RED]', 'red') + colored('[GREEN]', 'green') + colored('[RED]', 'red')) print(colored('[RED]', 'red') + colored('[UNDERLINE]', 'red', attrs=['underline']) + colored('[RED]', 'red')) ネストが深くなったり、長い文章同士を結合しようとすると毎回こういった指定をしなければならず、ちょっと大変です。そこで スタイリングライブラリ iro iroを使ってみましょう！ $ pip install iro GitHub さて、このライブラリを使うと、ネストされたスタイルをこのように表現することができます。 出力は上記の画像と全く同じです。 from iro import Iro, Color, Style print(Iro([ Color.RED, "[RED]", [ Color.GREEN, "[GREEN]" ], "[RED]" ])) print(Iro([ Color.RED, "[RED]", [ Style.UNDERLINE, "[UNDERLINE]" ], "[RED]" ])) ある深さにあるstrには、その深さまでの色・スタイルがすべて適用されます。深さが深い方が優先度が高くなるので、このように簡単に表すことができます。 また、これはネストがより深くなればなるほどわかりやすくなります。 以下の二つのコードは同等の出力をします。 print(Iro([ Color.RED, "[RED]", [ Style.UNDERLINE, "[RED/UNDERLINE]", [ Style.BOLD, Color.GREEN, "[GREEN/UNDERLINE/BOLD]", [ Style.INVERT, "[GREEN/UNDERLINE/BOLD/INVERT]", [ Style.OFF_BOLD, "[GREEN/UNDERLINE/INVERT]" ] ], "[GREEN/UNDERLINE/BOLD]" ], "[RED/UNDERLINE]" ] ])) print(colored("[RED]", 'red'), colored("[RED/UNDERLINE]", 'red', attrs=['underline']), colored("[GREEN/UNDERLINE/BOLD]", 'green', attrs=['underline', 'bold']), colored("[GREEN/UNDERLINE/BOLD/ITALIC]", 'green', attrs=['underline', 'bold', 'reverse']), colored("[GREEN/UNDERLINE/ITALIC]", 'green', attrs=['underline', 'reverse']), colored("[GREEN/UNDERLINE/BOLD]", 'green', attrs=['underline', 'bold']), colored("[RED/UNDERLINE]", 'red', attrs=['underline']), sep='' ) Iroのほうはあえて見やすさのために大量に改行などしてますが、非常に分かりやすくなっているのではないでしょうか？ これだけでなく、Iroは8-bitカラーや、RGBカラーを使用することもできます。 print(Iro([ [ Color256(135), Color256(217, bg=True), " SEXY COLOR " ], [ ColorRGB.from_color_code("EEA47F"), ColorRGB(0, 83, 156, bg=True), " AND THIS IS SICK " ], [ Color.BRIGHT_YELLOW, Color.BG_YELLOW, " BRIGHT!! " ] ], disable_rgb=False)) デフォルトではdisable_rgbはTrueで、ColorRGBをそれに最も近しいColor256に変換して出力します。 RGBカラーに対応していないコンソールも時々あるので、その時はdisable_rgbをTrueのままにしておきましょう。Trueにするとこのようになります。 確かによく見てみると、上の写真と下の写真、中央の色が違いますね！(その程度の差なので、結構精度は高い) Styleも豊富で RESET BOLD DIM ITALIC UNDERLINE SLOW_BLINK RAPID_BLINK INVERT HIDE STRIKE BLACKLETTER_FONT DOUBLY_UNDERLINE OFF_INTENSITY OFF_BOLD OFF_DIM OFF_ITALIC OFF_UNDERLINE OFF_BLINK OFF_INVERT OFF_HIDE OFF_STRIKE OFF_COLOR OFF_BG_COLOR OFF_OVERLINE などが使えます。 Iro([text...], optimize_level=1) とすることで最適化が行われ描画が高速になる可能性がありますが、端末によっては最適化によって動作が不安定になる場合があります。ただ、この効果は大きく、GitHubのQ&Aを見ると value = [ Color.RED, "[RED]", [ Style.UNDERLINE, "[RED/UNDERLINE]", [ Style.BOLD, Color.GREEN, "[GREEN/UNDERLINE/BOLD]", [ Style.INVERT, "[GREEN/UNDERLINE/BOLD/INVERT]", [ Style.OFF_BOLD, "[GREEN/UNDERLINE/INVERT]" ] ], "[GREEN/UNDERLINE/BOLD]" ], "[RED/UNDERLINE]" ] ] print(repr(Iro(value, optimize_level=1).text)) # '\x1b[31m[RED]\x1b[4m[RED/UNDERLINE]\x1b[1m\x1b[32m[GREEN/UNDERLINE/BOLD]\x1b[7m[GREEN/UNDERLINE/BOLD/INVERT] # \x1b[22m[GREEN/UNDERLINE/INVERT]\x1b[1m\x1b[27m[GREEN/UNDERLINE/BOLD]\x1b[22m\x1b[31m[RED/UNDERLINE] # \x1b[24m\x1b[39m\x1b[0m' print(repr(Iro(value, optimize_level=0).text)) # '\x1b[31m[RED]\x1b[31m\x1b[4m[RED/UNDERLINE]\x1b[31m\x1b[4m\x1b[1m\x1b[32m[GREEN/UNDERLINE/BOLD] # \x1b[31m\x1b[4m\x1b[1m\x1b[32m\x1b[7m[GREEN/UNDERLINE/BOLD/INVERT] # \x1b[31m\x1b[4m\x1b[1m\x1b[32m\x1b[7m\x1b[22m[GREEN/UNDERLINE/INVERT] # \x1b[39m\x1b[24m\x1b[22m\x1b[39m\x1b[27m\x1b[22m\x1b[31m\x1b[4m\x1b[1m\x1b[32m\x1b[7m\x1b[39m\x1b # [24m\x1b[22m\x1b[39m\x1b[27m\x1b[31m\x1b[4m\x1b[1m\x1b[32m[GREEN/UNDERLINE/BOLD] # \x1b[39m\x1b[24m\x1b[22m\x1b[39m\x1b[31m\x1b[4m[RED/UNDERLINE] # \x1b[39m\x1b[24m\x1b[31m\x1b[39m\x1b[0m' という程の変化があるようです。 おわり 是非便利にお使いください!(作者)

概要 scikit-learn の 分類用の決定木 DecisionTreeClassifier の各リーフノードでは、クラス確率が一番高いクラスに分類されます。 ではクラス確率が一番高いクラスが複数ある（同率首位）リーフノードではどのクラスに分類されるのか？ を確認しました。 （そもそもこのようなリーフノードは良くありませんが） 確認結果 クラス確率が一番高いクラスが複数あるリーフノードでは、 クラス確率が一番高いクラスのうち、クラスインデックスが一番小さいクラスに分類されます。 確認過程 実行環境 OS: Windows 10 Python 3.9.5 sklearn 0.24.2 matplotlib 3.4.2 公式ページの説明 sklearn.tree.DecisionTreeClassifier The predict method operates using the numpy.argmax function on the outputs of predict_proba. This means that in case the highest predicted probabilities are tied, the classifier will predict the tied class with the lowest index in classes_. 引用文の predict_proba は、各クラス確率を出力するメソッドです。（クラスAとクラスBの2値分類の問題の場合、[クラスAである確率, クラスBである確率]を出力） ※ predict_proba() メモ（DecisionTreeClassifier, RandomForestClassifier 対象） 確率値MAXのクラスが複数ある場合、各クラスに割り当てられたインデックスが一番小さいクラスに分類される、ということですね。 サンプルデータで確認 公式ページの説明で分かったつもりではありますが、念のためデータを見ながら確認しておきます。 お馴染みのirisデータセットを使います。 1. 決定木1: 生成とプロット まず学習まで。 from sklearn.datasets import load_iris from sklearn.tree import DecisionTreeClassifier, plot_tree iris = load_iris() # データ確認 print(iris.data[:5]) # 各説明変数の値（先頭5件） # -> [[5.1 3.5 1.4 0.2] # -> [4.9 3. 1.4 0.2] # -> [4.7 3.2 1.3 0.2] # -> [4.6 3.1 1.5 0.2] # -> [5. 3.6 1.4 0.2]] print(iris.feature_names) # 各説明変数名 # -> ['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)'] print(iris.target[:5]) # 正解ラベルの値（先頭5件） # -> [0 0 0 0 0] print(set(iris.target)) # 正解ラベルのユニーク値 # -> {0, 1, 2} print(iris.target_names) # -> ['setosa' 'versicolor' 'virginica'] ★ 0=setosa, 1=versicolor, 2=virginica # 学習 X, y = iris.data, iris.target tree_clf = DecisionTreeClassifier(max_depth=1, random_state=42) # 確認目的の都合で最大深度1 tree_clf.fit(X, y) 次に、クラスインデックス、クラス値の対応 を見ておきます。 for idx, class_val in enumerate(tree_clf.classes_): print(f"クラスインデックス: {idx}　クラス値: {class_val}") # -> クラスインデックス: 0　クラス値: 0 # -> クラスインデックス: 1　クラス値: 1 # -> クラスインデックス: 2　クラス値: 2 生成された決定木をプロットします。 ※ plot_tree()はsklearn Vr.0.21 で追加された機能です。 plot_tree( tree_clf, feature_names=iris.feature_names, class_names=iris.target_names, impurity=False, # 不純度は今回の話には不要なので表示しない node_ids=True # 各ノードのIDを表示 ) node #2 を見るとクラス1(versicolor), 2(virginica)のサンプル数は同じ（-> ノード内クラス確率が同じ）ですが、クラス1(versicolor)に分類するようです。 （クラスインデックスが小さい方 = クラス1 なので、ドキュメントの記載通り） 2. 決定木1: テストデータの各クラス確率と分類結果 node #2 にやってくるテストデータx_testを作り、各クラス確率と分類結果を確認します。 # テストデータ x_test = [[5.1, 3.5, 2.5, 1.7]] # petal length(cm)=3つ目の説明変数 > 2.45にする print(tree_clf.predict_proba(x_test)) # -> [[0. 0.5 0.5]] # [0(setosa)の確率, 1(versicolor)の確率, 2(virginica)の確率]　★1と2の確率が同じ print(tree_clf.predict(x_test)) # -> [1] # クラス1(versicolor)に分類された 3. 決定木2: 生成とプロット ここで、クラスインデックスが入れ替わるようにクラス値1を3に置換し、y2を作成します。 y2を正解ラベルとして新しいモデルを作成します。 y2 = y.copy() # クラス値1(versicolor) の値を 3 にする y2[y2==1] = 3 print(set(y2)) # -> {0, 2, 3} # 上記置換に合わせてクラス名のリストを作成 class_names = ["setosa", "virginica", "versicolor"] # 0=setosa, 2=virginica, 3=versicolor tree_clf2 = DecisionTreeClassifier(max_depth=1, random_state=42) tree_clf2.fit(X, y2) 新モデルもクラスインデックス、クラス値の対応を見ておきます。 # クラスインデックス、クラス値の対応 for idx, class_val in enumerate(tree_clf2.classes_): print(f"クラスインデックス: {idx}　クラス値: {class_val}") # -> クラスインデックス: 0　クラス値: 0 # -> クラスインデックス: 1　クラス値: 2 # -> クラスインデックス: 2　クラス値: 3 生成された決定木をプロットします。 plot_tree( tree_clf2, feature_names=iris.feature_names, class_names=class_names, impurity=False, node_ids=True ) 4. 決定木2: テストデータの各クラス確率と分類結果 上記3で使用したテストデータx_testの、新しいモデルによる各クラス確率と分類結果を確認します。 print(tree_clf2.predict_proba(x_test)) # -> [[0. 0.5 0.5]] # [0(setosa)の確率, 2(virginica)の確率, 3(versicolor)の確率] print(tree_clf2.predict(x_test)) # -> [2] # クラス2(virginica)に分類された クラスインデックスを変えただけで分類結果が変わりました。 感想 そもそもこのようなリーフノードでは、どちらのクラスに分類すべき、と言えません。 分類結果だけでなく決定木の状態確認も重要ですね。

概要 scikit-learn の RandomForestClassifier のメソッドpredict_proba()は各クラス確率を出力します。このクラス確率とは具体的に何か、メモを実行結果と共に残します。 まず、1本の決定木であるDecisionTreeClassifierのpredict_proba()を理解し、その後、複数の決定木で構成される RandomForestClassifierのpredict_proba()を確認しました。 確認結果 DecisionTreeClassifierのpredict_proba() \begin{align} クラスAの確率値&=\frac{そのリーフノード内のクラスAのサンプル数}{そのリーフノード内のすべてのサンプル数},\\ \\ クラスBの確率値&=\frac{そのリーフノード内のクラスBのサンプル数}{そのリーフノード内のすべてのサンプル数},\ \\ クラスC・・・ \end{align} RandomForestClassifierのpredict_proba() モデル内の各決定木の各クラス確率の平均 確認過程 実行環境 OS: Windows 10 Python 3.9.5 sklearn 0.24.2 matplotlib 3.4.2 DecisionTreeClassifier のpredict_proba() まず、決定木（分類用）1本の predict_proba() を確認します。 公式ページより sklearn.tree.DecisionTreeClassifier The predicted class probability is the fraction of samples of the same class in a leaf. ふむふむ、決定木のpredict_proba()は、そのリーフノード内のサンプルに対するそのクラスのサンプルの割合を表しているのですね。 サンプルデータで確認 お馴染みのirisデータセットを使います。 from sklearn.datasets import load_iris from sklearn.tree import DecisionTreeClassifier, plot_tree import matplotlib.pyplot as plt iris = load_iris() # データ確認 print(iris.data[:5]) # 各説明変数の値（先頭5件） # -> [[5.1 3.5 1.4 0.2] # -> [4.9 3. 1.4 0.2] # -> [4.7 3.2 1.3 0.2] # -> [4.6 3.1 1.5 0.2] # -> [5. 3.6 1.4 0.2]] print(iris.feature_names) # 各説明変数名 # -> ['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)'] print(iris.target[:5]) # 正解ラベルの値（先頭5件） # -> [0 0 0 0 0] print(set(iris.target)) # 正解ラベルのユニーク値 # -> {0, 1, 2} print(iris.target_names) # -> ['setosa' 'versicolor' 'virginica'] ★ 0=setosa, 1=versicolor, 2=virginica # 学習 X, y = iris.data, iris.target # 今回の目的のため、木の最大深さを2に設定してシンプルなモデルに。 tree_clf = DecisionTreeClassifier(max_depth=2, random_state=42) tree_clf.fit(X, y) 決定木をプロットします。 ここでは、最初にimportしておいた sklearn.tree.plot_tree()を使っています。 ※ sklearn Vr.0.21で追加されたプロット機能です。 plot_tree( tree_clf, feature_names=iris.feature_names, class_names=iris.target_names, impurity=False, # 不純度は今回の話には不要なので表示しない node_ids=True, # 各ノードのIDを表示 filled=True, # ノードを各クラスのサンプル割合に応じて色付け ) 【プロット1】 各ノード内の値の補足: - node #X: ノードID - samples: そのノードのサンプル数 - value : [クラス0(satosa)のサンプル数, クラス1(versicolor)のサンプル数、クラス2(virginica)のサンプル数] - classes: そのノードのクラス samples と value に関しては、以下のようなことがが分かります。 全サンプル数は150 ・・・node#0 の samples より node#3 のサンプル数は54 ・・・node #3 の samples より node#3 のサンプルの内訳は以下の通り・・・node#3 の value より クラス0(setosa)が0 クラス1(versicolor)が49 クラス2(virginica)が5 ここで、samples を%表示, value を比率で表示するように変更します。 # プロットの描画サイズ調整（今回は小数第8位までの出力するため見やすいサイズに） ax = plt.figure(figsize=(12, 6)).add_subplot() plot_tree( tree_clf, feature_names=iris.feature_names, class_names=iris.target_names, impurity=False, node_ids=True, filled=True, proportion=True, # 各ノードのsamplesを%表示, value を各クラスの比率で表示する precision=8, # 小数第8位まで表示 ax=ax # 事前準備しておいたaxを適用 ) 【プロット2】 オプション追加により、samples と value の表示が変わりました。 それぞれ、以下の通りであることが確認できます。 samples: そのノードのサンプル数が全サンプル数に占める割合(%) node#3 の場合: node#3内のサンプル数54 ÷ 全サンプル数150 * 100 = 36.0(%) value : ノード内の[クラス0(satosa)の比率, クラス1(versicolor)の比率、クラス2(virginica)の比率] node#3 の場合: クラス0の比率: node#3のクラス0のサンプル数0 ÷ node#3のサンプル数54 = 0.0 クラス1の比率: node#3のクラス1のサンプル数49 ÷ node#3のサンプル数54 = 0.90740741 クラス2の比率: node#3のクラス2のサンプル数5 ÷ node#3のサンプル数54 = 0.09259259 テストデータを作成し、predict_proba() の出力（各クラス確率）を確認します。 【プロット2】で確認した各クラスの比率ですね。 # 【テストデータA】node#3に入る tree_clf.predict_proba([[5.1, 3.5, 2.5, 1.7]]) # -> array([[0. , 0.90740741, 0.09259259]]) # 【テストデータB】node#1に入る tree_clf.predict_proba([[5.1, 3.5, 1.4, 0.2]]) # -> array([[1., 0., 0.]]) RandomForestClassifier のpredict_proba() 次に複数の決定木を並列に扱うランダムフォレスト（分類用）の predict_proba() を確認します。 公式ページより sklearn.ensemble.RandomForestClassifier The predicted class probabilities of an input sample are computed as the mean predicted class probabilities of the trees in the forest. ふむふむ、モデル内の決定木1本1本が出してくるクラス確率の平均だそうです。 サンプルデータで確認 引き続きirisデータセットで確認します。 まず、ランダムフォレストのモデルを作成します。 木の数2、と小さい値を指定しました(オプション n_estimators=2) from sklearn.ensemble import RandomForestClassifier rf_clf = RandomForestClassifier(n_estimators=2, max_depth=2, random_state=42) rf_clf.fit(X, y) 次に、決定木の確認で使用した【テストデータA】を使い ランダムフォレストの predict_proba() の出力（各クラス確率）を確認します。 rf_clf.predict_proba([[5.1, 3.5, 2.5, 1.7]]) # -> array([[0. , 0.48039216, 0.51960784]]) RandomForestClassifier の predict_proba() は、決定木1本1本が出してくるクラス確率の平均 ということなので、決定木ごとにクラス確率を見てみます。 print(f"ランダムフォレスト内の決定木の数: {len(rf_clf.estimators_)}", end="\n\n") fig = plt.figure(figsize=(25, 6)) for idx in range(len((rf_clf.estimators_))): print(f"tree_{idx}: サンプル数 {rf_clf.estimators_[idx].tree_.n_node_samples[0]}") print(f"tree_{idx}のpredict_proba()の出力: {rf_clf.estimators_[idx].predict_proba([[5.1, 3.5, 2.5, 1.7]])}", end="\n\n") ax = fig.add_subplot(1, 2, idx+1) plot_tree( rf_clf.estimators_[idx], feature_names=iris.feature_names, class_names=iris.target_names, impurity=False, node_ids=True, filled=True, proportion=True, precision=8, ax=ax, ) # -> ランダムフォレスト内の決定木の数: 2 # # -> tree_0: サンプル数 101 # -> tree_0のpredict_proba()の出力: [[0. 0.96078431 0.03921569]] # # -> tree_1: サンプル数 97 # -> tree_1のpredict_proba()の出力: [[0. 0. 1.]] 2つの決定木の predict_proba(【テストデータA】) の出力=クラス確率 を上記コード内のコメントに記載しました。 クラス毎に平均を計算してみます. クラス0(setosa) クラス1(versicolor) クラス2(virginica) tree_0 のクラス確率 0.0 0.96078431 0.03921569 tree_1 のクラス確率 0.0 0.0 1.0 平均 0.0 0.480392155 0.519607845 ここで、ランダムフォレストの predict_proba(【テストデータA】) の出力を再掲。 rf_clf.predict_proba([[5.1, 3.5, 2.5, 1.7]]) # -> array([[0. , 0.48039216, 0.51960784]]) 微細な誤差はありますが、上記で計算した各クラス確率の平均とほぼ一致しました。 データを見ながら具体的に predict_proba() について確認しました。

最近PythonのスクレイピングライブラリであるBeautifulSoupを勉強しているとかなり面白くて、何かに活用してみたくなったので、全上場銘柄の株価を取得するツールを作ってみました。取得したデータを使ってどういう分析をするかはまだ検討中ですが、色々なことができそうです。 概要 pandasを用いたデータの入出力 まず、こちらから東証上場銘柄一覧が記載されたCSVを保存します。私は今回ETFのデータは不要だったので事前に消しておきました。 pandasでは構造化データが扱いやすいDataFrameや、簡単にファイルの入出力ができるread_csv, to_csvメソッドが用意されており、非常に便利です。 requests, BeautifulSoupによるHTML取得・解析 requestsモジュールで指定したURLからHTMLを取得し、BeautifulSoupによって必要なデータを抽出します。 こちらに各銘柄の年度別株価データが載っているのですが、銘柄・年度ごとにページが分かれているので全銘柄の20年分のデータを手作業で集めるのはまず無理です。 そこで、事前に保存したコード一覧を読み込み、各コードに対して年度ごとに異なるURLにアクセスしてデータを取得します。 今回アクセスするページ数、取得するデータ量は非常に多いため、この自動化ツールを使っても私の環境では半日以上かかりました。笑 マルチプロセス等を使えばもう少し速くできるかもしれません。 ソースコード # ライブラリをインポート from bs4 import BeautifulSoup import requests import time import os import pandas as pd if __name__ == '__main__': # 開始時刻 start_time = time.time() # 銘柄・コード一覧を取得 df_input = pd.read_csv('data_j.csv', encoding='shift_jis') # 調べたいデータを指定 codes = df_input['コード'].to_list() # 20年間の株価を取得する years = range(2000, 2020) # 各銘柄（コード）に対して for code in codes: df = pd.DataFrame(columns=['日付', '終値']) # 日付と終値を格納するDataFrame row_start = 1 # 各年の最初のデータの行 flg_print = False file_name = '' # 保存するファイル名 # 各年に対して for year in years: try: # URLを取得 url = "https://kabuoji3.com/stock/" + str(code) + "/" + str(year) + "/" # headers情報は確認くん（https://www.ugtop.com/spill.shtml）で確認 # headers = {"User-Agent": "Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_4) AppleWebKit/605.1.15 (KHTML, like Gecko) Version/13.1 Safari/605.1.15 "} headers = { "User-Agent": "Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/89.0.4389.128 Safari/537.36" } # HTML解析用オブジェクト soup = BeautifulSoup(requests.get(url, headers=headers).content, 'html.parser') # タイトルを取得 title = soup.select_one("span.jp").text # タイトルを出力 if not flg_print: print(title) flg_print = True # ファイル名を指定 file_name = '{0}.csv'.format(str(title)) # 既にファイルが存在する場合はスキップ if os.path.exists(file_name): break # 株価を取り出す all_tr = soup.find_all('tr') # 先頭要素はヘッダなので削除 all_tr.remove(all_tr[0]) # その年のデータの最終行を取得 row_finish = row_start + len(all_tr) # 各行に日付・終値のデータを書き込む for i in range(row_start, row_finish): tr = all_tr[i - row_start].find_all('td') # リスト型 date = tr[0].text # 日付 fp = tr[4].text # 終値 # DataFrameの末尾に追加 df = df.append({'日付': date, '終値': fp}, ignore_index=True) # 次の年のデータの開始行を更新 row_start = row_finish except: # その年のデータが存在しなくてもエラー復帰せずに次の年のデータの取得に移る pass # データを保存 df.to_csv('data\\' + file_name, encoding='shift_jis', index=False) stop_time = time.time() # 終了時刻 elapsed = stop_time - start_time # 経過時間 # 経過時間の出力 print('処理が完了しました。経過時間: {0:.1f}秒'.format(elapsed))

はじめに 自分が通っている大学の工学部では、1年生時に多くの数学的処理に関して学びます。けれど、実際どのように役立っているのか不明なままそれらの処理を行っているのが実情だと感じています。 そして多くの学生は、これらの数学的処理がつながる先を考えずに問題を解ければいいと考えているかもしれません。けれど、学年が上がるにつれ1年生時に習った数学的処理が実世界を記述するのに大きな意味を持ち、それが分かったときに面白さが生じると思っています。今回はそんな面白さをを共有できたらと願いこの記事を執筆します。 1.凡授業で扱う一部の内容 工学部の授業でよく扱う問題として、2変数関数の最小値を扱う問題が出ますね。以下に例を示します。 　　　　　　　　　　　　　　$f(x_0,x_1) = x_0^2 + x_1^2$の最小値を求めよ。 　　 　　　　　　　　　　 学部1年生時代、この問題を解くのに$f(x,y)$を$x,y$で偏微分して、さらに二階の偏微分をして、判別式から最小値の吟味を行いますが、実際、この最小値を求める作業にどのような意味をもつのかは習わないと思います。 今回は、この2変数関数の最小値を求める行為が、ニューラルネットワーク構築にどのような場面で必要になるのかとともに解説できたらなと思います。 2.ニューラルネットワークにおける最小値を求める意味とは 　　　　　　　　　　 "損失関数の最小値を求めることです" もう少し解説したいと思います。 そもそも、ニューラルネットワークでの学習とは、重みとバイアスを訓練データに対応することができるようにする行為です。 そして、この学習という行為の中で大事な関数が損失関数というものです。教師データ（正解）に対しての訓練データは誤差を小さくすることであり、この誤差を最小とするのに損失関数が必要となってきます。 3.勾配降下法 上記で記したように、損失関数の最小値を求める方法として勾配降下法があります。では、具体的に勾配降下法について見ていきましょう。 ここで、勾配降下法の式を以下に示します。 x = x - μ\frac{∂f}{∂x}\\ y = y - μ\frac{∂f}{∂y}\\ と表すことができます（μは学習率）。ちなみに、この学習率とは上の式をみればわかるように1回の学習で、どれだけ重みやバイアスを更新するのかの値です。そして、この行為を関数f（損失関数）において、損失関数の最小値を求めていくことが大事になります。 f(x_0,x_1) = x_0^2 + x_1^2 の最小値を実際に求めてみましょう。 import numpy as np import sys def numercial_gradient(f, x): h = 1e-4 grad = np.zero_like(x) for idx in range(x.size): tmp_val = x[idx] x[idx] = tmp_val + h fxh1 = f(x) x[idx] = tmp_val - h fxh2 = f(x) grad[idx] = (fxh1 - fxh2)/(2*h) x[idx] = tmp_val return grad #勾配を定義 def gradient_descent(f, init_x, lr=0.01, step_num=100): ##init_xは初期値、学習率lrを適切な値に設定、step_numは繰り返し数を表す x = init_x for i in range(step_num): grad = numerical_gradient(f, x) x -= lr*grad #勾配降下法の式 return x これで、勾配法の式を表すことはできました。次に、ここで定義した勾配法の式を用いて$f(x_0,x_1) = x_0^2 + x_1^2$の最小値を求めていきます。 def function_2(x): return x[0]**2 + x[1]**2 init_x = np.array([-2.0, 3.0]) gradient_descent(function_2, init_x=init_x, lr=0.1, step_num=100) これを実行すると このような結果になりました。これは、($x_0,x_1$) = ($0,0$)を示していますね。 4.勾配法による更新のプロセスを追う この勾配法で追った過程を視覚的に捉えてみましょう import sys from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt import numpy as np from matplotlib import cm def func2(X, Y): return X ** 2 + Y ** 2 def main(args): x_0 = np.arange(-5, 5, 0.25) x_1 = np.arange(-5, 5, 0.25) X, Y = np.meshgrid(x, y) Z = func2(X, Y) fig = plt.figure() ax = Axes3D(fig) ax.set_xlabel("x") ax.set_ylabel("x") ax.set_zlabel("f(x, y)") contour = ax.contourf(X, Y, Z, levels=20, cmap=cm.coolwarm) fig.colorbar(contour, shrink=0.5, aspect=10) plt.show() if __name__ == "__main__": main(sys.argv) 今回の関数を等高線つきで見ると、以下のようになる。 たしかに$f(x,y)=(0,0)$で最小値をとっているようにみえますが、3次元で描かれたグラフではどこが最小値か見えにくいですね。次では、2次元に落とし込んで考えてみます。 コードは割愛します（長く冗長になるため）が、結果だけ示します。 これを見る限り$f(x_0,x_1)=(0,0)$で最小値をとっていることが分かりますね！ 5.まとめ 今回は、実際に大学で扱われる数学の問題が実際にディープラーニングという場に応用されている例を見ました。今回は基本的な関数で最小値を求めましたが、実際は損失関数を使います（実験レポートを作成する際によく使う二乗和誤差もその一部です。） 大学で習う数学が解くで終わるのではなく、応用できているんだという面白さが共有できていたら幸いです。

「PuLPで解く線形計画法と連立方程式とナップサック問題」と「PuLPで解く最適割当問題と最適輸送問題」の続編として、今度は最大流問題と最小費用流問題を解いてみます。 インストール !pip install pulp Collecting pulp [?25l Downloading https://files.pythonhosted.org/packages/14/c4/0eec14a0123209c261de6ff154ef3be5cad3fd557c084f468356662e0585/PuLP-2.4-py3-none-any.whl (40.6MB) [K |████████████████████████████████| 40.6MB 76kB/s [?25hCollecting amply>=0.1.2 Downloading https://files.pythonhosted.org/packages/f3/c5/dfa09dd2595a2ab2ab4e6fa7bebef9565812722e1980d04b0edce5032066/amply-0.1.4-py3-none-any.whl Requirement already satisfied: pyparsing in /usr/local/lib/python3.7/dist-packages (from amply>=0.1.2->pulp) (2.4.7) Requirement already satisfied: docutils>=0.3 in /usr/local/lib/python3.7/dist-packages (from amply>=0.1.2->pulp) (0.17.1) Installing collected packages: amply, pulp Successfully installed amply-0.1.4 pulp-2.4 必要なライブラリのインポート import pulp import numpy as np import matplotlib.pyplot as plt 問題設定 座標 (0, 0) から (7, 7) の 8 x 8 のグリッド状の頂点（ノード）を考えます。これらの頂点間に、次のように縦・横・斜めの辺（エッジ）が存在するものとします。 grid_size = 8 plt.figure(figsize=(12,12)) checkpoints = [] for x1 in range(grid_size): for y1 in range(grid_size): if x1 == y1: pass checkpoints.append((x1, y1)) plt.scatter([x1], [y1]) edges = {} for p1 in checkpoints: for p2 in checkpoints: if p1[0] == p2[0] and abs(p1[1] - p2[1]) == 1: plt.plot([p1[0], p2[0]], [p1[1], p2[1]]) elif p1[1] == p2[1] and abs(p1[0] - p2[0]) == 1: plt.plot([p1[0], p2[0]], [p1[1], p2[1]]) elif p1[1] - p2[1] == 1 and p1[0] - p2[0] == 1: plt.plot([p1[0], p2[0]], [p1[1], p2[1]]) else: continue if p1 not in edges.keys(): edges[p1] = {} if p2 not in edges.keys(): edges[p2] = {} weight = int(np.random.rand() * 50) if p2 not in edges[p1]: edges[p1][p2] = weight if p1 not in edges[p2]: edges[p2][p1] = weight plt.text((p1[0] + p2[0]) / 2, (p1[1] + p2[1]) / 2, edges[p1][p2]) 上の図で、辺に書かれている数字は、その辺の容量（その辺に沿って流すことができる「品物」の最大量）とします。 最大流問題 最大流問題（Maximum flow problem）とは、供給点（ソース) から需要点（シンク）まで流すことができる「品物」の最大量を求める問題です。 ここでは、原点 (0, 0) を供給点（ソース) とし、点 (7, 7) を需要点（シンク）とします。 まず PuLP を使って、問題設定するためのインスタンスを作成しましょう。 import pulp prob = pulp.LpProblem('max_flow_problem', sense = pulp.LpMaximize) 辺ごとの重み（品物の輸送量）の最小値がゼロ、最大値がその辺の容量となるようにします。 val = {} for k1, v1 in edges.items(): val[k1] = {} for k2, v2 in v1.items(): val[k1][k2] = pulp.LpVariable("({},{})-({},{})".format( k1[0], k1[1], k2[0], k2[1] ), lowBound = 0, upBound = v2) 最大化すべき値は、供給点（ソース) から需要点（シンク）まで流すことができる「品物」の最大量ですが、これは「供給点 (0, 0) から出てくる品物の総和」を最大化することと同義です。 prob += pulp.lpSum([v for v in val[(0, 0)].values()]) 「供給点（ソース) と需要点（シンク）以外の頂点（ノード）において、入ってくる品物の量の総和と出ていく品物の量の総和は同じである」という制約条件を加えます。 ここで、ノード (i, j) は (i - 1, j), (i, j - 1), (i - 1, j - 1) から品物が入ってきて、(i + 1, j), (i, j + 1), (i + 1, j + 1) に品物が出ていくものとし、逆流はしないものとします。 for x, y in checkpoints: if x == grid_size - 1 and y == grid_size - 1: continue elif x == 0 and y == 0: continue sahen = 0 if (x, y - 1) in val[(x, y)]: sahen += val[(x, y)][(x, y - 1)] if (x - 1, y) in val[(x, y)]: sahen += val[(x, y)][(x - 1, y)] if (x - 1, y - 1) in val[(x, y)]: sahen += val[(x, y)][(x - 1, y - 1)] uhen = 0 if (x, y + 1) in val[(x, y)]: uhen += val[(x, y)][(x, y + 1)] if (x + 1, y) in val[(x, y)]: uhen += val[(x, y)][(x + 1, y)] if (x + 1, y + 1) in val[(x, y)]: uhen += val[(x, y)][(x + 1, y + 1)] prob += sahen == uhen 今回の問題設定では、実装を簡単にするため、頂点 k1 から k2 に向かう辺と、 k2 から k1 に向かう辺の流量を同じであるとしておきます。 for k1, v1 in val.items(): for k2, v2 in v1.items(): prob += val[k1][k2] == val[k2][k1] これで目的関数の入力と制約条件の入力が完了したので、ソルバーを実行します。 result = prob.solve() pulp.LpStatus[result] 'Optimal' 結果は次のようになりました。 plt.figure(figsize=(12,12)) plt.title("Optimal value = {}".format(prob.objective.value())) for k1, v1 in val.items(): for k2, v2 in v1.items(): if v2.value() > 0: plt.plot([k1[0], k2[0]], [k1[1], k2[1]], linewidth=int(v2.value()), alpha=0.5) #print(k1, k2, v2.value()) plt.text( (k1[0] + k2[0]) / 2, (k1[1] + k2[1]) / 2, "{} / {}".format(int(v2.value()), edges[k1][k2]) ) for x1 in range(grid_size): for y1 in range(grid_size): plt.scatter([x1], [y1]) 品物の流れがゼロになった辺は省略し、たとえば容量が 33 の辺に品物が 2 だけ流れた時は 2 / 33 のように表示しています。流れの多い辺は太く表示しています。 最小費用流問題 最小費用流問題 (Minimum cost flow problem) は、最大流問題が最大の流量を求めるのに対して、流量を固定したときにそのコストの総量を最小化する問題です。 ここでは、上の最大流問題で求まった最大流量で固定することにしましょう。 total_flow = prob.objective.value() 今回新たに加えるのは、各辺において 1 の量の品物を流すためのコストです。 costs = {} for k1, v1 in edges.items(): for k2, v2 in v1.items(): cost = int(np.random.rand() * 100) if k1 not in costs.keys(): costs[k1] = {} if k2 not in costs.keys(): costs[k2] = {} costs[k1][k2] = cost costs[k2][k1] = cost 次のように、最大流問題のときとほぼ同じコードで解けます。違う点は、目的関数が「辺のコスト」と「流量」の積和であるという点です。 import pulp prob = pulp.LpProblem('minimum_flow_problem', sense = pulp.LpMinimize) sum_cost = 0 for k1, v1 in val.items(): for k2, v2 in v1.items(): sum_cost = costs[k1][k2] * v2 prob += sum_cost prob += pulp.lpSum([v for v in val[(0, 0)].values()]) == total_flow for x, y in checkpoints: if x == grid_size - 1 and y == grid_size - 1: continue elif x == 0 and y == 0: continue sahen = 0 if (x, y - 1) in val[(x, y)]: sahen += val[(x, y)][(x, y - 1)] if (x - 1, y) in val[(x, y)]: sahen += val[(x, y)][(x - 1, y)] if (x - 1, y - 1) in val[(x, y)]: sahen += val[(x, y)][(x - 1, y - 1)] uhen = 0 if (x, y + 1) in val[(x, y)]: uhen += val[(x, y)][(x, y + 1)] if (x + 1, y) in val[(x, y)]: uhen += val[(x, y)][(x + 1, y)] if (x + 1, y + 1) in val[(x, y)]: uhen += val[(x, y)][(x + 1, y + 1)] prob += sahen == uhen for k1, v1 in val.items(): for k2, v2 in v1.items(): prob += val[k1][k2] == val[k2][k1] result = prob.solve() pulp.LpStatus[result] 'Optimal' 結果は次のようになりました。 plt.figure(figsize=(12,12)) plt.title("Optimal value = {}".format(prob.objective.value())) for k1, v1 in val.items(): for k2, v2 in v1.items(): if v2.value() > 0: plt.plot([k1[0], k2[0]], [k1[1], k2[1]], linewidth=int(v2.value()), alpha=0.5) #print(k1, k2, v2.value()) plt.text((k1[0] + k2[0]) / 2, (k1[1] + k2[1]) / 2, "{} / {}".format(int(v2.value()), edges[k1][k2])) for x1 in range(grid_size): for y1 in range(grid_size): plt.scatter([x1], [y1]) 最大流問題で得られた結果と、最小費用流問題で得られた結果は、流量は同じですが微妙に異なりますね。どちらも、解が複数ある場合はその中の１つしか出力してくれないと思います。 最小費用流問題の総流量を別の数字に変更すると当然違う結果が得られます。試してみてください。

概要 PythonでWeb操作の自動化をする際に、ヘッドレス設定をそのまま使えるところがなく困ったのでメモしておく 環境 Selenium 4.0.0.a7 beta4 まで出てたのね・・nuget beta4 python 3.9.5 コード例 automateWeb.py # -*- coding: utf-8 -*- from selenium import webdriver from selenium.webdriver.common.by import By from selenium.webdriver.edge.options import Options options = Options() options.use_chromium = True options.headless = True service_args = ['--verbose'] browser = webdriver.Edge( executable_path='msedgedriver.exe', options=options, # service_args = service_args, # service_log_path=service_log_path, # verbose=True ) browser.implicitly_wait(10) browser.get('{target url}') # 省略 参考 Selenium 4 With Python: All You Need To Know テスト オートメーションに WebDriver (Chromium) を使用する Chromium-固有 のオプションを使用する

今回やること 前にJWT認証の記事を上げたのですがバックエンドの説明が皆無だったのでこれから複数回に分けてあげていきます（単に暇なだけというのもある）。使う言語はGoとPythonで迷ったのですが、Gin（Goの言語のAPIフレームワーク）は日本語の文献がほぼ皆無なので需要がありそうなPythonのDjangoを使います。今回からはサンプルとして簡単な記事を投稿できるシステムを作ります。 APIの設計 APIとはアプリケーションプログラミングインターフェースのことで、特にこの記事ではクライアント（Webブラウザなど）からのHTTPリクエストに応じてデータベースなどとやり取りして特定のデータフォーマット(JSONなど)のレスポンスを返すシステムを指すと思ってください。 今回のAPIエンドポイントの設計はこんな感じ。 path method params detail /api/users GET, POST null false /api/user/:id GET, PATCH, PUT, DELETE id: int true /api/posts GET, POST null false /api/posts/:slug GET, PATCH, PUT, DELETE slug: string true /api/token POST null true /api/token/refresh POST null true /api/token/logout POST null true slugとはurlによくくっついてるやつです(hello-worldみたいな)。5つ目からはJWTにかかわる話なのでここでは割愛します。 今回必要なもの まず今回必須になるものから説明します。 * Python 3 * Django * Django rest framework * pip * virtualenv virtualenvはpythonの仮想環境を作成＋管理するためのパッケージです。このおかげでバージョンの管理などで悩まなくて済みます。作業をする時は必ず仮想環境に入るようにしてください。 pipはpythonのパッケージ管理ソフトウェアです。Anacondaをインストールされている方はすでに入っていると思います。 必須でないものまでまとめるとこんな感じ。 requirements.txt appdirs==1.4.4 asgiref==3.2.10 black==21.6b0 certifi==2020.12.5 cffi==1.14.5 chardet==4.0.0 click==8.0.1 colorama==0.4.4 coreapi==2.3.3 coreschema==0.0.4 cryptography==3.4.7 defusedxml==0.7.1 Django==3.1 django-cors-headers==3.7.0 django-filter==2.4.0 django-templated-mail==1.1.1 django-versatileimagefield==2.0 djangorestframework==3.12.4 djangorestframework-jwt==1.11.0 djangorestframework-simplejwt==4.6.0 idna==2.10 itypes==1.2.0 Jinja2==3.0.0 jwt==1.2.0 MarkupSafe==2.0.0 mypy-extensions==0.4.3 oauthlib==3.1.0 pathspec==0.8.1 Pillow==8.2.0 pycparser==2.20 PyJWT==2.1.0 python-magic-bin==0.4.14 python-slugify==5.0.2 python3-openid==3.2.0 pytz==2021.1 requests==2.25.1 requests-oauthlib==1.3.0 six==1.16.0 social-auth-app-django==4.0.0 social-auth-core==4.1.0 sqlparse==0.4.1 text-unidecode==1.3 toml==0.10.2 uritemplate==3.0.1 urllib3==1.26.4 インストールするときはまとめて一つのファイルに記しておくと同時にインストールできます。コマンドは後程。 インストールと環境構築 まず仮想環境を作成します。virtualenvは以下のコマンドでインストールできます。 pip install virtualenv // pathが通っていない方はこちら python -m pip install virtualenv 次に仮想環境を作成します。 virtualenv myenv これでmyenvという仮想環境ができました（pythonのバージョンも指定できます)。 仮想環境を起動するには、myenvディレクトリに入って次のコマンドを打ってください。 D:/myenv> ../Scripts/activate.bat (myenv) D:/myenv> こんな感じで(name)みたいなのが出てくるのでこれが印です。出るときは同じ場所でdeactivateと打ってください。 ではここから先ほど紹介したものをインストールします。同ディレクトリに先ほどのrequirements.txtファイルを作成しインストールします。 (myenv) D:/myenv> pip install -r requirements.txt たまに信じられないほどインストールに時間がかかりますが仏の心で待ちましょう。 これですべてインストールできたので、いよいよプロジェクトを起動します。 (myenv) D:/myenv> django-admin startproject backend // backendのところはプロジェクト名 フォルダ移動。そしてアプリを登録。 (myenv) D:/myenv> cd backend (myenv) D:/myenv/backend> python manage.py startapp accounts (myenv) D:/myenv/backend> python manage.py startapp posts これでpostsとaccountsというフォルダがbackendディレクトリに作成されます。構成はこんな感じ。 backend ├─.pytest_cache │ └─v │ └─cache ├─.vscode ├─accounts │ ├─migrations │ │ └─__pycache__ │ └─__pycache__ ├─backend │ └─__pycache__ ├─posts │ ├─migrations │ │ └─__pycache__ │ └─__pycache__ └─utils └─__pycache__ これから話すのはVScodeを使っている人向けです。.vscodeフォルダ(上に書いてあるやつ)にlaunch.jsonというファイルを作成し以下を入力すると便利です。 launch.json { "version": "0.2.0", "configurations": [ { "name": "run server", "type": "python", "request": "launch", "program": "${workspaceFolder}\\manage.py", "args": ["runserver"], "console": "externalTerminal", "stopOnEntry": false, }, { "name": "make migration files", "request": "launch", "cwd": "${workspaceFolder}", "type": "python", "program": "${workspaceFolder}\\manage.py", "args": ["makemigrations"], "console": "externalTerminal" }, { "name": "migrate", "type": "python", "request": "launch", "cwd": "${workspaceFolder}", "program": "${workspaceFolder}\\manage.py", "args": ["migrate"], "console": "externalTerminal" } ] } これで毎回コマンド入力しなくてもF5で実行できます。 プロジェクトの設定 backendフォルダ内にsettings.pyというファイルがあるのでそこに先ほど追加したアプリを登録します。 settings.py # Application definition INSTALLED_APPS = [ 'django.contrib.admin', 'django.contrib.auth', 'django.contrib.contenttypes', 'django.contrib.sessions', 'django.contrib.messages', 'django.contrib.staticfiles', 'rest_framework', # 追加 'corsheaders', #追加 'posts', # 追加 'accounts', #追加 ] MIDDLEWARE = [ 'django.middleware.security.SecurityMiddleware', 'django.contrib.sessions.middleware.SessionMiddleware', 'corsheaders.middleware.CorsMiddleware', 'django.middleware.common.CommonMiddleware', 'django.middleware.csrf.CsrfViewMiddleware', 'django.contrib.auth.middleware.AuthenticationMiddleware', 'django.contrib.messages.middleware.MessageMiddleware', 'django.middleware.clickjacking.XFrameOptionsMiddleware', ] これから新しいアプリを登録するときには必ず同じ操作をしてください 言語の設定もしておきます settings.py # Internationalization # https://docs.djangoproject.com/en/3.2/topics/i18n/ LANGUAGE_CODE = 'ja' TIME_ZONE = 'Asia/Tokyo' USE_I18N = True USE_L10N = True USE_TZ = True 初回からいきなりハードですが、カスタムユーザーモデルを作ります。後で変更しようとして泥沼にはまったことがあるので先に作っておきます。デフォルトから大きく変更する必要はありません。ただ作っておいたほうがかなり便利です。認証にはemailとpasswordを使います。 accountsフォルダのmodels.pyファイルに以下を入力します。 accounts/models.py from django.db import models from django.contrib.auth.models import AbstractBaseUser, BaseUserManager class MyAccountManager(BaseUserManager): def create_user(self, email, username, password=None, **extra_kwargs): if not email: raise ValueError('Users must have an email address') if not username: raise ValueError('Users must have an username') user = self.model( email = self.normalize_email(email), username = username, **extra_kwargs ) user.set_password(password) user.save(using = self._db) return user def create_superuser(self, email, username, password): user = self.create_user( email = self.normalize_email(email), username = username, password = password ) user.is_admin = True user.is_staff = True user.is_superuser = True user.save(using=self._db) return user class User(AbstractBaseUser): email = models.EmailField(verbose_name="email", max_length=60, unique=True) username = models.CharField(max_length=30, unique=True) date_joined = models.DateTimeField( verbose_name='date joined', auto_now_add=True) last_login = models.DateTimeField(verbose_name='last login', auto_now=True) is_admin = models.BooleanField(default=False) is_active = models.BooleanField(default=True) is_staff = models.BooleanField(default=False) is_superuser = models.BooleanField(default=False) USERNAME_FIELD = 'email' REQUIRED_FIELDS = ['username'] objects = MyAccountManager() def __str__(self): return self.username def get_absolute_url(self): return '/users/{}'.format(self.id) # 絶対書くこと(必須) def has_perm(self, perm, obj=None): return self.is_admin # これも絶対書くこと(必須) def has_module_perms(self, app_label): return True class Meta: ordering = ['-date_joined'] 何が起こっているのか知りたくなる気持ちはわかりますがこれは水面下でかなり複雑なことが起こっているのでできればコピペしてください。 さらにこれよりもカスタマイズする方法は次回やります。 このモデルをプロジェクト設定に登録します。 settings.py from pathlib import Path import os from corsheaders.defaults import default_headers # Build paths inside the project like this: BASE_DIR / 'subdir'. BASE_DIR = Path(__file__).resolve().parent.parent # Quick-start development settings - unsuitable for production # See https://docs.djangoproject.com/en/3.2/howto/deployment/checklist/ # SECURITY WARNING: don't run with debug turned on in production! DEBUG = True ALLOWED_HOSTS = [] CORS_ALLOWED_ORIGINS = [ #frontend側の設定 ] CORS_ALLOW_HEADERS = default_headers + ('contenttype', 'Authorization') # Application definition INSTALLED_APPS = [ 'django.contrib.admin', 'django.contrib.auth', 'django.contrib.contenttypes', 'django.contrib.sessions', 'django.contrib.messages', 'django.contrib.staticfiles', 'rest_framework', 'corsheaders', 'posts', 'accounts', 'rest_framework_simplejwt.token_blacklist', #JWT認証を使わないなら要らない ] MIDDLEWARE = [ 'django.middleware.security.SecurityMiddleware', 'django.contrib.sessions.middleware.SessionMiddleware', 'corsheaders.middleware.CorsMiddleware', 'django.middleware.common.CommonMiddleware', 'django.middleware.csrf.CsrfViewMiddleware', 'django.contrib.auth.middleware.AuthenticationMiddleware', 'django.contrib.messages.middleware.MessageMiddleware', 'django.middleware.clickjacking.XFrameOptionsMiddleware', ] ROOT_URLCONF = 'backend.urls' TEMPLATES = [ { 'BACKEND': 'django.template.backends.django.DjangoTemplates', 'DIRS': [], 'APP_DIRS': True, 'OPTIONS': { 'context_processors': [ 'django.template.context_processors.debug', 'django.template.context_processors.request', 'django.contrib.auth.context_processors.auth', 'django.contrib.messages.context_processors.messages', ], }, }, ] WSGI_APPLICATION = 'backend.wsgi.application' # Database # https://docs.djangoproject.com/en/3.2/ref/settings/#databases DATABASES = { 'default': { 'ENGINE': 'django.db.backends.sqlite3', 'NAME': BASE_DIR / 'db.sqlite3', } } # Password validation # https://docs.djangoproject.com/en/3.2/ref/settings/#auth-password-validators AUTH_PASSWORD_VALIDATORS = [ { 'NAME': 'django.contrib.auth.password_validation.UserAttributeSimilarityValidator', }, { 'NAME': 'django.contrib.auth.password_validation.MinimumLengthValidator', }, { 'NAME': 'django.contrib.auth.password_validation.CommonPasswordValidator', }, { 'NAME': 'django.contrib.auth.password_validation.NumericPasswordValidator', }, ] # Internationalization # https://docs.djangoproject.com/en/3.2/topics/i18n/ LANGUAGE_CODE = 'ja' TIME_ZONE = 'Asia/Tokyo' USE_I18N = True USE_L10N = True USE_TZ = True # Static files (CSS, JavaScript, Images) # https://docs.djangoproject.com/en/3.2/howto/static-files/ STATIC_URL = '/static/' # Default primary key field type # https://docs.djangoproject.com/en/3.2/ref/settings/#default-auto-field DEFAULT_AUTO_FIELD = 'django.db.models.BigAutoField' AUTH_USER_MODEL = 'accounts.User' #この部分を新たに追加 MEDIA_URL = '/media/' MEDIA_ROOT = os.path.join(BASE_DIR, 'media') REST_FRAMEWORK = { 'DEFAULT_PERMISSION_CLASSES': ('rest_framework.permissions.AllowAny', ), 'DEFAULT_AUTHENTICATION_CLASSES': ('rest_framework_simplejwt.authentication.JWTAuthentication', ), } 次に管理サイトにこのモデルを追加します。同ディレクトリにあるadmin.pyファイルに移動し以下を入力します。 accounts/admin.py from django.contrib import admin from .models import User admin.site.register(User) 最後にマイグレーションしてサーバーを起動してみてください。 python manage.py makemigrations python manage.py migrate python manage.py runserver 管理サイトが開くはずです。 おまけ utilsなんてフォルダないぞ？と思ったあなた、よく見てますね。これは僕が勝手に作ったもので、ユーザーの権限とかの制限を管理するスクリプトを入れてます。あると便利なのでよかったらどうぞ。 utils/permissions.py from rest_framework import permissions class IsAdminOrReadOnly(permissions.BasePermission): def has_permission(self, request, view): if request.method in permissions.SAFE_METHODS: return True return request.user.is_staff class IsOwnerOrReadOnly(permissions.BasePermission): def has_object_permission(self, request, view, obj): if request.method in permissions.SAFE_METHODS: return True return obj.posts.author == request.user class IsSelfOrReadOnly(permissions.BasePermission): def has_object_permission(self, request, view, obj): if request.method in permissions.SAFE_METHODS: return True return obj == request.user utils/authentications.py from rest_framework_simplejwt.authentication import JWTAuthentication from rest_framework_simplejwt.exceptions import InvalidToken class JWTAuthenticationSafe(JWTAuthentication): def authenticate(self, request): try: return super().authenticate(request=request) except InvalidToken: return None 最後に 次回はいろいろやります。 お疲れさまでした。

Webメルカトルの各タイルのXY座標値とZoom値から端の緯度経度を下記で計算（https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames）。 num2deg.py import sys import math def num2deg(xtile, ytile, zoom): n = 2.0 ** zoom lon_deg = xtile / n * 360.0 - 180.0 lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n))) lat_deg = math.degrees(lat_rad) return (lat_deg, lon_deg) x = int(sys.argv[1]) y = int(sys.argv[2]) z = int(sys.argv[3]) ll = num2deg(x, y + 1, z) ur = num2deg(x + 1, y, z) print(ll[1], ll[0], ur[1], ur[0]) Zoom値=0のとき（1分割）。 $ python num2deg.py 0 0 0 -180.0 -85.0511287798066 180.0 85.0511287798066 Zoom値=1のとき（4分割）。 $ python num2deg.py 0 0 1 -180.0 0.0 0.0 85.0511287798066 $ python num2deg.py 1 0 1 0.0 0.0 180.0 85.0511287798066 $ python num2deg.py 0 1 1 -180.0 -85.0511287798066 0.0 0.0 $ python num2deg.py 1 1 1 0.0 -85.0511287798066 180.0 0.0 Zoom=2のとき（16分割）。 $ python num2deg.py 0 0 2 -180.0 66.51326044311186 -90.0 85.0511287798066 $ python num2deg.py 1 0 2 -90.0 66.51326044311186 0.0 85.0511287798066 $ python num2deg.py 2 0 2 0.0 66.51326044311186 90.0 85.0511287798066 $ python num2deg.py 3 0 2 90.0 66.51326044311186 180.0 85.0511287798066 ...

はじめに 前回に引き続き、物体検出に関する方法について忘備録としてまとめました。 Google Colaboratory上でYOLOv5とtorch hubを使用して、物体検出の推論結果の中身を見ていきます。 前回のおさらい 物体検出やYOLOv5の説明は以下の記事をご覧ください。 Pytorch Hub YOLOv5のPytorch Hubについて興味がある方は以下のリンクを参照してください。 YOLOv5の準備 まずはGoogle ColaboratoryでYOLOv5の準備をします。 !git clone https://github.com/ultralytics/yolov5 %cd yolov5/ !pip install -qr requirements.txt 以上で準備は完了です。 torch hubで物体検出してみる まずはモデルをダウンロードします。 import torch model = torch.hub.load('ultralytics/yolov5', 'yolov5x', pretrained=True) 次に物体検出する画像を指定します。 ローカルの画像を使用する場合はパスを指定します。 imgs = ['https://ultralytics.com/images/zidane.jpg'] results = model(imgs) 結果を保存します。 results.print() results.save() 前回と同様に物体検出ができました。 画像は/content/yolov5/runs/detect/exp/ に zidane.jpgとして保存されます。 Pytorch Hubで結果を確認する 推論結果の内容を確認していきます。 results.pandas().xyxy[0] 実行すると以下の出力を得ます。 物体検出の結果が一覧で表示されました。 各物体に対して、左から物体検出位置の座標、confidence（確信度）、番号、物体名となります。 なお、検出位置の座標は左上が原点(0,0)となることに注意してください。 画像と表を見比べてみると内容がわかりやすいかと思います。 モデルの中身を確認する 今回使用したモデルの物体名と番号は以下のように確認ができます。 print(model.names) ['person', 'bicycle', 'car', 'motorcycle', 'airplane', 'bus', 'train', 'truck', 'boat', 'traffic light', 'fire hydrant', 'stop sign', 'parking meter', 'bench', 'bird', 'cat', 'dog', 'horse', 'sheep', 'cow', 'elephant', 'bear', 'zebra', 'giraffe', 'backpack', 'umbrella', 'handbag', 'tie', 'suitcase', 'frisbee', 'skis', 'snowboard', 'sports ball', 'kite', 'baseball bat', 'baseball glove', 'skateboard', 'surfboard', 'tennis racket', 'bottle', 'wine glass', 'cup', 'fork', 'knife', 'spoon', 'bowl', 'banana', 'apple', 'sandwich', 'orange', 'broccoli', 'carrot', 'hot dog', 'pizza', 'donut', 'cake', 'chair', 'couch', 'potted plant', 'bed', 'dining table', 'toilet', 'tv', 'laptop', 'mouse', 'remote', 'keyboard', 'cell phone', 'microwave', 'oven', 'toaster', 'sink', 'refrigerator', 'book', 'clock', 'vase', 'scissors', 'teddy bear', 'hair drier', 'toothbrush'] personは０番目、tieは27番目であることがわかります。 YOLOv5のモデルを使えば、上記リストにある物体の検出が可能となります。 結果をcsvに出力する 先ほどで出力した表をcsvファイルとして出力します。 results.pandas().xyxy[0].to_csv('/content/yolov5/result.csv') 終わりに 最後までご覧いただきましてありがとうございました。

はじめに 前回に引き続き、物体検出に関する方法について忘備録としてまとめました。 Google Colaboratory上でYOLOv5とtorch hubを使用して、物体検出の推論結果の中身を見ていきます。 前回のおさらい 物体検出やYOLOv5の説明は以下の記事をご覧ください。 torch Hub YOLOv5のtorch Hubについて興味がある方は以下のリンクを参照してください。 YOLOv5の準備 まずはGoogle ColaboratoryでYOLOv5の準備をします。 !git clone https://github.com/ultralytics/yolov5 %cd yolov5/ !pip install -qr requirements.txt 以上で準備は完了です。 torch hubで物体検出してみる まずはモデルをダウンロードします。 import torch model = torch.hub.load('ultralytics/yolov5', 'yolov5x', pretrained=True) 次に物体検出する画像を指定します。 ローカルの画像を使用する場合はパスを指定します。 imgs = ['https://ultralytics.com/images/zidane.jpg'] results = model(imgs) 結果を保存します。 results.print() results.save() 前回と同様に物体検出ができました。 画像は/content/yolov5/runs/detect/exp/ に zidane.jpgとして保存されます。 torch Hubで結果を確認する 推論結果の内容を確認していきます。 results.pandas().xyxy[0] 実行すると以下の出力を得ます。 物体検出の結果が一覧で表示されました。 各物体に対して、左から物体検出位置の座標、confidence（確信度）、番号、物体名となります。 なお、検出位置の座標は左上が原点(0,0)となることに注意してください。 画像と表を見比べてみると内容がわかりやすいかと思います。 モデルの中身を確認する 今回使用したモデルの物体名と番号は以下のように確認ができます。 print(model.names) ['person', 'bicycle', 'car', 'motorcycle', 'airplane', 'bus', 'train', 'truck', 'boat', 'traffic light', 'fire hydrant', 'stop sign', 'parking meter', 'bench', 'bird', 'cat', 'dog', 'horse', 'sheep', 'cow', 'elephant', 'bear', 'zebra', 'giraffe', 'backpack', 'umbrella', 'handbag', 'tie', 'suitcase', 'frisbee', 'skis', 'snowboard', 'sports ball', 'kite', 'baseball bat', 'baseball glove', 'skateboard', 'surfboard', 'tennis racket', 'bottle', 'wine glass', 'cup', 'fork', 'knife', 'spoon', 'bowl', 'banana', 'apple', 'sandwich', 'orange', 'broccoli', 'carrot', 'hot dog', 'pizza', 'donut', 'cake', 'chair', 'couch', 'potted plant', 'bed', 'dining table', 'toilet', 'tv', 'laptop', 'mouse', 'remote', 'keyboard', 'cell phone', 'microwave', 'oven', 'toaster', 'sink', 'refrigerator', 'book', 'clock', 'vase', 'scissors', 'teddy bear', 'hair drier', 'toothbrush'] personは０番目、tieは27番目であることがわかります。 YOLOv5のモデルを使えば、上記リストにある物体の検出が可能となります。 結果をcsvに出力する 先ほどで出力した表をcsvファイルとして出力します。 results.pandas().xyxy[0].to_csv('/content/yolov5/result.csv') 終わりに 最後までご覧いただきましてありがとうございました。

目次 APIとWeb API MVCとMTV Model Tamplate View まとめ APIとWeb API API Application Programming Interfaceの頭文字を取った略称。 サービスの提供者がそのサービスを利用するために提供するインタフェースの事です。 開発者がAPIを利用すれば、同じ機能をもったサービスを開発する必要がないため 開発効率が向上し、開発費用を抑えることができる。 Web API API提供者とAPI利用者をHTTP／HTTPSベースでやりとりするAPI APIは提供者と利用者のプログラミング言語が同じであることが多いが、 WebAPIは、異なるプログラミング言語でアプリケーションを連携させることが可能であり、 汎用性が高い。 Web APIではREST,SOAPという代表的な実装方式が存在します。 REST REST: Representational State Transferの略称。 RESTの考え方に沿ってできたWeb APIをRESTful API,又は、REST APIと呼びます。 １、以下のHTTPメソッドでデータ操作を表します。 GET:取得 POST:作成 PUT,Patch:更新 Delete:削除 ２、ステートレス 前回のAPIの結果と2回目の結果で同じ値を返す。 ３、URIで操作対象のリソースを判別可能にする ユーザーというリソースを取得するため、usersというリソース名を付与したとする。 ユーザーIDで下記のようなURIフォーマットでURIを生成し、操作可能とする URI: http://example.com/api/users/userId/ 生成されたURI：　http://example.com/api/users/1/ ４、レスポンスとしてXML形式かJSON形式で操作結果を戻す データ操作できた場合、XMLもしくはJSONでデータ操作結果を記述し、HTTPのレスポンスボディーに含め、データを返す。 SOAP SOAPはSimpleObject Access Protocolの略称。 リクエスト、レスポンスをXML形式でやり取りする。 MVCとMTV 以下の表ような違いがあります。 MTV 説明 MVC Model データ、ロジック Model Template データの見せ方 View View どのようなデータを見せるか View URLディスパッチャ URLの振り分け、リクエスト制御 Controller Model アプリケーションに必要なデータ、そのデータに付随する処理を記述する場所 Djangoでアプリケーションを作成する際には保存するデータとそのデータに関連する処理を全てModelにまとめる Template Viewから渡ってきたデータを表示するhtmlファイルのこと View Django REST frameworkにおけるViewにおける役割はユーザーからのリクエストを元に、どのAPIを提供するか決定する処理をしている Viewには以下のようなクラスを継承することで簡単に開発することも可能。 FunctionView: Responseを返却する関数型View ClassBasedView: Httpメソッドごとにインスタンスメソッドを定義できる GenericView: クエリセットに紐づくため簡略化した書き方ができる ViewSet 取得、一覧、登録、更新、削除をまとめて管理するビュー 簡潔に素早くAPIロジックを作成したい場合 特別なカスタムロジックを実装しない場合 APIView 複雑なAPIロジックを作成したい場合 外部のAPIを呼び出す場合 localのファイルを扱う場合 最後に これからDRFを学習する方へ向けての投稿でしたが、次回以降、コード例を載せてより詳しく説明しようと思います。 私がDRFを学習する際、RESTって何？、そもそもAPIとWebAPIの違いがわからないと言った事がありましたので、そういった方のお役に立てれば幸いです。

はじめに こんにちは。Yuki | Kagglerです！ 先日、Shopeeコンペの順位が確定して銀メダルをいただき、晴れてCompetition Expertになることができました。区切りがいいのでここまで取り組んできたことをまとめてみました。 プログラミング&機械学習を始めて一年、ようやくKaggle Expertになることができました！！行列も正規分布も知らず、ターミナルなんて触ったこともない状態からのスタートでしたが、ようやくここまで来ました。ここまで来れたのは偏にこれまで関わってきた皆様のお陰です。これからも頑張ります！！ pic.twitter.com/kMkaFhqhU9— ユウキ | Kaggler (@Yuki_Kaggler) May 12, 2021 この記事の対象者 Kaggleをやってみたいが、何から始めればいいかわからない方 機械学習の勉強をしてみたいが、何から始めればいいかわからない方 機械学習の勉強をしているが、なかなか本屋に行けず、どの本を買えばいいかわからない方 Kaggleを頑張っているが、中々上位にいけなかったり、最後まで走りきれずに悩んでいる方 上記にきれいに当てはまらなくても、何か参考になることがあるかもしれません！！ 目次 1.自己紹介 2.注意事項 3.実際に書いていく 4.まとめ 5.今後の目標 1. 自己紹介 TwitterでYuki | Kagglerという名で活動している、東京大学理科一類2年生です。学部・学科は工学部物理工学科に進学予定で、量子機械学習を研究したいと考えています。元々は機械学習を学びたくて、色々なデータに触れて機械学習に対する興味を絞るためにKaggleを始めました。(と言いつつ最近は専ら画像コンペに出てます。) 2. 注意事項 この記事はあくまで僕がExpertになるまでにしたことを書いた記事であり、こうすればExpertになれる、といった類のものではありません。また、数学に関しては大学の授業に合わせて教科書で学んだため、機械学習を勉強するために学んだとは言い難いと判断して詳細は省きます。また、アフィリエイトは一切使っていません。 3. 実際に書いていく 基本的には1ヶ月単位で学んだこと、その時感じたことを書いていきます。 勉強開始前(~2020年3月) 受験が終わったのでとりあえず線型代数を勉強し始めました。この時点では微積・統計の知識は高校生レベルで、線型代数は何も知りませんでした。ちなみに、僕の世代は行列は高校生では習わず、しかも統計も2Bの範囲は無しで受験できたのでせいぜい相関係数は知ってるよレベルです。正規分布なんて聞いたこともありませんでした。もちろんプログラミングはほとんどやったことがなく、ターミナルなんて触ったこともありません。本当に0からのスタートです。 4月 「スッキリわかる Python入門」 この本でPythonを勉強し始めました。これ以上はないぐらい親切な内容で、1冊目として選んでよかったと感じています。 Amazonのリンクはこちら Udemyの「実践Pythonデータサイエンス」講座 ここでnumpy, pandas, matplotlib, seabornなどを一通り学びましたが、はっきり言ってあまりわかってはいませんでした。udemyの講座に関してはもっといい講座に出会えたので後に書きます。最後の機械学習のパートを楽しみにしていたものの、そのパートが全くわからず途方に暮れて手を出したのが次の本です。 リンクはこちら 「Pythonではじめる機械学習」 あっという間に挫折しました。教師あり学習、教師なし学習、特徴量エンジニアリング、モデルの評価、テキストデータの処理などが網羅的に解説されている一冊です。機械学習そのものに関する内容ももちろん豊富ですが、scikit-learnの使い方の解説という趣も割と強かったと感じています。前述の講座であまりnumpyやmatplotlibがわからなかった上に、機械学習そのものも何も知らないという状況だったので、言ってることもコードもわからない箇所が多くて辛かったです。また、日本語が若干不自然なところが割と多く、不自然な日本語を頭の中で補正ができなかったことも、今振り返れば致命的なポイントでした。 Amazonのリンクはこちら 5月 courseraの「Machine Learning」 腰を据えて機械学習を勉強しなければいけないと感じて受講しました。3ヶ月のコースでしたが、覚悟を決めて約10日で全て完了させました。機械学習について、初学者にもわかりやすく網羅的に解説されています。「Pythonではじめる機械学習」よりもっと詳しいので、本当の初心者でも最後までやっていけます。コーディング課題があるので知識が定着しやすいのもいいところです。courseraには有料の講座も多いですが、この講座は無料なのも嬉しいポイントです。全部英語でしたが、受験で英語はかなりやり込んだのもあり、だんだん慣れてきて後半は課題も全て翻訳なしで乗り切ることができました。英語が苦手でも、動画には日本語字幕をつけられますし、課題もDeeplに突っ込めばなんとかなると思うので、思い切って受講をしてしまうのがオススメです。この講座は本当に有名で、ネットにも解説記事が大量にあるので、困ったらそれらを活用するのも手ですね。使われている言語がPythonではなくoctave(またはMatlab)なのが少しネックではあります。ここが自分の中で大きな転換点になり、その後の勉強がうまくいったのでこれは本当にオススメです。 リンクはこちら 「ゼロから作る Deep Learning」(以下ゼロつくと略記) courseraで学ぶ中で深層学習が面白そうと感じて読みました。パーセプトロンから始まり、多層パーセプトロン、誤差逆伝播法、SGDやAdamといった最適化手法、CNNと豊富な内容で構成されています。この本は著者が日本人であることもあり、解説が非常にわかりやすくて日本語にも詰まることなく読めました。コードに関する説明も非常に詳しいので、Pythonに不安がある人でも十分読み切れると感じています。この本で深層学習(特にCNN)に感動して、この分野に人生をかけていいのではないかと思えたからこそ今があるので、僕の中で最高の一冊です。迷っている方はぜひ手に取ることを勧めます。 Amazonのリンクはこちら 「Pythonで始める機械学習」 前回挫折したリベンジをしました。読み切れた理由としては、機械学習に強くなったお陰で違和感のある日本語に悩まされなくなった、ゼロつくを読むうちにコーディング力が上がっていた、などが挙げられます。 中途半端にしか読んでいないので大きくは取り上げませんが、似た本に「Python 機械学習プログラミング」という本があります。この本の方が若干難易度は上がって数学力も求められますが、courseraの説明との親和性が高くてより詳しいので、こちらにチャレンジしてもいいかもしれないです。前半はDeepじゃない機械学習について、後半は深層学習がTensor Flowを用いて解説されています。僕は前半を読みながらコーディングしましたが、非常に勉強になりました。 Amazonのリンクはこちら 6月 「ゼロから作る Deep Learning 2」 ゼロから作るDeep Learningの続編で、自然言語処理編です。前作の復習から始まり、word2vec、シンプルなRNN、LSTM、LSTMを用いたSeq2Seq、Attentionなどがフルスクラッチで作り上げられています。こちらの方が量が多くて難易度も上がっており、前作ほどスラスラとはいきませんでした。しかし説明は非常にわかりやすく力がつくので、ゼロつく1を読み終えた方、深層学習での自然言語処理に興味がある方は是非手に取ってみてください。(全然関係ありませんが、帯の言葉が大好きです。痺れますね〜) Amazonのリンクはこちら 「ゼロから作る Deep Learning 3」 これも前作からの続編で、フレームワークを作っていくというものです。最後には、ゼロから作ったフレームワークでCNNやLSTMを実装して動かすというワクワクするシナリオになっています。1, 2と比べると人気は劣るような印象ですが、個人的には1, 2と同じぐらい大好きです。この本のお陰で何の違和感もなく、むしろ感動してPyTorchを使い始めることができたという点で、今振り返るとこの本も非常に重要でした。さらに分厚いですが、前作や前々作の復習もすることができ、さらにPythonそのもの(特にクラスについて)に対して強くなれました。 Amazonのリンクはこちら G検定の教科書、問題集 G検定というものがあるということを知り、どうせなら受けてみようということで勉強しました。今振り返って必要だったかと言われると必要はなかったと思います。自分が勉強していなかった各分野を浅く知れました。 教科書のAmazonのリンクはこちら 問題集のAmazonのリンクはこちら 7月 G検定合格 ここまでゼロつくなど様々な本を読んでいたので機械学習に関しては余裕でしたが、法律に関してははっきり言って無理ゲーという他なかったので、教科書にある分だけ抑えて、本番でわからないところは40回以上ググってやっつけました。体感は8割5部ぐらいでした。積極的に受ける必要はないかと思いますが、勉強したついでに受けてみるのもそれはそれでいいと思います。G検定とゼロつくなら間違いなくゼロつくです。少なくともKaggleに繋がるのは後者です。個人の所感なのであくまでご参考までに。 「直感Deep Learning」(Keras) 途中で挫折してしまいました。僕はフレームワークを学ぶ前にゼロつく3を読んでいたためにフレームワークはこういうものだというイメージが頭の中に既にあったので、それにそぐわないKerasを受け入れることがどうしてもできませんでした。全然直感的じゃなかったです...(人によるとは思いますが。) Amazonのリンクはこちら あと、本筋からはそれるのですが、GitやDockerもこのタイミングで一度勉強しました。ゼロつく2、3あたりは1人でGitを使う練習も兼ねていました。 8月 「実践Data Scienceシリーズ PythonではじめるKaggleスタートブック」 6月あたりに買ったものの、submitの仕方がわからないという理由で放置していたので引っ張り出してきました。Kaggleは全て英語で書かれているため心理的障壁が高くなりがちですが、この本を読めば取り組み方、コンペの種類など様々なことを知れるのでとてもよかったです。機械学習についての説明も優しく書かれていたので、この本から機械学習を始めるのもありなのかなと思いました。 Amazonのリンクはこちら Kaggle Titanicコンペ ・ House Prise コンペ Kaggleスタートブックを読んで、自分なりにEDAをしたりアンサンブルをしたりしてみました。ほとんどスコアは伸びなかったですが、これまでと違ってちゃんとデータにさわれてる感覚があってとても嬉しかったです。個人的には、Kaggleに取り組み始める前にしっかり勉強する、というよりはできると思ったら即参加してみる方がいいと思います。目標はそれぞれ違うとは思いますが、Kaggleで強くなりたいのであれば、Kaggleを肌で感じで目標を立てて進めた方が圧倒的に早いです(自戒です)。 「米国データサイエンティストがやさしく教えるデータサイエンスのためのPython講座」 udemyの講座に関してはもっといい講座に出会えたので後に書きます。 それがこの講座です。今まで何となく使っていたnumpyやpandas、matplotlibなどにしっかり向き合うことができたので非常によかったです。この講座は発売同時に買ったので4月時点では無かったですが、あの時にこの講座に出会えていたらもっとスラスラ他の本を進められていたんだろうなと考えてしまいます。この講座の作者の「かめさん」という方はこの講座の他にもDockerやGitの講座も出しており、全部わかりやすいです。また、最近、numpyなどのライブラリではなく純粋なPythonの講座を出したようなので、今から買うのであれば目的に応じて選ぶといいと思います。udemyのクーポンが取得できる、かめさんのブログ記事を貼っておきます。こちらからどうぞ。Pythonを学ぶならこの方の講座が一番だと思っています。 講座のリンクはこちら 「PyTorch for Deep Learning with Python Bootcamp」 udemyのPytorchの講座です。書き始めてすぐに、これが今まで求めてたものだと感動した覚えがあります。全て英語なので若干敷居が高いですが、発音がとても綺麗で聞き取りやすいですし、Pytorchがしっかり解説されている動画講座は意外と少ないのでとてもオススメです。 Udemyのリンクはこちら OSICコンペ参加、挫折 Pytorchを身につけた勢いで初めてメダルありのコンペに参加しました。しかし、実際のコンペのコードにとても苦戦し、おまけにセグメンテーションも必要かも知れないということに気づき、あっという間に撤退してしまいました。この後はこのギャップを何とか埋めるべく色々な本を漁ることになりますが、Kaggleで強くなりたいという目標からするとこの判断は非常にまずかったです。Kaggleで強くなるためにはKaggleに向き合うしかないです。1ヶ月ほどコーディングから離れますが、再び戻ってきて同じように絶望したので、ここからの1ヶ月はKaggleに関しての進歩はほぼゼロでした。 E資格の勉強 G検定を受けたのだからE資格も受けてみようと思い、AVILENさんの講座を受講しました。G検定と比べればだいぶ理論が詳しく説明されていて、何となく誤魔化してきた部分をしっかり学習することができたと感じています。僕は大学1年の途中だったので、内容としては特異値分解あたりの話、そして情報理論は改めて勉強する必要がありました。また、強化学習に関してはここまで全く勉強してこなかったのでとても苦しみました。正直今でもあんまりわかってないです。そのほかの内容に関してはほぼこれまでの教材で網羅され尽くされていると感じています。色々な手法を知ることができたので、Kaggleで、この単語知らないなぁ...となることがだいぶ減って楽になりました。受験は2月なので、それまで何度も復習していくことになります。 「これならわかる深層学習入門」 E資格のために、深層学習の理論面の勉強で用いた本です。E資格を受験する際にはよく「深層学習」というめちゃくちゃ厚い本を勧められますが、僕はあれは明らかにオーバーワークだと思ったのでやめて、そこで代わりに手に取りました。深層学習だけでなく、機械学習全般についてわかりやすく書かれており、coursera、はじパタで手に入れた理解をアップデートすることができました。個人的にはこの本は強化学習の章がとてもわかりやすく、E資格の勉強に際して何度も何度も読んで式変形を追いました。深層学習の理論を学ぶならイチオシの本です。 Amazonのリンクはこちら 9月 統計、機械学習に関する諸々の本 OSICで撃沈してから、若干現実逃避気味になり色々な本(完全独習ベイズ統計入門、統計学入門、データ解析のための統計モデリング入門、多変量解析、これなら分かる最適化数学、はじパタ、 ...etc)を読みました。統計にもともと興味があったのでそれはそれでよかったのですが、Kaggleで強くなりためにはだいぶ遠回りをしたと思ってます。 「つくりながら学ぶ! Pytorchによる発展Deep Learning」 物体検出やセグメンテーション、Transformerなどこれまで紹介した書籍には中々掲載されていない発展的な内容が詳しく解説されています。各章で独立しているため、必要に応じて章ごとに読むことができるのでとてもオススメです。ソースコードも豊富なので、基本的な書き方を一通り理解した後に実践的な書き方を学ぶのにちょうどよかったと感じています。多くの人がこの本でPytorchを学んでいる印象です。 Amazonのリンクはこちら ちなみに、現在はPytorchの入門書として有名であった洋書が翻訳されたPytorch実践入門が出版されており、Pytorchを深く学ぶならばこの本がいいのかもしれないです。(僕はもうKaggleに何度も参戦してそこで学んだので特に必要性をあまり感じておらず、買ったものの積読状態です。) Amazonのリンクはこちら MoAコンぺ参加 初心者でも取り組みやすい雰囲気のテーブルデータコンペ・かつニューラルネットが効きそうだったので参加しました。これが初めての本格参戦になります。上でも書きましたが、OSICで苦戦したことはそっくりそのままこちらでも苦戦しました。ここで、Kaggleと向き合っておくべきだったと後悔しました。今度はやり通せたのはタスクの難易度がこちらの方が少し低いのとOSICコンペから逃げた後悔のおかげです。本を読んで力がついたからではないと感じています。 10月 「Kaggleで勝つデータ分析の技術」(通称：Kaggle本) MoAコンペに参加していた時に読みました。コンペに対する取り組み方、様々な手法、バリデーションの仕方、などKaggleで重要な部分をたくさん学ぶことができました。現在はテーブルデータのコンペが減ったためにこの本に書いていることがそのまま全部通用することは少なくなりましたが、それでも上述のようなコンペにおける重要な部分を学ぶことは非常に意義深いことだと思うので、持っていて損はしないと思います。色々なことが網羅的に書かれているため、辞書的な使い方をしています。長く使える一冊なので、Kaggleに取り組むならばおすすめです。 Amazonのリンクはこちら 改めて書くことはしませんがMoAコンペを頑張り続けていました。 11月 MoAコンペ終了 Top13%で敗北しました。その時の振り返り記事があるのでもしよければご覧になってください。初めてコンペをやり切った感想、反省を詳しく書いています。こちらから飛べます。 Gitの勉強 詳しくは振り返り記事を参照していただきたいですが、MoAコンペの終盤でNotebook管理が崩壊して失速してしまったため、それを反省して管理体制を見直そうと考え、改めてGitを勉強しました。その時に使ったのは「Git: もう怖くないGit! チーム開発で必要なGitを完全マスター」というudemyの講座です。この直後に、Kaggle Notebookで、特定のVersionにボタンひとつで戻れる機能がついたので結局Gitは使わなかったです。 Udemyのリンクはこちら 12月~2月 いきなりまとめてすみません。Cassavaコンペにほぼ全振りしたので書くことが少ないです。 Cassavaコンペ もともと画像処理に一番興味があったので、MoAの雪辱を晴らすために参加しました。キャッサバの葉の画像を、病気や健康などのラベルに分類するタスクでした。ゼロつくで色々学んだとはいえやはり実際にコンペに取り組むとなるとハードルは高く、最初は非本質的でしょーもない実験も繰り返していました。この頃からKaggle日記(参考記事, fkubotaさん著)を使い始めました。これが自分の中でとてもヒットして、以降のコンペでもスムーズに実験を進めることができ、かつ結果をいつでも振り返れるので自分の経験値にしていくことができるようになり、成長が早まりました。結果はTop11%で終わってしまいました。終わった後はショックが大きすぎて本当にもうやめようかと思いましたが、腹いせに飲んだタピオカミルクティーがおいしかったので何とか元気を取り戻しました。振り返り記事はこちらで、Kaggle日記はこちらです。もしよろしければご覧ください。 E資格合格 2月下旬の試験で無事合格することができました。かなり高額ですが、受けてよかったと感じています。Cassavaコンペにも大きく生きました。ただ、満点を取れなかったことが悔やまれます。 E資格合格してました。なんで数学が満点じゃないんでしょうね...全体でも満点を取りたかったですが、力及びませんでした。これからも謙虚に頑張ります！#E資格 pic.twitter.com/eIPNDDKB8o— ユウキ | Kaggler (@Yuki_Kaggler) March 11, 2021 3月 RANZCRコンペ Cassavaコンペのリベンジを懸けて参加しました。CT画像を元に、カテーテルの種類などを分類するコンペでした。ここまで1人で戦ってきて全てだめだったので、思い切ってひらさんというかたにお声かけしてチームを組ませていただきました。チームを組むと自分で全部やる必要がなくなるため自分のやるべきことを深掘りできるようになったり、相手の戦い方を見たりして、1人では得られない学びをたくさん得ることができました。そして何より楽しかったです。画像処理の勉強をした甲斐があってTop6%で銅メダルを獲得することができました。 Kaggle日記はこちらです。もしよろしければご覧ください。 楽しさを伝えたいので、一番楽しかったときのslackでのやりとりを貼っておきます。 4月・5月 Shopeeコンペに参加 こちらもチームで、今回は4人で参加しました。shopeeという東南アジアのメルカリのようなサイトの商品の画像と文章を元にマッチングをするタスクでした。Top3%銀メダルを獲得することができ、Expertに昇格することができました。チームメイトの方々には本当に感謝です。 自然言語処理についてはBERTを動かすので精一杯だったので今後の課題です。Kaggle日記はこちらで、KaggleのDiscussionに掲載している解法はこちらです。もしよろしければご覧ください。 4. まとめ 僕が偉そうに言えることは何もありませんが、とにかく伝えたいのはKaggleで強くなりたいならKaggleをするのが最も近道であるということです。本で勉強した後にKaggleのコンペに実際に参加するとギャップに苦しむことになる可能性が高いですが、Kaggleには公開NotebookやDiscussionなどヒントがたくさん落ちているので、腰を据えて(英語が苦手であればDeeplを使って)向き合っていくことが大切だと思います。 そして、今回僕があえてExpertになった段階で勉強したことを全て公開したのは、Master以上の方々が書いた素晴らしい記事が既にたくさんあるからです。特に初心者の場合はMaster以上の方々との距離がだいぶ大きいので、距離が近い僕の記事にも価値があると判断しましたが、Master以上の方々の記事も是非ご覧になってください。僕自身も何度も見て参考にさせていただきました。 5. 今後の目標 言わずもがな、Kaggle Masterです。今年中に確実に昇格します。また、ここまでのメダルは全てチームでとっているため、どこかのタイミングでソロでメダルを取ります。 ここまでご覧いただきありがとうございました！少しでもお役に立てたのであれば嬉しい限りです。Happy Kaggling!!

はじめに この記事を見て、タイトルの動作ををしようとした時にハマったので、メモしておきます。 解決はしましたが、根本的なメカニズムを理解できていないので、詳しい方がいたらご教授ください。 やりたいこと MultiIndexのDataFramedfから、Index2 == 'abc'のものだけ取得する。 Index1 Index2 Content A abc A_abc B de B_de C abc C_abc D de D_de E abc E_abc F de F_de 解決方法 df.loc[(slice(None), 'abc' ), :] または df.loc[pd.IndexSlice[:, 'abc'], :] 詳細はこちらの記事を参考にしてください。 できなかった方法 Pandasのloc[]では後ろの:を省略できるため、通常この2つは同じ動作をします。 df.loc['C'] df.loc['C', :] ところが、今回:を省略したところエラーが発生。 df.loc[(slice(None), 'abc' )] # KeyError: 'abc' df.loc[pd.IndexSlice[:, 'abc']] # KeyError: 'abc' ちなみに、これは問題なく実行できました。 df.loc[(slice(None), 'abc' ), ] df.loc[pd.IndexSlice[:, 'abc'], ] エラーの理由がイマイチよく分かっていないので、知っている方がいれば教えていただけると助かります。 参考 pandasのMultiIndexから任意の行・列を選択、抽出(note.nkmk.me) https://note.nkmk.me/python-pandas-multiindex-indexing/

今回やる事 DRF、Vue.js、dockerでの環境構築をする上で、各ファイルの意味だったりコマンドだったりをよく忘れるので、とりあえずDRFとVueが連携出来るまでを、ほぼ自分用に残しておきます。 例として、ユーザー一覧を取得するだけの処理を実装します。 環境 docker $ docker version Version: 20.10.6 docker_compose $ docker-compose -v docker-compose version 1.29.1, build c34c88b2 python $ python --version Python 3.8.10 構築編 DBは開発段階ではsqliteを使い、DRF用にコンテナ1つ、Vue用にコンテナ1つを用意します。 サーバーはそれぞれの開発サーバーを使用するとします。 まずは、以下の様なディレクトリ構成でファイルを作成します。 アプリ名は「k2」です。 k2/ ├ api │ ├ Dockerfile │ └ requirements.txt ├ web │ └ Dockerfile └docker-compose.yml ※一部省略。 docker-composeを定義 まずは、docker-composeでDRF用のコンテナとVue用のコンテナを定義します。 docker-compose.yml version: '3' services: api: # Dockerfileの場所 build: ./api # どのimageを使うか。Dockerfileでimageは作成する。 image: k2-drf-image # コンテナ名を指定出来る。 container_name: k2-api # 公開用のポート。ホスト側:コンテナ側を指定 # コンテナ側のみ指定も可能だが、その時ホスト側はランダムになる。 ports: - '8000:8000' # マウントするパス volumes: - ./api/:/usr/src/k2/api/ # dockerコンテナを落とさず起動し続けるようにする tty: true web: build: ./web image: k2-vue-image container_name: k2-web ports: - '8080:8080' volumes: - ./web/:/usr/src/k2/web/ tty: true 次に、apiとwebのDockerfileを定義します。 apiのDockerfile # イメージを選択 FROM python:3.9.5-alpine # バイナリレイヤ下での標準出力とエラー出力を抑制 ENV PYTHONUNBUFFERED 1 # 開始ディレクトリ位置 WORKDIR /usr/src/k2/api # 必要なものをインストール # 今回はalpineなのでapkを使用する。 RUN apk update RUN apk add --no-cache --vertual=wow-ini-module \ net-tools \ sudo \ bzip2 \ curl \ gcc \ git \ python3-dev \ vim \ bash # ADDかCOPYを使うなら.dockerignoreを書いておく。 ADD . . RUN pip install --upgrade pip RUN pip install -r requirements.txt とりあえず最小限をrequirements.txtに記述する。 requirements.txt Django==3.2.4 django-cors-headers==3.4.0 django-environ==0.4.5 django-filter==2.3.0 djangorestframework==3.12.4 djangorestframework-jwt==1.11.0 gunicorn==20.0.4 requests==2.24.0 webのDockerfile FROM node:10.17.0 WORKDIR /usr/src/k2/web/ RUN apt-get update RUN apt-get install -y --no-install-recommends \ net-tools \ sudo \ bzip2 \ curl \ gcc \ git \ vim RUN apt-get clean RUN npm install -g yarn \ && yarn global add @vue/cli \ イメージ構築、コンテナ作成～立ち上げ $ docker-compose up -d --build コンテナに接続し、環境を整える ターミナルを2つ立ち上げてそれぞれで接続し、それぞれの最低限の環境を整える。 DRF側 $ docker-compose exec api bash $ django-admin startproject api $ django-admin startapp k2 $ python manage.py makemigrations $ python manage.py migrate $ python manage.py createsuperuser 上記で作ったファイル群を以下のフォルダ構成に直します。 k2/ └ api ├ api/ ├ k2/ ├ db.sqlite3 ├ Dockerfile ├ manage.py └ requirements.txt とりあえずこの段階で、python manage.py runserver 0.0.0.0:8000で開発サーバーは立ち上がる。 Vue側 $ docker-compose exec web bash vueプロジェクトの作成 ①現在のディレクトリにフォルダ群を展開する。 $ vue create . ? Generate project in current directory? (Y/n) y ②マニュアルで導入するものを選択する。 ? Please pick a preset: Default ([Vue 2] babel, eslint) Default (Vue 3) ([Vue 3] babel, eslint) ❯ Manually select features ③導入するものにスペースで☑を入れる。 ? Check the features needed for your project: ◉ Choose Vue version ◉ Babel ◯ TypeScript ◯ Progressive Web App (PWA) Support ◉ Router ◉ Vuex ◉ CSS Pre-processors ◉ Linter / Formatter ◉ Unit Testing ❯◉ E2E Testing ・Babel 　JSのコードを新しい書き方から古い書き方へと変換するツール。 　ブラウザによって古いバージョンしか対応していないため、このようなツールが必要になる。 ・PWA 　ネイティブアプリのような使い勝手を実現したWEBアプリの事。 　インストール不要で、ホーム画面やアイコン追加やプッシュ通知などが可能。 　そのほかに読み込み速度や表示の高速化、オフラインで閲覧化などのメリットもある。 ・CSS Pre-processor 　Sass, Less, Stylusなどのプリプロセッサ独自の構文でCSSを作成するプログラム。 　純粋なCSSにはない、ミックスイン、セレクターの入れ子、セレクター継承などが可能。 ・Linter 　ESLint等のコードフォーマッターなどを使用するか選択する。 ・Unit Testing（単体テスト） 　Mocha+Chainや、Jestなどを使用するか選択する。 ・E2E Testing（End to Endテスト） 　システム全体を通してテストを行う事。 　それらのライブラリなどを使うか選択する。 ④vueのバージョンを選択する。 3.x系だとvuetify Alphaしか対応してないみたいなので、今回は普通の2.xを使います。 ? Choose a version of Vue.js that you want to start the project with ❯ 2.x 3.x ⑤routerのhistoryモードを使うか選択。 ? Use history mode for router? (Requires proper server setup for index fallback in producti on) (Y/n) y ・historyモード 　ページのリロード無しにURL遷移を実現するhistory.pushState APIを利用したルーターのモード。シングルページのクライアントアプリなので、適切なサーバーの設定をしないと、404エラーが起きる。そのため、どのアセットにもURLがマッチしなかった場合に、index.htmlを返す設定が必要。 ・hashモード（vue-routerのデフォルト） 　完全なURLをhashを使ってシミュレートし、URLが変更された時にページのリロードが起きない。 ⑥ どのCSS pre-processorを使うか選択。 ? Pick a CSS pre-processor (PostCSS, Autoprefixer and CSS Modules are supported by default) : (Use arrow keys) ❯ Sass/SCSS (with dart-sass) Sass/SCSS (with node-sass) Less Stylus ⑦ linterを選択する。 ? Pick a linter / formatter config: (Use arrow keys) ❯ ESLint with error prevention only ESLint + Airbnb config ESLint + Standard config ESLint + Prettier ⑧ lintの追加の機能を設定 ? Pick additional lint features: (Press <space> to select, <a> to toggle all, <i> to invert selection) ❯◉ Lint on save ◯ Lint and fix on commit ⑨ unitテストに何を使用するか。 ? Pick a unit testing solution: (Use arrow keys) ❯ Mocha + Chai Jest ・Mochaは、BDD/TDDをするための枠組みを提供、chaiはテストを実施するための便利なメソッドを提供してくれる。TDD（テスト駆動開発）/BDD（振る舞い駆動開発） ・JestはFacebookが開発したテストランナ。 ⑩ E2Eテストに何を使用するか。 ? Pick an E2E testing solution: ❯ Cypress (Chrome only) Nightwatch (WebDriver-based) WebdriverIO (WebDriver/DevTools based) ・Cypress 　E2Eテストフレームワーク ・Nightwatch 　Seleniumを基盤として作られており、Node.jsでブラウザを操作してE2Eテストを行うツール ・WebdriverIO 　ウェブブラウザによるテストを自動的に実行するライブラリ ⑪ Babel, ESLintなどをどのファイルで管理するか。 ? Where do you prefer placing config for Babel, ESLint, etc.? (Use arrow keys) ❯ In dedicated config files In package.json ⑫ プリセットをセーブするか ? Save this as a preset for future projects? (y/N) n ⑬ どちらのパッケージマネージャーを使うか。 ? Pick the package manager to use when installing dependencies: (Use arrow keys) ❯ Use Yarn Use NPM これでvueのプロジェクトが出来上がり、yarn serveで開発サーバーが立ち上がる。 あと、必要に応じて下記のライブラリ等も導入しておく。 下の連携テストでは、 ・vuetify ・sass ・axios ・boxicons ・vue-loading-template を使っている。 # update yarn install # vuetify yarn add vuetify (vue add vuetify) # sass yarn add --dev sass # vuesax npm install vuesax # material-icons（アイコン） npm install material-icons --save # axios npm install --save axios vue-axios # boxicons（アイコン） npm i boxicons # lodash npm i --save lodash # vue-persistedstate npm i vuex-persistedstate # vue-session npm i vue-session # vue-loading-template npm install vue-loading-template --save 連携テスト とりあえず、ユーザー一覧をDRFから受け取ってそれを表示させるところまでをやります。 DRF側 やる事 初期設定 ルーティング View記述 DRF: 初期設定 k2/api/api/settings.py # 追加 import logging logging.basicConfig( level=logging.DEBUG, format='''%(levelname)s %(asctime)s %(pathname)s:%(funcName)s 行数:%(lineno)s:%(lineno)s %(message)s''' # 編集 ALLOWED_HOSTS = ['*'] # 追加 INSTALLED_APPS = [ ... 'rest_framework', # Json Web Token関連のライブラリ 'rest_framework_jwt', 'django_filters', 'k2.apps.K2Config', ] # 追加 REST_FRAMEWORK = { # デフォルトのパーミッションクラス 'DEFAULT_PERMISSION_CLASSES': [ 'rest_framework.permissions.DjangoModelPermissionsOrAnonReadOnly', ], # デフォルトの認証クラス 'DEFAULT_AUTHENTICATION_CLASSES': [ 'rest_framework_jwt.authentication.JSONWebTokenAuthentication', 'rest_framework.authentication.SessionAuthentication', 'rest_framework.authentication.BasicAuthentication', ], # デフォルトのフィルターバックエンド 'DEFAULT_FILTER_BACKENDS': [ 'django_filters.rest_framework.DjangoFilterBackend', ], # デフォルトのスロットルクラス 'DEFAULT_THROTTLE_CLASSES': [ 'rest_framework.throttling.AnonRateThrottle', 'rest_framework.throttling.UserRateThrottle', ], # デフォルトのスロットルの制限 'DEFAULT_THROTTLE_RATES': { 'anon': '100/day', 'user': '5/sec', } } # 追加 if DEBUG: # CORS関連の設定追加。 INSTALLED_APPS += ['corsheaders'] MIDDLEWARE = ['corsheaders.middleware.CorsMiddleware'] + MIDDLEWARE CORS_ORIGIN_ALLOW_ALL = True DRF: ルーティング k2/api/k2/urls.py from django.urls import path, include from . import views, viewsets from rest_framework.routers import DefaultRouter router = DefaultRouter() router.register('users', viewsets.UserViewSet) app_name = 'k2' urlpatterns = [ path('', include(router.urls)), ] k2/api/api/urls.py from django.contrib import admin from django.urls import path, include urlpatterns = [ path('admin/', admin.site.urls), path('api/', include('k2.urls')), ] DRF: View記述 今回はDjangoデフォルトのUserモデルを使用します。 k2/api/k2/viewsets.py from django.contrib.auth.models import User from serializers import UserSerializer from rest_framework import ( viewsets ) import logging logger = logging.getLogger(__name__) class BaseModelViewSet(viewsets.ModelViewSet): pass class UserViewSet(BaseModelViewSet): queryset = User.objects.all() serializer_class = UserSerializer def list(self, request, *args, **kwargs): return super().list(request, *args, **kwargs) k2/api/k2/serializers.py from rest_framework import serializers from django.contrib.auth.models import User import logging import logging logger = logging.getLogger(__name__) class UserSerializer(serializers.ModelSerializer): class Meta: model = User fields = '__all__' これで、localhost:8000/api/users/にアクセスすると、python manage.py createsuperuserコマンドなどで作成したユーザー一覧が返る。 Vue側 やる事 初期設定 ルーティング 画面実装 ajax投げて表示 Vue: 初期設定 main.jsでライブラリ等を読み込む。 k2/web/src/main.js import Vue from 'vue' import App from './App.vue' import router from './router' import store from './store' import http from '@/plugins/http' import vuetify from './plugins/vuetify' import eventHub from '@/plugins/eventHub' import 'boxicons/css/boxicons.min.css' Vue.config.productionTip = false Vue.use(http) Vue.use(eventHub) new Vue({ router, store, vuetify, render: h => h(App) }).$mount('#app') k2/web/src/App.vue <template> <v-app> <router-view/> </v-app> </template> <script> export default { name: 'App', } </script> axiosのデフォルト設定をしておく。 k2/web/src/plugins/http.js import Vue from 'vue' import axios from 'axios' import router from '@/router' export default { install: function (Vue, options) { const http = axios.create({ baseURL: 'http://localhost:8000/', xsrfCookieName: 'csrftoken', xsrfHeaderName: 'X-CSRFTOKEN', timeout: 10000, }) Vue.prototype.$axios = http } } k2/web/src/plugins/vuetify.js import Vue from 'vue' import Vuetify from 'vuetify/lib' Vue.use(Vuetify) export default new Vuetify({ }) 孫コンポーネントとか遠い関係のコンポーネントのイベントを検知するためのもの。 k2/web/src/plugins/eventHub.js import Vue from 'vue' const eventHub = { install: function (Vue, options) { Vue.prototype.$eventHub = new Vue() } } export default eventHub eslintがうるさいので最低限だまらせておく。 k2/web/.eslintrc.js ... rules: { 'no-console': process.env.NODE_ENV === 'production' ? 'warn' : 'off', 'no-debugger': process.env.NODE_ENV === 'production' ? 'warn' : 'off', 'no-tabs': 'off', 'indent': 'off', 'comma-dangle': 'off', 'eol-last': 'off', 'no-unused-vars': 'off', 'no-mixed-spaces-and-tabs': 'off', 'no-unneeded-ternary': 'off', 'vue/no-unused-components': 'off', 'no-multi-spaces': 'off' }, ... Vue: ルーティング 全体のルート k2/web/src/router/index.js import Vue from 'vue' import VueRouter from 'vue-router' import pages from './pages' Vue.use(VueRouter) const routes = [...pages] const router = new VueRouter({ mode: 'history', routes }) export default router ページ関連のルート k2/web/src/router/pages.js import { Home } from '@/views/index' const routes = [ { path: '/', name: 'Home', component: Home }, ] export default routes Vue: 画面実装 index.jsにコンポーネントを纏める。 k2/web/src/views/index.js export { default as Home } from './Home' Home画面のコンポーネント。今回はコンポーネント分けないで直書きしちゃいます。 k2/web/src/views/Home.vue <template> <div id="home_wrap"> <v-card flat tile > <div v-if="loading"> <Loading/> </div> <div v-else> <v-card-title> users </v-card-title> <v-list three-line> <template v-for="(user, i) in users"> <v-list-item :key="i" > <v-list-item-avatar> <i class='bx bxs-user'></i> </v-list-item-avatar> <v-list-item-content> <v-list-item-title> {{ user.username }} </v-list-item-title> <v-list-item-subtitle> id: {{ user.id }} </v-list-item-subtitle> </v-list-item-content> </v-list-item> </template> </v-list> </div> </v-card> </div> </template> <script> import Loading from '../components/Loading' export default { name: 'Home', components: { Loading, }, data: () => ({ loading: true, users: [], }), mounted () { this.getUsers() }, methods: { getUsers () { this.loading = true this.$axios({ url: '/api/users/', method: 'GET', }) .then(res => { console.log(res) this.users = res.data this.loading = false }) .catch(e => { console.log(e) this.loading = false }) } } } </script> <style lang="scss" scoped> #home_wrap { width: 1200px; margin: 0 auto; } </style> ロード画面のコンポーネント k2/web/src/components/Loading.vue <template> <div class="loading_wrap"> <vue-loading type="cylon" color="#aaa" :size="{ width: '50px', height: '50px' }" ></vue-loading> </div> </template> <script> import { VueLoading } from 'vue-loading-template' export default { name: 'Loading', components: { VueLoading, }, } </script> <style lang="scss" scoped> .loading_wrap { width: 50px; height: 50px; margin: 50px auto; } </style> まとめ こんな感じの画面になりました。

タスクスケジューラなど定期実行のツールを使用しました。 ~時間通りに生活するのはしんどいので対策です。~ Windowsユーザーです。タスクスケジューラの様な定期実行を行うプログラムを記録しました。 普段やっている作業が基本パソコン エクセルやパワポの様なツールを使うことが多い 普段から情報収集のためのサイトを多く活用している これがいつもの生活なんですけど、あんまり楽しくなさそうですよね( ；∀；) そんな毎日を送っていると　"この作業にかける時間無駄だな”　、とか思ってくるんです。 マイクロソフトのサービスを利用している人は聞いたことある　"vba"　ってありますよね。 あの言語を使って多少は改善してきてるんですけど、コーディングするまで大変ですよね。 でも言語さえ覚えれば便利なんでけど、ここで思い始めてきたのが　 ”もう、アプリ開かなくても動いてくれないかな”　 いやいや,それは素晴らしいけど...,甘えじゃないか…? 甘えだと分かってるんですけど、そう思い始めるとその方向に全力で舵を切ります!!!! て思ったらかなり常識的な話で、ただ私が知らないだけでした。 パソコンさえ起動させていれば、勝手にやってくれるサービスは普通にありました。 ツールが多すぎてどれにしようか迷ったくらいです。それでは参りましょう。 大定番。 というわけで,ここからはツールの詳細についてです。普段使っているパソコンのOSがWindowsであれば,必ず搭載されているツール,それがタスクスケジューラです。どういうものかというと, プログラムファイル,アプリケーションを起動する 起動する日時が設定できる この2つが主な目的です,基本的にタスクスケジューラというツールがツールバーだったり,デスクトップには表示されていないので,検索機能でそのままタスクスケジューラと入力してください。そうすると,タスクスケジューラを起動できます。ここから, タスクを作成 使用するプログラムの参照 時間と日付を決める テスト実行 こういった流れで,設定します。実際に使用したのが,ExcelやWordをタスク設定し,毎日10:00に起動し,報告書や集計表を作成していく様にしました。起動する前にパソコンがログインしている状態である事が条件です。便利というより,楽になるといった表現があってるかなと思います。ここまで如何でしょうか。ただ時間にファイル開くだけならそんなに便利でもないか…。と思われがちですが,このツールの最大の特性は, 他のアプリケーションを使用していようが,必ず起動するところです。 これに尽きます,実際プログラミングで賄えると思われがちですが,必ずそのファイルがアクティブな状態である事が条件になると思います。タスクスケジューラは,OSに付随されたツールなので,パソコン自体を開いていれば必ず作動します。Windowsユーザーには便利なので,是非お勧めです。続けてはこちら。 プログラムで動かす方法。 イメージはこんな感じです。 Pexels.com 続けて,プログラムで時間ごとに作動するフローを作る方法です。一応ツールを使用する前に,このやり方で話を聞いていたので,こっちも掲載します。使用するのはExcel VBAです。念のためにPythonもやり方を載せてみます。 private sub 自動実行用マクロ() msgbox "メッセージを表示します","メッセージ",vbyesnocancel end sub() private sub 自動実行用() application.ontime now + timevalue("00:00:15"), "自動実行用マクロ" end sub() こんな感じです。最初のマクロはメッセージボックスを表示するだけです。それを自動実行用というマクロが作動すると,15秒後に自動実行用マクロを実行します。この機能については,Excelブックを開いておかないといけません,それに今回のマクロはシンプルですが,複雑なマクロだとエラーが起きたり,アプリケーションが応答しない事があるので,都度確認する事をお勧めします。それでは次行きましょう。 一応Pythonも掲載されていたので,載せます。 Pythonは機械学習,深層学習に特化したライブラリを持つプログラミング言語です。 Pexels.com 次はPythonです。ここでは開発環境が必要になります,今回はコードだけ載せることになります。ご利用されていない場合は,スキップしてください。 #scheduleライブラリ #定期実行 pip install schedule import time import schedule #関数を設定 def job(): print("I'm working...") #schedule.every()で実行する日程を決める schedule.every(10).minutes.do(job) schedule.every().hour.do(job) schedule.every().day.at("10:30").do(job) schedule.every().monday.do(job) schedule.every().wednesday.at("13:15").do(job) schedule.every().minute.at(":17").do(job) while True: schedule.run_pending() time.sleep(1) scheduleライブラリをインストールします。外部パッケージなのでpipでインストールします。timeモジュールは時間プログラムの中の時間経過を把握できます,scheduleモジュールは実行する時間を時間軸や曜日で指定できます。こちらもプログラムが起動しなければ動かないのでご注意ください。最後はツールの紹介です。Google Apps Scriptは勉強中ですが,エラー続きなので今放置してます,再開したとき,報告します。MacOSに関してはユーザーではないので,割愛します。 タスクスケジューラ類似ツール。 猫の手も借りたい,そんな状況です。 Pexels.com 最後は類似ツールをいくつか掲載します。タスクスケジューラとほぼ変わらないものもあるので,ご興味あればどうぞ。 https://www.visualcron.com/scheduler.aspx https://zapier.com/ VisualCronはほとんどタスクスケジューラと変わらないツールみたいです,私がMicrosoftユーザーなのでタスクスケジューラを使っているので導入は控えています。zapierは定期実行より自動化ツールとして多くの人に使われているようです。私もサイト制作においてお世話になりました。その話はまた別の記事で,お話します。それでは,今回はこの辺で…。See you...

日本語で流行っているWikipediaカテゴリーと英語で流行っているWikipediaカテゴリーがどう違うか分析します。 前提 wikipedia記事のボリュームが正規化された分布上のどこにあるかによって流行を定義。 (後述のコードで具体的に書く。) 追記: 流行っているというよりは、「重視されている」という言葉のほうがニュアンスとして近いかもしれません。 前処理コード wikipedia dumpを展開 wikiextractorでdumpからtextを抽出します。 https://github.com/attardi/wikiextractor python -m wikiextractor.WikiExtractor <Wikipedia dump file> 記事ボリュームの抽出 import os from tqdm import tqdm def read_all_doc(root="../en_wiki/text"): out = [] for d in tqdm(os.listdir(root)): for fname in os.listdir(os.path.join(root, d)): if fname.startswith("wiki_"): with open(os.path.join(root, d, fname)) as f: lines = f.readlines() x = read_doc(lines) out += x return out def read_doc(doc): current = None out = [] for line in doc: if line.startswith("<doc"): line = line.split('"') title = line[-2] current = [title, 0] elif line.startswith("</doc"): out.append(tuple(current)) current = None else: current[1] += len(line) return out def format_data(data): out = {"title": [], "volume": []} for x1, x2 in data: out["title"].append(x1) out["volume"].append(x2) return out if __name__ == "__main__": import pandas as pd out = format_data(read_all_doc("../en_wiki/text")) pd.DataFrame(out).to_csv("data_en.csv", index=False) out = format_data(read_all_doc("../ja_wiki/text")) pd.DataFrame(out).to_csv("data_ja.csv", index=False) categorylinks, lanklinks, pageの抽出 sqlファイルをmysqlに入れるのは時間がかかりすぎるので正規表現で抽出します。 (本当は関数として抽象化しておきたいですが、とりあえず実行できるものを置いています。) # <scriptname>.py > categorylink.csv import re import mmap import os filename = "../jawiki-20210501-categorylinks.sql" file_size = os.stat(filename).st_size chunk_size = 2**32 regex = re.compile(rb"\(([0-9]+),(.+?),") max_length = 128 matches = [] f = open(filename, "rb") print("from_id,to_title") for i in range(0, file_size, chunk_size - max_length + 1): length = chunk_size if i + chunk_size <= file_size else file_size - i m = mmap.mmap(f.fileno(), length=length, offset=i, access=mmap.ACCESS_READ) for match in regex.finditer(m): if not (match.end() == len(m) and len(match.group()) < max_length and length == chunk_size): x1 = match.group(1).decode('utf-8') x2 = match.group(2).decode('utf-8') print("{},{}".format(x1, x2)) m.close() f.close() # <scriptname>.py > langlink_jaen.csv import re import mmap import os filename = "../jawiki-20210501-langlinks.sql" file_size = os.stat(filename).st_size chunk_size = 2**32 regex = re.compile(rb"\(([0-9]+),(.+?),(.+?),") max_length = 128 matches = [] f = open(filename, "rb") print("ja_id,en_title") for i in range(0, file_size, chunk_size - max_length + 1): length = chunk_size if i + chunk_size <= file_size else file_size - i m = mmap.mmap(f.fileno(), length=length, offset=i, access=mmap.ACCESS_READ) for match in regex.finditer(m): if not (match.end() == len(m) and len(match.group()) < max_length and length == chunk_size): x1 = match.group(1).decode('utf-8') x2 = match.group(2).decode('utf-8')[1:-1] try: x3 = match.group(3).decode('utf-8')[:-1] except: continue if x2 == "en": print("{},{}".format(x1, x3)) m.close() f.close() # <scriptname>.py > id_title_ja.csv import re import mmap import os filename = "../jawiki-20210501-page.sql" file_size = os.stat(filename).st_size chunk_size = 2**32 regex = re.compile(rb"\(([0-9]+),[0-9+],(.+?),") max_length = 128 matches = [] f = open(filename, "rb") print("id,title") for i in range(0, file_size, chunk_size - max_length + 1): length = chunk_size if i + chunk_size <= file_size else file_size - i m = mmap.mmap(f.fileno(), length=length, offset=i, access=mmap.ACCESS_READ) for match in regex.finditer(m): if not (match.end() == len(m) and len(match.group()) < max_length and length == chunk_size): print("{},{}".format( match.group(1).decode('utf-8'), match.group(2).decode('utf-8'))) m.close() f.close() データをマージする import os import pandas as pd dataroot = "./data" fnames = ["data_en.csv", "data_ja.csv", "id_title_ja.csv", "langlink_jaen.csv"] data = [pd.read_csv(os.path.join(dataroot, fname)) for fname in fnames] data[2]["title"] = [str(x)[1:-1].replace("_", " ") for x in data[2]["title"]] data[3]["en_title"] = [str(x)[1:-1] for x in data[3]["en_title"]] df = pd.merge(data[0], data[3], left_on='title', right_on='en_title') df2 = pd.merge(data[1], data[2], left_on='title', right_on='title') df = pd.merge(df, df2, left_on="ja_id", right_on="id") df.to_csv("data.csv", index=False) wikipediaページ番号を整形 ページ番号を0から順に対応付けなおします。これは処理の都合上のもので、本質的ではありません。 import json import pandas as pd df1 = pd.read_csv("../data/categorylink.csv") df3 = pd.read_csv("../data/data.csv") ids = set(df3["id"]) pdict = {} count = 0 for pid in df1["from_id"]: if pid in ids and pid not in pdict: pdict[pid] = count count += 1 with open("page2datanum.json", "w") as f: json.dump(pdict, f) categoryに番号を振る import pandas as pd df = pd.read_csv("./data/categorylink.csv") out = [] for i, x in enumerate(sorted(df["to_title"].unique())): out.append({"cid": i, "ctitle": x}) pd.DataFrame(out).to_csv("categories.csv", index=False) 特徴量設計する wikipediaカテゴリーをonehot特徴量とするとさすがにメモリが足らんのでスパース化して保存しておきます。 import numpy as np import pandas as pd from tqdm import tqdm import pickle from scipy.sparse import dok_matrix df1 = pd.read_csv("./data/categorylink.csv") df2 = pd.read_csv("./data/categories.csv") df3 = pd.read_csv("./data/data.csv") ids = set(df3["id"]) cdict = dict([(ctitle, cid) for cid, ctitle in zip(df2["cid"], df2["ctitle"])]) row = [0 for _ in range(df2.shape[0])] rowsize = 2004011 #記事数に合わせて修正する colsize = len(row) pdict = {} count = 0 for pid in df1["from_id"]: if pid in ids and pid not in pdict: pdict[pid] = count count += 1 out = dok_matrix((count, colsize), dtype=np.int8) for pid, title in tqdm(zip(df1["from_id"], df1["to_title"])): if int(pid) in ids: out[pdict[pid], cdict[title]] = 1 with open("features.pkl", "wb") as f: pickle.dump(out, f) データの整形 処理の都合上として変換したIDで整形します。 import json import pandas as pd df1 = pd.read_csv("../data/data.csv") with open("../data/page2datanum.json") as f: p2n = json.load(f) m = max(list(p2n.items()), key=lambda x: x[1])[1] out = [None for _ in range(m+1)] first_d = None for i, d in df1.iterrows(): try: key = p2n[str(d["id"])] except KeyError: continue if first_d is None: first_d = d out[key] = d first_d["title_x"] = None first_d["volume_x"] = 0, first_d["ja_id"] = 0, first_d["en_title"] = None, first_d["title_y"] = None first_d["volume_y"] = 0 first_d["id"] = 0 for i, x in enumerate(out): if x is None: out[i] = first_d pd.DataFrame(out).to_csv("data_fixed.csv", index=False) jupyterで実行 In [1]: import pickle import pandas as pd with open("../data/features.pkl", "rb") as f: X = pickle.load(f) df = pd.read_csv("../data/data_fixed.csv") In [2]: import numpy as np df["volume_x_ln"] = [np.log(x+1) for x in df["volume_x"]] df["volume_y_ln"] = [np.log(x+1) for x in df["volume_y"]] df[["volume_x_ln", "volume_y_ln"]].hist(bins=100) Out In [3]: from sklearn.preprocessing import StandardScaler data_scl = StandardScaler().fit_transform(df[["volume_x_ln", "volume_y_ln"]]) df["volume_x_ln_z"] = data_scl[:, 0] df["volume_y_ln_z"] = data_scl[:, 1] df[["volume_x_ln_z", "volume_y_ln_z"]].hist(bins=100) Out In [4]: import matplotlib.pyplot as plt plt.scatter(df["volume_x_ln_z"], df["volume_y_ln_z"]) Out このプロットから推測できるのは、 英語圏と日本語圏で共通して流行しているカテゴリー。 英語圏で日本語圏より流行しているカテゴリー。 日本語圏で英語圏より流行しているカテゴリー。 の3つに分割されていることです。 In [5]: from sklearn.linear_model import LogisticRegression y = df["volume_x_ln_z"] > df["volume_y_ln_z"] clf = LogisticRegression().fit(X, y) In [6]: df_cats = pd.read_csv("../data/categories.csv") cats = df_cats["ctitle"] In [7]: print("英語圏で日本語圏より重要視されるジャンル Top 100") sorted([(x, y) for x, y in zip(cats, clf.coef_[0])], key=lambda x: x[1], reverse=True)[:100] Out [("'化学式の曖昧さ回避'", 8.78965480123576), ("'2016年リオデジャネイロオリンピックの選手団'", 6.4868899626526755), ("'国別の空港の一覧'", 6.192994506435434), ("'ウィキデータ項目に接続されているリダイレクト'", 5.9852466820386425), ("'全米オープン_(テニス)'", 5.739147053420267), ("'ウィンブルドン選手権'", 5.494081617029156), ("'全豪オープン'", 5.044413803151733), ("'2012年ロンドンオリンピックの選手団'", 4.959988235376474), ("'朝鮮民主主義人民共和国の鉄道路線'", 4.953440370793841), ("'2018年平昌オリンピックの選手団'", 4.169210912094708), ("'NFLのシーズン'", 3.963965273829145), ("'2020年東京オリンピックの選手団'", 3.784828296567487), ("'医療翻訳プロジェクト'", 3.644218696183388), ("'ケニアのカウンティ'", 3.537512552672342), ("'各年の映画'", 3.525023577920066), ("'カンヌ国際映画祭'", 3.5091493203184814), ("'バンコクの区'", 3.4963174393787813), ("'多面体関連のスタブ項目'", 3.463447271219896), ("'モスクワ地下鉄の駅'", 3.430100714848111), ("'両江道の鉄道駅'", 3.4093101736218943), ("'各年の文学'", 3.3704473983916197), ("'ベトナムの市社'", 3.3415948871794687), ("'一様多面体'", 3.274475647260854), ("'長沙地下鉄の駅'", 3.2725114510074933), ("'ユダヤ系アメリカ人の野球選手'", 3.2180785734223054), ("'メトロノース鉄道の鉄道駅'", 3.1521912615497323), ("'2004年アテネオリンピックの選手団'", 3.1497131786928008), ("'江原道の鉄道駅_(北)'", 3.1162845268059995), ("'翻訳中'", 3.050924761475756), ("'アフリカのサッカーナショナルチーム'", 2.983501250381955), ("'ブルキナファソの県'", 2.9750502199478452), ("'全仏オープン'", 2.941470855135487), ("'歌手に関するサブスタブ'", 2.935080098470459), ("'イングランドのサッカー競技施設'", 2.920380086615955), ("'テレビ番組に関するサブスタブ'", 2.917601276840637), ("'各年のツール・ド・フランス'", 2.8523732643162436), ("'WWEのPPVイベント'", 2.817326794329752), ("'キリル文字'", 2.8055971062678924), ("'北中米カリブ海のサッカーナショナルチーム'", 2.80060747883836), ("'2012年の日本の各サッカークラブ'", 2.800395133712363), ("'ドイツの主要都市'", 2.799929355362397), ("'アメリカ合衆国の副知事の一覧'", 2.7370952721478043), ("'フランス陸軍の連隊'", 2.7297753302004155), ("'K._ラインムートが発見した天体'", 2.722889361157496), ("'アカデミー賞'", 2.7193328653915643), ("'DCコミックスの登場人物'", 2.711111261472102), ("'タイの県'", 2.7072267395449647), ("'銀行の一覧'", 2.6634085771235565), ("'LDH所属のアーティスト・歌手'", 2.6295816150066047), ("'咸鏡北道の鉄道駅'", 2.607147544748716), ("'台湾のテレビドラマ'", 2.5405701659634317), ("'都市の一覧'", 2.535954389587591), ("'イギリスの航空母艦'", 2.5038940096540063), ("'ルーマニアの県'", 2.4829255537462775), ("'マドンナの楽曲'", 2.461558238490475), ("'サッカーアルバニア代表選手'", 2.452063619405983), ("'ギリシャのサッカークラブ'", 2.447949256145987), ("'レッスルマニア'", 2.4401571457570794), ("'特筆性の基準を満たしていないおそれのある記事/2015年10月'", 2.4278335710965275), ("'咸鏡南道の鉄道駅'", 2.3879674843762806), ("'アルペンスキー世界選手権'", 2.377562580184821), ("'ジョンソンの立体'", 2.3649538983140053), ("'ドイツ語版ウィキペディアからの翻訳を必要とする選り抜き記事'", 2.35065373581786), ("'MLB所属チームの歴代監督一覧'", 2.34203026126281), ("'アメリカ合衆国のテレビ局'", 2.3408765641573455), ("'3等星'", 2.335035616106319), ("'シンガポールの鉄道駅'", 2.329221204684235), ("'アメリカ合衆国の列車'", 2.3244495438131842), ("'合気道家'", 2.3111974199141527), ("'ニューサウスウェールズ州の地方公共団体'", 2.308923757664836), ("'スタテンアイランド鉄道の駅'", 2.302326646100214), ("'北京地下鉄の路線'", 2.301391739237924), ("'英語版ウィキペディアからの翻訳を必要とする選り抜き記事'", 2.2588628019214263), ("'ネパールの郡'", 2.254855200070133), ("'コートジボワールの州'", 2.2505405415473425), ("'工藤静香の楽曲'", 2.243106209268534), ("'MLBオールスターゲーム'", 2.242207780995818), ("'UKPARL識別子が指定されている記事'", 2.2347560870938588), ("'バンベルク郡'", 2.2260537848125868), ("'日本の女子プロレスラー'", 2.2228866587228753), ("'イタリア系アメリカ人のボクサー'", 2.218979576333191), ("'平安北道の鉄道駅'", 2.216703613358715), ("'6世紀中国の女性'", 2.2052280279757004), ("'地域別野鳥一覧'", 2.1812595924783604), ("'聖書に登場する地名'", 2.173145325060096), ("'国際連合安全保障理事会決議の一覧'", 2.1721208013233246), ("'国・地域別のLGBTの権利'", 2.1721060993913395), ("'1929年発見の天体'", 2.1714758882796144), ("'母音'", 2.168575820861319), ("'ベトナムの省'", 2.168310637306939), ("'月の海'", 2.1670244380744546), ("'東亜プランのゲームソフト'", 2.1616478293778583), ("'アジアのサッカーナショナルチーム'", 2.160684396316022), ("'数学者関連のスタブ項目'", 2.1570669971347725), ("'ゴールデングローブ賞'", 2.1570660870654432), ("'アメリカ合衆国の野球場'", 2.156377367097728), ("'プエルトリコのボクサー'", 2.1497232942810847), ("'シンガポール関連のスタブ項目'", 2.135975452827605), ("'2012年生_(競走馬)'", 2.131448969218936), ("'ノヴァーラ県のコムーネ'", 2.1305873465559526)] In [8]: print("日本語圏で英語圏より重要視されるジャンル Top 100") sorted([(x, y) for x, y in zip(cats, clf.coef_[0])], key=lambda x: x[1], reverse=False)[:100] [("'ヘッセン州の行政区画'", -6.656503838669296), ("'日本の重力式コンクリートダム'", -4.677250334777994), ("'戦国武将'", -4.597537701941709), ("'中央競馬の競走'", -4.263524600449024), ("'ウィキデータにある日本相撲協会識別子'", -4.062260864597189), ("'良質な記事'", -4.000476177005945), ("'各年のメジャーリーグベースボール'", -3.773628832673474), ("'木星の衛星'", -3.7695251618472727), ("'日本の野球選手'", -3.7338732590073467), ("'秋元康が制作した楽曲'", -3.58052170022946), ("'土星の衛星'", -3.564344817840197), ("'サッカー日本代表'", -3.4810250620841643), ("'関東地方の鉄道路線'", -3.4658272948772377), ("'NBAのシーズン'", -3.4289618605720014), ("'日本の男子バスケットボール選手'", -3.4063404117333618), ("'Jリーグのシーズン'", -3.3893788224952917), ("'各年の鉄道'", -3.3530982470985373), ("'バラオ級潜水艦'", -3.3513589291698094), ("'フィニステール県のコミューン'", -3.3252031956331924), ("'カリフォルニア州の都市'", -3.285739146248556), ("'ノルトライン＝ヴェストファーレン州の行政区画'", -3.250677146009764), ("'戦国大名'", -3.2320592814495157), ("'日本民間放送連盟会員'", -3.189809433685342), ("'釜山交通公社の鉄道駅'", -3.1679351063913903), ("'ギリシャの県'", -3.1505736330497793), ("'X-ファイルのエピソード'", -3.143230057702075), ("'各年の日本のテレビ放送'", -3.12063972926185), ("'イシガメ科'", -3.1067526066382816), ("'ニーダーザクセン州の行政区画'", -3.1012895614535405), ("'各年の航空'", -3.0720082317291264), ("'日本の高速道路'", -3.0642048869297973), ("'オーストリアの郡'", -3.0625741007306297), ("'鹿児島本線'", -3.0484239636700425), ("'日本人の発見した小惑星'", -3.0396341664666595), ("'西日本旅客鉄道の鉄道駅'", -3.027178466090375), ("'Jリーグ公式戦開催競技場'", -3.006864940865679), ("'紀元前8世紀の各年'", -2.993774925054372), ("'近畿地方の鉄道路線'", -2.9823024867144476), ("'ガトー級潜水艦'", -2.981977664176823), ("'バーデン＝ヴュルテンベルク州の行政区画'", -2.9629892972302283), ("'ニーダーザクセン州のザムトゲマインデ'", -2.92071800294315), ("'日本車輌製造製の電車'", -2.858201187754197), ("'日本の女子サッカー選手'", -2.8499437513543278), ("'アメリカ合衆国のロマンス作家'", -2.8315523852717592), ("'日本貨物鉄道の鉄道駅'", -2.811665054295694), ("'ロードレース世界選手権のレースレポート'", -2.801282297215031), ("'日本の国立大学'", -2.787545917331532), ("'JFL'", -2.787134758305849), ("'日本の囲碁棋士'", -2.7832959823790113), ("'日本レコード大賞'", -2.76257516826387), ("'日本女子サッカーリーグ'", -2.753583222667426), ("'アメリカ合衆国の州'", -2.726731988835921), ("'日本国有鉄道の鉄道路線'", -2.693645938296784), ("'秀逸な記事'", -2.6845217634620964), ("'日本の公立大学'", -2.6762057562387387), ("'中部地方の鉄道路線'", -2.674953331834506), ("'トレント県のコムーネ'", -2.6686021539216647), ("'ロワール＝アトランティック県のコミューン'", -2.6616413161428096), ("'三菱重工業長崎造船所が建造した船舶'", -2.6283050174506575), ("'サッカー日本女子代表選手'", -2.6163726838460932), ("'ソンドリオ県のコムーネ'", -2.6008733553088375), ("'日本のサッカー競技施設'", -2.597480186206064), ("'EC_2.1.1'", -2.5793494034379476), ("'東急車輛製造製の電車'", -2.5422948670796703), ("'日本の廃止市町村一覧'", -2.531928482570778), ("'ウィキデータにあるJ.League_player_ID'", -2.5306233382664836), ("'中国の囲碁棋士'", -2.5232228084965076), ("'カルヴァドス県のコミューン'", -2.5174266924755795), ("'チャンパ王'", -2.5155379760759424), ("'韓国のサッカー選手'", -2.512717973478731), ("'地域リーグ_(サッカー)'", -2.5092659938331185), ("'西日本旅客鉄道の廃駅'", -2.498118817838307), ("'日本ゴルフツアー'", -2.491887432773248), ("'プログレス補給船'", -2.4850912198888517), ("'シロンスク・ピャスト家'", -2.4597853002112835), ("'Every_Little_Thingの楽曲'", -2.4473470561739035), ("'サンフレッチェ広島F.Cの選手'", -2.4368795566513466), ("'同名の条約'", -2.4122489600231236), ("'部首'", -2.408524332052575), ("'ノルウェーのスキージャンプ選手'", -2.3908719135063885), ("'日本の特急列車'", -2.3825918009949816), ("'古代エジプトの王朝'", -2.371403584625034), ("'モササウルス科'", -2.371160768526831), ("'プライムタイム・エミー賞'", -2.370694238733858), ("'川崎重工業製の電車'", -2.3702464677412576), ("'国際宇宙ステーションへのプログレス補給船'", -2.3610100085465904), ("'日本のサッカー選手'", -2.357607796365293), ("'ブルーグラス'", -2.3555303538295034), ("'Reflistで4列以上を指定しているページ'", -2.3552855477384127), ("'クイナ科'", -2.3533545440603167), ("'現在活動していない日本のサッカークラブ'", -2.318549520809771), ("'日本の町・字のスタブ項目'", -2.31738325991341), ("'日本のアニメスタジオ'", -2.3148947978511174), ("'PRIDEの大会'", -2.2920916445971407), ("'アメリカ合衆国のジャズ・バンド'", -2.287166077883642), ("'日立製作所製の電車'", -2.2841706145543657), ("'東京地下鉄の鉄道駅'", -2.283909048734933), ("'相撲に関するスタブ'", -2.2834021869957764), ("'日本の女子サッカークラブ'", -2.278460378256027), ("'各年のプロ野球日本シリーズ'", -2.2772900084762076 例えば、「ヘッセン州の行政区画」については、ドイツについて専門的に書いている日本語ユーザーの熱心さを表しています。 https://ja.wikipedia.org/wiki/%E5%88%A9%E7%94%A8%E8%80%85:%E6%B1%B2%E5%B9%B3 リンク https://github.com/sugiyamath/wikipedia_category_analysis

日本語で重視されているWikipediaカテゴリーと英語で重視されているWikipediaカテゴリーがどう違うか分析します。 前提 wikipedia記事のボリュームが正規化された分布上のどこにあるかによって重視されているかどうかを定義。(後述のコードで具体的に書く。) 前処理コード wikipedia dumpを展開 wikiextractorでdumpからtextを抽出します。 https://github.com/attardi/wikiextractor python -m wikiextractor.WikiExtractor <Wikipedia dump file> 記事ボリュームの抽出 import os from tqdm import tqdm def read_all_doc(root="../en_wiki/text"): out = [] for d in tqdm(os.listdir(root)): for fname in os.listdir(os.path.join(root, d)): if fname.startswith("wiki_"): with open(os.path.join(root, d, fname)) as f: lines = f.readlines() x = read_doc(lines) out += x return out def read_doc(doc): current = None out = [] for line in doc: if line.startswith("<doc"): line = line.split('"') title = line[-2] current = [title, 0] elif line.startswith("</doc"): out.append(tuple(current)) current = None else: current[1] += len(line) return out def format_data(data): out = {"title": [], "volume": []} for x1, x2 in data: out["title"].append(x1) out["volume"].append(x2) return out if __name__ == "__main__": import pandas as pd out = format_data(read_all_doc("../en_wiki/text")) pd.DataFrame(out).to_csv("data_en.csv", index=False) out = format_data(read_all_doc("../ja_wiki/text")) pd.DataFrame(out).to_csv("data_ja.csv", index=False) categorylinks, lanklinks, pageの抽出 sqlファイルをmysqlに入れるのは時間がかかりすぎるので正規表現で抽出します。 (本当は関数として抽象化しておきたいですが、とりあえず実行できるものを置いています。) # <scriptname>.py > categorylink.csv import re import mmap import os filename = "../jawiki-20210501-categorylinks.sql" file_size = os.stat(filename).st_size chunk_size = 2**32 regex = re.compile(rb"\(([0-9]+),(.+?),") max_length = 128 matches = [] f = open(filename, "rb") print("from_id,to_title") for i in range(0, file_size, chunk_size - max_length + 1): length = chunk_size if i + chunk_size <= file_size else file_size - i m = mmap.mmap(f.fileno(), length=length, offset=i, access=mmap.ACCESS_READ) for match in regex.finditer(m): if not (match.end() == len(m) and len(match.group()) < max_length and length == chunk_size): x1 = match.group(1).decode('utf-8') x2 = match.group(2).decode('utf-8') print("{},{}".format(x1, x2)) m.close() f.close() # <scriptname>.py > langlink_jaen.csv import re import mmap import os filename = "../jawiki-20210501-langlinks.sql" file_size = os.stat(filename).st_size chunk_size = 2**32 regex = re.compile(rb"\(([0-9]+),(.+?),(.+?),") max_length = 128 matches = [] f = open(filename, "rb") print("ja_id,en_title") for i in range(0, file_size, chunk_size - max_length + 1): length = chunk_size if i + chunk_size <= file_size else file_size - i m = mmap.mmap(f.fileno(), length=length, offset=i, access=mmap.ACCESS_READ) for match in regex.finditer(m): if not (match.end() == len(m) and len(match.group()) < max_length and length == chunk_size): x1 = match.group(1).decode('utf-8') x2 = match.group(2).decode('utf-8')[1:-1] try: x3 = match.group(3).decode('utf-8')[:-1] except: continue if x2 == "en": print("{},{}".format(x1, x3)) m.close() f.close() # <scriptname>.py > id_title_ja.csv import re import mmap import os filename = "../jawiki-20210501-page.sql" file_size = os.stat(filename).st_size chunk_size = 2**32 regex = re.compile(rb"\(([0-9]+),[0-9+],(.+?),") max_length = 128 matches = [] f = open(filename, "rb") print("id,title") for i in range(0, file_size, chunk_size - max_length + 1): length = chunk_size if i + chunk_size <= file_size else file_size - i m = mmap.mmap(f.fileno(), length=length, offset=i, access=mmap.ACCESS_READ) for match in regex.finditer(m): if not (match.end() == len(m) and len(match.group()) < max_length and length == chunk_size): print("{},{}".format( match.group(1).decode('utf-8'), match.group(2).decode('utf-8'))) m.close() f.close() データをマージする import os import pandas as pd dataroot = "./data" fnames = ["data_en.csv", "data_ja.csv", "id_title_ja.csv", "langlink_jaen.csv"] data = [pd.read_csv(os.path.join(dataroot, fname)) for fname in fnames] data[2]["title"] = [str(x)[1:-1].replace("_", " ") for x in data[2]["title"]] data[3]["en_title"] = [str(x)[1:-1] for x in data[3]["en_title"]] df = pd.merge(data[0], data[3], left_on='title', right_on='en_title') df2 = pd.merge(data[1], data[2], left_on='title', right_on='title') df = pd.merge(df, df2, left_on="ja_id", right_on="id") df.to_csv("data.csv", index=False) wikipediaページ番号を整形 ページ番号を0から順に対応付けなおします。これは処理の都合上のもので、本質的ではありません。 import json import pandas as pd df1 = pd.read_csv("../data/categorylink.csv") df3 = pd.read_csv("../data/data.csv") ids = set(df3["id"]) pdict = {} count = 0 for pid in df1["from_id"]: if pid in ids and pid not in pdict: pdict[pid] = count count += 1 with open("page2datanum.json", "w") as f: json.dump(pdict, f) categoryに番号を振る import pandas as pd df = pd.read_csv("./data/categorylink.csv") out = [] for i, x in enumerate(sorted(df["to_title"].unique())): out.append({"cid": i, "ctitle": x}) pd.DataFrame(out).to_csv("categories.csv", index=False) 特徴量設計する wikipediaカテゴリーをonehot特徴量とするとさすがにメモリが足らんのでスパース化して保存しておきます。 import numpy as np import pandas as pd from tqdm import tqdm import pickle from scipy.sparse import dok_matrix df1 = pd.read_csv("./data/categorylink.csv") df2 = pd.read_csv("./data/categories.csv") df3 = pd.read_csv("./data/data.csv") ids = set(df3["id"]) cdict = dict([(ctitle, cid) for cid, ctitle in zip(df2["cid"], df2["ctitle"])]) row = [0 for _ in range(df2.shape[0])] colsize = len(row) pdict = {} count = 0 for pid in df1["from_id"]: if pid in ids and pid not in pdict: pdict[pid] = count count += 1 out = dok_matrix((count, colsize), dtype=np.int8) for pid, title in tqdm(zip(df1["from_id"], df1["to_title"])): if int(pid) in ids: out[pdict[pid], cdict[title]] = 1 with open("features.pkl", "wb") as f: pickle.dump(out, f) データの整形 処理の都合上として変換したIDで整形します。 import json import pandas as pd df1 = pd.read_csv("../data/data.csv") with open("../data/page2datanum.json") as f: p2n = json.load(f) m = max(list(p2n.items()), key=lambda x: x[1])[1] out = [None for _ in range(m+1)] first_d = None for i, d in df1.iterrows(): try: key = p2n[str(d["id"])] except KeyError: continue if first_d is None: first_d = d out[key] = d first_d["title_x"] = None first_d["volume_x"] = 0, first_d["ja_id"] = 0, first_d["en_title"] = None, first_d["title_y"] = None first_d["volume_y"] = 0 first_d["id"] = 0 for i, x in enumerate(out): if x is None: out[i] = first_d pd.DataFrame(out).to_csv("data_fixed.csv", index=False) jupyterで実行 In [1]: import pickle import pandas as pd with open("../data/features.pkl", "rb") as f: X = pickle.load(f) df = pd.read_csv("../data/data_fixed.csv") In [2]: import numpy as np df["volume_x_ln"] = [np.log(x+1) for x in df["volume_x"]] df["volume_y_ln"] = [np.log(x+1) for x in df["volume_y"]] df[["volume_x_ln", "volume_y_ln"]].hist(bins=100) Out In [3]: from sklearn.preprocessing import StandardScaler data_scl = StandardScaler().fit_transform(df[["volume_x_ln", "volume_y_ln"]]) df["volume_x_ln_z"] = data_scl[:, 0] df["volume_y_ln_z"] = data_scl[:, 1] df[["volume_x_ln_z", "volume_y_ln_z"]].hist(bins=100) Out In [4]: import matplotlib.pyplot as plt plt.scatter(df["volume_x_ln_z"], df["volume_y_ln_z"]) Out このプロットから推測できるのは、 英語圏と日本語圏で共通して重視されているカテゴリー。 英語圏で日本語圏より重視されているカテゴリー。 日本語圏で英語圏より重視されているカテゴリー。 の軸に分かれていることです。 In [5]: from sklearn.linear_model import LogisticRegression y = df["volume_x_ln_z"] > df["volume_y_ln_z"] clf = LogisticRegression().fit(X, y) In [6]: df_cats = pd.read_csv("../data/categories.csv") cats = df_cats["ctitle"] In [7]: print("英語圏で日本語圏より重要視されるジャンル Top 100") sorted([(x, y) for x, y in zip(cats, clf.coef_[0])], key=lambda x: x[1], reverse=True)[:100] Out [("'化学式の曖昧さ回避'", 8.78965480123576), ("'2016年リオデジャネイロオリンピックの選手団'", 6.4868899626526755), ("'国別の空港の一覧'", 6.192994506435434), ("'ウィキデータ項目に接続されているリダイレクト'", 5.9852466820386425), ("'全米オープン_(テニス)'", 5.739147053420267), ("'ウィンブルドン選手権'", 5.494081617029156), ("'全豪オープン'", 5.044413803151733), ("'2012年ロンドンオリンピックの選手団'", 4.959988235376474), ("'朝鮮民主主義人民共和国の鉄道路線'", 4.953440370793841), ("'2018年平昌オリンピックの選手団'", 4.169210912094708), ("'NFLのシーズン'", 3.963965273829145), ("'2020年東京オリンピックの選手団'", 3.784828296567487), ("'医療翻訳プロジェクト'", 3.644218696183388), ("'ケニアのカウンティ'", 3.537512552672342), ("'各年の映画'", 3.525023577920066), ("'カンヌ国際映画祭'", 3.5091493203184814), ("'バンコクの区'", 3.4963174393787813), ("'多面体関連のスタブ項目'", 3.463447271219896), ("'モスクワ地下鉄の駅'", 3.430100714848111), ("'両江道の鉄道駅'", 3.4093101736218943), ("'各年の文学'", 3.3704473983916197), ("'ベトナムの市社'", 3.3415948871794687), ("'一様多面体'", 3.274475647260854), ("'長沙地下鉄の駅'", 3.2725114510074933), ("'ユダヤ系アメリカ人の野球選手'", 3.2180785734223054), ("'メトロノース鉄道の鉄道駅'", 3.1521912615497323), ("'2004年アテネオリンピックの選手団'", 3.1497131786928008), ("'江原道の鉄道駅_(北)'", 3.1162845268059995), ("'翻訳中'", 3.050924761475756), ("'アフリカのサッカーナショナルチーム'", 2.983501250381955), ("'ブルキナファソの県'", 2.9750502199478452), ("'全仏オープン'", 2.941470855135487), ("'歌手に関するサブスタブ'", 2.935080098470459), ("'イングランドのサッカー競技施設'", 2.920380086615955), ("'テレビ番組に関するサブスタブ'", 2.917601276840637), ("'各年のツール・ド・フランス'", 2.8523732643162436), ("'WWEのPPVイベント'", 2.817326794329752), ("'キリル文字'", 2.8055971062678924), ("'北中米カリブ海のサッカーナショナルチーム'", 2.80060747883836), ("'2012年の日本の各サッカークラブ'", 2.800395133712363), ("'ドイツの主要都市'", 2.799929355362397), ("'アメリカ合衆国の副知事の一覧'", 2.7370952721478043), ("'フランス陸軍の連隊'", 2.7297753302004155), ("'K._ラインムートが発見した天体'", 2.722889361157496), ("'アカデミー賞'", 2.7193328653915643), ("'DCコミックスの登場人物'", 2.711111261472102), ("'タイの県'", 2.7072267395449647), ("'銀行の一覧'", 2.6634085771235565), ("'LDH所属のアーティスト・歌手'", 2.6295816150066047), ("'咸鏡北道の鉄道駅'", 2.607147544748716), ("'台湾のテレビドラマ'", 2.5405701659634317), ("'都市の一覧'", 2.535954389587591), ("'イギリスの航空母艦'", 2.5038940096540063), ("'ルーマニアの県'", 2.4829255537462775), ("'マドンナの楽曲'", 2.461558238490475), ("'サッカーアルバニア代表選手'", 2.452063619405983), ("'ギリシャのサッカークラブ'", 2.447949256145987), ("'レッスルマニア'", 2.4401571457570794), ("'特筆性の基準を満たしていないおそれのある記事/2015年10月'", 2.4278335710965275), ("'咸鏡南道の鉄道駅'", 2.3879674843762806), ("'アルペンスキー世界選手権'", 2.377562580184821), ("'ジョンソンの立体'", 2.3649538983140053), ("'ドイツ語版ウィキペディアからの翻訳を必要とする選り抜き記事'", 2.35065373581786), ("'MLB所属チームの歴代監督一覧'", 2.34203026126281), ("'アメリカ合衆国のテレビ局'", 2.3408765641573455), ("'3等星'", 2.335035616106319), ("'シンガポールの鉄道駅'", 2.329221204684235), ("'アメリカ合衆国の列車'", 2.3244495438131842), ("'合気道家'", 2.3111974199141527), ("'ニューサウスウェールズ州の地方公共団体'", 2.308923757664836), ("'スタテンアイランド鉄道の駅'", 2.302326646100214), ("'北京地下鉄の路線'", 2.301391739237924), ("'英語版ウィキペディアからの翻訳を必要とする選り抜き記事'", 2.2588628019214263), ("'ネパールの郡'", 2.254855200070133), ("'コートジボワールの州'", 2.2505405415473425), ("'工藤静香の楽曲'", 2.243106209268534), ("'MLBオールスターゲーム'", 2.242207780995818), ("'UKPARL識別子が指定されている記事'", 2.2347560870938588), ("'バンベルク郡'", 2.2260537848125868), ("'日本の女子プロレスラー'", 2.2228866587228753), ("'イタリア系アメリカ人のボクサー'", 2.218979576333191), ("'平安北道の鉄道駅'", 2.216703613358715), ("'6世紀中国の女性'", 2.2052280279757004), ("'地域別野鳥一覧'", 2.1812595924783604), ("'聖書に登場する地名'", 2.173145325060096), ("'国際連合安全保障理事会決議の一覧'", 2.1721208013233246), ("'国・地域別のLGBTの権利'", 2.1721060993913395), ("'1929年発見の天体'", 2.1714758882796144), ("'母音'", 2.168575820861319), ("'ベトナムの省'", 2.168310637306939), ("'月の海'", 2.1670244380744546), ("'東亜プランのゲームソフト'", 2.1616478293778583), ("'アジアのサッカーナショナルチーム'", 2.160684396316022), ("'数学者関連のスタブ項目'", 2.1570669971347725), ("'ゴールデングローブ賞'", 2.1570660870654432), ("'アメリカ合衆国の野球場'", 2.156377367097728), ("'プエルトリコのボクサー'", 2.1497232942810847), ("'シンガポール関連のスタブ項目'", 2.135975452827605), ("'2012年生_(競走馬)'", 2.131448969218936), ("'ノヴァーラ県のコムーネ'", 2.1305873465559526)] In [8]: print("日本語圏で英語圏より重要視されるジャンル Top 100") sorted([(x, y) for x, y in zip(cats, clf.coef_[0])], key=lambda x: x[1], reverse=False)[:100] [("'ヘッセン州の行政区画'", -6.656503838669296), ("'日本の重力式コンクリートダム'", -4.677250334777994), ("'戦国武将'", -4.597537701941709), ("'中央競馬の競走'", -4.263524600449024), ("'ウィキデータにある日本相撲協会識別子'", -4.062260864597189), ("'良質な記事'", -4.000476177005945), ("'各年のメジャーリーグベースボール'", -3.773628832673474), ("'木星の衛星'", -3.7695251618472727), ("'日本の野球選手'", -3.7338732590073467), ("'秋元康が制作した楽曲'", -3.58052170022946), ("'土星の衛星'", -3.564344817840197), ("'サッカー日本代表'", -3.4810250620841643), ("'関東地方の鉄道路線'", -3.4658272948772377), ("'NBAのシーズン'", -3.4289618605720014), ("'日本の男子バスケットボール選手'", -3.4063404117333618), ("'Jリーグのシーズン'", -3.3893788224952917), ("'各年の鉄道'", -3.3530982470985373), ("'バラオ級潜水艦'", -3.3513589291698094), ("'フィニステール県のコミューン'", -3.3252031956331924), ("'カリフォルニア州の都市'", -3.285739146248556), ("'ノルトライン＝ヴェストファーレン州の行政区画'", -3.250677146009764), ("'戦国大名'", -3.2320592814495157), ("'日本民間放送連盟会員'", -3.189809433685342), ("'釜山交通公社の鉄道駅'", -3.1679351063913903), ("'ギリシャの県'", -3.1505736330497793), ("'X-ファイルのエピソード'", -3.143230057702075), ("'各年の日本のテレビ放送'", -3.12063972926185), ("'イシガメ科'", -3.1067526066382816), ("'ニーダーザクセン州の行政区画'", -3.1012895614535405), ("'各年の航空'", -3.0720082317291264), ("'日本の高速道路'", -3.0642048869297973), ("'オーストリアの郡'", -3.0625741007306297), ("'鹿児島本線'", -3.0484239636700425), ("'日本人の発見した小惑星'", -3.0396341664666595), ("'西日本旅客鉄道の鉄道駅'", -3.027178466090375), ("'Jリーグ公式戦開催競技場'", -3.006864940865679), ("'紀元前8世紀の各年'", -2.993774925054372), ("'近畿地方の鉄道路線'", -2.9823024867144476), ("'ガトー級潜水艦'", -2.981977664176823), ("'バーデン＝ヴュルテンベルク州の行政区画'", -2.9629892972302283), ("'ニーダーザクセン州のザムトゲマインデ'", -2.92071800294315), ("'日本車輌製造製の電車'", -2.858201187754197), ("'日本の女子サッカー選手'", -2.8499437513543278), ("'アメリカ合衆国のロマンス作家'", -2.8315523852717592), ("'日本貨物鉄道の鉄道駅'", -2.811665054295694), ("'ロードレース世界選手権のレースレポート'", -2.801282297215031), ("'日本の国立大学'", -2.787545917331532), ("'JFL'", -2.787134758305849), ("'日本の囲碁棋士'", -2.7832959823790113), ("'日本レコード大賞'", -2.76257516826387), ("'日本女子サッカーリーグ'", -2.753583222667426), ("'アメリカ合衆国の州'", -2.726731988835921), ("'日本国有鉄道の鉄道路線'", -2.693645938296784), ("'秀逸な記事'", -2.6845217634620964), ("'日本の公立大学'", -2.6762057562387387), ("'中部地方の鉄道路線'", -2.674953331834506), ("'トレント県のコムーネ'", -2.6686021539216647), ("'ロワール＝アトランティック県のコミューン'", -2.6616413161428096), ("'三菱重工業長崎造船所が建造した船舶'", -2.6283050174506575), ("'サッカー日本女子代表選手'", -2.6163726838460932), ("'ソンドリオ県のコムーネ'", -2.6008733553088375), ("'日本のサッカー競技施設'", -2.597480186206064), ("'EC_2.1.1'", -2.5793494034379476), ("'東急車輛製造製の電車'", -2.5422948670796703), ("'日本の廃止市町村一覧'", -2.531928482570778), ("'ウィキデータにあるJ.League_player_ID'", -2.5306233382664836), ("'中国の囲碁棋士'", -2.5232228084965076), ("'カルヴァドス県のコミューン'", -2.5174266924755795), ("'チャンパ王'", -2.5155379760759424), ("'韓国のサッカー選手'", -2.512717973478731), ("'地域リーグ_(サッカー)'", -2.5092659938331185), ("'西日本旅客鉄道の廃駅'", -2.498118817838307), ("'日本ゴルフツアー'", -2.491887432773248), ("'プログレス補給船'", -2.4850912198888517), ("'シロンスク・ピャスト家'", -2.4597853002112835), ("'Every_Little_Thingの楽曲'", -2.4473470561739035), ("'サンフレッチェ広島F.Cの選手'", -2.4368795566513466), ("'同名の条約'", -2.4122489600231236), ("'部首'", -2.408524332052575), ("'ノルウェーのスキージャンプ選手'", -2.3908719135063885), ("'日本の特急列車'", -2.3825918009949816), ("'古代エジプトの王朝'", -2.371403584625034), ("'モササウルス科'", -2.371160768526831), ("'プライムタイム・エミー賞'", -2.370694238733858), ("'川崎重工業製の電車'", -2.3702464677412576), ("'国際宇宙ステーションへのプログレス補給船'", -2.3610100085465904), ("'日本のサッカー選手'", -2.357607796365293), ("'ブルーグラス'", -2.3555303538295034), ("'Reflistで4列以上を指定しているページ'", -2.3552855477384127), ("'クイナ科'", -2.3533545440603167), ("'現在活動していない日本のサッカークラブ'", -2.318549520809771), ("'日本の町・字のスタブ項目'", -2.31738325991341), ("'日本のアニメスタジオ'", -2.3148947978511174), ("'PRIDEの大会'", -2.2920916445971407), ("'アメリカ合衆国のジャズ・バンド'", -2.287166077883642), ("'日立製作所製の電車'", -2.2841706145543657), ("'東京地下鉄の鉄道駅'", -2.283909048734933), ("'相撲に関するスタブ'", -2.2834021869957764), ("'日本の女子サッカークラブ'", -2.278460378256027), ("'各年のプロ野球日本シリーズ'", -2.2772900084762076 例えば、「ヘッセン州の行政区画」については、ドイツについて専門的に書いている日本語ユーザーの熱心さを表しています。 https://ja.wikipedia.org/wiki/%E5%88%A9%E7%94%A8%E8%80%85:%E6%B1%B2%E5%B9%B3 リンク https://github.com/sugiyamath/wikipedia_category_analysis

はじめに 移動平均線のゴールデンクロスよりもシグナルが早く出るということで人気のテクニカル指標macd。このmacdによる売買はテクニカル分析を紹介するサイトで数多く取り上げられています。 前回、macdには騙しが多く使えないのではないか ということを検証してみましたが、この検証における売買方針は、macdゴールデンクロスの翌日寄り買いし、10日後に売るというものでした。この10日という期間にあまり根拠はなく説得力がなかったので、今回はmacdゴールデンクロスの翌日寄り買いデッドクロスが出るまでホールドし、デッドクロスの翌日寄り売りの売買方針ではどうなるかを検証してみたいと思います。 チャートの薄赤の期間買い持ちしているようなかんじですね。この例では、ホールド中上がり続けていることもあるけど、売りシグナルが出るのが遅くて売り時を逃しているのが多い印象です。 macdゴールデンクロスで買いデッドクロスで売るとどうなるかの先行研究 実は、macdゴールデンクロスで買いデッドクロスで売るとどうなるかは、すでに検証されています。 MACDは有効か？【用語解説と統計データ公開】西村剛さんの記事にきれいにまとめられています。 以下、抜粋。 検証対象：日経平均採用銘柄（225種銘柄） 検証期間：2000/01/01～2020/09/30 1銘柄当たりの投資金額：20万円 【買い条件】・MACD（12日、26日）がシグナル（9日）を上抜けした銘柄 ⇒上記を満たした翌日に成行で買い 【売り条件】・MACD（12日、26日）がシグナル（9日）を下抜け ⇒上記を満たした翌日に成行で売り 勝率： 40.13 ％ 勝ち数： 18,561 回 負け数： 27,691 回 引き分け数： 676 回 平均損益（円）： 546 円　　平均損益（率）： 0.27 ％ 平均利益（円）： 13,762 円　　平均利益（率）： 6.88 ％ 平均損失（円）： -8,299 円　　平均損失（率）： -4.15 ％ 合計損益（円）： 25,635,459 円　　合計損益（率）： 12,815.52 ％ 合計利益（円）： 255,437,850 円　　合計利益（率）： 127,723.44 ％ 合計損失（円）： -229,802,391 円　　合計損失（率）： -114,907.92 ％ Profit Factor： 1.112 平均保持日数： 17.29 日 勝率は40％と低いものの、平均損益はプラスとなっています。 損益の推移を確認すると、一部期間を除いて、おおむね右肩上がりの推移となっています。 以上から、MACDは、統計的に有効なテクニカル指標と判断できるでしょう。 これとは別の記事、拙作「騙しに騙されるMACDテクニカル分析」 での検証結果。 検証対象：TOPIX500採用銘柄 検証期間：2015年1月1日～2020年1月20日 結論 ゴールデンクロスの60.3%が騙しであった。 平均損益率： 0.3795% 勝率＝100％ー騙しの割合 と解釈すると、西村剛さんの記事のゴールデンクロス買いの勝率40.13 ％と結果が一致しています。 拙作ではゴールデンクロスの10日後に反対売買、西村剛さんの記事ではデッドクロスで売りの違いはあるものの、平均損益率も近い値となっています。 macdゴールデンクロス売買の追試 あらためて、以下の条件でmacdゴールデンクロスで買いデッドクロスで売りの追試をしてみます。 検証対象：TOPIX500採用銘柄） 検証期間：2015/01/05～2020/01/20 【買い条件】・MACD（12日、26日）がシグナル（9日）を上抜けした銘柄 ⇒上記を満たした翌日に成行で買い 【売り条件】・MACD（12日、26日）がシグナル（9日）を下抜け ⇒上記を満たした翌日に成行で売り 検証データはコチラ。 検証コードは以下。本記事はpythonコードの巧みさでマウントをとるためのものではないので、コードは雑です。（というか、そもそも巧みなコードなど書けない。） test_macd_trading.py # -*- coding: utf-8 -*- """ macd-signalゴールデンクロスで買い、デッドクロスで売り """ import glob import os import math import numpy as np import pandas as pd import matplotlib.pyplot as plt DATA_PATH = "../data/" def calc_macd(df, es, el, sg): macd = pd.DataFrame() macd['ema_s'] = df['Close'].ewm(span=es).mean() macd['ema_l'] = df['Close'].ewm(span=el).mean() macd['macd'] = macd['ema_s'] - macd['ema_l'] macd['signal'] = macd['macd'].ewm(span=sg).mean() macd['diff'] = macd['macd'] - macd['signal'] f_plus = lambda x: x if x > 0 else 0 f_minus = lambda x: x if x < 0 else 0 macd['diff+'] = macd['diff'].map(f_plus) macd['diff-'] = macd['diff'].map(f_minus) return macd def calc_gc(shorts, longs): # ゴールデンクロス検出 gcs = [] short0 = shorts[0] long0 = longs[0] for i in range(len(shorts)): short = shorts[i] long = longs[i] if not math.isnan(short0) and not math.isnan(long0) and not math.isnan(short) and not math.isnan(long): if (short0-long0)<0 and (short-long)>0: gcs.append(shorts.index[i]) short0 = short long0 = long return gcs def calc_gains(df, trades): gains = [] if len(trades)==0: return gains for trade in trades: buy = df['Open'][trade[0]] sell = df['Open'][trade[1]] gain = sell/buy - 1.0 gains.append(gain) return gains def make_macd_trades(df, macd): # macd-signalゴールデンクロスで買い # macd-signalデッドクロスで売り macd_gcs = calc_gc(macd['macd'], macd['signal']) macd_dcs = calc_gc(macd['signal'], macd['macd']) trades = [] hold = False for i in range(len(macd)-1): if macd.index[i] in macd_gcs: # 買いシグナル hold = True buy_date = macd.index[i+1] elif macd.index[i] in macd_dcs: # 売りシグナル if hold: sell_date = macd.index[i+1] trades.append([buy_date, sell_date]) hold = False if hold: sell_date = macd.index[-1] trades.append([buy_date, sell_date]) return trades def walk_around(term=1300): all_gains = [] for i in range(1, 10): csv_files = glob.glob(DATA_PATH+'{}/*.csv'.format(i*1000)) for csv_file in csv_files: filename = os.path.basename(csv_file) code = int(filename[0:4]) code_dir = int(code/1000)*1000 input_path = DATA_PATH + str(code_dir) + "/" + str(code) + ".csv" df = pd.read_csv(input_path, header=None, names=['Date','Open','High','Low','Close','Volume'], encoding='UTF-8') macd = calc_macd(df, 12, 26, 9) trades = make_macd_trades(df, macd) gains = calc_gains(df, trades) all_gains.extend(gains) return all_gains def test1(term=100): all_gains = walk_around(term=term) wins = [ gain for gain in all_gains if gain>0 ] loses = [ gain for gain in all_gains if gain<0 ] draws = [ gain for gain in all_gains if gain==0 ] print("wins: {} {}".format(len(wins), np.mean(wins))) print("loses: {} {}".format(len(loses), np.mean(loses))) print("draws: {} {}".format(len(draws), np.mean(draws))) plt.hist(all_gains, bins=100) mean = np.mean(all_gains) print(mean) def main(): test1(term=1300) if __name__ == '__main__': main() 実行結果 wins: 9601 0.06514704638000195 loses: 14642 -0.03605581596854157 draws: 211 0.0 0.003989021627668806 先行研究の結果フォーマットにあわせると、 勝率： 9601/(9601+14642+211) = 39.26 ％ 勝ち数： 9601 回 負け数： 14642 回 引き分け数： 211 回 平均損益（率）： 0.3989 ％ 平均利益（率）： 6.514 ％ 平均損失（率）： -3.605 ％ Profit Factor： 6.514*9601/(3.605*14642)=62,502.51/52,784.41=1.1841 対象銘柄や検証期間は異なるものの、西村剛さんの記事の結果とほぼ同じになりました。 素朴な疑問 さて、ここからが本題です。 これまでの検証で、確かに、macdゴールデンクロスで買いデッドクロスで売りを繰り返すと、着実に利益が積み上がることが判りました。 ここで１つ素朴な疑問。過去20年では確かにそうなったのだろうけど、そもそも、検証期間2000/01/01～2020/09/30の20年で日経平均は、19002円から27444円に44.4%に上がっている訳ですから、単に日経平均連動投信を買ってその間寝かせているだけで相当の利益が出た訳です。macdゴールデンクロス売買はその利益より大きいのでしょうか？ そこで、なるべく条件を合わせてmacd売買と日経平均連動投信寝かせ売買のパフォーマンスを比較してみます。西村剛さんの記事の結果を使用します。 パフォーマンスを比較のため、総損益率を計算してみます。総損益率は、検証期間のすべての売買の損益率の合計で、全対象銘柄を一定額投資したとき、総損益率×投資額が総損益になります。 macd売買の場合は、毎回一定金額づつ（例えば100万円とか）売買するものとします（実際はこのような売買はできないけど）。macd売買で、対象期間中平均で何銘柄数分をホールドしていたかを計算し、その銘柄数分×100万円分の日経平均連動投信を購入したとして、パフォーマンスを比較します。 macd売買の総損益率 検証期間：2000/01/01～2020/09/30の7305日+270日=7575日間 1回のmacd売買における、1日当たりの損益率：損益率／検証期間＝0.27％／7575 売買回数：勝ち数＋負け数＋引き分け数＝18,561＋27,691＋676＝46,928回 のべ投資日数：平均保持日数×売買回数＝17.29日×46,928回＝811385日 1日当たりの売買回数＝売買回数／検証期間＝46,928回／7575日＝6.195回 (※1) 1日当たりの投資回数：のべ投資日数／検証期間＝811385／7575＝107.11回 (※2) 総損益率＝Σ(1日当たりの損益率)＝(Σ損益率)／7575＝0.27％×46928／7575＝126.705／7575 売買手数料（概算）を考慮すると、 SBI証券の例100万円で往復1280円（0.128％） 総損益率＝Σ(1日当たりの損益率)＝(Σ損益率)／7575 = (0.27％ー0.128％)×46928／7575＝66.638／7575 日経平均売買の総損益率 検証期間：2000/01/01～2020/09/30の7575日間 この間で、日経平均は、19002から27444に44.4%上昇している。 一日当たりの損益率＝0.444／7575 macd売買と投資規模を合わせる。macd売買では、1日平均107.11銘柄を (※2) ホールドしたのだから、 Σ(1日当たりの損益率) = 107.11×0.444／7575＝47.557／7575 売買手数料を考慮しなければ、、macd売買の方が126.705／47.557＝2.66倍パフォーマンスが良い 売買手数料を考慮すると、macd売買の方が66.638／47.557＝1.40倍パフォーマンスが良い 20年間せっせとmacdゴールデンクロス売買を繰り返すと、日経平均連動投信を買って20年間寝かせておくよりも1.40倍成績が良かった！この倍率を多いとみるか少ないとみるかは悩ましいところですが、macd売買は毎日平均6.195回 (※1)、往復だと12.39回 の売買をしなければならないので、その労力を考慮すると、そんなにおいしい訳ではないような。 まとめ 西村剛さんの記事では、平均損益はプラスで、損益の推移はおおむね右肩上がりだから、MACDは統計的に有効なテクニカル指標であると結論付けています。 それは事実ではありますが、検証期間で対象銘柄は指数で44%上がっている訳で、検証期間に日経平均連動投信を寝かせておいた場合とパフォーマンスを条件をあわせて比較すると、1.4倍良いだけにすぎない。または1.4倍も良いと見るべきなのか。 おまけ 最初の図を描いたpythonコードを載せておきます。macd信者様にとってはそれなりに有益かも。私はmacd信者ではなくなりましたので、もはや不要となりました。 コードをコピペしてすぐ動くようにするため、calc_macd(), calc_gc()のコードは先のソースと重複しています。 test_macd_chart.py # -*- coding: utf-8 -*- """ macdチャート """ import math import mplfinance as mpf import numpy as np import pandas as pd DATA_PATH = "../data/" def calc_macd(df, es, el, sg): macd = pd.DataFrame() macd['ema_s'] = df['Close'].ewm(span=es).mean() macd['ema_l'] = df['Close'].ewm(span=el).mean() macd['macd'] = macd['ema_s'] - macd['ema_l'] macd['signal'] = macd['macd'].ewm(span=sg).mean() macd['diff'] = macd['macd'] - macd['signal'] f_plus = lambda x: x if x > 0 else 0 f_minus = lambda x: x if x < 0 else 0 macd['diff+'] = macd['diff'].map(f_plus) macd['diff-'] = macd['diff'].map(f_minus) return macd def calc_gc(shorts, longs): # ゴールデンクロス検出 gcs = [] short0 = shorts[0] long0 = longs[0] for i in range(len(shorts)): short = shorts[i] long = longs[i] if not math.isnan(short0) and not math.isnan(long0) and not math.isnan(short) and not math.isnan(long): if (short0-long0)<0 and (short-long)>0: gcs.append(shorts.index[i]) short0 = short long0 = long return gcs def make_macd_trades(df, macd): # macd-signalゴールデンクロスで買い # macd-signalデッドクロスで売り macd_gcs = calc_gc(macd['macd'], macd['signal']) macd_dcs = calc_gc(macd['signal'], macd['macd']) trades = [] hold = False for i in range(len(macd)-1): if macd.index[i] in macd_gcs: # 買いシグナル hold = True buy_date = macd.index[i+1] elif macd.index[i] in macd_dcs: # 売りシグナル if hold: sell_date = macd.index[i+1] trades.append([buy_date, sell_date]) hold = False if hold: sell_date = macd.index[-1] trades.append([buy_date, sell_date]) return trades def calc_ylim(ser): max1 = ser.max() min1 = ser.min() if max1>0 and min1<0: max2 = max(max1, -min1) return max2, -max2 elif max1>0 and min1>=0: return max1, 0 else: return 0, min1 def draw_macd_chart(df): macd = calc_macd(df, 12, 26, 9) trades = make_macd_trades(df, macd) mh, ml = calc_ylim(macd['macd']) sh, sl = calc_ylim(macd['signal']) hold_line = [np.nan for i in range(len(df))] vlines = [] vcolors = [] for trade in trades: n = len(df['Close'][trade[0]:trade[1]].index) x0 = df.index.get_loc(trade[0]) x1 = x0 + n y0 = df['Open'][trade[0]] y1 = df['Open'][trade[1]] a = (y1-y0)/n b = (x1*y0-x0*y1)/n for j in range(n-1): vlines.append(df['Close'][trade[0]:trade[1]].index[j]) vcolors.append('r') y = a * (x0+j) + b hold_line[x0+j] = y macd['hold_line'] = hold_line vlines1 = dict(vlines=vlines,colors=vcolors,alpha=0.1, linewidths=2.5) add_plot = [ mpf.make_addplot(macd[['hold_line']], color='r'), mpf.make_addplot((macd[['macd', 'signal']]), panel=1, ylim=[min(ml, sl), max(mh, sh)]), mpf.make_addplot((macd['diff+']), type='bar', color='r', panel=1, ylim=[min(ml, sl), max(mh, sh)]), mpf.make_addplot((macd['diff-']), type='bar', color='b', panel=1, ylim=[min(ml, sl), max(mh, sh)]), ] mpf.plot(df, type='candle', addplot=add_plot, vlines=vlines1, volume=False) return def main(): code = 1414 term = 120 code_dir = int(code/1000)*1000 input_path = DATA_PATH + str(code_dir) + "/" + str(code) + ".csv" df = pd.read_csv(input_path, header=None, names=['Date','Open','High','Low','Close','Volume'], encoding='UTF-8') df['Date'] = pd.to_datetime(df['Date']) df = df.set_index("Date") df = df.tail(term) draw_macd_chart(df) if __name__ == '__main__': main()

はじめに データの基礎分析で質的変数の度数分布図（ヒストグラム）を作る必要がある場面になったとき、matplotlib の hist では少し物足りなくなった事情があったので、備忘録として記事に残しておく。 データセットの準備 データセットは質的変数が含まれていればなんでもいいので、iris の場合の例を示す。 #モジュールのインポート from sklearn import datasets import pandas as pd import matplotlib.pyplot as plt #（量的変数は使わないけど）データセットの読み込み iris = datasets.load_iris() df = pd.DataFrame(iris.data) #質的変数のラベルを追加する。 df['species'] = iris.target df.loc[df['species'] == 0, 'species'] = "setosa" df.loc[df['species'] == 1, 'species'] = "versicolor" df.loc[df['species'] == 2, 'species'] = "virginica" print(df) 出力 0 1 2 3 species 0 5.1 3.5 1.4 0.2 setosa 1 4.9 3.0 1.4 0.2 setosa 2 4.7 3.2 1.3 0.2 setosa 3 4.6 3.1 1.5 0.2 setosa 4 5.0 3.6 1.4 0.2 setosa .. ... ... ... ... ... 145 6.7 3.0 5.2 2.3 virginica 146 6.3 2.5 5.0 1.9 virginica 147 6.5 3.0 5.2 2.0 virginica 148 6.2 3.4 5.4 2.3 virginica 149 5.9 3.0 5.1 1.8 virginica histを使ったプロット 質的変数（ラベルデータ等…）でもhistでヒストグラムが作成できる。 fig1 = plt.figure() plt.hist(df["species"]) plt.show() fig1.savefig("iris_speacies_hist.png") しかし、完成したプロットを見て頂ければ分かる通り、setosa, versicolor は種類名がグラフの左側にあるのに対して、virginica はグラフの右側にある。種類の数が少なければ大きな問題には見えないが、種類が増えるとプロットが重なる原因になってしまう。何より不均一なのが気持ち悪い。 度数分布表を作成し、barを使ったプロット 上記の問題を解消するため、度数分布表を別途作成し、barによってヒストグラムを描く方法を取る。 #度数分布を辞書型で作成 col = set(df["species"]) data = df["species"].tolist() dosu = { key: data.count(key) for key in col} ind = range(len(dosu.keys())) #プロット部分 fig2 = plt.figure() plt.xticks(ind, dosu.keys()) plt.bar(ind, dosu.values()) plt.show() fig2.savefig("iris_speacies_bar.png") これによりきちんと種類名の真ん中にラベルが来て、見やすいヒストグラムとなった。 おわりに 質的変数を可視化する際にはちょっとした処理が必要だったので、そのメモを残すために記事を書いた。もし、histを使ってきれいに図がかける方法をご存じの方が居れば、是非コメントお願いします！

概要 数値積分による常微分方程式の数値的な求解を500×500個の初期値それぞれにたいして約6万ステップで行った。 通常のPythonは論外として、NumbaをつかったCPUレベルでの高速化をまず行ったところ、一晩かかった。しかしNumba CUDAを使って10秒ほどで計算が完了した。CUDAに入門してはや1週間ほど、めざましい成果が出たのはこれが初めてだったので記録がてら。 結果 CPUでひと晩かけてじっくりと計算した位相場 GPUで10秒で計算した位相場 結論 無数の初期値にたいして全く同じ演算を行う場合GPUの処理方式がベストマッチしていた。つかうべし。 ソースコード CUDAツールキットとかは一通り入っている前提 方程式のリミットサイクルを固定刻み幅の4段Runge-Kutta法でもとめる リミットサイクル上の点の値(x,y)が求まる。適当な点をえらび、またその点と(ほぼ)同じ点に帰ってくるまでの経過時間を位相とする。刻み幅は変えていないので点の個数が経過時間とそのまま対応する。 リミットサイクルを含んだ平面の格子点からn周期ぶんRunge-Kutta法をまわす それぞれの格子点がたどりついた点と、2で求めたリミットサイクル上の点とをコサイン類似度で比較 もっとも類似度が高い点と同じ位相をその格子点の位相として、色で対応付ける。 理論についてはよくわからん場合詳しくはこちら https://www.notion.so/4-1-9bead7cd44fb443f8da408e9bc11c0bc 3-4のところの処理 cpu.py import numpy as np import matplotlib.pyplot as plt from numba import jit r=500 @jit('f8(f8, f8)') def fhn_xdot(x, y): return x*(x+0.1)*(1-x)-y @jit('f8(f8, f8)') def fhn_ydot(x, y): return (x-0.5*y)/100 #Tは周期(今回は12648×0.01秒), lcvecotrはリミットサイクル上の点の座標がはいったndarray(T×2),lcnormはその点の原点からの距離(T×1) for i in range(r): x = np.linspace(-0.5, 1, r)[i] for j in range(r): y = np.linspace(-0.1, 0.3,r)[j] x2, y2, _ = rk4th_2d(fhn_xdot, fhn_ydot, x, y, T)#適宜積分する。 col = np.argmax(lcvector@np.array([x2[-1], y2[-1]])/(lcnorm*np.sqrt(x2[-1]**2+y2[-1]**2)))/T#コサイン類似度 pf[i][j] = col forは入れ子にしなくともいいのかもしれない。詳しい人教えて下さい gpu.py from numba import cuda #デバイス関数を定義しておく @cuda.jit('f8(f8, f8)', device=True) def fhn_xdot(x, y): return x*(x+0.1)*(1-x)-y @cuda.jit('f8(f8, f8)', device=True) def fhn_ydot(x, y): return (x-0.5*y)/100 @cuda.jit def rk4_gpu_fhn(pf, mx, my): #x→何番目の振動子か?、y→状態量はどれか? k, j = cuda.grid(2) x, y = mx[k], my[j] dt = 0.01 for _ in range(12648*5): k_1x = fhn_xdot(x,y) k_1y = fhn_ydot(x,y) k_2x = fhn_xdot(x+k_1x/2*dt,y+k_1y/2*dt) k_2y = fhn_ydot(x+k_1x/2*dt,y+k_1y/2*dt) k_3x = fhn_xdot(x+k_2x/2*dt,y+k_2y/2*dt) k_3y = fhn_ydot(x+k_2x/2*dt,y+k_2y/2*dt) k_4x = fhn_xdot(x+k_3x*dt,y+k_3y*dt) k_4y = fhn_ydot(x+k_3x*dt,y+k_3y*dt) x+=((k_1x+2*k_2x+2*k_3x+k_4x)/6*dt) y+=((k_1y+2*k_2y+2*k_3y+k_4y)/6*dt) pf[k, j] = x, y mx1 = np.linspace(-0.5, 1, r) my1 = np.linspace(-0.1, 0.3, r) griddim=r blockdim=1, r pf = np.empty((r, r,2)) rk4_gpu_fhn[griddim, blockdim](pf, mx1, my1) #いろつける pfc = np.empty((r, r)) for i in range(r): for j in range(r): col = np.argmax(lcvector@pf[i][j]/(lcnorm*np.linalg.norm(pf[i][j])))/T#コサイン類似度 pfc[i][j] = col cpuでは、添字をつかってじゅんぐりに積分していくようなことしかできなかったのだけど、CUDAのプログラミングモデルでは｢任意の添字にたいする一般的な処理｣という形で書くことができて、blockdim×griddimぶん同時にその処理を走らせることができる。｢任意の添字｣は｢cuda.grid(2)｣で取得している。 可視化も結構めんどくさかったので参考で書いておく visualize.py cm = plt.get_cmap("jet") fig=plt.figure(dpi=200) ax = fig.add_subplot(1, 1, 1) ax.plot(px1, py1) ax.set_xlim(-0.5, 1) ax.set_ylim(-0.1,0.3) ax.set_xlabel("u") ax.set_ylabel("v") xlim = ax.get_xlim() ylim = ax.get_ylim() mp=ax.imshow(pfc.T*2*np.pi, extent=[*xlim, *ylim], aspect='auto',vmin=0, vmax=2*np.pi, origin="lower", cmap=cm) cbar = fig.colorbar(mp, ticks=[0, np.pi, 2*np.pi]) cbar.ax.set_yticklabels([r'0', r'$\pi$', r'$2\pi$'])

1.はじめに pythonの学習を進める中でWebアクセスやスクレイピングなんてものにもチュートリアル程度に触れたのですが、 如何せん、しばらくは業務で使う機会もなさそうなので備忘録として記事を残すことにしました。 2.アプリを作成する アプリの作成に沿って使用したライブラリの使用例を挙げていきます。 今日の天気を気象庁のWebAPIから取得して表示してくれるアプリを作成します。 3.WebAPIにアクセス(urllib) まずは気象庁のサイトから本日の天気の情報を取得します。 今回はpythonに標準機能として用意されているurllibというライブラリを使用してWebAPIにアクセスを試みます。 下記のコードを実行します。 "urllib.request.urlopen(url)"で気象庁のWebApiにアクセス "response.read()"で内容を読出 #!/usr/bin/env python import urllib.request url = "https://www.jma.go.jp/bosai/forecast/data/forecast/170000.json" response = urllib.request.urlopen(url) print(response.read()) 実行結果として以下のレスポンスが得られました。 b'[{"publishingOffice":"\xe9\x87\x91\xe6\xb2\xa2\xe5\x9c\xb0\xe6\x96\xb9\xe6\xb0\x97\xe8\xb1\xa1\xe5\x8f\xb0","reportDatetime":"2021-06-21T05:00:00+09:00","timeSeries":[{"timeDefines":["2021-06-21T05:00:00+09:00","2021-06-22T00:00:00+09:00"],"areas":[{"area":{"name":"\xe5\x8a\xa0\xe8\xb3\x80","code":"170010"},"weatherCodes":["100","212"],"weathers":["\xe6\x99\xb4\xe3\x82\x8c","\xe3\x81\x8f\xe3\x82\x82\xe3\x82\x8a\xe3\x80\x80\xe5\xa4\x95\xe6\x96\xb9\xe3\x80\x80\xe3\x81\x8b\xe3\x82\x89\xe3\x80\x80\xe5\xa4\x9c\xe3\x81\xae\xe3\x81\xaf\xe3\x81\x98\xe3\x82\x81\xe9\xa0\x83\xe3\x80\x80\xe9\x9b\xa8\xe3\x80\x80\xe6\x89\x80\xe3\x81\xab\xe3\x82\x88\xe3\x82\x8a\xe3\x80\x80\xe6\x98\xbc\xe9\x81\x8e\xe3\x81\x8e\xe3\x80\x80\xe3\x81\x8b\xe3\x82\x89\xe3\x80\x80\xe5\xa4\x9c\xe3\x81\xae\xe3\x81\xaf\xe3\x81\x98\xe3\x82\x81\xe9\xa0\x83\xe3\x80\x80\xe9\x9b\xb7\xe3\x80\x80\xe3\x82\x92\xe4\xbc\xb4\xe3\x81\x86"],"winds":["\xe6\x9d\xb1\xe3\x81\xae\xe9\xa2\xa8\xe3\x80\x80\xe6\x97\xa5\xe4\xb8\xad\xe3\x80\x80\xe5\x8c\x97\xe3\x81\xae\xe9\xa2\xa8","\xe8\xa5\xbf\xe3\x81\xae\xe9\xa2\xa8"],"waves":["\xef\xbc\x91\xef\xbc\x8e以下略... 4.データを成型する(BeautifulSoup4) 上記で気象庁の公開している情報の取得は出来ましたが、 型がバイナリ型であったり、改行されていなかったりと正直データが見づらいです。 そこでBeautifulSoupというライブラリを使用してデータを成型します。 まず下記コマンドでインストールします。 pip install beautifulfulsoup4 インストールが完了したら下記のコードを実行します。 上記のコードからの変更点は "from bs4 import BeautifulSoup"でBeautifulSoupをインポート "soup = BeautifulSoup(response, 'html.parser')"でレスポンスを成型。 "html.parser"は構文を解析するツールを指定しています。 #!/usr/bin/env python import urllib.request from bs4 import BeautifulSoup # New!! url = "https://www.jma.go.jp/bosai/forecast/data/forecast/170000.json" response = urllib.request.urlopen(url) # print(response.read()) # New!! BeautifulSoupで変換 soup = BeautifulSoup(response, 'html.parser') print(soup) 実行結果は以下となります。 [{"publishingOffice":"金沢地方気象台","reportDatetime":"2021-06-21T05:00:00+09:00","timeSeries":[{"timeDefines":["2021-06-21T05:00:00+09:00","2021-06-22T00:00:00+09:00"],"areas":[{"area":{"name":"加賀","code":"170010"},"weatherCodes":["100","212"],"weathers":["晴れ","くもり　夕方　から　夜のはじめ頃　雨　所により　昼過ぎ　から　夜のはじめ頃　雷　を伴う"],"winds":["東の風　日中　 北の風","西の風"],"waves":["１．５メートル　後　１メートル　うねり　を伴う","１メートル"]},{"area":{"name":"能登","code":"170020"},"weatherCodes":["211","200"],"weathers":["くもり　昼前　から 　晴れ","くもり　所により　昼前　から　夕方　雨"],"winds":["南の風　日中　北東の風","南の風　日中　西の風"],"waves":["１．５メートル　後　１メートル　うねり　を伴う","１メートル"]}]},{"timeDefines":["2021-06-21T06:00:00+09:00","2021-06-21T12:00:00+09:00","2021-06-21T18:00:00+09:00","2021-06-22T00:00:00+09:00","2021-06-22T06:00:00+09:00","2021-06-22T12:00:00+09:00","2021-06-22T18:00:00+09:00"],"areas":[{"area":{"name":"加賀","code":"170010"},"pops":["0","0","0","10","10","50","50"]},{"area":{"name":"能登","code":"170020"},"pops":["0","0","0","10","10","20","20"]}]},{"timeDefines":["2021-06-21T09:00:00+09:00","2021-06-21T00:00:00+09:00","2021-06-22T00:00:00+09:00","2021-06-22T09:00:00+09:00"],"areas":[{"area":{"name":"金沢","code":"56227"},"temps":["30","30","20","26"]},{"area":{"name":"輪島","code":"56052"},"temps":["26","26","18","25"]}]}]},{"publishingOffice":"金沢地方気象台","reportDatetime":"2021-06-20T17:00:00+09:00","timeSeries":[{"timeDefines":["2021-06-21T00:00:00+09:00","2021-06-22T00:00:00+09:00","2021-06-23T00:00:00+09:00","2021-06-24T00:00:00+09:00","2021-06-25T00:00:00+09:00","2021-06-26T00:00:00+09:00","2021-06-27T00:00:00+09:00"],"areas":[{"area":{"name":"石川県","code":"170000"},"weatherCodes":["101","212","202","200","200","200","200"],"pops":["","50","50","40","30","40","40"],"reliabilities":["","","C","B","A","B","B"]}]},{"timeDefines":["2021-06-21T00:00:00+09:00","2021-06-22T00:00:00+09:00","2021-06-23T00:00:00+09:00","2021-06-24T00:00:00+09:00","2021-06-25T00:00:00+09:00","2021-06-26T00:00:00+09:00","2021-06-27T00:00:00+09:00"],"areas":[{"area":{"name":"金沢","code":"56227"},"tempsMin":["","20","19","19","20","20","20"],"tempsMinUpper":["","21","20","21","21","23","23"],"tempsMinLower":["","18","17","17","18","18","18"],"tempsMax":["","26","26","26","27","27","27"],"tempsMaxUpper":["","28","28","29","29","30","30"],"tempsMaxLower":["","25","24","24","24","23","23"]}]}],"tempAverage":{"areas":[{"area":{"name":"金沢","code":"56227"},"min":"19.7","max":"26.4"}]},"precipAverage":{"areas":[{"area":{"name":"金沢","code":"56227"},"min":"31.4","max":"70.1"}]}}] 5.更にデータを成型する(json) だいぶ正体が見えてきましたが、 jsonクラスを使用して更に人間側に分かりやすい形式に変換します。 #!/usr/bin/env python import urllib.request from bs4 import BeautifulSoup import json # New!! url = "https://www.jma.go.jp/bosai/forecast/data/forecast/170000.json" response = urllib.request.urlopen(url) # print(response.read()) soup = BeautifulSoup(response, 'html.parser') # print(soup) # New!!　見やすいJSON形式に変換 data = json.loads("{}".format(soup, 'utf-8')) data = json.dumps(data, indent=2, ensure_ascii=False) print(data) 実行結果は以下となります。 [ { "publishingOffice": "金沢地方気象台", "reportDatetime": "2021-06-21T11:00:00+09:00", "timeSeries": [ { "timeDefines": [ "2021-06-21T11:00:00+09:00", "2021-06-22T00:00:00+09:00", "2021-06-23T00:00:00+09:00" ], "areas": [ { "area": { "name": "加賀", "code": "170010" }, "weatherCodes": [ "100", "212", "202" ], "weathers": [ "晴れ", "くもり　夕方　から　夜のはじめ頃　雨　所により　昼過ぎ　から　夜のはじめ頃　雷　を伴う", "くもり　一時　雨" ], "winds": [ "北の風　後　東の風", "西の風", "南の風　後　北西の風" ], "waves": [ "１．５メートル　後　１メートル　うねり　を伴う", "１メートル", "０．５メートル" ] }, { "area": { "name": "能登", "code": "170020" }, "weatherCodes": [ "100", "200", "202" ], "weathers": [ "晴れ", "くもり　所により　昼前　から　夕方　雨", "くもり　一時　雨" ], "winds": [ "北東の風　後　南の風", "南の風　日中　西の風", "南の風　後　北の風" ], "waves": [ "１．５メートル　後　１メートル　うねり　を伴う", "１メートル", "０．５メートル" ] } ] }, { "timeDefines": [ "2021-06-21T12:00:00+09:00", "2021-06-21T18:00:00+09:00", "2021-06-22T00:00:00+09:00", "2021-06-22T06:00:00+09:00", "2021-06-22T12:00:00+09:00", "2021-06-22T18:00:00+09:00" ], "areas": [ { "area": { "name": "加賀", "code": "170010" }, "pops": [ "10", "0", "10", "10", "50", "50" ] }, { "area": { "name": "能登", "code": "170020" }, "pops": [ "0", "0", "10", "10", "20", "20" ] } ] }, { "timeDefines": [ "2021-06-21T09:00:00+09:00", "2021-06-21T00:00:00+09:00", "2021-06-22T00:00:00+09:00", "2021-06-22T09:00:00+09:00" ], "areas": [ { "area": { "name": "金沢", "code": "56227" }, "temps": [ "30", "30", "20", "26" ] }, { "area": { "name": "輪島", "code": "56052" }, "temps": [ "26", "26", "18", "25" ] } ] } ] }, { "publishingOffice": "金沢地方気象台", "reportDatetime": "2021-06-21T11:00:00+09:00", "timeSeries": [ { "timeDefines": [ "2021-06-22T00:00:00+09:00", "2021-06-23T00:00:00+09:00", "2021-06-24T00:00:00+09:00", "2021-06-25T00:00:00+09:00", "2021-06-26T00:00:00+09:00", "2021-06-27T00:00:00+09:00", "2021-06-28T00:00:00+09:00" ], "areas": [ { "area": { "name": "石川県", "code": "170000" }, "weatherCodes": [ "212", "202", "201", "201", "200", "202", "202" ], "pops": [ "", "50", "30", "30", "30", "50", "50" ], "reliabilities": [ "", "", "A", "A", "A", "C", "C" ] } ] }, { "timeDefines": [ "2021-06-22T00:00:00+09:00", "2021-06-23T00:00:00+09:00", "2021-06-24T00:00:00+09:00", "2021-06-25T00:00:00+09:00", "2021-06-26T00:00:00+09:00", "2021-06-27T00:00:00+09:00", "2021-06-28T00:00:00+09:00" ], "areas": [ { "area": { "name": "金沢", "code": "56227" }, "tempsMin": [ "", "19", "19", "19", "20", "21", "20" ], "tempsMinUpper": [ "", "21", "20", "21", "22", "22", "23" ], "tempsMinLower": [ "", "17", "17", "18", "18", "18", "18" ], "tempsMax": [ "", "25", "26", "27", "28", "27", "25" ], "tempsMaxUpper": [ "", "28", "27", "29", "30", "30", "29" ], "tempsMaxLower": [ "", "23", "23", "24", "25", "24", "22" ] } ] } ], "tempAverage": { "areas": [ { "area": { "name": "金沢", "code": "56227" }, "min": "19.8", "max": "26.5" } ] }, "precipAverage": { "areas": [ { "area": { "name": "金沢", "code": "56227" }, "min": "33.5", "max": "73.0" } ] } } ] 6.必要なデータを表示する 上記の取得値から"weathers"の最初のデータが当日の値であることが分かります。 JSONの構造が見えたので下記コードで天気を出力させます。 #!/usr/bin/env python import urllib.request from bs4 import BeautifulSoup import json url = "https://www.jma.go.jp/bosai/forecast/data/forecast/170000.json" response = urllib.request.urlopen(url) # print(response.read()) soup = BeautifulSoup(response, 'html.parser') # print(soup) data = json.loads("{}".format(soup, 'utf-8')) # data = json.dumps(data, indent=2, ensure_ascii=False) # print(data) # New!! 天気を表示 print("今日の天気は{}。".format(data[0]["timeSeries"][0]["areas"][0]["weathers"][0])) 実行結果は以下となります。 今日の天気は晴れ。 7.さいごに 今回使用した気象庁のWebAPI以外にも世の中には有料のものから無料のものまで様々な情報がWebAPIとして公開されているみたいです。 利用できるものは何でも利用して高機能なプログラムを素早く開発できると素敵だなと思いました。 8.参照 気象庁HP(https://www.jma.go.jp/jma/index.html)

以下の記事も参考にして、SPSS Modelerでタイムスタンプ関連のclem関数を使って、よく行う特徴量抽出を行ってみます。そしてその処理をPythonのpandasで書き換えてみます。 日付型については以下の記事で紹介しています。 https://qiita.com/kawada2017/items/c41f443e5246d331bea1 代表的なものとして以下の加工処理を行っていきます。 タイムスタンプデータの読み込み 文字列や数値からのタイムスタンプデータ生成 現在タイムスタンプの取得 タイムスタンプの年、月、日、時、分、秒の分解 タイムスタンプの差の計算 タイムスタンプの大小比較 タイムスタンプの加算減算 以下の、タイムスタンプ(TS),電力(Power),温度(Temperature),圧力(Pressure),起動時間(Uptime),状態(Status),エラーコード(Outcome)が記録された時系列のセンターデータを使います。毎分のセンサーの値が記録されています。今回はこの中のTSのみを使います。 1m．① 文字列や数値からのタイムスタンプデータ生成　Modeler版 文字列や数値からタイムスタンプ型データを生成してみます。タイムスタンプ演算をするために文字列データから型変換をすることはよくあります。また年、月、日、時、分、秒が別列に分かれて保存されている場合に、それらの列からタイムスタンプ型のデータを作成することもよくあります。 まずModelerのタイムスタンプや時間のフォーマットを確認しておきます。ストリームのオプションの「日付/時刻」の中にあります。YYYY-MM-DDとHH:MM:SSがデフォルトです。 CSVデータを読み込んだ際に、デフォルトでは「自動的に日付と時間を認識します」のオプションが有効になっており、ストリームオプションで指定されていたYYYY-MM-DD HH:MM:SSで解釈できる文字列データはタイムスタンプデータとして読み込みます。 1p．① 文字列や数値からのタイムスタンプデータ生成　pandas版 pythonには以下のような複数のTimestamp型があります。 * datetime.datetime * numpy.datetime64 * pandas._libs.tslibs.timestamps.Timestamp データ型確認 #datetime.datetime dtdt=datetime.datetime(2021, 6, 10, 3, 11, 0) print(dtdt) print(type(dtdt)) #numpy.datetime64 npdt=np.datetime64('2021-06-10 03:11:00') print(npdt) print(type(npdt)) #pandas._libs.tslibs.timestamps.Timestamp pddt=pd.to_datetime('2021-06-10 03:11:00') print(pddt) print(type(pddt)) 表示上はnumpy.datetime64は日付と時間の間に「T」が入っていました。 結果 2021-06-10 03:11:00 <class 'datetime.datetime'> 2021-06-10T03:11:00 <class 'numpy.datetime64'> 2021-06-10 03:11:00 <class 'pandas._libs.tslibs.timestamps.Timestamp'> この記事ではpandasのデータフレームを使うので、pandasのTimestamp型を使用していきます。 しかし、datetime.datetimeとnumpy.datetime64をpandasの列にいれると自動的にpandas._libs.tslibs.timestamps.Timestampに変換されていました。 dfdt.dtypesでみるといずれもdatetime64[ns]として表示されます。互換性があるようです。 dfdt =pd.DataFrame({'datetime':[dtdt], 'npdatetime':[npdt], 'pddatetime':[pddt]}) print(dfdt.dtypes) print(type(dfdt['datetime'][0])) print(type(dfdt['npdatetime'][0])) print(type(dfdt['pddatetime'][0])) 結果 datetime datetime64[ns] npdatetime datetime64[ns] pddatetime datetime64[ns] dtype: object <class 'pandas._libs.tslibs.timestamps.Timestamp'> <class 'pandas._libs.tslibs.timestamps.Timestamp'> <class 'pandas._libs.tslibs.timestamps.Timestamp'> csvファイルからpandasに読むこむときにはデフォルトでは文字列型で読み込まれてしまいますが、パースしてタイムスタンプ型として読み込むこともできます。 parse_datesに列を指定すると、その列がタイムスタンプとして解釈できればタイムスタンプ型データとして読み込まれます。 数値型データからタイムスタンプ型データ生成 #現在のタイムスタンプ df = pd.read_csv('COND2nts.csv', parse_dates=['TS']) print(df.dtypes) print(type(df['TS'][0]) 以下ではTS列が、文字列としてではなく、datetime64[ns]型で読み込まれています。 結果 TS datetime64[ns] Power int64 Temperature int64 Pressure int64 Uptime int64 Status int64 Outcome int64 dtype: object <class 'pandas._libs.tslibs.timestamps.Timestamp'> 参考 2m．② 文字列や数値からのタイムスタンプデータ生成　Modeler版 フィールド作成ノードで文字列からタイムスタンプ型データを生成してみます。datetime_dateという関数でストリームオプションで指定されていたYYYY-MM-DD形式で文字列を指定し、datetime_timeという関数でHH:MM:SS形式で文字列を指定します。そしてその日付型データと時間型データをdatetime_timestampでタイムスタンプ型として統合します。 datetime_timestamp(datetime_date('2021-04-01'),datetime_time('03:11:00')) 続いて、数値型データからタイムスタンプ型データを生成してみます。やはりdatetime_dateという関数を使いますが、年、月、日、時、分、秒の数値をカンマ区切りで指定します。 datetime_timestamp(2021,4,1,3,11,00) データ型ノードをつかって確認してみると、アイコンがカレンダーと時計になり、尺度が連続型になっています。正しくタイムスタンプ型が生成されていることが確認できました。 2p．② 文字列や数値からのタイムスタンプデータ生成　pandas版 あらためて文字列からタイムスタンプ型データを生成してみます。pd.to_datetimeという関数でYYYY-MM-DD HH:MM:SS形式で文字列を指定します。YYYY-MM-DD HH:MM:SS形式以外のフォーマットの場合はformatstrのオプションでフォーマットを指定します。 文字列からタイムスタンプ型データ生成 df1['REF_TS']=pd.to_datetime('2021-04-01 03:11:00') 続いて、数値型データからタイムスタンプ型データを生成してみます。いったんdatetime.dateで年、月、日、時、分、秒を数値で指定してdatetime.datetime型で生成します。pd.to_datetimeで変換せずとも自動的にpandas._libs.tslibs.timestamps.Timestampに変換されます。 数値からタイムスタンプ型データ生成 df1['REF_TS2']=datetime.datetime(2021,4,1,3,11,0) print(df1.dtypes) 結果 TS datetime64[ns] Power int64 Temperature int64 Pressure int64 Uptime int64 Status int64 Outcome int64 REF_TS datetime64[ns] REF_TS2 datetime64[ns] dtype: object 参考 3m．③現在のタイムスタンプの取得　Modeler版 「現在のタイムスタンプ」を取得してみます。現在のタイムスタンプはデータ生成日時の記録などによく使います。 フィールド作成ノードで「datetime_now」という関数を設定します。 現在のタイムスタンプが入ります。 3p．③ 現在のタイムスタンプの取得　pandas版 datetime.date.now()で現在のタイムスタンプを取得します。 なお、np.datetime64('now')でも取得はできましたが、マニュアルにこの記述方法がなかったので、datetimeを使う方法の方が堅いかと思います。 数値型データからタイムスタンプ型データ生成 #現在のタイムスタンプ df2['TODAYTS']=datetime.datetime.now() #以下でも動作するがマニュアルに見当たらなかった。 #df2['TODAYTS']=np.datetime64('now') print(df2.dtypes) 結果 TS datetime64[ns] Power int64 Temperature int64 Pressure int64 Uptime int64 Status int64 Outcome int64 TODAYTS datetime64[ns] dtype: object 参考 4m．④ タイムスタンプの年、月、日、時、分、秒の分解　Modeler版 タイムスタンプから年、月、日、時、分、秒を取得してみます。時間毎のグループ化集計などのために抜き出すことがよくあります。 フィールド作成ノードでdatetime_year,datetime_month、datetime_day、datetime_hour、datetime_minute、datetime_secondという関数をつかって抜き出すことができます。ここではdatetime_minuteで分を抜き出してみます。 4p．④ タイムスタンプの年、月、日、時、分、秒の分解　pandas版 dtアクセサをつかって年year、月month、日day、時hour、分minute、秒secondの分解が可能です。dtアクセサはdayofyear（年初からの日数）やquarter(四半期)なども抜き出すことができます。 ここではdt.minuteで月を抜き出します。結果はintで返ります。 数値型データからタイムスタンプ型データ生成 df3['MIN']=df3['TS'].dt.minute 結果 TS datetime64[ns] Power int64 Temperature int64 Pressure int64 Uptime int64 Status int64 Outcome int64 MIN int64 dtype: object *参考 5m．⑤ タイムスタンプの差の計算　Modeler版 タイムスタンプ同士の差を計算してみます。あるシグナルとシグナルの時間間隔や傾きをとる場合など、時間差は重要な特徴量になりえます。 フィールド作成ノードでdatetime_in_secondsで秒に変換してから差を取ります。 以下の例だとタイムスタンプ(TS)- 基準タイムスタンプ(REF_TS)の秒数を返します。基準タイムスタンプ(REF_TS)＞タイムスタンプ(TS)であればマイナスの値が返ります。 datetime_in_seconds(TS) - datetime_in_seconds(REF_TS) さらに、分の単位で差をだしたい場合は60秒で割ることで、 ( datetime_in_seconds(TS) - datetime_in_seconds(REF_TS))/60 時間の単位で差をだしたい場合は60秒*60分で割ることで求めることができます。 ( datetime_in_seconds(TS) - datetime_in_seconds(REF_TS))/(60*60) 2021-04-01 03:11:00という基準日とタイムスタンプ(TS)の差がTS_DIFF秒、TS_DIFF分、TS_DIFF、に入っています。 5p．⑤タイムスタンプの差の計算　pandas版 まず、以下のように単純にタイムスタンプの入った列の引算を行います。 df4["TS_DIFF_TD"]=(df4["TS"]-df4["REF_TS"]) これは整数ではなく、pandas._libs.tslibs.timedeltas.Timedeltaという特殊な型（df5.dtypesの結果だとtimedelta64[ns]）で値を返します。 ややこしいのは、Timedeltaはマイナスの時間を表さない仕様であることです。-1分の場合、-1 days +23:59:00というようなわかりにくい値が返ります（合計すると1分）。 これを.dt.secondsで秒数に変換してしまうと23:59:00を秒変換するので、86340秒（24*60*60-60）という期待しない値が返ってしまいます。 単なる秒数がほしいので、dtアクセサとをtotal_seconds()メソッドつかって.dt.total_seconds()で秒数に変換します。プロパティではなくメソッドなので（）が付くことに注意してください。 タイムスタンプの差の計算 df4["TS_DIFF_TD"]=(df4["TS"]-df4["REF_TS"]) print(type((df4["TS_DIFF_TD"])[0])) #以下は見た目も分かりやすい df4["TS_DIFF秒"]=(df4["TS"]-df4["REF_TS"]).dt.total_seconds().astype('int') df4["TS_DIFF分"]=(df4["TS"]-df4["REF_TS"]).dt.total_seconds()/60 df4["TS_DIFF時"]=(df4["TS"]-df4["REF_TS"]).dt.total_seconds()/(60*60) #以下だとマイナスの時に間違った値が戻る #df4["TS_DIFF"]=(df4["TS"]-df4["REF_TS"]).dt.seconds #以下でも計算はできる #df4["TS_DIFF秒"]=((df4["TS"]-df4["REF_TS"]) /np.timedelta64(1,'s')).astype('int') #df4["TS_DIFF分"]=((df4["TS"]-df4["REF_TS"]) /np.timedelta64(1,'m')) #df4["TS_DIFF時"]=((df4["TS"]-df4["REF_TS"]) /np.timedelta64(1,'h')) print(df4.dtypes) 結果 <class 'pandas._libs.tslibs.timedeltas.Timedelta'> TS datetime64[ns] Power int64 Temperature int64 Pressure int64 Uptime int64 ... REF_TS datetime64[ns] TS_DIFF_TD timedelta64[ns] TS_DIFF秒 int32 TS_DIFF分 float64 TS_DIFF時 float64 TS_DIFF_TDとTS_DIFF秒の結果は同じですが、TS_DIFF_TDにはtimedelta64[ns]型で格納され、TS_DIFF秒は整数型になっています。そしてマイナス値になるときにはTS_DIFF秒の方が理解しやすい値として返っています。例では-60秒です。 timedelta（Python）のマイナス値表現の謎 - Qiita https://qiita.com/__Attsun__/items/836bb3954ddcc6413dc4 pandas.Series.dt.total_seconds https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.dt.total_seconds.html# 6m．⑥タイムスタンプの大小比較　Modeler版 タイムスタンプの大小比較をしてみます。ある作業工程以降かというようなフラグの特徴量が作れます。 フィールド作成ノードで派生：フラグ型に設定し、条件に以下を入力します。基準日より後の購入日ならフラグを立てるという意味です。 TS>REF_TS 2021-04-01 03:11:00 という基準タイムスタンプ以降の場合にAFTER_REFにTが立っています。 6p．⑥ タイムスタンプの大小比較　pandas版 以下のように単純にタイムスタンプの比較式を代入します。結果はbool型で返ります。 タイムスタンプの差の計算 df5["AFTER_REF"]=(df5["TS"]>df5["REF_TS"]) print(df5.dtypes) 結果 TS datetime64[ns] Power int64 Temperature int64 Pressure int64 Uptime int64 Status int64 Outcome int64 REF_TS datetime64[ns] AFTER_REF bool dtype: object 2021-04-01 03:11:00という基準タイムスタンプ以降の売上の場合にAFTER_REFにTrueが立っています。 7m．⑦ タイムスタンプの加算減算　Modeler版 タイムスタンプに加算減算を計算してみます。起動からから3時間後というような計算です。 フィールド作成ノードでdatetime_in_secondsとdatetime_dateという関数をつかって、タイムスタンプの3時間後を生成してみます。 まずdatetime_in_seconds(TS)でタイムスタンプを基準日からの秒数に変換します。 基準日はデフォルトだと1900-01-01です。 次にこの秒数に3時間分の秒数を足します。3時間*60分*60秒です。 そしてこの秒数をdatetime_dateでタイムスタンプ型に変換しなおします。 つまり式としては以下になります。 datetime_timestamp(datetime_in_seconds(TS)+3*60*60) タイムスタンプ(TS)の3時間後のタイムスタンプが3H_AFTERに生成されました。 7p．⑦ タイムスタンプの加算減算　pandas版 datetime.timedeltaをつかって量と単位を指定して加算減算をします。以下の例は時間hoursを3個分足す、つまり3時間分足すという意味になります。hoursの部分を変えると分minutes、秒secondsなど別の単位でも加算減算が行えます。 タイムスタンプの加算減算 df5['3H_AFTER']=df5["TS"]+ datetime.timedelta(hours=3) #以下でも実現は可能 #df5['3H_AFTER']=df5["TS"]+ np.timedelta64(3,'h') print(df5.dtypes) 結果 TS datetime64[ns] Power int64 Temperature int64 Pressure int64 Uptime int64 Status int64 Outcome int64 3H_AFTER datetime64[ns] dtype: object タイムスタンプ(TS)の3時間後のタイムスタンプが3H_AFTERに生成されました。 8. サンプル サンプルは以下に置きました。 ストリーム https://github.com/hkwd/200611Modeler2Python/raw/master/datetime/ts.str notebook https://github.com/hkwd/200611Modeler2Python/blob/master/datetime/ts.ipynb データ https://raw.githubusercontent.com/hkwd/200611Modeler2Python/master/data/Cond4n_e.csv ■テスト環境 Modeler 18.2.2 Windows 10 64bit Python 3.8.5 pandas 1.0.5 numpy 1.19.2

Spectral Python Spectral Python(SPy)とは純度100%のpythonのモジュールで、ハイパースペクトル(HS)データの処理ができる。HSデータの読み込み、表示、計算そして分類まで一括して行える。MITライセンスのため、自由に使うことができる。 Introduction このユーザーガイドでは、下記の論文のデータを参照しチュートリアルを進める。またチュートリアルでは利便性のためfrom spectral import *としているが、モジュール内ネームスペースのコンテンツをすべてインポートすることは非推奨である。データはここから。 File Name Description 92AV3C.lan ERDAS/LanフォーマットのHS画像。AVRISでの撮影で得られた145*145ピクセル、220バンドのデータ。 92AV3GT.GIS 上のデータの土地利用図。 92AV3C.spc 上のデータのAVIRISフォーマットのバンドカリブレーションファイル Landgrebe, D. Multispectral data analysis: A signal theory perspective. School of Electr. Comput. Eng., Purdue Univ., West Lafayette, IN (1998). Installing Spy 基本的にはPythonの組み込みのものとNumpyがあればHSデータの処理は可能である。グラフィカルな機能を使用して解析を行いたい場合には別途モジュールが必要となる。下表にそれぞれの機能の要件を示す。 Dependency Notes Python 2.6+ or 3.3+ (1) Numpy Required Pillow or Python Imaging Library(PIL) Required if displaying or saving images wxPython (2) matplotlib Required if rendering raster displays or spectral plots IPython Required for interactive, non-blocking GUI windows PyOpenGL (2) Notes: (1): view_cube、view_nd以外はPython 3.3+で動作 (2): view_cube、view_ndを呼び出す際に必要 pip でのインストールは下記のコマンドで。 その他の方法はここから pip install spectral 基本的にHSデータは容量が大きいため、特定のディレクトリに保存されている。SPyでデータを開こうとするたびに絶対パスを入力しなくて済むようにSPECTRAL_DATAという環境変数を設定できる。SPyはこの環境変数に記されているディレクトリからイメージファイルを探す。またディレクトリパスのリストは:区切りで記載されなければならない。 Reading HSI Data Files 通常、SPyではHSイメージファイルにアクセスし開くためにimage関数を使用する。この関数はSpyFileオブジェクトのインスタンスを返す。 The SpyFile Interface SpyFileはHSデータを読み込むオブジェクトの基本的なクラスである。SpyFileオブジェクトの作成時、対応するファイルからデータを読み込むインターフェースが与えられる。画像を開いたとき、返される実際のオブジェクトはSpyFileのサブクラス(BipFile, BilFile, BsqFile)で、このサブクラスはデータの格納方法(interleaves)にそれぞれ対応している。 下記のコードではサンプル画像を開いている。 In [1]: from spectral import * In [2]: img = open_image('92AV3C.lan') In [3]: img.class Out[3]: spectral.io.bilfile.BilFile In [4]: print(img) Data Source: '/home/thomas/spectral_data/92AV3C.lan' # Rows: 145 # Samples: 145 # Bands: 220 Interleave: BIL Quantization: 16 bits Data format: int16 HS画像はワーキングディレクトリには存在しないが、SPECTRAL＿DATA環境変数のディレクトリから探している。また画像のピクセルデータがlineごとにinterleaveされているため、image関数はBilFileインスタンスを返す。(BILはデータの格納方式) HSデータファイルはサイズが大きいため、SpyFileオブジェクトが作られるときメタデータだけが読み込まれる。実際のイメージデータの値はSpyFileメソッドによってリクエストされた時のみ読み込まれる。Numpyのように配列に対して添字演算子(たぶんインデックスとスライシングとか)を適用することができる。SpyFileオブジェクトはM×N×Bの配列でMが画像の行数、Nが画像の列数、Bがバンド数にそれぞれなる。thenumberになってる。 In [5]: img.shape Out[5]: (145, 145, 220) In [6]: pixel = img[50,100] In [7]: pixel.shape Out[7]: (220,) In [8]: band6 = img[:,:,5] In [9]: band6.shape Out[9]: (145, 145, 1) 画像データの各ピクセルの値は添字演算子が適用されるまで(インデクシングやスライシングされるまで)読み込まれない。またPythonではインデックスは0から始まるため、img[50, 100]のような場合、51行、101列目のピクセルデータが読み込まれる。同様に、img[:,:,6]では画像全体での6バンド目の情報が読み込まれる。 SpyFileにはサブクラスのインスタンスが存在し、特定のイメージファイルを返すメソッドが実装されている。 Method Description read_band Reads a single band into an M×N array read_bands Reads multiple bands into an M×N×C array read_pixel Reads a single pixel into a length B array read_subregion Reads multiple bands from a rectangular sub-reagion of the image read_subimage Reads specified rows, columns, and bands SpyFileオブジェクトはbandsメンバーを持っており、このメンバーはイメージスペクトルバンドについての補足的な情報を持つBandInfoオブジェクトのインスタンスである。 Loading Entire Images 繰り返しにはなるが、画像データはSpyFileオブジェクトによって必要となったときに呼び出され、データはキャッシュされない。添字演算子やメソッドの呼び出しごとに対応する画像データファイルが読み込まれる。これは事前に同じデータが読み込まれていたとしても関係はない。都度読み込みを行うことで、巨大な画像データに対して処理を行う際に激しいメモリの消費を防ぐ。また画像データ中の一部だけを利用する際には無駄を省くことができる(画像表示のためにRGBだけ読み込むなど)。 しかし、データ全体に対して適用するようなアルゴリズムを走らせた場合には処理に時間がかかる。さらに何度も繰り返しデータにアクセスする必要がるアルゴリズムではよりひどい。 前述のような欠点を解消するために、loadメソッドを利用しデータ全体をメモリにロードすることが望ましい。このメソッドではImageArrayオブジェクトが返されるが、SpyFileインターフェースと同様にnumpy.ndarrayのインターフェース（インデクシングとかスライシングとか）をもつ。 In [1]: arr = img.load() In [2]: arr.class Out[2]: spectral.image.ImageArray In [3]: print(arr.info()) # Rows: 145 # Samples: 145 # Bands: 220 Data format: float32 In [4]: arr.shape Out[4]: (145, 145, 220) SPyは主にスペクトル分野での処理用途にデザインされているため、メモリ中のspectral.ImageArrayオブジェクトはもとのデータファイルのインターリーブスに関係なく常にピクセルごとのインターリーブスをもつ。言い換えるならば、numpy.ndarrayのshapeは(numRows, numCols, numBands)である。また配列の要素は常に32-bit floatsのデータ型である。 Note: loadメソッドを呼び出す前に、出力されうるImageArrayオブジェクトの消費メモリ量を考慮すべきである。spectral.ImageArrayオブジェクトは32-bit 浮動小数点の要素を持つため、メモリ消費量は4*numRows*numCols*numBandsバイトと概算される。 Numpy memmap Interface ファイル全体を読み込んでおく方法は他にもあり、遅いもののメモリ効率のよいnumpy memmapオブジェクトを利用した方法がある。ここではSpyFileオブジェクトをopen_memmapメソッドで開く。このmemmapオブジェクトはdate（dataの間違い？）をイメージファイルに書き込むことにも使われる。 File Formats Supported ENVI Headers ENVI*は地理空間画像情報を処理し解析するための商用ソフトウェアパッケージのひとつである。SPyはENVIヘッダーでまとめられた画像を読み込み、ENVIヘッダーでスペクトルライブラリの読み込みと書き込みができる。ENVIファイルはSPyのimage関数で自動的に開かれるが、明示的にENVIファイルと指定して開くことも可能である。ヘッダーとデータのディレクトリが異なる場合や、データがSPyで特定できない拡張子を持つ場合、ENVIファイルと明示的に指定しデータをひらく必要がある。 *ENVI is a registered trademark of Exelis Visual Information Solutions. In [5]: import spectral.io.envi as envi In [6]: img = envi.open('cup95eff.int.hdr', 'cup95eff.int') In [7]: import spectral.io.envi as envi In [8]: lib = envi.open('spectra.hdr') In [9]: lib.names[:5] Out[9]: ['construction asphalt', 'construction concrete', 'red smooth-faced brick', 'weathered red brick', 'bare red brick'] See also: ファイルに画像データを書き込むための関数: create_image: ディスク上の割り当てられたストレージに新しい画像ファイルを作成 save_image: ENVIヘッダーで既存の画像、またはndarrayを保存 AVIRIS SPyはAirborne Visible/Infrared Imaging Spectrometer (AVIRIS)**によって生成されたデータファイルをサポートしている。AVIRISファイルはopen_image関数によって自動的に認識されるが、スペクトルバンドのカリブレーションファイルは自動的には認識されない。そのため、それぞれのファイルを明示的に開くことになるかもしれない。 In [10]: img = aviris.open('f970619t01p02_r02_sc01.a.rfl', 'f970619t01p02_r02.a.spc') まあバンドカリブレーションファイルを別々に読み込むこともできる(AVIRISフォーマット中にバンドカリブレーションファイルはあるが、イメージファイルがない場合などに必要とりうる。) In [11]: img = open_image('92AV3C.lan') In [12]: img.bands = aviris.read_aviris_bands('92AV3C.spc') **AVIRIS ERDAS/Lan ERDAS/Lanファイルフォーマットはimage関数により自動的に認識される。ファイルをLanファイルとして明示的に開く必要があるとは考えにくいが、以下のようにして開くことができる。 In [13]: import spectral.io.erdas as erdas In [14]: img = erdas.open('92AV3C.lan')

Spectral Python Spectral Python(SPy)とは純度100%のpythonのモジュールで、ハイパースペクトル(HS)データの処理ができる。HSデータの読み込み、表示、計算そして分類まで一括して行える。MITライセンスのため、自由に使うことができる。 Introduction このユーザーガイドでは、下記の論文のデータを参照しチュートリアルを進める。またチュートリアルでは利便性のためfrom spectral import *としているが、モジュール内ネームスペースのコンテンツをすべてインポートすることは非推奨である。データはここから。 File Name Description 92AV3C.lan ERDAS/LanフォーマットのHS画像。AVRISでの撮影で得られた145*145ピクセル、220バンドのデータ。 92AV3GT.GIS 上のデータの土地利用図。 92AV3C.spc 上のデータのAVIRISフォーマットのバンドカリブレーションファイル Landgrebe, D. Multispectral data analysis: A signal theory perspective. School of Electr. Comput. Eng., Purdue Univ., West Lafayette, IN (1998). Installing Spy 基本的にはPythonの組み込みのものとNumpyがあればHSデータの処理は可能である。グラフィカルな機能を使用して解析を行いたい場合には別途モジュールが必要となる。下表にそれぞれの機能の要件を示す。 Dependency Notes Python 2.6+ or 3.3+ (1) Numpy Required Pillow or Python Imaging Library(PIL) Required if displaying or saving images wxPython (2) matplotlib Required if rendering raster displays or spectral plots IPython Required for interactive, non-blocking GUI windows PyOpenGL (2) Notes: (1): view_cube、view_nd以外はPython 3.3+で動作 (2): view_cube、view_ndを呼び出す際に必要 pip でのインストールは下記のコマンドで。 その他の方法はここから pip install spectral 基本的にHSデータは容量が大きいため、特定のディレクトリに保存されている。SPyでデータを開こうとするたびに絶対パスを入力しなくて済むようにSPECTRAL_DATAという環境変数を設定できる。SPyはこの環境変数に記されているディレクトリからイメージファイルを探す。またディレクトリパスのリストは:区切りで記載されなければならない。 Reading HSI Data Files 通常、SPyではHSイメージファイルにアクセスし開くためにimage関数を使用する。この関数はSpyFileオブジェクトのインスタンスを返す。 The SpyFile Interface SpyFileはHSデータを読み込むオブジェクトの基本的なクラスである。SpyFileオブジェクトの作成時、対応するファイルからデータを読み込むインターフェースが与えられる。画像を開いたとき、返される実際のオブジェクトはSpyFileのサブクラス(BipFile, BilFile, BsqFile)で、このサブクラスはデータの格納方法(interleaves)にそれぞれ対応している。 下記のコードではサンプル画像を開いている。 In [1]: from spectral import * In [2]: img = open_image('92AV3C.lan') In [3]: img.class Out[3]: spectral.io.bilfile.BilFile In [4]: print(img) Data Source: '/home/thomas/spectral_data/92AV3C.lan' # Rows: 145 # Samples: 145 # Bands: 220 Interleave: BIL Quantization: 16 bits Data format: int16 HS画像はワーキングディレクトリには存在しないが、SPECTRAL＿DATA環境変数のディレクトリから探している。また画像のピクセルデータがlineごとにinterleaveされているため、image関数はBilFileインスタンスを返す。(BILはデータの格納方式) HSデータファイルはサイズが大きいため、SpyFileオブジェクトが作られるときメタデータだけが読み込まれる。実際のイメージデータの値はSpyFileメソッドによってリクエストされた時のみ読み込まれる。Numpyのように配列に対して添字演算子(たぶんインデックスとスライシングとか)を適用することができる。SpyFileオブジェクトはM×N×Bの配列でMが画像の行数、Nが画像の列数、Bがバンド数にそれぞれなる。thenumberになってる。 In [5]: img.shape Out[5]: (145, 145, 220) In [6]: pixel = img[50,100] In [7]: pixel.shape Out[7]: (220,) In [8]: band6 = img[:,:,5] In [9]: band6.shape Out[9]: (145, 145, 1) 画像データの各ピクセルの値は添字演算子が適用されるまで(インデクシングやスライシングされるまで)読み込まれない。またPythonではインデックスは0から始まるため、img[50, 100]のような場合、51行、101列目のピクセルデータが読み込まれる。同様に、img[:,:,6]では画像全体での6バンド目の情報が読み込まれる。 SpyFileにはサブクラスのインスタンスが存在し、特定のイメージファイルを返すメソッドが実装されている。 Method Description read_band Reads a single band into an M×N array read_bands Reads multiple bands into an M×N×C array read_pixel Reads a single pixel into a length B array read_subregion Reads multiple bands from a rectangular sub-reagion of the image read_subimage Reads specified rows, columns, and bands SpyFileオブジェクトはbandsメンバーを持っており、このメンバーはイメージスペクトルバンドについての補足的な情報を持つBandInfoオブジェクトのインスタンスである。 Loading Entire Images 繰り返しにはなるが、画像データはSpyFileオブジェクトによって必要となったときに呼び出され、データはキャッシュされない。添字演算子やメソッドの呼び出しごとに対応する画像データファイルが読み込まれる。これは事前に同じデータが読み込まれていたとしても関係はない。都度読み込みを行うことで、巨大な画像データに対して処理を行う際に激しいメモリの消費を防ぐ。また画像データ中の一部だけを利用する際には無駄を省くことができる(画像表示のためにRGBだけ読み込むなど)。 しかし、データ全体に対して適用するようなアルゴリズムを走らせた場合には処理に時間がかかる。さらに何度も繰り返しデータにアクセスする必要がるアルゴリズムではよりひどい。 前述のような欠点を解消するために、loadメソッドを利用しデータ全体をメモリにロードすることが望ましい。このメソッドではImageArrayオブジェクトが返されるが、SpyFileインターフェースと同様にnumpy.ndarrayのインターフェース（インデクシングとかスライシングとか）をもつ。 In [1]: arr = img.load() In [2]: arr.class Out[2]: spectral.image.ImageArray In [3]: print(arr.info()) # Rows: 145 # Samples: 145 # Bands: 220 Data format: float32 In [4]: arr.shape Out[4]: (145, 145, 220) SPyは主にスペクトル分野での処理用途にデザインされているため、メモリ中のspectral.ImageArrayオブジェクトはもとのデータファイルのインターリーブスに関係なく常にピクセルごとのインターリーブスをもつ。言い換えるならば、numpy.ndarrayのshapeは(numRows, numCols, numBands)である。また配列の要素は常に32-bit floatsのデータ型である。 Note: loadメソッドを呼び出す前に、出力されうるImageArrayオブジェクトの消費メモリ量を考慮すべきである。spectral.ImageArrayオブジェクトは32-bit 浮動小数点の要素を持つため、メモリ消費量は4*numRows*numCols*numBandsバイトと概算される。 Numpy memmap Interface ファイル全体を読み込んでおく方法は他にもあり、遅いもののメモリ効率のよいnumpy memmapオブジェクトを利用した方法がある。ここではSpyFileオブジェクトをopen_memmapメソッドで開く。このmemmapオブジェクトはdate（dataの間違い？）をイメージファイルに書き込むことにも使われる。 File Formats Supported ENVI Headers ENVI*は地理空間画像情報を処理し解析するための商用ソフトウェアパッケージのひとつである。SPyはENVIヘッダーでまとめられた画像を読み込み、ENVIヘッダーでスペクトルライブラリの読み込みと書き込みができる。ENVIファイルはSPyのimage関数で自動的に開かれるが、明示的にENVIファイルと指定して開くことも可能である。ヘッダーとデータのディレクトリが異なる場合や、データがSPyで特定できない拡張子を持つ場合、ENVIファイルと明示的に指定しデータをひらく必要がある。 *ENVI is a registered trademark of Exelis Visual Information Solutions. In [5]: import spectral.io.envi as envi In [6]: img = envi.open('cup95eff.int.hdr', 'cup95eff.int') In [7]: import spectral.io.envi as envi In [8]: lib = envi.open('spectra.hdr') In [9]: lib.names[:5] Out[9]: ['construction asphalt', 'construction concrete', 'red smooth-faced brick', 'weathered red brick', 'bare red brick'] See also: ファイルに画像データを書き込むための関数: create_image: ディスク上の割り当てられたストレージに新しい画像ファイルを作成 save_image: ENVIヘッダーで既存の画像、またはndarrayを保存 AVIRIS SPyはAirborne Visible/Infrared Imaging Spectrometer (AVIRIS)**によって生成されたデータファイルをサポートしている。AVIRISファイルはopen_image関数によって自動的に認識されるが、スペクトルバンドのカリブレーションファイルは自動的には認識されない。そのため、それぞれのファイルを明示的に開くことになるかもしれない。 In [10]: img = aviris.open('f970619t01p02_r02_sc01.a.rfl', 'f970619t01p02_r02.a.spc') まあバンドカリブレーションファイルを別々に読み込むこともできる(AVIRISフォーマット中にバンドカリブレーションファイルはあるが、イメージファイルがない場合などに必要とりうる。) In [11]: img = open_image('92AV3C.lan') In [12]: img.bands = aviris.read_aviris_bands('92AV3C.spc') **AVIRIS ERDAS/Lan ERDAS/Lanファイルフォーマットはimage関数により自動的に認識される。ファイルをLanファイルとして明示的に開く必要があるとは考えにくいが、以下のようにして開くことができる。 In [13]: import spectral.io.erdas as erdas In [14]: img = erdas.open('92AV3C.lan')

備忘録：実行ファイルより上層に配置した自作モジュールをimport いつも躓くので書いておく 結論 moduleがあるディレクトリのpathをPYTHONPATHへ追加する # .pyの場合 import os import sys sys.path.append(os.path.dirname(__file__), '実行ファイルからの相対path') # .ipynbの場合 import os import sys sys.path.append(os.path.join(Path().resolve(), '実行ファイルからの相対path') ipynbではfile変数を参照できないのでPath()から参照する 通常pythonは環境変数PYTHONPATHの上位から順に探索してるらしい 感想 これで読み込めるが、sys.path.append("相対path")との違いがよくわかってない 首記のコードでは実行ファイルまでのpathに相対pathをつなげて(絶対pathとして)いるのだろうけど 備考 PYTHONPATHをコード内でいじらず、pthファイルに記述する方法もある

今日は久しぶりにパイソンを復習しようと思います。 その中でパイソンを利用してファイルをダウンロードする方法を調べてみようと思います。 方法 方法はrequestモジュールを活用することです。 pip install requests requestモジュールを使用するには、まずpip install requestsを使用してrequestsモジュールをインストールする必要があります。 import requests url = "http://google.com/favicon.ico"//..... 接続したいアドレス... req = requests.get(url, allow_redirects = True) // getを使ってもいいし、postを使ってもいいし、別の接続方式を使ってもいいです。 open('google.ico','wb').write(req.content) ここからより多様な使い方を確認することができます。

前回の投稿は「小学校の算数の授業で習った数直線=スカラーの推移のみを抽出した同心円集合」から出発しました。 【Pyrhon演算処理】同心集合(Concentric Set)について。①乗法的同心集合(Multiplicative Concentric Set)とは？ しかし残念ながら最初に着手した乗法的同心集合(Multiplicative Concentric Set)は、どれとしてこの形に到達できなかったばかりか、その本質に近づくほど同じものとは見えなくなっていったのです。 正方形が構築する同心円集合(演算$2^{\frac{n}{2}}$の結果集合) 球面の半径が構築する同心集合の一例(演算$\sqrt{3}*2^{\frac{n}{2}}$の結果集合) 正三角形が構築する同心円集合(演算$2^n$の結果集合) 演算$2^n$の結果集合 演算$e^n$の結果集合 この問題を解決するには加法群(Additive Group)概念の導入、あるいは(演算としてスカラー倍と加減算が保証された)ベクトル空間(Vector Space)への移行が必須となるのです。 加法的同心集合(Additive Concentric Set) さっそく同心集合の加法群定義に入りたいと思います。 【数理考古学】群論概念(Group Theory Concept)①基本定義 まずは元集合Aを用意する。例えばそれは乗法単位元(Multiplicative Identity)1を備えた乗法群(Multiplicative Group)であっても良い。 加法単位元(Additive Identity)0と逆元(Inverse Element)$A^{-1}$を準備する。 全てを結合する。 その結果として(規定演算によって球面座標系やトーラス座標系によって把握されてきた)元集合が、円柱座標系(Cylindrical Coordinate System)や円錐座標系(Conical Coordinate System)で扱える様になります。 円筒座標系 - Wikipedia 与えられた点Pの位置を以下の三つの座標成分(p,φ,z)で表す。ただし基準平面からの「距離」はその点が基準平面の（表または裏の）どちら側に面するかによって正または負の値を持つものとする。 軸距離(Axial Distance)または動径距離(Radial Distance)ρ…z-軸から点 P までのユークリッド距離(別表現では「特別に選ばれた基準軸からの距離」); 方位角(Azimuth)φ…基準平面上の基準方向と原点から点Pの基準平面への正射影へ結んだ直線との間の角(別表現では「特別に選ばれた基準方向に対する軸から測った方向」); 軸座標 (axial coordinate) または高さ (height) z は基準平面から点 P までの符号付き距離(基準軸に直交する特別に選ばれた基準平面からの距離)。 円筒座標系の原点(origin)とは、上記三つの座標成分がすべて0として与えることができるような点を言う。これは上記の基準平面と基準軸の交点になる。 この基準軸は円筒軸(Cylindrical Axis; 円柱軸)や緯線(Longitudinal Axis) などと呼ばれ、極軸(Polar Axis)とは区別される（極軸あるいは始線とは、基準平面上に載っている半直線で、原点から出てその平面上の基準方向を指し示すものを表すために用いられる）。緯線に垂直な任意の方向を射線(Radial Lines;放射状の直線)と総称する。 円筒軸からの距離は動径距離(Radial Distance;放射距離)や動径(Radius) と呼ばれ、円筒軸周りの偏角座標は角度位置(Angular Position)や方位角(Azimuth)などと言う。動径成分と方位角成分を併せて極座標成分(Polar Coordinates)と言い、これはその点を通り基準平面に平行な平面上の二次元の極座標系に対応する。第三座標は基準平面を水平面と見るとき高さ(Height)や高度(Altitude)と呼んだり、緯度(Longitudinal Position)や軸位置(Axial Position)などとも言う。 「第三座標は基準平面を水平面と見るとき高さ(Height)や高度(Altitude)と呼んだり、緯度(Longitudinal Position)や軸位置(Axial Position)などとも言う。」…率直に「時間軸(t)」を含めない理由が分かりません。 円筒座標系は緯線周りの何らかの回転対称性を持つ物体や現象（例えば、丸い断面を持つ直線パイプを流れる水流や、金属円柱の熱分布、長い真っ直ぐなワイヤー内の電荷から出る電場、天文学における降着円盤など）との関連で有意である。 一方、この射影(Projection)によって副次的に生じる円錐座標系は、例えば元元が乗法群であった場合、その写像(Map)となる。 射影（projection）とは - IT用語辞典 e-Words 円錐の体積を円柱座標系で真面目に計算してみた 「円筒座標の導入自体は群規定導入以前でも可能である」なる考え方も不可能ではないが、こういう議論もあるのでその場合にどういう制限を受ける事になるかについて十分な追加記述を必要とするものとする。 そういえば円錐座標系にも「HSV色空間モデル(Hue=色相,Saturation=彩度,Saturation=明度)」の様に逆元は備えていないが(あるいは元そのものとなるだけなので設定する意味がない)加法単位元(黒一色)と乗法単位元(基準色相環を実現する彩度と明度…あれ、次元それぞれに一つずつある？)なら備えている。こうした複雑な中間例が他にも幾らでも見つかりそうだったりする訳である。円錐座標系の意味・用法を知る彩度 - さいど | 武蔵野美術大学 造形ファイル そして、いわゆる「ヴェーバー‐フェヒナーの法則」に従うなら、例えば「HSV色空間モデル」も「音空間モデル」同様、その目盛りが均等尺(Even Scale)ではなく対数尺(Logarithmic Scale)かもしれないのである。【数理考古学】ヴェーバー‐フェヒナーの法則【数理考古学】「音階論」なる古くて新しい数理空間 かかる形での乗法群の加法群への写像がもたらす最大の恩恵は「傾き」の概念の追加というべきかもしれません。そう「(スカラー倍と加減算しか演算が存在しない)ベクトル空間」から「(角度や内積の概念が追加された)計量ベクトル空間」への拡張です。 傾き (数学) - Wikipedia 傾きmは傾斜角θによりm=tan(θ)と定義される。また平面上の直線の傾きは、垂直移動距離Δy($y_2-y_1$)を水平移動距離Δx($x_2-x_1$)で割ったm=Δy/Δxで定義される。 これらの等式が示す様に鉛直線（y軸に平行な直線）の傾きは、零除算となり、定義されない。 全体像としてはこういうイメージですね。 %matplotlib nbagg import math as m import cmath as c import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import matplotlib.animation as animation #円柱データ作成 c0=np.linspace(0,m.pi*120,1201,endpoint = True) s0=[] for nm in range(len(c0)): s0.append(complex(m.cos(c0[nm]),m.sin(c0[nm]))) s1=np.array(s0) z0=np.linspace(-1,1,1201,endpoint = True) #単位円データ作成 u0=np.linspace(0,m.pi*2,61,endpoint = True) u1=[] for nm in range(len(u0)): u1.append(complex(m.cos(u0[nm]),m.sin(u0[nm]))) uc=np.array(u1) uz0=np.repeat(-1,61) uz1=np.repeat(-0,61) uz2=np.repeat(1,61) #グラフ表示 plt.style.use('default') fig = plt.figure() ax = Axes3D(fig) #関数定義 def unit_cylinder(n): plt.cla() #円柱高計算 h1=2*n/len(uc)-1 hz=np.repeat(h1,len(uc)) #円柱描画 ax.plot(s1.real,s1.imag,z0,color="gray",lw=0.5) #スポーク描画 #for num in range(len(uc)): # ax.plot([0,uc[num].real],[0,uc[num].imag],[0,0],color="gray",lw=0.5) #for num in range(len(uc)): # ax.plot([0,uc[num].real],[0,uc[num].imag],[1,1],color="gray",lw=0.5) for num in range(len(uc)): ax.plot([0,uc[num].real],[0,uc[num].imag],[0,0],color="olivedrab",lw=0.5) #単位円描画 ax.plot(uc.real,uc.imag,uz0,color="red",lw=1) ax.plot(uc.real,uc.imag,uz1,color="green",lw=1) ax.plot(uc.real,uc.imag,uz2,color="blue",lw=1) #実数線追加 ax.plot([0,0],[0,0],[-1,1],color="black",lw=1) ax.plot([0,1],[0,0],[-1,-1],color="red",lw=1) ax.plot([0,1],[0,0],[0,0],color="black",lw=1) ax.plot([0,1],[0,0],[1,1],color="blue",lw=1) ax.plot([1,1],[0,0],[-1,0],color="red",lw=1) ax.plot([1,1],[0,0],[0,1],color="blue",lw=1) #移動円追加 for num in range(len(uc)): ax.plot([0,uc[num].real],[0,uc[num].imag],[h1,h1],color="black",lw=0.5) ax.plot(uc.real,uc.imag,hz,color="black",lw=1) ax.plot([0,1],[0,0],[0,h1],color="black",lw=1) #諸元追加 ax.set_ylim([-1.1,1.1]) ax.set_xlim([-1.1,1.1]) ax.set_zlim([-1.1,1.1]) ax.set_title("Unit Cylinder") ax.set_xlabel("Real") ax.set_ylabel("Imaginal") ax.set_zlabel("Cycle") # グラフを回転(elev=0にすると水平表示、90にすると垂直表示) ax.view_init(elev=45, azim=Time_code[n]) Time_code0=np.arange(0,360,6) Time_code=Time_code0[::-1] #unit_cylinder(len(s1)) #plt.show() ani = animation.FuncAnimation(fig, unit_cylinder, interval=50,frames=len(Time_code)) ani.save("cylinder003.gif", writer="pillow") 【Python演算処理】半径・直径・円周長・円の面積・球の表面積・球の体積の計算上の往復 1次元球の体積は2rとなり、既知の結果と一致する。 これまでこうした表現を無条件で受容してきましたが、実は加法群が成立する裏側には単位円が隠れていた(直交、すなわち内積0となる形で交わっていた)という訳です。だから「十進法導入により周期を無限に細かく刻んでいく事によって実数列を近似する」といった方便が必要となった訳ですね。ある意味これこそが群論でいう「環(Ring)」概念そのもの？ 【数理考古学】とある実数列(Real Sequance)の規定例①等差数列から加法整数群へ ではこの考え方にこの投稿の主題たる同心集合の概念を重ねてみましょう。 【Pyrhon演算処理】同心集合(Concentric Set)について。①乗法的同心集合(Multiplicative Concentric Set)とは？ %matplotlib nbagg import math as m import cmath as c import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import matplotlib.animation as animation #円柱データ作成 c0=np.linspace(0,m.pi*120,1201,endpoint = True) s0=[] for nm in range(len(c0)): s0.append(complex(m.cos(c0[nm]),m.sin(c0[nm]))) s1=np.array(s0) z0=np.linspace(-1,1,1201,endpoint = True) #単位円データ作成 u0=np.linspace(0,m.pi*2,61,endpoint = True) u1=[] for nm in range(len(u0)): u1.append(complex(m.cos(u0[nm]),m.sin(u0[nm]))) uc=np.array(u1) uz0=np.repeat(-1,61) uz1=np.repeat(-0,61) uz2=np.repeat(1,61) #グラフ表示 plt.style.use('default') fig = plt.figure() ax = Axes3D(fig) #関数定義 def unit_cylinder(n): plt.cla() #円柱高計算 h1=2*n/len(uc)-1 hz=np.repeat(h1,len(uc)) #円柱描画 ax.plot(s1.real,s1.imag,z0,color="gray",lw=0.5) #半円柱描画 ax.plot(s1.real/2,s1.imag/2,z0,color="gray",lw=0.5) #スポーク描画 #for num in range(len(uc)): # ax.plot([0,uc[num].real],[0,uc[num].imag],[0,0],color="gray",lw=0.5) #for num in range(len(uc)): # ax.plot([0,uc[num].real],[0,uc[num].imag],[1,1],color="gray",lw=0.5) for num in range(len(uc)): ax.plot([0,uc[num].real],[0,uc[num].imag],[0,0],color="olivedrab",lw=0.5) #単位円描画 ax.plot(uc.real,uc.imag,uz0,color="red",lw=1) ax.plot(uc.real,uc.imag,uz1,color="green",lw=1) ax.plot(uc.real,uc.imag,uz2,color="blue",lw=1) #半単位円描画 ax.plot(uc.real/2,uc.imag/2,uz0,color="red",lw=1) ax.plot(uc.real/2,uc.imag/2,uz1,color="green",lw=1) ax.plot(uc.real/2,uc.imag/2,uz2,color="blue",lw=1) #実数線追加 ax.plot([0,0],[0,0],[-1,1],color="black",lw=1) ax.plot([0,1],[0,0],[-1,-1],color="red",lw=1) ax.plot([0,1],[0,0],[0,0],color="black",lw=1) ax.plot([0,1],[0,0],[1,1],color="blue",lw=1) ax.plot([1,1],[0,0],[-1,0],color="red",lw=1) ax.plot([1,1],[0,0],[0,1],color="blue",lw=1) #半実数線追加 ax.plot([1/2,1/2],[0,0],[-1,0],color="red",lw=1) ax.plot([1/2,1/2],[0,0],[0,1],color="blue",lw=1) #移動円追加 for num in range(len(uc)): ax.plot([0,uc[num].real],[0,uc[num].imag],[h1,h1],color="black",lw=0.5) ax.plot(uc.real,uc.imag,hz,color="black",lw=1) ax.plot(uc.real/2,uc.imag/2,hz,color="black",lw=1) ax.plot([0,1],[0,0],[0,h1],color="black",lw=1) ax.plot([0,1/2],[0,0],[0,h1],color="black",lw=1) #諸元追加 ax.set_ylim([-1.1,1.1]) ax.set_xlim([-1.1,1.1]) ax.set_zlim([-1.1,1.1]) ax.set_title("Unit Cylinder") ax.set_xlabel("Real") ax.set_ylabel("Imaginal") ax.set_zlabel("Cycle") # グラフを回転(elev=0にすると水平表示、90にすると垂直表示) ax.view_init(elev=45, azim=Time_code[n]) Time_code0=np.arange(0,360,6) Time_code=Time_code0[::-1] #unit_cylinder(len(s1)) #plt.show() ani = animation.FuncAnimation(fig, unit_cylinder, interval=50,frames=len(Time_code)) ani.save("cylinder101.gif", writer="pillow") 1を単位元とする乗法同心群は、その半径が極限たる0と∞(トーラス座標系ではさらに両者を同一視し、連結させて消してしまう!!)に差し掛かった時に計算不可能となりますが「傾き」概念にその特徴が継承されている事が分かります。 円筒座標系 円柱の中心軸からの半径方向（r軸）、周方向（θ軸）および円柱の軸方向（z軸）の3軸から構成される座標系のことで「円柱座標系」と呼ばれることもあります。 直交座標系(x,y,z)と円筒座標系(r,φ,z)の間には以下の関係が成り立ちます。 x=r×cos(θ) y=r×sin(θ) z=z $r=\sqrt{x^2+y^2}$ $φ=sign(y)atan2(y,x)$ import sympy as sp import numpy as np from sympy import I, pi, E x,y,z,r,φ,m= sp.symbols('x,y,z,r,φ,m') #sp.Eq(z,z)と書いてもTrueという #論理演算結果が返るのみなので省略 #デカルト座標系(x,y,z)と極座標系(r,φ,θ)の相互変換 eq01=sp.Eq(r,sp.sqrt(x**2+y**2)) eq02=sp.Eq(φ,sp.sign(y)*sp.atan2(y,x)) eq03=sp.Eq(x,r*sp.cos(φ)) eq04=sp.Eq(y,r*sp.sin(φ)) eq05=sp.Eq(m,sp.tan(φ)) eq06=sp.Eq(m,z/r) #単位円(半径r=1)の場合 eq11=eq01.subs(r,1) eq12=eq02.subs(r,1) eq13=eq03.subs(r,1) eq14=eq04.subs(r,1) eq15=eq05.subs(r,1) eq16=eq06.subs(r,1) #tex sp.init_printing() print("デカルト座標系(x,y,z)と極座標系(r,φ,z)の相互変換") display(eq01) print(sp.latex(eq01)) display(eq02) print(sp.latex(eq02)) display(eq03) print(sp.latex(eq03)) display(eq04) print(sp.latex(eq04)) display(eq05) print(sp.latex(eq05)) display(eq06) print(sp.latex(eq06)) print("単位円(半径r=1)の場合") display(eq11) print(sp.latex(eq11)) display(eq12) print(sp.latex(eq12)) display(eq13) print(sp.latex(eq13)) display(eq14) print(sp.latex(eq14)) display(eq15) print(sp.latex(eq15)) display(eq16) print(sp.latex(eq16)) デカルト座標系(x,y,z)と極座標系(r,φ,z)の相互変換 \begin{align*} &r = \sqrt{x^{2} + y^{2}}\\ &φ = \operatorname{atan_{2}}{\left(y,x \right)} \operatorname{sign}{\left(y \right)}\\ &x = r \cos{\left(φ \right)}\\ &y = r \sin{\left(φ \right)}\\ &z = z\\ &m = \tan{\left(φ \right)}\\ &m = \frac{z}{r} \end{align*} 単位円(半径r=1)の場合(実数列$\mathbb{Z}$そのものとなる)。 \begin{align*} &r = 1\\ &φ = \operatorname{atan_{2}}{\left(y,x \right)} \operatorname{sign}{\left(y \right)}\\ &x = \cos{\left(φ \right)}\\ &y = \sin{\left(φ \right)}\\ &z = z\\ &m = \tan{\left(φ \right)}\\ &m = z \end{align*} こうして全体像を俯瞰した時、初めて「小学校の算数の授業で習った数直線=スカラーの推移のみを抽出した同心円集合」をより高次元の数学の領域で再掌握したといえるのかもしれません。そういう私自身、こうした考え方に辿り着く事によってやっと群論における環の概念を理解出来たのです（この理解で本当に数学的に合っているか、まだ全然自身がなかったりもする訳ったりもするのでしたすが。） とりあえず、そんな感じで以下続報…