pythonでffmpegを操作したい場合、ffmpeg-pythonを使うと思います。 ffmpeg.output()にはffmpegのオプションを渡すことができます。そこで、ffmpeg-pythonで動画にメタデータを追加するには、ffmpeg.output()に引数名metadataで"key=value"の文字列の形式でメタデータを渡します。 ffmpeg .input(input_video_path) .output(output_video_path, metadata="title=New Title") .run() ffmpeg.output()は可変長引数を取るので、複数のメタデータを設定したい場合はその分だけ引数を設定します。 ffmpeg .input(input_video_path) .output(output_video_path, metadata="title=New Title", metadata="title=New Artist", metadata="album=New Album") .run() 辞書型の引数を展開することもできます。 ffmpeg .input(input_video_path) .output(output_video_path, **{"metadata": "title=New Title", "metadata": "artist=New Artist", "metadata": "album=New Album"}) .run() 辞書型の引数を変数で保持しておきたい場合は少し工夫が必要です。キーをそのままmetadataにすると、キーが重複して上書きされてしまうためです。 キーをユニークにするためには、metadata:g:0のようにグローバルオプションgを付けて、更にその後にユニークなインデックスを振るようにします。グローバルオプションgを付けるとメタデータがファイル全体に適用されますが、メタデータはデフォルトでグローバルになっており、グローバルオプションを明示的に付けても省略した場合と同じ動作になるため、このようにしても無害です。また、グローバルオプションの後のインデックスは無視されます。 なお、グローバルオプションを付けずにインデックスだけを振ってmetadata:0とするとエラーになります。 metadata_dict = {"metadata:g:0": "title=New Title", "metadata:g:1": "artist=New Artist", "metadata:g:2": "album=New Album"} ffmpeg .input(input_video_path) .output(output_video_path, **metadata_dict) .run() このハックを使うと、次のように別のファイルからメタデータをコピーすることもできます。 tags = ffmpeg.probe(input_video_path)["format"]["tags"] metadata_dict = {} for i, key in enumerate(tags): metadata_dict[f"metadata:g:{i}"] = f"{key}={tags[key]}" ffmpeg .input(input_video_path) .output(output_video_path, **metadata_dict) .run() 参考

自分のようなど素人しかハマらない箇所かもしれませんが、うっかり他の人がハマった時の助けになるよう、記録しておきます。 発生したエラー １３６ページ、PyCharmを通じてrunserverを実行しサーバーを起動する操作でエラー発生。 django.core.exceptions.ImproperlyConfigured: Error loading psycopg2 module: No module named 'psycopg2' psycopg2がないという。確かに１０１ページでインストールしたのはpsycopg2-binaryだった。psycopg2を別途インストールしないといけないんだろうか？？？ 顛末 PyCharm上でrunserverを登録する際、pythonインタープリターが仮想環境内のものになっていなかったことが原因だった。このpythonにはpsycopg2-binaryがインストールされていないためエラーとなっていた。 なおこの設定は同書117ページにて案内されているが、上手くできていなかった。

はじめに MTCNNを使って顔検出するときにどんな情報が得られるのか知りたくて、ラズパイでサンプルコードを実行してみました。 全ての情報は以下にあります。 実行環境 ライブラリ バージョン python 3.7.3 mtcnn 0.1.1 opencv-python 4.5.3.56 tensorflow 1.14.0 やったこと example.pyの内容を少しいじって、以下のコードを実行しました。 import cv2 from mtcnn import MTCNN filename = "ivan.jpg" output = "drow_" + filename detector = MTCNN() image = cv2.cvtColor(cv2.imread(filename), cv2.COLOR_BGR2RGB) result = detector.detect_faces(image) # Result is an array with all the bounding boxes detected. We know that for 'ivan.jpg' there is only one. bounding_box = result[0]['box'] keypoints = result[0]['keypoints'] cv2.rectangle(image, (bounding_box[0], bounding_box[1]), (bounding_box[0]+bounding_box[2], bounding_box[1] + bounding_box[3]), (0,155,255), 2) cv2.circle(image,(keypoints['left_eye']), 2, (0,155,255), 2) cv2.circle(image,(keypoints['right_eye']), 2, (0,155,255), 2) cv2.circle(image,(keypoints['nose']), 2, (0,155,255), 2) cv2.circle(image,(keypoints['mouth_left']), 2, (0,155,255), 2) cv2.circle(image,(keypoints['mouth_right']), 2, (0,155,255), 2) cv2.imwrite(output, cv2.cvtColor(image, cv2.COLOR_RGB2BGR)) print(result) 実行結果 顔を検出した箇所にマーキングされました。画像のサイズは561×561のファイルです。 resultは以下のように出力されました。※見やすいように改行入れてます。 [{'box': [277, 93, 49, 62], 'confidence': 0.9999709725379944, 'keypoints': {'left_eye': (291, 117), 'right_eye': (314, 115), 'nose': (304, 130), 'mouth_left': (296, 143), 'mouth_right': (313, 142)}}, {'box': [307, 173, 37, 55], 'confidence': 0.8657229542732239, 'keypoints': {'left_eye': (327, 194), 'right_eye': (339, 191), 'nose': (341, 199), 'mouth_left': (334, 215), 'mouth_right': (342, 213)}}] なんだか2か所の顔を検出しているように見えるので、result[1]の方を書き込むようにして実行してみました。 どうやらマイクを持った右手で顔検出しているようです。 あんまり精度が良くないのかもしれません。サンプルの画像なのできちんと検出するものを置いておけばいいのに。 気になった点 コード実行時、以下のエラーとワーニングが出てました。このままでも動いたのでスルーしましたが、いつか気が向いたら調べようと思います。 E tensorflow/core/platform/hadoop/hadoop_file_system.cc:132] HadoopFileSystem load error: libhdfs.so: cannot open shared object file: No such file or directory WARNING:tensorflow:From /home/pi/.local/lib/python3.7/site-packages/tensorflow_core/python/ops/resource_variable_ops.py:1630: calling BaseResourceVariable.__init__ (from tensorflow.python.ops.resource_variable_ops) with constraint is deprecated and will be removed in a future version. Instructions for updating: If using Keras pass *_constraint arguments to layers.

https://atcoder.jp/contests/abc244/tasks/abc244_d 以下の３つの方法で解きます 全列挙で考察 転置数 乱択 1:全列挙で考察 R,G,Bの並び替えは6つしかないです。以下の解説のように図を書いて考えます。 https://atcoder.jp/contests/abc244/editorial/3594 コードは割愛します。(解説に書いてある通りです。 長さがNだとするとハッシュをO(N)でとれるとして 時間計算量はO(N) 空間計算量はO(N!)と極めて大きな数になります 2:転置数 これも上記の解説に書いてある通りです。ある1,2,3,4,...,nという数列が昇順で並んでいた時、異なるi, jを選んでswapするとその転置数は奇数になります。もう一度、適当なswapをすると転置数は偶数になります。 今回の場合、操作の回数が十分なので、偶数回の移動で両者の転置数の偶奇が一致するものに絶対遷移できます。なので、偶奇を判定します。 転置数の実装は典型通りにBITで頑張ります。 時間計算量: O(NlogN) 空間計算量: O(N) 実装(Python) # https://ikatakos.com/pot/programming_algorithm/data_structure/binary_indexed_tree class Bit: def __init__(self, n): self.size = n self.tree = [0] * (n + 1) def sum(self, i): s = 0 while i > 0: s += self.tree[i] i -= i & -i return s def add(self, i, x): while i <= self.size: self.tree[i] += x i += i & -i ############################################ # 転置数を求める # 1,2,3...のリスト。0はダメ版 def transposes(s): n = len(s) bit = Bit(n) ans = 0 for i in range(n): assert s[i] != 0 ans += i - bit.sum(s[i]) bit.add(s[i], 1) return ans a = input().replace("R", "1").replace("G", "2").replace("B", "3").split() b = input().replace("R", "1").replace("G", "2").replace("B", "3").split() a = list(map(int, a)) b = list(map(int, b)) if transposes(a)%2 == transposes(b)%2: print("Yes") else: print("No") 3: 乱択 実はsample-1を読むと、場面: 偶数回のSWAPの回数が残っていても、今、正解と同じ順番になっていればYesになることが分かります。 Nが3と小さいので以下のように考えます。 ある瞬間、あと偶数回swapできるなら2回適当にswapする もし、場面であればYes。Noならもう一度swapを繰り返す 十分な数これを繰り返してもダメならNo ここで、あり得る並び方は1の通り高々6通りのため、十分な数(例えば1000回など)swapすれば、まず、十分に判定できそうです。 このように、適当な試行回数で、全パターンを（ほぼ間違いなく=100%ではない）列挙しきれる場合、乱択は条件の考慮漏れなどの心配がなく、有効な手段です。 今回は愚直に比較をしているので時間計算量は試行回数をkとしてO(kN)ですが、工夫すれば(例えば、異なる数字の数を持っておけば)、O(k)になります。 実装(Python) from random import sample import sys s = input().split() t = input().split() for i in range(100000): if s == t and i%2==0: print("Yes") sys.exit(0) i, j = sample([0,1,2], 2) s[i], s[j] = s[j], s[i] # 適当な2つをswap print("No")

はじめに 本記事では、Microsoftのデータサイエンス初心者向けコースDay3の"Visualizing Proportion"のPython関連について Google Colaboratoryを利用して実際に動かしていく。エラーを特定して修正して動かすという体感を通してのみ学びがあると信じている。 対象読者 　Pythonとデータサイエンスに興味があって、英語が苦手な人（英語が得意な人は、参考文献を直接解いてね。 きのこについて知る。 import pandas as pd import matplotlib.pyplot as plt mushrooms = pd.read_csv('https://raw.githubusercontent.com/microsoft/Data-Science-For-Beginners/main/data/mushrooms.csv') mushrooms.head() きのこに関するデータセットの先頭５行を表示。 毒キノコかどうかの列があるので、極性判定とかできそうな感じ。 ラベルをカテゴリのdtypeとして扱い、食べれるキノコと毒キノコの属性を集計。 cols = mushrooms.select_dtypes(["object"]).columns mushrooms[cols] = mushrooms[cols].astype('category') edibleclass=mushrooms.groupby(['class']).count() edibleclass 食べれるキノコと毒キノコの割合をパイチャートで。 labels=['Edible','Poisonous'] plt.pie(edibleclass['population'],labels=labels,autopct='%.1f %%') plt.title('Edible?') plt.show() キノコはどこで育つのかをドーナッツチャートで。 habitat=mushrooms.groupby(['habitat']).count() habitat labels=['Grasses','Leaves','Meadows','Paths','Urban','Waste','Wood'] plt.pie(habitat['class'], labels=labels, autopct='%1.1f%%', pctdistance=0.85) center_circle = plt.Circle((0, 0), 0.40, fc='white') fig = plt.gcf() fig.gca().add_artist(center_circle) plt.title('Mushroom Habitats') plt.show() ワッフルチャート 「ワッフル」タイプのチャートは、数量を正方形の2次元配列として可視化する別の方法である。このデータセットで、キノコの傘の色の異なる数量を視覚化してみる。これを行うには、PyWaffleというライブラリをインストールし、Matplotlibを使う必要がある。 pip install pywaffle import pandas as pd import matplotlib.pyplot as plt from pywaffle import Waffle capcolor=mushrooms.groupby(['cap-color']).count() capcolor data ={'color': ['brown', 'buff', 'cinnamon', 'green', 'pink', 'purple', 'red', 'white', 'yellow'], 'amount': capcolor['class'] } df = pd.DataFrame(data) fig = plt.figure( FigureClass = Waffle, rows = 100, values = df.amount, labels = list(df.color), figsize = (30,30), colors=["brown", "tan", "maroon", "green", "pink", "purple", "red", "whitesmoke", "yellow"], ) 　おしゃれなWaffleグラフの出来上がり。ボナペティ！ 著者のTWITTERアカウント @keiji_dl ホームページ（pytorch/python/nlp） 参考文献

はじめに (作ろうと思ったきっかけ) 大学のWebサービスのログインパスコードをブラウザで自動生成した結果、非常に長く覚えらづらいパスコードが出来てしまいました。(それ自体は良いこと。) ところが、自分の端末で入力するときはブラウザによる自動入力があるのでいいものの、学内のパソコンでログインするとき、いちいち手打ちする必要がありました。 これがめちゃくちゃ面倒くさい。 てなわけで、自作キーボードの要領でパスコード自動入力機を作れないか構想していました。 Raspberry Pi Pico を選んだ理由 最初は自作キーボードといえばの、Pro Micro を使おうかとも思いましたが、主にマイコンの値段が高いことなどを理由にやめました。 (それのためだけに使うのはコスパ的に見合わないという意味) そこで目をつけたのが Raspberry Pi Pico (以下、ラズパイピコ) です。 ラズパイピコ自体は発売直後に手に入れていたのですが、使用用途が思いつかず、ずっと寝かしていました。 最近になり、ラズパイピコがHIDデバイス(ヒューマンインターフェースデバイス)として認識させることが出来ることを知りました。 すなわち、ラズパイピコをキーボードとしてパソコンに認識させることが出来るのです。 これは使うしか無い!! と、思い製作を始めました。 実行環境の選定 ラズパイピコのプログラミングにどのファームウェア(言語)を用いるかについてです。 今回の候補となったファームウェアについてまとめたいと思います。 ファームウェア 言語系 メリット デメリット C/C++ C/C++ 高速かつ安定 プログラミングが超面倒 MicroPython Python お手軽かつ公式で安心 HIDのライブラリが無い CircuitPython Python お手軽でありながら発展的 公式環境ではない Arduino C系 Arduinoライブラリの資産を活かせる 初期設定がちょい面倒 C/C++ まず、C/C++について、プログラム記述方法が非常に複雑かつ面倒で、挫折してしまいました。 コンパイルするので高速かつ安定して動作するので良いのですが、今回は断念。 MicroPython 続いてMicroPythonについて、当初は公式サポートのあるこの環境を用いようとしていました。 ところが、調べをすすめるうちにHIDのライブラリが無いことがわかり、やむなく使用を断念。 ただ、調査という意味でLチカやシリアル通信、ディスプレイ描画などのチュートリアルは行いました。 CircuitPython そこで代替手段として発見したのが、CircuitPythonです。 軽く調べたところ、このCircuitPythonというのはAdafruit社がMicroPythonを発展させたものとのことです。 そして、こちらのCircuitPython環境下であれば、HIDのライブラリがあります！ よって、今回はこの環境を採用です。 実環境調査 ラズパイピコが2台手元にあったので、以下の写真のようにブレッドボード上にならべて、MicroPythonとCircuitPythonの両方を試してみました。 MicroPythonでSSD1306のチュートリアル CircuitPythonでHIDデバイスの実装 Arduino 最後にArduino環境についてですが、未だ情報が少ないことや、初期設定が少し手間だということから今回はパスしました。 それ以外にも、この間までM5Stackシリーズの製品をArduinoを使ってプログラミングしまくっていたため、ちょっと飽きたというのもあります(笑)。 これとか↓ これとかも↓ コード では実際にコードを書いていきます。 パスコードは仮のものとしてABCDEfghij01234としています。 まずは全体像から。 code.py import usb_hid from time import sleep import digitalio from board import * from adafruit_hid.keyboard import Keyboard from adafruit_hid.keycode import Keycode keyboard = Keyboard(usb_hid.devices) led = digitalio.DigitalInOut(GP25) led.direction = digitalio.Direction.OUTPUT #パスコード(例) : ABCDEfghij01234 pass_list = [4,5,6,7,8,9,10,11,12,13,39,30,31,32,33] big_or_small = [1,1,1,1,1,0,0,0,0,0,0,0,0,0,0] led.value = True sleep(2.5) led.value = False for i in range(len(pass_list)): if big_or_small[i]==1: keyboard.press(225) sleep(0.01) keyboard.send(pass_list[i]) sleep(0.01) keyboard.release(225) led.value = True sleep(5) led.value = False それでは、コードの説明をしていきます。 ライブラリ 導入したライブラリは全部で6つありますが、前半の4つと、後半の2つに分けられます。 前半の4つ import usb_hid from time import sleep import digitalio from board import * これらはCircuitPythonの標準ライブラリです。上から順番に、 USB経由でHID認識させるもの スリープ(待機)させるもの GPIOを制御するもの ボードの定義を参照するもの です。特に設定等なく、記述するだけで利用できます。 後半の2つ from adafruit_hid.keyboard import Keyboard from adafruit_hid.keycode import Keycode 問題はこれらのライブラリについてです。 これらのライブラリは、ラズパイピコ内のライブラリ格納フォルダ(libという名前のフォルダ)に、ライブラリファイルをコピーしておかなければなりません。 CircuitPythonで利用可能なライブラリは、以下のサイトよりダウンロード出来ます。 CircuitPython Libraries 上の画像のように、adafruit-circuitpython-bundle-7.x-mpy-20220304というzipファイルをダウンロードします。 (バージョンは問いません。) フォルダを展開すると、まずはrequirementsフォルダを開いて、インポートしたいライブラリを探してください。 ライブラリのフォルダを開くと、requirements.txtというテキストファイルがあります。 これを参考に、インポートすべきライブラリを確認します。 ひとつ上の階層に戻り、libフォルダを開き、インポートしたいライブラリを探します。 探したら、それをコピーした後、ラズパイピコ本体のlibフォルダにペーストします。 今回の場合は、adafruit_hidフォルダをそのままコピーすればオッケーです。 これで、ライブラリを利用できます。 定義 keyboard = Keyboard(usb_hid.devices) led = digitalio.DigitalInOut(GP25) led.direction = digitalio.Direction.OUTPUT #パスコード(例) : ABCDEfghij01234 pass_list = [4,5,6,7,8,9,10,11,12,13,39,30,31,32,33] big_or_small = [1,1,1,1,1,0,0,0,0,0,0,0,0,0,0] キーボードとLEDの定義 1行目はキーボードの定義 3行目はLEDピンのピン番号の定義 4行目はLEDピンをアウトプットに定義 という流れです。 パスコードをキーコードに変換して変数に保存 adafruit_hid.keyboard.Keyboard Reference ↑のサイトにアクセスし、打ちたいパスコードの文字を1字ずつキーコードに変換していきます。 例えば最初の文字のAを打ちたいとするならば、A = 4 の記述より、pass_listの最初に4を入れます。 ただし、このままだと通常のキーボードと同様に小文字が入力されます。 すなわち、Shiftキーを押して大文字を入力する必要があります。 そこで、big_or_smallという、大文字・小文字を判別するためのリストを用意し、大文字なら1を、小文字なら0を入れることとします。 メインタスク led.value = True sleep(2.5) led.value = False for i in range(len(pass_list)): if big_or_small[i]==1: keyboard.press(225) sleep(0.01) keyboard.send(pass_list[i]) sleep(0.01) keyboard.release(225) led.value = True sleep(5) led.value = False プログラムの流れとしては、 LEDを点灯する (キー操作準備体制の表示) 2.5秒待機 (ホスト側に認識してもらうための待機) LEDを消灯 for文で全文字を入力する 大文字か小文字かの判定 大文字ならShiftキー押下(keyboard.press(225)) リストよりキーコードを取り出し、キー送信 Shiftキー押下を解除 LEDを点灯し、実行終了を知らせる 5秒待機した後、消灯 といった感じです。 実行風景 ラズパイピコをHIDデバイス(キーボード)と認識させて、任意の文字列を打ち込むテスト pic.twitter.com/YbUPXdT0Yr— Koushiro (@Koushiro_A) March 6, 2022 課題点 セキュリティ面が弱い 例えば誰かに無断で使われてしまった場合、モロにパスワードが流出します。 またCircuitPythonの場合、パソコンにストレージとしても認識されるため、ソースコードが丸見えです。 これらの点から、セキュリティ的にはあまりよろしくない状況です。 スタイリッシュじゃない このパスワード打ち込みのためだけにラズパイピコをパソコンとMicroUSBケーブルを使ってつなげるというのは、正直スマートでは無いと感じます。 そもそも、今回のプロダクトにおいてGPIOピンは一切必要ないので、もっと小さくて良いのです。 解決策 そこで、こんな物を発見しました。 こちらは、Adafruit Trinkey QT2040 - RP2040 USB Key with Stemma QT というもので、Adafruit社が製作しているRP2040搭載ボードです。 すなわち、今回のラズパイピコでの開発を、そのまま移植出来る(はず)なのです。 これなら、小型化でき、なおかつパソコンとの接続に必要なケーブルも要りません! セキュリティ面についても解決策が出来たので、 次回、こちらのボードを使ってデバイスを完成させたいと思います。 参考サイト Adafruit CircuitPython API Reference Adafruit HID Library

ラズパイにサーバーを構築し、metabaseでグラフを見る方法についてまとめました。 環境 環境は、下記の通りになります。 MariaDB version　 10.3.31-MariaDB-0+deb10u1 Python version: 3.7.3 OS: Raspbian 10.0 1.1 apacheのインストール apacheをインストールします。 $ sudo apt-get update $ sudo apt-get install apache2 1.2 apacheの確認 下記コマンドを実行し、inetに書かれているIPアドレスを確認します。 $ ifconfig wlan0 Chromeなどのブラウザからhttp://[、inetに書かれているIPアドレス]を入力してください。 成功すると、このようなwebページが出てきます。 1.3 apacheディレクトリ権限の変更 デフォルトの状態でApacheは、/var/www/html/ディレクトリ内にあるファイルを表示します。 下記を入力し、ディレクトの権限を変更します。 $ cd /var/www $ sudo chmod 777 html/ 1.4 PHPのインストール 下記コマンドでPHPをインストールしてください。 $ sudo apt install php libapache2-mod-php -y 2.1 データベース、metabaseの環境構築 データベースは、MariaDBを使用することにしました。MariaDBを見ると、mySQLよりも機能が高く、セキュリティ面が良いらしいので、試しに使ってみます。 MariaDB、セットアップ方法は、こちらを参考にさせていただきました。 2.2 metabaseの設定 metabaseのインストール サイトからmetabase.jarをダウンロードし、作成したフォルダに配置します。私は、windows端末からダウンロードし、winSCPを経由してラズパイのディレクトリに置きましたおきました。 $ sudo mkdir -p /apps/java $ cd /apps $ sudo chmod 777 java $ pwd /apps/java/metabase.jar 2.3 OSのユーザー設定 metabase実行用のユーザを定義して権限等の設定をします。 $ sudo groupadd -r metabase $ sudo useradd -r -s /bin/false -g metabase metabase $ sudo chown -R metabase:metabase /apps/java $ sudo touch /var/log/metabase.log $ sudo chown metabase:metabase /var/log/metabase.log $ sudo touch /etc/default/metabase $ sudo chmod 640 /etc/default/metabase 2.4 サービス作成 metabaseをサービスとして登録します。 $ sudo touch /etc/systemd/system/metabase.service $ sudo nano /etc/systemd/system/metabase.service metabase.service [Unit] Description=Metabase server After=syslog.target After=network.target [Service] WorkingDirectory=/apps/java/ ExecStart=/usr/bin/java -jar /apps/java/metabase.jar EnvironmentFile=/etc/default/metabase User=metabase Type=simple StandardOutput=syslog StandardError=syslog SyslogIdentifier=metabase SuccessExitStatus=143 TimeoutStopSec=120 Restart=always [Install] WantedBy=multi-user.target 2.5 syslogファイル作成とMetabaseサービス登録 $ sudo touch /etc/rsyslog.d/metabase.conf $ sudo nano /etc/rsyslog.d/metabase.conf 以下をmetabase.conf内に以下を追記します。 metabase.conf if $programname == 'metabase' then /var/log/metabase.log & stop syslogを再起動し、サービスを設定します。 $ sudo systemctl daemon-reload $ sudo systemctl start metabase.service $ sudo systemctl enable metabase.service 2.6 接続確認 下記コマンドを実行し、chromeでhttp://:3000でアクセスできるかどうかを確認します。 エラーが発生して場合、2.7を見てください。 $ java -jar /apps/java/metabase.jar 2.7 ポート重複によるエラー回避 私の環境では、3000のポートが使用していたため下記のようなエラーが出てしまいました。 WARNING: sun.reflect.Reflection.getCallerClass is not supported. This will impact performance. 2022-03-16 14:13:54,552 INFO metabase.util :: Maximum memory available to JVM: 948.0 MB 2022-03-16 14:14:32,072 INFO util.encryption :: Saved credentials encryption is DISABLED for this Metabase instance. ? For more information, see https://metabase.com/docs/latest/operations-guide/encrypting-database-details-at-rest.html 2022-03-16 14:14:44,643 INFO driver.impl :: Registered abstract driver :sql ? ⮦ Load driver :sql took 2.3 s 2022-03-16 14:14:44,716 INFO driver.impl :: Registered abstract driver :sql-jdbc (parents: [:sql]) ? Load driver :sql-jdbc took 2.4 s 2022-03-16 14:14:44,736 INFO driver.impl :: Registered driver :h2 (parents: [:sql-jdbc]) ? 2022-03-16 14:14:44,797 INFO driver.impl :: Registered driver :mysql (parents: [:sql-jdbc]) ? 2022-03-16 14:14:44,852 INFO driver.impl :: Registered driver :postgres (parents: [:sql-jdbc]) ? 2022-03-16 14:14:53,939 INFO metabase.core :: Metabase v0.42.2 (d6ff494 release-x.42.x) Copyright © 2022 Metabase, Inc. Metabase Enterprise Edition extensions are NOT PRESENT. 2022-03-16 14:14:53,968 INFO metabase.core :: Starting Metabase in STANDALONE mode 2022-03-16 14:14:54,103 INFO metabase.server :: Launching Embedded Jetty Webserver with config: {:port 3000} 2022-03-16 14:14:54,183 ERROR metabase.core :: Metabase Initialization FAILED java.io.IOException: Failed to bind to 0.0.0.0/0.0.0.0:3000 metabase.comを確認すると、3000番が使われてたら、jarを走らす前にMB_JETTY_PORTを使って別のところ使ってね。と記載があります。なので、下記のように変更し、走らせます。chromeでは、http://:3100でアクセスします。 $ java -DMB_JETTY_PORT=3100 -jar /apps/java/metabase.jar 3.1 データベースの追加 metabaseを起動すると、初期設定をするようなっています。それほど難しい項目はないので、この部分は割愛します。 metabaseへmySQLのデータを登録します。mySQL内でのデータベースの作成方法については、以前の記事に記述しています。 下記の記事を参照してください。 pythonでcsvのデータをデータベースへ書き込む。 3.2 データベースの登録 右上にある歯車マークをクリックし、管理者設定 -> データベースを追加する を選択する。 下記の項目に必要な情報を入力する。 データベースのタイプ: MySQL 表示名: 適当に入力(metabase内のタイトルになります) ホスト名: localhost ポート: 3306 データベース名: 適当に入力(設定したものを入力) ユーザー名: root(設定したものを入力) パスワード: xxx(設定したものを入力) ※ポート番号について mySQL(MariaDB)にアクセスし、下記コマンドを入力すると、mySQLのポート番号を確認することができます。 MariaDB[(none)]> show variables like 'port'; 3.3 グラフ作成 適当に操作して、グラフ、ダッシュボートを作成します。下記は、作成例です。 Link 下記サイトを参照させていただきました。 Raspberry PiにMetabaseをインストール mariadb mysql ポート番号を調べる

ABC001 メモ A - 積雪深差 問題文通り、答えは$H_1-H_2$。 001A.py H_1 = int(input()) H_2 = int(input()) ans = H_1 - H_2 print(ans) B - 視程の通報 入力は単位がメートル(m)で与えられるが、判定条件は単位がkmで与えられているため、少々複雑になる。 今回は判定条件をmに統一して計算した。 001B.py m = int(input()) if m < 100: VV = '00' elif m <= 5000: if len(str(m)) == 4: VV = str(m)[:2] else: VV = '0' + str(m)[0] elif m <= 30000: VV = str(m//1000+50) elif m <= 70000: VV = str((m//1000-30)//5+80) else: VV = '89' print(VV) C - 風力観測 問題文通り実装した。場合分けが煩雑なのと、小数を扱う場合には誤差に注意が必要。小数の扱いについては、整数で扱えるように条件式を加工した。 具体的には、風速は分速のまま扱えるようにした。今回の場合、風速$0.2$m/s以下という条件は、$0.25$m/s未満と同値なため、これを60倍して$15$m/min未満、というようにすることで、小数の計算を避けた。 公式の解説にあるように、場合分けを手動でたくさん書くのはあまり望ましくないのだが… 001C.py Deg, Dis = map(int, input().split()) if Deg <= 112: Dir = 'N' elif Deg <= 337: Dir = 'NNE' elif Deg <= 562: Dir = 'NE' elif Deg <= 787: Dir = 'ENE' elif Deg <= 1012: Dir = 'E' elif Deg <= 1237: Dir = 'ESE' elif Deg <= 1462: Dir = 'SE' elif Deg <= 1687: Dir = 'SSE' elif Deg <= 1912: Dir = 'S' elif Deg <= 2137: Dir = 'SSW' elif Deg <= 2362: Dir = 'SW' elif Deg <= 2587: Dir = 'WSW' elif Deg <= 2812: Dir = 'W' elif Deg <= 3037: Dir = 'WNW' elif Deg <= 3262: Dir = 'NW' elif Deg <= 3487: Dir = 'NNW' else: Dir = 'N' if Dis < 15: Dir = 'C' W = 0 elif Dis < 93: W = 1 elif Dis < 201: W = 2 elif Dis < 327: W = 3 elif Dis < 477: W = 4 elif Dis < 645: W = 5 elif Dis < 831: W = 6 elif Dis < 1029: W = 7 elif Dis < 1245: W = 8 elif Dis < 1467: W = 9 elif Dis < 1707: W = 10 elif Dis < 1959: W = 11 else: W = 12 print(Dir, W) D - 感雨時刻の整理 時間の区間を扱う問題だが、今回は最小の単位が分であり、最長でも1440分なので、配列で保持しても十分高速に計算できそう。 入力を00:00からの分に計算し直し、配列に保持する。この時、区間全てに対して計算を行う必要はなく、始まった時と終わった時に対して計算を行う、いわゆる「いもす法」で扱えば十分である。 計算量は$N$が大きい時でも$O(N)$程度であるため、十分高速。 001D.py def adj(num): h = str(num//60) m = str(num%60) time = h.zfill(2) + m.zfill(2) return time N = int(input()) lst = [0]*1442 for i in range(N): S, E = input().split('-') S = 60*int(S[:2]) + int(S[2:]) E = 60*int(E[:2]) + int(E[2:]) if S%10 < 5: S -= S%10 else: S -= S%10 S += 5 if 0 < E%10 < 6: E -= E%10 E += 5 else: E += 5 E -= E%10 lst[S] += 1 lst[E+1] -= 1 tmp = 0 S = -1 E = -1 for i in range(1442): tmp += lst[i] if tmp and S == -1: S = i if tmp == 0 and S != -1: E = i-1 S = adj(S) E = adj(E) print(S + '-' + E) S = -1

NGINX Unitについて PythonをNGINXで使用したいときに中継サーバーとして使用 Python,PHP,Go,Perl,Ruby,Node.js,Javaなどに対応 同じ言語の複数バージョンを同時に運用可能 サーバー停止せずに設定変更可能 公式ドキュメント https://unit.nginx.org/installation/ インストール $ sudo vi /etc/yum.repos.d/unit.repo [unit] name=unit repo baseurl=https://packages.nginx.org/unit/amzn2/$releasever/$basearch/ gpgcheck=0 enabled=1 $ sudo yum install unit $ sudo yum install unit-devel unit-go unit-jsc8 unit-perl \ unit-php unit-python27 unit-python37 自動再起動設定 sudo systemctl enable unit 再起動 sudo systemctl restart unit ログファイル /var/log/unit/unit.log

初めに 本記事は時系列解析手法の一つであるSARIMAモデルをpythonのstatsmodelsで実装するまでの流れをまとめたものになります。 学んだことをアウトプットすることを目的として作成しました。 目次 １. 定常性とは ２. ホワイトノイズ ３. SARIMAモデルとは 1.ARモデル 2.MAモデル 3.ARMAモデル 4.ARIMAモデル SARIMAモデル ４. SARIMAモデルの実装 1. 定常性とは 時系列解析を語る上で欠かせないものの一つに定常性があります。定常性とは、ある時系列データの値$Y$が、どの時点$t$においても平均、分散、自己共分散が一定であるという性質のことです。これは弱定常性とも言います。 簡単にいうと、各時点での値がある値を中心にばらついているということです。 2.ホワイトノイズ ホワイトノイズとは、平均が0、全ての時点で自己共分散が0である定常性を持った系列のことをいいます。 時系列分析では、モデルが説明できなかった誤差項がホワイトノイズの性質を持っていると仮定しています。 3.　SARIMAモデルとは SARIMAモデルとは、自己回帰モデル(ARモデル)と移動平均モデル(MAモデル)と和分モデル(I)を組み合わせた自己回帰和分移動平均モデル(ARIMAモデル)に季節的な周期変動を取り入れたモデルです。 SARIMAモデルの理解には上記のモデルの理解が必要不可欠なので、まずはそれらのモデルについて説明します。 1.　ARモデル ARモデルは、自身の過去の値で回帰を行うことで現在の値を予測するというモデルになります。式で表すと次のようになります。 $$ Y_t = c + φ_1Y_{t-1} + ε_t $$ ここで$Y_t$は時点$t$での値、$c$は定数項、$φ$は係数、　$ε_t$はホワイトノイズに従う誤差項を表しています。 上の式は、一つ前の時点の値から現在の値を予測しています。このようなARモデルをAR(1)過程といいます。 2.　MAモデル MAモデルは、過去の誤差の積み重ねにより現在の値を求めるというモデルになります。式で表すと次のようになります。 $$ Y_t = c + φ_1ε_{t-1} + ε_t $$ ARモデルと比較すると、係数$φ$が過去の値$Y_{t-1}$ではなく、過去の誤差項$ε_{t-1}$にかかっていることがわかります。この式はMA(1)過程ともいいます。 3. ARMAモデル ARMAモデルはARモデルとMAモデルを組み合わせたモデルになります。式で表すと次のようになります。 ARMA(1,1)過程ともいいます。 $$ Y_t = c + φ_1Y_{t-1} + ε_t + θ_1ε_{t-1} $$ 4. ARIMAモデル 時系列データによっては、その系列が定常性を満たさず、階差をとった系列が定常性を満たすものもあります。そのような時系列データについてはARIMAモデルを使用することができます。式で表すと次にようになります。 $$ ΔY_t = c + φ_1ΔY_{t-1} + ε_t + θ_1ε_{t-1} $$ ここでの$ΔY_t$は$Y_t - Y_{t-1}$のことです。このようなARIMAモデルをARIMA(1,1,1)過程といいます。(AR = 1、 I = 1、 MA =1 ) 5. SARIMAXモデル この記事では扱いませんが、SARIMAモデルに説明変数を加えたSARIMAXモデルも存在します。 4. SARIMAモデルの実装 pythonでは、statsmodelsというライブラリを使うことでSARIMAモデルを実装できます。 まず今回使用するライブラリをインポートし、データを読み込みます。この記事で使用するデータは気象庁ホームページからダウンロードした2017年1月~2022年3月15日までの東京都の気温データです。気温データ以外にも降水量や日照時間といったデータも含まれていますが、今回は使用しません。 import pandas as pd import matplotlib.pyplot as plt %matplotlib inline import japanize_matplotlib import statsmodels.api as sm from statsmodels.tsa.seasonal import seasonal_decompose df_tokyo = pd.read_csv('tokyo.csv',index_col="年月日", parse_dates=True) 気温についてグラフを作成し可視化します。 from matplotlib.pylab import rcParams rcParams['figure.figsize'] = 15, 6 y = df_tokyo["平均気温(℃)"] y.plot() 日本には季節があるので、グラフは周期的なものになります。 次に自己相関係数と偏自己相関係数を可視化します。 fig_1 = sm.graphics.tsa.plot_acf(y, lags=400) fig_2 = sm.graphics.tsa.plot_pacf(y, lags=50) 上のグラフは自己相関係数のコレログラムで、下の図は偏自己相関係数のコレログラムです。この図からラグ10で偏自己相関係数が0になることがわかります。 次に波形分解をしてみます。 result = seasonal_decompose(y, period=365, two_sided=True) result.plot() seasonal_decomposeの引数periodは、自分で周期を指定する必要があります。先ほどのコレログラムからも365日周期であることが考えられるので　period=365　と指定していいでしょう。 出力結果の一番上のグラフは元データ、二番目のグラフは周期性を除いたデータ、三番目は周期性を表すグラフ(今回は３６５と指定しているので365周期のグラフになっています)、四番目は残差のグラフになります。 波形分解をするとデータの周期性を除いた傾向がわかります。上の例で言うと、傾向のグラフを見ることで、2017年あたりはそれ以降の年に比べて全体的に気温が低かったことがわかります。 次にSARIMAモデルを実装していきます。 （と、言いたいところだったのですが、365周期で学習を行うと計算に時間がかかり過ぎてしまうので、月ごとにデータを取り直して実装していきたいと思います。） データの後ろ12ヶ月分のデータをテストデータとし、残りのデータで学習を行なっていきます。 df_tokyo_month = df_tokyo.resample("M").mean() y = df_tokyo_month["平均気温(℃)"] train = y[:-12] test = y[-12:] sarima_model = sm.tsa.SARIMAX(train, order=(1,1,1),seasonal_order=(1,1,1,12)) result = sarima_model.fit() result.summary() sm.tsa.SARIMAXでSARIMAXモデルを実装できますが、説明変数Xを指定しなければSARIMAモデルとして使用できます。 引数orderには左からAR、I、MAのそれぞれの差分の値を指定します。引数seasonal_orderにはAR、I、MAに加え周期性として適応する値を指定します。今回の例では、日ごとのデータを月ごとに取り直したので12を指定しています。これらのパラメータを調整してAICが低くなるようなモデルを作成します。 最後に予測を行います。 bestPred = result.predict('2021-4-30', '2022-03-31') plt.plot(df_tokyo_month["平均気温(℃)"]) plt.plot(bestPred, "red") 結果を見ると、おおよその予測はできていると考えていいと思います。パラメータをいじればもう少し精度が良くなるかもしれません。

以下は、複数スレッドから同時アクセスされたときに前回のアクセス時間よりも２秒以上経過していなければ待機して、経過後にアクセス時間を記録して、この値をreturnする処理 import time import threading lock = threading.Lock() access_time = time.time() def task(): while True: print(concurrent_access()) pass def concurrent_access(): """複数スレッドに同時にアクセスされるかもしれない処理""" # 前回のアクセス時間より2秒以上経過して "いる" # access_timeを現在時刻にしてから # access_timeをreturn # 前回のアクセス時間より2秒以上経過して 'いない" # 2秒経過するまで待機 # 2病経過したら上と同じ処理 global access_time try: # 現在のスレッド名を表示 print(threading.current_thread().name) lock.acquire() while True: if access_time + 2 < time.time(): access_time = time.time() return access_time finally: lock.release() ths = [] for i in range(3): th = threading.Thread(target=task, name=f'thread{i}') th.start() ths.append(th) for th in ths: th.join() 実行結果 thread0 thread1 thread2 1648015757.6682205 thread0 1648015759.6691394 thread1 1648015761.6700594 thread2 1648015763.6702552 thread0 1648015765.671093 thread1 1648015767.6731405 thread2

以下は、複数スレッドから同時アクセスされたときに前回のアクセス時間よりも２秒以上経過していなければ待機して、経過後にアクセス時間を記録して、この値をreturnする処理 import time import threading lock = threading.Lock() access_time = time.time() def task(): while True: print(concurrent_access()) pass def concurrent_access(): """複数スレッドに同時にアクセスされるかもしれない処理""" # 前回のアクセス時間より2秒以上経過して "いる" # access_timeを現在時刻にしてから # access_timeをreturn # 前回のアクセス時間より2秒以上経過して 'いない" # 2秒経過するまで待機 # 2病経過したら上と同じ処理 global access_time try: # 現在のスレッド名を表示 print(threading.current_thread().name) lock.acquire() while True: if access_time + 2 < time.time(): access_time = time.time() return access_time finally: lock.release() ths = [] for i in range(3): th = threading.Thread(target=task, name=f'thread{i}') th.start() ths.append(th) for th in ths: th.join() 実行結果 thread0 thread1 thread2 1648015757.6682205 thread0 1648015759.6691394 thread1 1648015761.6700594 thread2 1648015763.6702552 thread0 1648015765.671093 thread1 1648015767.6731405 thread2

前置き 2chで有名なコピペの「初カキコ…ども…」がありますね。 今週のオモコロの記事【第1回】初カキコ…ども…選手権の影響で初カキコが再び話題になっています。 このコピペを感情分析したら、ポジティブ・ネガティブになるか分析してみました。 初カキコの全文は以下になります。 初カキコ…ども… 俺みたいな中3でグロ見てる腐れ野郎、他に、いますかっていねーか、はは 今日のクラスの会話 あの流行りの曲かっこいい　とか　あの服ほしい　とか ま、それが普通ですわな かたや俺は電子の砂漠で死体を見て、呟くんすわ it’a true wolrd．狂ってる？それ、誉め言葉ね。 好きな音楽　eminem 尊敬する人間　アドルフ・ヒトラー（虐殺行為はNO） なんつってる間に4時っすよ（笑）　あ～あ、義務教育の辛いとこね、これ 実行結果 感情値が-0.488227のため、ネガティブという結果になりました。 ポジティブ・ネガティブの単語に色付け 全体的にネガティブが多い… 実行結果（全文表示） 一文ずつの感情値を出力しました。 全体的に感情値がネガティブになっていることが分かると思います。 全文の感情値の平均を求めると-0.488227となります。 -0.28484 初カキコ…ども… -0.98782 俺みたいな中3でグロ見てる腐れ野郎、他に、いますかっていねーか、はは -0.20401 今日のクラスの会話 -0.72328 あの流行りの曲かっこいい　とか　あの服ほしい　とか 0.0 ま、それが普通ですわな -0.77961 かたや俺は電子の砂漠で死体を見て、呟くんすわ -0.98429 it’a true wolrd．狂ってる？それ、誉め言葉ね。 0.12064 好きな音楽　eminem -0.63717 尊敬する人間　アドルフ・ヒトラー（虐殺行為はNO） -0.40189 なんつってる間に4時っすよ（笑）　あ～あ、義務教育の辛いとこね、これ Ave = -0.488227 感情分析とは？ 極性辞書 東工大が公開している単語感情極性対応表のこと。 単語に（ネガティブ）-1~（ポジティブ）1と感情値を振ってある表となります。 ポジティブな単語（一部抜粋） ネガティブな単語（一部抜粋） ポジティブ・ネガティブの出し方 「今日は良い天気だ」の感情値を求めると、0.47101という結果が出力され、ポジティブな文章と分析できます。 極性辞書より「今日」「良い」「天気」の単語の感情値を求めます。 今日：0.172375 良い：0.999995 天気：0.24065 感情値を全て足し合わせた後、単語の総数で割ることで0.47101となります。 (0.172375 + 0.999995 + 0.24065) / 3 = 0.47101 なぜポジティブ・ネガティブの結果になった？ 初カキコの文章がネガティブになった理由は、ネガティブの単語数が全体の9割を占めているため。 極性辞書の単語数は全55,125語 ポジティブの割合（0.0＜感情値≦1.0）:0.093(4882語) ネガティブの割合（-1.0≦感情値＜0.0）:0.907(47776語) ニュートラルの割合（感情値=1.0）:0.000(13語) コード説明（GoogleColabで実行） janomeをインストール JanomeはPythonの形態素解析エンジン 日本語のテキストを形態素ごとに分割して品詞を判定したり分かち書き（単語に分割）したりすることができる !pip install janome 極性辞書をアップロード 極性辞書（pn_ja.dic）を単語感情極性対応表からダウンロード GoogleColabにアップロードする from google.colab import files files.upload() テキスト（初カキコ）の文から空白などを削除 テキストを適当な変数に入れる original_text = """初カキコ…ども… 俺みたいな中3でグロ見てる腐れ野郎、他に、いますかっていねーか、はは 今日のクラスの会話 あの流行りの曲かっこいい　とか　あの服ほしい　とか ま、それが普通ですわな かたや俺は電子の砂漠で死体を見て、呟くんすわ it’a true wolrd．狂ってる？それ、誉め言葉ね。 好きな音楽　eminem 尊敬する人間　アドルフ・ヒトラー（虐殺行為はNO） なんつってる間に4時っすよ（笑）　あ～あ、義務教育の辛いとこね、これ""" 改行ごとに要素をリストに入れる 空白（\u3000）を削除 result = original_text.split('\n') result = [i for i in result if not i in ("","\u3000")] 感情値を求める コード sentiment_dic = {} with open('pn_ja.dic', 'r') as f: lines = f.readlines() for line in lines: line_components = line.split(':') sentiment_dic[line_components[0]] = line_components[3] from janome.tokenizer import Tokenizer keyword = [] def sentiment_analyse(text): t = Tokenizer() tokens = t.tokenize(text) sentiment_val = 0 word_cnt = 1e-6 #割り算で使うので0に近い値 for token in tokens: #単語の基本形を取り出し word = token.surface if word in sentiment_dic: keyword.append([word, float(sentiment_dic[word])]) sentiment_val = sentiment_val + float(sentiment_dic[word]) word_cnt += 1 return round(sentiment_val/word_cnt, 5) ave_list = [] cnt = 0 result_list = [] for i in result: print(sentiment_analyse(i),i) ave_list.append(sentiment_analyse(i)) result_list.append([i,sentiment_analyse(i)]) print("Ave = ",sum(ave_list)/len(ave_list)) 出力結果 -0.28484 初カキコ…ども… -0.98782 俺みたいな中3でグロ見てる腐れ野郎、他に、いますかっていねーか、はは -0.20401 今日のクラスの会話 -0.72328 あの流行りの曲かっこいい　とか　あの服ほしい　とか 0.0 ま、それが普通ですわな -0.77961 かたや俺は電子の砂漠で死体を見て、呟くんすわ -0.98429 it’a true wolrd．狂ってる？それ、誉め言葉ね。 0.12064 好きな音楽　eminem -0.63717 尊敬する人間　アドルフ・ヒトラー（虐殺行為はNO） -0.40189 なんつってる間に4時っすよ（笑）　あ～あ、義務教育の辛いとこね、これ Ave = -0.488227 ポジティブ・ネガティブの単語のリストを整理 重複している単語が多いため、重複単語を1つのみに変更 keyword = list(map(list, set(map(tuple, keyword)))) ポジティブ・ネガティブの単語に色付け ポジティブ単語は赤、ネガティブ単語は青で色付け import re color_dic = {'red':'\033[31m', 'blue':'\033[34m', 'end':'\033[0m'} def print_hl(text, keyword, color="yellow"): for kw in keyword: if float(kw[1]) < 0: bef = kw[0] aft = color_dic["blue"] + kw[0] + color_dic["end"] text = re.sub(bef, aft, text) elif float(kw[1]) > 0: bef = kw[0] aft = color_dic["red"] + kw[0] + color_dic["end"] text = re.sub(bef, aft, text) print(text) for i in result: print_hl(i, keyword) おまけ ユーザーローカル感情認識AIで感情分析してみました。 怒りの感情が多く含んでいるようです この文章の感情が怒り・好き・恐れが含まれているのは違和感ありますね、中二病の感情はわかりづらい 参考にしたサイト 【第1回】初カキコ…ども…選手権 単語感情極性対応表 ユーザーローカル感情認識AI 3. Pythonによる自然言語処理　5-3. 日本語文の感情値分析［単語感情極性値対応表］ 【Python】print出力で特定のキーワードの文字色や背景色を変える Python, Janomeで日本語の形態素解析、分かち書き（単語分割）

スクレイピングってなに？ スクレイピングとは、ウェブサイトから情報を抽出することです。 サイト上のデータ(商品、価格、天気などのデータ)を取ってくることが多いです。 pythonで書いてみよう 必要なもの(pip) ・requests ・BeautifulSoup ・調べたいサイトのURL 流れ requestsで取得 res = requests.get("URL") BeautifulSoupでサイトデータを抽出 soup = BeautifulSoup(res.text, "html.parser") 欲しい部分を決定 ex)クラス名:product_lists aaa の ul の中の liを取得 ※商品一覧などの構造はだいたいこんな感じ found_part = soup.find_all("ul", class_='product_lists aaa') for ul_tag in found_part: for li in ul_tag.find_all('li'): print(li.find('span').text) サンプルコード # coding: utf-8 import requests from bs4 import BeautifulSoup def get_soup(url): res = requests.get(url) soup = BeautifulSoup(res.text, "html.parser") tag_obj = soup.title print(tag_obj) if __name__ == '__main__': get_soup("ここにリンクのURL") ここ注意 サイトによってはスクレイピングを禁止しているものもあるため、確認する必要がある。

すべてのスレッドでevent.waitしている場合、各スレッドがロックされるためctrl+cが受け取れない。 そのため、event.waitを使わないスレッドを必ず１つ用意しておく。このスレッドでctrl+cを受け取り終了フラグをTrueにして子スレッドの条件に終了フラグを配置していおけばループから抜けることができる。 import time import threading e_task1 = threading.Event() e_task2 = threading.Event() f_sys_exit = False # =============================== # 制御不可能なchildスレッドの処理 # =============================== def can_not_control_task1(): while True: print(threading.current_thread().name, ' is waitting.') e_task1.wait() print(threading.current_thread().name, ' end.') def cant_not_control_task2(): while True: print(threading.current_thread().name, ' is waitting.') e_task2.wait() print(threading.current_thread().name, ' end.') # =============================== # 制御可能なchildスレッドの処理 # =============================== def can_control_task1(): while True: if f_sys_exit: break print(threading.current_thread().name, ' is waitting.') e_task1.wait() print(threading.current_thread().name, ' end.') def can_control_task2(): while True: if f_sys_exit: break print(threading.current_thread().name, ' is waitting.') e_task2.wait() print(threading.current_thread().name, ' end.') # =============================== # childスレッドの実行処理 # =============================== def run_can_not_control_child_thread(): """コントール不可能""" ths = [] for task in [can_not_control_task1, cant_not_control_task2]: th = threading.Thread(target=task, name=task.__name__, daemon=False) th.start() ths.append(th) for th in ths: th.join() def run_can_control_child_thread(): """コントール可能""" ths = [] for task in [can_control_task1, can_control_task2]: th = threading.Thread(target=task, name=task.__name__, daemon=False) th.start() ths.append(th) for th in ths: th.join() # =============================== # メインスレッドの実行処理 # =============================== def can_not_control(): """コントール不可能""" global f_sys_exit threading.Thread(target=run_can_not_control_child_thread, daemon=True).start() while True: try: time.sleep(1) print('parent thread is running.') except KeyboardInterrupt: print('pressed to Ctrl + C keys.') break def can_control(): """コントール可能""" global f_sys_exit threading.Thread(target=run_can_control_child_thread, daemon=True).start() while True: try: time.sleep(1) print('parent thread is running.') except KeyboardInterrupt: print('pressed to Ctrl + C keys.') # sys_exitフラグをTrueにしてから　 f_sys_exit = True # すべてのwait中のスレッドをsetする set_all_event() break # =============================== # イベント処理 # =============================== def set_all_event(): """すべてのイベントをset""" global e_task1, e_task2 e_task1.set() e_task2.set() if __name__ == '__main__': # can_not_control() can_control()

TL;DR 研究遂行経費報告書がまさかの手入力だったので、pythonで自動化した はじめに 今年度ももう終わり、様々な書類仕事をしなければならないです。博士課程にいる学生の中では、学振のしちめんどくさい色々な手続きをしている人もいるでしょう。学振には研究遂行経費というものがあり、所得の一部を研究経費として計上できる制度があります。これを利用すると所得税を節税できるので、この制度を利用している学振研究者は多いと思います。ExcelだとかGoogle spread sheetだとかを使って購入品目を管理している人もいるでしょう。私もその一人です。見返したら大体3~40品目ぐらい購入していました。 さて、そんな研究遂行経費を報告するWebフォームがこちら。 　 (　ﾟдﾟ) ＿(＿つ/￣￣￣/＿ 　　＼/　　　/ 　　　 ￣￣￣ 　 ( ﾟдﾟ ) ＿(＿つ/￣￣￣/＿ 　　＼/　　　/ 正気か？ まぁ色々とひどい。疑念だらけである。 手打ち？csvインポーターもない？ほとんどの人はExcelで管理してるって。品名は詳細を書けって言ってるくせに30文字以上でエラー出すのはなぜ？人によっては100品目超える場合もあるんやぞ。品名を適宜修正しながらの手入力じゃ1日かかるわい。予め規定のExcelファイルを配布しておけば互いに楽ではないのか？ etc... そして何より、 「こんなことに研究者の時間を奪っておいて「日本の研究力が～」とか言ってんじゃないだろうな？」 ...という腹立たしさが目覚めたとともに、他の人も絶対同じことを思い悲しんでいるはず。「せめて入力だけでも...」と思いpythonを使って自動化しました。今回の記事はその報告です。 ↓　成果物 日本学術振興会さん、学術を振興させたいなら来年度はまともなシステムに変えてください... 中身の説明 とはいえ、やっていることはpandasによるデータ整形とseleniumによる自動入力なのでn番煎じではあります。 githubにコードと専用のエクセルファイルを上げるので良かったら使ってあげてください。 python触れない人のためにもexe化も検討しています（やれるか怪しい）。 事前準備 まずgoogle chromeをインストールし、chromeのバージョンに合ったchromedriverをここからダウンロード。 私の環境ではChromeDriver 99.0.4844.51をダウンロードしました。 OSによってダウンロードするものが違うので注意。windowsはwin32, macはmac64 or mac64_m1を入れる。 解凍し、pythonを記述するファイルと同じディレクトリに中身を入れておく。これで準備はok スクリプトの記述 必要なもろもろをimport。足りないものはpip or condaでインストールしてください。 import.py import pandas as pd from selenium import webdriver from time import sleep import time import datetime from selenium.webdriver.support.select import Select そしてコツコツ入力しておいたExcel fileの読み込み read_excel.py df = pd.read_excel("研究遂行経費の支出報告書用入力シート.xlsx").dropna() Excel fileの中身はこんな感じになっています。 納入日 品名 単価 税抜・税込 個数 経費項目 領収書の有無 0 2021-04-21 00:00:00 foo 10000 税込 1 学術調査にかかる経費 あり 1 2021-05-22 00:00:00 bar 30000 税込 1 学術調査にかかる経費 あり 先のことを考え、列名を変更。ついでに各品目でかかった総額も出しておく。 rename_column.py df.columns = ["date", "item", "value", "is_tax","number", "expenditure", "is_receipt"] def add_tax(x): if x.is_tax =="税抜": return x.value*1.1 else: return x.value df['value'] = df.apply(lambda x: add_tax(x), axis = 1)*df["number"] ここからは、Web側の入力に合わせてデータを前処理する。 まず、プルダウンメニューについて。 こちらの記事によると、seleniumのselect.select_by_value()で指定ができるらしい。 引数はintであり、要はプルダウンメニュー内の順番。これに合わせてデータ側を改変していく。 「領収書の有無」はラジオボタンなので、それぞれのボタンのidを探してそのidを割り当てればok。 forum_impout.py df['date'] = df['date'].dt.strftime("%Y-%m") def convert_date_to_flag(x): if x.date == "2021-04": return 1 elif x.date == "2021-05": return 2 elif x.date == "2021-06": return 3 elif x.date == "2021-07": return 4 elif x.date == "2021-08": return 5 elif x.date == "2021-09": return 6 elif x.date == "2021-10": return 7 elif x.date == "2021-11": return 8 elif x.date == "2021-12": return 9 elif x.date == "2022-01": return 10 elif x.date == "2022-02": return 11 elif x.date == "2022-03": return 12 def is_receipt(x): if x.is_receipt == "あり": return "etr_needOrNotReceipt_01" else: return "etr_needOrNotReceipt_02" def assign_expenditure(x): if x.expenditure == "学会関係経費": return 1 elif x.expenditure == "各種研究集会等への参加費": return 2 elif x.expenditure =="学術調査にかかる経費": return 3 elif x.expenditure == "自宅での研究に必要な経費": return 4 elif x.expenditure == "所属・関連機関への交通費": return 5 df["date"] = df.apply(lambda x: convert_date_to_flag(x), axis = 1) df["is_receipt"] = df.apply(lambda x: is_receipt(x), axis = 1) df["expenditure"] = df.apply(lambda x: assign_expenditure(x), axis = 1) df_l = [list(row) for row in df.itertuples()] df = df[["date", "item", "expenditure", "is_receipt", "value"]] df["value"] = pd.Series(df["value"], dtype = 'int64') 変更後はこんな感じ。 date item expenditure is_receipt value 0 1 foo 3 etr_needOrNotReceipt_01 10000 1 2 bar 3 etr_needOrNotReceipt_01 30000 これで入力用のデータは作れた。あとはseleniumで自動入力させていけばいい。 seleniu.py driver = webdriver.Chrome(executable_path='./chromedriver.exe') JSPS_URL = "https://area31.smp.ne.jp/area/servlet/area.MyPageBundle?MyPageID=629871a6_lalj9lenala0oare4" JSPS_ID = "JSPSのID" JSPS_PASS = "JSPSのpass" J_SYSTEM_ID = "J_SYSTEMのID" J_SYSTEM_PASS = "J_SYSTEMのpass" driver.get(JSPS_URL) sleep(1) driver.find_element_by_name("SMPID").send_keys(JSPS_ID) driver.find_element_by_name("SMPPASSWORD").send_keys(JSPS_PASS) driver.find_element_by_class_name("btn").click() driver.find_element_by_link_text('⑤ 報告書等提出システム（Jシステム）').click() sleep(1) driver.switch_to.window(driver.window_handles[1]) driver.find_element_by_name("login_id").send_keys(J_SYSTEM_ID) driver.find_element_by_name("login_psw").send_keys(J_SYSTEM_PASS) driver.find_element_by_name("lg").click() sleep(1) driver.find_element_by_id("li_of_F003").click() sleep(1) def enter_forum(row): dropdown_expense = driver.find_element_by_id('etr_group') select_expense = Select(dropdown_expense) select_expense.select_by_index(row.expenditure) dropdown_date = driver.find_element_by_id('etr_dateOfComp') select_date = Select(dropdown_date) select_date.select_by_index(row.date) driver.find_element_by_id(str(row.is_receipt)).click() driver.find_element_by_name("etr_itemName").send_keys(row.item) driver.find_element_by_name("etr_billAmount").send_keys(row.value) driver.find_element_by_id("add_to_basetable").click() sleep(2) if driver.find_element_by_id("goto_F003_edit"): driver.find_element_by_id("goto_F003_edit").click() sleep(2) for row in df.itertuples(): enter_forum(row) else: sleep(2) for row in df.itertuples(): enter_forum(row) driver.find_element_by_id("additional_function_on_submitting").click() これで出来上がり。上手く行けば自動でchromeが立ち上がりガンガン入力されていくはず。 最後に「入力を完了しますか？」というダイアログが出るのでクリックしたら終わりです。 もし品目を追加したい場合は、追加した品目だけをExcelで入力して、再実行すれば大丈夫です。 注意点として、品名は30文字以内でないといけないので、もし超えている場合はExcel sheet上で修正しないといけない。 最後に コードの全文を載せます。ご自由にお使いください。 全文.py import pandas as pd from selenium import webdriver from time import sleep import time import datetime from selenium.webdriver.support.select import Select df = pd.read_excel("研究遂行経費の支出報告書用入力シート.xlsx").dropna() df.columns = ["date", "item", "value", "is_tax","number", "expenditure", "is_receipt"] df['date'] = df['date'].dt.strftime("%Y-%m") def add_tax(x): if x.is_tax =="税抜": return x.value*1.1 else: return x.value def convert_date_to_flag(x): if x.date == "2021-04": return 1 elif x.date == "2021-05": return 2 elif x.date == "2021-06": return 3 elif x.date == "2021-07": return 4 elif x.date == "2021-08": return 5 elif x.date == "2021-09": return 6 elif x.date == "2021-10": return 7 elif x.date == "2021-11": return 8 elif x.date == "2021-12": return 9 elif x.date == "2022-01": return 10 elif x.date == "2022-02": return 11 elif x.date == "2022-03": return 12 def is_receipt(x): if x.is_receipt == "あり": return "etr_needOrNotReceipt_01" else: return "etr_needOrNotReceipt_02" def assign_expenditure(x): if x.expenditure == "学会関係経費": return 1 elif x.expenditure == "各種研究集会等への参加費": return 2 elif x.expenditure =="学術調査にかかる経費": return 3 elif x.expenditure == "自宅での研究に必要な経費": return 4 elif x.expenditure == "所属・関連機関への交通費": return 5 df['value'] = df.apply(lambda x: add_tax(x), axis = 1)*df["number"] df["date"] = df.apply(lambda x: convert_date_to_flag(x), axis = 1) df["is_receipt"] = df.apply(lambda x: is_receipt(x), axis = 1) df["expenditure"] = df.apply(lambda x: assign_expenditure(x), axis = 1) df_l = [list(row) for row in df.itertuples()] df = df[["date", "item", "expenditure", "is_receipt", "value"]] df["value"] = pd.Series(df["value"], dtype = 'int64') driver = webdriver.Chrome(executable_path='./chromedriver.exe') JSPS_URL = "https://area31.smp.ne.jp/area/servlet/area.MyPageBundle?MyPageID=629871a6_lalj9lenala0oare4" JSPS_ID = "202112786" JSPS_PASS = "RCyumJGu9Sdm5wP" J_SYSTEM_ID = "JU34996" J_SYSTEM_PASS = "ccmvMe3rDYL" driver.get(JSPS_URL) sleep(1) driver.find_element_by_name("SMPID").send_keys(JSPS_ID) driver.find_element_by_name("SMPPASSWORD").send_keys(JSPS_PASS) driver.find_element_by_class_name("btn").click() driver.find_element_by_link_text('⑤ 報告書等提出システム（Jシステム）').click() sleep(1) driver.switch_to.window(driver.window_handles[1]) driver.find_element_by_name("login_id").send_keys(J_SYSTEM_ID) driver.find_element_by_name("login_psw").send_keys(J_SYSTEM_PASS) driver.find_element_by_name("lg").click() sleep(1) driver.find_element_by_id("li_of_F003").click() sleep(1) def enter_forum(row): dropdown_expense = driver.find_element_by_id('etr_group') select_expense = Select(dropdown_expense) select_expense.select_by_index(row.expenditure) dropdown_date = driver.find_element_by_id('etr_dateOfComp') select_date = Select(dropdown_date) select_date.select_by_index(row.date) driver.find_element_by_id(str(row.is_receipt)).click() driver.find_element_by_name("etr_itemName").send_keys(row.item) driver.find_element_by_name("etr_billAmount").send_keys(row.value) driver.find_element_by_id("add_to_basetable").click() sleep(2) if driver.find_element_by_id("goto_F003_edit"): driver.find_element_by_id("goto_F003_edit").click() sleep(2) for row in df.itertuples(): enter_forum(row) else: sleep(2) for row in df.itertuples(): enter_forum(row) driver.find_element_by_id("additional_function_on_submitting").click()

前提 使い方、導入の仕方などはこちらの記事がが非常に有用でしたのでご参照ください。 この記事では使い方ではなく困った時の解決法などをいくつか書いていきます。また、随時更新します。 競合が起こった時 gitでもよくあることですがチームでバージョン管理していると競合が起こります。 alembic history --verbose こちらのコマンドでheadリビジョンを確認することができます。 バージョンはRivision IDで管理しているので競合を起こしているリビジョンのpyファイルを退避させましょう。 alembic upgrade head でテーブルが最新リビジョンと相違ない状態になったらsqlalchemyで定義してあるテーブルと退避させたpyファイルの相違している部分を修正し、 alembic revision --autogenerate しましょう。 alter文でつくられたリビジョンが生成されるので再度 alembic upgrade head で目的のリビジョンを反映させることができます。 headが2つ存在してしまった場合 例えばですが alembic upgrade head をすることがないまま alembic revision --autogenerate -m "create tables" を2回してしまった場合、headリビジョンが2つ存在してしまうという状況になります。 もちろんこのように単純なケースであれば古い方のリビジョンのpyファイルを消してしまえば解決するのですが、それぞれのリビジョンから派生して複数のalter_table.pyがつくられてしまった場合、後に alembic upgrade head すると FAILED: Multiple head revisions are present for given argument 'head'; please specify a specific target revision, '<branchname>@head' to narrow to a specific head, or 'heads' for all heads というエラーメッセージが出ます。headを1つにしてくださいという意味ですね。 じゃあリビジョンを指定してあげて... alembic revision -m "hogehoge" としたところで FAILED: Multiple heads are present; please specify the head revision on which the new revision should be based, or perform a merge. と出てしまう。 実はheadの管理って revision = 'hoge' down_revision = None ここでしてるんですよね。 なので revision = 'hoge' down_revision = 'fuga' と指定してあげましょう。 そうするとhogeがheadに、fugaが1つ前のリビジョンになります。 あとは alembic upgrade head すればいいだけです。 autogenerateでカラムの型変更が検出されない時の解決方法 env.pyの context.configure に、 compare_type=True を追加してあげましょう。 その後 alembic revision --autogenerate すれば検出してくれます。

id text 1 メロスは激怒した。 2 必ず、かの邪智暴虐の王を除かなければならぬと決意した。 3 メロスには政治がわからぬ。 4 メロスは、村の牧人である。 5 笛を吹き、羊と遊んで暮して来た。 tokenizer_obj = dictionary.Dictionary(dict="full").create() mode = tokenizer.Tokenizer.SplitMode.C doc = [] for row in range(len(df)): t = tokenizer_obj.tokenize(df["text"][row], mode) d = [m.normalized_form() for m in t if m.part_of_speech()[0] in ["名詞", "動詞"]] doc.append(d) docs = pd.array([" ".join(doc[i]) for i in range(len(doc))]) print(docs) #<StringArray> #['メロス 激怒 為る', #'邪知 暴虐 王 除く 成る 決意 為る', #'メロス 政治 分かる', #'メロス 村 牧人 有る', #'笛 吹く 羊 遊ぶ 暮らす 来る'] #Length: 5, dtype: string

概要 Slackアプリケーションでチャンネル名からチャンネルIDを取得したい。 しかし、API一覧では、channel関係はすべてdeprecatedになっている。 解決策 conversations.list APIを使うらしい。 target_channelのIDをtarget_channel_idに入れる例: from slack_bolt import App app = App(token=bot_token) channels = app.client.conversations_list()["channels"] try: target_channel_id = next(ch["id"] for ch in channels if ch["name"] == target_channel) except StopIteration: raise ValueError(f"{target_channel} is not found") 参考 How to find a channel id from channel name using a slack api?

GitHub Trending をキャッチアップする習慣をつけて、強強エンジニアになろう。 この記事では、Python のGithubのトレンドデイリーランキングを25位まで紹介します。 トレンドデイリーランキング 【1 位】 microsoft/routeros-scanner ? 574 star 【2 位】 lxgr-linux/pokete ? 673 star 【3 位】 public-apis/public-apis ? 186,353 star 【4 位】 vnpy/vnpy ? 17,929 star 【5 位】 google-research/kubric ? 1,043 star 【6 位】 home-assistant/core ? 50,883 star 【7 位】 diphosphane/CodeCraft2022-PressureGenerator ? 47 star 【8 位】 josephmisiti/awesome-machine-learning ? 53,702 star 【9 位】 freqtrade/freqtrade ? 16,767 star 【10 位】 kholia/OSX-KVM ? 13,540 star 【11 位】 projectdiscovery/nuclei-templates ? 3,613 star 【12 位】 iperov/DeepFaceLive ? 3,837 star 【13 位】 3b1b/videos ? 3,320 star 【14 位】 TeamUltroid/Ultroid ? 1,603 star 【15 位】 edilgin/DeepForSpeed ? 181 star 【16 位】 ethereum/web3.py ? 3,114 star 【17 位】 sherlock-project/sherlock ? 30,120 star 【18 位】 a16z/nft-analyst-starter-pack ? 194 star 【19 位】 datvuthanh/HybridNets ? 122 star 【20 位】 python-telegram-bot/python-telegram-bot ? 18,003 star 【21 位】 TheAlgorithms/Python ? 132,689 star 【22 位】 MHProDev/PyRoxy ? 13 star 【23 位】 Unknow101/FuckThatPacker ? 417 star 【24 位】 dgtlmoon/changedetection.io ? 3,876 star 【25 位】 ermaozi/get_subscribe ? 1,373 star 最後に 最後まで見ていただきありがとうございました。 LGTMをもらえるととても励みになりますので、ぜひお願いします

はじめに はじめまして．株式会社音圧爆上げくんにプロKagglerとして所属していますAshmeと申します． 業務の一環としてKaggleの様々なコンペティションに参加し，そこで得られた知見などを記事にして投稿しております．よろしくお願いいたします． 今回は現在のKaggleコンペティションに直接関係のあるものではありません．Twitterで流れてきたものですが教師なし特徴量選択の手法であるFRUFS(Feature Relevance based Unsupervised Feature Selection)という手法に興味があったのと，コンペティションで特徴量選択の1つとして使用することもできると考えたためこちらについて紹介します． こちらの手法はDeepwizAIにて紹介されていた手法です．元の記事はこちらにありますので，興味がある方はぜひこちらを読んでみてください． 本記事ではこのFRUFSと通常の教師あり学習における特徴量選択の違い，FRUFSはどのような考え方に基づいた手法なのか，どのように実装するのかについて紹介していきます． 教師あり学習における特徴量選択 まずこれまでもよく用いられてきた教師あり学習における特徴量選択について説明します．教師あり学習での特徴量選択は主に3種類に分類され，フィルター法，ラッパー法，埋め込み法が存在します．まずはこれらの３つの手法について簡単に説明していきます． フィルター法 フィルター法は特徴量選択の中で一番コストが低い手法だと思います．これは統計的な考え方に基づいて重要な特徴量を選択する手法です．例えば多重共線性という言葉を聞いたことがあるのではないでしょうか．相関係数が高い特徴量があると学習の結果が悪くなるということを聞いたことがあると思います．そこでEDAなどの早い段階で各特徴量間の相関係数を計算し，相関が大きい特徴量に関しては片方を不要として学習に用いないとする方法があります．ただこの多重共線性においては機械学習では予測精度を高くすることを目的とするため，そこまで深く考えないこともあるそうです．多重共線性につきましては調べれば簡単に出てきますので興味のある方は調べてみてください． ラッパー法 ラッパー法と埋め込み法は機械学習モデルを学習させることで特徴量選択を行いますが，それぞれ方法が少し異なります．ラッパー法は使用する特徴量を徐々に増やす，もしくは少なくすることで特徴量選択をします．徐々に増やす場合は使用する特徴量を１つ増やしたときに最もスコアが増加する特徴量を追加するということを繰り返します．逆に少なくしていく場合は使用する特徴量を1つ減らしたときに最もスコアが減少しない特徴量を削除するということを繰り返します．使用する特徴量の数を事前にしてしておき，その数に達したら処理を終了します． 文字だけでは少し分かりづらいため，簡単な例を図にしてみました．例は使用する特徴量を増やす場合です．ほとんど意味はないですが，もともと特徴量が3つのものを2つにしたいと考えます．今回の性能はAccuracyを仮定し1に近いほど良いとします． Step1では3つの特徴量を1つずつ入力してモデルを学習し性能を評価します．この場合，Feature Bが最も性能が高いためFeature Bは1つ目の特徴量として使用することが決まります． Step2で2つ目の特徴量を決めます．このときStep1で最も性能が高かったFeature Bは必ず用いることにし，Feature AとFeature Bを入力した場合とFeature BとFeature Cを入力としてモデルを学習し性能を評価します．結果としてFeature BとFeature Cを入力としたほうが性能が良いためこの2つを特徴量として使用することが決まります． 今回の場合は使用する特徴量は2つにすればよいためここで終了となります．例えば100個の特徴量を50個にしたい場合はこれを特徴量が50個決まるまで繰り返します． この手法の考え方としては最もスコアが増加する特徴量は残っている特徴量の中で予測するのに最も重要である，最もスコアが減少しない特徴量は特徴量の中で予測に最も重要でないという考えに基づいていると考えられます．この手法は機械学習のどの手法でも利用することができますが，何度もモデルを学習するためコストがかかってしまいます．ただし使用するモデルに対して最も良い特徴量の組み合わせを得ることができます． 埋め込み法 埋め込み法はラッパー法と同様に機械学習モデルを学習させることで特徴量選択を行いますが，こちらは学習は1度のみになります．ただし，使用できるモデルが限られており，線形回帰モデルかXGBoostやLightGBMなどのGBDTを用いることができます．線形回帰モデルではモデル学習した後に，各特徴量の重みを重要度と考え，重みの絶対値が大きいものほどその特徴量の重要度が高いとします．またGBDT系のFeature Importanceは大きいものほど重要な特徴量としています．本記事で紹介するFRUFSはこの埋め込み法の考え方に近いです． また実際の予測にこれら以外のモデルを使用する場合は，上記のモデルで重要度が高い特徴量を選択し，それらの特徴量を用いてモデルを学習することができます． 本題とは外れますが，特徴量の重要度の考え方は，上記のようなものだけでなく局所的な近似モデルを再現することで解釈を可能にするLIMEや，ゲーム理論におけるシャープレイ値の考え方を用いるSHAPという手法も存在します．また画像に関してもGrad-Cam，系列データに対してはTransformerに用いられているAttentionによってどの時間の特徴量が重要なのかを知ることができるなど様々な手法が提案されています．これらを知っておくとモデルがどの特徴量に注目して予測しているかの解釈を得ることができるため興味ある方は調べてみると良いと思います． FRUFS 教師なし学習での特徴量選択 教師あり学習での特徴量選択法について知ることができたので，本題であるFRUFSに入っていきたいと思います．FRUFSの紹介については元の記事の順序に基づいて紹介していきたいと思います． FRUFSの考え方に似ているものは埋め込み法です．まず埋め込み法を線形回帰モデルで学習したときのことを考えてみます．線形回帰モデルを定式化すると以下のようになっています．線形回帰モデルで特徴量重要度を理解する際にバイアス（切片）は不要なため式から除いています． $$\hat{y} = \alpha_1 x_1 + \cdots + \alpha_D x_D = \sum_{d=1}^D \alpha_d x_d$$ $\hat{y}$はモデルの予測値，$\alpha_d$はd番目の特徴量に対する重み，$x_d$はd番目の入力を表しています．また入力次元はDとしています．このとき重み$\alpha_d$はd番目の特徴量が予測をするためにどれだけ寄与しているか，重要であるかを表していると考えることができます．教師なし学習では正解ラベルがないため同様の手法を用いることができません．そこで予測するものを正解ラベルではなく，他の入力特徴量に置き換えます．上記と同様に式にすると以下のようになります． $$x_k = \sum_{d \neq k}^D \alpha_d x_d$$ つまりk番目の特徴量をそれ以外の特徴量で表現するということになります． これを用いてどのように特徴量の重要度を考えるのかと言うことについて説明します． 例として1番目の特徴量に焦点を当てて他の特徴量，つまり1番目以外の特徴量を表すときのことを考えます． まずk番目の特徴量を表すときにj番目の特徴量にかけられる重みを$w_{jk}$とします．このとき1番目の特徴量について，2番目の特徴量を表すときの重み$w_{12}$，3番目の特徴量を表すときの重み$w_{13}$…というように1番目の特徴量に対してかけられる重みはD-1個存在します．この重みの平均が大きいほどその特徴量の重要度が大きいと考えます．上の式を今定義した$w_{jk}$を用いて書き直すと以下のようになります． $$x_k = \sum_{j \neq k}^D w_{jk}x_j$$ なぜなら$w_{jk}$はk番目の特徴量を表すときにj番目の特徴量をどれだけ重要とされているかを表しているためです．j番目の特徴量が他の特徴量を表すときに重要とされているということは，他の特徴量の情報をj番目の特徴量が多く含んでいると考えられます．したがって入力データ全体を入力データの一部の特徴量（重要度が大きい特徴量）で表現できるということになります． この考え方は他の式を見ると少しわかりやすいかと思います．上記の重みは以下のように決定されます．このとき入力データの計画行列を$X(N \times D)$，重み行列を$W(D \times D)$としており，$W_{kj}$はj番目の特徴量を表すときのk番目の特徴量の重みを表しています．つまりj列目にはj番目の特徴量にかける重みが並んでいることになります．また$W$の対角成分はすべて0になります．これはj番目の特徴量を表すのにj番目の特徴量を用いないためです． $$X = XW^T$$ $$X - XW^T = 0$$ 上記の式を完全に満たすことはほぼできないため，結果としては$X - XW^T$のフロベニウスノルム（行列の各要素の2乗和の平方根）を最小化することになると元の記事にはありますが，簡単に言えば元の入力データともとの入力データに$W^T$をかけて作成したデータの各要素の差が0に近づくようになればよいということであり，$XW^T$によって元データ$X$をできる限り表現できれば良いということです． ここで扱った線形回帰モデルを用いての教師なし学習での特徴量選択はすでに2015年の研究論文で発表されているようです[1]． これを元にFRUFSは少し手を加えたものになっています． FRUFSと前節の手法の違い ではFRUFSでは前節の手法と何が異なるのでしょうか？これは非常に単純なもので，線形モデルの代わりに非線形モデルを用いるという部分だけです．ここでいう非線形モデルとしてLightGBMを扱っています．線形モデルでは特徴量の重要度として各特徴量に対する重みを持ちていました．その代わりにGBDT系で計算することができるFeature Importancesを用いることになります．なのでおそらくですがLightGBMでなくXGBoostやCatBoostなどの特徴量重要度を計算できるモデルであれば用いることができると考えられます． FRUFSの実装 FRUFSをMNISTデータセットに対して適用したものが元の記事で実装されていますが，私自身でも実装してみました．そこまで難しい実装ではないためみなさんも一度実装してみると良いと思います．私が実装の全体（FRUFSの実装以外は結果を比較したいため同じです）はこちらにありますが，元の記事では並列処理をうまく使用していたりするため，元の記事のものを参考にするほうが良いと思います． 私の実装ではPandasのDataFrameを入力し，メソッドを呼び出すことで重要度の計算，可視化ができるようにクラス化しています．また比較用に線形モデルとLightGBMのどちらも利用できるようにしています． # 他の特徴量を計算するための重み/重要度を計算する def calc_coef(self): X = self.df.values for i in tqdm(range(X.shape[1]), total=X.shape[1], desc="calculate coefficient/importances"): # 重み/重要度を格納するインデックスの指定 indices = np.concatenate((np.arange(i).reshape(-1, 1), np.arange(i+1, X.shape[1]).reshape(-1, 1))) train_X = np.hstack((X[:, :i], X[:, i+1:])) # i番目の特徴量を外す train_y = X[:, i] # i番目の特徴量 # i番目の特徴量を他の特徴量で表現するための学習 model = self.model_func() model.fit(train_X, train_y) # モデルの重み/重要度をWに格納する if self.method == "linear": coef = model.coef_ coef = np.absolute(coef) elif self.method == "lgb": coef = model.feature_importances_ self.W[i, indices] = coef.reshape(-1, 1) # 各特徴量が他の特徴量を表すときの重み/重要度の平均を計算 self.W_average = self.W.mean(axis=0) # 列方向に平均を取る = ある特徴量にかけられる重みの平均 self.average_coef_df = pd.DataFrame({"columns": self.df.columns.values, "importances": self.W_average}) 実際にある特徴量を他の特徴量を用いて表現するときの特徴量重要度を計算し，データを保存する部分になります．実装の中でインデックスの部分が少し分かりづらいかと思うので説明します． # 重み/重要度を格納するインデックスの指定 indices = np.concatenate((np.arange(i).reshape(-1, 1), np.arange(i+1, X.shape[1]).reshape(-1, 1))) train_X = np.hstack((X[:, :i], X[:, i+1:])) # i番目の特徴量を外す train_y = X[:, i] # i番目の特徴量 この辺の部分は配列を扱う機会があればわかりやすいかと思いますが，あまりない方向けに説明します．ここには載っていませんが，重み$W$をインスタンスの初期化のときに作成しています．このとき$W$は説明の通り$D \times D$のゼロ行列として初期化しています．j列目はj番目の特徴量が他の特徴量を表現するときにかけられる重みが並んでいますが，k行目にはk番目の特徴量を表すときにk番目の特徴量以外にかける重みが並んでいます． 特徴量重要度の計算はk番目の特徴量を計算するときのそれ以外の特徴量重要度を得ることになるため，結果をk行名に格納することになります．またこのときk行k列目は0である必要があるためk番目はindicesに含まないようにしています． そして入力データとしてk番目以外の特徴量，出力の教師ラベルとしてk番目の特徴量をそれぞれtrain_X, train_yとしています．これを用いてモデルを学習させることになります． 結果 元の記事と同様の可視化やスコアの算出をしたところ以下のようになりました．データのとり方が少し異なっている（私の場合はscikit-learnを用いている），データの一部を用いており再現性確保をしていないため，多少結果がぶれますが，元の記事の実装とほとんど遜色ない結果を得ることができました．論文実装など様々なものがありますが，そういうものに慣れていない方でもこちらの実装はわかりやすいと思うので練習と思って一度実装してみると良いと思います． また比較用に線形モデルでの結果も載せておきます．こちらの結果だと特徴量すべてを用いたものが最もスコアが高いということでうまく特徴量選択できていないことが解ると思います．ただ流石に悪すぎる気もするので実装を一度見直してみようと思います． おわりに 今回は現在開催されているコンペティションに直接関係はしていませんが，教師なし学習での特徴量選択手法FRUFSについて紹介しました．結果を見ても教師なし学習のコンペティションのときに特徴量選択として使用できる可能性もあるのではないでしょうか．また教師あり学習のときにこの手法で特徴量選択をしてもすべての特徴量より良い性能を示すこともあるそうなので使用してみようと思っています． 最後に人材募集となりますが，株式会社音圧爆上げくんではプロKagglerを募集しています． 少しでも興味の有る方はぜひ以下のリンクをご覧の上ご応募ください． Wantedlyリンク 参考 [1] Zhu, Pengfei, Wangmeng Zuo, Lei Zhang, Qinghua Hu, and Simon CK Shiu. "Unsupervised feature selection by regularized self-representation." Pattern Recognition 48, no. 2 (2015): 438-446.

与えられた$n$について、素数かどうか判定したい。 アルゴリズムをいろいろ考えてみよう 2022/03/23 @StrawBerryMoon さん、ご指摘ありがとうございました。 1. ひたすら割る $2$以上の数でひたすら割っていき、割れなければ素数、割れれば合成数。 def trial(n: int): assert n > 1 i = 2 while i < n: if n%i == 0: return False i+=1 return True 2. 偶数は飛ばす $n\gt2$であれば、素数は全て奇数になる。 1よりも早く判定ができる。 def trial_odd(n: int): assert n > 1 if n == 2: return True elif n % 2 == 0: return False i = 3 while i < n: if n%i == 0: return False i += 2 return True 3. 探索範囲をsqrt(n)に絞る 実は$\sqrt{n}$まで調べれば、素数かどうかは判別できる。 $n=2^{32}-1$だとしても2と組み合わせてだいたい$32767$回で済める def trial_odd_sqrt(n: int): assert n > 1 if n == 2: return True elif n % 2 == 0: return False i = 3 while i**2 < n: if n%i == 0: return False i+=2 return True 4. Wheel factorizationを使う これを知っているひとは数学好きに違いない。 $\{2, 3, 5\}$が既知の素数であるとし、$i=7$から探索を開始する。このとき、$\{4, 2, 4, 4, 6, 2, 6\}$の周期で数値を加算しながら素数かどうかを確認できる。 偶数飛ばしよりもさらに(プログラム内では誤差かもしれないが)飛ばしていける。 def trial_wheel_factorize(n: int): assert n > 1 if n in [2, 3, 5]: return True elif n%2 == 0 or n%3 == 0 or n%5 == 0: return False i = 7 c = 0 wheel = [4, 2, 4, 2, 4, 6, 2, 6] while i < n: if n%i == 0: return False i += wheel[c%8] c += 1 return True ただし、既知の数を増やしても使えるというわけではない。理由はここを見るか、下の周期一覧を開いてもらえれば分かると思う。 周期一覧 $p = \{2, 3, 5\} \rightarrow \mathrm{wheel}=\{4, 2, 4, 2, 4, 6, 2, 6\}$ $p = \{2, 3, 5, 7\} \rightarrow \mathrm{wheel}=\{2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10\}$ $p = \{2, 3, 5, 7, 11\} \rightarrow \mathrm{wheel}=\{4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 14, 4, 6, 2, 10, 2, 6, 6, 4, 2, 4, 6, 2, 10, 2, 4, 2, 12, 10, 2, 4, 2, 4, 6, 2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 6, 8, 6, 10, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 6, 10, 2, 10, 2, 4, 2, 4, 6, 8, 4, 2, 4, 12, 2, 6, 4, 2, 6, 4, 6, 12, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 10, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 6, 6, 2, 6, 6, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 2, 10, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 8, 6, 4, 2, 4, 6, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 6, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 6, 4, 6, 2, 6, 4, 2, 4, 6, 6, 8, 4, 2, 6, 10, 8, 4, 2, 4, 2, 4, 8, 10, 6, 2, 4, 8, 6, 6, 4, 2, 4, 6, 2, 6, 4, 6, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 6, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 6, 4, 2, 4, 6, 8, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 10, 2, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 6, 6, 2, 6, 6, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 10, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 12, 6, 4, 6, 2, 4, 6, 2, 12, 4, 2, 4, 8, 6, 4, 2, 4, 2, 10, 2, 10, 6, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 10, 6, 8, 6, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 2, 4, 2, 10, 12, 2, 4, 2, 10, 2, 6, 4, 2, 4, 6, 6, 2, 10, 2, 6, 4, 14, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 12, 2, 12\}$ $p = \{2, 3, 5, 7, 11, 13\} \rightarrow　\mathrm{wheel}=$ {2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 14, 4, 6, 2, 10, 2, 6, 6, 4, 6, 6, 2, 10, 2, 4, 2, 12, 12, 4, 2, 4, 6, 2, 10, 6, 6, 6, 2, 6, 4, 2, 6, 4, 14, 4, 2, 4, 6, 8, 6, 10, 2, 4, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 8, 10, 2, 10, 2, 4, 2, 4, 6, 8, 4, 2, 4, 12, 8, 4, 2, 6, 4, 6, 12, 2, 4, 2, 12, 6, 4, 6, 6, 6, 2, 6, 10, 2, 4, 6, 2, 6, 6, 4, 2, 10, 2, 10, 2, 4, 6, 6, 2, 6, 6, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 8, 4, 6, 8, 10, 2, 10, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 14, 4, 2, 4, 6, 6, 2, 6, 4, 8, 10, 8, 4, 2, 4, 6, 8, 6, 4, 6, 6, 6, 2, 6, 6, 4, 2, 4, 6, 2, 10, 2, 4, 2, 10, 2, 10, 2, 6, 4, 8, 6, 4, 2, 4, 6, 6, 8, 4, 2, 6, 10, 8, 4, 2, 6, 4, 8, 10, 6, 2, 4, 8, 6, 6, 4, 2, 4, 6, 2, 6, 4, 6, 2, 10, 12, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 12, 2, 6, 6, 10, 6, 8, 4, 2, 4, 2, 4, 8, 6, 12, 4, 6, 2, 12, 4, 2, 4, 6, 8, 4, 2, 4, 2, 12, 10, 2, 4, 2, 4, 6, 2, 10, 2, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 6, 6, 8, 6, 4, 2, 10, 2, 10, 2, 4, 2, 10, 2, 6, 4, 2, 10, 6, 2, 6, 4, 2, 6, 4, 6, 8, 6, 4, 2, 12, 10, 6, 2, 4, 6, 2, 12, 4, 2, 4, 8, 6, 4, 2, 4, 2, 10, 2, 10, 6, 2, 4, 6, 2, 6, 4, 2, 10, 6, 2, 6, 4, 12, 6, 8, 6, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 6, 4, 6, 2, 6, 4, 2, 4, 2, 10, 12, 2, 4, 12, 2, 6, 4, 2, 4, 6, 6, 2, 12, 6, 4, 18, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 8, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 12, 2, 12, 6, 4, 6, 2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 14, 4, 6, 2, 10, 2, 6, 6, 4, 2, 10, 2, 10, 2, 4, 14, 10, 2, 4, 2, 4, 6, 2, 6, 10, 6, 6, 2, 10, 2, 6, 4, 6, 8, 4, 2, 4, 6, 8, 6, 10, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 6, 10, 2, 10, 6, 2, 4, 6, 8, 4, 2, 4, 12, 2, 6, 4, 2, 6, 4, 6, 12, 2, 4, 2, 4, 8, 6, 4, 6, 2, 10, 2, 6, 10, 6, 6, 2, 6, 4, 2, 4, 2, 10, 2, 12, 4, 6, 6, 2, 12, 4, 6, 6, 2, 6, 4, 2, 6, 4, 14, 4, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 6, 10, 2, 6, 4, 6, 2, 10, 2, 10, 2, 4, 2, 4, 8, 6, 4, 2, 4, 6, 6, 8, 4, 8, 4, 6, 8, 4, 2, 4, 2, 12, 6, 4, 6, 6, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 6, 4, 2, 10, 2, 10, 2, 6, 4, 6, 2, 6, 4, 2, 4, 6, 14, 4, 2, 6, 10, 8, 4, 2, 4, 2, 4, 8, 10, 8, 4, 8, 6, 10, 2, 4, 6, 2, 6, 4, 6, 2, 10, 2, 10, 2, 4, 2, 4, 6, 8, 4, 2, 4, 6, 6, 2, 6, 6, 6, 10, 8, 4, 2, 4, 6, 8, 6, 4, 8, 4, 6, 2, 6, 6, 4, 2, 4, 6, 12, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 8, 10, 2, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 8, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 6, 6, 2, 6, 6, 4, 2, 10, 12, 2, 4, 2, 4, 6, 2, 6, 4, 2, 16, 2, 6, 4, 2, 10, 6, 8, 4, 2, 4, 2, 12, 6, 10, 2, 4, 6, 2, 12, 4, 2, 4, 8, 6, 4, 2, 4, 2, 12, 10, 6, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 10, 6, 8, 10, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 6, 4, 6, 8, 4, 2, 4, 2, 10, 12, 2, 4, 2, 10, 2, 6, 4, 2, 4, 6, 6, 2, 10, 2, 6, 4, 14, 6, 4, 2, 4, 8, 10, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 8, 6, 4, 2, 4, 12, 2, 12, 4, 2, 4, 6, 2, 6, 4, 2, 10, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 4, 2, 4, 14, 4, 6, 2, 10, 8, 6, 4, 2, 4, 6, 2, 10, 2, 4, 2, 12, 10, 2, 4, 6, 6, 2, 6, 4, 6, 6, 6, 2, 6, 6, 6, 4, 6, 12, 2, 4, 6, 8, 6, 10, 2, 4, 8, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 6, 10, 2, 10, 2, 6, 4, 6, 8, 4, 6, 12, 2, 6, 4, 2, 6, 4, 6, 12, 2, 4, 2, 4, 14, 4, 6, 2, 4, 6, 2, 6, 10, 2, 10, 2, 6, 4, 2, 4, 12, 2, 10, 2, 4, 6, 6, 2, 6, 6, 4, 6, 6, 2, 10, 2, 6, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 2, 10, 2, 6, 4, 2, 6, 10, 2, 10, 6, 2, 4, 8, 6, 4, 2, 4, 6, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 12, 2, 6, 6, 4, 6, 6, 2, 6, 4, 2, 4, 2, 10, 2, 12, 6, 4, 6, 2, 10, 2, 4, 6, 6, 8, 4, 2, 6, 18, 4, 2, 4, 2, 4, 8, 10, 6, 2, 4, 8, 6, 6, 6, 4, 6, 2, 6, 4, 6, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 8, 6, 6, 4, 6, 8, 4, 2, 4, 2, 12, 6, 4, 12, 6, 2, 6, 6, 4, 2, 4, 6, 8, 6, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 10, 2, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 10, 6, 6, 2, 6, 6, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 8, 4, 2, 10, 6, 2, 6, 4, 2, 6, 10, 8, 4, 2, 4, 14, 6, 4, 6, 2, 4, 6, 2, 12, 4, 2, 4, 8, 10, 2, 4, 2, 10, 2, 10, 6, 2, 4, 8, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 10, 6, 8, 6, 6, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 2, 4, 2, 10, 12, 2, 4, 2, 10, 2, 6, 4, 2, 4, 12, 2, 10, 2, 10, 14, 4, 2, 4, 2, 4, 8, 6, 10, 2, 4, 6, 2, 12, 4, 2, 4, 6, 2, 6, 4, 2, 4, 14, 12, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 6, 2, 4, 14, 4, 6, 2, 10, 2, 6, 6, 4, 2, 4, 6, 12, 2, 4, 2, 12, 10, 2, 4, 2, 10, 2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 6, 4, 6, 8, 16, 2, 4, 6, 2, 6, 6, 4, 2, 4, 8, 6, 4, 2, 6, 10, 2, 10, 2, 4, 2, 4, 6, 8, 4, 2, 16, 2, 6, 4, 8, 4, 6, 12, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 8, 10, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 6, 6, 2, 6, 6, 4, 6, 6, 2, 6, 6, 6, 4, 6, 12, 2, 6, 4, 8, 6, 4, 6, 2, 4, 14, 6, 4, 2, 10, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 6, 4, 8, 6, 4, 6, 6, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 4, 2, 4, 14, 4, 6, 6, 6, 2, 6, 6, 4, 2, 10, 2, 6, 4, 2, 4, 12, 2, 10, 2, 6, 4, 6, 2, 6, 6, 4, 6, 6, 12, 2, 6, 10, 8, 4, 2, 4, 2, 4, 8, 10, 6, 2, 4, 8, 6, 6, 4, 2, 4, 6, 2, 6, 4, 8, 10, 2, 10, 6, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 6, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 8, 10, 2, 6, 6, 4, 6, 6, 8, 4, 2, 4, 2, 10, 2, 12, 4, 2, 4, 6, 2, 10, 2, 4, 6, 8, 6, 4, 2, 6, 4, 14, 4, 6, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 10, 6, 2, 6, 10, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 10, 6, 8, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 12, 6, 4, 6, 6, 6, 2, 12, 4, 2, 4, 8, 6, 6, 4, 2, 10, 2, 10, 6, 2, 4, 6, 2, 6, 4, 2, 4, 6, 8, 6, 4, 2, 10, 6, 8, 6, 4, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 12, 4, 6, 2, 6, 4, 2, 4, 2, 10, 12, 2, 4, 2, 10, 8, 4, 2, 4, 6, 6, 2, 10, 2, 6, 18, 4, 2, 4, 6, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 10, 2, 4, 12, 2, 12, 4, 2, 4, 8, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 6, 4, 14, 4, 6, 2, 10, 2, 6, 6, 4, 2, 4, 6, 2, 10, 2, 4, 2, 22, 2, 4, 2, 4, 6, 2, 6, 4, 6, 12, 2, 6, 4, 2, 10, 6, 8, 4, 2, 4, 6, 8, 6, 10, 2, 4, 6, 2, 12, 4, 2, 4, 6, 2, 6, 4, 2, 6, 12, 10, 2, 4, 2, 4, 6, 8, 4, 2, 4, 12, 2, 6, 4, 2, 6, 4, 6, 12, 6, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 10, 2, 4, 6, 8, 4, 2, 4, 2, 10, 2, 10, 2, 4, 12, 2, 6, 6, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 6, 6, 4, 8, 10, 6, 2, 4, 6, 8, 6, 4, 2, 12, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 8, 6, 4, 2, 10, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 6, 6, 8, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 10, 6, 2, 6, 4, 2, 4, 6, 6, 8, 6, 6, 10, 12, 2, 4, 2, 4, 8, 10, 6, 2, 4, 8, 6, 6, 4, 2, 4, 6, 2, 6, 4, 6, 2, 10, 2, 10, 2, 6, 4, 6, 2, 6, 4, 6, 6, 6, 2, 6, 6, 6, 4, 6, 8, 4, 2, 4, 2, 4, 14, 4, 8, 4, 6, 2, 6, 6, 4, 2, 10, 8, 4, 2, 4, 12, 2, 10, 2, 4, 2, 4, 6, 2, 12, 4, 6, 8, 10, 2, 6, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 6, 6, 2, 6, 6, 6, 10, 2, 10, 6, 2, 4, 6, 2, 6, 4, 2, 10, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 12, 6, 4, 6, 2, 10, 2, 12, 4, 6, 8, 6, 4, 2, 4, 2, 10, 2, 16, 2, 4, 6, 2, 10, 2, 4, 6, 6, 2, 6, 4, 2, 10, 14, 6, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 6, 2, 10, 12, 2, 4, 2, 10, 2, 6, 4, 2, 4, 6, 6, 12, 2, 6, 4, 14, 4, 2, 4, 2, 12, 6, 4, 6, 6, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 6, 4, 12, 2, 12, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 14, 4, 8, 10, 2, 6, 10, 2, 4, 6, 2, 10, 2, 4, 2, 12, 10, 2, 4, 2, 4, 6, 8, 4, 6, 6, 6, 2, 6, 4, 2, 6, 10, 8, 4, 2, 4, 6, 8, 6, 10, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 10, 2, 6, 10, 2, 10, 2, 4, 2, 4, 14, 4, 2, 4, 12, 2, 6, 4, 2, 6, 4, 6, 12, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 10, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 12, 2, 4, 6, 6, 2, 6, 6, 4, 12, 2, 6, 4, 2, 10, 6, 8, 4, 2, 6, 4, 8, 6, 10, 2, 4, 6, 14, 4, 2, 10, 2, 6, 4, 2, 4, 2, 12, 10, 2, 4, 2, 4, 8, 6, 4, 2, 4, 6, 6, 2, 6, 4, 8, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 6, 6, 6, 2, 6, 6, 4, 2, 4, 6, 8, 4, 2, 4, 2, 10, 2, 10, 2, 6, 10, 2, 6, 4, 2, 4, 6, 6, 8, 4, 2, 6, 10, 8, 6, 4, 2, 4, 8, 10, 6, 2, 4, 8, 6, 6, 4, 2, 4, 8, 6, 4, 6, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 10, 6, 2, 6, 12, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 8, 4, 6, 8, 6, 4, 2, 4, 6, 8, 4, 2, 4, 2, 10, 2, 10, 2, 4, 6, 6, 2, 10, 2, 4, 6, 8, 6, 6, 6, 4, 6, 12, 6, 2, 4, 8, 6, 4, 6, 2, 4, 8, 6, 6, 4, 6, 6, 2, 6, 6, 4, 2, 10, 2, 10, 2, 6, 4, 6, 2, 6, 4, 12, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 18, 4, 6, 2, 4, 6, 2, 12, 4, 2, 12, 6, 4, 2, 4, 12, 2, 10, 6, 2, 4, 6, 2, 6, 6, 4, 6, 6, 2, 10, 2, 10, 6, 8, 6, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 2, 6, 10, 12, 6, 2, 10, 2, 6, 4, 2, 4, 6, 6, 2, 10, 2, 6, 4, 14, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 10, 2, 6, 6, 4, 6, 6, 2, 6, 4, 2, 4, 12, 2, 12, 4, 2, 4, 6, 2, 10, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 14, 4, 2, 4, 2, 4, 14, 4, 6, 2, 10, 2, 6, 6, 6, 4, 6, 2, 10, 6, 2, 12, 10, 2, 4, 2, 4, 6, 2, 6, 4, 6, 6, 6, 8, 4, 2, 6, 4, 6, 8, 4, 2, 4, 14, 6, 10, 6, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 6, 6, 10, 2, 10, 2, 4, 2, 4, 6, 8, 4, 2, 4, 14, 6, 4, 2, 6, 4, 6, 12, 2, 4, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 10, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 6, 6, 8, 6, 4, 6, 6, 2, 6, 4, 2, 6, 10, 8, 4, 2, 10, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 2, 10, 2, 10, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 8, 6, 4, 2, 4, 6, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 6, 6, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 12, 2, 6, 4, 6, 2, 6, 4, 2, 4, 12, 8, 4, 2, 16, 8, 4, 2, 4, 2, 4, 8, 16, 2, 4, 8, 12, 4, 2, 4, 6, 2, 6, 4, 6, 2, 12, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 6, 6, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 6, 4, 2, 4, 6, 8, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 10, 2, 10, 2, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 10, 2, 4, 8, 10, 6, 2, 4, 6, 2, 6, 6, 4, 6, 8, 6, 6, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 10, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 4, 2, 12, 6, 4, 6, 2, 4, 6, 14, 4, 2, 4, 8, 6, 4, 2, 4, 2, 10, 2, 10, 6, 6, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 6, 10, 6, 14, 4, 2, 4, 8, 6, 4, 6, 2, 4, 8, 6, 6, 6, 4, 6, 2, 6, 4, 2, 4, 2, 10, 12, 2, 6, 10, 2, 6, 4, 6, 6, 6, 2, 10, 2, 6, 4, 14, 4, 2, 4, 2, 4, 14, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 10, 2, 6, 4, 2, 4, 12, 2, 12, 4, 2, 4, 6, 2, 6, 6, 4, 6, 6, 2, 10, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 14, 4, 6, 2, 10, 2, 6, 6, 4, 2, 4, 6, 2, 10, 2, 6, 12, 10, 6, 2, 4, 6, 2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 6, 8, 6, 10, 2, 10, 2, 6, 6, 4, 6, 6, 2, 6, 4, 2, 6, 10, 2, 12, 4, 2, 4, 6, 12, 2, 4, 12, 2, 6, 4, 2, 6, 4, 18, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 12, 4, 6, 2, 6, 4, 6, 2, 10, 2, 10, 2, 4, 6, 6, 2, 6, 6, 4, 6, 6, 8, 4, 2, 6, 4, 6, 8, 4, 2, 6, 12, 6, 4, 6, 6, 6, 8, 6, 4, 2, 10, 2, 6, 6, 4, 2, 10, 2, 10, 2, 4, 2, 4, 8, 6, 4, 2, 4, 6, 8, 6, 4, 8, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 12, 6, 2, 6, 10, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 6, 4, 6, 8, 4, 2, 4, 6, 6, 8, 4, 2, 6, 10, 8, 4, 2, 4, 6, 8, 10, 6, 2, 4, 8, 6, 6, 4, 2, 4, 6, 2, 10, 6, 2, 10, 2, 10, 2, 4, 2, 4, 8, 6, 4, 2, 4, 6, 6, 2, 6, 6, 6, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 8, 4, 6, 2, 6, 6, 4, 2, 4, 6, 8, 4, 2, 4, 2, 10, 12, 2, 4, 2, 4, 6, 2, 10, 2, 4, 14, 6, 4, 2, 10, 6, 8, 4, 6, 2, 4, 8, 6, 10, 2, 4, 6, 2, 12, 4, 6, 6, 2, 6, 6, 4, 2, 12, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 10, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 6, 2, 12, 6, 4, 6, 2, 4, 6, 2, 12, 4, 2, 4, 14, 4, 2, 4, 2, 10, 2, 10, 6, 2, 10, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 10, 6, 8, 6, 4, 2, 4, 8, 10, 6, 2, 4, 6, 2, 6, 6, 6, 4, 8, 6, 4, 2, 4, 2, 10, 12, 2, 4, 2, 10, 2, 6, 4, 2, 10, 6, 2, 10, 8, 4, 14, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 12, 2, 12, 4, 6, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 6, 6, 4, 6, 12, 2, 4, 2, 4, 14, 4, 6, 2, 12, 6, 6, 4, 2, 4, 6, 2, 10, 2, 4, 2, 12, 10, 2, 6, 4, 6, 2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 6, 14, 10, 2, 4, 6, 2, 6, 6, 4, 2, 10, 2, 6, 4, 2, 16, 2, 10, 2, 4, 2, 4, 6, 8, 6, 4, 12, 2, 10, 2, 6, 4, 6, 12, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 10, 2, 4, 6, 2, 6, 4, 2, 6, 10, 2, 10, 6, 6, 6, 2, 6, 6, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 6, 2, 10, 8, 6, 4, 12, 2, 6, 4, 2, 4, 2, 10, 2, 12, 4, 2, 4, 8, 10, 2, 4, 6, 6, 2, 6, 4, 8, 4, 14, 4, 2, 4, 2, 4, 8, 6, 4, 6, 6, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 6, 2, 10, 2, 10, 2, 6, 4, 6, 2, 6, 4, 2, 4, 6, 6, 8, 4, 2, 6, 10, 8, 4, 2, 4, 2, 12, 10, 6, 6, 8, 6, 6, 4, 2, 4, 6, 2, 6, 10, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 8, 6, 6, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 10, 2, 4, 6, 8, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 12, 2, 4, 6, 8, 6, 4, 2, 6, 10, 8, 4, 6, 6, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 6, 6, 2, 12, 4, 2, 10, 2, 10, 2, 4, 2, 4, 8, 6, 4, 2, 10, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 6, 12, 6, 4, 6, 2, 4, 6, 2, 12, 4, 2, 4, 8, 6, 4, 2, 4, 2, 10, 12, 6, 2, 4, 6, 2, 6, 4, 2, 4, 12, 2, 6, 4, 2, 10, 6, 8, 6, 4, 2, 4, 8, 6, 10, 2, 4, 6, 2, 12, 6, 4, 6, 2, 6, 4, 2, 4, 2, 22, 2, 4, 2, 10, 2, 6, 4, 2, 4, 6, 6, 2, 10, 2, 6, 4, 14, 4, 6, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 8, 4, 2, 4, 12, 2, 12, 4, 2, 10, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 6, 4, 2, 4, 18, 6, 2, 10, 2, 6, 6, 4, 2, 4, 8, 10, 2, 4, 2, 12, 10, 2, 4, 2, 4, 6, 2, 6, 4, 12, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 4, 6, 8, 6, 10, 2, 4, 6, 8, 6, 4, 2, 4, 6, 2, 6, 4, 2, 6, 10, 2, 10, 2, 4, 6, 6, 8, 4, 2, 4, 12, 2, 6, 6, 6, 4, 6, 12, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 8, 6, 10, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 10, 6, 2, 6, 10, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 6, 4, 14, 4, 6, 2, 4, 6, 8, 6, 4, 2, 10, 2, 6, 4, 2, 4, 12, 2, 10, 2, 4, 2, 4, 8, 6, 6, 4, 6, 6, 2, 10, 8, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 6, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 6, 10, 2, 10, 8, 4, 6, 2, 6, 4, 2, 4, 6, 6, 8, 4, 2, 6, 10, 8, 4, 2, 4, 2, 4, 8, 10, 6, 2, 12, 6, 6, 4, 6, 6, 2, 6, 4, 6, 2, 10, 2, 12, 4, 2, 4, 6, 2, 10, 2, 4, 6, 6, 2, 6, 6, 6, 4, 14, 4, 2, 4, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 6, 6, 4, 6, 8, 4, 6, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 10, 2, 4, 6, 14, 4, 2, 6, 4, 6, 8, 4, 6, 2, 12, 6, 4, 6, 6, 6, 2, 6, 6, 4, 6, 6, 2, 6, 6, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 10, 8, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 12, 6, 4, 8, 4, 6, 2, 16, 2, 4, 8, 6, 4, 2, 4, 2, 10, 2, 10, 6, 2, 4, 6, 8, 4, 2, 4, 6, 6, 2, 6, 4, 2, 16, 8, 6, 4, 6, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 6, 4, 6, 2, 10, 2, 4, 2, 10, 12, 2, 4, 2, 12, 6, 4, 2, 4, 6, 6, 2, 10, 2, 6, 4, 14, 4, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 12, 14, 4, 2, 4, 6, 2, 6, 4, 2, 4, 12, 2, 6, 4, 2, 10, 6, 8, 4, 2, 4, 2, 4, 14, 10, 2, 10, 2, 12, 4, 2, 4, 6, 2, 10, 2, 4, 2, 12, 10, 2, 4, 2, 4, 6, 2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 6, 6, 8, 6, 10, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 8, 4, 2, 6, 10, 2, 10, 2, 4, 2, 10, 8, 4, 2, 4, 12, 2, 6, 4, 2, 6, 4, 6, 14, 4, 2, 4, 8, 10, 6, 2, 4, 6, 2, 6, 10, 2, 4, 8, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 6, 6, 2, 6, 6, 10, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 2, 10, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 6, 8, 6, 4, 2, 4, 6, 6, 2, 6, 12, 4, 6, 12, 2, 4, 2, 4, 8, 6, 4, 6, 6, 8, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 6, 4, 6, 2, 6, 4, 6, 6, 6, 8, 4, 2, 6, 10, 8, 4, 2, 4, 2, 4, 18, 6, 2, 4, 8, 6, 6, 4, 2, 10, 2, 6, 4, 6, 12, 2, 10, 2, 4, 2, 4, 6, 2, 6, 6, 4, 6, 6, 2, 12, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 6, 4, 2, 4, 6, 8, 4, 2, 6, 10, 2, 10, 6, 2, 4, 6, 2, 10, 2, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 6, 2, 10, 2, 6, 6, 4, 6, 6, 2, 6, 6, 4, 2, 10, 2, 12, 4, 2, 4, 6, 2, 10, 2, 10, 6, 2, 6, 4, 2, 6, 4, 14, 4, 2, 4, 2, 12, 6, 4, 6, 2, 4, 6, 2, 12, 6, 4, 8, 6, 4, 6, 2, 10, 2, 10, 6, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 8, 4, 2, 10, 6, 8, 6, 4, 2, 12, 6, 4, 6, 6, 6, 2, 6, 6, 6, 4, 6, 2, 6, 6, 4, 2, 10, 12, 2, 4, 2, 10, 2, 6, 4, 2, 4, 6, 8, 10, 2, 6, 4, 14, 4, 2, 4, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 10, 2, 4, 6, 2, 6, 4, 2, 4, 12, 2, 12, 4, 2, 4, 6, 8, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 10, 8, 4, 2, 4, 6, 14, 4, 6, 2, 10, 2, 6, 6, 4, 2, 4, 6, 2, 10, 2, 4, 2, 12, 10, 2, 4, 2, 4, 8, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 10, 8, 6, 10, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 6, 10, 12, 2, 4, 2, 4, 6, 8, 4, 2, 4, 12, 2, 6, 4, 2, 10, 6, 12, 2, 4, 2, 4, 8, 6, 10, 2, 4, 6, 2, 16, 2, 4, 6, 2, 6, 4, 2, 4, 2, 12, 10, 2, 4, 6, 6, 2, 6, 6, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 8, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 2, 10, 8, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 12, 6, 4, 2, 4, 6, 6, 2, 6, 4, 8, 4, 6, 8, 6, 4, 2, 4, 8, 10, 6, 6, 6, 2, 6, 6, 4, 2, 4, 8, 6, 4, 2, 4, 2, 10, 2, 10, 2, 6, 4, 6, 2, 6, 4, 2, 10, 6, 8, 4, 8, 10, 8, 4, 2, 4, 2, 4, 8, 10, 6, 2, 4, 14, 6, 4, 2, 4, 6, 2, 6, 4, 6, 2, 10, 2, 10, 2, 4, 6, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 6, 6, 4, 6, 12, 2, 4, 2, 4, 8, 6, 4, 8, 4, 8, 6, 6, 4, 2, 4, 6, 8, 4, 2, 4, 2, 10, 2, 10, 2, 6, 4, 6, 2, 10, 6, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 6, 2, 4, 14, 4, 6, 2, 4, 6, 2, 6, 6, 4, 12, 2, 6, 6, 4, 12, 2, 10, 2, 4, 2, 4, 6, 2, 6, 6, 10, 6, 2, 10, 2, 6, 4, 6, 8, 4, 2, 4, 2, 12, 6, 4, 6, 2, 4, 6, 2, 12, 4, 2, 4, 8, 6, 4, 2, 6, 10, 2, 10, 6, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 10, 6, 8, 6, 4, 2, 4, 8, 6, 4, 6, 2, 10, 2, 6, 6, 10, 6, 2, 6, 4, 2, 4, 2, 10, 14, 4, 2, 10, 2, 10, 2, 4, 6, 6, 2, 10, 2, 6, 4, 14, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 6, 12, 2, 12, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 8, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 18, 4, 6, 12, 2, 6, 6, 4, 2, 4, 6, 2, 12, 4, 2, 12, 10, 2, 4, 2, 4, 6, 2, 6, 4, 6, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 6, 8, 6, 12, 4, 6, 2, 6, 10, 2, 4, 6, 2, 6, 4, 2, 6, 10, 2, 10, 2, 4, 2, 4, 6, 8, 4, 2, 4, 12, 2, 6, 4, 2, 6, 10, 12, 2, 4, 6, 8, 6, 4, 6, 2, 4, 6, 2, 6, 10, 2, 4, 6, 2, 10, 2, 4, 2, 10, 2, 10, 2, 4, 6, 8, 6, 6, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 2, 10, 2, 6, 4, 2, 4, 2, 10, 12, 2, 4, 2, 4, 8, 6, 4, 2, 4, 12, 2, 6, 4, 12, 6, 8, 4, 2, 4, 2, 4, 8, 6, 10, 6, 6, 2, 12, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 12, 10, 2, 6, 4, 6, 2, 6, 4, 2, 4, 6, 6, 8, 4, 2, 6, 10, 8, 4, 6, 2, 4, 8, 10, 6, 2, 4, 8, 6, 6, 4, 2, 4, 6, 8, 4, 6, 2, 10, 2, 10, 2, 4, 2, 10, 2, 6, 4, 2, 4, 6, 6, 2, 6, 6, 6, 4, 6, 8, 6, 4, 2, 4, 8, 10, 8, 4, 6, 2, 6, 6, 4, 2, 4, 14, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 10, 2, 10, 8, 6, 4, 8, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 6, 6, 2, 6, 6, 4, 2, 10, 2, 10, 2, 4, 6, 6, 2, 6, 4, 2, 10, 6, 2, 6, 6, 6, 4, 6, 12, 2, 4, 2, 12, 6, 4, 6, 2, 4, 8, 12, 4, 2, 4, 8, 6, 4, 2, 4, 2, 10, 2, 10, 8, 4, 6, 2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 10, 6, 8, 6, 4, 2, 4, 14, 4, 6, 2, 4, 6, 2, 6, 6, 6, 10, 2, 6, 4, 2, 4, 12, 12, 2, 4, 2, 10, 2, 6, 6, 4, 6, 6, 2, 10, 2, 6, 4, 14, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 16, 2, 16} $p = \{2, 3, 5, 7, 11, 13, 17\} \rightarrow 92160\quad\mathrm{cycle}$ $p = \{2, 3, 5, 7, 11, 13, 17, 19\} \rightarrow 1658880\quad\mathrm{cycle}$ 5. AKS素数判定法を使う 省略。こちらの記事が詳しく解説されています。(決定的アルゴリズム) 5. ミラーラビン素数判定法を使う 確率的アルゴリズム。100%ではないが素数かどうかを判定できる。約$n\lt 3.317\times10^{24}$までなら決定的なアルゴリズムであることが分かっている。 kを使う場合 誤差は$4^{-k}$。$k=20$であれば$\sim9.094\times 10^{-13}$ import random def mr_test_normal(n: int): if n == 2: return True elif n < 2 or n%2 == 0: return False d = n - 1 s = 0 while d & 1 == 0: s += 1; d >>= 1 k = 20 for _ in range(k): a = random.randint(2, n-1) x = pow(a, d, n) # a^d mod n if x == 1 or x == n-1: continue for _ in range(s): x = pow(x, 2, n) if x == n-1: # a^(d*2^r) mod n break else: return False return True kを使わない場合 def mr_test(n: int): if n == 2: return True elif n < 2 or n%2 == 0: return False assert n < int(3.317 * 10**24) d = n - 1 s = 0 while d & 1 == 0: s += 1; d >>= 1 it = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41) for a in it: if a == n: return False a %= n x = pow(a, d, n) # a^d mod n if x == 1 or x == n-1: continue for _ in range(s): x = pow(x, 2, n) if x == n-1: # a^(d*2^r) mod n break else: return False return True もうすこし厳密 def select_base(f): def wrap(n) -> bool: if n < 2047: bases = (2,) elif n < 1373653: bases = (2, 3) elif n < 9080191: bases = (31, 73) elif n < 25326001: bases = (2, 3, 5) elif n < 3215031751: bases = (2, 3, 5, 7) elif n < 350269456337: bases = (4230279247111683200, 14694767155120705706, 16641139526367750375) elif n < 55245642489451: bases = (2, 141889084524735, 1199124725622454117, 11096072698276303650) elif n < 7999252175582851: bases = (2, 4130806001517, 149795463772692060, 186635894390467037, 3967304179347715805) elif n < 585226005592931977: bases = (2, 123635709730000, 9233062284813009, 43835965440333360, 761179012939631437, 1263739024124850375) elif n < 2**64: bases = (2, 325, 9375, 28178, 450775, 9780504, 1795265022) elif n < 3317044064679887385961981: bases = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37) else: raise OverflowError return f(n, bases) return wrap @select_base def mr_test(n: int, mr_bases: tuple): if n == 2: return True elif n < 2 or n%2 == 0: return False d = n - 1 s = 0 while d & 1 == 0: s += 1; d >>= 1 for a in mr_bases: if a == n: return False a %= n x = pow(a, d, n) # a^d mod n if x == 1 or x == n-1: continue for _ in range(s): x = pow(x, 2, n) if x == n-1: # a^(d*2^r) mod n break else: return False return True (Wikipediaをもとに作成。) おわりに こんな感じでいろんなやり方があります。確実に高速で素数か判定したいのであれば、$10^{24}$までならミラーラビン、それ以上ならAKS、なんて使い方もいいかもしれません。