Prestoで1ヶ月後の時刻を求める時に気をつけること

普段の業務で利用しているPresto (treasuredata) で1ヶ月後の日付を求める機会があり、
こちらのドキュイメントになる、date_add という関数の挙動をテストした時に見つけた挙動のメモです。

6.13. Date and Time Functions and Operators

こちら、日付単位のデータのn日後やn日前を求める時には何も問題ないのですが、
nヶ月後を求める時に少し不思議な動きがありました。

まず前提として、元のデータの日付部分が、1〜28日の場合は、何も問題がありません。
5月16日の3ヶ月後は 8月16日ですし、
2ヶ月前は 3月16日です。


SELECT
    DATE_ADD('month', 3, timestamp '2019-05-16 12:00:00'),  -- 2019-08-16 12:00:00.000
    DATE_ADD('month', -2, timestamp '2019-05-16 12:00:00') -- 2019-03-16 12:00:00.000

問題は、月によって日数が違う点です。
例えば、8月は31日までありますが、9月は30日までしかありません。
そのため、8月31日の1ヶ月後をPrestoのDATE_ADDで算出すると、9月30日になります。
ちなみに、9月30日の1ヶ月前は8月30日です。
要するに、8月31日に1ヶ月足して1ヶ月引くと8月30日になり、元に戻らない。


SELECT
    DATE_ADD('month', 1, timestamp '2019-08-30'), -- 2019-09-30 00:00:00.000
    DATE_ADD('month', 1, timestamp '2019-08-31'),  -- 2019-09-30 00:00:00.000
    DATE_ADD('month', -1, timestamp '2019-09-30'),  -- 2019-08-30 00:00:00.000
    DATE_ADD('month', -1, timestamp '2019-10-01'),  -- 2019-09-01 00:00:00.000
    DATE_ADD('month', -1, DATE_ADD('month', 1, timestamp '2019-08-31')) -- 2019-08-30 00:00:00.000

こういう挙動を嫌って、いつもDATE_ADDは利用せず、60*60*24*30秒足したり引いたりするクエリを書いていたのですが、
30日と1ヶ月は厳密には違うので、その時に応じてよく考えて方法を選ぶ必要があります。

そして、問題はもう一点あります。
それは時刻まで含めて計算した時に、時刻の前後関係が前後することです。
8月30日と8月31日の1ヶ月後はどちらも9月30日と算出されますが、
これが時刻も含データの場合、時刻部分は1ヶ月足す前の値から変化しません。

要するに
(1)8月30日20時の1ヶ月後は9月30日20時で、
(2)8月31日7時の1ヶ月後は9月30日7時です。
元の時間は当然、(1)の方が前なのに、(1)と(2)の1ヶ月後の時刻は(2)の1ヶ月後の方が前になります。


SELECT
    DATE_ADD('month', 1, timestamp '2019-08-30 20:00:00'), -- 2019-09-30 20:00:00.000
    DATE_ADD('month', 1, timestamp '2019-08-31 07:00:00'),  -- 2019-09-30 07:00:00.000
    timestamp '2019-08-30 20:00:00' < timestamp '2019-08-31 07:00:00', -- true
    DATE_ADD('month', 1, timestamp '2019-08-30 20:00:00') < DATE_ADD('month', 1, timestamp '2019-08-31 07:00:00') --false

日付単位で行った時の不等号が等号になるのはまだ許容範囲かもしれませんが、
不等号の反転はちょっと困る。

幸い、nヶ月間以内の判定を時刻まで考慮して厳密に行う場面は少ない(これまでほぼなかった)ので
問題になることは少ないのですが、注意が必要です。
(timestampに変換する前に時刻を切り捨てるなどの処理を入れた方がいい。)

以上をまとめると、Prestoの DATE_ADDで月単位(month)の演算をする時の注意は次の3つです。

  1. 異なる日付の±nヶ月後が同じ日付になることがある
  2. ある日付のnヶ月後のnヶ月前が元の日付と異なることがある
  3. 2つの時間のnヶ月後を計算すると時間の前後関係が入れ替わることがある

とくに2月がからむと最悪で、
1月28,29,30,31日の1ヶ月後は全部2月28日になります。

requestsにタイムアウトを設定する

requestsを使ってurlのリストをチェックしていた時にredirectに加えて困ったのが、応答に時間がかかるサーバーへアクセスした時です。
数秒待って結果が戻って来ればまだいいのですが、そのままスクリプトが進まなくなってしまうことがありました。
これを防ぐには、timeoutを設定すると良いようです。
(デフォルトでは設定されていない)

ドキュメントはこちら。
クイックスタート – timeouts
日本語が少しおかしい気がします。
timeout パラメーターに秒数を指定すると、指定した秒数の間、Requestsのレスポンスの待機を止めることができます。
おそらく意図は「timeout パラメーターに秒数を指定すると、指定した秒数でRequestsのレスポンスの待機を止めることができます。」ではないかと。

早速ですがやってみましょう。
timeout になると例外が発生するので、キャッチできるようにします。
(本来はrequestsライブラリで定義されている専用の例外を使うべきなのですが、とりあえずException使います)


import requests

url = "https://httpstat.us/200?sleep=10000"  # 時間のかかるURL

try:
    response = requests.get(url, timeout=3)
    print(response.status_code)  # 実行されない
    print(response.url)  # 実行されない

except Exception as e:
    print(e.args)

# 出力
(ReadTimeoutError("HTTPSConnectionPool(host='httpstat.us', port=443): Read timed out. (read timeout=3)",),)

想定通りに動きました。

requestsを使って、GETでアクセスすると自動的にリダイレクトされる

日常的に使っていて、このブログでも紹介したことのあるrequestsの話です。
参考:requestsを使って、Webサイトのソースコードを取得する

これまで意識せずに使っていたのですが、requestsでgetすると、リダイレクトがあるページの場合、
自動的にリダイレクトされます。
ドキュメントにもはっきり書いてありますね。
リダイレクトと履歴

例えば、このブログはhttpでアクセスすると、httpsのurlにリダイレクトする設定になっています。
そのため、以下のコードは、”https://analytics-note.xyz/” ではなく、
そこからリダイレクトされて、”https://analytics-note.xyz/”にアクセスします。


import requests
url = "https://analytics-note.xyz/"
response = requests.get(url)
print(response.status_code)
print(response.url)

# 以下出力
200
https://analytics-note.xyz/

status_codeがリダイレクトの302ではなく、200になることや、
urlがhttpsの方に書き換えられていることがわかります。

ちなみにリダイレクトされたページへのアクセス結果は、Responseオブジェクトの、historyというプロパティに、
Responseオブジェクトの配列として格納されます。
今回リダイレクトは1回でしたが、複数回に及ぶ可能性もあるので配列で格納されているようです。


print(response.history)
# 出力
[<Response [302]>]

この自動的にリダイレクトしてくれる仕組みはデータ収集等では非常に便利なのですが、
作業の目的によっては逆に不便です。

リダイレクトして欲しくない時は、allow_redirectsという引数にFalseを渡すことでリダイレクトを禁止できます。


response = requests.get(url, allow_redirects=False)
print(response.status_code)
print(response.url)

# 以下出力
302
https://analytics-note.xyz/

graphvizで決定木を可視化

前回の記事でgraphvizをインストールしたので、早速決定木を可視化してみましょう。
サンプルなので、モデル自体は適当に作ります。(データもいつものirisです)


# ライブラリインポート
import graphviz
from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier
from sklearn.tree import export_graphviz
# データの読み込みとモデルの学習
iris = load_iris()
X = iris.data
y = iris.target
clf = DecisionTreeClassifier(min_samples_split=5)
clf.fit(X, y)

これで、モデルができたので可視化してみましょう。
参考にするコードはscikit-learnのサンプルがいいと思います。
ここではblogに挿入することを考慮して、pdfではなくpng画像に出力しました。


dot_data = export_graphviz(
                        clf,
                        class_names=iris.target_names,
                        feature_names=iris.feature_names,
                        filled=True,
                        rounded=True,
                        out_file=None
                    )
graph = graphviz.Source(dot_data)
graph.render("iris-tree", format="png")

結果がこちらです。
わかりやすく可視化できましたね。

macにgraphvizをインストールする

決定木の可視化などに使うgraphvizをMacにインストールする手順のメモです。
ついでにpythonライブラリもインストールします。

まず、本体の方はドキュメントにある通り、homebrewでインストールできます。


brew install graphviz

pythonライブラリの方はこちら。 pipコマンドでインストールします。


pip install graphviz

動作確認をかねて Quickstartのコードを動かしてみましょう。


>>> from graphviz import Digraph
>>> dot = Digraph(comment='The Round Table')
>>> dot  #doctest: +ELLIPSIS
>>> dot.node('A', 'King Arthur')
>>> dot.node('B', 'Sir Bedevere the Wise')
>>> dot.node('L', 'Sir Lancelot the Brave')
>>> dot.edges(['AB', 'AL'])
>>> dot.edge('B', 'L', constraint='false')
>>> print(dot.source)  # doctest: +NORMALIZE_WHITESPACE
>>> dot.render('test-output/round-table.gv', view=True)

これでpdfファイルを開いた時に3個の丸を含むグラフが表示されたらOKです。

(2020/3/21 追記)
このページへのアクセスが非常に多い理由が不明だったのですが、
conda でインストールして、うまく動かなかった人がいるのではないかと思ったので補足します。

condaの方にも、 graphviz と言うパッケージがあり、
conda install graphviz でインストールできますが、
これはgraphvizの本体が含まれないラッパーのみで、
しかも本体を別途 brew でインストールしてもパスが通らないことがあるようです。

conda で導入する場合は
conda install python-graphviz
を使います。

ほかのサイトでは、
brew で本体を入れた後に, graphviz と python-graphviz をそれぞれ入れている例を見かけますが、
試したところ、brewも無しで、
conda install python-graphviz
だけ、行えば動きました。

kerasのモデルの中間層の出力を可視化してみる

ディープラーニングのモデルを作成したとき、中間層の出力が気になることがよくあります。
きちんと活性化しているかとか、相関が高すぎて意味がないユニットが多くないかとか、
どんな条件の時に活性するのかなど、確認したい内容は時により様々です。

kerasの場合、学習済みのモデルの層を取り出して新しいモデルを作成することで中間層の出力を確認できます。

中間レイヤーの出力を得るには?

試しに以前下記の記事で作ったモデルでやってみましょう。
CNNで手書き数字文字の分類

公式ドキュメントに紹介されていたのと少し違う方法ですが、普通にSequentialモデルに学習済みの層を一個追加したら動いたので、
その方法で行います。
一層目には16ユニットあるのですが、そのうち2このユニットについて、出力を可視化しました。

# 学習済みモデルの1層目だけ取得してモデルを作成する


model_2 = Sequential()
model_2.add(model.layers[0])

# 元画像と1層目の出力2個を可視化
fig = plt.figure(figsize=(18, 30))
for i in range(5):
    # print(y_test[i].argmax())
    ax = fig.add_subplot(6, 3, 3*i+1)
    ax.imshow(X_test[i][:, :, 0], cmap='gray_r')
    ax = fig.add_subplot(6, 3, 3*i+2)
    ax.imshow(model_2.predict(X_test[i:i+1])[0][:, :, 0], cmap='gray_r')
    ax = fig.add_subplot(6, 3, 3*i+3)
    ax.imshow(model_2.predict(X_test[i:i+1])[0][:, :, 1], cmap='gray_r')
plt.show()

出力がこちらです。

真ん中の列の出力は横線の下辺に反応していることや、右側の列の結果は中抜き文字のような形で反応しているのがわかりますね。

ちなみに、それぞれのユニットのウェイト(バイアスは除く)を可視化すると次のようになります


fig = plt.figure(figsize=(5,10))
for i in range(2):
    w = model_2.get_weights()[0][:, :, 0, i].reshape(3, 3)
    ax = fig.add_subplot(2, 1, i+1)
    ax.imshow(w, cmap='gray_r')
plt.show()

イメージした通りのウェイトでした。

PrestoのURL関数

Webサービスのアクセスログデータを集計するとき、URLを扱う場面はよくあります。
WHERE句で特定の条件のURLを絞り込んだり、
SELECT句でパラメーターやパスごとに数えたりする場面はよくありますが、
Prestoにはそのようなときに便利な関数が準備されていいます。

ドキュメント 6.18. URL Functions

具体的には下記の関数群です。

  • url_extract_fragment(url)
  • url_extract_host(url)
  • url_extract_parameter(url, name)
  • url_extract_path(url)
  • url_extract_port(url)
  • url_extract_protocol(url)
  • url_extract_query(url)
  • URLを次の形式とすると、それぞれの関数で url_extract_xxxx の xxxx 部分を文字列として返してくれます。
    [protocol:][//host[:port]][path][?query][#fragment]

    これを知るまでは LIKE や正規表現を使うことが多く、稀に予期せぬパターンがマッチしてしまったりしていて困っていたのですが、
    これらのURL関数群を使うようになってからは、そのようなエラーは起こりにくくなりましたし、
    クエリも意図がわかりやすいものになりとても助かっています。

PrestoのCOALESCE関数でNULL値を別の値に置き換える

Prestoの便利な関数の紹介です。
Prestoに限らず、SQLデータを抽出するときにNULL値を別の値に置き換えたいケースは頻繁に発生します。
例えば、NULLならば0、その他の値はそのまま出力したいときは下記のような書き方になります。


CASE
    WHEN column_name1 IS NULL THEN 0
    ELSE column_name1
END AS column_name1

これを、下記のように少しだけスッキリ書くことがでます。


COALESCE(column_name1, 0) AS column_name1

ドキュメントはこちら

この他、次のように引数を複数並べて書くと、左から順番に評価して最初にNULLでなかった列の値を取得できます。
完全外部結合するときなどに便利ですね。


COALESCE(column_name1, column_name2, column_name3, 0) AS column_name1

pythonで編集距離(レーベンシュタイン距離)を求める

ごく稀にですが、文字列同士の編集距離を求める必要が発生するのでその時のメモです。

編集距離(レーベンシュタイン距離)とは、二つの文字列がどの程度異なっているかを表す距離の一種です。
Wikipediaにも解説があります。

一方の文字列に対して、1文字の挿入、削除、置換を最低何回施せばもう一方の文字列に等しくなるかで定まります。

pythonでこれを求めるときは、python-Levenshtein というライブラリが使えます。

インストール


pip install python-Levenshtein

使い方


>>> import Levenshtein
>>> text1 = 'Levenshtein'
>>> text2 = 'Lenvinsten'
>>> Levenshtein.distance(text1, text2)
4

CNNで手書き数字文字の分類

以前の記事で読み込んだ手書き数字文字データ(MINIST)を使って、0~9の数字を判定するモデルを作ってみます。

kerasのサンプルコードもあるのですが、せっかくなので少しだけパラメーターなどを変えてやってみましょう。

最初にライブラリをインポートしてデータを準備します。
前処理として配列の形をConv2Dのinputに合わせるのと、
0〜1への正規化、 ラベルの1hot化を行います。


# ライブラリの読み込み
from keras.datasets import mnist
from keras.utils import to_categorical
from keras.models import Sequential
from keras.layers import Dense, Dropout, Flatten
from keras.layers import Conv2D, MaxPooling2D
from keras.callbacks import EarlyStopping
import matplotlib.pyplot as plt
from sklearn.metrics import classification_report

# データの読み込み
(data_train, target_train), (data_test, target_test) = mnist.load_data()

# Conv2D の inputに合わせて変形
X_train = data_train.reshape(-1, 28, 28, 1)
X_test = data_test.reshape(-1, 28, 28, 1)

# 特徴量を0~1に正規化する
X_train = X_train / 255
X_test = X_test / 255

# ラベルを1 hot 表現に変換
y_train = to_categorical(target_train, 10)
y_test = to_categorical(target_test, 10)

続いてモデルの構築です


# モデルの構築
model = Sequential()
model.add(Conv2D(16, kernel_size=(3, 3),
                 activation='relu',
                 input_shape=(28, 28, 1)))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))
model.add(Conv2D(32, (3, 3), activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))
model.add(Flatten())
model.add(Dense(64, activation='relu'))
model.add(Dropout(0.5))
model.add(Dense(10, activation='softmax'))
model.compile(
    loss="categorical_crossentropy",
    optimizer="adam",
    metrics=['accuracy']
)
print(model.summary())

# 以下、出力
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
conv2d_1 (Conv2D)            (None, 26, 26, 16)        160       
_________________________________________________________________
max_pooling2d_1 (MaxPooling2 (None, 13, 13, 16)        0         
_________________________________________________________________
dropout_1 (Dropout)          (None, 13, 13, 16)        0         
_________________________________________________________________
conv2d_2 (Conv2D)            (None, 11, 11, 32)        4640      
_________________________________________________________________
max_pooling2d_2 (MaxPooling2 (None, 5, 5, 32)          0         
_________________________________________________________________
dropout_2 (Dropout)          (None, 5, 5, 32)          0         
_________________________________________________________________
flatten_1 (Flatten)          (None, 800)               0         
_________________________________________________________________
dense_1 (Dense)              (None, 64)                51264     
_________________________________________________________________
dropout_3 (Dropout)          (None, 64)                0         
_________________________________________________________________
dense_2 (Dense)              (None, 10)                650       
=================================================================
Total params: 56,714
Trainable params: 56,714
Non-trainable params: 0
_________________________________________________________________

そして学習です。


# 学習
early_stopping = EarlyStopping(
                        monitor='val_loss',
                        min_delta=0.0,
                        # patience=2,
                )

history = model.fit(X_train, y_train,
                    batch_size=128,
                    epochs=30,
                    verbose=2,
                    validation_data=(X_test, y_test),
                    callbacks=[early_stopping],
                    )

# 以下出力
Train on 60000 samples, validate on 10000 samples
Epoch 1/30
 - 18s - loss: 0.6387 - acc: 0.7918 - val_loss: 0.1158 - val_acc: 0.9651
Epoch 2/30
 - 18s - loss: 0.2342 - acc: 0.9294 - val_loss: 0.0727 - val_acc: 0.9772
Epoch 3/30
 - 17s - loss: 0.1827 - acc: 0.9464 - val_loss: 0.0571 - val_acc: 0.9815
Epoch 4/30
 - 18s - loss: 0.1541 - acc: 0.9552 - val_loss: 0.0519 - val_acc: 0.9826
Epoch 5/30
 - 18s - loss: 0.1359 - acc: 0.9598 - val_loss: 0.0420 - val_acc: 0.9862
Epoch 6/30
 - 17s - loss: 0.1260 - acc: 0.9620 - val_loss: 0.0392 - val_acc: 0.9880
Epoch 7/30
 - 18s - loss: 0.1157 - acc: 0.9657 - val_loss: 0.0381 - val_acc: 0.9885
Epoch 8/30
 - 19s - loss: 0.1106 - acc: 0.9673 - val_loss: 0.0349 - val_acc: 0.9889
Epoch 9/30
 - 17s - loss: 0.1035 - acc: 0.9694 - val_loss: 0.0359 - val_acc: 0.9885

学習の進み方をプロットしておきましょう。


# Epoch ごとの正解率と損失関数のプロット
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(2, 1, 1, title="loss")
ax.plot(history.epoch, history.history["loss"], label="train_loss")
ax.plot(history.epoch, history.history["val_loss"], linestyle="-.", label="val_loss")
ax.legend()
ax = fig.add_subplot(2, 1, 2, title="acc")
ax.plot(history.epoch, history.history["acc"], label="train_acc")
ax.plot(history.epoch, history.history["val_acc"], linestyle="-.", label="val_acc")
ax.legend()
plt.show()

val_acc は結構高い値を出していますが、一応クラスごとの成績も評価しておきましょう。


# 評価
y_predict = model.predict_classes(X_train)
print(classification_report(target_train, y_predict))

# 以下出力
             precision    recall  f1-score   support

          0       0.99      1.00      0.99      5923
          1       0.99      1.00      0.99      6742
          2       0.99      0.99      0.99      5958
          3       0.99      0.99      0.99      6131
          4       0.99      0.99      0.99      5842
          5       0.99      0.99      0.99      5421
          6       0.99      0.99      0.99      5918
          7       0.98      0.99      0.99      6265
          8       0.99      0.97      0.98      5851
          9       0.98      0.99      0.98      5949

avg / total       0.99      0.99      0.99     60000

ほとんど適当に作ったモデルでしたが、
ほぼほぼ正解できていますね。