概要

numpy初心者向けの記事です。Pythonを使い始めて間もない頃、1行列のshapeの出力について勘違いしてました。割と引っかかりやすいんじゃないかなと思い、実例をもとに共有いたします。

例

まず以下のarrayのshapeを見てみる。

a = np.array([[1, 2, 3], [4, 5, 6]])
print(a.shape)

出力は以下のようになります。

(2, 3)

左側に行列の行、右側に列が表示される。数学的にも2行3列と表現されるので納得。

\begin{bmatrix}
1 & 2 & 3 \\
4 & 5 & 6
\end{bmatrix}

行を一つ減らして、再びshapeを出力してみる。

b = np.array([1, 2, 3])
print(b.shape)

出力は以下。

(3, )

...? 1行3列で、(1, 3) が出力されると思ったけど違った。。

以下のように行列積を計算してみる

a = np.array([2, 2])
b = np.array([[1, 2], [3, 4]])
print(np.dot(a, b))
print(np.dot(b, a))

すると、結果は

[ 8 12]
[ 6 14]

行列積の計算の際は、もう一つの行列の shape に応じて、柔軟に行ベクトル、列ベクトルとして計算してくれることがわかる。

まとめ

一次元配列のshapeは要素が1個のタプルになり、(要素数, )として表現される。数学的には1つの行のみを持つ行列を行ベクトル、1つの列のみを持つ行列を列ベクトルと呼ぶが、ndarrayの一次元配列では行ベクトルと列ベクトルの区別はない。

はじめに

前回の記事に、感染症の拡散を説明するSIRモデルについて説明しました。その時に、SIRモデルを解く数値積分のPythonプログラムも掲載しました。今回の記事には、そのPythonプログラムをGUI化してみます。

前回の記事：感染病の数学予測モデル　(SIRモデル）：事例紹介(1) https://qiita.com/kotai2003/items/d74583b588841e6427a2

GUIの説明

自分が作ったプログラムをGUI化する理由は、「他人に使ってもらいたいから」です。将来的には、GUIのPyコードだけではなくて、setup.exeの形で配布したいと思います。

まず、今回に作成したGUIの画面構成です。SIRモデルの感染率（Beta)と除去率(Gamma)を右側の空欄に入力し、その下のDrawボタンを押すと、SIRモデルの計算結果が中央の画面にプロットされます。

\begin{align}

\frac{dS}{dt} &= -\beta SI \\
\frac{dI}{dt} &=  \beta SI -\gamma I \\
\frac{dR}{dt} &=  \gamma I \\

\end{align} 

\begin{align}
S &: 感染可能者 \quad \text{(Susceptible)}　\\
I &: 感染者 \quad \text{(Infectious)}　\\
R &: 感染後死亡者、もしくは免疫を獲得した者 \quad \text{(Removed)}　\\
\beta &: 感染率\quad \text{(The infectious rate)} \quad [1/day] \\
\gamma &:除去率\quad \text{(The Recovery rate)} \quad [1/day]　\\
\end{align}

初期の画面に、Typical Condition~と書いてありますが、これは感染率（Beta)と除去率(Gamma)を入力する際に参考にするためです。

次の画面が計算結果を示します。SIRモデルの数値積分結果がLine Plotで表示されます。新しく入力した感染率（Beta)と除去率(Gamma)の値が、タイトルに更新されます。

GUIの終了は、右下のQuitボタンを押すことで実行されます。

プログラムの説明

GUIのライブラリは、Tkinterを利用しました。
簡単にプログラムの動作を説明します。

(1) TkinterでCanvas Widgetを用意する。 TkinterのEntry Widgetで、感染率（Beta)と除去率(Gamma)を入力するテキストボックスを用意する。
(2) scipyでSIRモデルを計算し、その結果をmatplotlibでプロットする。
(3)　matplotlibのFigureを Canvas WidgetでDrawする。

必要なライブラリをインポートします。

import tkinter
import tkinter.messagebox as tkmsg
import matplotlib.pyplot as plt
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
from matplotlib.backends.backend_tkagg import NavigationToolbar2Tk

import numpy as np
from functools import partial
from scipy.integrate import odeint

SIRモデルの数値積分に必要な関数を用意します。

# SIR Differential Equation
def SIR_EQ(v, t, beta, gamma):
    '''
    dS/dt = -beta * S * I
    dI/dt = beta * S * I - gamma * I
    dR/dt = gamma * I

    [v[0], v[1], v[2]]=[S, I, R]

    dv[0]/dt = -beta * v[0] * v[1]
    dv[1]/dt = beta * v[0] * v[1] - gamma * v[1]
    dv[2]/dt = gamma * v[1]

    '''
    return [-beta*v[0]*v[1], beta * v[0] * v[1] - gamma * v[1], gamma * v[1]]

# Solving SIR Equation
def Calc_SIR(var_b, var_r):
    # parameters
    t_max = 14
    dt = 0.01
    #beta_const = 0.0026
    #gamma_const = 0.5

    beta_const = var_b
    gamma_const = var_r

    # initial_state
    S_0 = 762
    I_0 = 1
    R_0 = 0
    ini_state = [S_0, I_0, R_0]  # [S[0], I[0], R[0]]

    # numerical integration
    times = np.arange(1, t_max, dt)
    args = (beta_const, gamma_const)
    # Solver SIR model
    result = odeint(SIR_EQ, ini_state, times, args)
    return times,result

TkinterのボタンのEventの時に実行される関数を用意します。
DrawCanvas関数に、matplotlibのaxesインスタンス(=ax)を用意し、そのインスタンスをプロットする仕組みです。最後のcanvas.draw()で、matplotlibのaxesインスタンスが、canvas widgetに表示されることになります。

def Quit():
    tkmsg.showinfo("Tomomi Research Inc.","Thank you for running this program!")
    root.quit()
    root.destroy()

#Draw Button
def DrawCanvas(canvas, ax):
    value_beta = EditBox_beta.get()
    value_gamma = EditBox_gamma.get()
    if value_beta != '':
        EditBox_beta.delete(0, tkinter.END)
        EditBox_gamma.delete(0, tkinter.END)
        ax.cla()
        beta = float(value_beta)
        gamma = float(value_gamma)
        t_r, res_r = Calc_SIR(beta,gamma)

        ax.plot(t_r, res_r)
        ax.legend(['Susceptible', 'Infectious', 'Removed'])
        ax.set_title('Beta='+str(beta)+', Gamma='+str(gamma) )
        ax.set_xlabel('time(days)')
        ax.set_ylabel('population')
    canvas.draw()

この部分はメインのプログラムです。
Tkinterの各Widgetの設定は、この部分に書いてあります。

if __name__ == '__main__':
    try:
        # GUI generate
        root = tkinter.Tk()
        root.title("SIR model")

        # Graph setting
        fig, ax1 = plt.subplots()
        #fig.gca().set_aspect('equal', adjustable='box')  # グラフ領域の調整 #get current axes

        ax1.set_title('Typical Condition:  beta=0.0026, gamma=0.5, $S_0$=762, $I_0$=1')


        # Generate Canvas
        Canvas = FigureCanvasTkAgg(fig, master=root)
        Canvas.get_tk_widget().grid(row=0, column=0, rowspan=10)

        #Beta
        EditBox_beta = tkinter.Entry(width=5)  # テキストボックスの生成
        EditBox_beta.grid(row=1, column=2)

        GridLabel_beta = tkinter.Label(text="Beta")
        GridLabel_beta.grid(row=1, column=1)

        # Gamma
        EditBox_gamma = tkinter.Entry(width=5)  # テキストボックスの生成
        EditBox_gamma.grid(row=4, column=2)

        GridLabel_gamma = tkinter.Label(text="Gamma")
        GridLabel_gamma.grid(row=4, column=1)

        # ボタンに関する諸々の設定
        ReDrawButton = tkinter.Button(text="Draw", width=15, command=partial(DrawCanvas, Canvas, ax1))  # ボタンの生成
        ReDrawButton.grid(row=5, column=1, columnspan=2)  # 描画位置(テキトー)

        QuitButton = tkinter.Button(text="Quit", width=15, command=Quit)  # ボタンの生成
        QuitButton.grid(row=7, column=1, columnspan=2)  # 描画位置(テキトー)

        DrawCanvas(Canvas, ax1)
        root.mainloop()
    except:
        import traceback

        traceback.print_exc()
    finally:
        input(">>")  # エラー吐き出したときの表示待ち

まとめ

Python GUIのライブラリには、Tkinterと PyQTが有名だそうですが、Tkiknterが初心者にもプログラムしやすいと聞いたため、Tkinterを選んでみました。今だにもPyQTに興味がありますので、両方使ったことがある方は、コメントを残していただければと思います。

PythonのGUIのメリットは、Cross Platformであることだと思います。プログラムを一つ作成すれば、変更なしにWindows, Mac OS, Linuxに使えることが可能です。現在、WindowsでDeep Learningのプログラムを開発し、Nvidia製のJetson系(Linux）に実装することを考えています。その時、Python GUIが大変役に立つと思います。

参考資料

1. MatplotlibにGUI(Tkinter)を組み合わせる話 https://www.shtsno24.tokyo/2018/12/matplotlibguitkinter.html
2. Matplotlib&Seaborn実装ハンドブック https://amzn.to/2ujQ2CL
3. 感染病の数学予測モデルの紹介　(SIRモデル）https://qiita.com/kotai2003/items/3078f4095c3e94e5325c
4. 感染病の数学予測モデル　(SIRモデル）：事例紹介(1) https://qiita.com/kotai2003/items/d74583b588841e6427a2

AtCoder Beginner Contest 156 参戦記

ABC156A - Beginner

4分で突破. 書くだけ.

N, R = map(int, input().split())

if N >= 10:
    print(R)
else:
    print(R - 100 * (10 - N))

ABC156B - Digits

2分で突破. 書くだけ. K進数の 10 は K で、100 は K² だと分かっていれば、K で何回割れるか調べるだけと分かるはず.

N, K = map(int, input().split())

result = 0
while N != 0:
    result += 1
    N //= K
print(result)

ABC156C - Rally

3分半で突破. 書くだけ……なのは、この計算をしなれてるからなんですが. (X_i−P)² が最小になる P は重心. 減色の K-means で3次元のこの計算をしてます…….

N = int(input())
X = list(map(int, input().split()))

P = int(sum(X) / N + 0.5)
print(sum((x - P) * (x - P) for x in X))

ABC156D - Bouquet

70分半で突破. まず _nC₀+_nC₁+...+_nC_n = 2ⁿ を知っていないといけない. 私は知りませんでした orz. パスカルの三角形の Wikipedia 見てたら書いてあって知りました. これで過去に書いた mpow と mcomb を貼って終わりかと思ったら n! を求める部分があって 2≤n≤10⁹ に殺される. a, b はたかだか 10⁵ であることに気づき、_nC_k を求める別の式に組み替えてようやく解けた. 良問!

def mpow(x, n):
    result = 1
    while n != 0:
        if n & 1 == 1:
            result *= x
            result %= 1000000007
        x *= x
        x %= 1000000007
        n >>= 1
    return result


def mcomb(n, k):
    a = 1
    b = 1
    for i in range(k):
        a *= n - i
        a %= 1000000007
        b *= i + 1
        b %= 1000000007
    return a * mpow(b, 1000000005) % 1000000007


n, a, b = map(int, input().split())

result = mpow(2, n) - 1
result -= mcomb(n, a)
result + 1000000007
result %= 1000000007
result -= mcomb(n, b)
result + 1000000007
result %= 1000000007
print(result)

はじめに

基本的な機械学習の手順：②データを準備しようでは、BigQueryで作ったテーブルをPytohn環境にPandas Dataframe形式で取り込む処理をしてきました。

ただ、テーブルのサイズが大きくなると、結構時間がかかってしまいます。
たぶん、そんな悩みを持つ方が多かったのでしょう。そこで出てきたのが、BigQuery Storage APIという新たなサービス。

一説には７～８倍も速いと聞きましたが、どうなのでしょう。試してみたいと思います。

分析環境

Google BigQuery
Google Colaboratory

参考にしたサイト

BigQuery Storage API を使用して BigQuery データを pandas にダウンロードする

対象とするデータ

用いるテーブルは、myproject.mydataset.mytableという約100MBのテーブルです。
それを、下記の様に全件取ってくるだけというシンプルな処理で、Pandas Dataframe形式で取り込みます。

query="SELECT * FROM `myproject.mydataset.mytable`

1.BigQuery 標準API

まずは、これまでも使っているBigQueryの標準APIでやってみます。

import time
from google.cloud import bigquery
start = time.time()

client = bigquery.Client(project="myproject")
df = client.query(query).to_dataframe()

elapsed_time = time.time() - start

約120秒で処理できました。まあ、これくらいなら許容範囲ですが。

2.Pandas read_gbq

BigQueryのAPIを使わなくても、Pandasの機能でできるよね。ということで、そちらも試してみます。

import time
import pandas as pd

start = time.time()

df = pd.io.gbq.read_gbq(query, project_id="myproject", dialect="standard")

elapsed_time = time.time() - start

約135秒で処理が完了。BigQueryのAPIより少し遅くなっていますね。
BigQueryの標準APIでも、Pandasの機能に比べると、何か工夫がされているようです。

3.BigQuery Storage API

そして、今回のテーマであるBigQuery Storage APIの出番です。
Colabでライブラリーをimportしようとしたら、ライブラリーが無いよと言われてしまったので、まずはinstallから。

pip install --upgrade google-cloud-bigquery-storage

と、インストールするとランタイムの再起動を求めるメッセージが。時々他のライブラリーでも出ますが、ちょっと面倒ですね。

WARNING: The following packages were previously imported in this runtime:
  [google]
You must restart the runtime in order to use newly installed versions.

さて、ランタイムを再起動して、改めてライブラリーをimportして実行します。

import time
from google.cloud import bigquery
from google.cloud import bigquery_storage
start = time.time()

client = bigquery.Client(project="myproject")
bqstorageclient = bigquery_storage.BigQueryStorageClient()
df3 = (
    client.query(query)
    .result()
    .to_dataframe(bqstorage_client=bqstorageclient)
)

elapsed_time = time.time() - start

実行時間は、なんと驚異の約12秒。標準APIの7～8倍どころか、10倍出てます。
偶然かと思って、何度かやってみましたが、1～2秒程度の誤差はあるものの、ほぼこのスピードで完了しました。

おわりに

予想よりもかなり速い結果が返ってきて、ビックリしました。
普段の10倍速ければ、数GBとかのデータも短時間で取り込むことが可能ですね。（その後のPythonでの処理が重そうですが）

普通にBigQueryを回すのに加えて、$1.10 per TBと費用が掛かるので、乱発はできませんが、テーブルが大きすぎてデータ取り込みに何十分も待たないといけない、というときには使っていきたいサービスですね。

はじめに

これを作ろうと思ったきっかけについて。
Raspberry Piでちょっとしたバッチ処理を動かすことになったのですが、Raspberry Piに限らずサーバーを長期間稼働させていると、知らぬ間にサーバー停止していた。という事象があり得ます。
また、サーバー自体は動いているけど、何かしらの理由でバッチ処理がコケ始めるということも考えられますが、流石に毎日ログインして確認するのは手間になります。
なので、Raspberry Pi上でバッチ処理が正常か異常かの結果をスマホに通知できると嬉しいと思いました。通知がない場合はサーバーで異常が起きたと判断できますし。
それを実装する一つの方法として、
1. cronでGoogle Spreadsheetにデータ更新を行うpythonプログラムを定期実行させる
2. IFTTTでGoogle Spreadsheetの更新をトリガーにLINE通知させる
上記で実現できそうだったのでやってみた。という記事でございます。
※色々手探りでやってみた感じなので、拙い感じとなってますが。。。

(1) Google SpreadSheetへの投稿準備

• Python Quickstartにアクセスして、以下を実施
• "Enable the Google Sheets API"ボタンをクリックして、適当なプロジェクト名を入力。「credentials.json」をローカルに保存する
• 下記コマンドで専用のライブラリを導入する。

pip install --upgrade google-api-python-client google-auth-httplib2 google-auth-oauthlib

投稿準備について、過去の記事で分かりやすくまとまっているのでオススメです。
https://qiita.com/connvoi_tyou/items/7cd7ffd5a98f61855f5c

(2) Google SpreadSheetへの投稿用プログラムの作成

あらかじめ、Google Drive上にSpreadsheetを作成しておきます。
ファイルURLの一部をプログラム内に記述する必要があるので、控えておきます。

https://docs.google.com/spreadsheets/d/<XXXXXXXXXXXXXXXXXXX>/edit#gid=zzzzzzzz
# "<XXXXXXXXXXXXXXXXXXX>"を控えておく

下記のような、GoogleSpreadsheetへの
「server_check.py」という名前で作っておきました。

server_check.py

from __future__ import print_function
import pickle
import os.path
from googleapiclient.discovery import build
from google_auth_oauthlib.flow import InstalledAppFlow
from google.auth.transport.requests import Request

SCOPES = 'https://www.googleapis.com/auth/spreadsheets'
creds = None
if os.path.exists('token.pickle'):
    with open('token.pickle', 'rb') as token:
        creds = pickle.load(token)
# If there are no (valid) credentials available, let the user log in.
if not creds or not creds.valid:
    if creds and creds.expired and creds.refresh_token:
        creds.refresh(Request())
    else:
        flow = InstalledAppFlow.from_client_secrets_file('credentials.json', SCOPES)
        creds = flow.run_local_server(port=0)
    # Save the credentials for the next run
    with open('token.pickle', 'wb') as token:
        pickle.dump(creds, token)

service = build('sheets', 'v4', credentials=creds)

# ここまではgoogle提供のガイドに記載されているものとほぼ同じ
# 以下、Spreadsheet固有の情報を記入。
spreadsheet_id = '<XXXXXXXXXXXXXXXXXXX>'
sheetname='シート1'
range_ = sheetname

# 日付情報を取得する
import datetime
_str_dt=datetime.datetime.now().strftime("%Y/%m/%d %H:%M:%S")

# A列、B列、C列へ書き込むため、データを代入
# A列は日付を、B列、C列には任意の文字列を(今回はテスト的な文字列で代用)
values = [
  [_str_dt,"Status-1","Status-2"]
]
body = {
  'values' : values
}

# Spreadsheetへの書き込みを実施
result=service.spreadsheets().values().append(spreadsheetId=spreadsheet_id,valueInputOption='RAW',range=range_,body=body).execute()
print(result)

上記プログラムでは"Status-1"のようなテストメッセージにしています。
実際に運用させるときは、ログファイルの出力を抽出したりして、
values配列に追加する形になると思います。

(3) プログラムの実行とブラウザでのアクセス許可

必要な箇所が修正できたら、上記pythonプログラムを実行します。
その際、「credentials.json」ファイルが同じディレクトリに存在している必要があります。
また、プログラムの初回実行時はブラウザ認証する必要があるため、GUI環境で実行する必要があります。

実行すると、ブラウザが立ち上がりgoogleアカウントの確認ページが表示されます。画面に従ってアクセス権を付与します。
(途中、安全でないページへの警告が表示されますが続けます)
"The authentication flow has completed. You may close this window."みたいな文字がブラウザに出力されると設定完了ぽいです。

アクセス権の付与が完了すると、Spreadsheetにエントリーが追加(先頭行に値が入力)されているはずです。
また、プログラムと同じディレクトリに「token.pickle」ファイルが作成されます。
このファイルがあれば、次回以降はブラウザでの認証をSKIPして、行が追加されるようになります。

(4) IFTTTでの通知設定

1. IFTTTアプリが必要なので、ない場合はインストール＆初期セットアップ
2. "Make your own Applets from scratch"でアプレットを作る
3. This(トリガー)は「Google Sheets」->「New row addedto spreadsheet」
4. (2)で作成したシート情報を入力する
5. That(アクション)は「LINE」->「Send message」
6. Recipientは"1:1でLINE Notifyから通知を受け取る"
7. メッセージ本文をカスタマイズする

はまったところ

プログラム実行環境でブラウザが開けない場合

今回の最終的な目標はラズベリーパイから対象プログラムを動かすことでしたが、ラズベリーパイにはGUI環境は導入していませんでした（通常はSSH接続で作業）。
ラズベリーパイ上で"server_check.py"を実行してもブラウザが開けないため、アクセス許可の設定ができません。
＊コンソール上に表示されるURLを、PC上のブラウザにコピペすると最後の画面でエラーになってしまいます。

そういった場合ですが、別のGUI環境で同じプログラムを実行して、生成された「token.pickle」ファイルを対象のサーバーに配置すれば、CUI上でもプログラムを実行できると思われます。
自分は同じプログラムをMac上で実行して、生成されたtoken.pickleをラズベリーパイに送ることで、ラズベリーパイからも実行できるようになりました。

プログラム実行時になんかエラーが出てしまう

たまに、プログラムを実行するとエラーが出てしまうことがありました。

[root@localhost python]# python server_check.py 
Traceback (most recent call last):
  File "server_check.py", line 17, in <module>
    creds.refresh(Request())
  File "/usr/local/lib/python3.7/site-packages/google/oauth2/credentials.py", line 182, in refresh
    self._scopes,
  File "/usr/local/lib/python3.7/site-packages/google/oauth2/_client.py", line 248, in refresh_grant
    response_data = _token_endpoint_request(request, token_uri, body)
  File "/usr/local/lib/python3.7/site-packages/google/oauth2/_client.py", line 124, in _token_endpoint_request
    _handle_error_response(response_body)
  File "/usr/local/lib/python3.7/site-packages/google/oauth2/_client.py", line 60, in _handle_error_response
    raise exceptions.RefreshError(error_details, response_body)
google.auth.exceptions.RefreshError: ('invalid_scope: Some requested scopes were invalid. {invalid=[a, c, d, e, g, h, i, l, m, ., /, o, p, r, s, t, u, w, :]}', '{\n  "error": "invalid_scope",\n  "error_description": "Some requested scopes were invalid. {invalid\\u003d[a, c, d, e, g, h, i, l, m, ., /, o, p, r, s, t, u, w, :]}",\n  "error_uri": "http://code.google.com/apis/accounts/docs/OAuth2.html"\n}')
[root@localhost python]# 

どういったエラーなのかはっきりと分からないのですが、自分は「token.pickle」を再作成すると解消しました。
再作成するには、「token.pickle」を削除してプログラムを実行すれば良いです。
ただ、またブラウザからの設定からになるので、ブラウザを開けないCUI環境は注意してください。

ちなみに、どういった理由で起きるのかも不明なのですが、個人的な感覚だと下記のタイミングでエラーになった気がします。

• プログラムを大幅に改修した際
• 別のプログラムを実行して認証を行なった際

おわりに

これでバッチ処理の異常時や、サーバー停止時に早めに気づけるようになりましたが、毎日LINEがくるとちょっとうんざりですね。
本当は「便りが無いのは元気な証拠」スタイルで、異常を検知した際(サーバーから異常を検知、もしくは定時連絡がない)に、LINEへ通知させるようにしたかったのですが、それは次の課題ということで。

英語の文章を日本語に翻訳するスクリプトです。
以前のコードをクラスにしたものです。配布されているライブラリがうまく動かなかったので作成したものです。

GoogleCloudのAPI Keyは取得していただく必要があります。（すぐできます）

google_translate.py

import requests
import json
import time

class GoogleTranslate:
    def __init__(self):
        self.private_key = '<ここにAPIKeyを入れる>'

    def post_text(self,text):
        url_items = 'https://www.googleapis.com/language/translate/v2'
        item_data = {
            'target': 'ja',
            'source': 'en',
            'q':text
        }

        response = requests.post('https://www.googleapis.com/language/translate/v2?key={}'.format(self.private_key), data=item_data)
        return json.loads(response.text)["data"]["translations"][0]["translatedText"]


    def split_and_send_to_post(self,text):
        sen_list = text.split('.')

        to_google_sen = ""
        translated_text = ""

        for index, sen in enumerate(sen_list[:-1]):
            to_google_sen += sen + '. '
            if len(to_google_sen)>1000:
                #1000字を超えてたらgoogleに送信
                translated_text += self.post_text(to_google_sen)
                time.sleep(3)

                to_google_sen = ""

            if index == len(sen_list)-2:
                #最後の文章の翻訳
                translated_text += self.post_text(to_google_sen)
                time.sleep(3)

        return translated_text


    def main(self,text):
        original_text = text
        if original_text[-1] != '.':
            original_text+='.'
            #.で終わっていない場合に分割処理に支障が出るので
        text_translated = self.split_and_send_to_post(original_text)
        print(text_translated)
        return text_translated

test.py

from google_translate import *

input_text = input('英語に変換して検索したいワードを日本語で入力してください:')
google_tr = GoogleTranslate()
translated = google_tr.main(input_text)

print(translated)

はじめに

ふと思い立って勉強を始めた「ゼロから作るDeep LearningーーPythonで学ぶディープラーニングの理論と実装」の8章で私がつまずいたことのメモです。

実行環境はmacOS Mojave + Anaconda 2019.10、Pythonのバージョンは3.7.4です。詳細はこのメモの1章をご参照ください。

（このメモの他の章へ：1章 / 2章 / 3章 / 4章 / 5章 / 6章 / 7章 / 8章）

8章 ディープラーニング

この章は、層を深くしたディープなニューラルネットワークの説明です。

8.1 ネットワークをより深く

これまでに学んだことを使って、ディープなネットワークでMNISTの手書き文字認識の実装に挑戦します。残念ながらこの章はソースの解説が全くないので大変です。

前章までで学んだDropoutやAdamは実装をサボっていたのですが、今回使うのでここから片付けます。

(1)Dropoutレイヤーの実装

Dropoutレイヤーは本の「6.4.3 Droput」に実装の解説があるので、それを見ながら実装しました。

dropout.py

# coding: utf-8
import numpy as np


class Dropout:
    def __init__(self, dropout_ratio=0.5):
        """Dropoutレイヤー

        Args:
            dropout_ratio (float): 学習時のニューロンの消去割合、デフォルトは0.5。
        """
        self.dropout_ratio = dropout_ratio              # 学習時のニューロンの消去割合
        self.valid_ratio = 1.0 - self.dropout_ratio     # 学習時に生かしていた割合
        self.mask = None                                # 各ニューロンの消去有無を示すフラグの配列

    def forward(self, x, train_flg=True):
        """順伝播

        Args:
            x (numpy.ndarray): 入力
            train_flg (bool, optional): 学習中ならTrue、デフォルトはTrue。

        Returns:
            numpy.ndarray: 出力
        """
        if train_flg:
            # 学習時は消去するニューロンを決めるマスクを生成
            self.mask = np.random.rand(*x.shape) > self.dropout_ratio

            # 出力を算出
            return x * self.mask

        else:
            # 認識時はニューロンは消去しないが、学習時の消去割合を加味した出力に調整する
            return x * self.valid_ratio

    def backward(self, dout):
        """逆伝播

        Args:
            dout (numpy.ndarray): 右の層から伝わってくる微分値

        Returns:
            numpy.ndarray: 微分値（勾配）
        """
        # 消去しなかったニューロンのみ右の層の微分値を逆伝播
        assert self.mask is not None, '順伝播なしに逆伝播が呼ばれた'
        return dout * self.mask

(2)Adamの実装

最適化に使うAdamは、本の「6.1.6 Adam」に簡単な解説があるのですが、簡単すぎてこれだけでは実装できません。また、本のソースを見てもアルゴリズムが良く分かりませんでした。そこでまず、 @omiita さんの 【2020決定版】スーパーわかりやすい最適化アルゴリズム -損失関数からAdamとニュートン法- で大まかな仕組みを理解しました。そして、本で紹介されている 原著論文のPDF （本の参考文献[8]のサイト Adam: A Method for Stochastic Optimization の右上からダウンロードできます）のP.2の「Algorithm 1」の説明を見ながら実装しました。英語ですが擬似コードによる20行程度の説明なので、英語が苦手な私でもなんとかなりました。パラメータの初期値も、この論文の推奨値通りにしてみました。

adam.py

# coding: utf-8
import numpy as np


class Adam:

    def __init__(self, alpha=0.001, beta1=0.9, beta2=0.999):
        """Adamによるパラメーターの最適化

        Args:
            alpha (float, optional): 学習係数、デフォルトは0.001。
            beta1 (float, optional): Momentumにおける速度の過去と今の按分の係数、デフォルトは0.9。
            beta2 (float, optional): AdaGradにおける学習係数の過去と今の按分の係数、デフォルトは0.999。
        """
        self.alpha = alpha
        self.beta1 = beta1
        self.beta2 = beta2

        self.m = None   # Momentumにおける速度
        self.v = None   # AdaGradにおける学習係数
        self.t = 0      # タイムステップ

    def update(self, params, grads):
        """パラメーター更新

        Args:
            params (dict): 更新対象のパラメーターの辞書、keyは'W1'、'b1'など。
            grads (dict): paramsに対応する勾配の辞書
        """
        # mとvの初期化
        if self.m is None:
            self.m = {}
            self.v = {}
            for key, val in params.items():
                self.m[key] = np.zeros_like(val)
                self.v[key] = np.zeros_like(val)

        # 更新
        self.t += 1     # タイムステップ加算
        for key in params.keys():

            # mの更新、Momentumにおける速度の更新に相当
            # 過去と今の勾配を beta1 : 1 - beta1 で按分する
            self.m[key] = \
                self.beta1 * self.m[key] + (1 - self.beta1) * grads[key]

            # vの更新、AdaGradにおける学習係数の更新に相当
            # 過去と今の勾配を beta2 : 1 - beta2 で按分する
            self.v[key] = \
                self.beta2 * self.v[key] + (1 - self.beta2) * (grads[key] ** 2)

            # パラメーター更新のためのmとvの補正値算出
            hat_m = self.m[key] / (1.0 - self.beta1 ** self.t)
            hat_v = self.v[key] / (1.0 - self.beta2 ** self.t)

            # パラメーター更新、最後の1e-7は0除算回避
            params[key] -= self.alpha * hat_m / (np.sqrt(hat_v) + 1e-7)

(3)畳み込み層とプーリング層の出力サイズの計算

今回は層が多く、畳み込み層とプーリング層の出力サイズの計算が何度も出てきます。そのため、それぞれconv_output_sizeとpool_output_sizeという関数としてfunctions.pyに追加しました。他の関数は前章までのままです。

functions.py

# coding: utf-8
import numpy as np


def softmax(x):
    """ソフトマックス関数

    Args:
        x (numpy.ndarray): 入力

    Returns:
        numpy.ndarray: 出力
    """
    # バッチ処理の場合xは(バッチの数, 10)の2次元配列になる。
    # この場合、ブロードキャストを使ってうまく画像ごとに計算する必要がある。
    # ここでは1次元でも2次元でも共通化できるようnp.max()やnp.sum()はaxis=-1で算出し、
    # そのままブロードキャストできるようkeepdims=Trueで次元を維持する。
    c = np.max(x, axis=-1, keepdims=True)
    exp_a = np.exp(x - c)  # オーバーフロー対策
    sum_exp_a = np.sum(exp_a, axis=-1, keepdims=True)
    y = exp_a / sum_exp_a
    return y


def cross_entropy_error(y, t):
    """交差エントロピー誤差の算出

    Args:
        y (numpy.ndarray): ニューラルネットワークの出力
        t (numpy.ndarray): 正解のラベル

    Returns:
        float: 交差エントロピー誤差
    """

    # データ1つ場合は形状を整形（1データ1行にする）
    if y.ndim == 1:
        t = t.reshape(1, t.size)
        y = y.reshape(1, y.size)

    # 誤差を算出してバッチ数で正規化
    batch_size = y.shape[0]
    return -np.sum(t * np.log(y + 1e-7)) / batch_size


def conv_output_size(input_size, filter_size, pad, stride):
    """畳み込み層の出力サイズ算出

    Args:
        input_size (int): 入力の1辺のサイズ（縦横は同値の前提）
        filter_size (int): フィルターの1辺のサイズ（縦横は同値の前提）
        pad (int): パディングのサイズ（縦横は同値の前提）
        stride (int): ストライド幅（縦横は同値の前提）

    Returns:
        int: 出力の1辺のサイズ
    """
    assert (input_size + 2 * pad - filter_size) \
        % stride == 0, '畳み込み層の出力サイズが割り切れない！'
    return int((input_size + 2 * pad - filter_size) / stride + 1)


def pool_output_size(input_size, pool_size, stride):
    """プーリング層の出力サイズ算出

    Args:
        input_size (int): 入力の1辺のサイズ（縦横は同値の前提）
        pool_size (int): プーリングのウインドウサイズ（縦横は同値の前提）
        stride (int): ストライド幅（縦横は同値の前提）

    Returns:
        int: 出力の1辺のサイズ
    """
    assert (input_size - pool_size) % stride == 0, 'プーリング層の出力サイズが割り切れない！'
    return int((input_size - pool_size) / stride + 1)

(4)ディープなCNNの実装

これで必要なパーツの実装が終わったので、いよいよネットワークの実装です。

まず、今回のネットワークにおける入出力の整理です。

レイヤー	入出力の形状	実装時の形状
	$ (バッチサイズN, チャンネル数CH, 画像の高さH, 幅W) $	$ (100, 1, 28, 28) $
[1] Convolution #1	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 16, 28, 28) $
[2] ReLU #1	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 16, 28, 28) $
[3] Convolution #2	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 16, 28, 28) $
[4] ReLU #2	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 16, 28, 28) $
[5] Pooling #1	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 16, 14, 14) $
[6] Convolution #3	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 32, 14, 14) $
[7] ReLU #3	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 32, 14, 14) $
[8] Convolution #4	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 32, 16, 16) $
[9] ReLU #4	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 32, 16, 16) $
[10] Pooling #2	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 32, 8, 8) $
[11] Convolution #5	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 64, 8, 8) $
[12] ReLU #5	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 64, 8, 8) $
[13] Convolution #6	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 64, 8, 8) $
[14] ReLU #6	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 64, 8, 8) $
[15] Pooling #3	↓
	$ (バッチサイズN, フィルター数FN, 出力の高さOH, 幅OW) $	$ (100, 64, 4, 4) $
[16] Affine #1	↓
	$ (バッチサイズN, 隠れ層のサイズ) $	$ (100, 50) $
[17] ReLU #7	↓
	$ (バッチサイズN, 隠れ層のサイズ) $	$ (100, 50) $
[18] Dropout #1	↓
	$ (バッチサイズN, 隠れ層のサイズ) $	$ (100, 50) $
[19] Affine #2	↓
	$ (バッチサイズN, 隠れ層のサイズ) $	$ (100, 10) $
[20] Dropout #2	↓
	$ (バッチサイズN, 隠れ層のサイズ) $	$ (100, 10) $
[21] Softmax	↓
	$ (バッチサイズN, 最終出力サイズ) $	$ (100, 10) $

壮大な表になりましたが、これを１層ずつ実装します。

本のコードはループを使ってシンプルにまとめられていますが、各層の入出力サイズの計算が混乱しそうなので、私は1層ずつパラメーターの初期化とレイヤーの生成を実装しました。かなり泥臭いコードになっています。なお、パラメーターは「Heの初期値」で初期化しています。

deep_conv_net.py

# coding: utf-8
import numpy as np
from affine import Affine
from convolution import Convolution
from dropout import Dropout
from functions import conv_output_size, pool_output_size
from pooling import Pooling
from relu import ReLU
from softmax_with_loss import SoftmaxWithLoss


class DeepConvNet:

    def __init__(
        self, input_dim=(1, 28, 28),
        conv_param_1={
            'filter_num': 16, 'filter_size': 3, 'pad': 1, 'stride': 1
        },
        conv_param_2={
            'filter_num': 16, 'filter_size': 3, 'pad': 1, 'stride': 1
        },
        conv_param_3={
            'filter_num': 32, 'filter_size': 3, 'pad': 1, 'stride': 1
        },
        conv_param_4={
            'filter_num': 32, 'filter_size': 3, 'pad': 2, 'stride': 1
        },
        conv_param_5={
            'filter_num': 64, 'filter_size': 3, 'pad': 1, 'stride': 1
        },
        conv_param_6={
            'filter_num': 64, 'filter_size': 3, 'pad': 1, 'stride': 1
        },
        hidden_size=50, output_size=10
    ):
        """ディープな畳み込みニューラルネットワーク

        Args:
            input_dim (tuple, optional): 入力データの形状、デフォルトは(1, 28, 28)。
            conv_param_1 (dict, optional): 畳み込み層1のハイパーパラメーター、
                デフォルトは{'filter_num':16, 'filter_size':3, 'pad':1, 'stride':1}。
            conv_param_2 (dict, optional): 畳み込み層2のハイパーパラメーター、
                デフォルトは{'filter_num':16, 'filter_size':3, 'pad':1, 'stride':1}。
            conv_param_3 (dict, optional): 畳み込み層3のハイパーパラメーター、
                デフォルトは{'filter_num':32, 'filter_size':3, 'pad':1, 'stride':1}。
            conv_param_4 (dict, optional): 畳み込み層4のハイパーパラメーター、
                デフォルトは{'filter_num':32, 'filter_size':3, 'pad':2, 'stride':1}。
            conv_param_5 (dict, optional): 畳み込み層5のハイパーパラメーター、
                デフォルトは{'filter_num':64, 'filter_size':3, 'pad':1, 'stride':1}。
            conv_param_6 (dict, optional): 畳み込み層6のハイパーパラメーター、
                デフォルトは{'filter_num':64, 'filter_size':3, 'pad':1, 'stride':1}。
            hidden_size (int, optional): 隠れ層のニューロンの数、デフォルトは50。
            output_size (int, optional): 出力層のニューロンの数、デフォルトは10。
        """
        assert input_dim[1] == input_dim[2], '入力データは高さと幅が同じ前提！'

        # パラメーターの初期化とレイヤー生成
        self.params = {}    # パラメーター
        self.layers = {}    # レイヤー（Python 3.7からは辞書の格納順が保持されるので、OrderedDictは不要）

        # 入力サイズ
        channel_num = input_dim[0]                          # 入力のチャンネル数
        input_size = input_dim[1]                           # 入力サイズ

        # [1] 畳み込み層#1 : パラメーター初期化とレイヤー生成
        filter_num, filter_size, pad, stride = list(conv_param_1.values())
        pre_node_num = channel_num * (filter_size ** 2)     # 1ノードに対する前層の接続ノード数
        key_w, key_b = 'W1', 'b1'                           # 辞書格納時のkey
        self.params[key_w] = np.random.normal(
            scale=np.sqrt(2.0 / pre_node_num),              # Heの初期値の標準偏差
            size=(filter_num, channel_num, filter_size, filter_size)
        )
        self.params[key_b] = np.zeros(filter_num)

        self.layers['Conv1'] = Convolution(
            self.params[key_w], self.params[key_b], stride, pad
        )

        # 次の層の入力サイズ算出
        channel_num = filter_num
        input_size = conv_output_size(input_size, filter_size, pad, stride)

        # [2] ReLU層#1 : レイヤー生成
        self.layers['ReLU1'] = ReLU()

        # [3] 畳み込み層#2 : パラメーター初期化とレイヤー生成
        filter_num, filter_size, pad, stride = list(conv_param_2.values())
        pre_node_num = channel_num * (filter_size ** 2)     # 1ノードに対する前層の接続ノード数
        key_w, key_b = 'W2', 'b2'                           # 辞書格納時のkey
        self.params[key_w] = np.random.normal(
            scale=np.sqrt(2.0 / pre_node_num),              # Heの初期値の標準偏差
            size=(filter_num, channel_num, filter_size, filter_size)
        )
        self.params[key_b] = np.zeros(filter_num)

        self.layers['Conv2'] = Convolution(
            self.params[key_w], self.params[key_b], stride, pad
        )

        # 次の層の入力サイズ算出
        channel_num = filter_num
        input_size = conv_output_size(input_size, filter_size, pad, stride)

        # [4] ReLU層#2 : レイヤー生成
        self.layers['ReLU2'] = ReLU()

        # [5] プーリング層#1 : レイヤー生成
        self.layers['Pool1'] = Pooling(pool_h=2, pool_w=2, stride=2)

        # 次の層の入力サイズ算出
        input_size = pool_output_size(input_size, pool_size=2, stride=2)

        # [6] 畳み込み層#3 : パラメーター初期化とレイヤー生成
        filter_num, filter_size, pad, stride = list(conv_param_3.values())
        pre_node_num = channel_num * (filter_size ** 2)     # 1ノードに対する前層の接続ノード数
        key_w, key_b = 'W3', 'b3'                           # 辞書格納時のkey
        self.params[key_w] = np.random.normal(
            scale=np.sqrt(2.0 / pre_node_num),              # Heの初期値の標準偏差
            size=(filter_num, channel_num, filter_size, filter_size)
        )
        self.params[key_b] = np.zeros(filter_num)

        self.layers['Conv3'] = Convolution(
            self.params[key_w], self.params[key_b], stride, pad
        )

        # 次の層の入力サイズ算出
        channel_num = filter_num
        input_size = conv_output_size(input_size, filter_size, pad, stride)

        # [7] ReLU層#3 : レイヤー生成
        self.layers['ReLU3'] = ReLU()

        # [8] 畳み込み層#4 : パラメーター初期化とレイヤー生成
        filter_num, filter_size, pad, stride = list(conv_param_4.values())
        pre_node_num = channel_num * (filter_size ** 2)     # 1ノードに対する前層の接続ノード数
        key_w, key_b = 'W4', 'b4'                           # 辞書格納時のkey
        self.params[key_w] = np.random.normal(
            scale=np.sqrt(2.0 / pre_node_num),              # Heの初期値の標準偏差
            size=(filter_num, channel_num, filter_size, filter_size)
        )
        self.params[key_b] = np.zeros(filter_num)

        self.layers['Conv4'] = Convolution(
            self.params[key_w], self.params[key_b], stride, pad
        )

        # 次の層の入力サイズ算出
        channel_num = filter_num
        input_size = conv_output_size(input_size, filter_size, pad, stride)

        # [9] ReLU層#4 : レイヤー生成
        self.layers['ReLU4'] = ReLU()

        # [10] プーリング層#2 : レイヤー生成
        self.layers['Pool2'] = Pooling(pool_h=2, pool_w=2, stride=2)

        # 次の層の入力サイズ算出
        input_size = pool_output_size(input_size, pool_size=2, stride=2)

        # [11] 畳み込み層#5 : パラメーター初期化とレイヤー生成
        filter_num, filter_size, pad, stride = list(conv_param_5.values())
        pre_node_num = channel_num * (filter_size ** 2)     # 1ノードに対する前層の接続ノード数
        key_w, key_b = 'W5', 'b5'                           # 辞書格納時のkey
        self.params[key_w] = np.random.normal(
            scale=np.sqrt(2.0 / pre_node_num),              # Heの初期値の標準偏差
            size=(filter_num, channel_num, filter_size, filter_size)
        )
        self.params[key_b] = np.zeros(filter_num)

        self.layers['Conv5'] = Convolution(
            self.params[key_w], self.params[key_b], stride, pad
        )

        # 次の層の入力サイズ算出
        channel_num = filter_num
        input_size = conv_output_size(input_size, filter_size, pad, stride)

        # [12] ReLU層#5 : レイヤー生成
        self.layers['ReLU5'] = ReLU()

        # [13] 畳み込み層#6 : パラメーター初期化とレイヤー生成
        filter_num, filter_size, pad, stride = list(conv_param_6.values())
        pre_node_num = channel_num * (filter_size ** 2)     # 1ノードに対する前層の接続ノード数
        key_w, key_b = 'W6', 'b6'                           # 辞書格納時のkey
        self.params[key_w] = np.random.normal(
            scale=np.sqrt(2.0 / pre_node_num),              # Heの初期値の標準偏差
            size=(filter_num, channel_num, filter_size, filter_size)
        )
        self.params[key_b] = np.zeros(filter_num)

        self.layers['Conv6'] = Convolution(
            self.params[key_w], self.params[key_b], stride, pad
        )

        # 次の層の入力サイズ算出
        channel_num = filter_num
        input_size = conv_output_size(input_size, filter_size, pad, stride)

        # [14] ReLU層#6 : レイヤー生成
        self.layers['ReLU6'] = ReLU()

        # [15] プーリング層#3 : レイヤー生成
        self.layers['Pool3'] = Pooling(pool_h=2, pool_w=2, stride=2)

        # 次の層の入力サイズ算出
        input_size = pool_output_size(input_size, pool_size=2, stride=2)

        # [16] Affine層#1　: パラメーター初期化とレイヤー生成
        pre_node_num = channel_num * (input_size ** 2)      # 1ノードに対する前層の接続ノード数
        key_w, key_b = 'W7', 'b7'                           # 辞書格納時のkey
        self.params[key_w] = np.random.normal(
            scale=np.sqrt(2.0 / pre_node_num),              # Heの初期値の標準偏差
            size=(channel_num * (input_size ** 2), hidden_size)
        )
        self.params[key_b] = np.zeros(hidden_size)

        self.layers['Affine1'] = Affine(self.params[key_w], self.params[key_b])

        # 次の層の入力サイズ算出
        input_size = hidden_size

        # [17] ReLU層#7 : レイヤー生成
        self.layers['ReLU7'] = ReLU()

        # [18] Dropout層#１ : レイヤー生成
        self.layers['Drop1'] = Dropout(dropout_ratio=0.5)

        # [19] Affine層#2　: パラメーター初期化とレイヤー生成
        pre_node_num = input_size                           # 1ノードに対する前層の接続ノード数
        key_w, key_b = 'W8', 'b8'                           # 辞書格納時のkey
        self.params[key_w] = np.random.normal(
            scale=np.sqrt(2.0 / pre_node_num),              # Heの初期値の標準偏差
            size=(input_size, output_size)
        )
        self.params[key_b] = np.zeros(output_size)

        self.layers['Affine2'] = Affine(self.params[key_w], self.params[key_b])

        # [20] Dropout層#2 : レイヤー生成
        self.layers['Drop2'] = Dropout(dropout_ratio=0.5)

        # [21] Softmax層 : レイヤー生成
        self.lastLayer = SoftmaxWithLoss()

    def predict(self, x, train_flg=False):
        """ニューラルネットワークによる推論

        Args:
            x (numpy.ndarray): ニューラルネットワークへの入力
            train_flg (Boolean): 学習中ならTrue（Dropout層でニューロンの消去を実施）

        Returns:
            numpy.ndarray: ニューラルネットワークの出力
        """
        # レイヤーを順伝播
        for layer in self.layers.values():
            if isinstance(layer, Dropout):
                x = layer.forward(x, train_flg)  # Dropout層の場合は、学習中かどうかを伝える
            else:
                x = layer.forward(x)
        return x

    def loss(self, x, t):
        """損失関数の値算出

        Args:
            x (numpy.ndarray): ニューラルネットワークへの入力
            t (numpy.ndarray): 正解のラベル

        Returns:
            float: 損失関数の値
        """
        # 推論
        y = self.predict(x, True)   # 損失は学習中しか算出しないので常にTrue

        # Softmax-with-Lossレイヤーの順伝播で算出
        loss = self.lastLayer.forward(y, t)

        return loss

    def accuracy(self, x, t, batch_size=100):
        """認識精度算出
        batch_sizeは算出時のバッチサイズ。一度に大量データを算出しようとすると
        im2colでメモリを食い過ぎてスラッシングが起きてしまい動かなくなるため、
        その回避のためのもの。

        Args:
            x (numpy.ndarray): ニューラルネットワークへの入力
            t (numpy.ndarray): 正解のラベル（one-hot）
            batch_size (int), optional): 算出時のバッチサイズ、デフォルトは100。

        Returns:
            float: 認識精度
        """
        # 分割数算出
        batch_num = max(int(x.shape[0] / batch_size), 1)

        # 分割
        x_list = np.array_split(x, batch_num, 0)
        t_list = np.array_split(t, batch_num, 0)

        # 分割した単位で処理
        correct_num = 0  # 正答数の合計
        for (sub_x, sub_t) in zip(x_list, t_list):
            assert sub_x.shape[0] == sub_t.shape[0], '分割境界がずれた？'
            y = self.predict(sub_x, False)  # 認識精度は学習中は算出しないので常にFalse
            y = np.argmax(y, axis=1)
            t = np.argmax(sub_t, axis=1)
            correct_num += np.sum(y == t)

        # 認識精度の算出
        return correct_num / x.shape[0]

    def gradient(self, x, t):
        """重みパラメーターに対する勾配を誤差逆伝播法で算出

         Args:
            x (numpy.ndarray): ニューラルネットワークへの入力
            t (numpy.ndarray): 正解のラベル

        Returns:
            dictionary: 勾配を格納した辞書
        """
        # 順伝播
        self.loss(x, t)     # 損失値算出のために順伝播する

        # 逆伝播
        dout = self.lastLayer.backward()
        for layer in reversed(list(self.layers.values())):
            dout = layer.backward(dout)

        # 各レイヤーの微分値を取り出し
        grads = {}
        layer = self.layers['Conv1']
        grads['W1'], grads['b1'] = layer.dW, layer.db
        layer = self.layers['Conv2']
        grads['W2'], grads['b2'] = layer.dW, layer.db
        layer = self.layers['Conv3']
        grads['W3'], grads['b3'] = layer.dW, layer.db
        layer = self.layers['Conv4']
        grads['W4'], grads['b4'] = layer.dW, layer.db
        layer = self.layers['Conv5']
        grads['W5'], grads['b5'] = layer.dW, layer.db
        layer = self.layers['Conv6']
        grads['W6'], grads['b6'] = layer.dW, layer.db
        layer = self.layers['Affine1']
        grads['W7'], grads['b7'] = layer.dW, layer.db
        layer = self.layers['Affine2']
        grads['W8'], grads['b8'] = layer.dW, layer.db

        return grads

(5)学習の実装

学習は前の章のコードとほとんど変わりません。本のコードに合わせて Trainer クラスを実装しようかと思っていたのですが、もう最後の章になってしまって今回の実装で終わりなので従来のままにしています。

更新回数は12,000（20エポック）にしてみました。

mnist.py

# coding: utf-8
import os
import sys
import matplotlib.pylab as plt
import numpy as np
from adam import Adam
from deep_conv_net import DeepConvNet
sys.path.append(os.pardir)  # パスに親ディレクトリ追加
from dataset.mnist import load_mnist


# MNISTの訓練データとテストデータ読み込み
(x_train, t_train), (x_test, t_test) = \
    load_mnist(normalize=True, flatten=False, one_hot_label=True)

# ハイパーパラメーター設定
iters_num = 12000           # 更新回数
batch_size = 100            # バッチサイズ
adam_param_alpha = 0.001    # Adamのパラメーター
adam_param_beta1 = 0.9      # Adamのパラメーター
adam_param_beta2 = 0.999    # Adamのパラメーター

train_size = x_train.shape[0]  # 訓練データのサイズ
iter_per_epoch = max(int(train_size / batch_size), 1)    # 1エポック当たりの繰り返し数

# ディープな畳み込みニューラルネットワーク生成
network = DeepConvNet()

# オプティマイザー生成、Adamを使用
optimizer = Adam(adam_param_alpha, adam_param_beta1, adam_param_beta2)

# 学習前の認識精度の確認
train_acc = network.accuracy(x_train, t_train)
test_acc = network.accuracy(x_test, t_test)
train_loss_list = []            # 損失関数の値の推移の格納先
train_acc_list = [train_acc]    # 訓練データに対する認識精度の推移の格納先
test_acc_list = [test_acc]      # テストデータに対する認識精度の推移の格納先
print(f'学習前 [訓練データの認識精度]{train_acc:.4f} [テストデータの認識精度]{test_acc:.4f}')

# 学習開始
for i in range(iters_num):

    # ミニバッチ生成
    batch_mask = np.random.choice(train_size, batch_size, replace=False)
    x_batch = x_train[batch_mask]
    t_batch = t_train[batch_mask]

    # 勾配の計算
    grads = network.gradient(x_batch, t_batch)

    # 重みパラメーター更新
    optimizer.update(network.params, grads)

    # 損失関数の値算出
    loss = network.loss(x_batch, t_batch)
    train_loss_list.append(loss)

    # 1エポックごとに認識精度算出
    if (i + 1) % iter_per_epoch == 0:
        train_acc = network.accuracy(x_train, t_train)
        test_acc = network.accuracy(x_test, t_test)
        train_acc_list.append(train_acc)
        test_acc_list.append(test_acc)

        # 経過表示
        print(
            f'[エポック]{(i + 1) // iter_per_epoch:>2} '
            f'[更新数]{i + 1:>5} [損失関数の値]{loss:.4f} '
            f'[訓練データの認識精度]{train_acc:.4f} [テストデータの認識精度]{test_acc:.4f}'
        )

# 損失関数の値の推移を描画
x = np.arange(len(train_loss_list))
plt.plot(x, train_loss_list, label='loss')
plt.xlabel('iteration')
plt.ylabel('loss')
plt.xlim(left=0)
plt.ylim(0, 2.5)
plt.show()

# 訓練データとテストデータの認識精度の推移を描画
x2 = np.arange(len(train_acc_list))
plt.plot(x2, train_acc_list, label='train acc')
plt.plot(x2, test_acc_list, label='test acc', linestyle='--')
plt.xlabel('epochs')
plt.ylabel('accuracy')
plt.xlim(left=0)
plt.ylim(0, 1.0)
plt.legend(loc='lower right')
plt.show()

(6)実行結果

以下、実行結果です。私の環境では半日くらいかかりました。

学習前 [訓練データの認識精度]0.0975 [テストデータの認識精度]0.0974
[エポック] 1 [更新数]  600 [損失関数の値]1.0798 [訓練データの認識精度]0.9798 [テストデータの認識精度]0.9811
[エポック] 2 [更新数] 1200 [損失関数の値]0.8792 [訓練データの認識精度]0.9881 [テストデータの認識精度]0.9872
[エポック] 3 [更新数] 1800 [損失関数の値]0.9032 [訓練データの認識精度]0.9884 [テストデータの認識精度]0.9890
[エポック] 4 [更新数] 2400 [損失関数の値]0.8012 [訓練データの認識精度]0.9914 [テストデータの認識精度]0.9906
[エポック] 5 [更新数] 3000 [損失関数の値]0.9475 [訓練データの認識精度]0.9932 [テストデータの認識精度]0.9907
[エポック] 6 [更新数] 3600 [損失関数の値]0.8105 [訓練データの認識精度]0.9939 [テストデータの認識精度]0.9910
[エポック] 7 [更新数] 4200 [損失関数の値]0.8369 [訓練データの認識精度]0.9920 [テストデータの認識精度]0.9915
[エポック] 8 [更新数] 4800 [損失関数の値]0.8727 [訓練データの認識精度]0.9954 [テストデータの認識精度]0.9939
[エポック] 9 [更新数] 5400 [損失関数の値]0.9640 [訓練データの認識精度]0.9958 [テストデータの認識精度]0.9935
[エポック]10 [更新数] 6000 [損失関数の値]0.8375 [訓練データの認識精度]0.9953 [テストデータの認識精度]0.9925
[エポック]11 [更新数] 6600 [損失関数の値]0.8500 [訓練データの認識精度]0.9955 [テストデータの認識精度]0.9915
[エポック]12 [更新数] 7200 [損失関数の値]0.7959 [訓練データの認識精度]0.9966 [テストデータの認識精度]0.9932
[エポック]13 [更新数] 7800 [損失関数の値]0.7778 [訓練データの認識精度]0.9946 [テストデータの認識精度]0.9919
[エポック]14 [更新数] 8400 [損失関数の値]0.9212 [訓練データの認識精度]0.9973 [テストデータの認識精度]0.9929
[エポック]15 [更新数] 9000 [損失関数の値]0.9046 [訓練データの認識精度]0.9974 [テストデータの認識精度]0.9934
[エポック]16 [更新数] 9600 [損失関数の値]0.9806 [訓練データの認識精度]0.9970 [テストデータの認識精度]0.9924
[エポック]17 [更新数]10200 [損失関数の値]0.7837 [訓練データの認識精度]0.9975 [テストデータの認識精度]0.9931
[エポック]18 [更新数]10800 [損失関数の値]0.8948 [訓練データの認識精度]0.9976 [テストデータの認識精度]0.9928
[エポック]19 [更新数]11400 [損失関数の値]0.7936 [訓練データの認識精度]0.9980 [テストデータの認識精度]0.9932
[エポック]20 [更新数]12000 [損失関数の値]0.8072 [訓練データの認識精度]0.9984 [テストデータの認識精度]0.9939

最終的な認識精度は99.39%でした。前の章のCNNは98.60%だったので、0.79ポイント アップです。層を深くすることの可能性を感じさせる結果になりました。

なお、前の章の結果と比べて認識精度に対する損失関数の値が大きいのですが、これはDropoutの影響かと思います。認識精度はすべてのニューロンを使っていますが、損失関数の算出時は半分（Dropoutのレートを0.5で実行したので）のニューロンが削除状態のためです。

この本で実装するのはここまでですが、さらに認識精度を高めるために、アンサンブル学習やData Augmentationなどの手法が紹介されています。また、層を深くすることのメリットについてもまとめられています。

8.2 ディープラーニングの小歴史

ディープラーニングのトレンドの紹介です。いずれも、ここまでに学んだCNNが基本であることが理解できました。

8.3 ディープラーニングの高速化

高速化についての説明です。興味深かったのは、ディープラーニングだと単精度浮動小数点では精度が高すぎてもったいないので、半精度の浮動小数点が注目されているという点です。これまで私が使ってきた開発言語では半精度浮動小数点型というのを聞いたことがなかったのですが、NumPyにはfloat16という型があることを知りました。

8.4 ディープラーニングの実用例

物体検出、セグメンテーション、画像のキャプション生成と、面白そうなことがすでに実現されていることが分かりました。ただ、その仕組みについては、これまでに学んだレベルではまだまだ理解しきれません。

8.5 ディープラーニングの未来

画像の生成や自動運転、強化学習など、研究中の分野の紹介です。ディープラーニングの可能性を感じますね。

8.6 まとめ

なんとか最後の実装も終えました。本の通りの精度が出せて一安心です。ディープラーニングの可能性についても学ぶことができました。

この章は以上です。誤りなどありましたら、ご指摘いただけますとうれしいです。

最後に

この本のおかげで、ディープラーニングへの理解が深まりました。著者の斎藤康毅さんに感謝しています。また、みなさんのからのいいね、ストック、はてなブックマークなどに励まされ、途中であきらめずに読破することができました。ありがとうございます。
このメモが、後に続く方の参考になることがあれば幸いです。

（このメモの他の章へ：1章 / 2章 / 3章 / 4章 / 5章 / 6章 / 7章 / 8章）

AtCorder Beginner Contest 153 にバーチャル参加しました。
ABCDの4問ACでした。
Python3を使用しています。

A問題

サーバルはモンスターと戦っています。
モンスターの体力は H です。
サーバルが攻撃を 1 回行うとモンスターの体力を A 減らすことができます。攻撃以外の方法でモンスターの体力を減らすことはできません。
モンスターの体力を 0 以下にすればサーバルの勝ちです。
サーバルがモンスターに勝つために必要な攻撃の回数を求めてください。

モンスターの体力が何回の攻撃で 0 以下になるかを考えるので、割り算の商を求める必要があると考えます。
しかし、回数は整数でないとおかしいため、割り算の商が小数を取った場合は繰り上げた整数を答える必要があります。

H, A = map(int,input().split())
if H % A == 0:
    print(H // A)
else:
    print(H // A + 1)

小数点切り上げだとif文も使わなくて良かった。

import math
H, A = map(int,input().split())
print(math.ceil( H // A ))

B問題

アライグマはモンスターと戦っています。
モンスターの体力は H です。
アライグマは N 種類の必殺技を使うことができ、i 番目の必殺技を使うとモンスターの体力を Ai 減らすことができます。 必殺技を使う以外の方法でモンスターの体力を減らすことはできません。
モンスターの体力を 0 以下にすればアライグマの勝ちです。
アライグマが同じ必殺技を 2 度以上使うことなくモンスターに勝つことができるなら Yes を、できないなら No を出力してください。
必殺技を 1 回しか使えないという条件から、全ての必殺技を使ったときにモンスターを倒せるのか考えます。

全ての必殺技で減らせる体力がモンスターの体力より大きければ勝つことができます。

H, N = map(int,input().split())
A = list(map(int,input().split()))
if sum(A) >= H:
    print("Yes")
else:
    print("No")

C問題

フェネックは N 体のモンスターと戦っています。
i 番目のモンスターの体力は Hi です。
フェネックは次の 2 種類の行動を行うことができます。
・攻撃：モンスターを 1 体選んで攻撃することで、そのモンスターの体力を 1 減らす
・必殺技：モンスターを 1 体選んで必殺技を使うことで、そのモンスターの体力を 0 にする
攻撃と必殺技以外の方法でモンスターの体力を減らすことはできません。
全てのモンスターの体力を 0 以下にすればフェネックの勝ちです。
フェネックが K 回まで必殺技を使えるとき、モンスターに勝つまでに行う攻撃の回数 (必殺技は数えません) の最小値を求めてください。

モンスターの体力を 0 にする必殺技は、体力の大きいモンスターに使いたいと考えます。
そこでモンスターを体力の大きい順(降順)に並べ、 K 体のモンスターを必殺技で倒します。
残った N - K 体のモンスターは攻撃で倒すため、この N - K 体のモンスターの総合体力の回数だけ攻撃を行う必要があります。

N, K = map(int,input().split())
A = list(map(int,input().split()))
A.sort(reverse = True)
del A[:K]
print(sum(A))

D問題

カラカルはモンスターと戦っています。
モンスターの体力は H です。
カラカルはモンスターを 1 体選んで攻撃することができます。モンスターを攻撃したとき、攻撃対象のモンスターの体力に応じて、次のどちらかが起こります。
・モンスターの体力が 1 なら、そのモンスターの体力は 0 になる
・モンスターの体力が X > 1 なら、そのモンスターは消滅し、体力が ⌊X/2⌋ のモンスターが新たに 2 体現れる
（⌊r⌋ は r を超えない最大の整数を表す）
全てのモンスターの体力を 0 以下にすればカラカルの勝ちです。
カラカルがモンスターに勝つまでに行う攻撃の回数の最小値を求めてください。

最初は体力 H のモンスターが 1 体いる、と問題文から読み取ります。
1 回の攻撃で体力が半分になったモンスターが 2 体になる(小数点以下の値は切り捨て)、という操作を繰り返します。
モンスターの体力 H が 2 の累乗でない場合、 H は約数に奇数を持ちます。
そのため体力を半分にする過程で小数点の切り捨てが起こり、モンスターが分裂するだけで、総合体力が減ります。
そこで、体力 H が 2 の累乗で攻撃回数に変化が起きるのではと考えました。

ここからは実際に数値計算し、一般論を導きます。

H 攻撃回数
1 1 回
2 ~ 3 3回
4 ~ 7 7回
8 ~ 15 15回
H が 2ⁿ から 2ⁿ⁺¹ - 1 のとき、攻撃回数は 2⁰ から 2ⁿ までの和で表せるということが分かります。

そこで、今度は H が 2 の n 乗と 2 の n + 1 乗 の間にあるかを考えます。
10 進数を 2 進数に変換して文字数を数えることで求めます。

H = int(input())
n = len(bin(H)) - 3
N = [2 ** i for i in range(n + 1)]
print(sum(N))

H が 2ⁿ から 2ⁿ⁺¹ - 1 のとき、攻撃回数は 2ⁿ⁺¹ - 1 で表せた。
さらに、底が 2 の log を取って小数点以下を切り捨てることで n が求められた。

import math
H = int(input())
n = math.floor(math.log2(H))
print(2 ** (n + 1) -1)

E問題

トキはモンスターと戦っています。
モンスターの体力は H です。
トキは N 種類の魔法が使え、i 番目の魔法を使うと、モンスターの体力を Ai 減らすことができますが、トキの魔力を Bi 消耗します。
同じ魔法は何度でも使うことができます。魔法以外の方法でモンスターの体力を減らすことはできません。
モンスターの体力を 0 以下にすればトキの勝ちです。
トキがモンスターに勝つまでに消耗する魔力の合計の最小値を求めてください。

どうやら、個数制限ナップサックという有名な問題だったらしい。
DP(動的計画法)もさっぱり分からなかったので、こちらを参考にさせてもらいました。
分かりやすい。

https://qiita.com/drken/items/dc53c683d6de8aeacf5a

まずTLEになった結果です。(Pypy3だとAC)
使える N 種類の魔法のうち、 i 番目の魔法について考えます。
i 番目の魔法を使って体力を合計で j 減らせる場合の消耗魔力は
dp[j] = dp[j - a] + b
i 番目の魔法を使わずに体力を合計で j 減らせる場合の消耗魔力は
dp[j] = dp[j]
となります。
これらを比較して、最小値を求めdp[j]とします。
モンスターの体力ぴったりで倒す必要はないため、dp[]の要素数はH + max(減らせる体力)となっています。
初期値はdp[0] = 0
また、0が最小値になってしまうためif文で場合分けしています。

H, N= map(int,input().split())
A = [0] * N
B = [0] * N
for i in range(N):
    A[i], B[i] = map(int, input().split())
ma = max(A)
dp = [0] * (H + ma)
dp[0] = 0
for i in range(N):
    a = A[i]
    b = B[i]
    for j in range(1, H + ma):
        if dp[j] == 0:
            dp[j] = dp[j - a] + b
        else:
            dp[j] = min(dp[j - a] + b, dp[j])
print(min(dp[H:]))

二重のforループの中でif文は遅くなりそうだなと反省し改善したのですが、TLEです。(Pypy3だとAC)
dp[]の要素に無限大を入れて作っているので場合分けして判定しなくて良くなった。
見やすくなったけど、1 つ目の解答より時間はかかっている。

H, N= map(int,input().split())
A = [0] * N
B = [0] * N
for i in range(N):
    A[i], B[i] = map(int, input().split())
ma = max(A)
dp = [float("inf")] * (H + ma)
dp[0] = 0
for i in range(N):
    a = A[i]
    b = B[i]
    for j in range(1, H + ma):
        dp[j] = min(dp[j - a] + b, dp[j])
print(min(dp[H:]))

PythonでACしたい。

F問題

ギンギツネは N 体のモンスターと戦っています。
モンスターは 1 列に並んでおり、数直線上にいるとみなすことができます。
i 番目のモンスターは座標 Xi にいて、体力は Hi です。
ギンギツネは爆弾を使ってモンスターを攻撃することができます。 座標 x で爆弾を使うと、座標が x−D 以上 x+D 以下の範囲にいる全てのモンスターの体力を A 減らすことができます。 爆弾を使う以外の方法でモンスターの体力を減らすことはできません。
全てのモンスターの体力を 0 以下にすればギンギツネの勝ちです。
ギンギツネがモンスターに勝つまでに爆弾を使う回数の最小値を求めてください。

サンプルだと正しい結果が出るんだけど、WAだった。
距離を昇順にならべて、一番小さい値が0になるまで攻撃する。
それで同時に攻撃される範囲のモンスターの体力を減らして、新たにリストを作り直すと考えた。

import math
N, D, A= map(int,input().split())
X = [0] * N
H = [0] * N
Y = []
for i in range(N):
    X[i], H[i] = map(int, input().split())
    Y.append([X[i] ,math.ceil(H[i] / A)])
Y.sort()
d = 0

while len(Y) > 0:
    y_min = Y[0][1]
    c = 0
    d += y_min
    for i in range(len(Y)):
        if Y[i][0] < Y[0][0] + 2 * D + 1:
            Y[i][1] = Y[i][1] - y_min
            if Y[i][1] <= 0:
                c += 1
    del Y[0:c]
print(d)

1

import tarfile

with tarfile.open('test.tar.gz', 'w:gz') as tr:
    tr.add('test_dir')
#test_dirというディレクトリをtest.tar.gzというファイル名でtarfileに圧縮

with tarfile.open('test.tar.gz', 'r:gz') as tr:
    tr.extractall(path='test_tar')
#test_tarというディレクトリにtest.tar.gzを展開

2

import tarfile

with tarfile.open('test.tar.gz', 'r:gz') as tr:
    with tr.extractfile('test_dir/sub_dir/sub_test.txt') as f:
        print(f.read())
#tarfileに圧縮されたtest.tar.gzを展開することなく、sub_test.txtを出力

はじめに

仕事で Python を使っていて、
「あー、ここのリスト作る処理、リスト内包表記にしたら簡潔に書けるのになー。」
と思うことがあります。
でも、私がリスト内包表記を書きたいと思う場面が意外と少ないため、いざ書こうと思ったときに、
「あれ、リスト内包表記って、どうやって書くんだっけ？うーん、時間もないし、いつも通りループにするか。」
と、残念な感じになるので、私自身のための覚え書きの意味も含めて。

リスト内包表記とは？

初めてこの言葉を聞いたときは、頭の中が？マークだらけに。。。
すごーく簡単に言うと、リスト内包表記とは、
リストデータを１行で作る方法
です。

ただし、何でも1行で書けるわけではなく、（若干語弊がありますが、）ループでリストを作成する処理を1行で書けるものです。
私自身、「仕事で書くコードは複雑なんだよ！」と思っていたのですが、意外と使えます。

書き方

通常のループ（これはリスト内包表記ではありません。）

data = []
for i in range(10):
    data.append(i)

print(data)
# [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

リスト内包表記

print([i for i in range(10)])
# [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

解説

通常のループと、リスト内包表記の結果は同じです。
しかし、リスト内包表記の方が、記述量が少なくてすみます。
そして、コードが簡潔で読みやすくなります。（慣れないと読みにくいですが。）

リスト内包表記のコードを見ていきます。

[...]
ご存知の通り、これでリストができるため、通常のループの data = [] が不要になります。

for i in range(10)
これは通常のループでも、リスト内包表記でも同じです。
ただし、リスト内包表記の方には、 : は不要です。

i for i in range(10)
最初にでてくる i は、ループで取り出した i （2つめの i ）の値が格納されます。
ループのたびに最初にでてくる i の値が書き換わります。
その i が、ループの回数分、リストに追記（append）されます。

もう少し具体的に

もう少し具体的な方が理解しやすいと思うので、次のケースで説明します。

ケース1

リスト内のデータを加工したい場合。
例えば、 ['aa', 'b', 'ccc'] というリストの各データに、 @ を付加したい場合。

items = ['aa', 'b', 'ccc']
print([item + '@' for item in items])
# ['aa@', 'b@', 'ccc@']

items のデータに @ をつけて、別のリストを作成しています。
リストデータから、別のリストデータを１行で作ったというわけです。

ケース2

データ加工処理など、決まったルールで文字列の途中に、別の文字を付加したい場合。
例えば、abcdefgh というデータに1文字ごとに、 - を付加したい場合。

items = 'abcdefgh'
print('-'.join([item for item in items]))
# a-b-c-d-e-f-g-h

こちらは、文字列からリストデータを１行で作っています。
そのリストのデータを - で join して文字列に戻しています。

ちなみに、２文字ごとに、 - を付加したい場合。

items = 'abcdefgh'
print('-'.join([items[i*2: i*2+2] for i in range(int(len(items)/2))]))
# ab-cd-ef-gh

ちょっと詰め込みすぎな感じがするので、こんな感じでしょうか。

items = 'abcdefgh'
loop_count = int(len(items)/2)
print('-'.join([items[i*2: i*2+2] for i in range(loop_count)]))
# ab-cd-ef-gh

いかがでしょうか？

ケース1では、リストのデータを加工して、別のリストを作る、
ケース2では、文字列を加工して、リストを作る、というサンプルです。

どちらのケースも、ループを使って実現できますが、
リスト内包表記を使えば、シンプルかつ簡潔にコードが書けます。

はじめに

開発のテストではライブラリやフレームワークを使用して自動化するのが一般的になっています。そこでpytestを使用してflaskの単体テストを自動化しようとしましたが、シンプルな例が見つけられなかったのでシンプルな例と簡単な説明をまとめました。

環境

• python:3.6.5
• flask:1.0.2
• pytest:5.3.5

インストール

pip install pytestでインストールするだけです。

pytestで自動化するのに必要なもの

pytestで単体テストを自動化するために必要なものは、テスト対象のソース(テストされる開発物)とテスト方法を記載したソースが必要になります。テスト方法のソースは、テスト対象の引数と関数の結果を与えてどのように比較するかを記載しています。

簡単な関数の単体テスト自動化

flaskの単体テストの自動化の前に、簡単な関数を通じてpytestの使い方を見ていきます。

テスト対象のソース

テスト対象のソースが無ければテストはできないため、テスト対象のソースを用意します。
例では、引数を加算して返却する関数を用意しましたが、本来の開発であれば開発物が相当します。

testing_mod.py

def add_calc(a, b):
    return a + b

テスト方法を書いたソース

テスト対象のソースの呼び方とソースの結果を書いたソースを作ります。このソースはテスト対象のソースの関数を呼び出して、テスト対象の関数が返した結果とこちらが想定した結果を比較して正しければOK、誤っていればNGになります。
例では、テスト対象のtesting_modをimport testing_mod でimportしてtesting_mod.add_calc()に1と2を渡して返却される結果が3であればOKとなっています。

py_test_main.py

import pytest
import testing_mod

def test_ok_sample():
    result = testing_mod.add_calc(1, 2)
    assert 3 == result

単体テスト実行

テスト対象とテスト方法のソースができたため、関数ごとに結果を見たいため-vオプションをつけて実行します。

# pytest -v py_test_main.py

py_test_main.py .         [100%]                                                                                                                                       
====== 1 passed in 0.05s ======
PS C:\Users\xxxx\program\python> pytest -v py_test_main.py                                   
====== test session starts ======
platform win32 -- Python 3.6.5, pytest-5.3.5, py-1.8.1, pluggy-0.13.1 -- c:\users\xxxx\appdata\local\programs\python\python36-32\python.exe
cachedir: .pytest_cache
rootdir: C:\Users\xxxx\program\python
collected 1 item                                                                                                                     
py_test_main.py::test_ok_sample PASSED     [100%] 

====== 1 passed in 0.02s ======

結果を見ると先ほど作成したtest_ok_sampleがPASSEDとなり正常に終わったため、テストOKになります。
関数をたくさん作成すると表示される関数の数が増えていきます。

簡単なflaskの単体テスト自動化

flaskの単体テストを自動化します。上の簡単な関数の自動化とは異なり、flaskはクライアントからの通信を必要としますが、単体テストではflaskの機能を使用して単体テストの自動化をします。

テスト対象のソース

flaskのソースを作成します。例では/にアクセスするとroot文字列を返却するものを作成します。本来の開発であれば開発物が相当します。
flaskについては以前のflaskについてまとめたものを参照してください。

flask_mod.py

from flask import Flask, jsonify

app = Flask(__name__)

@app.route('/')
def root():
    return "root"

テスト方法を書いたソース

flaskのテスト方法を書いたソースは、関数のソースと異なりflaskのクライアントを生成してからそのクライアントを使用してリクエストを発行して結果を確認する必要があります。

テスト用flaskクライアントの生成

まずは、テスト用のクライアントの生成を行います。テスト対象ソースのappをimportしてappのテスト用configをtrueに変更します。その後appのtest_client()を使用してクライアントを生成します。
下の例でいうとテスト対象のソースのimportはfrom flask_mod import appになります。

py_test_main.py

import pytest
from flask_mod import app

def test_flask_simple():
    app.config['TESTING'] = True
    client = app.test_client() 

テスト対象の関数の実行

上で生成したクライアントを使用して、テスト対象のURLに向けてget関数やpost関数を使用してリクエストを発行します。その結果がflaskからのレスポンスになるため、期待した答えになるかをpytestのassertでチェックします。
下の例でいうとresult = client.get('/') で/にgetリクエストを発行してその結果がresultに格納されるのでdata(body)とrootを比較しています。

py_test_main.py

import pytest
from flask_mod import app

def test_flask_simple():
    app.config['TESTING'] = True
    client = app.test_client() 
    result = client.get('/')
    assert b'root' == result.data

単体テスト実行

テスト対象とテスト方法のソースができたため、実行します。

# pytest -v py_test_main.py

====== 1 passed in 0.22s =======
PS C:\Users\xxxx\program\python\flask> pytest -v .\pytest_flask.py
====== test session starts ======
platform win32 -- Python 3.6.5, pytest-5.3.5, py-1.8.1, pluggy-0.13.1 -- c:\users\xxxx\appdata\local\programs\python\python36-32\python.exe
cachedir: .pytest_cache
rootdir: C:\Users\xxxx\program\python\flask
collected 1 item

pytest_flask.py::test_flask_simple PASSED  [100%]

====== 1 passed in 0.20s =======

結果を見ると先ほど作成したtest_flask_simpleがPASSEDとなり正常に終わったため、テストOKになります。
関数をたくさん作成するとここの関数の数が増えていきます。

単体テストがエラーの時の例

試しに、flaskが返却してくる文字列と比較する文字列をsampleにしたときの結果を見てみます。

# pytest -v pytest_flask.py
======= test session starts =======
platform win32 -- Python 3.6.5, pytest-5.3.5, py-1.8.1, pluggy-0.13.1 -- c:\users\xxxx\appdata\local\programs\python\python36-32\python.exe
cachedir: .pytest_cache
rootdir: C:\Users\xxxx\program\python\flask
collected 1 item

pytest_flask.py::test_flask_simple FAILED                                                                                                                          [100%]

============ FAILURES ============= 
____________ test_flask_simple ____________

    def test_flask_simple():
        app.config['TESTING'] = True
        client = app.test_client()
        result = client.get('/')
>       assert b'sample' == result.data
E       AssertionError: assert b'sample' == b'root'
E         At index 0 diff: b's' != b'r'
E         Full diff:
E         - b'sample'
E         + b'root'

pytest_flask.py:8: AssertionError
======== 1 failed in 0.26s ========

ちゃんとsampleとrootが異なるのでAssertionError: assert b'sample' == b'root'と表示されました。

おわりに

単体テストの自動化は、自動化スクリプトの作成の労力が少なければとても便利な仕組みになります。
その労力を少なくする方法がフレームワークですが、上記以外にも、テストの前処理と後処理をする、同じテスト方法で複数のパラメータを試すなど便利な機能があります。次はその便利な方法をまとめていきます。

中国語形態素解析エンジンjieba

jiebaのPython版で使ってみました。
他のプログラミング言語のバージョンもあります。

インストール

$ pip install jieba

テキストセグメンテーション

>>> import jieba
>>> text = "我明天去东京大学上课。早上十点开始。"
#"私は明日東京大学の授業に出ます。朝の十時から。"

jieba.cutの戻り値はジェネレータ
jieba.lcutの戻り値はリスト
jieba.cut_for_searchの戻り値はジェネレータ
jieba.lcut_for_searchの戻り値はリスト

Accurate Mode

>>> segments = jieba.cut(text)
>>> list(segments)
['我', '明天', '去', '东京大学', '上课', '。', '早上', '十点', '开始', '。']

>>> segments = jieba.lcut(text)
>>> segments
['我', '明天', '去', '东京大学', '上课', '。', '早上', '十点', '开始', '。']

东京大学がひと単語になっています、いいですね！

Full Mode

cut_all=Trueにする。

>>> segments = jieba.cut(text, cut_all=True)
>>> list(segments)
['我', '明天', '去', '东京', '东京大学', '大学', '学上', '上课', '。', '早上', '十点', '开始', '。']

>>> segments = jieba.lcut(text, cut_all=True)
>>> segments
['我', '明天', '去', '东京', '东京大学', '大学', '学上', '上课', '。', '早上', '十点', '开始', '。']

Search Engine Mode

>>> segments = jieba.cut_for_search(text)
>>> list(segments)
['我', '明天', '去', '东京', '大学', '东京大学', '上课', '。', '早上', '十点', '开始', '。']

>>> segments = jieba.lcut_for_search(text)
>>> segments
['我', '明天', '去', '东京', '大学', '东京大学', '上课', '。', '早上', '十点', '开始', '。']

キーワード抽出

>>> import jieba.analyse
>>> text = '''
... 伴随全球化进程不断加速，人类正面临着日益严峻的挑战，能源危机、环境污染、金融环境不稳定以及贫穷等各种问题不断显现在我们面前。应对这些挑战，各类人才需要通力协作，献力献策，共同攻克这些全球化难题。在这种背景下，作为领导型人才的摇篮，东京大学肩负着义不容辞的责任。我们将会以无限的勇气、智慧与责任感迎面直对这些挑战。
... 学贵精专，攀登学术巅峰毋庸置疑成为治学的原动力。对于在治学道路上苦苦求索的年轻学者和广大学生而言，他们不断前进、不断成长的动力源自于获取最精尖知识的一种兴奋感与喜悦感。东京大学不负国民所托，为这些学生、学者提供了成长的乐土，使其更好地造福社会。
... 东京大学创立至今，一直保持着东西方文化相结合的学术观点，不断发展，放眼世界，形成一面独特的旗帜。秉承传统，展望未来，东京大学吸引了各式人才汇聚于此，探索求知，灵变创新，将逐步成为“全球性知识创新与协作的据点”。在东大，国界、文化、年龄的壁垒通通被打破，新领域学术研究超越文理界限，产官学合作在此通力展开。为实现这一目标，首先需要创立卓越性与国际性二者兼备的研究生院，并积极开展跨学科的新学术领域拓展。... 东京大学宪章中规定，东京大学致力于为世界和平与人类福祉作出不朽的贡献。而现代社会正飞速发展，我们需要顺应时代发展的要求，将学术研究赋予新时代的意义。在坚守传统的同时，体制改革亦不可或缺。在深化本科生教育改革的同时，不断推进研究生院的根本变革，使其创造知识价值主体的意义显现出来。此外，亟需推进人事制度的改革创新，大力倡导男女平等，多给年轻人创造学以致用的机会， 并实现兼顾人力资源的流动性与稳定性。 不可忽视的一个问题是，推进上述改革的前提是提升整个社会对科研学术的信赖感，而这种信赖感建立在强化科研道德水准，推动科研成果转化的基础之上。
... 通过不断推陈出新，我们致力于将东京大学发展成为为日本民众所爱戴，乃至为整个世界所爱戴的东京大学。
... '''

テキストは東京大学総長談論中国語版になります。

tf-idf値による抽出

>>> keywords = jieba.analyse.extract_tags(text, topK=20, withWeight=False, allowPOS=())
>>> keywords
['东京大学', '不断', '信赖感', '学术', '挑战', '人才', '治学', '知识', '研究生院', '爱戴', '学术研究', '创新', '推进', '全球化', '改革', '科研', '这些', '致力于', '喜悦感', '创立']

良さそうですね。日本の漢字と少し違いますが、だいたい読めます。

TextRankに基づく抽出

>>> keywords = jieba.analyse.textrank(text, topK=20, withWeight=False, allowPOS=('ns', 'n', 'vn', 'v'))
>>> keywords
['成为', '知识', '创新', '学术', '发展', '需要', '改革', '人才', '推进', '科研', '挑战', '实现', '领域', '意义', '社会', '学术研究', '人类', '文化', '治学', '勇气']

その他

他にもたくさんの機能を持っていて、
辞書をいじったり、品詞タグ付けしたり、などの事ができるので、詳しくは公式をみた方が良さそうです。README.mdの前半は中国語ですが、後半は英語訳になっています。

筆者は東京大学と全く関係ありあません。

解決したいこと

・DjangoでHttpステータスコードを指定してreturnさせる。
・jQueryでAjax通信を行なっており、特定条件で強制的にfailを実行させたい。

環境

・Windows10 premium
・VScode
・jQuery

結果

Ajaxのコードは省略
仮に適当な文字列が送信され、送信された文字列内にテストという文字が含まれていない場合の処理とする。

test.py

import HttpResponse

get_post = request.POST["testData"]

if not "テスト" in get_post:
    response = HttpResponse(status=500) # ステータスコードに500を指定
    return response

参考

HTTPステータスコード一覧

解決したいこと

・Djangoでアップロードしたファイルの読み込みを行い、何らかの操作を行うという処理をcronで回そうとした際、ファイルが読めないとパーミッションエラーが発生。
・ファイルアップロード時の権限は「600」

試したこと

前提としてchmodを行う際、管理者権限が必要な環境です。

os.py

import os
os.system('sudo shmod 644 path')

python os.py -> 権限切り替わる。
cronで実行 -> 権限切り替わらず。

subprocess.py

import subprocess
subprocess.call('sudo shmod 644 path')

python subprocess.py -> 権限切り替わる。
cronで実行 -> 権限切り替わらず。

解決策

settings.pyに「FILE_UPLOAD_PERMISSIONS = 0o644」と記載しておくことでアップロードしたファイルの権限が「644」になった。
セキュリティ考慮してデフォルト「600」にしているっぽいですが、時間かかった。。

cron回したときにsubprocessやos.systemで切り替わらなかった理由は分かってないです。

はじめに

古典的な計算機上で素因数分解をするには愚直には$O(\sqrt{N})$の計算量がかかり、現在開発されている最も効率的なアルゴリズムでも桁数の指数時間程度はかかってしまうようです。でも300桁レベルの素因数分解をしようと思うと、死ぬほど時間がかかってしまいます。え？多項式時間で素因数分解を？できらぁ！というのがショアのアルゴリズムです。この記事ではそのショアのアルゴリズムをqiskitで実装して、実際の量子コンピュータで実行するのがこの記事の目標です。今回は15を因数分解していきます。この記事にはおそらく（というか確実に）間違いや、誤植があると思うので、もし見つけたら優しく教えてください...

流れ

ショアのアルゴリズムを説明するには、いろいろ準備が必要なのですがまず全体の流れを書いておきたいと思います。

まず、素因数分解したい自然数$N$と互いに素な自然数$x$を考えます（互いに素かどうかはユークリッドの互除法で判定できます。）。以下modNの世界で考えることにします。$x^r = 1$となるような最小のrを位数と呼びます。この位数を求めることが素因数分解に有効であることが知られています。しかし古典的な計算機ではこの計算にも指数時間がかかってしまいます。これを解決するのが位相推定アルゴリズムです。位相推定アルゴリズムはユニタリ行列の固有値（絶対値は1）の位相を近似的に与えてくれます。$\mod N$の世界で$U^r = I$となるようなUを考えるとその固有値からrを求めることができるようになります。この位相推定アルゴリズムのサブルーチンとしてこのあと説明するアダマールテストと量子フーリエ変換が必要です。

アダマールテスト

まずアダマールテストと呼ばれるアルゴリズムを説明します。

アダマールテストは以下のような回路で実現されます（$U$は任意のゲートです）。

この回路に$q_0 = \left|0\right>, q_1 = \left|\varphi\right>$ を通すとどうなるか考えてみましょう。ただし$\left|\varphi\right>$は$U$の固有ベクトルであり、その固有値を$e^{i\lambda}$とします。
$$
\left|0\right>\left|\varphi\right>\rightarrow \frac{1}{\sqrt{2}}(\left|0\right>\left|\varphi\right>+\left|1\right>\left|\varphi\right>)\rightarrow
\frac{1}{\sqrt{2}}(\left|0\right>\left|\varphi\right>+e^{i\lambda}\left|1\right>\left|\varphi\right>)\rightarrow
\left(\frac{1+e^{i\lambda}}{2}\left|0\right>+\frac{1-e^{i\lambda}}{2}\left|1\right>\right)\otimes\left|\varphi\right>
$$

control-Uゲートを通過した後の量子ビットは$\frac{1}{\sqrt{2}}(\left|0\right>+e^{i\lambda}\left|1\right>)\otimes \left|\varphi\right>$という状態になります。Uの固有値の位相が前（第一量子ビット）に位相として出てきているのがわかると思います。そこにさらにアダマールゲートを作用させることで、第一量子ビットの観測確率としてUの固有値の位相を得ることができます。コントロールビットにはなんの操作もしていないのに、まるで２番目の量子ビットから位相が乗り移ったかのように観測確率が変化するのは不思議ですね。ただし、この方法で固有値の位相を知ろうとすると何度も測定しなければなりません。この問題点を解決するアルゴリズムとして位相推定アルゴリズムが出てくることになります。

量子フーリエ変換

量子フーリエ変換は離散フーリエ変換を行うアルゴリズムです。数列の長さを$2^n$とすると、古典計算機ではみんな大好き高速フーリエ変換で離散フーリエ変換を$O(n2^n)$で行うことができますが、量子フーリエ変換では$O(n^2)$という$n$の多項式時間で解くことができます！

まず、離散フーリエ変換の定義式を思い出してみましょう。
$$
y_k = \frac{1}{\sqrt{2^n}} \sum_{j=0}^{2^n-1} x_j e^{i\frac{2\pi kj}{2^n}}
$$
これをqubit上で表現すると、
$$
|x\rangle = \sum_{j=0}^{2^n-1} x_j |j\rangle \rightarrow |y\rangle = \sum_{k = 0}^{2^n-1} y_k |k\rangle \\
|j\rangle \rightarrow \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} e^{i\frac{2\pi jk}{2^n}}|k\rangle
$$
となります。

実はこの変換はunitaryであり、量子コンピュータはユニバーサルゲートセットによって任意のユニタリ変換を近似できるので、この変換は量子コンピュータ上で実現可能です。しかし、このままではどのようなゲートを噛ませればいいのかわからないのでもっと使いやすい形に変形したいと思います。kを二進数で展開していろいろ頑張ると以下のようになります。
$$
|j\rangle \rightarrow\otimes_{l = 1}^{2^n} (|0\rangle+e^{\frac{2\pi i}{2^l}}|1\rangle)
$$
この形にすると、それぞれのq-bitのテンソル積としてかけているので、フーリエ変換がアダマールゲートと回転ゲートにより、実現できることがわかると思います（下図）。

実はショアのアルゴリズムではその逆変換である逆量子フーリエ変換を使うのでこれも下に示しておきます。これはただゲートを反対からかけていっただけです(ただし回転は逆回転）。

ところで量子ビットの読む方向には注意が必要です。ここではq0q1q2q3の順に読むことにします。つまり、$q0 = 1, q1 = q2 = q3 = 0$なら$x = 1000(2) = 8$です。
一応ゲートの説明もしておくと、Hはアダマールゲート、U1はRzゲートと同じで、$|1\rangle$にのみ$e^{i\theta}$の位相をつけるゲートです。詳しくはこちら

位相推定アルゴリズム

位相アルゴリズムはアダマールテストを改良したようなアルゴリズムとなっています。

サイズ$2^n$のユニタリ行列$U$の固有値を$e^{2\pi i\lambda}$, それに対応する固有ベクトルを$|\varphi\rangle$とします。もし、$\lambda$が二進法で$m$桁の小数で表すことができ、
$$
|y\rangle = \frac{1}{\sqrt{2^n}}\sum_{k = 0}^{2^m-1}e^{2\pi ij\lambda}|j\rangle
$$
という状態を作り出すことができれば、逆離散フーリエ展開で$\lambda$がわかります（前章の定義式と見比べてみてください）。

位相推定アルゴリズムにはm+n個の量子ビットを使い、入力には$|0\rangle|\varphi\rangle$を入れることにします。まずアダマールゲートを最初のm個の量子ビットに適用することで$|0\rangle$から$|2^m-1\rangle$までの重ね合わせ状態を作ります。次に$i$の$k$桁目が1ならば後ろの$|\varphi\rangle$に$U^k$をかけます。するとアダマールテストと同じ原理で$|i\rangle\rightarrow e^{2\pi ik\lambda}|i\rangle$になります。これを$k = 1, 2, \cdots m$で繰り返すことで、上位mビットに上の状態を実現することができます！（jを二進数で考えてみてください）回路図は下のようになります。回路図を見た方がわかりやすいかもしれません。

素因数分解への応用

なぜ、位数rが求まると、素因数分解できるのでしょうか？まず、rは偶数であると仮定します（ランダムに取るといい確率で偶数になってくれるらしいです）。すると以下のような変形ができます。
$$
x^r = 1 \mod N\
(x^{r/2}-1)(x^{r/2}+1) = 0\mod N
$$
つまり、$(x^{r/2}-1)$と$(x^{r/2}+1)$のどちらかはNと非自明な公約数を持つことになります。どちらもがN　の倍数となってしまう確率は低いです。公約数はユークリッドの互助法によって、古典計算機で高速に計算できます。そこで$U^r = I$なる行列を作りたいのですが、これは簡単で、$U|i\rangle = |i\cdot x \mod N\rangle$となるように定義すれば良いです。ただし$i\geq N$については何もしないことにします。実はこれもユニタリ行列であることが示せます。これに位相推定アルゴリズムを用いれば位数rが十分な確率で求まります。一つ問題があるのは、位相推定アルゴリズムでは、固有ベクトルを利用していましたが、今回はわかりません。しかし、入力としては$\left|1\right>$を用いれば十分です。なぜなら、固有ベクトルで展開すれば、位相推定アルゴリズムで得られる値は固有値のどれかであり、そのどれであっても$s/r (s = 0, 1, \dots r-1)$という形で表せるため、連分数アルゴリズムによってrを高速に求められるからです($s = 0$の場合は無視します)。連分数アルゴリズムは小数を連分数で近似することでもっともらしい分数としての形を推定します。詳しいことはcontinued fraction algorithmでググりましょう。一応実装は後に載せてあります。

Qiskitでの実装

さて、ここからが本番です。上で見てきたアルゴリズムをシュミレーターで実行して$N = 15$を$x = 4$で因数分解してみましょう。今回はqiskitという量子計算ライブラリを使います。qiskitについて詳しくは公式ドキュメントを参照してください。実行は全てjupyter notebook上で行なっています。

量子コンピュター部分は以下のようになります。前章で出てきたUの実装ですが、一般の剰余演算を実装するのは困難なので、$N = 15, x = 4$という事実を活用してしまっています（${U}^{2^i}$については位数が2であることも使ってしまってい省略しています...）が、まあ今回は雰囲気を掴むだけなのでいいでしょう。

from qiskit import *
from qiskit.providers.ibmq import least_busy
from qiskit.visualization import plot_histogram
from qiskit.tools.monitor import job_monitor
import math

q = QuantumRegister(8, "q")
c = ClassicalRegister(4, "c")
circuit = QuantumCircuit(q, c)
circuit.x(7)
#Hadamard transform
for i in range(4):
    circuit.h(i)
#U4^1 gate
circuit.cswap(3, 4, 6)
circuit.cswap(3, 5, 7)
circuit.barrier()
#U4^2 gate
#U4^4 gate
#U4^8 gate
circuit.barrier()
#inverse quantum Fourier transform
circuit.swap(0, 3)
circuit.swap(1, 2)
circuit.h(3)
circuit.cu1(-np.pi/2, 3, 2)
circuit.h(2)
circuit.barrier()
circuit.cu1(-np.pi/4, 3, 1)
circuit.cu1(-np.pi/2, 2, 1)
circuit.h(1)
circuit.barrier()
circuit.cu1(-np.pi/8, 3, 0)
circuit.cu1(-np.pi/4, 2, 0)
circuit.cu1(-np.pi/2, 1, 0)
circuit.h(0)
circuit.barrier()
#measure
for i in range(4):
    circuit.measure(q[i], c[3-i])#二進法でx = q1q2q3q4になるようにしている
circuit.draw(output='mpl')

さて、あとは古典計算の部分です。連分数展開は以下のように実装しました（結構適当に書いてちょっとサンプル試しただけなのでかなり怪しいですが...）。相対誤差がeps以下になったら終了するようにしています。

eps = 0.01
def continued_fractions_algorithm(x):
    res = [int(x)]
    x-=int(x)
    if x/(res[0]+0.1)>eps:
        a = continued_fractions_algorithm(1/x)
        res+=a
    return res

さて、ここで（量子計算機としての）計算量を見積もりましょう。素因数分解したい自然数が$n$bitであるとして、離散フーリエ変換パートが $O(n^2)$、位相推定アルゴリズムパートが$O(n^3)$です。古典パートはボトルネックとならないので、全体で$O(n^3)$の多項式時間で解くことができました！
あとは以下のコードを実行すると、

def shor_algorithm(use_simulator):
    if use_simulator:
        backend = BasicAer.get_backend('qasm_simulator')
    else:
        backend = IBMQ.get_provider().get_backend('ibmq_16_melbourne')
    flag = False
    job = execute(circuit, backend, shots = N)
    job_monitor(job)
    result = job.result()
    counts = result.get_counts(circuit)
    measures = np.array(list(map(lambda x:int(x, 2), counts.keys())), dtype = np.int64)
    probs = np.array(list(counts.values()))/N
    for i in range(5):
        output = np.random.choice(measures, p = probs)
        a = continued_fractions_algorithm(output/2**4)
        r , s =1,  a[-1]
        for i in range(len(a)-1)[::-1]:
            r, s = s, a[i]*s+r 
        if r % 15 == 0 or s == 0:
            continue
        d = math.gcd(15, 4**(r-1)-1)
        if d != 1:
            flag = True
            break
    plot_histogram(result.get_counts())
    if flag:
        print('{0} devides 15 with r = {1}'.format(d, r))
    else:
        print('the algorithm failed')    
    return result

%matplotlib inline
use_simulator = True
result = shor_algorithm(use_simulator)
plot_histogram(result.get_counts())

となり、因数分解に成功しました！

実機での計算

qiskitから無料で使える量子コンピューターであるIBMQを操作することができます。公式ページはこちら。操作にはアカウント作成が必要です。

まず以下のコードを実行しておきます。

from qiskit import IBMQ
my_token = ""
IBMQ.save_account(my_token)
provider = IBMQ.load_account()

my_tokenのところには公式ページのMy Accountのところから取得したトークンを入れてください。

今回は8qubit使っているのでibmq_16_melbourneを使います。少し時間がかかりますが、上のコードでuse_simulator = Falseとして実行すると、

あれっ？

まとめ

いかがでしたか？今回量子コンピュータによる因数分解を試したわけですが、前の章でみてきた通り実機だとかなり誤差があります。もちろん今回は誰でも使えるもの簡易的なもので試しているわけですが、15の因数分解にすら手こずるとなると量子コンピュータが実用化されるには時間がかかりそうですね...。

参考

Quantum Native Dojo https://dojo.qulacs.org/ja/latest/
宮野健二郎,古澤明, 量子コンピュータ入門　日本評論社 2016
Michael A. Nilsen,Isaac L. Chuang Quantum Computation and Quantum Information 10th Anniversary Edition 2010