無心になったが、とりあえず サンプルも読む。 なるほど。 例えば** 1 2 3 4 が与えられたとすると **1x2+(1+2)x3+(1+2+3)x4とならないか？ 後は表現方法が問題だ。 だがしかし、tle & wa ans.py N = int(input()) A = list(map(int,input().split())) score = A[0]*A[1] for i in range(2,N): #計算量 O(N) score += sum(A[:i])*A[i] #計算量 O(N)...合算してほぼ O(N^2) score %= 7+10**9 print(score) 配列の sum の計算量は配列長さに依存する。 なので最終的には、それだけで O(N) となりうる。 っということで書き換えたら、通った。 ans_r1.py N = int(input()) A = list(map(int,input().split())) X = sum(A) score = 0 for i in range(N-1,0,-1): score += (X-A[i])*A[i] score %= 7+10**9 X -= A[i] print(score)

1.はじめに 久しぶりにロータリーエンコーダを使ってみたけど、構造を明確に理解していなかったことに気付いたので分解して勉強。 2.分解してみた 今回分解したロータリーエンコーダは、回転軸がプッシュスイッチになっているものです。エンコーダー検出用の3pin(A,B,GND)の他に、スイッチ用の2pin(SW,GND)がついています。 分解すると3階建て構造になっていました。 1階層目はプッシュスイッチの接点。2階層目はロータリーエンコーダの接点、3階層目はクリック感を出すためのばねがしこまれていました。 2-1.1階層目 一階層目をばらすと、プッシュスイッチ接点を、上部から貫通してきているエンコーダ軸が押し込む構造になっていました。スイッチは接点が同時にばねになっている一般的なタクトスイッチと同じ構造でした。 2-2.2階層目 2階層目は、エンコーダーの接点があります。放射状に延びているのがAstateとBstateの接点です。扇状に配置されているのがGNDにつながる接点です。これらの接点を接続するのが上面に存在するコンタクトです。コンタクトは一枚のリングから4本の接点が伸びている構造で、A,B各接点とGNDを適宜スイッチする構造です。 AsteteとBstateは1/4周、GNDは1/2周の範囲にしか配置されていませんが、コンタクトの接点は4ヶ所（赤矢印部）存在するため、全周にわたって接点が存在しているのと同等の働きをします。 2-3.3階層目 3階層目には板ばねが仕込まれていて、回転する際にコンタクト裏面の溝と接触してクリック感を出しています。 3.まとめ ロータリーエンコーダーの仕組みについては何となく知っていましたが、分解してより理解が深まりました。特に。AstateとBstateを同じ面に1/4周ずつ配置し、コンタクトの接点4個を使って各ステートが一周配置されているのと同等の処理を行うという構造には目からウロコの思いでした。 コンタクト自体は薄い鉄板状の端子の弾性を使って接点に押し付ける構造だったため、チャタリングは避けられない仕組みです。高級品はこのあたりの部品が良くなってチャタリングを防止しているのでしょうか。 前回の記事ではソフトでチャタリング対策を行いましたが、適当なコンデンサを端子間に挟んでハード面でのノイズ対策する方がスマートだなと思えたことに電子工作歴2年の成長を感じています。

ここでは、csvファイルの読み込みを説明します。 Pythonを学習していく上でデータインポートは必須であり、print()文やfor文を覚えたら、次に知りたいのが、データファイルのインポートのはずです。 トピック 1. pandasでcsvファイルを読み込む: pd.read_csv() 2. pandasで行、列の範囲を指定して csvを読み込む 2.1 日本語を含んだファイルの場合: pd.read_csv(encoding="shift_jis") 2.2 範囲を指定する: pd.read_csv(skiprows=, skipfooter=, usecols=) 3. 日付をindexとして読み込む: pd.read_csv(parse_dates=True) まとめ ポイント整理 参考ページ https://shilabo.com/python/web_self/ https://note.nkmk.me/python-pandas-read-csv-tsv/ https://qiita.com/sakabe/items/ae1fa47a58c796006627 読み込むファイルはcsv形式(.csv)が好んで用いられます。csvは書式がなく、扱いやすいのです。 Excel(.xlsx)を読み込む事も多いですが、まずはcsv形式を習得する事が重要です。 csvが読み込めれば、Excelもほぼ同じ手順でできます。 後で見るように、Pythonは日本語がひどく苦手です。インポートするデータは英語だけにしておくべきです。 データのインポートが出来ないと、何も出来ないに等しいです。 尚、当サイトではよく使う部分のみを説明し、詳細は他のHPに譲ることにします。 例えば、以下のサイトが非常に役に立ちます。極めて質の高い無料のサイトです。 1. pandasでcsvファイルを読み込む csvの読み込みは、pandasのpd.read_csv() を使います。 これを使う為に最初にimport pandas as pdと書きます。 読み込んだデータは基本DataFrameとなります。 ここでは、以下のようなcsvファイルを読み込みます。 データはA列が行の名前、B～D列がカラム名（列名）です。 通常、このようにデータには、行の名前(11, 12, 13)と列名（clm1, clm2, clm3）が付いています。列名はヘッダー(header)と呼ばれます。 これを、行の名前が(11, 12, 13)となるように、読み込もうと思います。DataFrameとして読み込んだ、行の名前のことをindexといいます。 以下のリンク（GitHub : 外部サイト）を右クリックして「名前を付けてリンク先を保存」か、一度開いてメモ帳にコピペして保存してください。その際に、拡張子が自動で「.txt」となりますが、「.csv」とタイプして保存してください。 a001_001a.csv ID,clm1,clm2,clm3 11,1863,7911,2634 12,1850,8000,2623 13,1853,7980,2578 Pythonでは pd.read_csv(‘ファイルのパス’)で読み込みます。 ファイルのパスは、例えば、Windowsのドキュメントに保存した場合、 C:\Users\ユーザ名\Documents となります。以下はユーザー名が “shilabo” で、ドキュメントフォルダの中に “SHiLABO_python” というフォルダを作った場合です。 パスを指定する場合”r”を頭に書くと\を\と書く必要がなくなり便利です。Windowsのエクスプローラーからパスをコピーして、そのままペーストするだけよくなります。 import pandas as pd df1a = pd.read_csv(r'C:\users\shilabo\documents\price1_01.csv') print(df1a) # ID clm1 clm2 clm3 #0 11 1863 7911 2634 #1 12 1850 8000 2623 #2 13 1853 7980 2578 あるいは、パスとファイル名を分けた方がいい場合があります。例えば、ループ(for文)で同じパスの複数のファイルを読み込みたい場合です。 その為、実務では以下のように書くことがあります。 p1 = r'C:\Users\shilabo\Documents\SHiLABO_python' f1 = 'a001_001a.csv' p1f1 = p1 +'\\'+f1 df1b = pd.read_csv(p1f1) print(df1b) # ID clm1 clm2 clm3 #0 11 1863 7911 2634 #1 12 1850 8000 2623 #2 13 1853 7980 2578 結果はどちらも同じです。 ただし、これでは目的が達成されていません。indexが自動で0,1,2と振られ、A列の11, 12, 13はデータとして扱われています。本来はA列の11, 12, 13をindexにしたかったはずです。 そこで、オプションのindex_col=0, header=0を使います。デフォルトはheader=0なので、これは省略できます。 df1c =pd.read_csv(p1f1, index_col=0, header=0) print(df1c) # clm1 clm2 clm3 #ID #11 1863 7911 2634 #12 1850 8000 2623 #13 1853 7980 2578 indexはカラム名(列名)で指定することもできます。以下は、’clm1’をindexにした場合です。 df1d =pd.read_csv(p1f1, index_col='clm1', header=0) print(df1d) # ID clm2 clm3 #clm1 #1863 11 7911 2634 #1850 12 8000 2623 #1853 13 7980 2578 2. pandasで行、列の範囲を指定してcsvを読み込む 2.1 日本語を含んだファイルの場合 しかし、現実には、このように元がきれいなデータである事は少なく、多くの場合、次のように、余計なものが含まれています。（GitHubのリンクa001_001b.csvをクリックすると形が見えます） しかも、日本語の嵐。Pythonは日本語が苦手です。 欲しいのは、セル「B4:E7」の範囲だけです。 部署,金融商品開発部,,, 作成日,2021/3/24,,, ,,,, ,ID,A,B,C ,11,1863,7911,2634 ,12,1850,8000,2623 ,13,1853,7980,2578 ,,,, *1: Aの参照期間は…。,,,, *2: Bは〇〇の配当あり。,,,, これもpandasのpd.read_csv()で読み込みますが、日本語を含む場合、オプションにencoding='shift_jis' を指定します。 “shift”と”jis”の間はアンダースコア(_)であり、ハイフン(-)ではありません。 これを間違い、データが読めないと悩む人が多いです。因みに、他に英語を読み取るencoding = 'utf-8'がありますが、こちらはハイフン(-)です。 とりあずencoding='shift_jis'を指定して実行すると、以下のようになります。 f2 = 'a001_001b.csv' p1f2 = p1 + '\\' + f2 # p1 = r'C:\Users\shilabo\Documents\SHiLABO_python' df2a = pd.read_csv(p1f2, encoding='shift_jis') # encoding='shift_jis'　が必要 print(df2a) # 部署 金融商品開発部 Unnamed: 2 Unnamed: 3 Unnamed: 4 #0 作成日 2021/3/24 NaN NaN NaN #1 NaN NaN NaN NaN NaN #2 NaN ID A B C #3 NaN 11 1863 7911 2634 #4 NaN 12 1850 8000 2623 #5 NaN 13 1853 7980 2578 #6 NaN NaN NaN NaN NaN #7 *1: Aの参照期間は…。 NaN NaN NaN NaN #8 *2: Bは〇〇の配当あり。 NaN NaN NaN NaN 読み込めていますが、csv上で空白も読み込み(NaN)、不必要なデータ（日本語部分）も読み込んでいます。これは意図とは異なります。データだけが必要なのです。むしろ、データ以外は取り込んではいけない、と言った方がいいでしょう。 2.2 範囲を指定する そこで範囲を指定する為に、以下のオプションを指定します。 skiprows=3 : 最初の3行はスキップ(1～3行目) skipfooter=3 : 下から3行はスキップ(8～10行目） usecols=range(1,5) : 列の指定(Python上での1列～4列目 = csv上では2列目(B列)から5列目(F列)) skipfooter=でWarningが出た場合、更にオプションで engine='python'を指定します。 ここで、usecols=について説明します。 range(1,5)と書くと、0始まりのPythonから見て、1, 2, 3, 4列目となります。 つまり0列目はスキップし、1, 2, 3, 4列目を読み、5列目は含みません。 csv上では普通1, 2, 3…列と数えるし、終わりの5も含むか含まないか、この感覚が難しいです。 実際には、print()で意図通りか必ず確認します。 また、index_col=0はusecolsで切り取った後の、Pythonの0列目をindexにする、という意味です。 同様に、header=0は、skiprowsで切り取った後の、Pythonの0行目をヘッダーにする、という意味です。 Pythonは基本0で始まります(1ではない）。 df2b= pd.read_csv(p1f2, encoding='shift_jis', skiprows=3, skipfooter=3, usecols=range(1,5), index_col=0, header=0, engine='python') #engine='python'で実行できる # print(df2b) # A B C #ID #11 1863 7911 2634 #12 1850 8000 2623 #13 1853 7980 2578 usecols の範囲はリストで指定することもできます。 この例の場合、usecols=[1,2,3,4]。ここでも、0始まりを意識しましょう。 df2c= pd.read_csv(p1f2, encoding='shift_jis', skiprows=3, skipfooter=3, usecols=[1,2,3,4], index_col=0, header=0, engine='python') print(df2c) # A B C #ID #11 1863 7911 2634 #12 1850 8000 2623 #13 1853 7980 2578 3. 日付をindexとして読み込む 日付がある時系列データを読み込む場合、”型”をdatetime型として読み込むと便利です。 index_col=でindexとするカラム名（列名）を指定し、parse_dates=Trueのオプションを指定します。 以下のcsvファイルを読み込んでみます。(GitHubのa001_001c.csvのリンクをクリックすると形が見えます。) ValDate,clm1,clm2,clm3 1999/1/1,1863,7911,2634 1999/2/12,1850,8000,2623 1999/3/23,1853,7980,2578 2005/1/23,1858,7905,2578 2005/2/12,1873,7952,2578 2005/3/1,1844,7820,2599 2021/1/12,1820,7480,2556 2021/2/23,1885,7349,2557 2021/3/1,1879,7199,2569 f3 = 'a001_001c.csv' p1f3 = p1 + '\\' +f3 #p1　=r'C:\Users\shilabo\Documents\SHiLABO_python' df3 = pd.read_csv(p1f3, index_col='ValDate',parse_dates=True) print(df3) # clm1 clm2 clm3 #ValDate #1999-01-01 1863 7911 2634 #1999-02-12 1850 8000 2623 #1999-03-23 1853 7980 2578 #2005-01-23 1858 7905 2578 #2005-02-12 1873 7952 2578 #2005-03-01 1844 7820 2599 #2021-01-12 1820 7480 2556 #2021-02-23 1885 7349 2557 #2021-03-01 1879 7199 2569 indexをdatetime型にすると、df.index.year == , df.index.month == , df.index.day ==と書くことで、年、月、日でフィルターをかけることができます。 以下は、1999年だけ、2月だけ、23日だけ、の例です。 print(df3[df3.index.year == 1999]) # clm1 clm2 clm3 #ValDate #1999-01-01 1863 7911 2634 #1999-02-12 1850 8000 2623 #1999-03-23 1853 7980 2578 print(df3[df3.index.month == 2]) # clm1 clm2 clm3 #ValDate #1999-02-12 1850 8000 2623 #2005-02-12 1873 7952 2578 #2021-02-23 1885 7349 2557 print(df3[df3.index.day == 23]) # clm1 clm2 clm3 #ValDate #1999-03-23 1853 7980 2578 #2005-01-23 1858 7905 2578 #2021-02-23 1885 7349 2557 年、年月の指定ならば、以下のようにdf['1999']やdf['2021-02']でも可能です。 print(df3['1999']) # clm1 clm2 clm3 #ValDate #1999-01-01 1863 7911 2634 #1999-02-12 1850 8000 2623 #1999-03-23 1853 7980 2578 print(df3['2021-02']) # clm1 clm2 clm3 #ValDate #2021-02-23 1885 7349 2557 indexをdatetime型にした場合の使い方は、次のサイトが詳しいです。 まとめ pd.read_csv()にいくつもオプションを指定すれば、不要なものを含むcsvファイルもデータとして読み込めます。しかし、とても面倒です。 現実には難しいですが、出来れば、読み込むcsvファイルの不要部分を削除し、シンプルなデータ形式にしておくのがいいです。その方が時間が省略できます。 また、日本語が含まれると、エラーの原因になり、その解決に時間が掛かるので、可能ならば日本語は避けた方がいいです。 ポイント整理 pd.read_csv()を使う。 日本語を含む場合、encoding='shift_jis'をオプションに指定する。 skiprows=, skipfooter=で上下の範囲（行）を指定できる。 skipfooter=でエラーが出たらengine='python'を指定する。 usecols=で左右（列）の範囲を指定できる。 index_col=でindexにする列を指定できる。 indexをdatetime型で読み込む場合、parse_dates=Trueを指定する。 参考ページ 実践で使っている技術が掲載されています。書籍のように構成されています。 ここに掲載したコードが全て掲載されています。 記事内で紹介したページです。

投稿日：2021/09/15 更新日：2021/09/16 (drawCircle関数、Mysimulation、Mysimulation2関数の修正) はじめに Pythonで結晶を描画するプログラムを作るシリーズの第2回目です。今回はCIFデータ内で使われている分極座標 (fractional coordinate system) (文献によっては原子座標、規格化座標などと呼ばれる)と直交座標系の関係性とその分極座標上のアフィン変換について説明していきます。今回は座学多めで、感覚的に対称操作というものを理解する事を目的とします。 Pythonで結晶を描画するプログラムを作る~単位格子を描く~ 1. 分極座標 (fractional coordinate system) 第1回目の記事でみたように6つの結晶パラメータから大きさや形が異なる単位格子が生成されます。その内、結晶をなすものは7つの晶系に分類され、結晶はそのどれかに必ず分類されます。これらを統一の規格で評価する座標系の1つが分極座標です。分極座標とは単位格子を単位立方体に変換した座標系で単位格子内の相対位置を表します。つまり、単位格子頂点の位置ベクトル$\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}$を変換$F$を用いて F\boldsymbol{A}=\boldsymbol{a}=\left(\begin{array}{c} 1\\ 0\\ 0 \end{array}\right), F\boldsymbol{B}=\boldsymbol{b}=\left(\begin{array}{c} 0\\ 1\\ 0 \end{array}\right),F\boldsymbol{C}=\boldsymbol{c}=\left(\begin{array}{c} 0\\ 0\\ 1 \end{array}\right) とした変換した座標系が分極座標系です。ここで、$\boldsymbol{a},\boldsymbol{b},\boldsymbol{c}$は分極座標系における単位格子頂点の位置ベクトルです。この変換$F$を用いて、単位格子の実際の形を描画する直交座標系(実格子座標系)上の点$\boldsymbol{X}=(X,Y,Z)'$から分極座標上の点$\boldsymbol{x}=(x,y,z)'$への変換は以下のように与えられます： \boldsymbol{x}=F\boldsymbol{X}=\left(\begin{array}{ccc} 1/A_{x} & -B_{x}C_{z}/V & (B_{x}C_{y}-B_{y}C_{x})/V \\ 0 & 1/B_{y} & -A_{x}C_{y}/V \\ 0 & 0 & 1/C_{z} \end{array}\right)\left(\begin{array}{c} X\\ Y\\ Z \end{array}\right)\tag{1} $A_{x},B_{x},\dots,V$の具体的な形は第一回目をご覧ください。この変換$F$は扱う結晶の晶系が立方晶($a=b=c,\alpha=\beta=\gamma=90^\circ$)のとき、 F_{cubic}=\left(\begin{array}{ccc} 1/a & 0 & 0 \\ 0 & 1/a & 0\\ 0 & 0 & 1/a \end{array}\right)=\frac{1}{a}I_{3} です。立方晶の角度は全て$90^\circ$なので単純に各軸に対して規格化を施すというような処理になっているのが見て取れるかと思います。CIFデータに記述される座標情報は全てこの分極座標系で記述されています。つまり、Pythonで結晶を描画するプログラムを作るためにはこの変換の逆変換$F^{-1}$が必要であり、下記のように表されます。 \boldsymbol{X}=F^{-1}\boldsymbol{x}=\left(\begin{array}{ccc} A_{x} &B_{x} & C_{x} \\ 0 & B_{y} & C_{y} \\ 0 & 0 & C_{z} \end{array}\right)\left(\begin{array}{c} x\\ y\\ z \end{array}\right)\tag{2} この変換$F^{-1}$は立方晶のとき、 F_{cubic}^{-1}=\left(\begin{array}{ccc} a & 0 & 0 \\ 0 & a & 0\\ 0 & 0 & a \end{array}\right)=aI_{3} です。復元するなら各軸を$a$倍しようという事です。これらの変換$F$と$F^{-1}$の間に必ず成り立つ性質として以下があります。 FF^{-1}=F^{-1}F=I_{3} 2. ザイツの記法によるアフィン変換の表現 CIFデータの"_symmetry_equiv_pos_as_xyz"は簡単に言えば分極座標系(xyz)における原子位置が記述されています。厳密な定義については次回取り扱います。SrTiO3では下記のように書かれています。 _symmetry_equiv_pos_as_xyz 1 'x, y, z' 2 '-x, -y, -z' 3 '-y, x, z' 4 'y, -x, -z' 5 '-x, -y, z' 6 'x, y, -z' 7 'y, -x, z' 8 '-y, x, -z' 分極座標上の原子位置は上記のように記述されていますが、これはアフィン変換と密接な関わりがあります。アフィン変換とは位置$\boldsymbol{x}=(x,y,z)'$にダミー次元を加えたアフィン空間上の点$\tilde{x}=(x,y,z,1)'$に対する変換であり、この変換をザイツの記法で表記すると以下のように記述されます。 (A|\boldsymbol{b})=\left(\begin{array}{cc} A & \boldsymbol{b}\\ 0 & 1 \end{array}\right)\tag{3} この行列は拡大行列とも呼ばれます。ここで、$A=(a_{ij})\quad (i,j=1,2,3)$は恒等操作$E$、反転$i$、$n$回回転$C_{n}$、鏡映$\sigma$、回映$S$、回反$C_{\bar{n}}$といった対称操作を表す$3\times 3$行列であり、$\boldsymbol{b}=(b_{i})\quad (i=1,2,3)$は平行移動を表し$3\times 1$の列ベクトルです。$A$と$\boldsymbol{b}$について整理すると、$A$は原点を中心とした変換であり、$\boldsymbol{b}$は原点を移す変換です。例として、$\boldsymbol{x}=(x,y,z)'$を$z$軸の周りに$90^\circ$回転($C_{4}(z)$)させ、$\boldsymbol{b}=(1,1,1)'$と平行移動させる事を考えてみます。まずはザイツの記法を使わすに計算してみます。 \begin{eqnarray} C_{4}(z)\boldsymbol{x}+\boldsymbol{b}&=&\left(\begin{array}{ccc} \cos 90^\circ & -\sin 90^\circ & 0\\ \sin90^\circ & \cos 90^\circ & 0 \\ 0 & 0 & 1 \end{array}\right)\left(\begin{array}{c} x\\ y\\ z \end{array}\right)+\left(\begin{array}{c} 1\\ 1\\ 1 \end{array}\right)\\ &=& \left(\begin{array}{ccc} 0 & -1 & 0\\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right)\left(\begin{array}{c} x\\ y\\ z \end{array}\right)+\left(\begin{array}{c} 1\\ 1\\ 1 \end{array}\right)\\ &=& \left(\begin{array}{c} -y+1\\ x+1\\ z+1 \end{array}\right) \end{eqnarray} これより、$\boldsymbol{x}=(x,y,z)'$は$\boldsymbol{x}'=(-y+1,x+1,z+1)'$に移動する事がわかります。ザイツの記法では計算すると以下が得られます。 \begin{eqnarray} \left(\begin{array}{cc} C_{4}(z) & \boldsymbol{b}\\ 0 & 1 \end{array}\right)\tilde{x}&=&\left(\begin{array}{cccc} \cos 90^\circ & -\sin 90^\circ & 0 & 1\\ \sin90^\circ & \cos 90^\circ & 0 & 1\\ 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 1 \end{array}\right)\left(\begin{array}{c} x\\ y\\ z\\ 1 \end{array}\right)\\ &=& \left(\begin{array}{c} -y+1\\ x+1\\ z+1\\ 1 \end{array}\right) \end{eqnarray} この記法でも同様の点に移動している事がわかります。ザイツの記法のおいしい点は対称操作と平行移動を含めた操作を1つの線形変換のように取り扱える事にあります。さらにダミー次元はどんな変換を考えても必ず$1$になるので計算結果からこの次元を除けば直ぐにどの点に移動したかわかります。余談ですがニューラルネットワークなどはアフィン変換や活性化関数などを組み合わせたものであるため、ザイツの記法を使うとよりシンプルに表現できます。計算結果の解釈について考えていきます。分極座標系における単位格子の頂点は整数$n_{1},n_{2},n_{3}$を用いて、 \boldsymbol{r} = n_{1}\boldsymbol{a}+n_{2}\boldsymbol{b}+n_{3}\boldsymbol{c}:= (n_{1},n_{2},n_{3}) と表現されます。$(0,0,0)$,$(0,0,1)$,$(0,1,0)$,$(1,0,0)$,$(1,1,0)$,$(0,1,1)$,$(1,0,1)$,$(1,1,1)$の8つの頂点に囲まれた単位格子の領域を$U_{000}$としたとき、実格子空間上の$U_{000}$に対応した単位格子内の原子位置は必ず$0\leq x,y,z\leq 1$になければいけません。$(x+1,y+1,z+1)'$という点を考えたとき、原点を$(1,1,1)$に移した単位格子$U_{111}$上の$(x,y,z)'$という点にあるという事が言えます。つまり、$\boldsymbol{x'}$は$x$方向に$1$格子移動、$y$方向に$1$格子移動、$z$方向に$1$格子移動した単位格子$U_{111}$を中心とした$(-y,x,z)'$の点を表している事と同義と言えます。結晶とは長期的な周期構造を持つ原子の集まりであるため、$(-y,x,z)'$が他の単位格子に存在するのであれば当然ながら$U_{000}$にもこの点が存在する必要があります。つまり、$(x,y,z)'=(x\pm n_{1},y\pm n_{2},z\pm n_{3})'$という関係を持ちます。もし、これをCIFデータの"_symmetry_equiv_pos_as_xyz"に書くのであれば下記になります(これを要素に持つ空間群はないのであくまでも例だと考えて下さい)。 _symmetry_equiv_pos_as_xyz 1 'x, y, z' 2 '-y, x, z' 3. 任意の軸周りの回転 ザイツの記法による拡大行列を構成する要素は前述したように対称操作$A$と平行移動$b$です。ここでは対称操作$A$における回転$C_{n}$について取り扱っていきます。まずは有名な以下の定理を紹介します。 ロドリゲスの回転公式 (Rodrigues' rotation formula) 任意の回転軸$p=(p_1,p_2,p_3)'$まわりに$\theta$回転させる回転行列$R_{p}(\theta)$は R_{p}(\theta)= \left(\begin{array}{ccc} c + p_{1}^2(1-c) & p_{1}p_{2}(1-c)-p_{3}s & p_{1}p_{3}(1-c)+ p_{2}s \\ p_1 p_2 (1-c)+p_3 s & c+p_2^2 (1-c) & p_2 p_3 (1-c) -p_1 s \\ p_1 p_3(1-c)-p_2 s & p_2 p_3 (1-c)+ p_1 s & c + p_3^2 (1-c) \end{array}\right) である。ここで、$s=\sin\theta,c=\cos\theta$を表す。 特に、$p=(0,0,1)'$かつ$\theta = 90^\circ$のとき、$R_{p}(90^\circ)=C_{4}(z)$です。 ここまでの勉強した所を使って点の移動の軌跡をみるプログラムを作っていきましょう。まず、注意点としてmathライブラリにおける三角関数は誤差が出てしまうという事が挙げられます。例えば、$\cos(\pi/2)=0$という関係を素直にコードで書けば"math.cos(math.pi/2)"なのですが、出力結果は6.123233995736766e-17です。非常に小さな桁で誤差が出るため、厳密な$0$になりません。簡単な対処法は出力結果を四捨五入して$0$とするという方法です。そもそも、小数点十数位の変化がピクセルをまたいで変化して見える事がないので適当な桁で処理してしまうのが良いかと思います。また、numpyは1次元配列はベクトルの転置をとったものになっているので、行列積を考えるときは転置".T"をとる必要があります。これらを踏まえて、ロドリゲスの回転公式に代入しこの回転行列と平行移動ベクトルから式(3)の拡大行列を出力、この行列と移動させたい点の行列積を計算し、移動先の点を求め、図示という手順で処理していきたいと思います。 import numpy as np import math from decimal import Decimal, ROUND_HALF_UP import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 関数 def deg2rad(deg): """ deg→rad """ return deg*(math.pi/180) def Round(vec,digit): """ 入力された値をdigitsの桁で四捨五入する """ num = len(vec) Rounding = [] for i in range(num): a = Decimal(str(vec[i])).quantize(Decimal(digit),rounding=ROUND_HALF_UP) a = float(a) Rounding.append(a) return Rounding def Rodrigues(vecP,theta,digit='0.01'): """ ロドリゲスの回転行列。 ------------------------------ vecP : array-like 回転軸を表すベクトル theta : float 反時計周りにθ回転させる。θはdegreeで入力する。 """ theta = deg2rad(theta) p1,p2,p3 = vecP c = math.cos(theta) s = math.sin(theta) # 回転行列の成分 R11 = c+(1-c)*p1**2 R12 = p1*p2*(1-c)-p3*s R13 = p1*p3*(1-c)+p2*s R21 = p1*p2*(1-c)+p3*s R22 = c+(1-c)*p2**2 R23 = p2*p3*(1-c)-p1*s R31 = p1*p3*(1-c)-p2*s R32 = p2*p3*(1-c)+p1*s R33 = c+(1-c)*p3**2 return np.array([Round([R11,R12,R13],digit) ,Round([R21,R22,R23],digit) ,Round([R31,R32,R33],digit)])+0 def seitzMtx(mtx,vecB): """ ザイツの記法による拡大行列。 --------------------------------- mtx：array 3×3の行列 vecB：array-like 平行移動を表すベクトル """ seitzMtx = np.insert(mtx,3,vecB,axis=1) # 4列目にvecBを加える seitzMtx = np.insert(seitzMtx,3,[0,0,0,1],axis=0) # 4行目に[0,0,0,1]を加える return seitzMtx def movePos(seitzMtx,x): """ ザイツ記法の拡大行列とxの行列積。 ----------------------------------- x：array-like xは実格子上の点とする """ xtilde = np.insert(x,3,1,axis=0).T return np.dot(seitzMtx,xtilde) def drawCircle(x_bef,vecP): """ rotMtxが回転する平面の図示。 """ p1,p2,p3 = vecP Theta = [i for i in range(0,360,1)] x = np.zeros(len(Theta)) y = np.zeros(len(Theta)) z = np.zeros(len(Theta)) vecP_unit = np.array([p1,p2,p3])/math.sqrt(p1**2+p2**2+p3**2) # vecPを単位ベクトルに変換 xc,yc,zc = np.dot(np.array(x_bef).T,vecP_unit.T)*vecP_unit # x_befからvecPの線上に垂直におろした点(xc,yc,zc) for i,theta in enumerate(Theta): x[i],y[i],z[i] = np.dot(Rodrigues(vecP_unit,theta),np.array(x_bef).T) return x,y,z,xc,yc,zc def Mysimulation(x_bef,vecP,vecB,theta,size=40,fontsize=15,fig_size1=8,fig_size2=8,color='r',arrow_length_ratio=0.1,lim=8): """ pos_bftからpos_aftに向かう矢印、x,y,z軸の描画 """ # ●計算 rotMtx = Rodrigues(vecP,theta) seitz = seitzMtx(rotMtx,vecB) x_aft = movePos(seitz,x_bef)[:-1] # 末尾はダミー次元なので使わない # ●描画 fig = plt.figure(figsize=(fig_size1,fig_size2)) ax = fig.add_subplot(111,projection='3d') # 移動前と移動後の点のプロットおよびその間に矢印を引く x,y,z = x_bef u,v,w = x_aft p1,p2,p3 = vecP ax.quiver(x,y,z,u-x,v-y,w-z,color=color,arrow_length_ratio=arrow_length_ratio) ax.scatter(x,y,z,color='b',s=size) ax.text(x,y,z,s=r'$x_{bft}$',fontsize=fontsize) ax.scatter(u,v,w,color='r',s=size) ax.text(u,v,w,s=r'$x_{aft}$',fontsize=fontsize) # x,y,z軸の追加 ax.quiver(-lim,0,0,2*lim,0,0,linestyle='-', arrow_length_ratio=0.1,color='black') # x軸 ax.quiver(0,-lim,0,0,2*lim,0,linestyle='-', arrow_length_ratio=0.1,color='black') # y軸 ax.quiver(0,0,-lim,0,0,2*lim,linestyle='-', arrow_length_ratio=0.1,color='black') # z軸 # x_befが回転する平面 ax.quiver(-lim*p1,-lim*p2,-lim*p3,2*lim*p1,2*lim*p2,2*lim*p3,linestyle='--', arrow_length_ratio=0.1,color='b') # 回転軸 r1,r2,r3,xc,yc,zc = drawCircle(x_bef,vecP) # x_befが回転する平面 ax.plot(r1,r2,r3,color='g',ls='--') ax.quiver(xc,yc,zc,x-xc,y-yc,z-zc,linestyle='-', arrow_length_ratio=0.1,color='g') # x_befの回転先 xr,yr,zr = np.dot(rotMtx,np.array(x_bef).T) ax.scatter(xr,yr,zr,color='g',s=size) ax.quiver(xc,yc,zc,xr-xc,yr-yc,zr-zc,linestyle='-', arrow_length_ratio=0.1,color='g') ax.quiver(x,y,z,xr-x,yr-y,zr-z,linestyle='--', arrow_length_ratio=0.1,color='r') ax.text(xr,yr,zr,s=r'$x_{rot}$',fontsize=fontsize) # 回転先から移動後のプロットまでに矢印を引く ax.quiver(xr,yr,zr,u-xr,v-yr,w-zr,linestyle='--', arrow_length_ratio=0.1,color='r') # 出力範囲 ax.set_xlim(-lim,lim) ax.set_ylim(-lim,lim) ax.set_zlim(-lim,lim) # 軸ラベル ax.set_xlabel('x',size=14) ax.set_ylabel('y',size=14) ax.set_zlabel('z',size=14) return rotMtx,seitz,x_aft # 回転軸ベクトルP, 点x, 平行移動vecB, 回転角thetaの指定 vecP = [0,0,1] # z軸 x_bef = [-0.5,-0.5,-1] vecB = [1,1,1] # b=(1,1,1)'の移動 theta = 90 _,_,_ = Mysimulation(x_bef,vecP,vecB,theta,lim=2) # plt.savefig('test.png',transparent=True,dpi=300) ここで、Rodrigues関数の戻り値に$+0$を付けているのですが、これは$-0$が出力されるのが避けられないのでこれを$0$に直す処理として加えています。上記を実行すると以下の図が得られます。今回は対称操作の1つである回転に限ってのプログラムを作成しましたが、seitzMtx関数の引数である"mtx"に恒等操作$E$、反転$i$、鏡映$\sigma$、回映$S$、回反$C_{\bar{n}}$を表す$3\times 3$行列を代入する事で様々な操作による変化を楽しむことができます。 最後に複数の点が対称操作と平行移動によってどう変化するのかを描画するプログラムを作成して終了したいと思います。 def movePoss(seitzMtx,pSet_bef): """ 点集合pSet_befをseitzMtxによってpSet_aftに移す ----------------------------------- pSet：array 点集合の座標 """ num = pSet_bef.shape[1] pSet_aft = np.zeros_like(pSet_bef) for i in range(num): xtilde = np.insert(pSet_bef[:,i],3,1,axis=0).T pSet_aft[:,i] = np.dot(seitzMtx,xtilde)[:-1] return pSet_aft def Mysimulation2(vecP,vecB,pSet_bef,theta,size=40,fontsize=15,fig_size1=8,fig_size2=8,color='r',arrow_length_ratio=0.1,lim=8): """ pSet_bef, pSet_aftの描画 """ # ●計算 rotMtx = Rodrigues(vecP,theta) seitz = seitzMtx(rotMtx,vecB) seitz_rot = seitzMtx(rotMtx,[0,0,0]) # 回転のみ seitz_trans = seitzMtx(np.eye(3),vecB) # 平行移動のみ pSet_rot = movePoss(seitz_rot,pSet_bef) pSet_trans = movePoss(seitz_trans,pSet_bef) pSet_aft = movePoss(seitz,pSet_bef) # ●描画 fig = plt.figure(figsize=(fig_size1,fig_size2)) ax = fig.add_subplot(111,projection='3d') # 移動前と移動後の点のプロット ax.scatter(pSet_bef[0],pSet_bef[1],pSet_bef[2],color='b',s=size,label=r'$x_{bef}$') ax.scatter(pSet_trans[0],pSet_trans[1],pSet_trans[2],color='g',s=size,label=r'$x_{trans}$') ax.scatter(pSet_rot[0],pSet_rot[1],pSet_rot[2],color='c',s=size,label=r'$x_{rot}$') ax.scatter(pSet_aft[0],pSet_aft[1],pSet_aft[2],color='r',s=size,label=r'$x_{aft}$') # x,y,z軸の追加 ax.quiver(-lim,0,0,2*lim,0,0,linestyle='-', arrow_length_ratio=0.1,color='black') # x軸 ax.quiver(0,-lim,0,0,2*lim,0,linestyle='-', arrow_length_ratio=0.1,color='black') # y軸 ax.quiver(0,0,-lim,0,0,2*lim,linestyle='-', arrow_length_ratio=0.1,color='black') # z軸 # 回転軸 p1,p2,p3 = vecP ax.quiver(-lim*p1,-lim*p2,-lim*p3,2*lim*p1,2*lim*p2,2*lim*p3,linestyle='--', arrow_length_ratio=0.1,color='b') # 出力範囲 ax.set_xlim(-lim,lim) ax.set_ylim(-lim,lim) ax.set_zlim(-lim,lim) # 軸ラベル ax.set_xlabel('x',size=14) ax.set_ylabel('y',size=14) ax.set_zlabel('z',size=14) plt.legend() return rotMtx,seitz,pSet_aft # 単位格子内の点集合を乱数で生成 np.random.seed(1) # シード値の固定 points = 20 # 生成する点の数 x = np.random.rand(points) # x座標 y = np.random.rand(points) # y座標 z = np.random.rand(points) # z座標 pSet_bef = np.array([x,y,z]) # 回転軸ベクトルP,平行移動vecB, 回転角thetaの指定 vecP = [0,0,1] vecB = [-1,-1,-1] theta = 90 # 描画 _,_,_=Mysimulation2(vecP,vecB,pSet_bef,theta,lim=2) #plt.savefig('test2.png',transparent=True,dpi=300) 上記のプログラムを実行する下図が得られます。この例では点の集合を乱数で作成していますが、前の記事扱った平行六面体の頂点に実行してみて下さい。

タイトルの通り、LeetCode 344の解き方についてシェアしてきます。 問題文はそのままコピーします、 Write a function that reverses a string. The input string is given as an array of characters s. Example 1: Input: s = ["h","e","l","l","o"] Output: ["o","l","l","e","h"] Example 2: Input: s = ["H","a","n","n","a","h"] Output: ["h","a","n","n","a","H"] Constraints: 1 <= s.length <= 105 s[i] is a printable ascii character. ま、簡単に言うと、string型の'hello'を'olleh'で変換することです。 各言語内の関数やモジュールを使えばええじゃんと思うかもしれないが、実はこれを実現するためにはいろいろな方法があります。 今回は2つの方法を紹介します、一つはタイトルの通り再帰で実現します。もう一つはTwo pointersで実現します。 Recursion solution Recursionの概念はわかると思いますが、ざっくりで言うと、 1. 関数内部にもう一度今所属している関数を呼び出す。 2．永遠に呼び出されたらキリがないので、必ずなにか終わりの条件を追加する。 実はこの問題に関して、Two pointersのやり方は多分理解しやすいが、Recursionで実現するのも面白いです。 では、大体のイメージ図からどのように再帰で実現したかをみてみましょう！ 解説： まず、1回目の再帰はhとoを同時に再帰します(多分理解にくいが、コード見ればわかると思います)。 2回目に入って、また再帰し、lのところで再帰が止まります。 止まったら、今回は元に戻します。 ここで、2回目の再帰の部分はeとlを交換し、1回目はhとoを交換します。 それで終了となります。 コード： def recursion(s, left, right): if left >= right: return recursion(s, left+1, right-1) s[left], s[right] = s[right], s[left] left >= rightが再帰の終了条件。 4行目の部分は再帰の部分で、毎回の再帰は左側を右に1 step,　右側は左に1 stepの形ですすみます。 再帰が終了した後(left < rightの時)、各文字を交換します。 Two pointers def recursion(s, left, right): while left < right: s[left], s[right] = s[right], s[left] left += 1 right -= 1 return s 考え方としては、再帰の解き方と似ていますが、実は違いますね。 参考

概要 自作モジュールを相対パスでimportしたく、sys.path.append(os.path.dirname(__file__)を入れたのに何故かimportができなかった。 結論 pipで入れていたパッケージと自作モジュールで利用していた名前が被っていたため、pipで入れたパッケージを読み込みにいってしまっていた。 環境 Windows10 miniconda 詳細 以下のような構成で自作モジュールを呼び出そうとしたところ、libs.utilsが見つからないとModuleNotFoundErrorが発生。 parent_dir/ |- my_modules/ | |- libs/ | |- utils.py | |- my_module.py |- controller.py controller.py ⇒ my_module.py ⇒ utils.pyと呼び出す流れ。 各モジュールのコード controller.py from my_modules.my_module import * module_func() my_modules\my_module.py from libs.utils import * def module_func(): print("This is module_func") utils_func() my_modules\libs\utils.py def utils_func(): print("This is utils_func") やりたいこと parent_dirにいる状態でpython controller.pyを実行 ⇒以下のような結果が欲しい。 AnacondaPrompt [自分のフォルダ]\parent_dir>python controller.py This is module_func This is utils_func しかし、結果は以下の通り。 AnacondaPrompt File "[自分のフォルダ]\parent_dir\my_modules\my_module.py", line 1, in <module> from libs.utils import * ModuleNotFoundError: No module named 'libs.utils' 原因と解決策 原因：コマンドから実行する場合、controller.pyを実行しているparent_dirが基準になるため、parent_dirからlibs.utilsを探しているが見つからない。 解決策：こちらの記事を参考に「モジュール探索パスを追加して絶対インポート」する方法を実践。 つまり、my_modules\my_module.pyにsys.path.append(os.path.dirname(__file__)を追加する。 my_modules\my_module.py # 追加分 import os import sys sys.path.append(os.path.dirname(__file__)) # 追加分ここまで from libs.utils import * def module_func(): print("This is module_func") utils_func() これにより、my_moduleフォルダもパスに追加されるため、相対パスとして読み込んでくれるはず！ 結果 AnacondaPrompt File "[自分のフォルダ]\parent_dir\my_modules\my_module.py", line 1, in <module> from libs.utils import * ModuleNotFoundError: No module named 'libs.utils' 変わらない！なぜ！ 現象の調査 いろいろ試してみると、自分が使っているminicondaの仮想環境でこの事象が発生することが判明。(新規作成したまっさらな仮想環境では事象が発生しない。) 試したこと 仮想環境にinstallした何かしらのパッケージが原因と考え、問題環境のinstallパッケージの一部を削除しながら原因のパッケージを特定することに。 おそらく他に効率が良いやり方があるのだろうが、知識のない私は地道な作業で調査。 ※以下のコマンドはこちらを参考にさせていただきました。 ①仮想環境の情報をyamlファイルに書き出す。 AnacondaPrompt conda env export > ファイル名.yml ②yamlファイルのうちdependenciesを2分探索の要領で分割 ③分割したyamlファイルを基に新規環境作成 AnacondaPrompt conda env create -n 新たな環境名 -f ファイル名.yml 特定 そしてついに突き止めた。 env.yaml dependencies: - pip=21.2.4=pyhd8ed1ab_0 - pip: - libs==0.0.10 このlibsが存在していると問題が発生する。 ここまでくれば原因は明確。 libs.utilsを見に行く際に、pipでinstallしているlibsを優先的に探していたということであった。 教訓 自作モジュールを使う場合は汎用的な命名を避ける。(my_libsのような名前にする。) condaとpip：混ぜるな危険←今回は直接関係なかったが、改めて確認するとpip installしているパッケージが多かったと反省。

うーん。さっぱり分からん。 20分悩んで降参。 なるほど。 a^2 <= N となるから a <= (N)^(0.5)が言える。。か。 想定される最大値 N と、想定される最小値 a^b(=a^2) の関係(不等式)から a の条件を導き出している。 Unexpressed.py N = int(input()) lis = [] for a in range(2,1+10**5): b = 2 while True: if a**b <= N: lis.append(a**b) b += 1 else: break lis = set(lis) #二重の値を省く処理 print(N - len(lis)) 勉強になりました。m(_ _)m

herokuの問題点 herokuは画像がアップロードできない仕様。 AWSなら良いのだけど、herokuの方が管理が楽で安い！ そして無料プランもある！ ってことでCloudinaryを使って実装してみる。 herokuにアドオンを追加 herokuにログインして自分のアプリに移動する。 「Configure Add-ons」→「Find more add-ons」 command + FでCloudinaryを検索。 Image Processingの「Cloudinary - Image and Video Management」内でアドオンを追加できます。 pip 画像を扱うにはPillowが必要なので入ってなければ入れておく。 terminal $ pip install Pillow 次にCloudinaryのパッケージを入れる。 terminal $ pip install django-cloudinary-storage #動画の場合は $ pip install django-cloudinary-storage[video] requirements.txtに書き込む。 terminal $ pip freeze > requirements.txt Djangoに設定する まずはINSTALLED_APPSに追加する。 setting.py #メディアファイルのみの場合 INSTALLED_APPS = [ #... 'django.contrib.staticfiles', 'cloudinary_storage', 'cloudinary', #... ] #メディアファイル以外がある場合 INSTALLED_APPS = [ #... 'cloudinary_storage', 'django.contrib.staticfiles', 'cloudinary', #... ] メディアファイルのみの場合 cloudinary_storageはdjango.contrib.staticfilesより下に置く。 メディアファイル以外の静的ファイルがある場合 django.contrib.staticfilesはcloudinary_storageより下に置く。 ※ 公式サイトを参照して下さい 次にCloudinaryのシークレットキーなどを設定。 キー等はCloudinaryのダッシュボードに書いてあります。 setting.py #setting.pyの最後に追加でOK CLOUDINARY_STORAGE = { 'CLOUD_NAME':'your_cloud_name', 'API_KEY':'env('CLOUDINARY_API_KEY')', 'API_SECRET':'env('CLOUDINARY_API_SECRET')' } シークレットキーなど公開するべきでないものは.envファイルに書いとく。 あとはパスの設定 setting.py MEDIA_URL = '/media/' DEFAULT_FILE_STORAGE = 'cloudinary_storage.storage.MediaCloudinaryStorage' これで画像がアップできるようになります。 モデルでのフィールド ImageFieldを使った通常通りの書き方でOKです！ models.py class TestModel(models.Model): name = models.CharField(max_length=10) image = models.ImageField(upload_to="image/", blank=True) テンプレート 他同様.image.urlでインスタンス化できます。 html <img src = "{{ test_model_instance.image.url }}"> ユーザーがアップロードした元のサイズで画像がダウンロードされるのでコントロールしたい場合は、 html {% load cloudinary %} <!-- 画像を挿入したい場所に --> {% cloudinary test_model_instance.image.name width=100 height=100 %} .image.urlでなく.image.nameなので注意。 これで100 × 100のサイズが読み込まれるようになります。 さらに、<img>を使いたい場合は絶対パスでURL内に設定を書き込めば良いです。 html <img src="https://res.cloudinary.com/<your_cloud_name>/image/upload/c_fill,h_100,w_100/v1/{{ test_model_instance.image.name }}"> 絶対パスの取得手順 適当な画像のURLを取得 画像のパスを取得 絶対パスの作成 オプションを設定 html <!-- CloudinaryのURLの構造 --> https://res.cloudinary.com/<cloud_name>/<asset_type>/<delivery_type>/<transformations>/<version>/<public_id>.<extension> 1,実際のURLはこんな感じかと https://res.cloudinary.com/<your_cloud_name>/image/upload/v1/<ディレクトリ名>/hoge_svqqsy 2,画像のパスを取得 {{ test_model_instance.image.name }}で<ディレクトリ>以降が取得できます。 3,絶対パスの作成 １と２を組み合わせる。 https://res.cloudinary.com/<your_cloud_name>/image/upload/v1/{{ test_model_instance.image.name }} 4,サイズ等オプションの設定 <transformations>の場所に画像のサイズ等の設定を書きます。 https://res.cloudinary.com/<your_cloud_name>/image/upload/h_150,w_150/v1/{{ test_model_instance.image.name }} Cloudinaryの便利機能 Cloudinaryにはサイズの変更だけでなく、様々な機能があります。 顔認証などもあり、全てを書くと膨大な量になるので一部だけ紹介。 /c_fill,h_150,w_150/ 150×150のサイズでいい感じにフィットさせる。 /c_fill,g_face,h_150,w_150/ 顔認証機能で顔にフォーカス 詳しく知りたい方は公式ページでご確認ください。 使われていない画像の削除 Cloudinaryには、未使用画像の一括削除という非常に便利な機能があります。 django-cloudinary-storageの場合管理コマンドは次の３つのみ。 collectstatic deleteorphanedmedia deleteredundantstatic 不要なメディアファイルの削除はdeleteorphanedmedia 使えるようにするためsetting.pyに記述。 setting.py CLOUDINARY_STORAGE = { #その他の設定 'EXCLUDE_DELETE_ORPHANED_MEDIA_PATHS':('path/', 'second-path/') } 引数に除外するパスを書きます。 上記の場合、path/とsecond-path/のメディアファイルは削除しない設定です。 全てを対象とする場合は引数無しでOK。 管理コマンドの実行 本番環境で実行するためheroku runを最初につけます。 terminal $ heroku run python manage.py deleteorphanedmedia 5 files will be deleted: - media/hoge1_g5bg8j - media/hoge2_h8roax - media/hoge3_uwhmgw - media/hoge4_svqqsy - media/hoge5_nqyljq If you are sure to delete them, please type 'yes': yesで実行すると使われていない画像が削除されます。楽ちんだ！ まとめると setting.py INSTALLED_APPS = [ #... 'django.contrib.staticfiles', 'cloudinary_storage', 'cloudinary', #... ] #色んな設定 MEDIA_URL = '/media/' DEFAULT_FILE_STORAGE = 'cloudinary_storage.storage.MediaCloudinaryStorage' CLOUDINARY_STORAGE = { 'CLOUD_NAME':'your_cloud_name', 'API_KEY':'env('CLOUDINARY_API_KEY')', 'API_SECRET':'env('CLOUDINARY_API_SECRET')' 'EXCLUDE_DELETE_ORPHANED_MEDIA_PATHS':(), } herokuはサーバー管理が楽なので本当に良いですね！ django-cloudinary-storageのREADME.mdで最新の情報を見れます。 こんなの作ってます オンラインでの仲間づくりに活用できる、心理学を利用したレジュメサービスです。 自分や色んな人の性格を見ることができますよ！ 無料なので使ってみてもらえると嬉しいです！（フィードバックも歓迎です！）

やりたいこと アクセスログをもとに機械学習で不正アクセスを予測する仕組み作る。 1.Splunkを使い、access.logをCSVに変換 2.CSVファイルをもとにJupyterLab上で教師データを作る 3.ディープラーニングを用いて、予測精度を確かめる　　⇦今回の内容 環境 ・JupyterLab　2.1.4 @ AWS Cloud9 ・Python 3 ・Splunk Enterprise　8.1 @ AWS EC2 別ノートで作業するので初期設定を行う。 初期設定 import numpy as np import matplotlib.pyplot as plt import pandas as pd # Keras from tensorflow import keras from tensorflow.keras.layers import Dense, Activation # データの分割 from sklearn.model_selection import train_test_split # MSE from sklearn.metrics import mean_squared_error # JupyterNotebook上でグラフを表示する設定 %matplotlib inline # DataFrameで全ての列を表示する設定 pd.options.display.max_columns = None データ取り込み 前回作成した特徴量データの確認 dataset = pd.read_csv('dataset.csv') dataset.head() 機械学習構築 目的変数(Y),説明変数(X)を設定する。 Y = np.array(dataset['clf']) X = np.array(dataset.iloc[:,:28]) 形状を確認する。 print(Y.shape) print(X.shape) # データの分割 X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.3, random_state=0) X_train, X_valid, Y_train, Y_valid = train_test_split(X_train, Y_train, test_size=0.3, random_state=0) それぞれの形状を確認 print("Y_train=", Y_train.shape, ", X_train=", X_train.shape) print("Y_valid=", Y_valid.shape, ", X_valid=", X_valid.shape) print("Y_test=", Y_test.shape, ", X_test=", X_test.shape) 機械学習モデルの設定と構築を行う。 model = keras.Sequential() model.add(Dense(8, activation='relu', input_shape=(28,))) model.add(Dense(32, activation='relu')) model.add(Dense(1, activation='sigmoid')) model.compile(optimizer = "rmsprop", loss='binary_crossentropy', metrics=['accuracy']) 学習前に、モデルの構造を確認する。 ※EC2のスペックに応じて変更 model.summary() 学習開始 モデルを用いて、機械学習開始。epochs=200もいらなかった。 %%time log = model.fit(X_train, Y_train, epochs=200, batch_size=128, verbose=True, callbacks=[keras.callbacks.EarlyStopping(monitor='val_loss', min_delta=0, patience=100, verbose=1)], validation_data=(X_valid, Y_valid)) 完了後、グラフで表示して確認。２００回もいらなかったかも。 plt.plot(log.history['loss'], label='loss') plt.plot(log.history['val_loss'], label='val_loss') plt.legend(frameon=False) # 凡例の表示 plt.xlabel("epochs") plt.ylabel("crossentropy") plt.show() 構築したモデルを使い予測を行う。 Y_pred = model.predict(X_test) Y_pred_cls = (Y_pred > 0.5).astype("int32") 形状を目的変数に合わせる。 Y_pred_ = Y_pred_cls.reshape(-1) モデルの評価 モデルの評価を行う。スコアは良好。 from sklearn.metrics import classification_report print(classification_report(Y_test, Y_pred_)) 以上。 参考 Scutum技術ブログ 異常検出アルゴリズム Isolation Forest

概要 こちらは前回Flask+スクレイピングで記事一覧画面を作ろう②の続きです。 テンプレートを作ろう まずは、前回のFlaskScraping/app.pyに必要なモジュールのインポートを行います。 差分は(-)と(+)で表します。 FlaskScraping/app.py # モジュールのインポート from flask import Flask, (+)render_template (+)import requests (+)from bs4 import BeautifulSoup # Flaskアプリのインスタンスを作成 app = Flask(__name__) ... requestsとBeautifulSoupはスクレイピングを行うためのモジュールですね。 render_templateは、簡単に言うと、対応するhtmlファイルを表示することができるという関数です。 後ほど使うので、頭の片隅に置いておいてください。 続いて、テンプレートを作っていきます。 FlaskScraipingフォルダ内にtemplatesフォルダを作成してください。 FlaskScraping 　| - templates そして、templatesフォルダ内にbase.htmlとindex.htmlを作りましょう。 FlaskScraping 　| - templates 　|　| - base.html 　|　| - index.html まずはbase.htmlから描いていきましょう。このファイルはhtmlファイルのベースとなるファイルです。つまり、htmlファイルに共通したルールなどをここでまとめて書いてしまおうということです。 コードは以下です。 templates/base.html <!DOCTYPE html> <html lang="ja"> <head> <meta charset="UTF-8"> <link href="https://cdn.jsdelivr.net/npm/bootstrap@5.0.0-beta1/dist/css/bootstrap.min.css" rel="stylesheet" integrity="sha384-giJF6kkoqNQ00vy+HMDP7azOuL0xtbfIcaT9wjKHr8RbDVddVHyTfAAsrekwKmP1" crossorigin="anonymous"> <title>Flask+スクレイピングで記事一覧画面を作ろう</title> </head> <body> {% block content %} {% endblock %} </body> </html> 今回、使用するUI(デザイン)はBootStrapです。 BootStrapは、CSSで良い感じの見た目にされたページ機能を提供してくれているものです。 これで、サイトが一気にオシャレになることでしょう。 チートシートに色々なコード例が載っています。 {% block content %}については、index.htmlで説明します。 とりあえず書いておきましょう。 では、index.htmlを書いていきましょう。 とりあえず、HelloFlask!と表示させてみましょう。 templates/index.html {% extends 'base.html' %} {% block content %} <p>Hello Flask!</p> {% endblock %} ここで、index.htmlにも{% block content %}というのがありますね。 これは、まず1行目の{% extends 'base.html' %}でbase.htmlを読み込みます。 そして、base.htmlタグのbodyタグの中にある{% block content %}をindex.htmlから呼び出すことで、base.htmlを継承することができます。 これにより、どのページからでも統一されたテンプレートを呼び出すことができるようになりました。 テンプレートを表示させよう それでは、作成したテンプレートを表示させていきたいと思います。 app.pyからhtmlファイルを呼び出して表示させるのですが、ここでrender_templateの出番です。 render_template('index.html')と指定することで、pythonファイルからindex.htmlを呼び出すことができます。 実際にコードを書いてみましょう。 FlaskScraping/app.py ... # ルーティングの設定 @app.route('/') def index(): (-)return '<h1>Hello World!</h1>' (+)return render_template('index.html') # このコードによって、Flaskがweb上で走ります if __name__ == '__main__': # debug=Trueとすることで、保存するだけで変更が適用されます # エラーも見やすくなります app.run(debug=True) では、実行してみましょう！ (flaskenv)FlaskScraping$ python app.py これで、HelloFlask!と表示されれば成功です！ 今回実装したコード 以上で、テンプレートの作成と表示ができるようになりました。 今回実装したコードの全体像を表示しておきます。確認してください！ FlaskScraping/app.py # モジュールのインポート from flask import Flask, render_template import requests from bs4 import BeautifulSoup # Flaskアプリのインスタンスを作成 app = Flask(__name__) # ルーティングの設定 @app.route('/') def index(): return render_template('index.html') # このコードによって、Flaskがweb上で走ります if __name__ == '__main__': # debug=Trueとすることで、保存するだけで変更が適用されます # エラーも見やすくなります app.run(debug=True) templates/base.html <!DOCTYPE html> <html lang="ja"> <head> <meta charset="UTF-8"> <link href="https://cdn.jsdelivr.net/npm/bootstrap@5.0.0-beta1/dist/css/bootstrap.min.css" rel="stylesheet" integrity="sha384-giJF6kkoqNQ00vy+HMDP7azOuL0xtbfIcaT9wjKHr8RbDVddVHyTfAAsrekwKmP1" crossorigin="anonymous"> <title>Flask+スクレイピングで記事一覧画面を作ろう</title> </head> <body> {% block content %} {% endblock %} </body> </html> templates/index.html {% extends 'base.html' %} {% block content %} <p>Hello Flask!</p> {% endblock %} 次回は、スクレイピングしたデータを表示できるようにします。 ありがとうございました?‍♂️

本記事では「レッドコーダーが教える、競プロ・AtCoder上達のガイドライン【初級編：競プロを始めよう】」で紹介されているAOJの「Introduction To Programming I」の40問をPythonで解説します。 ITP1_1_A Hello World ITP1_1_A print("Hello world") 【解説】 print関数を使用し、Hello Worldを出力します。 詳しい解説はこちら ITP1_1_B X:Cubic ITP1_1_B x = int(input()) print(x ** 3) 【解説】 input()はキーボードからの入力を文字列型として受け取るのでint()で数値型にしてから処理します。 詳しい解説はこちら ITP1_1_C Rectangle ITP1_1_C a, b = map(int, input().split()) print(a*b, 2*(a+b)) 【解説】 input().split()を使用し空白区切りで複数の入力を行います。ただし文字列として受け取るのでmap関数の第一引数にint関数を設定し、数値にしてから変数a,bに設定します。入力については「AtCoderで始めるPython入門」が参考になります。 詳しい解説はこちら ITP1_1_D Watch ITP1_1_D s = int(input()) h = s // 3600 m = s % 3600 // 60 s = s % 60 print("{0}:{1}:{2}".format(h, m, s)) 【解説】秒数を3600で割った商がhに、秒数を3600で割った余りを60で割った商がmに、秒数を60で割った余りがsになります。h,m,sの表示にはformat関数を使用します。中括弧のに挟まれた番号に従ってformat関数の引数が設定されます。 詳しい解説はこちら ITP1_2_A Small, Large, or Equal ITP1_2_A a, b = map(int, input().split()) if a > b: print('a > b') elif a < b: print('a < b') elif a == b: print('a == b') else: pass 【解説】if文を使用し、条件ごとに表示する文字を設定します。 詳しい解説はこちら ITP1_2_B Range ITP1_2_B a, b, c = map(int, input().split()) if a < b < c: print('Yes') else: print('No') 【解説】a　<　b　and　b　<　cと記述しても問題ありません。 詳しい解説はこちら ITP1_2_C Sorting Three Numbers ITP1_2_C get_lst = list(map(int, input().split())) get_lst.sort() print(*get_lst) 【解説】listは複数の要素を含むオブジェクトです。map関数を使用し、数値型（int型）として受け取った複数の値をlistとして変数に設定します。リスト名.sort()とすることで昇順に並べ替え、print関数での表示の際には、リスト名の前にアスタリスクをつけてリストの中身を表示します。 詳しい解説はこちら ITP1_2_D Circle in a Rectangle ITP1_2_D W, H, x, y, r = map(int, input().split()) if x - r < 0: print('No') elif y - r < 0: print('No') elif x + r > W: print('No') elif y + r > H: print('No') else: print('Yes') 【解説】横と縦でそれぞれ円が長方形からはみ出さないかを確認します。 下記の場合は長方形から円がはみ出るので、”No”を表示します。 条件式がx – r < 0の場合、円が左にはみ出ます。 条件式がy – r < 0の場合、円が下にはみ出ます。 条件式がx + r > Wの場合、円が右にはみ出ます。 条件式がy + r > Hの場合、円が上にはみ出ます。 上記以外の場合、円は長方形の中にあるので、”Yes”を表示します。 詳しい解説はこちら ITP1_3_A Print Many Hello World ITP1_3_A for _ in range(1000): print('Hello World') 【解説】for文を使用し、繰り返し処理を行います。今回は1000回処理を繰り返すのでfor文のイテラブルにrange(1000)を設定します。そうすると0〜999までの1000回処理が行われます。 詳しい解説はこちら ITP1_3_B Print Test Cases ITP1_3_B num_lst = [] idx = [] for i in range(10000): num_lst.append(int(input())) idx.append(i+1) if num_lst[i] == 0: num_lst.pop() break for j in range(len(num_lst)): print('Case {}: {}'.format(idx[j], num_lst[j])) 【解説】入力用のfor文と出力用のfor文で分けています。数字入力用の配列num_lstとインデックス用の配列idxを用意し、0が入力されるまで入力用のfor文を回します。0が入力されたら、popメソッドよ使用し、最後に入力した数字を削除します。出力用のfor文で入力された数字を順番通りに出力します。 詳しい解説はこちら ITP1_3_C Swapping Two Numbers ITP1_3_C while True: x, y = map(int, input().split()) if x == 0 and y == 0: break elif x < y: print(x, y) else: print(y, x) 【解説】whileで0 0が入力されるまで実行します。あとはx,yの大小関係を確認して出力します。 詳しい解説はこちら ITP1_3_D How Many Divisors? ITP1_3_D a, b, c = map(int, input().split()) cnt = 0 for i in range(a, b + 1): if c % i == 0: cnt += 1 print(cnt) 【解説】aからbまでfor文を回し、cの約数かどうか調べます。約数であればcnt（カウント）を1つ大きくします。 詳しい解説はこちら ITP1_4_A A/B Problem ITP1_4_A a, b = map(int, input().split()) print('{0} {1} {2:.5f}'.format(a//b, a % b, a/b)) 【解説】//は商、%は余りを求める演算子です。表示にはformat関数を使用し、{0} {1} {2}の番号通りに引数が表示されます。 詳しい解説はこちら ITP1_4_B Circle ITP1_4_B import math r = float(input()) print('{0:.6f} {1:.6f}'.format(r*r*math.pi, 2*r*math.pi)) ITP1_4_B（別解） from math import pi r = float(input()) print('{0:.6f} {1:.6f}'.format(r*r*pi, 2*r*pi)) 【解説】input関数をfloat関数で囲み、実数として値を受け取ります。mathモジュールをインポートする事で円周率を使用する事ができるのでそれを使用し、円の面積、円周の長さを表示します。 詳しい解説はこちら ITP1_4_C Simple Calculator ITP1_4_C while True: a, op, b = map(str, input().split()) a = int(a) b = int(b) if op == "?": break elif op == "+": print(int(a + b)) elif op == "-": print(int(a - b)) elif op == "*": print(int(a * b)) elif op == "/": print(int(a / b)) else: pass ITP1_4_C（別解） while True: table = input().split() a = int(table[0]) op = table[1] b = int(table[2]) if op == "?": break elif op == "+": print(int(a + b)) elif op == "-": print(int(a - b)) elif op == "*": print(int(a * b)) elif op == "/": print(int(a / b)) else: pass 【解説】whileで?がでるまで実行します。あとは足し算、引き算、掛け算、割り算の演算子の入力に従って計算結果を出力します。 詳しい解説はこちら ITP1_4_D Min, Max and Sum ITP1_4_D n = int(input()) get_lst = list(map(int, input().split())) print(min(get_lst), max(get_lst), sum(get_lst)) 【解説】Pythonではlistの最大や最小、合計をそれぞれmax(list名)min(list名)sum(list名)で求めることができます。 詳しい解説はこちら ITP1_5_A Print a Rectangle ITP1_5_A while True: H, W = map(int, input().split()) if H == 0 and W == 0: break for i in range(H): print("#"*W) print() 【解説】whileループを使用し、H, Wがともに0が入力されるまで入力を続けます。　そして、forループでW(横の長さ)分"#"を表示します。 詳しい解説はこちら ITP1_5_B Print a Frame ITP1_5_B while True: H, W = map(int, input().split()) if H == 0 and W == 0: break for i in range(H): for k in range(W): if i == 0 or i == H-1 or k == 0 or k == W-1: print("#", end="") else: print(".", end="") print() print() 【解説】以下のときは描画する四角形の枠と判断し"#"を描画します。 range(H)で0~H-1までの値を取り出す際の0は上の枠 range(H)で0~H-1までの値を取り出す際のH-1は下の枠 range(W)で0~W-1までの値を取り出す際の0は左の枠 range(W)で0~W-1までの値を取り出す際のW-1は右の枠 それ以外は"."を描画します。 詳しい解説はこちら ITP1_5_C Print a Chessboard ITP1_5_C while True: H, W = map(int, input().split()) if H == 0 and W == 0: break for i in range(H): for j in range(W): if (i+j) % 2 == 0: print("#", end="") else: print(".", end="") print() print() 【解説】forのなかにforを書くことができます。二重ループと言われます。（参考:「Pythonで多重ループ（ネストしたforループ）からbreak」)出力したいチェック柄をじっと見ると、行数と列数を足した値が偶数の時は#、奇数の時は.であることに気づきます（Pythonでは0行目、1行目…と数えることに注意）。各行で順に出力していき、最後の列で改行します。ただprint("#")のようにすると#のあとに改行されてしまいます。それを避けるためにend=""で出力時に改行しないようにします。 詳しい解説はこちら ITP1_5_D Structured Programming ITP1_5_D n = int(input()) for i in range(1, n + 1): if i % 3 == 0 or "3" in str(i): print(" {}".format(i), end="") print() 【解説】問題の意味をまずは理解します。C++のコードを読めないと辛い問題ですが、競技プログラミングの解説ではC++がよく使われるため、読めるようにしておくといいかもしれません。例えば通称「蟻本」と呼ばれる本もコードはC++で書かれています。今回の問題では「n以下の自然数のうち、3の倍数もしくは3がつく数を小さいものから順に出力しなさい」という問題です。3がつくことの確認は今回は文字列としてinを用いて行いました。任意の文字列を含むかの判定は「Pythonで文字列を検索（〜を含むか判定、位置取得、カウント）」が参考になります。 詳しい解説はこちら ITP1_6_A Reversing Numbers ITP1_6_A n = int(input()) get_lst = list(map(int, input().split())) get_lst.reverse() print(*get_lst) 【解説】input().split()で空白区切りの入力をint型に変換したmap関数をlist関数の引数に設定することにより、リスト型として変数get_lstに設定します。リストの順を逆順にするには、リスト型のメソッドreverse、組み込み関数のreversed、スライスによって並び替えることができます。今回はリスト型のメソッドreverseを使用しました。リストの逆順は「Pythonでリストや文字列を逆順に並べ替え」が、スライスは「Pythonのスライスによるリストや文字列の部分選択・代入」が参考になります。 詳しい解説はこちら ITP1_6_B Finding Missing Cards ITP1_6_B all_cards = [(s, n) for s in ['S', 'H', 'C', 'D'] for n in range(1, 14)] n = int(input()) hold_cards = [] for i in range(n): suit, num = input().split() num = int(num) hold_cards.append((suit, num)) for card in all_cards: if card not in hold_cards: print(*card) 【解説】all_cards = [(s, n) for s in ['S', 'H', 'C', 'D'] for n in range(1, 14)]でトランプの全てのパターンを作成し、一つ目のfor文でn組のカードの取得、二つ目のfor文で不足しているカードの出力を行います。二つ目のfor文ではnot inを使用し、全ての組み合わせからn組のカード以外のカードの出力を行います。 詳しい解説はこちら ITP1_6_C Official House ITP1_6_C houses = [[[0 for r in range(10)] for f in range(3)] for b in range(4)] n = int(input()) for i in range(n): b, f, r, v = map(int, input().split()) houses[b - 1][f - 1][r - 1] += v for b in range(4): for f in houses[b]: print('', *f) if b != 3: print('#' * 20) 【解説】この問題は多次元配列で解くことができます。「Pythonで多次元配列を扱う方法【初心者向け】」が参考になります。出力の形式を合わせるのが少し大変な問題です。 houses = [[[0 for r in range(10)] for f in range(3)] for b in range(4)]で10部屋、3階建ての公舎4棟のグループの作成を行います。一つ目のfor文で4つの整数b(棟), f(階), r(番目), v(人)の読み込み、二つ目のfor文で出力を行います。 詳しい解説はこちら ITP1_6_D Matrix Vector Multiplication ITP1_6_D n, m = map(int, input().split()) A = [list(map(int, input().split())) for _ in range(n)] b = [int(input()) for _ in range(m)] for i in range(n): sum = 0 for j in range(m): sum += A[i][j] * b[j] print(sum) 【解説】この問題も多次元配列で処理することで解くことができます。 まず、n×mの行列の読み込み変数Aに設定。続いてm×1のベクトルを読み込み変数bに設定。 for文で入れ子にし、行列の計算を行います。 詳しい解説はこちら ITP1_7_A Grading ITP1_7_A while True: m, f, r = map(int, input().split()) if m == f == r == -1: break if m == -1 or f == -1 or m + f < 30: print('F') elif m + f >= 80: print('A') elif m + f >= 65: print('B') elif m + f >= 50: print('C') elif m + f >= 30: if r >= 50: print('C') else: print('D') else: pass 【解説】while Trueで無限ループにし、整数値が全て-1が入力されるまで入力を続けます。続いて条件分岐ですが、まず成績がFランクの生徒かどうか判断、続いてelifで点数の大きい順にAからDまでの条件分岐を行います。プログラムは上から順に実行されるため、点数順にif文を記述すると条件式の簡略化につながり効率が良いです。 詳しい解説はこちら ITP1_7_B How many ways? ITP1_7_C while True: n, x = map(int, input().split()) if n == 0 and x == 0: break cnt = 0 for i in range(1, n-1): for j in range(i+1, n): if j < x - i - j <= n: cnt += 1 print(cnt) ITP1_7_C（別解） while True: n, x = map(int, input().split()) if n == 0 and x == 0: break cnt = 0 for i in range(1, n + 1): for j in range(i+1, n + 1): for k in range(j+1, n + 1): if i + j + k == x: cnt += 1 print(cnt) 【解説】while文とif文を使用し、空白区切りで0が入力されるまで入力を続けます。 はじめに外側のfor文で1からn-1までを選びiの設定します、内側のfor文ではそれより大きいi+1からnを選びjに設定します。あとはx-i-jがjよりも大きくn以下であれば前に選んだ数よりも大きくなり、組み合わせとして成立します。 別解では、for文を3つ使用し、range関数の開始値をずらして設定し、合計値がxと同じ時に組み合わせの一つとしてカウントします。 詳しい解説はこちら ITP1_7_C Spreadsheet ITP1_7_C r, c = map(int, input().split()) sheet = [list(map(int, input().split())) for i in range(r)] for i in range(r): sheet[i].append(sum(sheet[i])) Column_sum = [0]*(c+1) for j in range(c+1): for k in range(r): Column_sum[j] += sheet[k][j] for i in range(r): print(*sheet[i]) print(*Column_sum) 【解説】はじめにr*cの2次元配列の入力を行います。続いて、appendメソッドとsum関数を使用し、行の合計値を追加します。続いてのfor文で列の合計値をもとめ、最後の表示用のfor文で行の合計値を追加した行列の表示、続いて、列の合計値の表示を行います。 詳しい解説はこちら ITP1_7_D Matrix Multiplication ITP1_7_D n, m, l = map(int, input().split()) A = [list(map(int, input().split())) for _ in range(n)] B = [list(map(int, input().split())) for _ in range(m)] for i in range(n): C = [] for j in range(l): tmp = 0 for k in range(m): tmp += A[i][k] * B[k][j] C.append(tmp) print(*C) 【解説】ITP1_6_Dと同様、行列の計算です。行列の計算の性質上、一方の行数と他方の列数（またはその逆）が同じでなければなりません。なので一番内側のfor文で一方の列と他方の行のカウントアップを行い、合計値を算出します。 詳しい解説はこちら ITP1_8_A Toggling Cases ITP1_8_A words = input() print(words.swapcase()) 【解説】input関数を使い、キーボードから文字列を受け取り変数に代入した後、swapcaseメソッドを使用し、大文字と小文字を変換します。 詳しい解説はこちら ITP1_8_B Sum of Numbers ITP1_8_B while True: str_x = input() if str_x == "0": break ans = 0 for n in str_x: ans += int(n) print(ans) 【解説】ある入力までその処理を続ける場合はwhile Trueで無限ループを回し、その中で入力を行い、if文で終了条件を設定、break文でループ処理を終了します。今回は123が入力されると、1+2+3という処理をしたいので、input関数をint関数で囲まず、文字列として入力を受け取ります。これは後で使用するfor文のイテラブル（繰り返し可能なオブジェクト）として文字列を設定すると文字列が一つずつ取り出され、処理を行う事ができるからです。今回はその合計値を出力します。 詳しい解説はこちら ITP1_8_C Counting Characters ITP1_8_C import sys texts=sys.stdin.read() texts=texts.lower() cnt=[0]*26 letters=list('abcdefghijklmnopqrstuvwxyz') for x in texts: i=0 for y in letters: if x==y: cnt[i]+=1 i+=1 for i in range(26): print(letters[i]+" : "+str(cnt[i])) 【解説】今回は複数行かつ行数がわからない英文を読み取らなければならないのでsysモジュールのsys.stdin.readを利用します。複数行の文字列を複数行のまま受け取ることができます。入力後は、lowerメソッドですべて小文字にして、あらかじめ用意したaからzまでの文字列と比較、一致した場合に配列のカウンタを+1します。 詳しい解説はこちら ITP1_8_D Ring ITP1_8_D s = input() p = input() s = s * 2 if s.find(p) != -1: print('Yes') else: print('No') 【解説】リング状の文字列からある文字列を作成できるかを考えます。今回の問題ではpはsよりも短いという条件がついています。そのためリング2周分で考えれば十分です。リング2周に該当する文字列を作成し、そこでpを作ることができるか調べます。pがあるかを確かめるために文字列のメソッドfindを利用します。findは特定の文字列の位置を取得しますが、文字列がなければ-1を返します。 詳しい解説はこちら ITP1_9_A Finding a Word ITP1_9_A import sys word = input() text = sys.stdin.read() print(text.lower().split().count(word)) 【解説】まずinput関数で単語の入力を行います。続いて文章の入力はsys.stdin.readで行います。text.lower().split().count(word)について解説します。lower()メソッドで小文字にします。split()メソッドでスペースごとに区切りリストにします。リストのなかに含まれる個数を調べるにはcountメソッドを利用します。 詳しい解説はこちら ITP1_9_B Shuffle ITP1_9_B while True: cards = input() if cards == "-": break m = int(input()) for i in range(m): sh = int(input()) former = cards[:sh] later = cards[sh:] cards = later+former print(cards) 【解説】リストのスライスを利用して解くことができます。指定の文字位置より前と後に分けて変数に設定し、交換したものを表示します。 詳しい解説はこちら ITP1_9_C Card Game ITP1_9_C n = int(input()) T = 0 H = 0 for i in range(n): card_t, card_h = input().split() if card_t == card_h: T += 1 H += 1 else: if card_h > card_t: H += 3 else: T += 3 print(T, H) 【解説】文字列に>や<を用いると、辞書順で比較した結果が返されます。これを利用して解くことができます。 詳しい解説はこちら ITP1_9_D Transformation ITP1_9_D text = input() q = int(input()) for i in range(q): order = input().split() a, b = map(int, order[1:3]) if order[0] == "print": print(text[a:b+1]) elif order[0] == "reverse": re_text = text[a:b+1] text = text[:a]+re_text[::-1]+text[b+1:] else: text = text[:a]+order[3]+text[b+1:] 【解説】先ほどと同じく、リストのスライスを利用して解くことができます。ただしreverseのときに少し注意が必要です。スライスでstepを-1にするときstartとendは文字通りの解釈がされません。こちらの質問が参考になります。今回はreverseしたい範囲のリストを作成してから逆順にしました。 詳しい解説はこちら ITP1_10_A Distance ITP1_10_A import math x1, y1, x2, y2 = map(float, input().split()) print(math.hypot(x2 - x1, y2 - y1)) 【解説】mathモジュールの1つであるmath.hypotを使用することで原点から点 (x, y)のベクトルの長さを求める事ができます。 詳しい解説はこちら ITP1_10_B Triangle ITP1_10_B import math a, b, C = map(int, input().split()) rad = math.radians(C) S = 0.5*a*b*math.sin(rad) h = 2*S/a c = math.sqrt(a**2+b**2-2*a*b*math.cos(rad)) L = a+b+c print(S, L, h, sep="\n") 【解説】三角形の面積、周の長さ、高さの求め方はたくさんありますので、mathモジュールをうまく利用して解きます。 詳しい解説はこちら ITP1_10_C Standard Deviation ITP1_10_C while True: n = int(input()) if n == 0: break score = list(map(int, input().split())) mean = sum(score)/n var_sum = 0 for i in range(n): var_sum += (score[i]-mean)**2 print((var_sum/n)**0.5) 【解説】while Trueとif文を使用し、終了条件の0が入力された際に処理が終了するようにします。続いて、得点をリストとして読み込み、その平均値を求めます。続いて、for文で各得点から平均値を引いたものを2乗し、その合計値を求めます。この合計値（偏差の合計）をデータの総数で割ったもの(分散)を0.5乗して、標準偏差を算出します。 詳しい解説はこちら ITP1_10_D Distance II ITP1_10_D n = int(input()) X = list(map(int, input().split())) Y = list(map(int, input().split())) # p=1,2,3 for p in range(1, 4): print("{0:.6f}".format(sum([abs(a-b)**p for a, b in zip(X, Y)])**(1/p))) # p=infinity print("{0:.6f}".format(max([abs(a-b) for a, b in zip(X, Y)]))) 【解説】p=1,2,3のミンコフスキー距離とp=infinityのチェビシェフの距離は別に計算をして表示します。 詳しい解説はこちら

ABC218のA-E問題の解説。 目次 A - Weather Forecast B - qwerty C - Shapes D - Rectangles E - Destruction A - Weather Forecast 解説 $N$日後の天気予報が晴れ、つまり、oかどうかを判定する問題。 $S$の$N-1$文字目がoだったらYes。 コード def main(): n = int(input()) s = input() if s[n-1] == 'o': print('Yes') else: print('No') if __name__ == '__main__': main() B - qwerty 解説 Pの入力に対して、対応する英小文字をつなげていくという問題。 数字を文字列に変換するときは、chr()を用いる。chr(97)がaにあたるので、$+96$して文字を連結してあげれば良い。 コード def main(): plist = list(map(int, input().split())) s = '' for p in plist: s += chr(p+96) print(s) if __name__ == '__main__': main() C - Shapes 解説 $S$のグリッドの中で、$T$の#でつくられた形があるかどうかを判定する問題。 まずは、SとTの#の数が等しいかをみる。そもそも#の数が等しくなければ、この2つを一致させることができない。 一致していれば、ループに入る。 記録したsの#の座標をもとに、tの座標をひとつずつみていく。 一致しなかった場合は、ただちにループを抜けることで計算量を削減できる。 また、一致しなかった場合は、sの座標自体を90度回転させ、同じように#の位置を見ていく。 コード def main(): n = int(input()) s = list(list(input()) for _ in range(n)) t = list(list(input()) for _ in range(n)) # sの'#'の座標を記録 slist = list((i, j) for i in range(n) for j in range(n) if s[i][j] == '#') # tの'#'の数を記録 tcnt = sum(1 for i in range(n) for j in range(n) if t[i][j] == '#') # '#'の数が合うかを確認し、あっていなければ'No' if len(slist) != tcnt: print("No") exit() # 90度x3回分回す for _ in range(4): for i in range(-n+1, n): for j in range(-n+1, n): flag = True for x, y in slist: if not (0 <= x+i < n and 0 <= y+j < n and t[x+i][y+j] == '#'): flag = False break if flag: print("Yes") exit() # 90度回転 slist = [(y, n-1-x) for x, y in slist] print("No") if __name__ == '__main__': main() D - Rectangles 解説 $2$次元平面上にある$N$個の点のうちから4つを選んで、何個の長方形が作れるかという問題。 まずxを固定して、それに対応するyを配列に保存する。つまり、各xがもつyを配列として保存する。 それらxを2つずつみていき、そのx_i, x_jがもつyの個数が2以上の場合は、x_i, x_jの共通するyを見つける。それらから2つをとった組み合わせの個数を数えてあげると、x_i, x_jとの長方形の数がわかるので、これを足していったものが答え。 コード def main(): from itertools import combinations n = int(input()) d = dict() for _ in range(n): x, y = map(int, input().split()) if x not in d: d[x] = [y] else: d[x].append(y) # xの座標 xaxis = [k for k in d.keys()] cnt = 0 for i in range(len(xaxis)-1): for j in range(i+1, len(xaxis)): # そのx_i, x_jがもつyの個数が2以上の場合 if len(d[xaxis[i]]) >= 2 and len(d[xaxis[j]]) >= 2: # x_i, x_jの共通するyを見つける test = set(d[xaxis[i]]) & set(d[xaxis[j]]) # それらの組み合わせの個数を数える cnt += len(list(combinations(test, 2))) print(cnt) if __name__ == '__main__': main() E - Destruction 解説 N頂点とM辺の連結無向グラフが存在するとき、重みのついた辺を取り除くと、その値に応じて報酬がもらえる。また、重みが負の値であれば、その値の絶対値分の罰金を支払う。 このとき、$0$本以上の辺を取り除き、すべての頂点がつながっている場合に得られる、最大の報酬を求める問題。 先に結論をいってしまうと、これは最小全域木を求めるアルゴリズムを用いると解くことができる。 このアルゴリズムには、プリム法とクラスカル法がある。こちらの仕組みついては解説記事をご覧いただきたい。どちらのアルゴリズムでも最小全域木を求めることにはかわりはない。 そもそも最小全域木というのは、存在するすべての頂点がつながっていて、辺の重みの総和が最小になるデータ構造のことをいう。 この問題で、なぜ最小全域木を求める必要があるのか。それは、罰金があるすべての辺、または、報酬が小さい辺はできるだけ残しておいた全域木を作成する必要があるためである。 ここからはコードを使って説明していく。 まず、重み（報酬）が正の値を取るとき、costに値を足し、先にすべての報酬を受けとっておく。これが辺をすべて取り除いた場合の報酬の最大値である。しかし、これではグラフが連結であるという条件を満たしていない。つまり、ここから最小全域木を作成する必要がある。 最小全域木をつくるにあたって、最小全域木の辺の値が負（罰金）の場合、辺を取り除く必要がない、つまり、罰金は取る必要がないため、そのままにしておく。 最小全域木の辺の値が正の場合、先程の操作でcostが報酬の最大値なので、報酬を受け取れない分、costから報酬分を引く。 これをすべての辺に対して行うことで、今回の問題の報酬の最大値を求めることができる。 コード1（プリム法） プリム法ではheapq（優先度付きキュー）を用いているが、より詳しく知りたい方はこちらを参考にしていただきたい。 def main(): from heapq import heappush, heappop n, m = map(int, input().split()) graph = [[] for _ in range(n)] cost = 0 for _ in range(m): u, v, c = map(int, input().split()) u, v = u-1, v-1 graph[u].append((v, c)) # u->vの辺 graph[v].append((u, c)) # v->uの辺 if c >= 0: cost += c # プリム法 # 頂点がマークされているか確認する配列 marked = [False for _ in range(n)] # マークされている頂点数を数える marked_cnt = 0 # はじめに0頂点をマーク marked[0] = True marked_cnt += 1 # heap q = [] # 頂点0に隣接する辺を保存 for j, c in graph[0]: heappush(q, (c, j)) # すべての頂点をマークするまでwhileループ while marked_cnt < n: # 最小の重みの辺をheapで取り出す c, i = heappop(q) # マークされているならスキップ if marked[i]: continue # 頂点をマーク marked[i] = True marked_cnt += 1 # 罰金でなければ報酬を # 受け取りすぎているので報酬分の値を引く if c >= 0: cost -= c # 頂点iに隣接する辺を保存 for j, c in graph[i]: # マークされていればスキップ if marked[j]: continue heappush(q, (c, j)) print(cost) if __name__ == '__main__': main() コード2（クラスカル法） クラスカル法ではUnion-Findを用いているが、より詳しく知りたい方はこちらを参考にしていただきたい。 # Union-Find from collections import defaultdict class UnionFind(): def __init__(self, n): self.n = n self.parents = [-1] * n def find(self, x): if self.parents[x] < 0: return x else: self.parents[x] = self.find(self.parents[x]) return self.parents[x] def union(self, x, y): x = self.find(x) y = self.find(y) if x == y: return if self.parents[x] > self.parents[y]: x, y = y, x self.parents[x] += self.parents[y] self.parents[y] = x def size(self, x): return -self.parents[self.find(x)] def same(self, x, y): return self.find(x) == self.find(y) def members(self, x): root = self.find(x) return [i for i in range(self.n) if self.find(i) == root] def roots(self): return [i for i, x in enumerate(self.parents) if x < 0] def group_count(self): return len(self.roots()) def all_group_members(self): group_members = defaultdict(list) for member in range(self.n): group_members[self.find(member)].append(member) return group_members def __str__(self): return '\n'.join(f'{r}: {m}' for r, m in self.all_group_members().items()) # Code def main(): n, m = map(int, input().split()) uf = UnionFind(n) # 辺を保存する配列 edges = [] # 報酬、つまり、辺の重みを保存する変数 cost = 0 for _ in range(m): u, v, c = map(int, input().split()) u, v = u-1, v-1 # インデックスずらし edges.append((c, u, v)) # c_iが0以上のとき、c_iの報酬を得る if c >= 0: cost += c # 重み（報酬）が小さい順に辺をソート edges.sort() # できるだけ辺（報酬）をとりたいので、 # 負（罰金）であるか値が小さい辺（報酬）を残した最小全域木を作りたい # （それ以外の辺は報酬として受け取れる） for edge in edges: c, u, v = edge # 罰金なら、それらの頂点をつなげる if c <= 0: uf.union(u, v) else: # 頂点がつながっていなければ if not uf.same(u, v): cost -= c # 報酬を引き uf.union(u, v) # 頂点同士をつなげる print(cost) if __name__ == '__main__': main() 編集後記 最小全域木結構面白いけど、例題が少ないのが難点。

MathⅠ RSAで暗号化された文を復号する問題 RSA暗号 https://ja.wikipedia.org/wiki/RSA暗号 桁数が大きい合成数の素因数分解問題が困難であることを安全性の根拠とした公開鍵暗号の一つ。 暗号化 平分を暗号化するには以下の式を計算する必要がある。 c = m^e\; mod\; n cは暗号化文、mは平分、eは適当に決める値、nは２つの大きな素数p,qの積であり、これらは自分で設定する。問題では既にこれらの値はわかっている。 暗号化する時には、復号に必要な秘密鍵dを一緒に作成しているが、詳しいことは後述する。 復号化 暗号文を復号するには以下の式を計算する必要がある。 m = c^d\; mod\; n 今回の問題は秘密鍵dを求めることが肝になっている。 解く cを復号して平分mを求める問題になっている。 問題をみると、e,n,c,p,qは既に決まっている。 復号化で述べたように復号するには秘密鍵dが必要になる。dは以下の計算で求められる。 d = e^{-1}\; mod\; lcm((p-1),(q-1)) lcm((p-1),(q-1))は、p-1とq-1の最小公倍数、a mod bはaをbで割ったあまりのこと。 これでdが求められるので実際にpythonで計算してみる。 import math e=65537 n=1517330236262917595314610888889322115651087080826711948897066340883208205571592392362650858571076247939805436226544833224526137582834770402681005343930059463684528957271778199162575053306238099823295117697031968370690372250916935800738698142103275969223264184374648246277564306900886005299731265812255274723175925185522344831066577166867786835955092059346244885587228196357297758371381557924676260190209536670230008561217008649261974735505203813478978893582292682827884118215872470401293272325715864815977064075988643101088355047954735427424641386870772845440782632933485165110172437511822736907550777817722248753671107339823410418938404382732079381329288400012929311347390423061254658780185245562668131009832293474920208834795460061115101364091252176594144096675899952570380792978037217747311595899301451192342027799533264325948876556110474850761538179748318187805312451895898751337975457949549497666542175077894987697085521882531938339334715190663665300179658557458036053188152532948734992896239950564081581184284728802682982779186068791931259198917308153082917381616147108543673346682338045309449569430550618884202465809290850964525390539782080230737593560891353558335337408957948041667929154230334506735825418239563481028126435029 c=225549592628492616152632265482125315868911125659971085929712296366214355608049224179339757637982541542745010822022226409126123627804953064072055667012172681551500780763483172914389813057444669314726404135978565446282309019729994976815925850916487257699707478206132474710963752590399332920672607440793116387051071191919835316845827838287954541558777355864714782464299278036910958484272003656702623646042688124964364376687297742060363382322519436200343894901785951095760714894439233966409337996138592489997024933882003852590408577812535049335652212448474376457015077047529818315877549614859586475504070051201054704954654093482056493092930700787890579346065916834434739980791402216175555075896066616519150164831990626727591876115821219941268309678240872298029611746575376322733311657394502859852213595389607239431585120943268774679785316133478171225719729917877009624611286702010936951705160870997184123775488592130586606070277173392647225589257616518666852404878425355285270687131724258281902727717116041282358028398978152480549468694659695121115046850718180640407034795656480263573773381753855724693739080045739160297875306923958599742379878734638341856117533253251168244471273520476474579680250862738227337561115160603373096699944163 q=44481453884385518268018625442920628989497457642625668259648790876723318635861137128631112417617317160816537010595885992856520476731882382742220627466006460645416066646852266992087386855491152795237153901319521506429873434336969666536995399866125781057768075533560120399184566956433129854995464893265403724034960689938351450709950699740508459206785093693277541785285699733873530541918483842122691276322286810422297015782658645129421043160749040846216892671031156465364652681036828461619272427318758098538927727392459501761203842363017121432657534770898181975532066012149902177196510416802134121754859407938165610800223 p=34111525225922333955113751419357677129436029651245533697825114748126342624744832960936498161825269430327019858323450578875242014583535842110912370431931233957939950911741013017595977471949767235426490850284286661592357779825212265055931705799916913817655743434497422993498931394618832741336247426815710164342599150990608143637331068220244525541794855651643135012846039439355101027994945120698530177329829213208761057392236875366458197098507252851244132455996468628957560178868724310000317011912994632328371761486669358065577269198065792981537378448324923622959249447066754504943097391628716371245206444816309511381323 n=p*q lcm = ((p-1)*(q-1) // math.gcd(p-1,q-1)) print(lcm) d=pow(e,-1) % lcm print(d) int too large to convert to floatというエラーが出てきた。 値が大きすぎてオーバーフローしてしまった。 他のやり方でdを計算してみる。 F=lcm((p-1),(q-1)) とする。dはFを法としたeの逆数であることから、以下の式が成り立つ。 de \equiv 1\; mod\; F この計算は拡張されたユークリッド互除法で容易に求められる。 ちなみに上の式の意味は、deをFで割ったあまり＝1をFで割ったあまりであることを意味している。当然1をFで割ったあまりは1であるため、deをFで割ったあまりは1であり、0 = de-1 mod Fとも書ける。 拡張されたユークリッド互除法を使ってpythonで解いてみる。 参考：https://tech.camph.net/rsa-public-key-encryption/ import math import binascii def ex_euclid(x, y): c0, c1 = x, y a0, a1 = 1, 0 b0, b1 = 0, 1 while c1 != 0: m = c0 % c1 q = c0 // c1 c0, c1 = c1, m a0, a1 = a1, (a0 - q * a1) b0, b1 = b1, (b0 - q * b1) return a0 e=65537 n=1517330236262917595314610888889322115651087080826711948897066340883208205571592392362650858571076247939805436226544833224526137582834770402681005343930059463684528957271778199162575053306238099823295117697031968370690372250916935800738698142103275969223264184374648246277564306900886005299731265812255274723175925185522344831066577166867786835955092059346244885587228196357297758371381557924676260190209536670230008561217008649261974735505203813478978893582292682827884118215872470401293272325715864815977064075988643101088355047954735427424641386870772845440782632933485165110172437511822736907550777817722248753671107339823410418938404382732079381329288400012929311347390423061254658780185245562668131009832293474920208834795460061115101364091252176594144096675899952570380792978037217747311595899301451192342027799533264325948876556110474850761538179748318187805312451895898751337975457949549497666542175077894987697085521882531938339334715190663665300179658557458036053188152532948734992896239950564081581184284728802682982779186068791931259198917308153082917381616147108543673346682338045309449569430550618884202465809290850964525390539782080230737593560891353558335337408957948041667929154230334506735825418239563481028126435029 c=225549592628492616152632265482125315868911125659971085929712296366214355608049224179339757637982541542745010822022226409126123627804953064072055667012172681551500780763483172914389813057444669314726404135978565446282309019729994976815925850916487257699707478206132474710963752590399332920672607440793116387051071191919835316845827838287954541558777355864714782464299278036910958484272003656702623646042688124964364376687297742060363382322519436200343894901785951095760714894439233966409337996138592489997024933882003852590408577812535049335652212448474376457015077047529818315877549614859586475504070051201054704954654093482056493092930700787890579346065916834434739980791402216175555075896066616519150164831990626727591876115821219941268309678240872298029611746575376322733311657394502859852213595389607239431585120943268774679785316133478171225719729917877009624611286702010936951705160870997184123775488592130586606070277173392647225589257616518666852404878425355285270687131724258281902727717116041282358028398978152480549468694659695121115046850718180640407034795656480263573773381753855724693739080045739160297875306923958599742379878734638341856117533253251168244471273520476474579680250862738227337561115160603373096699944163 q=44481453884385518268018625442920628989497457642625668259648790876723318635861137128631112417617317160816537010595885992856520476731882382742220627466006460645416066646852266992087386855491152795237153901319521506429873434336969666536995399866125781057768075533560120399184566956433129854995464893265403724034960689938351450709950699740508459206785093693277541785285699733873530541918483842122691276322286810422297015782658645129421043160749040846216892671031156465364652681036828461619272427318758098538927727392459501761203842363017121432657534770898181975532066012149902177196510416802134121754859407938165610800223 p=34111525225922333955113751419357677129436029651245533697825114748126342624744832960936498161825269430327019858323450578875242014583535842110912370431931233957939950911741013017595977471949767235426490850284286661592357779825212265055931705799916913817655743434497422993498931394618832741336247426815710164342599150990608143637331068220244525541794855651643135012846039439355101027994945120698530177329829213208761057392236875366458197098507252851244132455996468628957560178868724310000317011912994632328371761486669358065577269198065792981537378448324923622959249447066754504943097391628716371245206444816309511381323 n=p*q F = ((p-1) * (q-1)) // math.gcd(p-1, q-1) d=ex_euclid(e,F) m = pow(c,d,n) print("%0512x"%m) binascii.unhexlify('546869732069732034303936626974205253412e0a54686520666c616720697320464c41475f4e625a366751654b444d56464a4d61550a50414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050414444494e472050') FLAGが得られた。 print ("%0512x"%m).decode("hex")でうまくデコードできなかったためbinasciiを使ってデコードした。ここからでもできる。 まとめ RSAの復号ができた。 p、qが知られると秘密鍵が計算できてしまい第3者に復号されてしまう危険性があることがわかった。

python,Qiita初心者なので、備忘録として記載していきます。 なにせ初心者なので、知識不足はご理解ください。 知識を深めながら追記していきたいと思います。 ブール型とは 基本 bool型はint型の派生クラス Trueと1は等価 Falseと0は等価 値 意味 True 真 False 偽 >>> True == 1 True >>> False == 0 True Pythonにおいて、偽(False)とみなされるのは以下のオブジェクト 定数：NoneとFalse 数値型におけるゼロ：「0」, 「0.0,」 「0j」(複素数), Decimal(0), Fraction(0, 1) 空のシーケンスまたはコレクション：'', (), [], {}, set(), range(0) 上記以外はすべて真(True)と判定される bool関数：ブール型に変換する bool関数で真偽値の確認 >>> bool(10) True >>> bool(0) False >>> bool('0') True >>> bool('False') True >>> bool('') False >>> bool([False]) True >>> bool([]) False >>> bool(1) True >>> bool(0) False int型の派生クラスなので、普通に演算することも可能 >>> issubclass(bool, int) # intの派生クラスという事が分かる True >>> +True 1 >>> -True -1 >>> -False 0 >>> True + True # 1 * 1 として演算 2 >>> 10 * True # 10 * 1 として演算 10 >>> 10 / False # 10 / 0 として演算 # 0除算エラーになる >>> 'はい' * True # 'はい' * 1 として演算 'はい' >>> 'はい' * False # 'はい' * 0として演算 '' >>> ('はい','いいえ')[True] # []はブール型になる式でも良い 'いいえ' >>> ('はい','いいえ')[False] # []はブール型になる式でも良い 'はい'

注意事項 本記事は、2021年9月15日時点におけるDatabricksでの調査結果です。 概要 Delta Lake 1.0から利用できるようになったGenerated Column（生成列）利用時におけるMerge操作に制約があるのですが、Merge内のInsert処理時にエラーとなるようです。 Delta Lakeのドキュメントにて、Generated Columnには下記のようにMerge処理がサポートされていないことが記載されております。 MERGE operations do not support generated columns. 引用元：Table batch reads and writes — Delta Lake Documentation Merge処理時のどの操作をサポートしていないのを検証したところ、Update処理をサポートしておりましたが、Insert処理がサポートされていないことがわかりました。下記が、Insert処理時のエラーメッセージです。チェック制約違反となることが原因であることから、Generated Columns列の制約を満たすカラムをソースとすることでMerege操作時にはエラーとなりませんでした。 Error in SQL statement: InvariantViolationException: CHECK constraint Generated Column (count_string <=> > CAST(count AS STRING)) violated by row with values: - count : 1 - count_string : null Databricks Runtime 9.1以降では、Insert処理にも対応するでき、実質的にMerge操作に制約がなくなるようです。 下記が実行結果です。 検証環境 databricks runtime : 9.0.x-scala2.12 Python version : 3.8.10 pyspark version : 3.1.2 検証手順 下記のGithub pagesのページをご確認ください。 Databricks Runtime 9.1での実行結果は下記のページです。 コードを実行したい方は、下記のdbcファイルを取り込んでください。 https://github.com/manabian-/databricks_tecks_for_qiita/blob/main/tecks/merge_test_generate_columns/merge_test_generate_columns.dbc

Lexとは AWSのbotのAIサービス。Amazon Lex を使うと、すべてのデベロッパーが Amazon Alexa に採用されている深層学習技術と同じ技術を利用できるため、自然言語での高度な対話ボット (チャットボット) を短時間で簡単に構築できます。 今回はlexを用いてアイスクリームのフレーバーや容器について答えるボットを作成して行きます。 コードの量があるので順番に解説して行きます。 ボット作成ファイルの全体像 lex_create_bot.py import boto3 import time iam = boto3.client('iam') iam.create_service_linked_role(AWSServiceName='lex.amazonaws.com') lex = boto3.client('lex-models', 'us-east-1') #フレーバーのスロットタイプの作成 flavor_slot_type = lex.put_slot_type( name='FlavorSlotType', enumerationValues=[ {'value': 'vanilla'}, {'value': 'chocolate', 'synonyms': ['choc']}, {'value': 'strawberry', 'synonyms': ['berry']} ], valueSelectionStrategy='TOP_RESOLUTION', createVersion=True) print('slot type:', flavor_slot_type['name']) #容器のスロットタイプを作成 container_slot_type = lex.put_slot_type( name='ContainerSlotType', enumerationValues=[ {'value': 'corn'}, {'value': 'cup'} ], valueSelectionStrategy='TOP_RESOLUTION', createVersion=True) print('slot type:', container_slot_type['name']) #インテントの作成 intent = lex.put_intent( name='OrderIntent', #インテント内のスロット slots=[ { 'name': 'Flavor', 'slotConstraint':'Required', 'slotType':'FlavorSlotType', 'slotTypeVersion': '1', 'valueElicitationPrompt': { 'messages': [{ 'contentType': 'PlainText', 'content': 'Vanilla, chocolate or strawberry?' }], 'maxAttempts': 3 } }, { 'name': 'Container', 'slotConstraint': 'Required', 'slotType': 'ContainerSlotType', 'slotTypeVersion': '1', 'valueElicitationPrompt': { 'messages': [{ 'contentType': 'PlainText', 'content': 'Corn or cup?' }], 'maxAttempts': 3 } } ], #発話例 sampleUtterances=[ 'I want {Flavor} ice cream in {Container}', '{Flavor} ice cream {Container}', 'ice create' ], #完了時のセリフ conclusionStatement={ 'messages': [{ 'contentType': 'PlainText', 'content': 'OK, {Flavor} ice cream in {Container}' }], }, #完了の動作 fulfillmentActivity={'type': 'ReturnIntent'}, createVersion=True) print('intent:',intent['name']) #ボットの作成 bot = lex.put_bot( name ='MyBot', locale='en-US', childDirected=False, #インテント intents=[ { 'intentName': 'OrderIntent', 'intentVersion': '1' } ], #中止時のセリフ abortStatement={ 'messages':[ { 'contentType': 'PlainText', 'content': 'Please try again.' } ] }, voiceId='Joanna', createVersion=True) print('bot:', bot['name']) #ボット作成の進捗表示 start = time.time() status = '' while status not in ['READY', 'FAILED']: #ボットを取得 status = lex.get_bot(name='MyBot', versionOrAlias='1')['status'] time.sleep(10) print('{:7.2f} {}'.format(time.time()-start, status)) #作成に失敗した場合は理由を表示 if status == 'FAILED': print(lex.get_bot( name='Mybot', versionOrAlias='1')['failureReason']) #ボットエイリアスを作成 bot_alias = lex.put_bot_alias( name='MyBotAlias', botName='MyBot', botVersion='1') print('bot alias:', bot_alias['name']) ①importとサービスクライアント作成 lex_create_bot.py import boto3 import time iam = boto3.client('iam') iam.create_service_linked_role(AWSServiceName='lex.amazonaws.com') lex = boto3.client('lex-models', 'us-east-1') boto3とtimeのインポートとiamのロールとlexのサービスクライアントを作成 ②スロットタイプ作成 スロットとはユーザーから聞き出す必要のあるパラメーターです。 今回はフレーバーと容器をスロットに設定しています。 lex_create_bot.py #フレーバーのスロットタイプ flavor_slot_type = lex.put_slot_type( name='FlavorSlotType', enumerationValues=[ {'value': 'vanilla'}, {'value': 'chocolate', 'synonyms': ['choc']}, {'value': 'strawberry', 'synonyms': ['berry']} ], valueSelectionStrategy='TOP_RESOLUTION', createVersion=True) print('slot type:', flavor_slot_type['name']) #容器のスロットタイプ container_slot_type = lex.put_slot_type( name='ContainerSlotType', enumerationValues=[ {'value': 'corn'}, {'value': 'cup'} ], valueSelectionStrategy='TOP_RESOLUTION', createVersion=True) print('slot type:', container_slot_type['name']) put_slot_typeメソッドでスロットを作成します。 enumerationValuesでスロットが取りうる値を辞書型で指定 valueSelectionStrategyの'TOP_RESOLUTION'はスロットタイプの中からユーザーの言葉に近いものを返します。今回でいえばchocoと入力すると本来の値のchocolateが返ってくる。 ③インテントの作成 intentとは意思や意図という意味がありますが、lexの中ではユーザーと会話する台本のイメージになります。 今回をアイスクリーム注文するコマンドを受け取るインテントを作成します。パラメーターとしてアイスクリームのフレーバーと容器を受け取リます。 put_intentメソッドでインテントを作成します。 引数slotsについて lex_create_bot.py #インテントの作成 intent = lex.put_intent( name='OrderIntent', #インテント内のスロット slots=[ { 'name': 'Flavor', 'slotConstraint':'Required', 'slotType':'FlavorSlotType', 'slotTypeVersion': '1', 'valueElicitationPrompt': { 'messages': [{ 'contentType': 'PlainText', 'content': 'Vanilla, chocolate or strawberry?' }], 'maxAttempts': 3 } }, { 'name': 'Container', 'slotConstraint': 'Required', 'slotType': 'ContainerSlotType', 'slotTypeVersion': '1', 'valueElicitationPrompt': { 'messages': [{ 'contentType': 'PlainText', 'content': 'Corn or cup?' }], 'maxAttempts': 3 } } ], 引数slotsはインテント内にスロットを作成するために使います。 ・slotConstraint:スロットの制約です。'Required'または'Optinal'で選択。 ・slotTypeVersion:スロットのバージョン。最初はバージョン'1' ・valueElicitationPrompt:スロットの値をユーザーから促すためのセリフ →valueElicitationPromptに指定する辞書の項目 ・message:セリフを表す辞書リスト。その中にcontentTypeとcontentが含まれる ・maxAttempts:聞き返す回数の上限 →→messageに指定する辞書の項目 ・contentType:セリフのタイプ。'PlainText','SSML','CustomPayload'から指定 ・content:セリフの内容 引数sampleUtterances lex_create_bot.py #発話例 sampleUtterances=[ 'I want {Flavor} ice cream in {Container}', '{Flavor} ice cream {Container}', 'ice create' ], 引数sampleUtterancesはインテントに対するユーザーの発話例です。 {スロット名}でスロットを表します。{Flavor}の部分にはフレーバー(chocolate,strawberry,vanilla)が入り、{Container}には(cupなど)が入ります。 今回ユーザーの発話を見ることはなさそうですね。 引数conclusionStatement lex_create_bot.py #完了時のセリフ conclusionStatement={ 'messages': [{ 'contentType': 'PlainText', 'content': 'OK, {Flavor} ice cream in {Container}' }], }, 引数conclusionStatementはボットがユーザーから必要なスロットの情報を聞き出すことができて、インテントが完了したときのセリフです。辞書型の項目として引数slotsのvalueElicitationPromptに指定する辞書の項目と同じです。 引数fulfillmentActivity lex_create_bot.py #完了の動作 fulfillmentActivity={'type': 'ReturnIntent'}, createVersion=True) print('intent:',intent['name']) 引数fulfillmentActivityはインテントが完了したときの動作です。 引数にはキー'type'を含む辞書を指定します。 typeの種類 ReruenIntent：呼び出し元のプログラムにインテントの情報を返します。 CodeHook：Lambda呼び出してLambdaに登録されたプログラムを実行します。 ④botの作成 最後にボットとボットエイリアスを作成します。Lexでbotを公開するにはボットエイリアスを作成する必要があります。 lex_create_bot.py #ボットの作成 bot = lex.put_bot( name ='MyBot', locale='en-US', childDirected=False, #インテント intents=[ { 'intentName': 'OrderIntent', 'intentVersion': '1' } ], #中止時のセリフ abortStatement={ 'messages':[ { 'contentType': 'PlainText', 'content': 'Please try again.' } ] }, voiceId='Joanna', createVersion=True) print('bot:', bot['name']) put_botメソッドを呼び出してボットを作成します。 引数のnameは名前、localeは今回'en-US'、childDirected(子供向け)は子供向けの場合にはTrue、子供向けではない場合はFalseに指定。 引数intentはnameとintentVersionを辞書型で渡します。 引数abortStatementには中止時のセリフを辞書型のリストにして指定します。 voiceIdは'en-US'の中から'Joanna'に指定します。 作成時の進捗時間表示 lex_create_bot.py #ボット作成の進捗表示 start = time.time() status = '' while status not in ['READY', 'FAILED']: #ボットを取得 status = lex.get_bot(name='MyBot', versionOrAlias='1')['status'] time.sleep(10) print('{:7.2f} {}'.format(time.time()-start, status)) #作成に失敗した場合は理由を表示 if status == 'FAILED': print(lex.get_bot( name='Mybot', versionOrAlias='1')['failureReason']) #ボットエイリアスを作成 bot_alias = lex.put_bot_alias( name='MyBotAlias', botName='MyBot', botVersion='1') print('bot alias:', bot_alias['name']) get_botメソッドでボットの情報を取得しstatusへ、10秒間時間を待ち開始からに時間と取得したボットの情報をstatusに代入して表示します。 作成に失敗した場合は、失敗した理由を取得し表示。 put_bot_aliasメソッドでボットエイリアスを作成し、画面に表示します。 まとめ ここまで長くなりましたが、lexでのbot作成のコードと解説でした。 次回はbotにテキストを与えて動かしてみます。 参考文献 この記事は以下の情報を参考にして執筆しました AWSでつくるAIプログラミング入門

Python って型書けるの？ 書けます！！！ Pythonは別にコンパイルするわけではないので、単体で実行前に静的解析ができるわけではないです...... が、サードパーティ制のソフトウェアを活用することでいい感じ(※要検証)にすることができます。 (あと、実は実行時にも書いた型は影響しません) typing --- 型ヒントのサポート Pythonの公式ドキュメントを読んでください。 完 ......いやまあこれが真理ではあるんですが、簡単に説明しつつ個人的なハマりポイントを紹介していきます。 typing のつけかた hoge: str = "Hello, World!" 変数の場合、変数名の後に : <型名> と記述することで型を書くことができます。 def func(hoge: str) -> str: return "Hello, World!" 関数の戻り値は、後ろに -> <型名> をつけることで指定できます。 使える型 List や Callable など、 typing.py に記載された型 https://github.com/python/cpython/blob/3.9/Lib/typing.py#L1631-L1709 int や str などの組み込みクラス 自作のクラス 型変数 など (ちなみに上2つは下2つに含まれるので、分ける必要があるかといえばない) ジェネリクス ドキュメントを読め 意訳：別の回で話でもしようかな Ellipsis ググラビリティがカスな悪い組み込みオブジェクトに ... というものがあります。 t: Tuple[str, ...] = ("H", "e", "l", "l", "o", "w") 型を省略するときに使う記号で、以下のように実装に(内部的に)利用されています。 if len(params) == 2 and params[1] is ...: msg = "Tuple[t, ...]: t must be a type." p = _type_check(params[0], msg) return self.copy_with((p, _TypingEllipsis)) ハマりどころ 循環インポート 返り値を自分自身にする (メソッドチェーン) class A: def func(self) -> A: # エラー return self なんかこんな感じのやつを書きたくなることってあるじゃないですか？ fluent interfaceみたいなことしたくなった時とか (こっちはこっちで長くなりそうだから別で作ろう) ただ、クラス定義内で自身の名前を使うことはできないので、 NameError: name '<type>' is not defined みたいなエラーになります。(もちろんIDEにも怒られます) こんな場合には forward-reference という機能を使います。 class A: def func(self) -> 'A': return self 文字列として型名を記述することで、解決を後回しにすることができます。 (IDEでは依存関係が解決できた場合、実際に型を記述したときと同じ挙動をします) 相互的に参照を持つ これはそもそも設計が悪いから見直せ みたいな記事が散見されたので、良くないんだろうなあとは思いつつ、必要なものは必要なので使わせてくれ！みたいな人向け (代替案的な設計があればぜひコメントください) main.py from a import A from b import B b: B = B() a: A = A(b) b.set_a(a) a.py from b import B class A: def __init__(self, b: B): self.b: B = b b.py from typing import Optional from a import A class B: def __init__(self): self.a: Optional[A] = None def set_a(self, a: A): self.a = a このように、A, Bの両方で相手を参照するようなインポートを行おうとすると以下のようなエラーが出ます ImportError: cannot import name 'A' from partially initialized module 'a' (most likely due to a circular import) (a.py) これを解決するためには、 typing.TYPE_CHECKING を利用します https://docs.python.org/ja/3/library/typing.html#typing.TYPE_CHECKING a.py from typing import TYPE_CHECKING if TYPE_CHECKING: from b import B class A: def __init__(self, b: 'B'): self.b: B = b b.py from typing import Optional, TYPE_CHECKING if TYPE_CHECKING: from a import A class B: def __init__(self): self.a: Optional[A] = None def set_a(self, a: 'A'): self.a = a TYPE_CHECKING は型検査中のみに True になる (※非保証) 値。 if文などでインポートをチェックすることにより、エラーを回避できる。 ただ、この状態だと Type Hinting に使っている型名が存在しないのでエラーとなる。 NameError: name 'B' is not defined ため、前述した forward-refarence と組み合わせることにより評価も遅延して解決することができる。 複雑な(再帰的)データを型付けする { "type": "or", "filters": [{ "type": "or", "filters": [...] }, { "type": "range", "from": 1, "to": 10 }] } なんかこんな感じで再帰的にデータが入りうるような型を作りたい！ → 型変数を使いましょう TypeOp = Dict[str, Union[str, int, List['TypeOp']]] json_: TypeOp = { "type": "or", "filters": [{ "type": "or", "filters": [...] }, { "type": "range", "from": 1, "to": 10 }] } 戻り値がvoid Python に void なんてものはない。 def func() pass print(func()) # None Python では return を省略すると None が返るので、型も同様に指定すると良い 上限つきワイルドカード型 親クラスを直接書いてやれば動きます Type.TypeVar を使うことで無理矢理 Java のときみたいにも書ける (書いてた) んですが、必要はなさそうです。 キャストしたい isinstance でチェックすればよしなにしてくれます。 ただ、直接触っても予測変換がでないので一旦正しい型の変数に代入しましょう。 a: Any = 123 if isinstance(a, A): b: int = a 参考文献 PythonとType Hintsで書くバックエンド PyCharm の型ヒント Pythonではじまる、型のある世界 今から始める型安全 Python

Python って型書けるの？ 書けます！！！ Pythonは別にコンパイルするわけではないので、単体で実行前に静的解析ができるわけではないです...... が、サードパーティ制のソフトウェアを活用することでいい感じ(※要検証)にすることができます。 (あと、実は実行時にも書いた型は影響しません) typing --- 型ヒントのサポート Pythonの公式ドキュメントを読んでください。 完 ......いやまあこれが真理ではあるんですが、簡単に説明しつつ個人的なハマりポイントを紹介していきます。 typing のつけかた hoge: str = "Hello, World!" 変数の場合、変数名の後に : <型名> と記述することで型を書くことができます。 def func(hoge: str) -> str: return "Hello, World!" 関数の戻り値は、後ろに -> <型名> をつけることで指定できます。 使える型 List や Callable など、 typing.py に記載された型 https://github.com/python/cpython/blob/3.9/Lib/typing.py#L1631-L1709 int や str などの組み込みクラス 自作のクラス 型変数 など (ちなみに上2つは下2つに含まれるので、分ける必要があるかといえばない) ジェネリクス ドキュメントを読め 意訳：別の回で話でもしようかな Ellipsis ググラビリティがカスな悪い組み込みオブジェクトに ... というものがあります。 t: Tuple[str, ...] = ("H", "e", "l", "l", "o", "w") 型を省略するときに使う記号で、以下のように実装に(内部的に)利用されています。 if len(params) == 2 and params[1] is ...: msg = "Tuple[t, ...]: t must be a type." p = _type_check(params[0], msg) return self.copy_with((p, _TypingEllipsis)) ハマりどころ 循環インポート 返り値を自分自身にする (メソッドチェーン) class A: def func(self) -> A: # エラー return self なんかこんな感じのやつを書きたくなることってあるじゃないですか？ fluent interfaceみたいなことしたくなった時とか (こっちはこっちで長くなりそうだから別で作ろう) ただ、クラス定義内で自身の名前を使うことはできないので、 NameError: name '<type>' is not defined みたいなエラーになります。(もちろんIDEにも怒られます) こんな場合には forward-reference という機能を使います。 class A: def func(self) -> 'A': return self 文字列として型名を記述することで、解決を後回しにすることができます。 (IDEでは依存関係が解決できた場合、実際に型を記述したときと同じ挙動をします) 相互的に参照を持つ これはそもそも設計が悪いから見直せ みたいな記事が散見されたので、良くないんだろうなあとは思いつつ、必要なものは必要なので使わせてくれ！みたいな人向け (代替案的な設計があればぜひコメントください) main.py from a import A from b import B b: B = B() a: A = A(b) b.set_a(a) a.py from b import B class A: def __init__(self, b: B): self.b: B = b b.py from typing import Optional from a import A class B: def __init__(self): self.a: Optional[A] = None def set_a(self, a: A): self.a = a このように、A, Bの両方で相手を参照するようなインポートを行おうとすると以下のようなエラーが出ます ImportError: cannot import name 'A' from partially initialized module 'a' (most likely due to a circular import) (a.py) これを解決するためには、 typing.TYPE_CHECKING を利用します https://docs.python.org/ja/3/library/typing.html#typing.TYPE_CHECKING a.py from typing import TYPE_CHECKING if TYPE_CHECKING: from b import B class A: def __init__(self, b: 'B'): self.b: B = b b.py from typing import Optional, TYPE_CHECKING if TYPE_CHECKING: from a import A class B: def __init__(self): self.a: Optional[A] = None def set_a(self, a: 'A'): self.a = a TYPE_CHECKING は型検査中のみに True になる (※非保証) 値。 if文などでインポートをチェックすることにより、エラーを回避できる。 ただ、この状態だと Type Hinting に使っている型名が存在しないのでエラーとなる。 NameError: name 'B' is not defined ため、前述した forward-refarence と組み合わせることにより評価も遅延して解決することができる。 複雑な(再帰的)データを型付けする { "type": "or", "filters": [{ "type": "or", "filters": [...] }, { "type": "range", "from": 1, "to": 10 }] } なんかこんな感じで再帰的にデータが入りうるような型を作りたい！ → 型変数を使いましょう TypeOp = Dict[str, Union[str, int, List['TypeOp']]] json_: TypeOp = { "type": "or", "filters": [{ "type": "or", "filters": [...] }, { "type": "range", "from": 1, "to": 10 }] } 戻り値がvoid Python に void なんてものはない。 def func() pass print(func()) # None Python では return を省略すると None が返るので、型も同様に指定すると良い 上限つきワイルドカード型 親クラスを直接書いてやれば動きます Type.TypeVar を使うことで無理矢理 Java のときみたいにも書ける (書いてた) んですが、必要はなさそうです。 キャストしたい isinstance でチェックすればよしなにしてくれます。 ただ、直接触っても予測変換がでないので一旦正しい型の変数に代入しましょう。 a: Any = 123 if isinstance(a, A): b: int = a 参考文献 PythonとType Hintsで書くバックエンド PyCharm の型ヒント Pythonではじまる、型のある世界 今から始める型安全 Python

概要 公式のドキュメント「クイックスタート: コマンド ラインから Azure に Python 関数を作成する 」を整理・補足したものになります。 最後に、matplotlibを使ってグラフを生成してレスポンスするウェブAPIのサンプルコードを付けています。 動作確認環境 Windows Windows 10 Pro 20H2 準備 Azure Functions Core Tools のインストール Azure Functions Core Tools の最新版をインストールする。コマンドラインからバージョンを確認する。3.0.3734(x64) AzureFunctionsCoreToolsのバージョン確認 func --version Azure CLI のインストール と Azureログイン Azure CLI の最新版をインストールする。コマンドラインからバージョンを確認する。2.28.0 AzureCLIのバージョン確認 az --version Azure CLI で Azure にログインする。コマンドラインから下記を実行すると、ウェブブラウザが立ち上がるので指示に従って Azure アカウントにログインする。 Azurログイン az login Pythonのバージョン確認と仮想環境の設定 Pythonのバージョンを確認する。Function App（関数アプリ）で、Python（ランタイムスタック）は 3.6 から 3.9 に対応している。 ここでは、3.8.5 を想定します。 Pythonバージョン確認 python -V Pythonの仮想環境 azure-func を作成して有効化する（virtualenv + virtualenvwrapper を想定）。pipのバージョンを最新にする。 mkvirtualenv azure-func workon azure-func pip install --upgrade pip プロジェクトの作成 関数アプリのプロジェクトを管理するためのフォルダ AzureFunc02 を作成して移動する。 cd "C:\Users\xxxx\Documents\Projects\AzureFunc02" プロジェクトの初期化とHTTPトリガの関数 HttpTrigger01（テンプレートとなるものを自動生成）する。 func new --name HttpTrigger01 --template "HTTP trigger" --authlevel "anonymous" 実行すると Use the up/down arrow keys to select a worker runtime: と尋ねられるので上下キーで Python を選択する。 生成された関数フォルダ HttpTrigger01 に移動する。__init__.py と function.json の2ファイルがある。 cd HttpTrigger01 __init__.py を確認して、必要なモジュールをインストールする。また、編集を行なう。 pip install azure.functions ローカルでの動作確認 プロジェクト AzureFunc02 のトップフォルダに戻って、ローカルで動作テストをする。 ローカルで動作テスト func start --python --port 8000 クラウドに発行して動作確認 Azure に xxxx-p-func2 という名前の Function App（関数アプリ）を新規作成する。既にAzureに、XXXX-JapanEast-RG というリソースグループ、storage01xxxx というストレージアカウントが作成済みであることを前提にしています。 Azureにリソース追加 az functionapp create --consumption-plan-location japaneast ^ --runtime python ^ --runtime-version 3.8 ^ --functions-version 3 ^ --name xxxx-p-func2 ^ --os-type linux ^ --resource-group XXXX-JapanEast-RG ^ --storage-account storage01xxxx ^ --disable-app-insights true ローカルで開発してものを Azure に発行する。 発行する func azure functionapp publish xxxx-p-func2 補足 requirements.txt に matplotlib を追加する。 __init__.py（Matplotlibで画像を生成してレスポンス） import logging import io import numpy as np import azure.functions as func import matplotlib import matplotlib.pyplot as plt matplotlib.use('Agg') def main(req: func.HttpRequest) -> func.HttpResponse: logging.info('Python HTTP trigger function processed a request.') image = io.BytesIO() x = np.linspace(0, 10) y = np.sin(x) plt.plot(x, y) plt.savefig(image, format='png') image.seek(0) return func.HttpResponse(image.read(),mimetype='image/png') 参考資料 VS Code の拡張機能で同様のことを行なう。

始めに 知り合いがSTM32L4A6を使ったボードを自作し、これにMicroPythonを入れて動かして欲しいと頼まれたので、動かすまでにやったことの備忘録です。といっても、ほとんどこちらを参考にさせていただきました。 MicroPythonビルド環境構築 こちらを参考 自作ボードにMicroPythonを対応させる 今回、自作ボード用にMicroPythonに追加したフォルダ、ファイルは以下の通り。 ports/stm32/ ├ boards/ │ └ STM32L4a6GDISC/ -> new folder │ ├ mpconfigboard.h -> new file │ ├ mpconfigboard.mk -> new file │ ├ pins.csv -> new file │ └ stm32l4xx_hal_conf.h -> new file ├ stm32l4a6_af.csv 　　　　　　-> new file └ stm32l4a6xg.ld 　　　　　　　-> new file 自作ボードの回路図・ピンアサ、チップのデータシート、すでにMicroPythonでサポートされている自分と似たチップのデータを参考にこれらのファイルを作成、編集します。 stm32l4a6xg.ld (リンカスクリプト) チップのデータシートのMemory mappingを参考に記述します。 今回は、同チップシリーズのSTM32L496のMemory mappingと見比べたところ、同じメモリ配置であったため、すでにMicroPythonプロジェクトに存在するstm32l496xg.ldをそのままコピーしました。 stm32l4a6xg.ld /* GNU linker script for STM32L4A6XG */ /* Specify the memory areas */ MEMORY { FLASH (rx) : ORIGIN = 0x08000000, LENGTH = 612K /* sectors 0-305 */ FLASH_FS (r) : ORIGIN = 0x08099000, LENGTH = 412K /* sectors 306-511 412 KiB */ RAM (xrw) : ORIGIN = 0x20000000, LENGTH = 320K /* SRAM1, 256K + SRAM2, 64K */ } /* produce a link error if there is not this amount of RAM for these sections */ _minimum_stack_size = 2K; _minimum_heap_size = 16K; /* Define the stack. The stack is full descending so begins just above last byte of RAM, or bottom of FS cache.. Note that EABI requires the stack to be 8-byte aligned for a call. */ /* RAM extents for the garbage collector */ _ram_start = ORIGIN(RAM); _ram_end = ORIGIN(RAM) + LENGTH(RAM); _ram_fs_cache_end = _ram_end; _ram_fs_cache_start = _ram_fs_cache_end - 2K; /* fs cache = 2K */ _estack = _ram_fs_cache_start - _estack_reserve; _sstack = _estack - 16K; /* stack = 16K */ _heap_start = _ebss; /* heap starts just after statically allocated memory */ _heap_end = _sstack; /* bss + heap = 302K, tunable by adjusting stack size */ _flash_fs_start = ORIGIN(FLASH_FS); _flash_fs_end = ORIGIN(FLASH_FS) + LENGTH(FLASH_FS); stm32l4a6_af.csv (Alternate Functionsリスト) チップのデータシートのAlternate functionテーブルとpin definitionsテーブルを参考に記述します。 今回は、同チップシリーズのSTM32L496のデータシートと見比べたところ、テーブルの内容が同一のものであったため、既に存在するstm32l496_af.csvをそのままコピーしました。 pins.csv 自作ボードの回路図を見ながら、各ピンの名前と信号線の対応を記述します。 後ほど、MicroPython上からは以下のように使うことができます。 from machine import Pin cs = Pin("SPI2_NSS", Pin.OUT) cs.on() mpconfigboard.h (ボードのコンフィグ) 自作ボードの回路図、他のボードのmpconfigboard.hを参考にしながら記述します。 今回はこのようになりました。 mpconfigboard.h #define MICROPY_HW_BOARD_NAME "HOGE_BOARD" #define MICROPY_HW_MCU_NAME "STM32L4A6" #define MICROPY_HW_ENABLE_INTERNAL_FLASH_STORAGE (0) #define MICROPY_HW_ENABLE_RNG (1) #define MICROPY_HW_ENABLE_RTC (1) #define MICROPY_HW_ENABLE_TIMER (1) #define MICROPY_HW_ENABLE_SDCARD (1) // MSI is used and is 4MHz, // Resulting core frequency is 80MHz: #define MICROPY_HW_CLK_PLLM (1) #define MICROPY_HW_CLK_PLLN (40) #define MICROPY_HW_CLK_PLLP (RCC_PLLP_DIV7) #define MICROPY_HW_CLK_PLLR (RCC_PLLR_DIV2) #define MICROPY_HW_CLK_PLLQ (RCC_PLLQ_DIV2) #define MICROPY_HW_FLASH_LATENCY FLASH_LATENCY_4 // The board has an external 32kHz crystal #define MICROPY_HW_RTC_USE_LSE (1) // USART config #define MICROPY_HW_UART1_TX (pin_A9) #define MICROPY_HW_UART1_RX (pin_A10) #define MICROPY_HW_UART2_TX (pin_A2) #define MICROPY_HW_UART2_RX (pin_A3) // USART 1 is connected to the virtual com port on the ST-LINK #define MICROPY_HW_UART_REPL PYB_UART_1 #define MICROPY_HW_UART_REPL_BAUD 115200 // I2C buses #define MICROPY_HW_I2C1_SCL (pin_B6) #define MICROPY_HW_I2C1_SDA (pin_B7) #define MICROPY_HW_I2C2_SCL (pin_B10) #define MICROPY_HW_I2C2_SDA (pin_B11) // SPI buses #define MICROPY_HW_SPI1_NSS (pin_A4) #define MICROPY_HW_SPI1_SCK (pin_A5) #define MICROPY_HW_SPI1_MISO (pin_A6) #define MICROPY_HW_SPI1_MOSI (pin_A7) #define MICROPY_HW_SPI2_NSS (pin_B12) #define MICROPY_HW_SPI2_SCK (pin_B13) #define MICROPY_HW_SPI2_MISO (pin_B14) #define MICROPY_HW_SPI2_MOSI (pin_B15) // LED (The orange LED is controlled over MFX) #define MICROPY_HW_LED1 (pin_B0) #define MICROPY_HW_LED2 (pin_B1) #define MICROPY_HW_LED_ON(pin) (mp_hal_pin_low(pin)) #define MICROPY_HW_LED_OFF(pin) (mp_hal_pin_high(pin)) // SD card #define MICROPY_HW_SDCARD_SDMMC (1) #define MICROPY_HW_SDCARD_CK (pin_C12) #define MICROPY_HW_SDCARD_CMD (pin_D2) #define MICROPY_HW_SDCARD_D0 (pin_C8) #define MICROPY_HW_SDCARD_D1 (pin_C9) #define MICROPY_HW_SDCARD_D2 (pin_C10) #define MICROPY_HW_SDCARD_D3 (pin_C11) #define MICROPY_HW_SDCARD_DETECT_PIN (pin_A8) #define MICROPY_HW_SDCARD_DETECT_PULL (GPIO_PULLUP) #define MICROPY_HW_SDCARD_DETECT_PRESENT (GPIO_PIN_RESET) mpconfigboard.mk mpconfigboard.mk MCU_SERIES = l4 CMSIS_MCU = STM32L4A6xx AF_FILE = boards/stm32l4a6_af.csv LD_FILES = boards/stm32l4a6xg.ld boards/common_basic.ld OPENOCD_CONFIG = boards/openocd_stm32l4.cfg stm32l4xx_hal_conf.h こちらもSTM32L496用のものをそのまま流用しました。 stm32l4xx_hal_conf.h /* This file is part of the MicroPython project, http://micropython.org/ * The MIT License (MIT) * Copyright (c) 2019 Damien P. George */ #ifndef MICROPY_INCLUDED_STM32L4XX_HAL_CONF_H #define MICROPY_INCLUDED_STM32L4XX_HAL_CONF_H #include "boards/stm32l4xx_hal_conf_base.h" // Oscillator values in Hz #define HSE_VALUE (8000000) #define LSE_VALUE (32768) #define EXTERNAL_SAI1_CLOCK_VALUE (48000) #define EXTERNAL_SAI2_CLOCK_VALUE (48000) // Oscillator timeouts in ms #define HSE_STARTUP_TIMEOUT (100) #define LSE_STARTUP_TIMEOUT (5000) #endif // MICROPY_INCLUDED_STM32L4XX_HAL_CONF_H その他　ports/stm32/adc.c 上記ファイルを追加し、ビルドを実行したところ、adc.cでエラーが発生したので、必要箇所に、 defined(STM32L4A6xx) を追記しました。 ビルド ports/stm32ディレクトリに移動し、 $ make BOARD=STM32L4a6GDISC 作成されたbuild-STM32L4a6GDISC/firmware.hexをチップに書き込むと、無事MicroPythonが起動しました。 おまけ SDカードへのアクセス SD内のファイルを確認 >>> import os >>> os.listdir("/sd") ['appl.mot', 'kernel_config.txt', 'System Volume Information', 'hello.py', 'tmp', 'spi.py', 'spi_bmx160.py'] >>> SD内のPythonスクリプトを実行 >>> exec(open("/sd/hello.py").read()) Hello World >>> SPI通信でセンサの情報を取得する 自作ボードにはボッシュ社のBMX160 9軸センサが載っており、チップとはSPIでつながっています。MicroPythonからセンサの情報を取得するスクリプトは以下のようになります。 spi_bmx160.py from struct import unpack import time from machine import SPI, Pin WRITE_BIT = 0 << 7 READ_BIT = 1 << 7 CHIP_ID_ADDRESS = 0x00 MAG_CONF_ADDRESS = 0x44 MAG_IF_0_ADDRESS = 0x4C MAG_IF_1_ADDRESS = 0x4D MAG_IF_2_ADDRESS = 0x4E MAG_IF_3_ADDRESS = 0x4F CMD_REG_ADDRESS = 0x7E DATA_START_ADDRESS = 0x04 DATA_END_ADDRESS = 0x1A DATA_NUM = DATA_END_ADDRESS - DATA_START_ADDRESS + 1 spi = SPI(2) cs = Pin("SPI2_NSS", Pin.OUT) cs.on() time.sleep(0.5) def write_spi(reg_address: int, write_data: int) -> None: send_data = WRITE_BIT | reg_address buf = bytearray([send_data, write_data]) print(f"Write: {buf}") cs.off() spi.write(buf) cs.on() time.sleep(0.1) def read_spi_1byte(reg_address: int) -> None: send_data = READ_BIT | reg_address print(f"Send: 0x{send_data:x}") cs.off() recv = spi.read(2, send_data) cs.on() print(f"Recv: 0x{recv[1]:x}\n") def read_spi_multiple(start_address: int, nbytes: int) -> bytearray: write_data = READ_BIT | start_address cs.off() read_buf = spi.read(nbytes + 1, write_data) cs.on() return read_buf[1:] def convert_mag_data(mag_data_byte: bytearray) -> float: mag_data_raw = unpack("<h", mag_data_byte)[0] return mag_data_raw / 16.0 def convert_gyro_data(gyro_data_byte: bytearray) -> float: gyro_data_raw = unpack("<h", gyro_data_byte)[0] return gyro_data_raw * 61.0 / 1000 def convert_acc_data(acc_data_byte: bytearray) -> float: acc_data_raw = unpack("<h", acc_data_byte)[0] return acc_data_raw * 2 / 0x7FFF def convert_sensor_time(sensor_time_byte: bytearray) -> float: sensor_time_raw = int.from_bytes(sensor_time_byte, "little") return sensor_time_raw * 39 / 1000000 print("CHIP ID") read_spi_1byte(CHIP_ID_ADDRESS) write_spi(CMD_REG_ADDRESS, 0x11) write_spi(CMD_REG_ADDRESS, 0x15) write_spi(CMD_REG_ADDRESS, 0x19) write_spi(MAG_IF_0_ADDRESS, 0x80) write_spi(MAG_IF_3_ADDRESS, 0x01) write_spi(MAG_IF_2_ADDRESS, 0x4B) write_spi(MAG_IF_3_ADDRESS, 0x01) write_spi(MAG_IF_2_ADDRESS, 0x51) write_spi(MAG_IF_3_ADDRESS, 0x0E) write_spi(MAG_IF_2_ADDRESS, 0x52) write_spi(MAG_IF_3_ADDRESS, 0x02) write_spi(MAG_IF_2_ADDRESS, 0x4C) write_spi(MAG_IF_1_ADDRESS, 0x42) write_spi(MAG_CONF_ADDRESS, 0x05) write_spi(MAG_IF_0_ADDRESS, 0x00) write_spi(CMD_REG_ADDRESS, 0x1A) while True: data_bytearray = read_spi_multiple(DATA_START_ADDRESS, DATA_NUM) mag_x = convert_mag_data(data_bytearray[0:2]) mag_y = convert_mag_data(data_bytearray[2:4]) mag_z = convert_mag_data(data_bytearray[4:6]) gyro_x = convert_gyro_data(data_bytearray[8:10]) gyro_y = convert_gyro_data(data_bytearray[10:12]) gyro_z = convert_gyro_data(data_bytearray[12:14]) acc_x = convert_acc_data(data_bytearray[14:16]) acc_y = convert_acc_data(data_bytearray[16:18]) acc_z = convert_acc_data(data_bytearray[18:20]) sensor_time = convert_sensor_time(data_bytearray[20:23]) print(f"Mag X: {mag_x:.2f}, Mag Y: {mag_y:.2f}, Mag Z: {mag_z:.2f} [μT]") print( f"Gyro X: {gyro_x:.2f}, Gyro Y: {gyro_y:.2f}, Gyro Z: {gyro_z:.2f} [deg/s]" ) print(f"Acc X: {acc_x:.2f}, Acc Y: {acc_y:.2f}, Acc Z: {acc_z:.2f} [g]") print(f"Sensor time: {sensor_time:.3f} [sec]") print(data_bytearray) print("") time.sleep(1) 実行すると1秒間隔でセンサの情報を取ってくれました。 Mag X: 29.56, Mag Y: -4.94, Mag Z: -4.19 [μT] Gyro X: -0.12, Gyro Y: 0.18, Gyro Z: 0.06 [deg/s] Acc X: 0.99, Acc Y: -0.02, Acc Z: 0.05 [g] Sensor time: 156.778 [sec] b'\xd9\x01\xb1\xff\xbd\xff\rZ\xfe\xff\x03\x00\x01\x00)?\xc5\xfe\x03\x03\xecV=' Mag X: 30.56, Mag Y: -5.94, Mag Z: -3.81 [μT] Gyro X: -0.06, Gyro Y: 0.18, Gyro Z: 0.18 [deg/s] Acc X: 0.99, Acc Y: -0.02, Acc Z: 0.05 [g] Sensor time: 157.816 [sec] b'\xe9\x01\xa1\xff\xc3\xff\rZ\xff\xff\x03\x00\x03\x00$?\xc4\xfe\xe5\x02\xe7\xbe=' 参考URL https://blog.boochow.com/article/mp-stm-nucleo-l432kc.html https://www.st.com/ja/microcontrollers-microprocessors/stm32l4x6.html https://ja.wikipedia.org/wiki/MicroPython https://github.com/micropython/micropython https://docs.micropython.org/en/latest/develop/gettingstarted.html https://www.st.com/resource/en/datasheet/DM00284207.pdf https://www.mouser.jp/new/bosch/bosch-sensortec-bmx160/

最適化について一般的なことを説明し、その中でも勾配降下法なども説明しつつ、ガウス・ニュートン法の導出について簡単に説明する。ガウス・ニュートン法によるモデルのパラメータ推定についても説明し、具体的な問題に適用した場合も実演する。 この記事を書いた動機 書籍やネット上で最適化についての資料はたくさんあるが、ガウス・ニュートン法についてもっと簡単に説明できないかと思ったため。 間違い等あれば指摘お願いいたします。 最適化とは 最適化とはいろいろな制約の中で、複数の選択肢から最適なものを選ぶことである。 最適化は様々な分野で用いられる。 企業の経営戦略:どの事業に投資すれば売り上げ、利益を最大化できるか？ 設計:強度があって軽い製品にするための材料はどうするべきか？ 生産計画:限られた時間とコストの中で、できるだけ多くの製品を作りたい 制御:システムが目的通りの挙動するように最適化する モデル化:現実の現象を最もよく表すモデルのパラメータは何か？ 機械学習:目的関数を最小化するモデルのパラメータは何か？ 最適化問題 これら最適化問題について数学的に考えてみる。 ある変数$x$がある範囲$S$の中で動くとき、$x$についての関数$f(x)$を最小化する$x$を探索する問題は \begin{align} f(x)\rightarrow min\\ x\in S \end{align} と書ける。$x$を「決定変数」、$f(x)$を目的の達成度を表す「目的関数」、決定変数に課される条件$S$を「制約条件」と呼ぶ。変数が1次元の時の場合のイメージを図にした。 最適化問題の種類 最適化問題はいろいろな種類の問題があり、分類の方法もいろいろある。詳しくはWikipediaの最適化問題のページを参照 決定変数について分類 連続最適化問題:連続的な値を取る場合 離散最適化問題:0か1、離散的な値を取る場合 目的関数についての分類 線形計画問題:目的関数が線形関数 非線形計画問題:目的関数が線形でない 最適化問題の解法の種類 またその最適化問題を解く手法(アルゴリズム)にも様々な種類がある。 導関数を使うアルゴリズム 最急降下法(勾配降下法) ニュートン法 準ニュートン法 ガウス・ニュートン法(最小二乗法でパラメータを求める場合) 導関数を使わないアルゴリズム アニーリング法 タブー探索法 遺伝アルゴリズム 粒子群最適法 ベイズ最適化 各アルゴリズムのPythonライブラリ 有名なアルゴリズムについてはPythonのモジュールやサンプルコードあるので自分で実装する必要はないが、そのアルゴリズムについて基礎的なことは理解していた方がいい。 ガウス・ニュートン法(に似た、より効率の良いアルゴリズム):scipy.optimize.least_squares 遺伝アルゴリズム:DEAP ベイズ最適化:GpyOpt ガウス・ニュートン法の導出 ここから各種最適化アルゴリズムの中でも、導関数を使った最適化の1つであるガウス・ニュートン法を導出し、その後ガウス・ニュートン法を使ったモデルのパラメータ推定について実験してみる。 ガウス・ニュートン法はある現象をパラメータを含んだモデル関数によってモデル化するとき、効率的にパラメータの探索するアルゴリズムである。 ここからその手法を導出していく。本来は、最適化問題についてはもっと厳密な議論が必要だが、今回はイメージ優先で厳密な議論はしない。正確な議論を知りたい場合は参考文献や論文など参照してほしい。 1変数関数の最小化 まず、1変数関数の最小化について説明する。 最急降下法 まずもっとも単純な最急降下法(勾配降下法)について述べる。 まず上図のような関数$f(x)$があったとする。 反復数のカウンターを$k=0$とする。下図のように初期値を$x_0$とし適当に置く。 以下3~4を反復する。 $x_k$での勾配$f'(x_k)$を計算する 下図のように$x_{k+1}=x_k-\alpha f'(x_k)$と更新する。 $||f'(x_k)||$が十分小さくなるまで続ける。 $\alpha$は更新の度合いを決めるパラメータである。 ここで、効率的に最小値を求める方法を考える。最小値ではないかもしれないが、少なくとも$f(x)$が極小値となる点を探す。それはつまり下図のように$f(x)$の微分が0、つまり f'(x)=0 となる$x$が分かればいい。 そこで、4で$x_{k+1}$を求める際に、$f'(x_{k+1})$が0に近くなるようにαを定めると、効率的に最小値を求められる。つまり$f'(x_{k+1})=0$という式を解くことになる。 一般に$g(y)=0$を解くことを方程式を解くという。数値的に(数値的：プログラムで計算して)方程式を解く方法として、ニュートン法というものがある。つまり最急降下法はニュートン法を用いることで、効率的に最小値が求められるようになる。 ニュートン法 ここで数値的に方程式を解くニュートン法を、下図のような関数$g(y)$について、$g(y)=0$を解いてみて説明する。 1.$y_0$を適当に置く 2.以下3~4を繰り返す 3.y_kの曲線上での接線を求める 4.その接線と横軸との交点を求めて$y_{k+1}$とする 接線は2点$(y_k,g(y_k))$、$(y_{k+1},0)$を通り、その傾きは$g(y_k)/(y_k-y_{k+1})$と表せられる。また同時に点$y_k$の微分$g'(y_k)$もその点での接線の傾きなので、$g(y_k)/(y_k-y_{k+1})=g'(y_k)$となる。よって$y_{k+1}$を求めるには、 y_{k+1}=y_k-\frac{g(y_k)}{g'(y_k)} を計算すればいい。これを反復し、$g(y_d)=0$となる$y_d$を探索する。 最急降下法にニュートン法を用いる ここまでを振り返ってみると、 勾配法では、 x_{k+1}=x_k-\alpha f'(x_k) と$x_{k+1}$を更新する。 $f'(x)=0$となる$x$を探せば効率的に最小値が探索できる。 また、$\alpha$は更新の度合いを決めるパラメータで自由に決められる。 ニュートン法では、 $g(y)=0$となる$y$を探す y_{k+1}=y_k-\frac{g(y_k)}{g'(y_k)} と$y_{k+1}$を更新する。 この式を見比べれば、$g(y)$にあたるものが$f'(x)$なので、$g'(y)$にあたるものは$f''(x)$である。よって効率的に最小値を探索するには、$\alpha=1/f''(x_k)$とし、 x_{k+1}=x_k- \frac{f'(x_k)}{f''(x_k)} と$x_k$を更新すれば$f'(x)=0$となる$x$を探せ、効率的に最小値が探索できることが分かる。(しかしこれでは極小値ではなく極大値を探索してしまう可能性もあるが、それについては無視することとする。) 多変数関数の最小化 今まで、1変数関数の最小化を考えてきたが、ここでは多変数関数の最小化について考える。多変数関数の場合、決定変数は$(x_1,x_2,\cdots)$と複数の変数を持つことになり、最適化する関数も$f(x_1,x_2,\cdots)$と複数の値に対し、値を返す多変数関数となる。複数の変数をまとめて決定変数ベクトル$\pmb{x}=(x_1,x_2,\cdots)$と書く。最適化問題は \begin{align} f(\pmb{x})\rightarrow \min \\ \pmb{x}\in S \end{align} となる。 多変数関数の最急降下法 1変数関数の最急降下法の場合 x_{k+1}=x_k-\alpha f'(x_k) を繰り返し、その時に f'(x)=0 となる$x$が分かれば効率的に最小化できたのであった。 多変数関数を最急降下法で最小化する場合も、同様に \pmb{x}_{k+1}=\pmb{x}_k-\alpha \nabla f(\pmb{x}_k) を繰り返し、その時に \nabla f(\pmb{x})=0 となる$\pmb{x}$が分かれば効率的に最小化できる。ただし$\nabla f(\pmb{x})$は$\pmb{x}$に関する$f(\pmb{x})$の勾配ベクトルである。勾配ベクトルは$f(\pmb{x})$を$x_1,x_2,\cdots$のそれぞれで偏微分したベクトルで、 \nabla f(\pmb{x})=\left(\frac{\partial f}{\partial x_1},\frac{\partial f}{\partial x_2},\cdots\right)^T と表せる。 多変数ベクトル値関数のニュートン法 多変数関数の最急降下法を効率的に行うために、 \nabla f(\pmb{x})=0 を解く必要がある。この$\nabla f(\pmb{x})$は複数の変数を持ち、複数の値を返す、多変数ベクトル値関数なので、多変数ベクトル値関数に関するニュートン法について考えなければならない。 多変数ベクトル値関数を$\pmb{g}(\pmb{y})$とすると、その全ての値が$0$、つまり、 \begin{align} \pmb{g}(\pmb{y})&=(g_1(y_1,y_2,...),g_2(y_1,y_2,...),...)^T\\ &=(0,0,...)^T \end{align} \tag{1} を数値的に解く場合について考える。 前に述べたように1変数関数のニュートン法の場合、 y_{k+1}=y_k-\frac{g(y_k)}{g'(y_k)} を繰り返せばよかったから、多変数ベクトル値関数の場合も、正確ではないが、気持ちとしては \pmb{y}_{k+1}=\pmb{y}_k-\frac{\pmb{g}(\pmb{y}_k)}{\pmb{g}'(\pmb{y}_k)} を繰り返せばいいと考えられる。 実際、正確には、 \nabla \pmb{g}(\pmb{y}_k)\Delta\pmb{y}_k=-\pmb{g}(\pmb{y}_k) \tag{2} を満たす$\Delta\pmb{y}_k$を計算し、そして、 \pmb{y}_{k+1}=\pmb{y}_k+\Delta\pmb{y}_k を計算し、これを反復すれば式(1)を満たす$\pmb{y}$を見つけられる。 ただし、ここで$\nabla \pmb{g}(\pmb{y})$はベクトル$\pmb{y}$に対する$\pmb{g}(\pmb{y})$の勾配行列で \nabla \pmb{g}(\pmb{y})= \begin{pmatrix} \frac{\partial g_1}{\partial y_1} & \frac{\partial g_1}{\partial y_2} & \cdots \\ \frac{\partial g_2}{\partial y_1} & \cdots & \cdots\\ \cdots & \cdots & \cdots \end{pmatrix} である。さきほどの$\nabla f(\pmb{x})$については、$f(\pmb{x})$が多変数に対する1値関数だったので、$\nabla f(\pmb{x})$は勾配ベクトルとなったが、$\pmb{g}(\pmb{y})$は多変数ベクトル値関数なので、$\nabla \pmb{g}(\pmb{y})$は勾配行列となる。 一般的に、式(2)を解くことは、線形連立方程式を解くことになり、$LU$分解などの連立方程式を数値的に解くアルゴリズムが使える。 多変数関数の最急降下法にニュートン法を用いる 1変数関数と同様に考えて、$\pmb{g}(\pmb{y})$にあたるものが$\nabla f(\pmb{x})$なので、 \nabla^2 f(\pmb{x}_k)\Delta\pmb{x}=-\nabla f(\pmb{x}_k) を満たす$\Delta\pmb{x}$を用いて、 \pmb{x}_{k+1}=\pmb{x}_k+\Delta\pmb{x}_k と更新することを繰り返せば、$f(\pmb{x})$を最小化する$\pmb{x}$を見つけることができる。(ここでも極小値ではなく極大値を探索してしまう可能性もあるが、それについては無視することとする。) ただし、ここで$\nabla^2 f(\pmb{x})$はベクトル$\pmb{x}$に対する$\pmb{f}(\pmb{x})$の2回微分の行列でヘッシアン$H$といい、 \begin{align} \nabla^2 f(\pmb{x})&=H\\ &= \begin{pmatrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1\partial x_2} & \cdots \\ \frac{\partial f}{\partial x_2\partial x_1} & \cdots & \cdots\\ \cdots & \cdots & \cdots \end{pmatrix} \end{align} である。 テイラー展開を用いた多変数関数の最小化 今まで述べた多変数関数の最適化について、テイラー展開を用いて説明することもできる。数理最適化の教科書はこの説明を使うことが多い。 放物線の最小化 まず、下に凸な放物線の式($a>0$) g(x)=\frac{a}{2}x^2+bx+c\tag{3} を最小化する$x$は、この式を平方完成し、軸の位置を求めることで求まる。 \begin{align} g(x)&=\frac{a}{2}x^2+bx+c\\ &=\frac{a}{2}\left(x+\frac{b}{a}\right)^2+c-\frac{b^2}{2a} \end{align} よって$g(x)$は$x=-b/a$で最小となる。後の為に少し式変形すると ax=-b\tag{4} を満たす$x$で$g(x)$は最小になるとも言える。 テイラー展開を用いた多変数関数の最小化 多変数関数$f(\pmb{x})$を最小化したいとする。 \pmb{x}_{k+1}=\pmb{x}_k+\Delta\pmb{x}_k のように反復的に最適化することを考える。 この$f(\pmb{x})$について、更新前の点$\pmb{x}_k$に近い場所の点$\pmb{x}$での$f(\pmb{x})$の値は、多変数関数に対するテイラー展開より、微小なベクトル$\Delta\pmb{x}=\pmb{x}-\pmb{x}_k$を用いて、 f(\pmb{x})\approx f(\pmb{x}_k)+\nabla f(\pmb{x}_k)^T\Delta\pmb{x}+\frac{1}{2}\Delta\pmb{x}^T\nabla^2 f(\pmb{x}_k)\Delta\pmb{x} のように表すことができる。$\pmb{x}_k$を固定した時、この右辺を最小化したい場合、 \frac{1}{2}\Delta\pmb{x}^T\nabla^2 f(\pmb{x}_k)\Delta\pmb{x}+\nabla f(\pmb{x}_k)^T\Delta\pmb{x}\tag{5} を最小化すればいい。この式は多変数関数の2次関数で、もしこれが下に凸である場合、前節の放物線の最小化と同じように最小化できる。つまりこの式(5)を式(3)と見比べて、式(4)と同じ形の式を作ればいい。よって \nabla^2 f(\pmb{x}_k)\Delta\pmb{x}=-\nabla f(\pmb{x}_k) を満たす$\Delta\pmb{x}$で式(5)は最小となる。よって、この$\Delta\pmb{x}_k$を用いて、 \pmb{x}_{k+1}=\pmb{x}_k+\Delta\pmb{x}_k と更新すれば、$f(\pmb{x})$を最小化する$\pmb{x}$を見つけることができる。このように上で述べた結論と同じことが、テイラー展開を用いた方法で説明することができる。 多変数関数の最小二乗法 多変数関数の最小二乗法について考える。 最小二乗法とは 最小二乗法とは、測定したデータに合いそうな関数を当てはめる手法で、モデル関数による計算値と実験での計測などで得られたデータの残差の二乗和が最小になるように当てめる。これも最適化問題である。 例えば、フックの法則によれば、ばねに加える力$x$とそのときのばねの伸び$y$は、ばね係数を$k$とすると y=kx となる。様々な大きさ$x_i$の力をばねに加え、そのときのばねの伸び$y_i$を測定することで、最小二乗法でこのばね係数$k$を求めることができる。 モデル関数を$f(x)$、計測で得られた各データ点を$(x_i,y_i)$とすると、 モデル関数に変数である$x_i$を入れた計算値$f(x_i)$と実際の測定値$y_i$の差(残差)は r_i=f(x_i)-y_i である。残差はひとまとめにして残差ベクトル、$\pmb{r}$と表す。そして、二乗誤差は S=\frac{1}{2}\sum_i r_i^2 である。二乗誤差$S$が最小になる$f(x)$のパラメータを探せば、それが最も妥当なモデル関数だと思われる。 例えばモデル関数が$f(x)=ax+b$の場合、パラメータは$(a,b)$である。これは一次関数の最小二乗法で、$S$を最小にするパラメータは \begin{align} \frac{\partial S}{\partial a}=0\\ \frac{\partial S}{\partial b}=0 \end{align} の解である。これは平均と分散を用いて \begin{align} \hat{a}=\frac{cov(x,y)}{\sigma^2(x)} \\ \hat{b}=\mu_y-\hat{a}\mu_x \end{align} と計算できる。詳しくは「高校数学の美しい物語」の最小二乗法（直線）の簡単な説明のページを参照 より一般的には最小二乗法とは、パラメータ$\pmb{\lambda}$が含まれるモデル関数$f(x)$について S(\pmb{\lambda})=\frac{1}{2}\sum_i (f(x_i)-y_i)^2 となる$S(\pmb{\lambda})$の最小化問題である。ここで、モデル関数$f(x)$の変数は$x$だが、最小化の対象の$S(\pmb{\lambda})$の変数はパラメータ$\pmb{\lambda}$である点に注意すること。 多変数関数の最小二乗法に最急降下法とニュートン法を用いる 複数のパラメータ$\pmb{\lambda}$が含まれる多変数関数$f(x)$について、 計測で得られた各データ点を$(x_i,y_i)$について、 S=\frac{1}{2}\sum_i (f(x_i) - y_i)^2 を最小にするパラメータ$\pmb{\lambda}$を求めることを考える。もしこの関数が非線形だと、微分=0という式は簡単には解けない。そこで、最急降下法とニュートン法を組み合わせてパラメータを求めることができる。 前の節で述べた結果から、勾配$\nabla S(\pmb{\lambda})$とヘッシアン$H=\nabla^2 S(\pmb{\lambda})$を用いて、 H\Delta\pmb{\lambda}_k=-\nabla S(\pmb{\lambda}) を満たす$\Delta\pmb{\lambda_k}$を計算し、そして、 \pmb{\lambda}_{k+1}=\pmb{\lambda}_k+\Delta\pmb{\lambda}_k を計算し、これを反復すれば$S(\pmb{\lambda})$を最小にするパラメータ$\pmb{\lambda}$を求められる。つまりモデル式$f(x)$の最適なパラメータ$\pmb{\lambda}$が求めることができる。 ガウス・ニュートン法 ガウス・ニュートン法 - Wikipedia ただ、ヘッシアンは2回微分が必要で計算コスト高く、またその式を導出するのも手間がかかる。 ヘッシアンを近似することで計算量と手間を減らすのがガウス・ニュートン法である。 まず、勾配$\nabla S$を考えてみると、 \nabla S=\frac{\partial S}{\partial \pmb{\lambda}}=\left( \sum_i (f(x_i)-y_i)\frac{\partial f(x_i)}{\partial \lambda_1},\sum_i (f(x_i)-y_i)\frac{\partial f(x_i)}{\partial \lambda_2}, ...\right)^T で、成分表示をし、さらに残差$r$を用いれば、次のように表せる。 \begin{align} (\nabla S)_j&=\sum_i (f(x_i)-y_i)\frac{\partial f(x_i)}{\partial \lambda_j}\\ &=\sum_i r_i\frac{\partial f(x_i)}{\partial \lambda_j}\tag{6} \end{align} ここでヤコビアン$J$を導入する。ヤコビアンも行列でその成分$J_{ij}$は、 J_{ij}=\frac{\partial f(x_i)}{\partial \lambda_j} である。ヤコビアンは縦が測定数、横がパラメータ数の行列となる。 よって式(6)はヤコビアンを用いて (\nabla S)_j= \frac{\partial S}{\partial \lambda_j}=\sum_i r_i J_{ij} つまり、 \nabla S=J^T\pmb{r} と表せる。 次にヘッシアン$H=\nabla^2 S$は \begin{align} H_{jk}&=(\nabla^2 S)_{jk}\\ &=\sum_i \left(\frac{\partial f_i}{\partial \lambda_k} \frac{\partial f_i}{\partial \lambda_j}+ (f_i-y_i)\frac{\partial^2 f_i}{\partial \lambda_j\lambda_k}\right)\\ &=\sum_i \left( J_{ik}J_{ij}+r_i\frac{\partial^2 f_i}{\partial \lambda_j\lambda_k}\right) \end{align}\tag{7} となる。 しかし式(7)の最後の式の右辺第2項目の2回微分は計算コストが高く、プログラムで式を実装するのも手間がかかる。ガウス・ニュートン法ではこの2回微分を無視し、ヘッシアンを H_{jk}\approx\sum_i J_{ik}J_{ij} つまり H\approx J^TJ と近似する。 以上からガウス・ニュートン法では J^TJ\Delta\pmb{\lambda}_k=-J^T\pmb{r} を満たす$\Delta\pmb{\lambda_k}$を計算し、そして、 \pmb{\lambda}_{k+1}=\pmb{\lambda}_k+\Delta\pmb{\lambda}_k を計算する。これを反復することで、$S(\pmb{\lambda})$を最小にするパラメータ$\pmb{\lambda}$を求められる。つまりモデル式$f(x)$の最適なパラメータ$\pmb{\lambda}$が求めることができる。 ガウス・ニュートン法によるモデルのパラメータ推定 ガウス・ニュートン法について具体的な実装と図示を見ながら確認する。 ガウス・ニュートン法 - Wikipediaに記載されているミカエリス・メンテン式の例を借りる。詳しい内容はリンク先を参照。 (出典：［ガウス・ニュートン法」『フリー百科事典　ウィキペディア日本語版』。2020年10月29日 (木) 10:09　UTC、URL: https://ja.wikipedia.org) 酵素の基質濃度と反応速度の関係のデータがある。 i 1 2 3 4 5 6 7 [S] 0.038 0.194 0.425 0.626 1.253 2.500 3.740 v 0.050 0.127 0.094 0.2122 0.2729 0.2665 0.3317 これに対し、ミカエリス・メンテン式をモデル関数としてフィッティングする。 v=\frac{V_{max}}{K_M+[S]} つまり二乗誤差 S=\frac{1}{2}\sum_i (v([S]_i) - v_i)^2 を最小化するようなミカエリス・メンテン式のパラメータ$V_{max}$と$K_M$を求める。 初期値 (V_{max},K_M)=(0.2,0.8) から出発し、収束させる。収束の推移を等高線図とプロットと矢印で描写する。 収束はしているが、等高線と垂直に降下していない、とくに最初の降下で顕著である。これはヘッシアンに近似を用いているためと考えられる。 また、二乗誤差の収束の様子を図示する。 急速に、誤差関数が収束していることが分かる。 また、測定点と算出したV_maxとK_Mでミカエリス・メンテンのモデル式を図示する。 パラメータフィッティングと図を算出するため作ったPythonコードをこのページの最後の捕捉に示す。 ライブラリの利用 実際にガウス・ニュートン法を用いて多変数モデル式のパラメータを決定したい場合、自分で実装するより、上のことを分かった上で、ライブラリを使う方がいい。 scipy.optimize.least_squares Scipyのoptimeze.least_squaresではガウス・ニュートン法よりより効果的な手法が使え、さらに制約条件についても考慮してくれるなど、機能が充実している。 この関数には手法(method)とヤコビアン(jac)を指定することができ、いくつかの手法を使う場合は、ヤコビアンを算出する関数を代入することで高速に探索ができる。もしヤコビアンを指定しない場合は、この関数内部で自動的に有限差分法でヤコビアンを算出する。自動的にヤコビアンを計算してくれるので便利に思えるが、有限差分でヤコビアンを計算すると、パラメータの数に比例して計算時間が増えてしまう欠点がある。 ここで有限差分法でのヤコビアンの計算方法について説明する。 パラメータ$\pmb{\lambda}=(\lambda_1,\lambda_2,\cdots)$でのモデル関数の計算値を$f(x)(\lambda_1,\lambda_2,\cdots)$、微小な$\lambda$の変化を$\Delta\lambda$とする。そうすると、有限差分でヤコビアンを計算する場合、最初のパラメータ$\lambda_1$に対する偏微分は \begin{align} J_{i1}&=\frac{\partial f(x_i)}{\partial \lambda_1}\\ &\approx\frac{f(x_i)(\lambda_1+\Delta\lambda,\lambda_2,\cdots)-f(x_i)(\lambda_1-\Delta\lambda,\lambda_2,\cdots)}{2\Delta\lambda}\\ \end{align} と計算できる。これを各パラメータ$(\lambda_1,\lambda_2,\cdots)$について行わなければならない。よって有限差分法でヤコビアンを計算すると、計算時間はパラメータの数に比例して増えていく。 もし、$\frac{\partial f(x_i)}{\partial \lambda_j}$が解析的に計算できるのであれば、その解析式をヤコビアンの関数としてプログラムで実装し、このoptimeze.least_squaresに引数として代入すれば高速化できることになる。そしてその高速化の効果はパラメータが増えれば増えるほど大きくなることになる。 さらにこの関数の詳しい利用方法はドキュメントなどを参照。 補足 ミカエリス・メンテン式のパラメータフィッティングと図示のコード。 import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import LogNorm from matplotlib import ticker,cm def michaelis(conce,V_max,K_m): """ ミカエリス・メンテン式。反応率を算出する。 """ return V_max*conce/(K_m+conce) def compute_j(matrix_con,p_x=np.array([[0.5],[0.5]])): """ モデル関数について微分し、ヤコビアンを計算する。 Args: matrix_con : 基質濃度の1次元配列 """ # 要素が0の縦2、横1の行列を用意 dx_base=np.zeros([2,1]) # 微小距離を定義 dx=1e-3 # ヤコビアンのリストを設定 jacobian=[] # パラメータの数分ループ for i in range(2): dx_array=dx_base.copy() dx_array[i,0]=dx # p_xから微小距離だけ移動したベクトルを作る。 p_x_d=p_x+dx_array # ミカエリス・メンテン式で反応率を計算する。 v_p=michaelis(matrix_con,p_x_d[0,0],p_x_d[1,0]) p_x_d=p_x-dx_array v_m=michaelis(matrix_con,p_x_d[0,0],p_x_d[1,0]) jacobian.append((v_p-v_m)/dx/2) return np.array(jacobian).T # 測定データを用意 matrix_con=np.array([0.038,0.194,0.425,0.626,1.253,2.500,3.740]) cons_rate=np.array([0.050,0.127,0.094,0.2122,0.2729,0.2665,0.3317]) # 初期値を与える。 p_x=np.array([[0.8],[0.2]]) # 最小値に向かう様子の推移を記録するリストの用意 p_x_tr=[p_x.copy()] s_tr=[] for i in range(10): # ヤコビアンの算出 jacobian=compute_j(matrix_con,p_x) # 残差の算出 res=michaelis(matrix_con,p_x[0],p_x[1]) - cons_rate # この後の行列積のため、縦測定数、横1の行列に変換 res=res[:,np.newaxis] # 疑似ヘッシアンを算出 hessian = jacobian.T @ jacobian # 最小値に向かう移動ベクトルを算出 d_p_x= - np.linalg.pinv(hessian) @ (jacobian.T @ res) # 移動 p_x = p_x+d_p_x p_x_tr.append(p_x.copy()) s_tr.append(np.sum(res**2)/2) # 図示の準備 array_1=np.arange(0,1.2,0.01) array_2=np.arange(0,1.2,0.01) mat_1,mat_2=np.meshgrid(array_1,array_2) length=matrix_con.shape[0] sse=np.zeros_like(mat_1) for i in range(length): calc_cons=michaelis(matrix_con[i],mat_1,mat_2) a_res=cons_rate[i]-calc_cons sse += a_res**2 # 図示 plt.rcParams["font.size"] = 18 fig = plt.figure() ax1 = fig.add_subplot(111) lev_exp = np.arange(np.floor(np.log10(sse.min())-1), np.ceil(np.log10(sse.max())+1)) levs = np.power(10, lev_exp) # 等高線図を図示する。カラーマップは対数スケールとする。 cs=ax1.contourf(mat_1,mat_2,sse,levs, location=ticker.LogLocator(),norm=LogNorm()) c_bar=fig.colorbar(cs) # 収束の推移をプロットと矢印で描写する。 for i in range(9): ax1.plot(p_x_tr[i][0],p_x_tr[i][1],'om') ax1.annotate('',xy=p_x_tr[i+1],xytext=p_x_tr[i], arrowprops=dict(shrink=0, width=1, headwidth=8, headlength=10, connectionstyle='arc3', facecolor='k', edgecolor='k')) ax1.plot(p_x_tr[10][0],p_x_tr[10][1],'om') ax1.set_xlabel('V_max') ax1.set_ylabel('K_M') c_bar.set_label('reaction rate') plt.show() # 二乗誤差の収束の様子を図示 fig=plt.figure() ax2=fig.add_subplot(111) ax2.plot(s_tr) ax2.set_yscale('log') ax2.set_xlabel('step') ax2.set_ylabel('sum of square error') plt.show() # 測定点と # 算出したV_maxとK_Mでミカエリス・メンテンのモデル式を図示 conce_array=np.arange(0,4,0.1) rr_array=michaelis(conce_array,p_x[0,0],p_x[1,0]) fig=plt.figure() ax3=fig.add_subplot(111) ax3.plot(conce_array,rr_array) ax3.plot(matrix_con,cons_rate,'o') ax3.set_xlabel('matrix concentration') ax3.set_ylabel('reaction rate') plt.show()

1. はじめに Python三次元ポアソン方程式による静電場解析により，基本的課題として空間上の点電荷の電場解析が可能なことを先の投稿において報告した(1)． 本稿では，応用課題として円板状電荷を対象に下記を実施する． (1)三次元ポアソン方程式を差分化（再掲） (2)繰り返し法によるプログラムを作成． ①メッシュ幅は1.0と0.5を選択できるようにする． ②電位を計算 ③得られた電位から電場を計算 (3)計算結果を理論値と比較することにより解析の妥当性を検証 2. 例題の問題設定 各軸のメッシュ数は50とし，中心に帯電した円板状電荷を配置する． 　　　　　図1　問題設定-帯電した円板状電荷の配置 3. ポアソン方程式の離散化(再掲) 下記に手順を記す． 　ここで注意が必要な点はρは電荷密度であることである． 4. プログラム import pprint from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt import numpy as np import matplotlib.cm as cm delta_L=1.0 #delta_L=0.5 LN=50 HLN=int(LN/2) nx = LN ny = LN nz = LN xmin = 0 xmax = LN*delta_L ymin = 0 ymax = LN*delta_L zmin = 0 zmax = LN*delta_L p = np.zeros((nz, ny, nx)) pd = np.zeros((nz, ny, nx)) charge = np.zeros((nz, ny, nx)) x = np.linspace(xmin, xmax, nx) y = np.linspace(ymin, ymax, ny) z = np.linspace(zmin, zmax, nz) eps0=1 y00=1 Q1=1 q=Q1/(delta_L**3) for i in range(0, LN): for j in range(0,LN): rr=np.sqrt((i-HLN)**2+(j-HLN)**2) if rr <= 3: charge[i,j,HLN] = q#/(delta_L**3)#*(delta_L**2) print("Riron") q1=Q1/(delta_L**2) for i in range(0,8): y00=i*delta_L x2=q1*(np.sqrt((3*delta_L)**2+y00**2)-y00)/2 print('z=',y00+25,' V=',x2) for i in range(0,8): y00=i*delta_L x2=q1*(1-y00/np.sqrt((3*delta_L)**2+y00**2))/2 print('Ez=',y00+25,'->',x2) for I in range(500): #2000 pd = p.copy() p[1:-1,1:-1,1:-1] = (pd[2:,1:-1,1:-1] + pd[:-2,1:-1,1:-1] + pd[1:-1,2:,1:-1] + pd[1:-1,:-2,1:-1] + pd[1:-1,1:-1,2:] + pd[1:-1,1:-1,:-2] + (delta_L**2)*charge[1:-1, 1:-1, 1:-1]) /6 p[0][:][:]=0 p[nx-1][:][:]=0 p[:][0][:]=0 p[:][ny-1][:]=0 p[:][:][0]=0 p[:][:][nz-1]=0 def plot2D(x, y, p): fig = plt.figure(figsize=(11, 7), dpi=100) ax = fig.gca(projection='3d') X, Y = np.meshgrid(x, y) surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False) ax.view_init(30, 225) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('p') plot2D(x,y,p[:,:,HLN]) def plot_surface(pp): nx = LN ny = LN fig = plt.figure(figsize=(11,7), dpi=100) fig xmin = 0 xmax = LN-1 ymin = 0 ymax = LN-1 x = np.linspace(xmin, xmax, nx) y = np.linspace(ymin, ymax, ny) X, Y = np.meshgrid(x, y) plt.contourf(X, Y, pp, alpha=0.1, cmap=cm.viridis) plt.colorbar() plt.xlabel('X') plt.ylabel('Y') plot_surface(p[:,:,HLN]) plt.xlim([HLN-5,HLN+5]) plt.ylim([HLN-5,HLN+5]) plt.style.use('ggplot') plt.rcParams["axes.facecolor"] = 'white' fig = plt.figure() ax = fig.gca(projection='3d') L=9 NLS=HLN#25 NLY=2 X, Y, Z= np.meshgrid(np.arange(0, L,1), np.arange(0, L,1), np.arange(0, NLY,1)) Ex = np.zeros([L,L,abs(NLY)]) Ey = np.zeros([L,L,abs(NLY)]) Ez = np.zeros([L,L,abs(NLY)]) for i in range(L): for j in range(L): for k in range(0,NLY): ii=i+HLN-4#21 jj=j+HLN-4#21 kk=k+NLS Ex[i,j,k]=-(p[ii,jj+1,kk]-p[ii,jj-1,kk])/2/delta_L Ey[i,j,k]=-(p[ii+1,jj,kk]-p[ii-1,jj,kk])/2/delta_L # Ez[i,j,k]=-(p[ii,jj,kk+1]-p[ii,jj,kk-1])/2 Ez[i,j,k]=-(p[ii,jj,kk+1]-p[ii,jj,kk])/delta_L print("Keisan") for i in range(20,33): print('Vz=',25+(i-25)*delta_L,' ',p[HLN,HLN,i]) i=HLN j=HLN for k in range(HLN,HLN+6): print("Ez=",k,"-> ",-(p[i,j,k+1]-p[i,j,k-1])/2) print("Ez=",25+(k-25)*delta_L,"-> ",-(p[i,j,k+1]-p[i,j,k])/delta_L) ax.quiver(X,Y,Z, Ex, Ey, Ez, color='red',length=1.0*delta_L**2, normalize=False) for i in range(NLY): X1,Y1, Z1=4,4,i ax.scatter3D(X1,Y1,Z1,"o", color='blue') plt.grid() plt.draw() plt.show() 5. 計算結果-delta_T=1.0の場合 図2(a)にz=25面の電位分布の2Dﾌﾟﾛｯﾄ図を，図2(b)にz=25（電荷を配置）とz=26の電場図を示す． 　　　　　　　　　　　　　　図2(a)　電位分布　　　　　　　　　　　　　　　 　　　　　　　　　　　　　　図2(b)　電場 　また，図3にz=23,24,25,26層の電場を示す．上下対象になっていることが分かる． 　　　　　　　　　図3 z=23,24,25,26層の電場 6. 解析の妥当性検証 6.1円板状電荷周りの電位，電場の理論的計算 図4に理論計算のためのモデルを示す． 　　図4 理論計算のためのモデル 　理論値は下記で計算できる． 6.2 電位，電場の理論値と計算値の比較 表1にdelta_L=1.0の場合の電位，電場の理論値と計算値の比較を示す． 計算値は理論値にほぼ一致している． 　　　表1 delta_L=1.0の場合の電位，電場の理論値と計算値の比較 参考文献 (1) Python三次元ポアソン方程式による静電場解析 (1)点電荷の電場解析 Qiita

はじめに 2021年9月15日（火）に、データサイエンティスト検定を受験しました。この検定は、2021年9月11日（土）からスタートした新たな試験です。本記事では検定の概要と、私が実際に受験して感じた所感を書いていきたいと思います。 私は実務経験者ですが、初学者や実務未経験の方にも参考になるように書きました。 今後、受検を考えている方に少しでもお役立てできれば幸いです。 なぜ受検したか？ 私は30歳手前で未経験からデータ分析者になり、業務で必要な知識をその都度気合で習得してきました。 体系的に勉強した経験がないため、どの適度の知識があるか客観的に計測したかったのが一番の動機です。 データサイエンティスト検定って何？ 「そんなのは知ってるよ！」という方は読み飛ばしてください。 3行でまとめると ・データサイエンティストとしての実務能力と知識を証明するための試験だよ ・難易度は見習いレベルの基本的な内容が多いよ ・とはいっても、カバーすべき領域は結構広いよ 『データサイエンティスト検定™ リテラシーレベル』（略称：DS検定™ ★）とは、アシスタント・データサイエンティスト（見習いレベル：★）と数理・データサイエンス教育強化拠点コンソーシアムが公開している数理・データサイエンス・AI（リテラシーレベル）におけるモデルカリキュラムを総合し、実務能力と知識を有することを証明する試験です。 参考：https://www.datascientist.or.jp/dskentei/ 検定概要 受験資格 なし 実施概要 選択式問題 全国の試験会場で開催(CBT) 問題数80問程度 試験時間90分 出題範囲 スキルチェックリストの3カテゴリ（データサイエンス力、データエンジニアリング力、ビジネス力）の★１（見習いレベル）相当と数理・データサイエンス・AI（リテラシーレベル）におけるモデルカリキュラムを総合した範囲 参考：https://www.datascientist.or.jp/dskentei/ 試験範囲 参考：https://www.datascientist.or.jp/dskentei/ 検定の対象 本検定の受験対象者は以下のような方を想定しています。 ・データサイエンティスト初学者 ・これからデータサイエンティストを目指すビジネスパーソン ・データサイエンティストに興味を持つ大学生や専門学校生など 参考：https://www.datascientist.or.jp/dskentei/ 私のバックグラウンド ・文系未経験からデータサイエンティストを目指す ・データ分析歴4年 ・データサイエンティスト系の資格は、特に持っていません 詳しくは下記を参照ください 28歳、一念発起して未経験からデータサイエンティストを志した2年間の軌跡 受検に受けた勉強 本検定に向けて下記を1周しました。 ・公式テキスト：5時間 ・対策アプリ：4時間 公式テキストの模試（45問）とアプリのオリジナル問題（90問）を解き、不正解だった分野を中心に勉強しました。 検定を終えての感想 試験内容の詳細は載せることができませんのでざっくり。あくまでも私個人の感想です。 データサイエンス力 ・公式テキストの語句の理解がきちんとできれば問題なし ・基礎統計や高校レベルの数学力があれば問題なし ・統計検定やE資格の勉強した経験がある人にとっては楽勝かも。 偉そうに書いてますが、私は全然楽勝ではなかったです。 データエンジニア力 ・SQL・BI・機械学習・DLの問題は即答 ・データインフラ系に苦戦 ・他の二つに比べると難易度が高め 実務で経験してきたところは割とできましたが、未経験＆勉強不足の箇所は何問か落としてしまいました。 ビジネス力 ・ビジネスというより、国語力や常識力が必要 ・データサイエンス関係なくねって問題もちらほら ・公式テキストやアプリのオリジナル問題の方が良問 実務未経験者への配慮がある一方で、もう少し実務よりな出題があっても良かったのではという印象です。 で、何点取れたの？ 内緒です。合格していることを祈ります。 やっておいた方がよさそうなこと 実務経験者 私のような実務からスタートした人間は、教科書的な基礎統計の数式理解や実務で扱ったことがない分野がどうしても弱いので、その部分は素直にがっつり勉強しておいた方がよさそうです。 実務未経験者・初学者 実務未経験者の方は、公式テキストやアプリの対策問題をがっつりやり込むのが良さそうです。 統計、SQL、機械学習・DL、サーバー、データベース、セキュリティ等々についてやや広めの範囲から浅めの内容が出題されました。専門書でがっつりやるよりも試験に特化したコンテンツに取り組むほうが効率が良さそうです。 実務未経験者に不利な問題は多少捨てても合格ラインには達すると思われます。 これからデータサイエンティストを目指す方へ Twitter等では「超簡単で余裕だったw」「ノー勉だったけど9割取れた」みたいな感想がちらほらみかけますが、「これからデータサイエンティストを目指す初学者」にとっては、決して簡単ではありません。 理由としては、扱う範囲がとても広いためです。 統計、SQL、機械学習・DL、サーバー、データベース、セキュリティ等々を初学者が体系的に勉強しようとすると膨大な時間がかかります。 本検定のために勉強し合格を目標にすることは、データサイエンティストになるためにプラスにはなりますが、近道ではないです。 職業としてデータサイエンティストを目指される方は、受検勉強に膨大な時間を費やすより、一度実務の現場に飛び込んで経験した方が、圧倒的に成長できる思います。 理想論ですが、受検勉強と実務を並行してどっちも頑張るのがおすすめです。 おわりに 幅広い分野についてバランスよく出題されていて、データサイエンティストとしてどの部分に課題があるかを確認する良い機会でした。 公式テキスト等できちんと勉強していれば十分対応可能で受検のハードルは高くないの興味がある方は是非チャレンジしてほしいです。

はじめに 2021年9月14日（火）に、データサイエンティスト検定を受検しました。この検定は、2021年9月11日（土）からスタートした新たな試験です。本記事では検定の概要と、私が実際に受検して感じた所感を書いていきます。 私は実務経験者ですが、初学者や実務未経験の方にも参考になるかと思います。 今後、受検を考えている方に少しでもお役立てできれば幸いです。 なぜ受検したか？ 私は30歳手前で未経験からデータサイエンティストを目指し、業務で必要な知識をその都度気合で習得してきました。 体系的に勉強した経験がないため、どの程度の知識があるか客観的に知りたかったのが一番の動機です。 データサイエンティスト検定って何？ 「そんなのは知ってるよ！」という方は読み飛ばしてください。 3行でまとめると ・データサイエンティストとしての実務能力と知識を証明するための試験だよ ・難易度は見習いレベルの基本的な内容が多いよ ・とはいっても、カバーすべき領域は結構広いよ 『データサイエンティスト検定™ リテラシーレベル』（略称：DS検定™ ★）とは、アシスタント・データサイエンティスト（見習いレベル：★）と数理・データサイエンス教育強化拠点コンソーシアムが公開している数理・データサイエンス・AI（リテラシーレベル）におけるモデルカリキュラムを総合し、実務能力と知識を有することを証明する試験です。 参考：https://www.datascientist.or.jp/dskentei/ 検定概要 受験資格 なし 実施概要 選択式問題 全国の試験会場で開催(CBT) 問題数80問程度 試験時間90分 出題範囲 スキルチェックリストの3カテゴリ（データサイエンス力、データエンジニアリング力、ビジネス力）の★１（見習いレベル）相当と数理・データサイエンス・AI（リテラシーレベル）におけるモデルカリキュラムを総合した範囲 参考：https://www.datascientist.or.jp/dskentei/ 試験範囲 参考：https://www.datascientist.or.jp/dskentei/ 検定の対象 本検定の受験対象者は以下のような方を想定しています。 ・データサイエンティスト初学者 ・これからデータサイエンティストを目指すビジネスパーソン ・データサイエンティストに興味を持つ大学生や専門学校生など 参考：https://www.datascientist.or.jp/dskentei/ 私のバックグラウンド ・文系未経験からデータサイエンティストを目指す ・データ分析歴4年 ・データサイエンティスト系の資格は、特に持っていません 詳しくは下記を参照ください 28歳、一念発起して未経験からデータサイエンティストを志した2年間の軌跡 受検に受けた勉強 本検定に向けて下記を1周しました。 ・公式テキスト：5時間 ・対策アプリ：4時間 公式テキストの模試（45問）とアプリのオリジナル問題（90問）を解き、不正解だった分野を中心に勉強しました。 検定を終えての感想 試験内容の詳細は載せることができませんのでざっくり。あくまでも私個人の感想です。 データサイエンス力 ・公式テキストの語句の理解がきちんとできていれば問題なし ・基礎統計や高校レベルの数学力があれば問題なし ・統計検定やE資格の勉強した経験がある人にとっては楽勝かも。 偉そうに書いてますが、私は全然楽勝ではなかったです。 データエンジニア力 ・SQL,BI,機械学習,深層学習,自然言語処理の問題は即答 ・データインフラ系に苦戦 ・他の二つに比べると難易度が高め 実務で経験してきたところは割とできましたが、未経験＆勉強不足の箇所は何問か落としてしまいました。 ビジネス力 ・ビジネスというより、国語力や常識力が必要 ・データサイエンス関係なくねって問題もちらほら ・公式テキストやアプリのオリジナル問題の方が良問 実務未経験者への配慮がある一方で、もう少し実務よりな出題があっても良かったのではという印象です。 で、何点取れたの？ 内緒です。合格していることを祈ります。 やっておいた方がよさそうなこと 実務経験者 私のような実務からスタートした人間は、教科書的な基礎統計の数式理解や実務で扱ったことがない分野がどうしても弱いので、その部分は素直にがっつり勉強しておいた方がよさそうです。 実務未経験者・初学者 実務未経験者の方は、公式テキストやアプリの対策問題をがっつりやり込むのが良さそうです。 統計、SQL、機械学習・DL、サーバー、データベース、セキュリティ等々についてやや広めの範囲から浅めの内容が出題されました。専門書でがっつりやるよりも試験に特化した公式テキストやアプリを中心に取り組むほうが効率が良さそうです。 実務未経験者に不利な問題は多少捨てても合格ラインには達すると思われます。 これからデータサイエンティストを目指す方へ Twitter等では「超簡単で余裕だったw」「ノー勉だったけど9割取れた」みたいな感想をちらほら見かけますが、「これからデータサイエンティストを目指す初学者」にとっては、決して簡単ではありません。 理由としては、扱う範囲がとても広いためです。 統計、SQL、機械学習・DL、サーバー、データベース、セキュリティ等々を初学者が体系的に勉強しようとすると膨大な時間がかかります。 本検定のために勉強し合格を目標にすることは、データサイエンティストになるためにプラスにはなりますが、近道ではないです。 職業としてデータサイエンティストを目指される方は、受検勉強に膨大な時間を費やすより、一度実務の現場に飛び込んで経験した方が、圧倒的に成長できる思います。 理想論ですが、受検勉強と実務を並行してどっちも頑張るのがおすすめです。 おわりに 幅広い分野についてバランスよく出題されていて、データサイエンティストとしてどの部分に課題があるかを確認する良い機会でした。 公式テキスト等できちんと勉強していれば十分対応可能で受検のハードルは高くないの興味がある方は是非チャレンジしてみてください。

AtCoder Beginner Contest 217 バーチャル参戦記 87:39でABCDE完. とはいえ、本番だったらさすがにA問題のRE2はしなかったと思うので、実質は77:39と思うと、参加できてたらパフォ1225だった模様. 参加できなくてよかったですね(笑). ABC217A - Signed Difficulty 3分で突破、RE2. 書くだけ. 全く動作確認をせず突っ込んで死んだ(笑). S, T = input().split() if S < T: print('Yes') else: print('No') ABC217B - AtCoder Quiz 2分で突破. どうとでも書けるだろうけど、set が一番楽かなと思ってこうなった. S1 = input() S2 = input() S3 = input() result = {'ABC' , 'ARC' , 'AGC' , 'AHC'} result.remove(S1) result.remove(S2) result.remove(S3) print(*result) ABC217C - Inverse of Permutation 3分半で突破. 簡単すぎて、自分が題意を勘違いしているんじゃないかと心配になった. N, *p = map(int, open(0).read().split()) Q = [None] * N for i in range(N): Q[p[i] - 1] = i + 1 print(*Q) ABC217E - Sorting Queries 23分?で突破. 考察に時間がかかったけど、実装はすぐでしたね. 一つのデータ構造でどう実現するかを考えていたけど、無理そうだなと思って2つならと考えたら後はすぐでしたね. from sys import stdin from collections import deque from heapq import heappop, heappush readline = stdin.readline Q = int(readline()) h = [] d = deque([]) result = [] for _ in range(Q): query = readline() if query[0] == '1': _, x = map(int, query.split()) d.append(x) elif query[0] == '2': if len(h) == 0: result.append(d.popleft()) else: result.append(heappop(h)) elif query[0] == '3': for x in d: heappush(h, x) d = deque([]) print(*result, sep='\n') ABC217D - Cutting Woods 40分半?で突破. 座標圧縮と2つのセグ木で書いたけど、D問題にしては牛刀すぎるように感じるので、出題者の想定解とは違うだろうなと思っていたけど、Python における想定解そのものだったとは orz. 平衡二分探索木がビルトインな言語ってどれくらいあるんだろう. 少なくとも Python と Go は入ってないし、Java も無かった気がするよなあ. C++ だけ簡単な問題になっている感. 逆順と Union Find で出来ないかなとも考えてみて思いつかなかったけど、想定解にあるってことは書きようがあるんだなあ. どう書くんだろ. from sys import stdin class SegmentTree: def __init__(self, size, op, e): self._op = op self._e = e self._size = size t = 1 while t < size: t *= 2 self._offset = t - 1 self._data = [e] * (t * 2 - 1) def __getitem__(self, key): return self._data[self._offset + key] def __setitem__(self, key, value): op = self._op data = self._data i = self._offset + key data[i] = value while i >= 1: i = (i - 1) // 2 data[i] = op(data[i * 2 + 1], data[i * 2 + 2]) def build(self, iterable): op = self._op data = self._data data[self._offset:self._offset + self._size] = iterable for i in range(self._offset - 1, -1, -1): data[i] = op(data[i * 2 + 1], data[i * 2 + 2]) def query(self, start, stop): def iter_segments(data, l, r): while l < r: if l & 1 == 0: yield data[l] if r & 1 == 0: yield data[r - 1] l = l // 2 r = (r - 1) // 2 op = self._op it = iter_segments(self._data, start + self._offset, stop + self._offset) result = self._e for v in it: result = op(result, v) return result readline = stdin.readline L, Q = map(int, readline().split()) cx = [tuple(map(int, readline().split())) for _ in range(Q)] a = sorted([x for _, x in cx] + [0, L]) b = {} for i in range(len(a)): b[a[i]] = i st1 = SegmentTree(len(a), max, -1) st2 = SegmentTree(len(a), min, 10 ** 20) st1[b[0]] = b[0] st1[b[L]] = b[L] st2[b[0]] = b[0] st2[b[L]] = b[L] result = [] for c, x in cx: if c == 1: st1[b[x]] = b[x] st2[b[x]] = b[x] elif c == 2: result.append(a[st2.query(b[x], b[L] + 1)] - a[st1.query(0, b[x])]) print(*result, sep='\n')

はじめに データ管理で下図のように、商品No.の検索から別シートにある情報を転記できたらいいのになと思った。 数が少ないうちはコピペなどで入力したらいいけど、データ数が膨大になったらめんどくさいのでpythonで自動化できないか試してみた。 環境(Colabratory) windows 10 Python 3.7.10 pandas 1.1.5 内容 #データの読み込み import pandas as pd import numpy as np from google.colab import files files.upload() df = pd.read_csv('/content/dict.csv',encoding='SHIFT-JIS',index_col=0) df.head() こんな感じのデータ。 #辞書型に変更 d_records = df.to_dict('index') d_records キー(インデックス)と値の紐づけ。 #選ぶインデックスのリスト selects = [3,4,5,4,5,2,10,15] #空のリスト l_data = [] #ループを回して値を取得し、リストに入れていく for i in range(0,len(selects)): select = selects[i] data = d_records[select] l_d = list(data.values()) l_data.append(l_d) #リストをデータフレーム化 files_df = pd.DataFrame(l_data,columns=['code','number']) #選んだインデックスと選ばれたデータを結合して保存 selects = pd.DataFrame(selects) c_df = pd.concat((selects,files_df),axis=1) c_df.to_csv('tenki.csv', encoding='shift_jis') 想定しているものはできた。体裁は適宜自分でそろえる必要がありそう。 まとめ 今回は試しのためデータ数を少なくしているが、多くなってくるとより便利に使えそう。 データは下記にアップ。 https://github.com/dem-kk/File/tree/main/dict 参考 https://techacademy.jp/magazine/15639 https://www.delftstack.com/ja/howto/python/python-convert-dictionary-to-list/

はじめに 昨日、macOSのcronで自動実行するためにmacOSの省電力モードをOFFにすると1ヶ月で1,200円程度電気代がかかるのでは!?という記事を書きました(あくまで試算です) やはり、オンプレ運用はよくないと。これからはクラウドっしょ！ ただ、AWS Lambdaってよくわからん。 って思っている方いますね？(私) ってことで本日はAWS Lambda へPython実行環境を移行にチャレンジしてみました。 少し長くなるので、事前準備編とAWS編で分けて投稿します。 環境 ・既にAWSのアカウントを持っている(お持ちでない方でもすぐ発行できますよ。) ・Python3.9 ・macOS Catalina10.15.7 ・定期実行したいpythonがある 事前解説 今回はサンプルとして、ccxtという複数の仮想通貨取引所のAPI操作が集約されたライブラリをAWS Lambdaで実行したいと思います。 AWS Lambdaにはデフォルトではccxtがインストールされていません。 またmacOSならターミナルで簡単にpip3コマンドでccxtをインストールできますが、残念ながらAWS Lambdaにはインストールできません。なので、みなさんが一生懸命作ったpythonプログラムの1行目のimportコマンドでいきなりエラーになります。 import ccxt　→エラー発生 { "errorMessage": "Unable to import module 'lambda_function': No module named 'ccxt'", "errorType": "Runtime.ImportModuleError" } この時点でAWS Lambda へPython実行環境を移行を諦める人多いと思います。(私) でも今日は諦めません。だって電気代高いんだもん！ ではどうすればいいのか。簡単に書きます。 ◆事前準備編◆ 　1.macOSのデスクトップにtempフォルダを作る 　2.tempフォルダにccxtをpip3でインストールする 　3.既に作ったpythonプログラムを lambda_function.py へリネームしてtempフォルダに入れる 　4.lambda_function.py の先頭に指定プログラムを追記する 　5.ターミナルでzipでまとめる ◆AWS編◆ 　6.AWSへログインし、5をアップロードする 　7.EventBridge (CloudWatch Events)に定期実行ルールを登録する 今回は1-5までを説明します。(次回は6-7) 手順 1.macOSのデスクトップにtempフォルダを作る 　1-1.ターミナルを起動します。pwdコマンドで自分がどこにいるか確認します。 % pwd /Users/Hiroki 　1-2.デスクトップに移動します。 % cd /Users/Hiroki/Desktop Desktop % 　1-3.デスクトップにtempフォルダを作ります。 Desktop % mkdir temp 　1-4.tempフォルダに移動します。 Desktop % cd temp temp % 2.tempフォルダにccxtをpip3でインストールする 　2-1.ccxtをフォルダ指定してインストールします。 temp % pip3 install ccxt -t ./ 　2-2.デスクトップにあるtempフォルダにccxtフォルダがあるか確認します。 3.既に作ったpythonプログラムを lambda_function.pyへリネームしてtempフォルダに入れる 4.lambda_function.py の先頭に指定のプログラムを追記する lambda_function.py def lambda_handler(event, context): return { } おまじないです。 5.ターミナルでzipでまとめる temp % zip -r zip_file ./* tempフォルダの中に　zip_file.zip　というファイルがあれば成功です。 最後に うまく事前準備できましたでしょうか。 次回はAWSへアップロードして、定期実行を設定していきます。 AWSアカウントをお持ちでない方はご準備をお願いします。 ご参考になれば幸いです。

Pythonバージョンの確認 python -V 環境の一覧表示 workon 環境「azure」のアクティブ化 workon azure pip自体をアップグレード pip install --upgrade pip インストール済みのパッケージの一覧表示 pip list 環境「azure」の非アクティブ化 deactivate 環境「azure2」の作成 mkvirtualenv azure2 環境「azure2」の削除 rmvirtualenv azure2 導入 sudo pip3 install virtualenv sudo pip3 install virtualenvwrapper vim ~/.profile .profileに追記 source /usr/local/bin/virtualenvwrapper.sh export WORKON_HOME=~/.virtualenvs

Sphinx拡張で使おうと思っている正規表現のメモ。 思い付く限りのパターンを試して、問題ないことを確認したサンプル。 サンプル #!/usr/bin/python import re regex = re.compile(r'([ 　]*)((.*?)\|)*([^_]*)((_)([0-9a-z]*)?)?') ps = [] ss = [ "読み１|用語１", "読み２|用語２_120a3", "読み３|用語３_", "読み４|用語４_120a3あいうえお", "読み５|用語５_あいうえお", "用語６", "用語７_120a3", "用語８_", "用語９_120a3あいうえお", "用語Ａ_あいうえお", " 読み１|用語１", " 読み２|用語２_120a3", " 読み３|用語３_", " 読み４|用語４_120a3あいうえお", " 読み５|用語５_あいうえお", " 用語６", " 用語７_120a3", " 用語８_", " 用語９_120a3あいうえお", " 用語Ａ_あいうえお", " 　読み１|用語１", " 　読み２|用語２_120a3", " 　読み３|用語３_", " 　読み４|用語４_120a3あいうえお", " 　読み５|用語５_あいうえお", " 　用語６", " 　用語７_120a3", " 　用語８_", " 　用語９_120a3あいうえお", " 　用語Ａ_あいうえお", ] for s in ss: p = regex.match(s) print(p.groups()) 以上

環境 VS Code(ターミナルGit Bash) Python3.9.1 仮想環境作成 アプリ作成するフォルダに移動 env: 任意 python -m venv env アクティブ化 windows source env/Scripts/activate Mac source env/bin/activate VS Codeインタープリンタ選択 コマンドパレットで「Python: Select Interpreter」入力し、 作成した仮想環境のPythonを選択