はじめに

例によってまたpaizaのBランクレベルアップセット。
一応公式でコード公開オッケーされてる問題なのであしからず。

今回はBランクの本ちゃんのスキルチェックにもあった鼠小僧とほぼ同じ問題。
どっち向いてる時に座標がどう動くかだけ間違わなければ簡単に解ける。

zahyo_muki.py

x,y,n = [int(i) for i in input().split()]

direction = "N"


for i in range(n):
    a = input()

    if direction == "N":
        if a == "R":
            x += 1
            direction = "E"

        else:
            x -= 1
            direction = "W"


    elif direction == "E":
        if a == "R":
            y += 1
            direction = "S"

        else:
            y -= 1
            direction= "N"

    elif direction == "S":
        if a == "R":
            x -= 1
            direction = "W"

        else:
            x += 1
            direction = "E"


    else:
        if a == "R":
            y -= 1
            direction = "N"

        else:
            y += 1
            direction = "S"
    print(x,y)

感想

特になし。
強いて言えば初期位置でどっちの方角向いてるのかが明記されておらず、テストケースの出力で判断するしかなかったのが困った。
この問題解けれるレベル感ならBランクの鼠小僧は解けると思う。

それよりも、、、

この問題の一個手前（これより簡単なはずの位置付け）の問題に手も足も出せず完封された。。。

この時計回りの動きをどうしてもコードに落とし込めず。
とりあえず現実逃避してタイトルの問題に挑んだものの簡単すぎて拍子抜け。
まじでわからん。。。白旗

実行環境

Google Colaboratory

Google Colaboratoryで画像を読み込む為の準備

from google.colab import files
from google.colab import drive
drive.mount('/content/drive')

必要なライブラリの読み込み

import cv2 #opencv
import matplotlib.pyplot as plt 
%matplotlib inline

画像準備

img = plt.imread("/content/drive/My Drive/Colab Notebooks/img/Lenna.bmp")
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)

アダプティブスレッショルド処理

近傍から閾値を求めて変換する。

plt.figure(figsize=(9, 6), dpi=100,
           facecolor='w', linewidth=0, edgecolor='w')

plt.gray()

#オリジナル画像
plt.subplot(2,2,1)
plt.axis('off') 
plt.imshow(gray)

henkango = 255 #閾値を超えたものの数値をどう変換するか
blocksize = 11 #閾値算出の為の近傍領域サイズ(3以降の奇数)
c = 16 #減算値


#アダプティブスレッショナルド処理：MEAN
plt.subplot(2,2,3)
plt.title("MEAN", fontsize=10)
dst = cv2.adaptiveThreshold(gray, henkango, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, blocksize, c)
plt.axis('off') 
plt.imshow(dst)

#アダプティブスレッショナルド処理：GAUSSIAN_C
plt.subplot(2,2,4)
plt.title("GAUSSIAN", fontsize=10)
dst = cv2.adaptiveThreshold(gray, henkango, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, blocksize, c)
plt.axis('off') 
plt.imshow(dst)

Pythonで静的ファイルを使ってみる

PythonでCSSやJsファイルなどのファイルを使ってみます。

Nginxのファイルを設定

まず、マイサーバーはNginx+Gunicornの組み合わせなので、
Nginxの振り分け設定を行います

/staticのアドレスエイリアスを作成します

   location /static {
       alias /app/sampleProj/sampleProj/static;
   }

DjangoのStaticルートを作成

次にDjango内のstaticフォルダを設定します
アプリ内のsetting.pyを編集します

setting.py

STATIC_ROOT = os.path.join(BASE_DIR, 'sampleProj/static')
STATIC_URL = '/static/'

※BASE_DIRはDjangoが設定している定数です。manage.pyの存在するディレクトリになります
ベースディレクトリから、相対パスで、設定できるわけですね

Djangoのstatic設定コマンドを設定

Djangoはデフォルト状態だと、pipによってインストールされたDjangoのインストールフォルダに静的ファイルを設定しています。
よって、現在のフォルダを静的ファイルにするコマンドを打ちます

python manage.py collectstatic

コレで完了！
CSSやJS等の性的ファイルがDjangoに設定されました

オマケ

追加された静的ファイルディレクトリは以下のコマンドで確認できます！

python manage.py findstatic .

参考URL

Django入門｜Nginxでstaticファイルを読み込む設定
https://dot-blog.jp/news/django-nginx-static-load/
BASE_DIRを分かりやすく解説【これで場所と構造がわかります】
https://codor.co.jp/django/about-basedir

Raspberry PiでTensorflowの環境構築【2020年】

情報が少なく、Version管理の難しいRaspberry Piでの環境構築を以下に記します。

1. Raspberry Pi OS (32-bit) LiteをSDに書き込む
   

   僕はmobileNetV2を用いてリアルタイムの画像処理をするのにOSのサイズをできる限り軽量化する必要がありました。
   

2. 自分のパソコンをsshでつなぐ。同じwifi環境下でのみ
   

   Wifiの管理ファイルにご自身のWifi設定を書き加え、sshで接続します。
   

3. 下記のコマンドでOS関係を最新の状態にします。
   sudo apt update
   sudo apt upgrade -y
   sudo reboot

4. この記事を参考にスワップ領域を拡張します。
   

   少しでもメモリ容量を持たせるために拡張します。
   

   テキストファイル開くときおすすめはnanoです。Nanoは使い方は基本的にメモ帳と一緒ですが、セーブがctrl+Xでそのあと保存ファイルの確認が出ますがエンター押してください
   使用例：sudo nano /etc/dphys-swapfile (違うかも)
   ファイルの中身は上のサイトを参考に変更してください。
   

5. いよいよこの記事を参考にTensorflowの環境を構築します。
   上記の記事を参考にインストール
   ただし一部変更があるため以下のコマンドが絶対です。
   必ず記事と見比べながらお願いします。
   ファイルの中身の変更はこの記事を参考に...

   sudo nano /etc/sysctl.conf
   sudo sysctl -p
   sudo nano /etc/rc.local
   sudo reboot　#再起動
   sudo apt install -y libhdf5-dev libqtwebkit4 libqt4-test libatlas-base-dev libjasper-dev
   sudo apt install python3-pip
   sudo apt install python3-dev -y
   
   sudo pip3 install pip -U
   sudo pip3 install setuptools -U
   

   numpy

   sudo pip3 install numpy==1.16.4

   sklearn

   sudo pip3 install scipy
   sudo pip3 install scikit-learn==0.21.3

   matplotlib

   sudo pip3 install matplotlib

   pandas

   sudo pip3 install pandas==0.24.2

   seaborn

   sudo pip3 install seaborn

   Tensorflow

   sudo pip3 install tensorflow==1.14.0

   keras

   sudo pip3 install keras

   flask

   sudo pip3 install flask flask_cors -U

6. 最後にOpenCVのインストールをこの記事を参考に行います。

   
   全部入れた後はpip listで確認してください。

Raspberry PiでTensorflowの環境構築【2020年】

情報が少なく、Version管理の難しいRaspberry Piでの環境構築を以下に記します。

1. Raspberry Pi OS (32-bit) LiteをSDに書き込む
   

   僕はmobileNetV2を用いてリアルタイムの画像処理をするのにOSのサイズをできる限り軽量化する必要がありました。
   

2. 自分のパソコンをsshでつなぐ。同じwifi環境下でのみ
   

   Wifiの管理ファイルにご自身のWifi設定を書き加え、sshで接続します。
   

3. 下記のコマンドでOS関係を最新の状態にします。
   sudo apt update
   sudo apt upgrade -y
   sudo reboot

4. この記事を参考にスワップ領域を拡張します。
   

   少しでもメモリ容量を持たせるために拡張します。
   

   テキストファイル開くときおすすめはnanoです。Nanoは使い方は基本的にメモ帳と一緒ですが、セーブがctrl+Xでそのあと保存ファイルの確認が出ますがエンター押してください
   使用例：sudo nano /etc/dphys-swapfile (違うかも)
   ファイルの中身は上のサイトを参考に変更してください。
   

5. いよいよこの記事を参考にTensorflowの環境を構築します。
   上記の記事を参考にインストール
   ただし一部変更があるため以下のコマンドが絶対です。
   必ず記事と見比べながらお願いします。
   ファイルの中身の変更はこの記事を参考に...

   sudo nano /etc/sysctl.conf
   sudo sysctl -p
   sudo nano /etc/rc.local
   sudo reboot　#再起動
   sudo apt install -y libhdf5-dev libqtwebkit4 libqt4-test libatlas-base-dev libjasper-dev
   sudo apt install python3-pip
   sudo apt install python3-dev -y
   
   sudo pip3 install pip -U
   sudo pip3 install setuptools -U
   

   numpy

   sudo pip3 install numpy==1.16.4

   sklearn

   sudo pip3 install scipy
   sudo pip3 install scikit-learn==0.21.3

   matplotlib

   sudo pip3 install matplotlib

   pandas

   sudo pip3 install pandas==0.24.2

   seaborn

   sudo pip3 install seaborn

   Tensorflow

   sudo pip3 install tensorflow==1.14.0

   keras

   sudo pip3 install keras

   flask

   sudo pip3 install flask flask_cors -U

6. 最後にOpenCVのインストールをこの記事を参考に行います。

   
   全部入れた後はpip listで確認してください。

はじめまして、柴犬のぽん太です

　飼い主がプログラミングを勉強しろとうるさいんです。犬にもプログラミングできるのか挑戦です。ちなみに、Markdownも初めての経験です。よろしくワン！

飼い主の指令によりPython3をインストールします

犬なのにパソコンを使います。キーボードを打つのも大変です。飼い主がMacをくれたのでかじってみましたが、食べれませんでした。プログラミングなんてなんのことやらさっぱりわからないです。途方に暮れていたら、飼い主が「Pythonでもやったら？」とアドバイスをくれました。また、「Python3だよ、間違えるなよ！」と怖い顔をするのでとりあえず飼い主としばらくは目を合わせないようにしたいと思います。
ええと、何の話でしたっけ？
ああ、思い出した！Python3をインストールするのでしたね。

https://www.python.org/downloads

このサイトから、python-3.8.5-macosx10.9.pkgというインストーラーパッケージをダウンロードしました。ダブルクリックでインストール開始しましたワン！

あ、終わりました。

Ponta@dog # python -V
Python 3.8.5
Ponta@dog # 

犬にしては上出来だワン！
きょうはこのくらいにしてあげる！
明日は何しようかな？バイバイ！

要望　加工品のつまり、終了を見抜きたい

機械工作工場にて自動レーンでの製品のつまり、最終加工品の完了した時を知りたい。

バルブの赤いハンドルの位置（開閉表示）を知りたい。

製品は黄色、背景（レーン）は白、青、又は黒の原色。
既製品のセンサーが高すぎるので安く作って欲しいとの要望。

完成品

ラズパイをローカルWIFIにつなげて。
http://"ラズパイのIPアドレス"/home.php

カメラ画像から当該箇所を指定して抜き出す。
色範囲を指定して、色範囲に入った場合GPIOを出力する。
GPIOの出力にて警報ランプを起動させる。
警報解除は警報リセットにて行う。

ラズパイ環境設定

python OpenCV　ダウンロード
参考サイト　計算機科学を学んだソフトウェア技術者の経験やメモなどの記録様

sudo apt-get install libopencv-dev
sudo apt-get install python-opencv

アパッチダウンロード
参考サイト　@mono_taro様
sudo apt-get update
sudo apt-get upgrade
sudo apt-get install apache2
apachectl -v

phpダウンロード
参考サイト　The modern stone age.様
sudo apt install php7.2 libapache2-mod-php7.2

php及びhtml作成

/var/www/html/home.php

<html>
<head>
<meta http-equiv="content-type" charset="utf-8">
<meta http-equiv="Cache-Control" content="no-cache">
<meta http-equiv="refresh" content="5" >
</head>
</html>
<form action="set.php" method="get">
<body>
    <h1>カメラ画像</h1>
    <?php
        $lines_reset = file('/home/pi/Documents/reset.csv');
        foreach($lines_reset as $line_reset){
            $data_reset = explode(',',$line_reset);
        if ($data_reset[0] == 1){
        echo "<FONT COLOR=\"RED\"> 現在警報動作中 解除する場合は下ボタンにて </FONT>";}
        }
    ?>
        <img src="1.jpg?1">
        <h2>抜出画像</h2>
        <img src="2.jpg?1">
</body>
<?php
$lines = file('/home/pi/Documents/color.csv');
foreach($lines as $line){
    $data = explode(',',$line); 
    echo '<p>';
    echo ' 赤色R',$data[0];
    echo ' 緑色G：',$data[1];
    echo ' 青色B：',$data[2];    
    echo '</p>';
}
?>
<a href="./set.html">設定画面</a>
<a href="./reset.html">警報リセット</a>
</form>

/var/www/html/set.html

<form action="set.php" method="get">
<p>
<a href="./home.php">戻る</a>
</p>
<img src="1.jpg">
<p>色の範囲指定</p>
<p>赤下限: <input type="number" name="r_2_1" value="1" min="0" max="255"></p>
<p>赤上限: <input type="number" name="r_2_2" value="1" min="0" max="255"></p>
<p>緑下限: <input type="number" name="g_2_1" value="1" min="0" max="255"></p>
<p>緑上限: <input type="number" name="g_2_2" value="1" min="0" max="255"></p>
<p>青下限: <input type="number" name="b_2_1" value="1" min="0" max="255"></p>
<p>青上限: <input type="number" name="b_2_2" value="1" min="0" max="255"></p>
<p>画像抜出ファイル　上下は上からの距離、左右は左からの距離</p>
<p>上top: <input type="number" name="top" value="1" min="0" max="2000"></p>
<p>下bot: <input type="number" name="bot" value="100" min="0" max="2000"></p>
<p>左left: <input type="number" name="left" value="1" min="0" max="2000"></p>
<p>右right: <input type="number" name="right" value="100" min="0" max="2000"></p>
<p>シャッター間隔、規定回数</p>
<p>シャッター間隔　秒（ｓ）: <input type="number" name="interval" value="60" min="10" max="6000"></p>
<p>回数: <input type="number" name="nutime" value="1" min="1" max="100"></p>
<input type="submit" />
<p>　</p>
<p>補足色説明</p>
<p>赤色　R255 G0 B0</p>
<p>黒　R0 G0 B0</p>
<p>白　R255 G255 B255</p>
<p>黄色　R255 G255 B0</p>
<p>灰色　R128 G128 B128</p>
<p>その他の色はRGB色で検索願います。</p>
<p>規定回数とは連続で指定範囲内に入った回数の時、リレーが作動する回数</p>
</form>

/var/www/html/set.php

<meta http-equiv="content-type" charset="utf-8">
<?php 
     $r_2_1=htmlspecialchars($_GET["r_2_1"],ENT_QUOTES);
     $r_2_2=htmlspecialchars($_GET["r_2_2"],ENT_QUOTES);
     $g_2_1=htmlspecialchars($_GET["g_2_1"],ENT_QUOTES);
     $g_2_2=htmlspecialchars($_GET["g_2_2"],ENT_QUOTES);
     $b_2_1=htmlspecialchars($_GET["b_2_1"],ENT_QUOTES);
     $b_2_2=htmlspecialchars($_GET["b_2_2"],ENT_QUOTES);
     $top=htmlspecialchars($_GET["top"],ENT_QUOTES);
     $bot=htmlspecialchars($_GET["bot"],ENT_QUOTES);
     $left=htmlspecialchars($_GET["left"],ENT_QUOTES);
     $right=htmlspecialchars($_GET["right"],ENT_QUOTES);
     $interval=htmlspecialchars($_GET["interval"],ENT_QUOTES);
     $nutime=htmlspecialchars($_GET["nutime"],ENT_QUOTES);
?>
<body>
<?php echo $_GET["r_2_1"]; ?>
<?php 
$data = [
    [$_GET["r_2_1"], $_GET["r_2_2"], $_GET["g_2_1"], $_GET["g_2_2"], $_GET["b_2_1"], $_GET["b_2_2"],
    $_GET["top"], $_GET["bot"], $_GET["left"], $_GET["right"],
    $_GET["interval"], $_GET["nutime"]],
];
$fp = fopen('/home/pi/Documents/position.csv', 'w');
foreach ($data as $line) {
    fputcsv($fp, $line);
}
fclose($fp);
?>
<p>
<a href="./home.php">戻る</a>
</p>

/var/www/html/reset.html

<meta http-equiv="content-type" charset="utf-8">
<meta http-equiv="Cache-Control" content="no-cache">
<form action="reset.php" method="get">
<input type="submit" />
</form>

/var/www/html/reset.php

<meta http-equiv="content-type" charset="utf-8">
<?php 
$data = [["0"]];
$fp = fopen('/home/pi/Documents/reset.csv', 'w');
foreach ($data as $line) {
    fputcsv($fp, $line);}
fclose($fp);
?>
<p>
<a href="./home.php">戻る</a>
</p>

python 及びCSV作成

/home/pi/Documents/gazou.py

#!/usr/bin/python
# -*- coding: utf-8 -*-
import cv2, os
import numpy as np
import csv
import time
import RPi.GPIO as GPIO
time.sleep(60)#ラズパイ起動時の待機時間
#GPIO設定
PIN = 14
GPIO.setmode(GPIO.BCM)
GPIO.setup(PIN, GPIO.OUT)
#イメージファイルの保存場所
image_file = '/var/www/html/1.jpg' #写真撮影
image_file_2 = '/var/www/html/2.jpg' #検証１
#読み込み用CSVの場所
poss = '/home/pi/Documents/position.csv'
colo_0 = '/home/pi/Documents/color0.csv'
colo_1 = '/home/pi/Documents/color.csv'
poss_reset = '/home/pi/Documents/reset.csv'
while True:
    try :
        with open(poss) as data0:#切り出し場所を読み取る
            reader = csv.reader(data0)
            data1 = [row for row in reader]
        #範囲読み取り
        r_2_1 = int(data1[0][0]) #赤下限 img2
        r_2_2 = int(data1[0][1]) #赤上限 img2
        g_2_1 = int(data1[0][2]) #緑下限 img2
        g_2_2 = int(data1[0][3]) #緑上限 img2
        b_2_1 = int(data1[0][4]) #青下限 img2
        b_2_2 = int(data1[0][5]) #青上限 img2
        #img2.jpgの切り出し部分
        top_2   = int(data1[0][6])
        bot_2   = int(data1[0][7])
        left_2  = int(data1[0][8])
        right_2 = int(data1[0][9])
        #繰り返し秒数及び動作確認数
        jobs_s = int(data1[0][10])
        time.sleep(jobs_s)
        with open(poss_reset) as data_reset:#リセットCSV読み取り
            reader_reset = csv.reader(data_reset)
            data_reset1 = [row for row in reader_reset]
        data_resset2 = int(data_reset1[0][0])
        with open(colo_0) as data2:#書き込み用のデフォルトデータを読み込む
            reader = csv.reader(data2)
            data3 = [row for row in reader]
        cam = cv2.VideoCapture(0)
        # カメラから映像を読み込む
        _, img = cam.read()
        # 画像ファイルとして書き出す
        img2 = img[top_2 : bot_2, left_2: right_2]
        #イメージ保存
        cv2.imwrite(image_file, img)
        cv2.imwrite(image_file_2, img2)
        # RGB平均値を出力
        # flattenで一次元化しmeanで平均を取得 
        r2 = img2.T[2].flatten().mean()
        g2 = img2.T[1].flatten().mean()
        b2 = img2.T[0].flatten().mean()     
        if data_resset2 == 0:
            GPIO.output(PIN, GPIO.LOW) #とりあえず切る
            if r_2_1 < r2 < r_2_2 and g_2_1 < g2 < g_2_2 and b_2_1 < b2 < b_2_2:
                GPIO.output(PIN, GPIO.HIGH)
                data_resset2 = 1
        # RGB平均値を取得
        data3[0][0] = '%.0f' % r2
        data3[0][1] = '%.0f' % g2
        data3[0][2] = '%.0f' % b2
        #CSVを保存
        with open(colo_1, 'w') as f:
            writer = csv.writer(f)
            for data4 in data3:
                writer.writerow(data4)
        data_reset1[0][0] = data_resset2
        with open(poss_reset, 'w') as f:
            writer = csv.writer(f)
            for data_reset2 in data_reset1:
                writer.writerow(data_reset2)
        # 事後処理
        cam.release()
        cv2.destroyAllWindows()
    except:
        time.sleep(30)
        pass

/home/pi/Documents/position.csv

1,1,1,1,1,1,1,100,1,100,60,1

/home/pi/Documents/color0.csv

100,100,100,100,100,100

/home/pi/Documents/color.csv

100,100,100,100,100,100

/home/pi/Documents/reset.csv

0

起動時にpythonを起動する設定をする。

sudo にて gazou.pyを自動起動する。
sudo でないと/var/www/html/にjpgの画像を書き込まない
参考サイト@yuru-camper様

使用したカメラ

C270N HD WEBCAM
logicool様サイト

将来やりたい事

設定画面が下記のように設定しにくいので抜き出し箇所はマウスで選択できるようにしたい。
色設定もマウスで決定できるようにしたい。
現在の状態だと設定がトライアンドエラーをたくさん繰り返さなくてはいけない。

TL;DR

Pythonで3分かかっていた数値計算・約800万行のテキストファイル書き出し処理をRustに置き換えたら約5倍速くなった

はじめに

https://github.com/MIERUNE/japan-mesh-tool

上記はこちらの記事でも紹介した、日本の地域メッシュを.geojsonl形式で生成するPythonモジュールです。標準モジュールのみで書けるくらい純粋な数値計算ですが、たとえば3次メッシュでは800万弱のポリゴン数となり、後述しますが所要時間は3分くらいになります。

このPythonの処理をRustに置き換えたらどれくらい速くなるのか、Rustの勉強を兼ねて試してみました。ちなみにRustのコードも上記のリポジトリに含まれています。

参考：Python上での高速化を検討

アルゴリズムをいじっても劇的な速度改善とはなりませんでした。Python自体インタプリタ言語ですし、速度を求めるべきではないでしょう。また、numpyを使えばC言語の力を借りて速くなるでしょうが、外部モジュールへの依存なしの状態を保ちたかったためボツとしました。

Python vs Rust 速度比較

テスト条件

• macOS: Core i7 9750H RAM:16GB
• 3次メッシュを日本全域で生成する
• Python 3.8.3
• Rust 1.43.0
• Rustはreleaseビルドした実行ファイルで測定する

ファイッ！

Python

time python main.py 3
initializing...
making meshes...
writing file...
done

real    2m46.227s
user    2m31.846s
sys     0m11.249s

166.227秒

Rust

time ./target/release/japan-mesh 3

real    0m35.955s
user    0m34.233s
sys     0m1.272s

35.955秒

結果

35.955/166.227=21.63006%

約5倍の高速化！

Rustについて

なぜRust？

• 速いらしい
• 今アツいらしい

これくらいの軽い理解で始めました。今考えると、何より魅力的なのは、環境構築の容易さだと思います。コマンド一発でコンパイラやパッケージマネージャとか全部入ります。

Rustって難しいんでしょう？

多少かじったくらいですが、Pythonよりは難しいです。ただRustを難しくしているのは所有権システムやライフタイムやらの概念で、今回のようなただの数値計算なら、コンパイラに怒られながらも動くものは作れます。

感想

• アレなコードになってるはずのRustでもちゃんと速いのがうれしい
• 当然ながら計算部分の速度差は歴然だが、ファイル書き出しはそんなに変わらなそう
• Rustはコンパイラが優秀でおかしなコードを書くと怒ってくれる
• いずれもコードに最適化の余地があると思う
• 当初はgeorust/geojsonを使っていたが、to_string()がとても遅かった

はじめに

初投稿になります。プログラミング歴１年半のPHPエンジニアです。

普段の業務ではAWSを使っているのですが、まだ本番環境に関わったことがないので、今回のHerokuへのデプロイを通していろいろと勉強になったことを初心者の視点でまとめていきたいと思います。

今回作ったアプリについて簡単にご紹介します。

• Twitterで特定のワードのリプライに反応して自動で返信

• flaskで自動返信用テキストの追加・編集を行う管理画面を作成

• Twitter Streaming API (Filter realtime Tweets) を利用して、特定のワードのリプライを24時間監視

• BotはPythonのtwitterライブラリを用いて実装

ディレクトリ構成

app/
├── web.py           # flaskアプリの実行ファイル
├── web/             # flaskアプリ
│   ├── models/
│   ├── static/
│   ├── templates/
│   ├── views/
│   ├── __init__.py
│   ├── config.py
│   └── database.py
├── twitter.py       # 自動返信Botの実行ファイル
├── twitter/         # 自動返信Bot（一部省略）
│   ├── __init__.py
│   ├── config.py
│   └── database.py
├── Procfile
├── requirements.txt
├── runtime.txt
├── .env
├── .gitignore
├── migrations/      # Flask-Migrateで作成
└── venv/            # ローカルの仮想環境

flaskについては、以下の記事を参考にさせていただきました。

• Flask + SQLAlchemyプロジェクトを始める手順

• Flaskのディレクトリ構成について考えてみたらBlueprintsを試していた話

Herokuにデプロイする上で必要なファイル

Herokuに環境、必要なライブラリや実行ファイルなどを伝えるためのファイルを準備する必要があるみたいです。

requirements.txt

必要なライブラリをバージョンも含めて記載します。
venvで開発を行ったので、仮想環境にインストールしたパッケージを以下のコマンドでファイルに出力しました。

(venv)$ pip freeze > requirements.txt

インストールしたライブラリ一覧

requirements.txt

alembic==1.4.2
autopep8==1.5.3
cffi==1.14.0
click==7.1.2
cryptography==2.9.2
Flask==1.1.2
Flask-Login==0.5.0
Flask-Migrate==2.5.3
Flask-SQLAlchemy==2.4.3
gunicorn==20.0.4
itsdangerous==1.1.0
Jinja2==2.11.2
Mako==1.1.3
MarkupSafe==1.1.1
pycodestyle==2.6.0
pycparser==2.20
PyMySQL==0.9.3
python-dateutil==2.8.1
python-dotenv==0.14.0
python-editor==1.0.4
six==1.15.0
SQLAlchemy==1.3.18
toml==0.10.1
twitter==1.18.0
Werkzeug==1.0.1

runtime.txt

これは別になくても良いみたいですが、pythonのバージョンを指定するために準備します。
また、Herokuがサポートしているバージョンを指定するよう注意が必要です。

runtime.txt

python-3.7.8

Procfile

アプリの起動方法を指定します。
flaskアプリは本番環境ではgunicornなどのWSGIサーバーを使う必要があります。

Procfile

web: gunicorn web:app --log-file=-
worker: python twitter.py

webでは、web.py中のappというFlaskインスタンスを起動しています。
workerでは、自動返信Botの実行ファイルtwitter.pyを実行しています。

環境変数の設定

ローカル環境と本番環境で異なる環境変数（ex. DB情報）をコードを書き換えずにどう設定するかにつまづいたので、まとめます。

Herokuの環境変数

以下のコマンドで環境変数を設定できます。
（コマンドを使用するためには、Heroku CLI のインストールが必要です）

$ heroku config:set DB_HOST=xxxxxxxxxxx

環境変数の一覧を表示して、設定できているか確認できます。

$ heroku config

MySQL DB の環境変数の設定は、こちらを参考にさせていただきました。
・HerokuでMySQLを利用する方法

ローカル環境の環境変数

.envファイルを用意して、そこに環境変数を記載していきます。
機密情報なので、gitignoreするのを忘れないようにしましょう。

.env

ENV = 'LOCAL'

# DB
DB_HOST = 'xxxxx'
DB_NAME = 'xxxxx'
DB_USER = 'xxxxx'
DB_PASSWORD = 'xxxxx'

# Session
SESSION_SECRET_KEY = 'xxxxx'

アプリ側で環境変数を読み込む

app/web/config.py で以下のように環境変数を読み込んで、アプリ側の設定を行います。

config.py

"""FlaskのConfigを提供する"""
import os
from os.path import join, dirname
from dotenv import load_dotenv

dotenv_path = join(dirname(__file__), '../.env')
load_dotenv(dotenv_path)


class Config:

    # Flask
    if (os.environ.get('ENV') == 'LOCAL'):
        DEBUG = True
    else:
        DEBUG = False

    # Session
    SECRET_KEY = os.environ.get('SESSION_SECRET_KEY')

    # SQLAlchemy
    SQLALCHEMY_DATABASE_URI = 'mysql+pymysql://{user}:{password}@{host}/{db_name}?charset=utf8'.format(**{
        'user': os.environ.get('DB_USER'),
        'password': os.environ.get('DB_PASSWORD'),
        'host': os.environ.get('DB_HOST'),
        'db_name': os.environ.get('DB_NAME')
    })
    SQLALCHEMY_TRACK_MODIFICATIONS = False
    SQLALCHEMY_ECHO = False


Config = Config

このように設定すると、Heroku環境でもローカル環境でも環境変数を読み込むことができます。
Pythonにおける環境変数の扱いについては、こちらを参考にさせていただきました。
・【GitHub】に載せたくない環境変数の書き方 Python

1アカウントにつき1日1リプライをどう実装するか

短時間に何度も自動返信するのは、Twitter API的によろしくないので、1アカウントにつき1日1リプライまでという仕様にしました。

一度リプライしたユーザをどうやって管理するか迷ったんですけど、

また、Dynosは、Herokuで実行されているアプリケーションの正常性を維持するために、少なくとも1日1回再起動（循環）されます。ローカルファイルシステムへの変更はすべて削除されます。サイクリングは24時間ごとに発生します（さらに、アプリケーションのすべてのdynoが同時に再起動しないように、最大​​216のランダムな分）。

Herokuの自動再起動が24時間周期で行われるみたいなので、今回はworkerに指定したPythonファイル上（常時ループし続けている）のリストでリプライしたユーザを管理することにしました。
こうすることで、自動再起動されるたびにリストは初期化されるので、24時間ごとに再び自動返信を行うことができるようになります。

twitter.py

replied_user_list = [] # 返信したユーザをリストで管理

# 24時間監視し続ける（プログラムが走り続ける）
for tweet in twitter_stream.statuses.filter(language='ja', track=tracking_text):
    # 以下、特定のワード（tracking_text）が含まれるツイートを検知したときの処理
    # ...
    # 自動返信を行うと、リストに追加する
    replied_user_list.append(user_id)
    # ...

webがsleepするとworkerもsleepする問題

Heroku の無料プラン (Free Dyno) では、30分間Webアプリにアクセスがないと自動的にsleepします。
これに関しては、全く問題がないのですが、どうやらwebがsleepするとworkerもsleepしてしまうみたいです。[https://devcenter.heroku.com/articles/free-dyno-hours#dyno-sleeping]

Apps that only utilise a free worker dyno do not sleep, because they do not respond to web requests. Be mindful of this as they may run 24/7 and consume from your pool of hours.

workerがsleepしてしまうと、自動返信Botが反応しなくなってしまうので、これはかなり問題です。

そこで、webをsleepさせないために以下のようにリクエストを定期的に送るなどの対処法があるみたいですが、
・Herokuの無料dynoをスリープさせないで24時間稼働させる4つの方法

これだと、webとworkerのdynoがともに24時間稼働し続けることになってしまい、無料のdyno時間である1000時間/月をオーバーしてしまいます。

結果、webは必要なとき以外は停止しておくことにしました。（Heroku ダッシュボードで簡単に起動できる）
worker単体ではsleepすることがないので、これで自動返信Botが24時間稼働し続けられるようになりました。

まとめ

pythonでのWebアプリ開発は初めてで、そういうことも書こうと思ったんですが、今回はHerokuへデプロイする上でつまづいたことを中心に書きました。

コーディング以外の知識をもっと身に付ける必要があるなと痛感しました。
今後はAWSの勉強をして、AWSソリューションアーキテクトのアソシエイトを取りたいと思います。

初めての技術ブログですので、読みにくい等の技術以外のアドバイスもどんどんよろしくお願いします。もちろん、技術的な指摘もよろしくお願いします。

使用するライブラリ

Requests

Requestsは、 Python の HTTP 通信ライブラリです。 Requests を使うとWebサイトの情報取得や画像の収集などを簡単に行うことができます。

Beautiful Soup

HTML や XML から狙ったデータを抽出するためのライブラリです。
公式ドキュメントの冒頭の説明を見るとこれは HTML や XML のパーサーそのものではなく、パーサーをラップして扱いやすくするライブラリのようです。

データの格納

まずrequestsをインポート
欲しいURLを取得する

qiita.py

import requests

def get_rakuten():
    url = "https://search.rakuten.co.jp/search/mall/pythonチュートリアル/"
    responce = requests.get(url)
    print(responce)

get_rakuten()

実行結果

HTTPステータスコードは、HTTPにおいてWebサーバからのレスポンスの意味を表現する3桁の数字からなるコード
200の場合、リクエストは成功し、レスポンスとともに要求に応じた情報が返される。

今回の場合取得したいデータがresponceに格納できたことを示す。

<Response [200]>

取得したいページのhtmlを取得

.textとするとhtmlの情報が取得できる

qiita.py

import requests

def get_rakuten():
    url = "https://search.rakuten.co.jp/search/mall/pythonチュートリアル/"
    responce = requests.get(url)
    print(responce.text)

get_rakuten()

実行結果の一部

<!DOCTYPE html>
<html>
<head>
  <meta http-equiv="X-UA-Compatible" content="IE=edge">
  <meta charset="utf-8">
  <meta name="viewport" content="width=device-width, initial-scale=1.0">

    <script> window.__INITIAL_STATE__ = {};
    window.__REWIRED_SCHEMAS__ = ["pc"];</script>
    <title data-duiid="7q54vmlfe">【楽天市場】pythonチュートリアルの通販</title><meta name="description" content="楽天市場-「pythonチュートリアル」22件 人気の商品を価格比較・ランキング･レビュー・口コミで検討できます。ご購入でポイント取得がお得。セール商品・送料無料商品も多数。「あす楽」なら翌日お届けも可能です。"/><meta name="keywords" content="pythonチュートリアル,通販,インターネット,ショッピング,楽天市場"/><meta name="robots" content="noodp"/><meta property="og:image" content="https://r.r10s.jp/com/img/logo/rakuten_group_red.svg"/><meta property="dui:trace" content="b9cb1ccc-0db5-4869-b005-516a67cf13c0"/><link rel="canonical" href="https://search.rakuten.co.jp/search/mall/python%E3%83%81%E3%83%A5%E3%83%BC%E3%83%88%E3%83%AA%E3%82%A2%E3%83%AB/"/><link rel="alternate" href="android-app://jp.co.rakuten.android/rakuten-ichiba-app/www.rakuten.co.jp/search?source=appindexing_search&amp;genreId=0&amp;keyword=python%E3%83%81%E3%83%A5%E3%83%BC%E3%83%88%E3%83%AA%E3%82%A2%E3%83%AB"/><link rel="icon" href="https://www.rakuten.co.jp/favicon.ico"/>
  <script type="text/javascript" src="https://r.r10s.jp/com/js/d/shisa/shisa-1.0.2.min.js" id="script-shisa"
          data-group="ichiba-search"></script>

.textで取得したhtml情報を解析する
その際 BeautifulSoupを使用する

qiita.py

import requests
from bs4 import BeautifulSoup

def get_rakuten():
    url = "https://search.rakuten.co.jp/search/mall/pythonチュートリアル/"
    responce = requests.get(url)
    html = responce.text
    soup = BeautifulSoup(html, "html.parser")
    print(soup)

get_rakuten()

実行結果の一部
divクラスが整理されているのがわかる

</div>
<div class="item">
<div class="dui-container conditionfilter open" data-duiid="cif84ng3s"><div class="dui-filter-menu"><div class="header" title="コンディション"><span>コンディション</span></div><div class="content"><div class="dui-list"><a class="item" href="https://search.rakuten.co.jp/search/mall/python%E3%83%81%E3%83%A5%E3%83%BC%E3%83%88%E3%83%AA%E3%82%A2%E3%83%AB/?used=0">新品</a><a class="item" href="https://search.rakuten.co.jp/search/mall/python%E3%83%81%E3%83%A5%E3%83%BC%E3%83%88%E3%83%AA%E3%82%A2%E3%83%AB/?used=1">中古品</a><a class="item" href="https://search.rakuten.co.jp/search/mall/python%E3%83%81%E3%83%A5%E3%83%BC%E3%83%88%E3%83%AA%E3%82%A2%E3%83%AB/-/tg5002023/">レンタル</a></div></div></div></div>
</div>

URLから商品タイトルを取得

デベロッパーツールで確認するとtitleタグに商品名があったので.selectで取得したいタグを指定
.selectではtitleとして見つかった複数のものがリスト型で出力されるのでfor文で回す

qiita.py

import requests
from bs4 import BeautifulSoup

def get_rakuten():
    url = "https://search.rakuten.co.jp/search/mall/pythonチュートリアル/"
    responce = requests.get(url)
    html = responce.text
    soup = BeautifulSoup(html, "html.parser")
    titles = soup.select(".title")
    for title in titles:
        print(title.text)

get_rakuten()

出力結果の一部

中古 【中古】 Excel　2010　＆　SQL　Server　2008　R2による企業データ分 PDCAサイクルに沿った最新BIの実 / / [単行本（ソフトカバー）]【ネコポス発送】
中古 【中古】 Vocaloidをたのしもう vol．2 / ヤマハミュージックメディア / ヤマハミュージックメディア [ムック]【ネコポス発送】
中古 【中古】 mule／viスーパーリファレンス / 伊藤 真人, 田谷 文彦, 小巻 賢二郎, 前田 雄一郎 / ソフトバンククリエイティブ [単行本]【ネコポス発送】

お届け先で設定された都道府県（離島除く）への最も安い配送方法での送料が各商品に表示されます。



翌日配送（あす楽）対応あり



楽天プレミアム特典あり

余分な空白が出力されたので除去

.searchresultitemの変数をリスト型で格納し、その中でtitleを一つ一つ出力していく
select_oneとすれば見つかった最初の情報のみが保持されるためリスト型にはならない

qiita.py

import requests
from bs4 import BeautifulSoup

def get_rakuten():
    url = "https://search.rakuten.co.jp/search/mall/pythonチュートリアル/"
    responce = requests.get(url)
    html = responce.text
    soup = BeautifulSoup(html, "html.parser")
    items = soup.select(".searchresultitem")

    for item in items:
        title = item.select_one(".title")
        print(title.text)

get_rakuten()

実行結果の一部
空白が消えていることがわかる

中古 【中古】 mule／viスーパーリファレンス / 伊藤 真人, 田谷 文彦, 小巻 賢二郎, 前田 雄一郎 / ソフトバンククリエイティブ [単行本]【宅配便出荷】
中古 【中古】 Excel　2010　＆　SQL　Server　2008　R2による企業データ分 PDCAサイクルに沿った最新BIの実 / / [単行本（ソフトカバー）]【ネコポス発送】
中古 【中古】 Vocaloidをたのしもう vol．2 / ヤマハミュージックメディア / ヤマハミュージックメディア [ムック]【ネコポス発送】
中古 【中古】 mule／viスーパーリファレンス / 伊藤 真人, 田谷 文彦, 小巻 賢二郎, 前田 雄一郎 / ソフトバンククリエイティブ [単行本]【ネコポス発送】

価格を表示

価格はimportantセレクタに記載されているため
select_oneで取得できる

qiita.py

import requests
from bs4 import BeautifulSoup

def get_rakuten():
    url = "https://search.rakuten.co.jp/search/mall/pythonチュートリアル/"
    responce = requests.get(url)
    html = responce.text
    soup = BeautifulSoup(html, "html.parser")
    items = soup.select(".searchresultitem")

    for item in items:
        title = item.select_one(".title")
        print(title.text)
        price = item.select_one(".important")
        print(price.text)

get_rakuten()

実行結果の一部

Pythonチュートリアル第3版 Python　3．5対応 [ グイド・ファン・ロッサム ]
1,980円
Pythonチュートリアル／GuidovanRossum／鴨澤眞夫【1000円以上送料無料】
1,980円
中古 【中古】 Pythonチュートリアル　第3版 ／Guido　van　Rossum(著者),鴨澤眞夫(訳者) 【中古】afb
1,870円

titleを.textとしてテキスト化し、replaceメソッドを使えるようにして改行をなくす
さらに、価格を取得したあとは改行を入れて見やすくする

qiita.py

import requests
from bs4 import BeautifulSoup

def get_rakuten():
    url = "https://search.rakuten.co.jp/search/mall/pythonチュートリアル/"
    responce = requests.get(url)
    html = responce.text
    soup = BeautifulSoup(html, "html.parser")
    items = soup.select(".searchresultitem")

    for item in items:
        title = item.select_one(".title").text.replace("\n", "")
        price = item.select_one(".important").text
        print(title)
        print(price + "\n")

get_rakuten()

商品データを入力, 値段の安い順に表示

?s=2をURLの語尾につけると価格の安い順に表示
?f=2の場合、送料無料の物のみ表示
?f=2&s=2、なら送料無料、価格の安い順にソートした結果が返される

qiita.py

import requests
from bs4 import BeautifulSoup

keyward = input("比較したい商品を入力してください:\n")

def get_rakuten():
    url = "https://search.rakuten.co.jp/search/mall/" + keyward + ?f=2&s=2
    responce = requests.get(url)
    html = responce.text
    soup = BeautifulSoup(html, "html.parser")
    items = soup.select(".searchresultitem")

    for item in items:
        title = item.select_one(".title").text.replace("\n", "")
        price = item.select_one(".important").text
        print(title)
        print(price + "\n")

get_rakuten()

課題

スクレイピングした結果が欲しい商品ではない可能性がある
そのため自分で、選ぶ必要がある
item_numberに初期値１を与えて、for文で回して1加算
さらに、リスト型に取得した価格を格納して選べるようにする

1,192円をint型にするためにreplaceで『,』と『円』を消した

qiita.py

import requests
from bs4 import BeautifulSoup

keyward = input("比較したい商品を入力してください:\n")

def get_rakuten():
    url = "https://search.rakuten.co.jp/search/mall/" + keyward + "?f=2&s=2"
    responce = requests.get(url)
    html = responce.text
    soup = BeautifulSoup(html, "html.parser")
    items = soup.select(".searchresultitem")

    item_number = 1
    price_list = []
    for item in items:
        title = item.select_one(".title").text.replace("\n", "")
        price = item.select_one(".important").text.replace(",","").replace("円","")
        price_list.append(price)
        print(item_number)
        print(title)
        print(price + "\n")
        item_number+=1

    selected_item_number = int(input("楽天：商品番号を入力してください：\n"))
    selected_price = int(price_list[selected_item_number-1])
    return selected_price

rakuten_price = get_rakuten()
print(rakuten_price)

蛍光X線プロファイル

最近は，蛍光X線の中心エネルギー，強度比，タイプ，自然幅は，xraylibを使えば，簡単に取得できる．しかし，実際上は，金属の効果でラインプロファイルが広がる効果を加味する必要があり，例えば，MnのKalpha線はKalpha1,2の強度比が1:2で終わりのはずが，7つや8つのVoigt関数をいれてフィットしないと，正しく装置の分解能を評価できない．いくつか代表的なラインプロファイルのパラメータと簡単な生成方法をpythonでまとめておいた．

実測スペクトルをモデルフィットする場合は，pythonで複数のvoigt関数にレスポンスを入れてフィットする方法 を参照してほしい．

コード

Lineprofile クラスが，ラインのエネルギー，強度，自然幅，装置の分解能を加味して，Voigt関数を生成する．mymodel クラスの中では，それをrawfuncという関数にして，必要なら低エネルギー側のテイルの計算も行えるようにしている．

plot_manylines_wletail_qiita.py

#!/usr/bin/env python
__version__= '1.0'

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'serif'
import numpy as np
import scipy.special

def mymodel(x,params, consts=[], tailonly = False):
    norm,gw,gain,P_tailfrac,P_tailtau,bkg1,bkg2 = params    
    # norm : nomarlizaion
    # gw : sigma of gaussian
    # gain : gain of the spectrum 
    # P_tailfrac : fraction of tail 
    # P_tailtau : width of the low energy tail
    # bkg1 : constant of background
    # bkg2 : linearity of background    
    initparams = [norm,gw,gain,bkg1,bkg2]
    def rawfunc(x): # local function, updated when mymodel is called 
        return Lineprofile(x,initparams,consts=consts)               
    model_y = smear(rawfunc, x, P_tailfrac, P_tailtau, tailonly=tailonly)
    return model_y

def Lineprofile(xval,params,consts=[]):
    norm,gw,gain,bkg1,bkg2 = params
    # norm : normalization 
    # gw : sigma of the gaussian 
    # gain : if gain changes
    # consttant facter if needed 
    prob = (amp * lgamma) / np.sum(amp * lgamma) # probabilites for each lines. 
    model_y = 0 
    if len(consts) == 0:
        consts = np.ones(len(energy))
    else:
        consts = consts

    for i, (ene,lg,pr,con) in enumerate(zip(energy,lgamma,prob,consts)):
        voi = voigt(xval,[ene*gain,lg*0.5,gw])
        model_y += norm * con * pr * voi

    background = bkg1 * np.ones(len(xval)) + (xval - np.mean(xval)) * bkg2
    model_y = model_y + background
    # print "bkg1,bkg2 = ", bkg1,bkg2, background
    return model_y

def voigt(xval,params):
    center,lw,gw = params
    # center : center of Lorentzian line
    # lw : HWFM of Lorentzian (half-width at half-maximum (HWHM))
    # gw : sigma of the gaussian 
    z = (xval - center + 1j*lw)/(gw * np.sqrt(2.0))
    w = scipy.special.wofz(z)
    model_y = (w.real)/(gw * np.sqrt(2.0*np.pi))
    return model_y

def smear(rawfunc, x, P_tailfrac, P_tailtau, tailonly = False):
    if P_tailfrac <= 1e-5:
        return rawfunc(x)

    dx = x[1] - x[0]
    freq = np.fft.rfftfreq(len(x), d=dx)
    rawspectrum = rawfunc(x)
    ft = np.fft.rfft(rawspectrum)
    if tailonly:
        ft *= P_tailfrac * (1.0 / (1 - 2j * np.pi * freq * P_tailtau) - 0)
    else:
        ft += ft * P_tailfrac * (1.0 / (1 - 2j * np.pi * freq * P_tailtau) - 1)

    smoothspectrum = np.fft.irfft(ft, n=len(x))
    if tailonly:
        pass
    else:
        smoothspectrum[smoothspectrum < 0] = 0
    return smoothspectrum    

class aline:
    def __init__(self,x,y,name):
        self.x = x
        self.y = y
        self.name = name

# global variables
gfwhm = 2
gw = gfwhm / 2.35
norm = 500000.0
gain = 1.0
bkg1 = 1.0
bkg2 = 0.0
P_tailfrac = 1e-6 # no tail 
P_tailtau = 10
nbin=1000
ewidth=500
init_params=[norm,gw,gain,P_tailfrac,P_tailtau,bkg1,bkg2]

linelist = []


#################################################################
name="Ti Kalpha"
energy = np.array((4510.918, 4509.954, 4507.763, 4514.002, 4504.910, 4503.088))
lgamma = np.array((1.37, 2.22, 3.75, 1.70, 1.88, 4.49))
amp = np.array((4549, 626, 236, 143, 2034, 54))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Ti KBeta"
energy = np.array((25.37, 30.096, 31.967, 35.59)) + 4900
lgamma = np.array((16.3, 4.25, 0.42, 0.47))
amp = np.array((199, 455, 326, 19.2))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="V Kalpha"
energy = np.array((4952.237, 4950.656, 4948.266, 4955.269, 4944.672, 4943.014))
lgamma = np.array((1.45, 2.00, 1.81, 1.76, 2.94, 3.09))
amp = np.array((25832, 5410, 1536, 956, 12971, 603))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="V KBeta"
energy = np.array((18.19, 24.50, 26.992)) + 5400
lgamma = np.array((18.86, 5.48, 2.499))
amp = np.array((258, 236, 507))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Cr Kalpha"
energy = 5400 + np.array([14.874, 14.099, 12.745, 10.583, 18.304, 5.551, 3.986])
lgamma = np.array([1.457, 1.760, 3.138, 5.149, 1.988, 2.224, 4.4740])
amp = np.array([882, 237, 85, 45, 15, 386, 36])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Cr KBeta"
energy = 5900 + np.array((47.00, 35.31, 46.24, 42.04, 44.93))
lgamma = np.array([1.70, 15.98, 1.90, 6.69, 3.37])
amp = np.array([670, 55, 337, 82, 151])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Mn Kalpha"
energy = 5800 + np.array((98.853, 97.867, 94.829, 96.532, 99.417, 102.712, 87.743, 86.495))
lgamma = np.array([1.715, 2.043, 4.499, 2.663, 0.969, 1.553, 2.361, 4.216])
amp = np.array([790, 264, 68, 96, 71, 10, 372, 100])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Mn KBeta"
energy = 6400 + np.array((90.89, 86.31, 77.73, 90.06, 88.83))
lgamma = np.array((1.83, 9.40, 13.22, 1.81, 2.81))
amp = np.array([608, 109, 77, 397, 176])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Fe Kalpha"
energy = np.array((6404.148, 6403.295, 6400.653, 6402.077, 6391.190, 6389.106, 6390.275))
lgamma = np.array((1.613, 1.965, 4.833, 2.803, 2.487, 2.339, 4.433))
amp = np.array([697, 376, 88, 136, 339, 60, 102])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Fe Kbeta"
energy = np.array((7046.90, 7057.21, 7058.36, 7054.75))
lgamma = np.array((14.17, 3.12, 1.97, 6.38))
amp = np.array([107, 448, 615, 141])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Co Kalpha"
energy = np.array((6930.425, 6929.388, 6927.676, 6930.941, 6915.713, 6914.659, 6913.078))
lgamma = np.array((1.795, 2.695, 4.555, 0.808, 2.406, 2.773, 4.463))
amp = np.array((809, 205, 107, 41, 314, 131, 43))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Co Kbeta"
energy = np.array((7649.60, 7647.83, 7639.87, 7645.49, 7636.21, 7654.13))
lgamma = np.array((3.05, 3.58, 9.78, 4.89, 13.59, 3.79))
amp = np.array((798, 286, 85, 114, 33, 35))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Ni Kalpha"
energy = np.array((7478.281, 7476.529, 7461.131, 7459.874, 7458.029))
lgamma = np.array((2.013, 4.711, 2.674, 3.039, 4.476))
amp = np.array((909, 136, 351, 79, 24))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Ni Kbeta"
energy = np.array((8265.01, 8263.01, 8256.67, 8268.70))
lgamma = np.array((3.76, 4.34, 13.70, 5.18))
amp = np.array((722, 358, 89, 104))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Cu Kalpha"
energy = np.array([8047.8372, 8045.3672, 8027.9935, 8026.5041])
lgamma = np.array([2.285, 3.358, 2.667, 3.571])
amp = np.array([957, 90, 334, 111])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Cu KBeta"
energy = np.array([8905.532, 8903.109, 8908.462, 8897.387, 8911.39])
lgamma = np.array([3.52, 3.52, 3.55, 8.08, 5.31])
amp = np.array([757, 388, 171, 68, 55])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="As Kalpha"
energy = np.array((10543.2674,10507.50))
lgamma = np.array((3.08, 3.17))
amp = np.array((1.00, 0.51))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="As KBeta"
energy = np.array((11725.73,11719.86))
lgamma = np.array((2.09+2.25, 2.09+2.25))
amp = np.array((0.13, 0.06))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


plt.figure(figsize=(10,8))
plt.title("Line Profiles")

for oneline in linelist:

    plt.xlabel("Energy (eV)")
    plt.plot(oneline.x, oneline.y, '-', label = oneline.name)
    plt.legend(numpoints=1, frameon=False, loc="best")
    plt.grid(linestyle='dotted',alpha=0.5)

plt.savefig("linelist.png")
plt.show()

使い方

単純に実行すると，このような図が生成される．

もし，装置の低エネルギー側のテイルを入れたい場合は，P_tailfrac = 1e-6 # no tail の部分を大きな数字にすればよい．

energy，lgamma，amp，という global 変数を用意して，それを書き換えるという邪道な書き方をしているが，正しくは，データベースとなる情報はそもそもコードと一体管理しないほうが良いし，ラインごとにちゃんとしたオブジェクトを生成した方がよい．

• データサイエンス100本ノックをやりきる会を作りました?
• こちらのSlack招待URLからご参加ください！！
• ぜひ一緒に励まし合いながら、データサイエンス100本ノックをやり切りたいと思っています！

Youtube

動画解説もしています。

問題

P-031: レシート明細データフレーム（df_receipt）に対し、店舗コード（store_cd）ごとに売上金額（amount）の標本標準偏差を計算し、降順でTOP5を表示せよ。

解答

コード

df_receipt.groupby('store_cd').amount.std(ddof=0).reset_index().sort_values('amount', ascending=False).head(5)

出力

	store_cd	amount
28	S13052	663.391816
31	S14011	553.456916
42	S14034	544.903736
5	S13001	543.536561
12	S13015	543.409938

解説

・PandasのDataFrame/Seriesにて、同じ値を持つデータをまとめて処理し、同じ値を持つデータの合計や平均などを確認したい時に使用します。
・'groupby'は、同じ値や文字列を持つデータをまとめて、それぞれの同じ値や文字列に対して、共通の操作（合計や平均など）を行いたい時に使います。
・'<列名>.std(ddof=0)'は、標本標準偏差を算出する際に'N-ddof(=0)'で割ることを指します。
・'.reset_index()'は、'groupby'によってバラバラになったインデックス番号を0始まりの連番に振り直す操作を行いたい時に使います。
・'.sort_values('amount', ascending=False)'で'amount'を降順に表示しています。**

※こちらのコードでは不偏標準偏差を出力します。不偏標準偏差の場合は、'ddof=1'となります。

コード

df_receipt.groupby('store_cd').agg({'amount':'std'}).reset_index().sort_values('amount', ascending=False).head(5)

活性化関数のトレンドとして最近SwishやMishが有名です。
Geluに関して説明している日本語の記事がないので自分なりの理解をまとめたいと思います。
これに関する説明がおそらくSwishの理解につながると思われます。
論文：https://arxiv.org/abs/1606.08415

Gelu定義

Gelu = x \Phi (x) = x * \frac{1}{2}[1+erf(\frac{x}{\sqrt{2}})]

Geluは以上のように定義されます。
ここで$\Phi (x)$は正規分布（ガウス分布）の累積分布関数です。$erf()$は誤差関数です。
これが理想的な関数ですが一般に誤差関数は初等関数で計算できないので近似が行われます。

\begin{align}
Gelu &= 0.5x(1+tanh[\sqrt{\frac{2}{\pi}}(x+0.044715x^3)])\\
Gelu &= x\sigma (1.702x)\\
Swish &= x\sigma (\beta x)\\
\end{align}

このGelu、およびGeluの近似をプロットしてみて、relu関数との差分を考えてみます。
ここでGelu idealとGelu approximate 1は大体等しいです。
一方、Gelu approximate 2はGelu idealと若干誤差があることが分かります。
ここで$\beta = 1.702$のSwish関数はGelu approximate 2と等価であることが分かります。
従ってGeluが有利である理由を掴めれば、それがSwishが有利な理由になります。
なお、この論文中ではβ=1のSwish関数はSiLUと呼称しています。

ただし$\beta = 1.702$のSwish関数とGelu approximate 2と等しいとはいえ、厳密なGelu関数はSwish関数とは異なることに注意する必要があります。

gelu.py

import matplotlib.pyplot as plt
import numpy as np
from scipy import special

x = np.arange(-6, 6, 0.05)
y1 = np.array([i if i>0 else 0 for i in x])                  # Relu
y2 = x*0.5*(1+special.erf(np.sqrt(0.5)*x))                   # Gelu ideal
y3 = x*0.5*(1+np.tanh(np.sqrt(2/np.pi)*(x+0.044715*x*x*x)))  # Gelu approximate 1
y4 = x/(1+np.exp(-1.702*x))                                  # Gelu approximate 2

plt.plot(x, y1, label="Relu")
plt.plot(x, y2, label="Gelu ideal")
plt.plot(x, y3, label="Gelu approximate 1")
plt.plot(x, y4, label="Gelu approximate 2")
plt.legend()
plt.show()

plt.plot(x, y2-y1, label="Gelu ideal - Relu")
plt.plot(x, y3-y1, label="Gelu approximate 1 - Relu")
plt.plot(x, y4-y1, label="Gelu approximate 2 - Relu")
plt.legend()
plt.show()

Geluが有利である理由の推測

Gelu-Reluのグラフにもう一度注目します。この時Gelu-Reluは二つの谷があるグラフになります。そして勾配が0になるのは$x=±1$付近です。またGeluはGAUSSIAN ERROR LINEAR UNITSの略です。
以上から推測を行うと、この活性化関数は関数の入力を1か-1に近づける勾配成分を持つのではないかと思われます。(いわゆる正則化を推進する項が活性化関数に含まれる)

activity_regularizer(中間層の出力に関する正則化)のL1,L2正則化との比較を考えてみたい。この場合の追加する損失関数成分は$\lambda |x|, \lambda x^2$である。
この損失関数を加えると出力はどの値でもゼロの方向に近づくが、出力が小さくなりすぎると逆に追加の損失関数よりも本来の損失関数が増大するので、出力の大きさは結局どこかで釣り合うようになる。
Gelu-Reluの差分を追加の損失関数と見なせば出力がゼロに近すぎるとき、逆にこの関数は出力を大きくしようとする成分を持つ。これはactivity_regularizerのL1,L2正則化にはない項である。

Mishとの比較

Geluと$\beta = 1$のSwishとMishのReluの差分を比べたところ以下の通りでした。
Mishは正と負で谷の大きさと場所が異なります。
Mishは負の領域で$\beta = 1.0$のSwishに近いですが、正の領域でGelu approximate 2($\beta = 1.702$のSwish)に近い気がします。
（ついでに何故、Mishが負の領域の方に偏っている理由を妄想しますと、負の値の収束値を大きくした方が正則化した際に正の値を取る入力の割合が大きくなるからとか考えました。）

gelu.py

y5 = x/(1+np.exp(-1.0*x))                                    # Swish (beta=1.0)
y6 = x*np.tanh(np.log(1+np.exp(x)))                          # Mish

plt.plot(x, y2-y1, label="Gelu ideal - Relu")
plt.plot(x, y4-y1, label="Gelu approximate 2 - Relu")
plt.plot(x, y5-y1, label="Swish (beta=1.0) - Relu")
plt.plot(x, y6-y1, label="Mish - Relu")
plt.legend()
plt.show()

まとめ：

GeluをRelu関数からの変分成分を表示させた。
Relu関数からの変分という解釈でプロットする記事が見当たらないので良い解釈かどうか分からないが、この場合Geluは入力を1か-1に近づける勾配成分を持つ、というのが自分の理解である。
これはMish論文にself regularized(自己正則化)という単語で現れるがそれではないかと思います。

L1正則化とLeakyReluの誤差逆伝播における影響を比べてみました。
L1正則化は機械学習では不要な説明変数をそぎ落とす次元圧縮のために用いるという説明がなされます。一方で、LeakyReluはxが負の時も僅かな傾きをもたせることで多層でも学習が止まってしまうことを防ぐ活性化関数です。
何故こんな全く関係ないものを比較してしまったのか、自分でも意味が分かりません。

通常の全結合の誤差逆伝播

x1 = Input(shape=(124,))
x2 = Dense(100)(x1)
x3 = Activation('relu')(x2)
y = Dense(10)(x3)

model.compile(loss='mean_squared_error')

\begin{align}
x_1 &= input\\
x_2 &= Ax_1 + B\\
x_3 &= relu(x_2)\\
y &= Cx_3+D\\
L_{mse} &= \frac{1}{2}(t-y)^2\\
\end{align}

損失関数、活性化関数、全結合計算が上記のように行われたとする。
この場合、各偏微分式は

\begin{align}
\frac{\partial L_{mse}}{\partial y} &= -(t - y)\\
\frac{\partial y}{\partial C} &= x_3\\
\frac{\partial y}{\partial D} &= 1\\
\frac{\partial y}{\partial x_3} &= C\\
\frac{\partial x_3}{\partial x_2} &= H(x_2)\\
 H(x)&=\left\{
\begin{array}{ll}
1 & (x \geq 0) \\
0 & (x \lt 0)
\end{array}
\right.\\
\frac{\partial x_2}{\partial A} &= x_1\\
\frac{\partial x_2}{\partial B} &= 1\\
\end{align}

からモデル重み$A,B,C,D$の更新量は連鎖律を使って

\begin{align}
\frac{\partial L_{mse}}{\partial A} &= \frac{\partial L_{mse}}{\partial y}\frac{\partial y}{\partial x_3}\frac{\partial x_3}{\partial x_2}\frac{\partial x_2}{\partial A}\\
&=-(t - y)\cdot C \cdot H(x_2) \cdot x_1\\ 
\frac{\partial L_{mse}}{\partial B} &= \frac{\partial L_{mse}}{\partial y}\frac{\partial y}{\partial x_3}\frac{\partial x_3}{\partial x_2}\frac{\partial x_2}{\partial B}\\
&=-(t - y)\cdot C \cdot H(x_2) \\ 
\frac{\partial L_{mse}}{\partial C} &= \frac{\partial L_{mse}}{\partial y}\frac{\partial y}{\partial C}\\
&=-(t - y)\cdot x_3 \\
\frac{\partial L_{mse}}{\partial D} &= \frac{\partial L_{mse}}{\partial y}\frac{\partial y}{\partial D}\\
&=-(t - y)\\
\end{align}

と示すことができます。

L1正則化の全結合層の誤差逆伝播

Kerasにおいて正則化は以下の三つを持ちます。
・kernel_regularizer
・bias_regularizer
・activity_regularizer
https://keras.io/ja/regularizers/

これらを前節における一回目の全結合層$x_2 = Ax_1 + B$での係数でいえばkernel_regularizerは重み$A$の、bias_regularizerは重み$B$の、activity_regularizerは出力$x_2$の正則化をそれぞれ行う事になります。
この場合、L1正則化とは損失関数に追加で小さい係数$\lambda$を掛けた変数の絶対値を追加します。
それぞれ$\lambda |A|$、$\lambda |B|$、$\lambda |x_2|$。
一方、L2正則化では損失関数に追加で小さい係数$\lambda$を掛けた変数の二乗を追加します。
それぞれ$\lambda A^2$、$\lambda B^2$、$\lambda x_2^2$。

ここではL1正則化のactivity_regularizer、$\lambda |x_2|$について考えます。
追加の損失関数を下記のように考えれば$|x|$の微分は$x>0$で$1$、$x<0$で$-1$だから

L_{L1} = \lambda |x_2| \\
\frac{\partial |x|}{\partial x} = \left\{
\begin{array}{ll}
1 & (x \geq 0) \\
-1 & (x \lt 0)
\end{array}
\right.\\
 H(x)=\left\{
\begin{array}{ll}
1 & (x \geq 0) \\
0 & (x \lt 0)
\end{array}
\right.\\
\frac{\partial L_{L1}}{\partial x_2} = \lambda(2H(x_2)-1)\\

従ってモデル重み$A,B,C,D$の更新量のL1正則化のactivity_regularizerによる変化分は

\begin{align}
\frac{\partial L_{L1}}{\partial A} &= \frac{\partial L_{L1}}{\partial x_2}\frac{\partial x_2}{\partial A}\\
&=\lambda(2H(x_2)-1) \cdot x_1\\ 
\frac{\partial L_{L1}}{\partial B} &= \frac{\partial L_{L1}}{\partial x_2}\frac{\partial x_2}{\partial B}\\
&=\lambda(2H(x_2)-1) \\ 
\frac{\partial L_{L1}}{\partial C} &= 0\\
\frac{\partial L_{L1}}{\partial D} &= 0\\
\end{align}

なお、L1正則化が次元圧縮のために用いられるという説明はactivity_regularizerではなくkernel_regularizerの$L_{L1} = \lambda |A|$のL1正則化損失項を考えた場合の説明です。

x1 = Input(shape=(124,))
x2 = Dense(100, activity_regularizer=regularizers.l1(0.01))(x1)
x3 = Activation('relu')(x2)
y = Dense(10)(x3)

model.compile(loss='mean_squared_error')

LeakyReluの全結合層の誤差逆伝播

LeakyReluは以下の様な活性化関数です。
一般に$\alpha ≪1$で

 LeakyRelu(x)=\left\{
\begin{array}{ll}
x & (x \geq 0) \\
\alpha x & (x \lt 0)
\end{array}
\right.\\
\frac{\partial LeakyRelu(x)}{\partial x}
 =\left\{
\begin{array}{ll}
1 & (x \geq 0) \\
\alpha  & (x \lt 0)
\end{array}
\right.\\

から

 H(x)=\left\{
\begin{array}{ll}
1 & (x \geq 0) \\
0 & (x \lt 0)
\end{array}
\right.\\

を使ってLeakyReluと通常のRelu勾配差分は

\begin{align}
\frac{\partial LeakyRelu(x)}{\partial x} &=(1-\alpha )H(x) + \alpha \\
\frac{\partial LeakyRelu(x)}{\partial x}-\frac{\partial Relu(x)}{\partial x} 
 &=((1-\alpha )H(x) + \alpha )-H(x)\\
&=-\alpha (H(x) -1)\\
\end{align}

と書けます。
従ってモデル重み$A,B,C,D$の更新量のRelu=>LeakyReluによる変化分は

\begin{align}
\frac{\partial L_{mse}}{\partial A} &=(t - y)\cdot C \cdot \alpha (H(x_2)-1) \cdot x_1\\ 
\frac{\partial L_{mse}}{\partial B} &=(t - y)\cdot C \cdot \alpha (H(x_2)-1) \\ 
\frac{\partial L_{mse}}{\partial C} &= 0\\
\frac{\partial L_{mse}}{\partial D} &= 0\\
\end{align}

と示すことができます。
またこれをあえて変形させると以下の様になります。
ここで二項目は通常の誤差逆伝播に$\alpha $を掛けた分に過ぎませんので以降無視します。

\begin{align}
\frac{\partial L_{mse}}{\partial A} &=(t - y)\cdot C \cdot \alpha (2H(x_2)-1) \cdot x_1-[\alpha (t - y)\cdot C \cdot H(x_2) \cdot x_1]\\ 
\frac{\partial L_{mse}}{\partial B} &=(t - y)\cdot C \cdot \alpha (2H(x_2)-1)-[\alpha (t - y)\cdot C \cdot H(x_2)] \\ 
\frac{\partial L_{mse}}{\partial C} &= 0\\
\frac{\partial L_{mse}}{\partial D} &= 0\\
\end{align}

x1 = Input(shape=(124,))
x2 = Dense(100)(x1)
x3 = LeakyReLU(alpha=0.01)(x2)
y = Dense(10)(x3)

model.compile(loss='mean_squared_error')

比較

通常の逆伝播からの変化分を比較してみます。
ここで$\lambda≪1,\alpha ≪1$で共に1より十分小さな値です。

• L1正則化のactivity_regularizerによる変化分

\begin{align}
\frac{\partial L_{L1}}{\partial A} &=\lambda(2H(x_2)-1) \cdot x_1\\ 
\frac{\partial L_{L1}}{\partial B} &=\lambda(2H(x_2)-1) \\ 
\frac{\partial L_{L1}}{\partial C} &= 0\\
\frac{\partial L_{L1}}{\partial D} &= 0\\
\end{align}

• Relu=>LeakyReluによる変化分

\begin{align}
\frac{\partial L_{mse}}{\partial A} &=(t - y)\cdot C \cdot \alpha (2H(x_2)-1) \cdot x_1\\ 
\frac{\partial L_{mse}}{\partial B} &=(t - y)\cdot C \cdot \alpha (2H(x_2)-1) \\ 
\frac{\partial L_{mse}}{\partial C} &= 0\\
\frac{\partial L_{mse}}{\partial D} &= 0\\
\end{align}

ここでRelu=>LeakyReluによる変化分は$(t - y)\cdot C$がゼロと見なせない学習の初期においてactivity_regularizerのL1正則化による変化分と近いように見なせないでしょうか。
一方、学習が十分進み、$(t - y)\cdot C$がゼロに近づいてしまうとL1正則化のactivity_regularizerが弱まるように見なせないでしょうか。……などと考えてみました。
ただし、L1正則化のactivity_regularizerの変化分はあくまで正則化を掛けた層にしか関係しませんが、Relu=>LeakyReluによる変化分は活性化関数を掛けた以前の層にも伝達します。

また、Relu=>LeakyReluによる変化分の二項目のみの場合、通常の誤差逆伝播に$\alpha $を掛けた分になります。

• Relu=>LeakyReluによる二項目の変化分

\begin{align}
\frac{\partial L_{mse}}{\partial A} &=-\alpha (t - y)\cdot C \cdot H(x_2) \cdot x_1\\ 
\frac{\partial L_{mse}}{\partial B} &=-\alpha (t - y)\cdot C \cdot H(x_2) \\ 
\frac{\partial L_{mse}}{\partial C} &= 0\\
\frac{\partial L_{mse}}{\partial D} &= 0\\
\end{align}

元の勾配と合わせLeakyReluの活性化関数を掛けた以前の全部の勾配が$(1+\alpha )$ 倍に大きくなります。
モデルが深い場合において、すべてのRelu関数をLeakyRelu関数に置き換えると入力に近いほうの浅い層の勾配が$(1+\alpha )^n$の対数的に大きくできるのが示唆されます。

まとめ：

L1正則化とLeakyReluを比較しました。
Relu=>LeakyReluによる変化分はactivity_regularizerのL1正則化による変化と、活性化関数以前の勾配を$(1+\alpha )$ 倍にするように思いました。
この解釈はLeakyRelを以下の様な関数の和であるとします。

\begin{align}
 LeakyRelu(x)&=\left\{
\begin{array}{ll}
x & (x \geq 0) \\
\alpha x & (x \lt 0)
\end{array}
\right.\\
&=\left\{
\begin{array}{ll}
(1+\alpha )x & (x \geq 0) \\
0 & (x \lt 0)
\end{array}
\right.  - \alpha|x|\\
&=(1+\alpha )Relu(x)- \alpha|x|
\end{align}

はじめに

時系列データでAIモデルを構築する際、スライディングウィンドウによる特徴量計算を行った上で各種アルゴリズムで学習を行う、という処理フローはよく実施されます。この時系列データの特徴量抽出ですが、メモリに乗りきらないような容量の大きいデータ（いわゆるビッグデータ）では、pandasやnumpyなどの便利なパッケージが使えないため処理の難易度が上がります。このような場合、どのように実現していますか？
本記事ではその一例としてPySparkを用いる方法をご紹介します。文量が多くなってしまいますので、２回に分けてご紹介していきます。

今回は基礎編ということで、一般的な「PySparkでスライディングウィンドウによる特徴量計算する方法」を紹介します。

検証環境

PySparkの実行環境としてAzure Synapse Analyticsを使いました。
主要パッケージのバージョンは以下の通りです。（2020年8月20日時点のデフォルト設定）

Apache Spark 2.4
Python version 3.6.1

PySparkでスライディングウィンドウによる特徴量計算する方法

検証環境作成方法は割愛させていただきます。
Azure Synapse AnalyticsでどのようにSpark実行環境を構築するかについては、ご要望がありましたら記事にしたいと思います。コメントでご要望いただけると幸いです。

1. データの準備

適当なデータをPySparkのデータフレームとして定義します。

df = sqlContext.createDataFrame([
    (1, 2.65,2.42,6.90,4.93),
    (2, 2.57,8.50,2.40,5.37),
    (3, 2.13,3.76,7.52,7.67),
    (4, 3.09,7.28,3.59,6.34),
    (5, 5.75,4.69,5.26,3.11),
    (6, 6.91,4.04,2.03,6.28),
    (7, 5.44,3.22,2.87,7.14),
    (8, 4.86,7.47,3.68,0.32),
    (9, 9.70,7.43,4.43,7.74),
    (10,6.30,7.72,7.78,7.91),
    ], 
    ["time", "data1", "data2", "data3", "data4"])

df.show()

# +----+-----+-----+-----+-----+
# |time|data1|data2|data3|data4|
# +----+-----+-----+-----+-----+
# |   1| 2.65| 2.42|  6.9| 4.93|
# |   2| 2.57|  8.5|  2.4| 5.37|
# |   3| 2.13| 3.76| 7.52| 7.67|
# |   4| 3.09| 7.28| 3.59| 6.34|
# |   5| 5.75| 4.69| 5.26| 3.11|
# |   6| 6.91| 4.04| 2.03| 6.28|
# |   7| 5.44| 3.22| 2.87| 7.14|
# |   8| 4.86| 7.47| 3.68| 0.32|
# |   9|  9.7| 7.43| 4.43| 7.74|
# |  10|  6.3| 7.72| 7.78| 7.91|
# +----+-----+-----+-----+-----+

仮に以下のようなデータが取れたものと考えてください。

カラム名	意味
time	記録時間（秒）
data1~data6	計測データ

2. スライディングウィンドウの定義

PySparkのスライディングウィンドウはpyspark.sql.windowのWindowsで定義します。ここではウィンドウ幅５秒のスライディングウィンドウを定義しています。

from pyspark.sql.window import Window
import pyspark.sql.functions as F

# sliding-windowの設定
window_size = 5
sliding_window = Window.orderBy(F.col("time")).rowsBetween(Window.currentRow, window_size-1)

orderBy("カラム名")でソートキーを指定しています。Sparkでは処理順序が保証されないため、ソートキーを指定することはとても重要です。本例では記録時間を示すtimeの順、つまりtimeを昇順に並べた上で１レコード目から順に処理して欲しいため、orderBy(F.col("time"))という指定をしています。ちなみにデフォルトではASC(昇順)で処理されます。DESC(降順)で処理したい場合は以下のように記述します。

sliding_window = Window.orderBy(F.col("time").desc()).rowsBetween(Window.currentRow, window_size-1)

F.col("time")に.desc()を付けると降順扱いとなります。

次にウィンドウ幅をrowsBetween(Window.currentRow, window_size-1)で定義しています。第１引数は開始位置の定義であり、ここではWindow.currentRowと現在の行を指定しています。第２引数は終了位置の定義であり、ここではwindow_size-1と現在の行から４行先（４秒先）を指定しています。これで現在の行を含む４行先までの５行分（５秒分）のデータを一つのウィンドウとして定義できたことになります。

3. 特徴量計算

先ほど設定したスライディングウィンドウの定義を使った特徴量抽出を実施します。data1に対してウィンドウ幅内のmax(最大値)、min(最小値)、avg(平均値)を取得してみます。

df.withColumn('feat_max_data1', F.max('data1').over(sliding_window))\
  .withColumn('feat_min_data1', F.min('data1').over(sliding_window))\
  .withColumn('feat_avg_data1', F.avg('data1').over(sliding_window))\
  .select('time', 'data1', 'feat_max_data1', 'feat_min_data1', 'feat_avg_data1')\
  .show()

# +----+-----+--------------+--------------+------------------+
# |time|data1|feat_max_data1|feat_min_data1|    feat_avg_data1|
# +----+-----+--------------+--------------+------------------+
# |   1| 2.65|          5.75|          2.13|3.2379999999999995|
# |   2| 2.57|          6.91|          2.13|              4.09|
# |   3| 2.13|          6.91|          2.13|             4.664|
# |   4| 3.09|          6.91|          3.09|              5.21|
# |   5| 5.75|           9.7|          4.86| 6.531999999999999|
# |   6| 6.91|           9.7|          4.86|             6.642|
# |   7| 5.44|           9.7|          4.86|             6.575|
# |   8| 4.86|           9.7|          4.86| 6.953333333333333|
# |   9|  9.7|           9.7|           6.3|               8.0|
# |  10|  6.3|           6.3|           6.3|               6.3|
# +----+-----+--------------+--------------+------------------+

withColumn("カラム名", 処理内容)でデータフレームの新規カラムとして、指定したカラム名で指定した処理内容の結果が追加されます。maxを計算する処理コードを見るとwithColumn('feat_max_data1', F.max('data1').over(sliding_window))となっており、「data1のmaxを取ります。over(sliding_window)の条件で。結果はfeat_max_data1カラムとして追加します。」という解釈になります。PySparkでのスライディングウィンドウ指定はover()で定義します。
PySparkでは１つずつ処理を定義していくため、この例の様に１つのカラムから複数の特徴量を取得する場合は複数の処理コードを羅列していく必要があります。

まとめ

以上、「PySparkでスライディングウィンドウによる特徴量計算する方法」の基礎編を紹介しました。
今回紹介した処理方法は一般的な手法であり、データが少ない場合や抽出する特徴量が少ない場合は、十分事足りると思います。しかし、データ量が多い場合、処理するカラムが多い場合、抽出する特徴量が多い場合など、この処理方法では処理効率が悪く処理コストが大きくなってしまいます。ですが、処理方法を工夫することで処理コストを劇的に改善することが可能です。
次回は応用編として、どのような工夫を行えばより効率的に処理を行えるのかをご紹介します。

ご一読いただき、ありがとうございました。

以下からどうぞ
↓↓↓↓↓↓
PYTHON製WEBフレームワークDJANGOのススメ

TL;DR

np.clipで最小値（最大値）を指定する！

実行環境

• python: 3.7.7
• numpy: 1.19.1

なにがおきるの?

pythonのfloat型ではおおよそ1e-323未満の値は0.0になってしまいます．
pythonのfloat型は多くの環境では,-1.7976931348623157e+308から1.7976931348623157e+308の間で範囲を取りうる．

$ python
>>> a = 1e-323
>>> b = 1e-324
>>> print(a)
　1e-323
>>> print(b)
　0.0

bが0.0になってしまうため，このbに対してlogをとってしまうとエラーがでて-infになってしまいます．
また，0割の警告も出てきます．

>>> import numpy as np
>>> b = 1e-324
>>> np.log(b)
　__main__:1: RuntimeWarning: divide by zero encountered in log
　-inf

解決策

np.clipで最小値(最大値)を指定してあげる．

>>> import numpy as np
>>> b = 1e-324
>>> np.log(np.clip(a=b, a_min=1e-323, a_max=1e+10))
　-743.7469247408213

最小値(最大値)を指定してあげることで，floatの範囲外の値が入った時でも一応計算できるようになります．ただ，最小値（最大値）を超えた値はそれぞれ設定した値に丸め込まれてしまうのでその分誤差が大きくなってしまうデメリットもあります．

昨夜は、【データサイエンティスト入門】科学計算、データ加工、グラフ描画ライブラリの使い方の基礎としてNumpyの基礎をまとめましたが、今夜はScipyの基礎をまとめます。ただし、Scipyは無暗に写経しても理解は難しいのと膨大なので、本書に乗った解説を補足することとします。
【注意】
「東京大学のデータサイエンティスト育成講座」を読んで、ちょっと疑問を持ったところや有用だと感じた部分をまとめて行こうと思う。
したがって、あらすじはまんまになると思うが内容は本書とは関係ないと思って読んでいただきたい。

Scipyの基礎

Scipyは、科学技術計算をするためのライブラリで、多様な数学処理（線形代数の計算、フーリエ変換など）が出来ます。
ここでは、線形代数の逆行列や固有値、方程式の解などを求めてみましょう。
※固有値問題はnumpyでも出来ましたが計算速度はScipyの方が速いです

2-3-1 Scipyライブラリのインポート

>>> import scipy as sp
>>> import scipy.linalg as linalg
>>> from scipy.optimize import minimize_scalar

2-3-2 行列計算

2-3-2-1 行列式と逆行列の計算

>>> matrix = np.array([[1,-1,-1],[-1,1,-1],[-1,-1,1]])
>>> matrix
array([[ 1, -1, -1],
       [-1,  1, -1],
       [-1, -1,  1]])
>>> print('行列式', linalg.det(matrix))
行列式 -4.0
>>> print('逆行列','\n', linalg.inv(matrix))
逆行列
 [[ 0.  -0.5 -0.5]
 [-0.5 -0.  -0.5]
 [-0.5 -0.5  0. ]]
>>> print('行列の積','\n',matrix.dot(linalg.inv(matrix)))
行列の積
 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
>>> matrix*linalg.inv(matrix)   #numpyの要素積
array([[ 0. ,  0.5,  0.5],
       [ 0.5, -0. ,  0.5],
       [ 0.5,  0.5,  0. ]])
>>> matrix@linalg.inv(matrix)   #numpyの積
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

2-3-2-2 固有値と固有ベクトル

>>> eig_value, eig_vector = linalg.eig(matrix)
>>> print('固有値','\n',eig_value)
固有値
 [-1.+0.j  2.+0.j  2.+0.j]
>>> print('固有ベクトル','\n',eig_vector)
固有ベクトル
 [[ 0.57735027 -0.81649658  0.42755853]
 [ 0.57735027  0.40824829 -0.81618716]
 [ 0.57735027  0.40824829  0.38862863]]
>>> eig_value, eig_vector = np.linalg.eig(matrix)  #numpyでも同じ
>>> print('固有ベクトル','\n',eig_vector)
固有ベクトル
 [[ 0.57735027 -0.81649658  0.42755853]
 [ 0.57735027  0.40824829 -0.81618716]
 [ 0.57735027  0.40824829  0.38862863]]
>>> print('固有値','\n',eig_value)   #固有値は実数表記
固有値
 [-1.  2.  2.]

scipy.linalg vs numpy.linalgについて
【参考】
Linear Algebra (scipy.linalg)
「scipy.linalgには、numpy.linalgのすべての関数が含まれています。加えて、numpy.linalgに含まれていない他のいくつかのより高度なものが含まれています」
「numpy.linalgよりもscipy.linalgを使用するもう1つの利点は、常にBLAS / LAPACKサポートを使用してコンパイルされるのに対し、numpyの場合はオプションです。したがって、numpyのインストール方法によっては、scipyのバージョンの方が速い場合があります。」
「したがって、numpyプログラムへの依存関係としてscipyを追加したくない場合を除いて、numpy.linalgの代わりにscipy.linalgを使用します。」

2-3-3 ニュートン法

最適化計算の例

2-3-3-1 方程式の解

>>> def my_function(x):
...     return (x**2 + 2*x + 1)
...
>>> from scipy.optimize import newton
>>> newton(my_function,0)
-0.9999999852953906

2-3-3-2 最小値を求める

Brent法を利用して、my_function()の最小値を求める。

>>> minimize_scalar(my_function, method = 'Brent')
     fun: 0.0
    nfev: 9
     nit: 4
 success: True
       x: -1.0000000000000002

練習問題2-6

>>> def my_function2(x):
...     return (x**3 + 2*x + 1)
...
>>> newton(my_function2,0)
-0.45339765151640365

SciPy Tutorial

ちょっと、これだけだと単にほんとの入り口の使い方しか書けていないので、以下そのうち、写経したいと思います。
SciPy Tutorial
Introduction
Basic functions
Special functions (scipy.special)
Integration (scipy.integrate)
Optimization (scipy.optimize)
Interpolation (scipy.interpolate)
Fourier Transforms (scipy.fft)
Signal Processing (scipy.signal)
Linear Algebra (scipy.linalg)
Sparse eigenvalue problems with ARPACK
Compressed Sparse Graph Routines (scipy.sparse.csgraph)
Spatial data structures and algorithms (scipy.spatial)
Statistics (scipy.stats)
Multidimensional image processing (scipy.ndimage)
File IO (scipy.io)

まとめ

・本書のScipyの基礎をまとめた
・Scipyはnumpyと被るものも多いが、基本Scipyの方が高速であるようだ

・SciPy Tutorialはボリュームはあるが面白い例が多く勉強になりそうなのでまとめたいと思う