「DiscordBotを作ってみたい、だけどPCを持っていない」という人向けにAndroid端末1つでDiscordBotを作って実際に動作させるまでの方法を紹介します。
なるべく初心者向けに説明していますが、Pythonやdiscord.py, Linuxコマンド, Vimコマンドなどの詳しい説明は省略しています。

動作環境

• Android 9
• Termux 0.90
• Python 3.8.1
• discord.py 1.3.1

1. Termuxのインストール

TermuxはGoogle Playストアからインストールできます。
https://play.google.com/store/apps/details?id=com.termux

Termuxとはroot化不要でLinux環境を用意することができるエミュレータです。
ここではPythonとVimを使ってDiscordBotを書く方法を説明します。

2. Termuxでの準備

まずターミナルで以下のコマンドを実行し、aptの更新とPython, Vimのインストールを行います。($はコマンドラインであることを表すもので入力は不要です)

$ apt update
$ apt upgrade
$ apt install python
$ apt install vim

次にPythonで簡単にDiscordBotを動かすことができるdiscord.pyというライブラリをインストールします。

$ pip install discord.py

※もしインストールが上手くいかない場合はpip install yarlで先にyarlをインストールしておくことで解決する場合があります。

次に実行するPythonファイルを置くディレクトリを作り、移動します。
※{}で囲まれた部分は説明してる内容で置き換えてください

$ mkdir {好きなディレクトリ名}
$ cd {上で作ったディレクトリ名}

3. Pythonファイルの作成

$ vim {好きなファイル名}.py

上のコマンドを実行すると新規ファイルがVimが開かれるので、キーボードでiと入力し、INSERTモードにします。
ここで下のコードをコピペしてください。(後述のサイトにあるものを一部改変し拝借させていただきました)

import discord

client = discord.Client()

@client.event
async def on_ready():
    print('Logged in as')
    print(client.user.name)
    print(client.user.id)
    print('------')

@client.event
async def on_message(message):
    # 送り主がBotだった場合反応したくないので
    if message.author.bot:
        return

    # 「おはよう」で始まるか調べる
    if message.content.startswith("おはよう"):
        # メッセージを書きます
        m = "おはようございます" + message.author.name + "さん！"
        # メッセージが送られてきたチャンネルへメッセージを送ります
        await message.channel.send(m)

client.run("token")

4. Botアカウントを作成する

Chromeなどのブラウザで「https://qiita.com/PinappleHunter/items/af4ccdbb04727437477f 」のサイトのサイトを開き、「Bot用のトークンを手に入れる」の項に従ってBotアカウントを作成して使用したいサーバーにBotを追加しておきます。
ここで「後で使います」と言われているトークンをコピーしておいてください。
※このトークンはBotを動かすためのパスワードのようなものなので扱いには十分気を付けてください。

5. Pythonファイルを実行する

Termuxを開き、cliant.run("token")のtokenの部分に、先程コピーしてきたトークンをペーストし、ESC→:wqでファイルを保存しVimを閉じます。

下のコマンドでPythonファイルを実行します。

$ python {先程作成したファイル名}.py

これでLogged in as...のように表示されればBotが動作していることになります。

Botを追加したサーバーでおはようなどを送信して返答が来るか試してみてください。
終了したいときはTermuxを終了すればBotも動かなくなります。もしBotを常時動かしたい場合は、Herokuなどを利用することでできます。
ちなみにTermuxではGitも使えるので、AndroidのみでHerokuにデプロイすることも可能です。

参考サイト

Pythonで簡単なDiscord Botの作り方
Android端末上でDiscordのbotを動かす(Termux)
AndroidでTermuxを使いroot化なしでLinux環境を作る！

こんにちは、こんばんは、おはようございます。
K.S.ロジャース の島袋です。

以前PythonのWebフレームワークの次期覇権はコレ みたいな記事を書きましたが、それからしばらくPythonから離れたエンジニア生活を過ごしていたところ、なんと FastAPI なるフレームワークが登場しているではないですか！

正直、何番煎じかわからない遅れ方してますが、ぜひとも紹介させてください。

FastAPIの特徴

公式トップページ に書いてあるには、

1. 高速 - NodeJSやGoなみに早い
2. コード書くの早い - 2倍3倍で書ける！
3. バグが減る - 40％くらい減る！
4. 直感的 - デバッグが簡単
5. 簡単 - 使いやすい。ドキュメント読むのも簡単！
6. 短い - 重複が少ない。パラメータ宣言で複数の機能を呼べる
7. 堅牢 - 自動生成あるからほんと楽
8. スタンダード - OpenAPI標準だから互換性あるよ（Swaggerと互換性あるよ）

とありますが一番下の OpenAPI標準だから互換性あるよ（Swaggerと互換性あるよ） これが素晴らしすぎる。
FastAPI、実はSwagger（ドキュメント）の自動生成ができるのです！

ドキュメント自動生成

pipで fastapi と uvicorn をインストールしたら早速触ってみます。
ちなみに標準のpipには環境切り分けなどがないので、pipenv を利用することをおすすめします。

from fastapi import FastAPI

app = FastAPI()


@app.get("/")
async def root():
    return {"message": "Hello World"}

上記を main.py で作成後に uvicorn main:app --reload で起動。
http://127.0.0.1:8000/docs を確認すると…

はい！素敵！！
試しに path parameter を設定すると…

Path parameters

@app.get("/{name}")
async def root():
    return {"message": "Hello " + name}

自動で反映されます。
ドキュメントが自動で作られていくとか最強です。
とくにドキュメントを重要視している弊社では便利すぎる機能です。

Query Parameters

@app.get("/")
async def read_item(name: str = ""):
    return {"message": "Hello " + name}

Requests

Modelを使います。

from fastapi import FastAPI
from pydantic import BaseModel

class Item(BaseModel):
    name: str
    description: str = None

app = FastAPI()


@app.post("/items/")
async def create_item(item: Item):
    return item

さらにModelクラス内でのプロパティにバリデーションを指定できるので、取り回しが凄く良いです。
ここらへんは Django のModelと似たような構成ですね。

class Item(BaseModel):
    name: str = Field(None, min_length=2, max_length=5)
    description: str = Field(None, max_length=100)

CRUD (Crate, Read, Update, Delete)

基本のCRUDをざっくり書いてみます。
一部抜粋して書いているので、詳しくは書きリンクをご参考ください。

https://github.com/shuto-S/sample-fast-api

DB Connection & Migration

DBへの接続とマイグレーションは databases と sqlalchemy を使います。
今回は取り回しのよい SQLite をDBに使います。

さくっとコネクション周りを書いて、

import databases
import sqlalchemy

DATABASE_URL = "sqlite:///./db.sqlite3"

database = databases.Database(DATABASE_URL)
engine = sqlalchemy.create_engine(DATABASE_URL, echo=False)
metadata = sqlalchemy.MetaData()

Schemas

テーブル定義を記述して、同ファイルを実行しておきます。

import sqlalchemy
from db import metadata, engine


items = sqlalchemy.Table(
    "items",
    metadata,
    sqlalchemy.Column("id", sqlalchemy.Integer, primary_key=True),
    sqlalchemy.Column("name", sqlalchemy.String),
    sqlalchemy.Column("description", sqlalchemy.String),
)

metadata.create_all(engine)

Main

main.py も改修

from fastapi import FastAPI
from db import database
from starlette.requests import Request
from routers import items

app = FastAPI()

@app.on_event("startup")
async def startup():
    # DBコネクション開始
    await database.connect()

@app.on_event("shutdown")
async def shutdown():
    # DBコネクション切断
    await database.disconnect()

# routersを登録
app.include_router(items.router)

# ミドルウェアでDBコネクション埋め込んでおく（routerで取得できるように）
@app.middleware("http")
async def db_session_middleware(request: Request, call_next):
    request.state.connection = database
    response = await call_next(request)
    return response

Routers

router もファイル分割してエンドポイント毎に記述できるようにします。

from fastapi import APIRouter, Depends, HTTPException
from typing import List
from databases import Database
from starlette.status import HTTP_204_NO_CONTENT

from utils import get_db_connection
from schemas import items
from models.item import ItemModel


router = APIRouter()

@router.get("/items", tags=["items"], response_model=List[ItemModel])
async def list_item(database: Database = Depends(get_db_connection)):
    query = items.select()
    return await database.fetch_all(query)

@router.post("/items", tags=["items"], response_model=ItemModel)
async def create_item(data: ItemModel, database: Database = Depends(get_db_connection)):
    query = items.insert()
    await database.execute(query, data.dict())
    return {**data.dict()}

@router.patch("/items/{item_id}", tags=["items"], response_model=ItemModel)
async def update_item(item_id: int, data: ItemModel, database: Database = Depends(get_db_connection)):
    query = items.update().where(items.columns.id==item_id)
    ret = await database.execute(query, data.dict())
    if not ret:
        raise HTTPException(status_code=404, detail="Not Found")
    return {**data.dict()}

@router.delete("/items/{item_id}", tags=["items"], status_code=HTTP_204_NO_CONTENT)
async def delete_item(item_id: int, database: Database = Depends(get_db_connection)):
    query = items.delete().where(items.columns.id==item_id)
    ret = await database.execute(query)
    if not ret:
        raise HTTPException(status_code=404, detail="Not Found")

まとめ

どうでしょう？ドキュメント自動生成はかなり魅力的な機能ではないでしょうか？
Responder もかなり良いフレームワークだと思いますが、FastAPI のほうが癖が少なく書きやすく感じました。

あとづけ

ちなみに弊社、Tech系以外にも会社ブログも掲載してますので、気になった方は是非どうぞ。
https://www.wantedly.com/companies/ks-rogers

1

class Person(object):
    def talk(self):
        print('talk')
    def run(self):
        print('person run')

class Car(object):
    def run(self):
        print('car run')

class PersonCarRobot(Person, Car):
    def fly(self):
        print('fly')

person_car_robot = PersonCarRobot()
person_car_robot.talk()
person_car_robot.run()
person_car_robot.fly()

1の実行結果

talk
person run
fly

PersonCarRobotクラスは、
PersonクラスとCarクラスの両方を継承している。
なので、
両方のクラスのメソッドを持っている。

ここで、
PersonクラスもCarクラスも runというメソッドを持っている。

この場合、
class PersonCarRobot(Person, Car)となっており、
先にPersonクラスを継承して、ないものをCarクラスから継承するイメージ。
なので、Carクラスのrunメソッドは継承されない。

class PersonCarRobot(Car, Person)と入れ替えると
先にCarクラスを継承して、ないものをPersonクラスから継承する。
Personクラスのrunメソッドは継承されない。

2

class Person(object):
    def talk(self):
        print('talk')
    def run(self):
        print('person run')

class Car(object):
    def run(self):
        print('car run')

class PersonCarRobot(Car, Person):
    def fly(self):
        print('fly')

person_car_robot = PersonCarRobot()
person_car_robot.talk()
person_car_robot.run()
person_car_robot.fly()

2の実行結果

talk
car run
fly

pandasとは

pandasとは、Pythonのライブラリです。
pandasは、主にデータを操作する時に使われます。AIや機械学習にpandasは必須です。
ですが、Excelのファイルを操作する機能も搭載されています。
今回は、pandasを使ったExcelファイルの読み込みの仕方について説明していきます。

Excelファイル読み込み

最初のシートのみの読み込み

sample.xlsxの中身が

だった場合

main.py

import pandas as pd
df = pd.read_excel('sample.xlsx')
print(df)

上記のmain.pyを実行すると

となり、無事Excelファイルを読み込む事に成功しました。

全シートから読み込み

因みに、main.pyの二行目を

df = pd.read_excel('sample.xlsx', sheet_name =None)

とすれば、sample.xlsxの全てのシートからデータを取ってくる事ができます。

指定したシート読み込み

番号で指定

sheet_nameをNoneにすると、全シート読み込みとなりました。
sheet_nameに0から始まる番号を指定することで、[指定した数字+1]番目のシートを読み込むことができます。

df = pd.read_excel('sample.xlsx', sheet_name =0) #1枚目のシート読み込み
df = pd.read_excel('sample.xlsx', sheet_name =2) #3枚目のシート読み込み

名前で指定

sheet_nameに文字列を指定することで、その文字列をシート名に持つシートを読み込むことができます。

df = pd.read_excel('sample.xlsx', sheet_name ='sheet1') #名前が'sheet1'のシート読み込み

まとめ

pandasを使ってExcelファイルを読み込む方法を紹介しました。
この機能はpandasの全機能の氷山の一角に過ぎません。

非常に様々な場面で使われるライブラリなので、気になった方は是非詳しいサイトなどを調べてみてください。

動画を連結したかったので作ってみた．

シチュエーション

複数の動画が入ったディレクトリあり，その中の動画を連結した一つの動画を作成するプログラム．

ソース

import cv2
import glob

def combine_movie(movie_files,out_path):

    # 形式はmp4
    fourcc = cv2.VideoWriter_fourcc('m','p','4','v')

    #動画情報の取得
    movie = cv2.VideoCapture(movie_files[0])
    fps    = movie.get(cv2.CAP_PROP_FPS)
    height = movie.get(cv2.CAP_PROP_FRAME_HEIGHT)
    width  = movie.get(cv2.CAP_PROP_FRAME_WIDTH)


    # 出力先のファイルを開く
    out = cv2.VideoWriter(out_path, int(fourcc), fps, (int(width), int(height)))


    for movies in (movie_files):
        print(movies)
        movie = cv2.VideoCapture(movies)

        if movie.isOpened() == True:
            ret,frame = movie.read()
        else:
            ret = False

        while ret:
            cv2.putText(frame,' ',(400,400),cv2.FONT_HERSHEY_SIMPLEX,0.75,(0,0,255),2)

            # 読み込んだフレームを書き込み
            out.write(frame)

            # 次のフレーム読み込み
            ret,frame = movie.read()

# ディレクトリ内の動画をリストで取り出す
files = sorted(glob.glob("./movie_dir/*.mp4"))

# 出力ファイル名
out_path = "movie_out1.mp4"

combine_movie(files,out_path)

強化学習結果などの沢山の動画を一発で連結したかったので作ってみた．

シチュエーション

複数の動画が入ったディレクトリあり，その中の動画を連結した一つの動画を作成するプログラム．

ソース

import cv2
import glob

def comb_movie(movie_files,out_path):

    # 形式はmp4
    fourcc = cv2.VideoWriter_fourcc('m','p','4','v')

    #動画情報の取得
    movie = cv2.VideoCapture(movie_files[0])
    fps    = movie.get(cv2.CAP_PROP_FPS)
    height = movie.get(cv2.CAP_PROP_FRAME_HEIGHT)
    width  = movie.get(cv2.CAP_PROP_FRAME_WIDTH)


    # 出力先のファイルを開く
    out = cv2.VideoWriter(out_path, int(fourcc), fps, (int(width), int(height)))


    for movies in (movie_files):
        print(movies)
        movie = cv2.VideoCapture(movies)

        if movie.isOpened() == True:
            ret,frame = movie.read()
        else:
            ret = False

        while ret:
            cv2.putText(frame,' ',(400,400),cv2.FONT_HERSHEY_SIMPLEX,0.75,(0,0,255),2)

            # 読み込んだフレームを書き込み
            out.write(frame)

            # 次のフレーム読み込み
            ret,frame = movie.read()

# ディレクトリ内の動画をリストで取り出す
files = sorted(glob.glob("./movie_dir/*.mp4"))

# 出力ファイル名
out_path = "movie_out1.mp4"

comb_movie(files,out_path)

12.1 コードのテスト

12.1.1 pylint、pyflakes、pep8によるチェック

style1

a=1
b=2
print(a)
print(b)
print(c)

実行結果

$ pylint style1.py
************* Module style1
style1.py:1:1: C0326: Exactly one space required around assignment
a=1
 ^ (bad-whitespace)
style1.py:2:1: C0326: Exactly one space required around assignment
b=2
 ^ (bad-whitespace)
style1.py:1:0: C0114: Missing module docstring (missing-module-docstring)
style1.py:1:0: C0103: Constant name "a" doesn't conform to UPPER_CASE naming style (invalid-name)
style1.py:2:0: C0103: Constant name "b" doesn't conform to UPPER_CASE naming style (invalid-name)
#先頭がEとなっているのはエラーを示す。
style1.py:5:6: E0602: Undefined variable 'c' (undefined-variable)

-------------------------------------
Your code has been rated at -10.00/10

エラー解消のため、cの値を追加する。

style2.py

a=1
b=2
c=3
print(a)
print(b)
print(c)

実行結果

$ pylint style2.py
************* Module style2
style2.py:1:1: C0326: Exactly one space required around assignment
a=1
 ^ (bad-whitespace)
style2.py:2:1: C0326: Exactly one space required around assignment
b=2
 ^ (bad-whitespace)
style2.py:3:1: C0326: Exactly one space required around assignment
c=3
 ^ (bad-whitespace)
style2.py:1:0: C0114: Missing module docstring (missing-module-docstring)
style2.py:1:0: C0103: Constant name "a" doesn't conform to UPPER_CASE naming style (invalid-name)
style2.py:2:0: C0103: Constant name "b" doesn't conform to UPPER_CASE naming style (invalid-name)
style2.py:3:0: C0103: Constant name "c" doesn't conform to UPPER_CASE naming style (invalid-name)

------------------------------------
Your code has been rated at -1.67/10

変数名を長くしてみる。

style3.py

def func():

    first=1
    second=2
    third=3
    print(first)
    print(second)
    print(third)

func()

実行結果

$ pylist style3.py
-bash: pylist: command not found
uemuratntonoAir:bin uemura$ pylint style3.py
************* Module style3
style3.py:3:9: C0326: Exactly one space required around assignment
    first=1
         ^ (bad-whitespace)
style3.py:4:10: C0326: Exactly one space required around assignment
    second=2
          ^ (bad-whitespace)
style3.py:5:9: C0326: Exactly one space required around assignment
    third=3
         ^ (bad-whitespace)
style3.py:1:0: C0114: Missing module docstring (missing-module-docstring)
style3.py:1:0: C0116: Missing function or method docstring (missing-function-docstring)

-----------------------------------
Your code has been rated at 3.75/10

12.1.2 unittestによるテスト

• コードをソース管理システムにコミットする前に、まずテストプログラムを書いて、コードが全てのテストに合格することを確認する作業はグッドプラクティス。

cap.py

def just_do_it(text):
     return text.capitalize()

test_cap.py

import unittest
import cap

class TestCap(unittest.TestCase):

    def setUp(self):
        pass

    def tearDown(self):
        pass

    def test_one_word(self):
        text="duck"
        result=cap.just_do_it(text)
        self.assertEqual(result, "Duck")

    def test_multiple_words(self):
        text="a veritable flock of ducks"
        result=cap.just_do_it(text)
        self.assertEqual(result, "A Veritable Flock Of Ducks")

#assertで始まる名前を持つメソッドで結果をチェックする。

if __name__=="__main__":
    unittest.main()

実行結果

$ python test_cap.py
..
----------------------------------------------------------------------
Ran 2 tests in 0.000s

OK
$ python test_cap.py
..
----------------------------------------------------------------------
Ran 2 tests in 0.000s

OK
$ python test_cap.py
..F
======================================================================
FAIL: test_words_with_apostrophes (__main__.TestCap)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "test_cap.py", line 25, in test_words_with_apostrophes
    self.assertEqual(result, "I'm Fresh Out Of Ideas")
AssertionError: "I'M Fresh Out Of Ideas" != "I'm Fresh Out Of Ideas"

#^は実際に文字列が異なっている箇所を示している。
- I'M Fresh Out Of Ideas
?   ^
+ I'm Fresh Out Of Ideas
?   ^


----------------------------------------------------------------------
Ran 3 tests in 0.001s

FAILED (failures=1)

capitalize()は先頭の単語のみタイトルケースにする。

cap.py

#capitalize()をtitle()に書き換え
def just_do_it(text):
     return text.title()

実行結果

$ python test_cap.py
...
----------------------------------------------------------------------
Ran 3 tests in 0.000s

OK

test_cap.py

import unittest
import cap

class TestCap(unittest.TestCase):

    def setUp(self):
        pass

    def tearDown(self):
        pass

    def test_one_word(self):
        text="duck"
        result=cap.just_do_it(text)
        self.assertEqual(result, "Duck")

    def test_multiple_words(self):
        text="a veritable flock of ducks"
        result=cap.just_do_it(text)
        self.assertEqual(result, "A Veritable Flock Of Ducks")

#assertで始まる名前を持つメソッドで結果をチェックする。

    def test_words_with_apostrophes(self):
        text="I'm fresh out of ideas"
        result=cap.just_do_it(text)
        self.assertEqual(result, "I'm Fresh Out Of Ideas")

    def test_words_with_quotes(self):
        text="\"You're despicable,\" said Daffy Duck"
        result=cap.just_do_it(text)
        self.assertEqual(result, "\"You're Despicable,\" Said Daffy Duck")


if __name__=="__main__":
    unittest.main()

実行結果

$ python test_cap.py
...F
======================================================================
FAIL: test_words_with_quotes (__main__.TestCap)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "test_cap.py", line 30, in test_words_with_quotes
    self.assertEqual(result, "\"You're Despicable,\" Said Daffy Duck")
AssertionError: '"you\'re Despicable," Said Daffy Duck' != '"You\'re Despicable," Said Daffy Duck'
- "you're Despicable," Said Daffy Duck
?  ^
+ "You're Despicable," Said Daffy Duck
?  ^


----------------------------------------------------------------------
Ran 4 tests in 0.002s

FAILED (failures=1)

12.1.3 noseによるテスト

• noseはサードパーティパッケージ。

test_cap_nose.py

import cap
from nose.tools import eq_

    def test_one_word(self):
        text="duck"
        result=cap.just_do_it(text)
        eq_(result, "Duck")

    def test_multiple_words(self):
        text="a veritable flock of ducks"
        result=cap.just_do_it(text)
        eq_(result, "A Veritable Flock Of Ducks")

    def test_words_with_apostrophes(self):
        text="I'm fresh out of ideas"
        result=cap.just_do_it(text)
        eq_(result, "I'm Fresh Out Of Ideas")

    def test_words_with_quotes(self):
        text="\"You're despicable,\" said Daffy Duck"
        result=cap.just_do_it(text)
        eq_(result, "\"You're Despicable,\" Said Daffy Duck")

unittestnでテストしたときに見つけたのと同じバグ。

実行結果

$ nosetests test_cap_nose.py
E
======================================================================
ERROR: Failure: IndentationError (unexpected indent (test_cap_nose.py, line 4))
----------------------------------------------------------------------
...(省略)
    def test_one_word(self):
    ^
IndentationError: unexpected indent

----------------------------------------------------------------------
Ran 1 test in 0.001s

FAILED (errors=1)

12.2 Pythonコードのデバッグ

test_dump.py

from dump1 import dump

@dump
def double(*args, **kwargs):
    "Double every argument"
    output_list=[2*x for x in args]
    output_dict={k:2*v for k,v in kwargs.items()}
    return output_list, output_dict

if __name__=="__main__":
    output=double(3, 5, first=100,next=98.6, last=-40)

実行結果

$ python test_dump.py
Function name: double
Input arguments: 3 5
Input keyword arguments: dict_items([('first', 100), ('next', 98.6), ('last', -40)])
Output: ([6, 10], {'first': 200, 'next': 197.2, 'last': -80})

12.3 pdbによるデバッグ

capitals.py

def process_cities(filename):
    with open(filename, "rt") as file:
        for line in file:
            line=line.strip()
            if "quit" in line.lower():
                return
            country, city=line.split(",")
            city=city.strip()
            country=country.strip()
            print(city.title(),country.title(), sep=",")

if __name__=="__main__":
    import sys
    process_cities(sys.argv[1])

実行結果

$ python capitals.py cities1.csv
Paris,France
Caracas,Venuzuela
Vilnius,Lithunia

cities2.csvを実行中、15行中わずか5行出力したところで停止。

cities2.csv

$ python capitals.py cities2.csv
Buenos Aires,Argentinam
La Paz,Bolivia
Brasilia,Brazil
Santiago,Chile
Bogota,Colombia

実行結果

$ python -m pdb capitals.py cities2.csv

#cと入力すると正常終了もしくはエラーで止まる。
(Pdb) c
Buenos Aires,Argentinam
La Paz,Bolivia
Brasilia,Brazil
Santiago,Chile
Bogota,Colombia
The program finished and will be restarted
> /Users/practice/bin/capitals.py(1)<module>()
-> def process_cities(filename):

#sと入力すると行を順に実行する。
(Pdb) s
> /Users/practice/bin/capitals.py(13)<module>()
-> import sys
(Pdb) s
> /Users/practice/bin/capitals.py(14)<module>()
-> process_cities(sys.argv[1])
(Pdb) s
--Call--
> /Users/practice/bin/capitals.py(1)process_cities()
-> def process_cities(filename):
(Pdb) s
> /Users/practice/bin/capitals.py(2)process_cities()
-> with open(filename, "rt") as file:

#第6行にブレークポイントセットした。
(Pdb) b 6
Breakpoint 1 at /Users/uemura/practice/bin/capitals.py:6

#ブレークポイントに当たるか、入力行を全て読み出して終了するまでプログラム実行する。
(Pdb) c
Buenos Aires,Argentinam
La Paz,Bolivia
Brasilia,Brazil
Santiago,Chile
Bogota,Colombia
> /Users/practice/bin/capitals.py(6)process_cities()
-> return
#何を読み出したのか確認する。
(Pdb) p line
'ecuador,quito'
(Pdb) b
Num Type         Disp Enb   Where
1   breakpoint   keep yes   at /Users/uemura/practice/bin/capitals.py:6
    breakpoint already hit 1 time
#1はコード行、現在位置(->)、ブレークポイント(B)を示す。
#オプションなしのlは前回のlの表示の末尾から表示を始めるのでオプションで指定する。
(Pdb) l 1
  1     def process_cities(filename):
  2         with open(filename, "rt") as file:
  3             for line in file:
  4                 line=line.strip()
  5                 if "quit" in line.lower():
  6 B->                 return
  7                 country, city=line.split(",")
  8                 city=city.strip()
  9                 country=country.strip()
 10                 print(city.title(),country.title(), sep=",")
 11     

quitテストが行全体でquitになっているものだけマッチするように書き換える。

capitals2.py

def process_cities(filename):
    with open(filename, "rt") as file:
        for line in file:
            line=line.strip()
            if "quit" ==line.lower():
                return
            country, city=line.split(",")
            city=city.strip()
            country=country.strip()
            print(city.title(),country.title(), sep=",")

if __name__=="__main__":
    import sys
    process_cities(sys.argv[1])

実行結果

$ python capitals2.py cities2.csv
Buenos Aires,Argentinam
La Paz,Bolivia
Brasilia,Brazil
Santiago,Chile
Bogota,Colombia
Quito,Ecuador
Stanley,Falkland Islands
Cayenne,French Guiana
Asuncion,Paraguay
Lima,Peru
Paramaribo,Suriname
Montevideo,Uruguay
Caracas,Venezuela

12.4 エラーメッセージのロギング

• ロギングのPython標準ライブラリはloggingモジュール。
  • ログに保存したいメッセージ
  • ランク付けのための優先順位レベルとそれに対応する関数
  • モジュールとの主要な通信経路となる一つ以上のロガーオブジェクト
  • メッセージを端末、ファイル等そのほかの場所に送るハンドラ
  • 出力を作成するフォーマッタ
  • 入力に基づいて判断を下すフィルタ

>>> import logging
>>> logging.debug("Looks like rain")
>>> logging.info("And hill")
#デフォルトの優先順位レベルはWARNIGである。
>>> logging.warn("Did I hear thunder?")
__main__:1: DeprecationWarning: The 'warn' function is deprecated, use 'warning' instead
WARNING:root:Did I hear thunder?
>>> logging.critical("Stop fencing and get inside")
CRITICAL:root:Stop fencing and get inside

#デフォルトレベルはbasicConfig()で設定できる。
#DEBUG以上の全てのレベルがログに書かれる。
>>> import logging
>>> logging.basicConfig(level="DEBUG")
>>> logger.debug("Timber")
DEBUG:bunyan:Timber
>>> logging.info("And hill")
INFO:root:And hill

12.5 コードの最適化

12.5.1 実行時間の計測

time1.py

#現在の時刻から処理終了時刻をひく。
from time import time

t1=time()
num=5
num*=2
print(time()-t1)

$ python time1.py
3.0994415283203125e-06
$ python time1.py
1.9073486328125e-06
$ python time1.py
1.6689300537109375e-06

time2.py

from time import time, sleep

t1=time()
sleep(1.0)
print(time()-t1)

実行結果

$ python time2.py
1.0051441192626953
uemuratntonoAir:bin uemura$ python time2.py
1.0001447200775146

timeitモジュールのtimeit()関数を使うと便利。

timeit1.py

#time.timeit("code",number,count)の構文
#timeitの内部で実行するためcodeは""で囲む必要がある。

from timeit import timeit
print(timeit("num=5;num*=2", number=1))

実行結果

$ python timeit1.py
1.63200000000141e-06
uemuratntonoAir:bin uemura$ python timeit1.py
1.252999999999671e-06

time2.py

#repeat()関数のrepeat引数を使うことで実行回数が増やせる。
from timeit import repeat
print(repeat("num=5;num*=2", number=1,repeat=3))

実行結果

$ python timeit2.py
[8.809999999977169e-07, 2.57000000000035e-07, 1.659999999993611e-07]

12.5.2 アルゴリズムとデータ構造

リスト内包表記とforループのどちらが早いか。

time_lists.py

from timeit import timeit

def make1():
     result=[]
     for value in range(1000):
         result.append(value)
     return result

def make2():
    result=[value for value in range(1000)]
    return result

print("make1 takes", timeit(make1, number=1000), "seconds")
print("make2 takes", timeit(make2, number=1000), "seconds")

リスト内包表記の方が高速である。

実行結果

$ python time_lists.py
make1 takes 0.07954489899999999 seconds
make2 takes 0.035908797000000006 seconds

12.6 Gitでのソース管理

#新しいディレクトリ作成後、移動する。
$ mkdir newdir
$ cd newdir

#カレントディレクトリのnewdirにローカルGitリポジトリを作る。
$ git init
Initialized empty Git repository in /Users/practice/bin/newdir/.git/

#Gitリポジトリにファイルを追加する。
$ git add test.py

#Gitの状態確認。
$ git status
On branch master

No commits yet

Changes to be committed:
  (use "git rm --cached <file>..." to unstage)

    new file:   test.py

#-m "my first commit"というコミットメッセージをつけてコミットする。
$ git commit -m "my first commit"
[master (root-commit) e03cd1c] my first commit
 1 file changed, 1 insertion(+)
 create mode 100644 test.py

#変更はコミット済み。
$ git status
On branch master
nothing to commit, working tree clean

感想

1周目ですが約1ヶ月で完走できた。ただ、アウトプットができていないのでアウトプットをメインに次の1ヶ月はやっていく。

参考文献

「Bill Lubanovic著 『入門 Python3』(オライリージャパン発行)」

クラスタリングの手法のひとつに、混合ガウスモデルの推定があります。これは、データセットを複数個のガウス分布で表現、及びクラスタ分けする手法です。通常、これらのガウス分布は有限個を仮定しますが、もし無限個あったとしたら、どんなクラスタリングになるのでしょうか。

本記事では、機械学習プロフェッショナルシリーズ『ノンパラメトリックベイズ　点過程と統計的機械学習の数理』を参考にして、「無限混合ガウスモデル」及び「ノンパラメトリック」の概要と実装を紹介します。この本、非常に分かりやすいのでおすすめです。

詳細な理論の話は省略し、概要とストーリーを重視したいと思います。
また、実装はソースコードを記載する程度に留めます。

※前回の記事
【Python】周辺化ギブスサンプリングを実装してみた
ただし、本記事だけ読んでも流れが伝わるように意識しています。

無限混合ガウスモデルの導入

通常の混合ガウスモデルでは、データセットが有限個のガウス分布で表現され、各データがその中のガウス分布と紐付いていると考えます。クラスラ数(=ガウス分布数)はあらかじめ決めておく必要がありますが、適切なクラスタ数はデータセットによって変化し、それを把握するのは容易ではありません。

もしクラスタリング中にクラスタ数も自動的に決まれば、上記の課題に悩まずに済みますよね。それを実現するのが、無限混合ガウスモデルです。

無限混合ガウスモデルでは、「ガウス分布は無限個存在し、データセットのデータは、そのうちの有限個と結びついている」という考え方に基づいています。
次の図は、後ほどお見せするクラスタリングの結果です。（めんどうなので図を使いまわします。ご了承ください。）

このデータセットは3個のガウス分布で表現されています。
通常の混合ガウスモデルでは、この世界にちょうど3個のガウス分布があると考えるのに対して、無限混合ガウスモデルでは、データセットと結びつかない潜在的なガウス分布が無限個あると考えます。
このように考えることで、データが既存の3個のクラスタに結びつく確率に加え、平均がどこにあるのか分からない未知のクラスタと結びつく確率も計算できるのです。

つまり、無限混合ガウスモデルによるクラスタリングでは、データが新しいクラスタに割り当てられたり、一方で既存のクラスタがどのデータとの結びつかず消滅したりと、クラスタ数を増減させながら、最後は適切な数に収束させることができます。

確率分布の算出

ここから計算式を掲載します。
通常の混合ガウス分布では、（$x_{1:N}$¹を観測しないときの）$z_i$の確率分布にカテゴリカル分布を適用します。

P(z_i = k | \pi) = \pi_k

ここで、$\pi := (\pi_1, \cdots, \pi_K)$は$(K-1)$次元単体上の点です。²
$\pi$の値の決め方として、一様分布を仮定したり、事前分布を仮定してベイズ推定する方法³があります。

一方、無限混合ガウスモデルでは、中華料理過程と呼ばれる確率過程を適用します。これは、$i$番目以外のデータのクラスタ$z_{1:N}^{\backslash i}$から、$i$番目のデータのクラスタ$z_i$の確率分布を算出する仕組みです。

P(z_i = k | z_{1:N}^{\backslash i}, \alpha) =
\left\{
\begin{array}{ll}
\frac{n_k}{N - 1 + \alpha} & (k = 1, \cdots, K) \\
\frac{\alpha}{N - 1 + \alpha} & (k = K + 1)
\end{array}
\right.

$k=1, \cdots, K$が既存のクラスタ、$k=K+1$が新しいクラスタに相当します。ここで、$n_k$は$i \neq j, z_j=k$となる$j$の個数、$\alpha$は正の値を取るハイパパラメータです。式から分かる通り、$\alpha$が大きほど新しいクラスタに割当される確率が大きくなります。今後、この確率分布を${\rm CRP}(k | z_{1:N}^{\backslash i}, \alpha)$と書くことにします。

次に、$x_{1:N}$を観測したときの確率分布を考えます。通常の混合ガウスモデルでは、次のようにして求めることができます。

\begin{align}
P(z_i = k | x_{1:N}, \mu_{1:K}, \Sigma_{1:K}, \pi) 
&= \frac{P(z_i = k | \pi) N(x_i | \mu_k, \Sigma_k)}{\sum_{l=1}^{K}P(z_i = l | \pi) N(x_i | \mu_l, \Sigma_l)} \\
&= \frac{\pi_k N(x_i | \mu_k, \Sigma_k)}{\sum_{l=1}^{K}\pi_l N(x_i | \mu_l, \Sigma_l)}
\end{align}

$N(x | \mu, \Sigma)$は多変量ガウス分布の確率密度関数です。
無限混合ガウスモデルでは、上式の$P(z_i = k | \pi)$を中華料理過程の式に変えることで確率分布を求めることができます。

\begin{align}
P(z_i = k | x_{1:N}, \mu_{1:K}, \Sigma_{1:K}, z_{1:N}^{\backslash i}, \alpha) 
&= \frac{{\rm CRP}(k | z_{1:N}^{\backslash i}, \alpha) N(x_i | \mu_k, \Sigma_k)}{\sum_{l=1}^{K+1}{\rm CRP}(l | z_{1:N}^{\backslash i}, \alpha) N(x_i | \mu_l, \Sigma_l)}
\end{align}

ここで、ひとつ注意があります。既存のクラスタ$k=1, \cdots, K$については、各クラスタに所属するデータから平均$\mu_k$と分散共分散行列$\Sigma_k$を推定することが可能です。
ところが、新しいクラスタ$k=K+1$については、クラスタに所属するデータがないため、直接的な平均と分散共分散行列の推定ができません。
そこで、今回は次の方法を採用します。

• $\mu_k$に事前分布を仮定し、同時分布$P(x_i, \mu_k| \cdots)$を周辺化して、$\mu_k$を陽に推定せずに確率分布を求める
• $\Sigma_k$を全てのクラスタで共通の行列$\sigma^2 I$とする。（$\sigma^2$はハイパパラメータ）

平均に関して、ベイズ推定と周辺化によって、事前分布のみから確率分布を計算します。
分散共分散行列に関して、上記のように固定してもクラスタリングには十分なので、今回は推定なしにします。⁴もちろん、分散を推定してクラスタリングすることも可能です。

場合分けするのはナンセンスなので、この手法を$k=1, \cdots, K+1$のときにも適用します。
なお、周辺化については前回の記事をご参照ください。

ノンパラメトリックベイズによるクラスタリング

$N$個のデータのクラスタ$z_{1 :N}$を推定するため、$z_1$から$z_N$まで一つずつ順次確率的にサンプリングすることにします。この方法をギブスサンプリングといい、今回のように、無限次元の確率分布をもとにした方法を特にノンパラメトリックベイズといいます。

この方法を実行するため、初期のクラスタリング$z_{1:N}$を何かしら仮定する必要があります。今回は全て同じクラスタ$k=1$を初期クラスタリングとして、$z_1$から$z_N$までのサンプリングを複数回反復することにします。最初は1つしかクラスタがありませんが、平均から遠いデータが新しいクラスタと結びついて、クラスタ数を増やしながらクラスリングが進むイメージです。

実装

無限混合ガウスモデルによるクラスタリングを実装しました。
前回の実装をベースにしていますが、全体的にブラッシュアップしており、差分が大きいかもしれません。

ソースコード

from collections import Counter
from functools import partial
from typing import Sequence, Tuple
import numpy as np
from scipy.special import logsumexp


def _log_gaussian(x: np.ndarray, mu: np.ndarray, var: float) -> np.ndarray:
    # 前回記事では2*np.piを省略していましたが、分かりやすさのため今回は計算に入れます
    d = x.shape[0]
    norm_squared = ((x - mu) * (x - mu)).sum(axis=1)
    return - d * np.log(2 * np.pi * var) / 2 - norm_squared / (2 * var)


# クラスタを格納する配列クラス
class InfiniteClusterArray(object):
    def __init__(self, array: Sequence[int]):
        self._array = np.array(array, dtype=int)
        self._counter = Counter(array)
        self._n_clusters = max(array) + 1

    def update(self, i: int, k: int):
        assert k <= self._n_clusters

        pre_value = self._array[i]
        if pre_value == k:
            return

        if self._counter[pre_value] > 0:
            self._counter[pre_value] -= 1
        self._array[i] = k
        self._counter[k] += 1

        if k == self._n_clusters:
            self._n_clusters += 1

        if not self._counter[pre_value]:
            self._array -= 1 * (self._array > pre_value)
            self._counter = Counter(self._array)
            self._n_clusters -= 1

    @property
    def array(self) -> np.ndarray:
        return self._array.copy()

    @property
    def n_clusters(self) -> int:
        return self._n_clusters

    def count(self, k: int, excluded=None) -> int:
        assert (excluded is None) or isinstance(excluded, int)
        if excluded is None:
            return self._counter[k]

        counted = self._counter[k] - bool(self._array[excluded] == k)
        return counted

    def __getitem__(self, i: int) -> int:
        return self._array[i]



# データセットを格納する配列クラス
class DataArray(object):
    def __init__(self, array: np.ndarray):
        assert array.ndim == 2
        self._array = array

    @property
    def n(self) -> int:
        return self._array.shape[0]

    @property
    def d(self) -> int:
        return self._array.shape[1]

    def mean(self, excluded=None) -> np.ndarray:
        assert (excluded is None) or isinstance(excluded, Sequence)

        if excluded is None:
            return self._array.mean(axis=0)
        else:
            return self._array[excluded].mean(axis=0)

    def __getitem__(self, i):
        return self._array.__getitem__(i)

    def __iter__(self):
        for i in range(self.n):
            yield self._array[i]

    def __str__(self) -> str:
        return f'DataArray({self._array})'


# 無限混合ガウスのギブスサンプリングクラス
class IGMMGibbsSampling(object):
    def __init__(self, var: float, var_pri: float, alpha: float, seed=None):
        # 事前分布の平均は、データセット全体の平均に設定されます。
        assert (var > 0) and (var_pri > 0)
        assert alpha > 0

        self.var = var
        self.var_pri = var_pri
        self.alpha = alpha
        self._random = np.random.RandomState(seed)

    def predict_clusters(self, X: np.ndarray, n_iter: int, init_z=None) -> np.ndarray:
        # クラスタを推定するメソッド
        assert X.ndim == 2
        X = DataArray(X)
        mu_pri = X.mean()

        if init_z is None:
            init_z = np.zeros(X.n, dtype=int)
        z = InfiniteClusterArray(init_z)

        for _ in range(n_iter):
            compute_probs = partial(self._compute_pmf_zi, X, z, mu_pri)
            for i in range(X.n):
                probs = compute_probs(i)
                z_i = self._random.multinomial(1, probs)
                z_i = np.where(z_i)[0][0]  # One-hot to index
                z.update(i, z_i)

        return z.array

    def _compute_pmf_zi(self, X: DataArray, z: InfiniteClusterArray, mu_pri: np.ndarray, i: int) -> np.ndarray:
        mu, var = self._compute_marginal_distribution(X, z, mu_pri, i)
        pi = np.array([
            z.count(k, excluded=i) / (X.n - 1 + self.alpha)
            if k < z.n_clusters
            else self.alpha / (X.n - 1 + self.alpha)
            for k in range(z.n_clusters + 1)
        ])

        log_pdf_xi = _log_gaussian(X[i], mu, var)
        with np.errstate(divide='ignore'):
            pmf_zi = np.exp(np.log(pi) + log_pdf_xi - logsumexp(np.log(pi) + log_pdf_xi)[np.newaxis])

        return pmf_zi

    def _compute_marginal_distribution(self, X: DataArray, z: InfiniteClusterArray,
                                       mu_pri: np.ndarray, i: int) -> Tuple[np.ndarray, float]:

        n_vec = np.array([
            z.count(k, excluded=i)
            for k in range(z.n_clusters + 1)
        ], dtype=np.int)

        tmp_vec = 1 / (n_vec / self.var + 1 / self.var_pri)

        mu = np.tile(mu_pri / self.var_pri, (z.n_clusters + 1, 1))  # Shape (z.n_clusters + 1, X.d)
        for k in range(z.n_clusters):
            excluded_list = [
                j for j in range(X.n)
                if (j != i) and (z[j] == k)
            ]
            if excluded_list:
                mean = X.mean(excluded=excluded_list)
                mu[k] += n_vec[k] * mean / self.var
        mu *= tmp_vec[:, np.newaxis]
        var = tmp_vec + self.var

        return mu, var

クラスタリングの実践

試しに、実装した自作クラスを使ってクラスタリングしてみます。
今回は、次のデータセットに適用します。

ソースコード

Jupyter Notebookを使用しています。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
%matplotlib inline

np.random.seed(0)
cmap = plt.rcParams['axes.prop_cycle'].by_key()['color']

pi = [0.6, 0.2, 0.2]
mu = np.array([[1, 2], [-2, 0], [3, -2]])
var = [0.8, 0.4, 0.2]
size = 100
ndim = 2

def make_data():
    z = np.random.multinomial(1, pvals=pi)
    z = np.where(z==1)[0][0]
    cov = var[z] * np.eye(2)

    return np.random.multivariate_normal(mu[z], cov)

X = np.array([make_data() for _ in range(100)])

fig = plt.figure()
ax = plt.axes()
ax.scatter(X[:, 0], X[:, 1], color='gray')
ax.set_aspect('equal')
plt.show()

序盤でお見せしたデータセットの、まだクラスタリングされていない状態のものです。
データの個数は100個で、3つのガウス分布からランダムに選んで作成しています。
したがって、3つのクラスタに分けるのが適切です。

それでは、クラスタリングしてみましょう。

model = IGMMGibbsSampling(var=0.5, var_pri=4, alpha=0.1, seed=0)  # ハイパパラメータをコンストラクタにわたす
clusters = model.predict_clusters(X, n_iter=10)  #クラスタリングの実行

# 可視化
fig = plt.figure()
ax = plt.axes()
for k in range(max(clusters) + 1):
    data = np.array([
        X[i, :] for i in range(X.shape[0])
        if clusters[i] == k
    ])
    ax.scatter(data[:, 0], data[:, 1], label=f'cluster {k + 1}', color=cmap[k])

ax.set_aspect('equal')
plt.legend()
plt.show()

クラスタで共通の分散varはデータセット小さめの$0.5$、平均の事前分布の分散var_priは、一様分布っぽくするために大きめの$4$、中華料理過程のパラメータalphaは$0.1$にしています。
ちなみに、平均の事前分布の平均は、データセット全体の平均になるように実装されています。
クラスタリングメソッドpredict_clustersに、データセットと反復回数n_iterを渡しています。ここで、反復回数を10としました。

結果は、次のようになりました。

クラスタ数を指定せず、全て同じクラスタの状態からスタートしたのにもかかわらず、3個のクラスタに分かれてくれました。

せっかくなので、反復回数ごとのクラスタリングの推移を調べてみます。
n_iter=2, 4, 6, 8, 10におけるクラスタリング結果を図示しました。

• n_iter = 2

• n_iter = 4

• n_iter = 6

• n_iter = 8

• n_iter = 10

初期状態から多めにクラスタができ、そこから徐々にクラスタが減って3個になったようです。
ちなみに、ここからn_iterを増やしても、クラスタ数は基本的に3個の状態を維持していました。
※たまにデータが1個だけのクラスタができて、クラスタ数が4個になることもありました。

ハイパパラメータの調整

中華料理過程のハイパパラメータ$\alpha$は、上述した通り、大きほど新しいクラスタが出来やすくなります。
もともとクラスタ数の推定を省略するという動機があり、この$\alpha$の調整に苦戦してしまったら導入した意義が薄れてしまいます。
$\alpha$を変えながら試行したところ、基本的には小さめにすれば問題ないという印象を受けました。
次の図は、$\alpha=0.0001$、反復回数10回の結果です。

しっかり3つのクラスタに分かれています。
$\alpha$をかなり小さめにしても、反復回数に大きなインパクトなく、適切にクラスタリングできるようです。

ちなみに、$\alpha=1$にすると、10回反復時点で5個のクラスタができ、その後も3個に減ることはありませんでした。$\alpha$が大きすぎるとうまくクラスタ数が決まらないようです。
機械学習で使われる最適化アルゴリズムのステップ幅と同じ性質を持っていますね。

おわりに

以上、無限混合ガウスモデルによるクラスタリングを試してみました。
狙い通りのクラスタ数に分かれてくれたのでおもしろかったです。

---

1. $N$個の変数$x_1, \cdots, x_N$を、省略して$x_{1:N}$と書くことにします。また、その中から$i$番目を除いた$N-1$個の変数を$x_{1:N}^{\backslash i}$と書くことにします。 ↩

2. $\pi_1, \cdots, \pi_K \ge 0$及び$\sum_{k=1}^K\pi_k=1$を満たす$(\pi_1, \cdots, \pi_K)$の集合を$(K-1)$次元単体と言います。もともとは幾何学で使われる集合で、$0$次元単体は1点、$1$次元単体は線分、$2$次元単体は三角形を表しています。ベイズ推定の本を読んでると、確率の文脈でこの集合が使われるのを目にします。 ↩

3. $\pi$を$z_{1:N}^{\backslash i}$の観測から最尤推定する方法もありますが、この方法ではクラスタ$k$が消滅するすると常に$\pi_k=0$となり、クラスタが復活しなくなるため注意が必要です。 ↩

4. $\Sigma_k$をクラスタ共通の固定値にしても、k-meansと同様の精度でクラスタリングができると思います。なぜなら、k-means法は、上記のように分散を固定した混合ガウスモデルにおいて、各$\mu_k$や$z_i$を最尤推定した場合と同等だからです。 ↩

問題

https://atcoder.jp/contests/abc134/tasks/abc134_d

難しかった。というか問題文を読むことが難しかった。

今回の学び

1. 任意の整数＝全ての整数のようなニュアンス
2. 後ろから順に解いていくと一意に確定する
3. $O(\sum(N/i)) = O(NlogN)$

計算量の求め方

発生した疑問：$O(\sum(N/i))$って間に合うの？

今回の問題では任意のiにおいて$N/i$を求める必要がある。これを考えるとき、積分を行うと良い。$\int(1/N) = logN$であることから最悪でもlogNで押さえつけることができる。

まとめ

1. 数学的な文章になれる
2. 配列を前もしくは後ろから判定するのはよくあるし気づかなければいけない
3. 積分を考えると面積から計算時間を求めることができる。

コード

N = int(input())
a = list(map(int, input().split()))
a.insert(0, 0)
b = [0] * N
M = 0
for i in range(N, 0, -1):
    ball = 0
    for j in range(i + i, N + 1, i):
        ball ^= b[j - 1]
    b[i - 1] = ball ^ a[i]
M = sum(b)
print(M)
for i in range(N):
    if b[i]:
        print(i + 1, end = ' ')

嘘つき族と正直族の問題をPythonで解いてみた。
答え合わせはしていない。Python楽。

import itertools

# 嘘つきなら反対を答える
def ans( person, liars, b):
    if person in liars:
        b = not b
    return b

def print_res(liars, cards):
    print('-----')
    for i in range(4):
        print( '{}:{} {}'.format('ABCD'[i], '嘘' if i in liars else '正', cards[i]))

# 嘘つきな人 と カードの組み合わせ
liars_g = itertools.combinations(list(range(4)),2) # 嘘つきな人
cards_g = itertools.permutations(list(range(1,5))) # カード
for liars, cards in itertools.product(liars_g,cards_g):
    # Ａさん：私のカードは，偶数です．
    # Ｂさん：私のカードは，３か４のどちらかです．
    # Ｃさん：Ｂさんは，正直族です．
    # Ｄさん：私のカードは，１です．
    if  ans( 0, liars, cards[0] % 2 == 0) and \
        ans( 1, liars, cards[1] in (3,4)) and \
        ans( 2, liars, 1 not in liars) and \
        ans( 3, liars, cards[3] == 1):
        print_res(liars, cards)

答え

-----
A:嘘 1
B:正 3
C:正 2
D:嘘 4
-----
A:嘘 1
B:正 3
C:正 4
D:嘘 2
-----
A:嘘 1
B:正 4
C:正 2
D:嘘 3
-----
A:嘘 1
B:正 4
C:正 3
D:嘘 2
-----
A:嘘 3
B:正 4
C:正 1
D:嘘 2
-----
A:正 4
B:嘘 2
C:嘘 3
D:正 1

当記事の記載範囲

この記事ではPythonでmongodbに接続してから、aggregate（SQLで言うところの集計関数）の使い方について記載します。
mongodbの起動やpymongoのインストール方法については以下の記事をご覧いただければ幸いです。
https://qiita.com/bc_yuuuuuki/items/2b92598434f6cc320112

準備データ

準備データは以下の記事でmongoDBに突っ込んだQiitaの記事情報を使用します
[Python]Qiitaの記事情報をmongoDBに突っ込んだ

aggregateの使い方

mongoDBのaggregateの使い方はSQLに慣れているとイマイチピンと来ないです。
下表がSQLとaggregateの対比表です。

SQL	aggregate
WHERE	$match
GROUP BY	$group
HAVING	$match
SELECT	$project
ORDER BY	$sort
LIMIT	$limit
SUM()	$sum
COUNT()	$sum

mongoDBの操作クラス

pymongoを使用して色々とmongoDBを使用するクラスを作ってます。

mongo_sample.py

from pymongo import MongoClient

class MongoSample(object):

    def __init__(self, dbName, collectionName):
        self.client = MongoClient()
        self.db = self.client[dbName] #DB名を設定
        self.collection = self.db.get_collection(collectionName)

    def aggregate(self, filter, **keyword):
        return self.collection.aggregate(filter, keyword)

aggregate呼び出し用の関数を作っているだけです。

mongoDBからデータを取得する

まずは、コードです。

aggregate_sample.py

from mongo_sample import MongoSample
import pprint
# arg1:DB Name
# arg2:Collection Name
mongo = MongoSample("db", "qiita")

# 最大値
pipeline = [
    {"$group":{ "_id":"title","page_max_view":{"$max":"$page_views_count"}}}
]
results = mongo.aggregate(pipeline)
print("------------------------最大値-----------------------------")
pprint.pprint(list(results))

# 最小値
pipeline = [
    {"$group":{ "_id":"title","page_min_view":{"$min":"$page_views_count"}}}
]
results = mongo.aggregate(pipeline)
print("------------------------最小値-----------------------------")
pprint.pprint(list(results))

# 平均値
pipeline = [
    {"$group":{ "_id":"average","page_average_view":{"$avg":"$page_views_count"}}}
]

# 合計
pipeline = [
    {"$group":{"_id":"page_total_count","total":{"$sum":"$page_views_count"}}}
]
results = mongo.aggregate(pipeline)
print("------------------------平均値-----------------------------")
pprint.pprint(list(results))

# tag毎の出現回数カウント
pipeline = [
    { "$unwind": "$tag_list"}, 
    { "$group": { "_id": "$tag_list", "count": { "$sum":1}}},
    { "$sort": {"count": -1, "_id":1}}
]

results = mongo.aggregate(pipeline)
print("------------------------集計値-----------------------------")
pprint.pprint(list(results))

やっていることは大したことはありません。
最大値、最小値、平均値、タグ毎のカウントを取得しています。

pprintはinstallする必要があります。

pip install pprint

それぞれ、mongoDBの操作方法と見比べていきます。

最大値／最小値／平均値／合計値

まずは、mongoDBのコマンド
例は最大値のみとします。maxをminやavg、sumに変えれば最小値／平均値／最大値になります。

db.qiita.aggregate([{$group:{_id:"page_max_views",total:{$max:"$page_views_count"}}}])

aggregate_sample.py

pipeline = [
    {"$group":{ "_id":"title","page_max_view":{"$max":"$page_views_count"}}}
]

実行結果

[{'_id': 'title', 'page_max_view': 2461}]

このやり方だと、"_id"を"title"に固定して、全レコード中の最大値を取得するという結果になりました。

ただ、どの記事が一番読まれているのかを知りたいので記事のタイトルを表示したいです。

mongoDBコマンド

> db.qiita.aggregate([{$project:{title:1,page_views_count:1}},{$group:{_id:"$title", total:{$max:"$page_views_count"}}},{$sort:{total:-1}}])
{ "_id" : "Pythonでmongodbを操作する～その２：find編～", "total" : 2461 }
{ "_id" : "Pythonでmongodbを操作する～その３：update編～", "total" : 1137 }
{ "_id" : "Pythonでmongodbを操作する～その４：insert編～", "total" : 1102 }
{ "_id" : "pymongoを使った様々な検索条件（AND／OR／部分一致／範囲検索）", "total" : 1019 }
（略）

このコマンドだと、記事のタイトルとページの閲覧回数を見れました。
当たり前ですが、記事名でグループ化されているのであまり意味のない集計ですね。。
グループ化が必要ない最大値であれば、findでsortしてlimitをかけるのが良さそうです。

tag1毎の最大値を取得してみます。

> db.qiita.aggregate([{$group:{_id:"$tag1", total:{$max:"$page_views_count"}}},{$sort:{total:-1}}])
{ "_id" : "Python", "total" : 2461 }
{ "_id" : "Vagrant", "total" : 946 }
{ "_id" : "Java", "total" : 617 }
{ "_id" : "Hyperledger", "total" : 598 }
{ "_id" : "solidity", "total" : 363 }
{ "_id" : "Ethereum", "total" : 347 }
{ "_id" : "ブロックチェーン", "total" : 232 }
{ "_id" : "Blockchain", "total" : 201 }
{ "_id" : "coverage", "total" : 199 }

はい。良い感じで取得出来ました。

一応、pythonのコードも変えてみます。

aggregate_sample.py

# 最大値
pipeline = [
    {"$group":{ "_id":"$tag1","page_max_view":{"$max":"$page_views_count"}}}
]
results = mongo.aggregate(pipeline)
print("------------------------最大値-----------------------------")
pprint.pprint(list(results))

タグ毎の集計

タグ毎に何個記事を書いているか集計したいと思います。
集計はtag_listという項目を使用しますが、このデータは以下のようになっています。

> db.qiita.find({},{_id:0,tag_list:1})
{ "tag_list" : [ "Python", "MongoDB", "Python3", "pymongo" ] }
{ "tag_list" : [ "Python", "Python3" ] }
{ "tag_list" : [ "Python", "Python3", "Blockchain", "ブロックチェーン", "Hyperledger-Iroha" ] }
{ "tag_list" : [ "Blockchain", "ブロックチェーン", "Hyperledger-Iroha" ] }
{ "tag_list" : [ "Blockchain", "Ethereum", "Hyperledger", "Hyperledger-sawtooth" ] }
{ "tag_list" : [ "ブロックチェーン", "Hyperledger", "Hyperledger-sawtooth" ] }
{ "tag_list" : [ "Java", "ブロックチェーン", "Hyperledger", "Hyperledger-Iroha" ] }
{ "tag_list" : [ "ブロックチェーン", "Hyperledger", "Hyperledger-Iroha" ] }
{ "tag_list" : [ "Java", "Ethereum", "ブロックチェーン", "Hyperledger", "Hyperledger-Iroha" ] }
{ "tag_list" : [ "Java", "ブロックチェーン", "Hyperledger", "Hyperledger-Iroha" ] }
{ "tag_list" : [ "Hyperledger", "Hyperledger-Iroha", "Hyperledger-burrow", "Hyperledger-sawtooth", "Hyperledger-besu" ] }
{ "tag_list" : [ "Vagrant", "VirtualBox", "Hyper-V" ] }
{ "tag_list" : [ "Java", "Ethereum", "solidity", "ブロックチェーン", "web3j" ] }
{ "tag_list" : [ "Java", "Ethereum", "ブロックチェーン", "web3j" ] }
{ "tag_list" : [ "Java", "Ethereum", "ブロックチェーン", "web3j" ] }
{ "tag_list" : [ "Java", "Ethereum", "ブロックチェーン", "web3j" ] }
{ "tag_list" : [ "Java", "Ethereum", "solidity", "ブロックチェーン", "web3j" ] }
{ "tag_list" : [ "Java", "Ethereum", "ブロックチェーン", "web3j" ] }
{ "tag_list" : [ "Java", "Ethereum", "ブロックチェーン", "web3j" ] }
{ "tag_list" : [ "Ethereum", "ブロックチェーン" ] }

SQLでこんな形式で格納されているデータを集計するとかなり面倒くさいですね。。

mongoDBではunwindというものを使うことでLIST形式のデータを分割して集計可能になります。

> db.qiita.aggregate( { $project:{tag_list:1}}, { $unwind: "$tag_list"}, { $group: { _id: "$tag_list", count: { $sum:1}}},{ $sort: {"count": -1, "_id":1}} )
{ "_id" : "ブロックチェーン", "count" : 16 }
{ "_id" : "Ethereum", "count" : 11 }
{ "_id" : "Java", "count" : 10 }
{ "_id" : "Python", "count" : 9 }
{ "_id" : "Python3", "count" : 9 }
{ "_id" : "Hyperledger", "count" : 7 }
{ "_id" : "Hyperledger-Iroha", "count" : 7 }
{ "_id" : "MongoDB", "count" : 7 }
{ "_id" : "web3j", "count" : 7 }
{ "_id" : "solidity", "count" : 4 }
{ "_id" : "Blockchain", "count" : 3 }
{ "_id" : "Hyperledger-sawtooth", "count" : 3 }
{ "_id" : "Hyper-V", "count" : 1 }
{ "_id" : "Hyperledger-besu", "count" : 1 }
{ "_id" : "Hyperledger-burrow", "count" : 1 }
{ "_id" : "Vagrant", "count" : 1 }
{ "_id" : "VirtualBox", "count" : 1 }
{ "_id" : "coverage", "count" : 1 }
{ "_id" : "pymongo", "count" : 1 }
{ "_id" : "truffle", "count" : 1 }

pythonのコードは「{ "$project": {"tag_list": 1}}」を入れていません。
あっても無くても、結果が変わりませんでした。
イマイチ、このprojectの使い方がよくわからない。

感想

SQLに慣れていると分かり辛い部分が多々あるが、unwindなどを使用することで柔軟な集計が出来そう

関連記事

• Pythonでmongodbを操作する～その１～
• Pythonでmongodbを操作する～その２：find編～
• Pythonでmongodbを操作する～その３：update編～
• Pythonでmongodbを操作する～その４：insert編～
• Pythonでmongodbを操作する～その５：delete編～
• pymongoを使った様々な検索条件（AND／OR／部分一致／範囲検索）

はじめに

SparkのDataFrameで作ったデータをPythonの各モジュールで使いたい時、toPandas()メソッドを利用してPandasのDataFrameに変換するが、その際にメモリエラーが起きる場合がしばしばある。
メモリに載せられるように試行錯誤した中で、有効そうなものをまとめた。

もっといい手段がありそうなので、ご存じの方は是非教えてください！

やり方

daskを利用して変換する

sparkによる変換は、spark.driver.memoryやspark.driver.maxResultSizeの影響を受けるが、daskの場合はそれがないため、エラーを回避しやすい。

daskを利用した変換

import dask.dataframe as dd

df.write.parquet(parquet_path)
dask_df = dd.read_parquet(parquet_path)
pandas_df = dask_df.compute()

データ型を変更する

バイト数を削減するように、変数のデータ型を変更する。

※勿論、バイト数を削減することで、その後の演算の精度は損なわれる

データ型を変更

# 例えば、int32型(4バイト)をint8型(1バイト)へ変換する
dask_df = dask_dt.astype({k: 'int8' for k in dask_df.dtypes[dask_df.dtypes == 'int32'].index})

最近，新型肺炎が流行っていますね．
以前にSIRモデルという感染症の数理モデルを使った研究をしていたので，新型肺炎のモデルに当てはめてみました．数理モデルを使うと，感染症の未来を予測できます．

今回は，発生の中心である 武漢 に注目して，感染症のモデル化と感染症の未来を予測します．

keyword

疫学, SIR model, nCoV-2019, 新型肺炎, 新型コロナウィルス

SIRモデル

SIRモデルは，感染した人数の推移を微分方程式として表したモデルです（Wikipediaでもわかりやすく解析されています）． SIRモデルでは，感染症に対して人は3つの状態を取ると考えられます．

1. 感染する可能性がある人 : S
2. 感染している人：I
3. 感染から治って免疫を獲得した人，または死亡した人 : R

時刻tにおける感染する可能性のある人，感染している人と感染から治った人をそれぞれS(t), I(t)とR(t)表記すると，SIRモデルは，

$$
\dot{S}(t) = -\beta S(t)I(t),\
\dot{I}(t) = \beta S(t)I(t) - \gamma I(t)\
\dot{R}(t) = \gamma I(t)
$$

と記述されます．
ここで，βは感染率を表し，γは回復率(+死亡率)を表しています．
感染者の増え方は，感染率β，感染する可能性がある人S(t)と感染している人I(t)に比例します．

死亡した人も感染を引き起こさないので，感染から治った人と同一視することに注意してください．

ここで，
$$
S(t) + I(t) + R(t) = N
$$
は一定となり，その地域の人口と一致します．
今回は，武漢の人口を使います．

武漢における新型肝炎の感染データをつかって，感染率βと回復率γを学習し，武漢の未来を予測します．

使ったデータ

感染データはkaggleに公開されているこちらから引用しています．
また，武漢の人口データは，こちらのに記述されている2017年の人口統計データを用いています．

SIRモデルを使ったパラメーターの学習

武漢における，2020年1月22日から2月4日までの感染の推移は，以下のようになります．

青線が感染者数，オレンジ線が死者数，緑線が回復した人数です．
回復した人数と死者数がだいたい同じであることは，少しおかしいかな？と思います．

このデータを表現できるように，SIRモデルでフィッティングしてみました．

Recovered peopleは回復した人と死者数の和としています．
また，青点が実際に観測されている感染者数，青線がSIRモデルで近似した結果です．
オレンジ点とオレンジ線は回復した人数の実測値と予測値です．

十分に近似できていそうです．

武漢の感染の未来を予測

十分に近似できていそうなので，学習したパラメータを使って，武漢の感染の未来を予測してみました．

次の図が，2月4日から10日間の予測になります．

点が実測値，線が予測値です．
SIRモデルよれば，まだまだ増加するようです．

次が2月4日から一年間の予測になります．

感染が全然止まらいようです．

考察

• SIRモデルを用いて武漢における新型肺炎の感染の拡大を予測してみましたが，拡大が止まらないという予測結果になりました．
• 武漢において回復した人数が正確に計量されていないことが原因だと思われます．
• 他のデータを用いて，回復率γを正確に測定できれば，より良い予測ができると考えられます．

今後の課題

次は，交通量をもとに，中国全土の感染の広がりを予測しようと思います．

コード

https://github.com/yuji0001/202002nCoV_analysis

はじめに：

CNNモデルにおいてよくPooling処理は畳み込み処理(Conv2D)と区別して解説が行われます。
しかし、勘の良い人なら(2,2)サイズのAveragePooling2D()は畳み込みフィルタ((0.25,0.25),(0.25,0.25))をstrides=2で畳み込み処理させたものに等しいことに既に気付いているでしょう。
一方でMaxPooling2DをConv2Dの畳み込み処理だけで表現することが可能かと問いかければ一見不可能に思えます。しかし、よくよく考えてみると実はConv2DのみでMaxpooling2Dを表現可能でしたので、それについて書いてみます。

max関数とreluの関係：

$max(a,b)$を$relu(x)$を使って表現することを考えます。

max(a,b) = relu(a-b) + b

上記の時、$a>b$で$max(a,b)=a$で、$a\leqq b$で$max(a,b)=b$であることが分かります。
つまり、$max$関数は$relu$関数を利用して表現可能ということです。

さて、MaxPooling2Dを考えるにあたり、pool_size(2,2)として4つの要素から最大値を取り出す作業を考えます。$max(a,b,c,d)$を考えた場合、$relu$関数を利用して下記のように記述可能です。

\begin{align}
max(a,b,c,d) &= max(max(a,b), max(c,d))\\
&=max(relu(a-b)+b,relu(c-d)+d)\\
&=relu((relu(a-b)+b)-(relu(c-d)+d))+relu(c-d)+d\\
&=relu(relu(a-b) - relu(c-d) + (b - d)) + relu(c-d) + d
\end{align}

さて、ここで$relu(a-b)$を畳み込みフィルタが((1,-1),(0,0))をstrides=2掛けた結果に活性化関数$relu$を掛けたもの見なすとしましょう。
同様に$relu(c-d)$は畳み込みフィルタが((0,0),(1,-1))に活性化関数$relu$を、
$b$は畳み込みフィルタが((0,1),(0,0))に活性化関数に恒等関数を、
$d$は畳み込みフィルタが((0,0),(0,1))に活性化関数に恒等関数を掛けたものとみなします。
これらの足し合わせ結果に$relu$と恒等関数を掛けた結果を足し合わせればMaxpoolingを再現できるのです。

実証コード：

下記のように式から求まった畳み込みフィルタ重みを定義してやり、その重みを.set_weightsで設定しました。

import numpy as np
from keras.layers import Input, Conv2D, Add, Concatenate
from keras.models import Model

inputs = Input(shape=(16,16,1))
x1 = Conv2D(1, (2, 2), strides=2, padding='same', activation='relu', use_bias=False)(inputs)
x2 = Conv2D(1, (2, 2), strides=2, padding='same', activation='relu', use_bias=False)(inputs)
x3 = Conv2D(1, (2, 2), strides=2, padding='same', activation='linear', use_bias=False)(inputs)
x4 = Conv2D(1, (2, 2), strides=2, padding='same', activation='linear', use_bias=False)(inputs)
x5 = Concatenate()([x1, x2, x3, x4])
x6 = Conv2D(1, (1, 1), activation='relu', use_bias=False)(x5)
x7 = Conv2D(1, (1, 1), activation='linear', use_bias=False)(x5)
outputs = Add()([x6, x7])

model = Model(inputs=inputs, outputs=outputs)
model.summary()

weight1 = np.array([[[[1]],[[-1]]],[[[0]],[[0]]]]) # relu(a-b)
weight2 = np.array([[[[0]],[[0]]],[[[1]],[[-1]]]]) # relu(c-d)
weight3 = np.array([[[[0]],[[1]]],[[[0]],[[0]]]])  # b
weight4 = np.array([[[[0]],[[0]]],[[[0]],[[1]]]])  # d
weight5 = np.array([[[[1],[-1],[1],[-1]]]])        # relu(a-b) - relu(c-d) + b - d
weight6 = np.array([[[[0],[1],[0],[1]]]])          # relu(c-d) + d

model.get_layer(name='conv2d_1').set_weights([weight1])
model.get_layer(name='conv2d_2').set_weights([weight2])
model.get_layer(name='conv2d_3').set_weights([weight3])
model.get_layer(name='conv2d_4').set_weights([weight4])
model.get_layer(name='conv2d_5').set_weights([weight5])
model.get_layer(name='conv2d_6').set_weights([weight6])

X = np.random.randint(-10,11,(1,16,16,1))
Y = model.predict(X)
print('X=\n',X[0,:,:,0])
print('Y=\n',Y[0,:,:,0])

適当な入力を与えた時、その出力が(2,2)サイズのMaxpooling2D()と等価であることを確認しました。
モデルにはConv2D()とConcatenate()とAdd()しか使っていません。

X=
 [[ -7   7   0  -8   8  -3  -1   7  -6   9   4 -10   8   7  -6  10]
 [ -4  -5  -5   0 -10   7   1   8   1  -9  10  -3   5 -10   5  -9]
 [ -7   9   6  -9   0  -7   3   0   4   9  -6  -1   9   1   0   0]
 [  1  -3  -7  -5   7   3   6   7  -4  -2   6  -8   7  -6   0  -2]
 [ -2  -6   9   4   4   3  10   3   9   9  -5   2   0   2   9  -3]
 [  2   7   5  -3   9  -7  -1 -10   7  -5  -4  -6   0   7   8 -10]
 [  1  -3  -3   9  -5  -6  -7  -7  -4   9  -7  -9  -6   2   1  -9]
 [ -1  -5  -3   1  -2   9   0  10 -10   5  -9  -8  -2   8  -4   3]
 [  1  -4  -2  -5  -2   3   5   4  -5   3  -6   9   0   2  -3   6]
 [  6   1   4  -8  -6   7  -8   4 -10 -10  -5   7  -8  -7  -1   5]
 [  8  -2   4   9   6   9 -10  -4  -3  -9   7   1  -7   4   7   0]
 [ -6   5   6  -1  -8  -2   0   0   6   3  10  -3   3   9   1  -2]
 [  2   3  -6   6  -1   1   9  -2  -3   2   4   5 -10  -7   5   4]
 [ -5   5   0   9   4   2 -10  -8   7   4  -7   2  -8   7  -3   3]
 [  5   0   3   2  -4   2  -3  10   1  -7  -7   2   7   5  -4   2]
 [  0  -9   6   2   1  -2  -4   3  -4   7   9  -9   7  -5   4  -1]]
Y=
 [[ 7.  0.  8.  8.  9. 10.  8. 10.]
 [ 9.  6.  7.  7.  9.  6.  9.  0.]
 [ 7.  9.  9. 10.  9.  2.  7.  9.]
 [ 1.  9.  9. 10.  9. -7.  8.  3.]
 [ 6.  4.  7.  5.  3.  9.  2.  6.]
 [ 8.  9.  9.  0.  6. 10.  9.  7.]
 [ 5.  9.  4.  9.  7.  5.  7.  5.]
 [ 5.  6.  2. 10.  7.  9.  7.  4.]]

まとめ：

Maxpooling2D()とてConv2d()を用いて等価なモデルを記述可能である。
a,b,c,dの最大値の取り方は$max(max(a,b), max(c,d))$の代わりに$max(max(a,c), max(b,d))$でもいいので他の係数でもMaxpooling2D()を再現できるかもしれない。
とはいえ、活性化関数に恒等関数を選ぶ必要があるので(あまり見掛けない)、普通のモデルでConv2d()のみで偶然にMaxpooling2D()と等価なレイヤーを作ってしまうことはあまりないかと思われます。
(identityの結合があればあり得るのでしょうか？)