yukicoder contest 278 参戦記

A 1335 1337

書くだけ.

A = input()
S = input()

result = ''
for c in S:
    if c in '0123456789':
        result += A[int(c)]
    else:
        result += c
print(result)

B 1336 Union on Blackboard

問題文を見た瞬間に★1.5!?ってなった. 2 3 0 や 2 3 1 を手計算してみたところ、どの順序でも同じ数値になったので、順番に関係ないのかと踏んで前から順に処理するコードを書いて AC.

T = int(input())

m = 1000000007

for _ in range(T):
    N = int(input())
    A = list(map(int, input().split()))
    t = 0
    for a in A:
        t = t * a + (t + a)
        t %= m
    print(t)

C 1337 Fair Otoshidama

最初に100円も持っているんだったら数字は自由に変えれそうなので、結局 X + Y + Z を3人で分けて同じ数にできるかだけだよなということでこうなった.

X, Y, Z = map(int, input().split())

if (X + Y + Z) % 3 == 0:
    print('Yes')
else:
    print('No')

D 1338 Giant Class

長さ W の配列を初期値 H に初期化して、Y - 1 の最小値を取り、毎回減った数を H * W から引いていけば答えが求まる……と思ったら W≦10⁹ で死んだ(爆). 辞書で記録するようにして AC.

from sys import stdin

readline = stdin.readline

H, W, Q = map(int, readline().split())

c = H * W
t = {}
result = []
for _ in range(Q):
    Y, X = map(int, readline().split())
    t.setdefault(X - 1, H)
    p = t[X - 1]
    t[X - 1] = min(t[X - 1], Y - 1)
    c -= p - t[X - 1]
    result.append(c)
print(*result, sep='\n')

はじめに

Twitterで一時期流行していた 100 Days Of Code なるものを先日知りました。本記事は、初学者である私が100日の学習を通してどの程度成長できるか記録を残すこと、アウトプットすることを目的とします。誤っている点、読みにくい点多々あると思います。ご指摘いただけると幸いです！

今回学習する教材

• Effective Python

  • 8章構成
  • 本章216ページ

今日の進捗

• 進行状況：87-91ページ
• 第4章：メタクラスと属性
• 本日学んだことの中で、よく忘れるところ、知らなかったところを書いていきます。

getやsetメソッドよりも素のままの属性を使う

getやsetメソッドを定義すると、計算をするときにぎこちないコードになってしまいます。

class OldResistor(object):
    def __init__(self, ohms):
        self._ohms = ohms

    def get_ohms(self):
        return self._ohms

    def set_ohms(self, ohms):
        self._ohms = ohms

r0 = OldResistor(100)
r0.set_ohms(r0.get_ohms() + 500)
print(r0.get_ohms())
# 600

それよりも素の属性を使ったほうが計算時にすっきりしたコードになります。

class Resistor(object):
    def __init__(self, ohms):
        self.ohms = ohms
        self.voltage = 0
        self.current = 0

r1 = Resistor(50e3)
r1.ohms = 400
r1.ohms += 200
print(r1.ohms)
# 600

後に、属性が設定された時に特別な振る舞いが必要になる場合は、@propertyデコレータとそれに対応する、setter属性を付与すれば良いです。
また、正しく動作するためには、セッターメソッドとゲッターメソッドの両方の名前が、意図しているプロパティ名と一致しないとダメです。

class VoltageResistance(Resistor):
    def __init__(self, ohms):    # 初期化
        super().__init__(ohms)   # ohmsを引数として、親クラスの__init__を呼び出す
        self._voltage = 0        # インスタンス変数 _voltageを宣言

    @property                    # _voltageに入っている値を返す
    def voltage(self): 
        return self._voltage

    @voltage.setter              # voltageに新たな値を代入するたびにセッターメソッドが実行される
    def voltage(self, voltage):  # プロパティ名と一致
        self._voltage = voltage  
        self.current = self._voltage / self.ohms

r2 = VoltageResistance(100)
print(r2.current)
# 0
r2.voltage = 10                 # voltage.setterが呼ばれ、voltageに10を渡す
print(r2.current)
# 0.1

このように、素の属性を使い、@propertyを使うことですっきりしたコードを書きながら機能を拡張することができます。@propertyメソッドを使ってセッターやゲッターを実装する際に、可読性を落とさないために以下のポイントを抑えるべしです。
関連するオブジェクト状態だけを変更し、オブジェクトを超えて、モジュールを動的にインポートする、遅いヘルパー庵数を実行するなどの、呼び出し元が予期しない動作をしないように注意！
また、複雑であったり遅くなるような場合は通常のメソッドを使って実装すべきなようです。

はじめに

データ分析をする際に初めにやることといえば、データがどのような特徴を持っているのか把握することからだと思います。そんなときにpandas-profilingを使うとEDAを一括して行ってくれるのでとても便利です。しかし、Jupyternotebookで試したところデータの日本語のカラムなどが□□□のように文字化け(豆腐)となってしまったので解決策をまとめていきたいと思います。

環境

• Python3.8
• Docker
• Jupyter notebook

原因

pandas-profilingで文字化けが起きる原因はmatplotlibとseabornが日本語化対応できていないことに起因します。日本語化対応できれば、matplotlibとseabornを利用したpandas-profilingも日本語対応されます。

この記事ではmatplotlibとseabornの日本語化の手順について説明することになります。

Matplotlibとseabornの日本語化

この記事ではDockerを利用したJupyternotebook環境の日本語対応について説明します。
さらに効率の良い方法があるかもしれませんので、その際にはコメントを頂けると嬉しいです。

1. 日本語のフォントをダウンロードする

こちらのサイトからipaexg00401.zip(4.0MB)をダウンロードして解凍します。
ipaexg00401フォルダの中にあるipaexg.ttfをDockerfileがあるディレクトリに移動します。

2. seabornの日本語対応のためコンテナ上のファイルをホストにコピーする

ここで行う作業は、seabornを日本語化させるために必要なrcmod.pyをローカルにダウンロードして内容を書き換え、docker-compose upの度に、コンテナ上のrcmod.pyをホストの書き換え済みのrcmod.pyで上書きするように設定します。このような流れをとることでdocker-compose upの度にrcmod.pyを書き換えずに済みます。

(本当はDockerfileでコンテナ上のを書き換えたいのですがわかりませんでした)

日本語対応していない状態でdocker-compose upをします。
別のターミナルを開いて、コンテナIDを確認します。

# コンテナIDを確認する
$ docker ps

次にコンテナ上のrcmod.pyをホスト(ローカル上)に保存します。

$ docker cp [コンテナID]:opt/conda/lib/python3.8/site-packages/seaborn/rcmod.py [保存先(C:\Users\....など)]

最後に保存したrcmod.pyをDockerfileがあるディレクトリにコピーします。

3. rcmod.pyを書き換える

rcmod.pyを開いて以下の内容を変更します。

86-87行目のdef set(context="notebook", ...)のfontの部分をfont="IPAexGothic"に変更します。

def set_theme(context="notebook", style="darkgrid", palette="deep",
              font="IPAexGothic", font_scale=1, color_codes=True, rc=None):

次に205行目の"font.family": ["sans-serif"]を以下に変更します

"font.family": ["IPAexGothic"]

これでseabornの日本語対応のための書き換えは終了です。

4. Dockerfileに以下を追加する

# matplotlibとscipyの日本語化
# 日本語のフォントをコピーする
COPY ipaexg.ttf /opt/conda/lib/python3.8/site-packages/matplotlib/mpl-data/fonts/ttf/ipaexg.ttf
# 書き換えたrcmod.pyでコンテナ上のrcmod.pyを上書きする
COPY settings/localize_ja/rcmod.py /opt/conda/lib/python3.8/site-packages/seaborn/rcmod.py
# matplotlibの設定ファイルの最後にfont.family : IPAexGothicを追加する
RUN echo "font.family : IPAexGothic" >>  /opt/conda/lib/python3.8/site-packages/matplotlib/mpl-data/matplotlibrc
# キャッシュを削除する
RUN rm -r ./.cache

これでmatplotlibとseabornの日本語対応が行えます。
以下のコードで文字化けしないか確かめることができます。

# matplotlibが日本対応できているか確認する
import matplotlib.pyplot as plt
plt.plot([1, 2, 3, 4])
plt.xlabel('日本語化')
plt.ylabel('matplotlibの')
plt.show()

# seabornが日本語対応できているか確認する
import seaborn as sns
sns.set(style="whitegrid")

# Load the example Titanic dataset
titanic = sns.load_dataset("titanic")

# Draw a nested barplot to show survival for class and sex
g = sns.catplot(x="class", y="survived", hue="sex", data=titanic,
                height=6, kind="bar", palette="muted")
g.despine(left=True)
g.set_ylabels("seabornの日本語化")

ラベルに日本語が使われているので、ラベルが文字化け(豆腐)してなければ成功です。

matplotlibとseabornの日本語対応が確認出来たら、pandas-profilingでも日本語対応されているはずです。

終わりに

EDAを行う前にとりあえずpandas-profilingするのが定番になりそうな気がしています。

この記事を読んでみて、とても便利だなーと思ったので共有したいと思います。

モチベーション

以下のコードを想定しましょう。

num1 = 30
num2 = 40

print(num1)
print(num2)

出力は

30
40

となります。さて、上の出力を見て30と40がnum1とnum2のどちらに対応しているかすぐに判別がつくでしょうか？
出力変数が少ないうちはこれでも問題はなさそうですが5個とかになると少し面倒くさくなりますね。

とりあえずの解決策として次のようにすれば良いでしょう。

num1 = 30
num2 = 40

print('num1', num1)
print('num2', num2)

出力は

num1 30
num2 40

確かに変数名と変数値の対応は元のコードよりは分かりやすくなりました。しかし、問題はコードの記述が面倒という点です。

Icecreamとは？

上のような問題への解決策としてpythonライブラリのIcecreamがあります。
インストールはpipで行えます。

pip install icecream

from icecream import ic

num1 = 30
num2 = 40

ic(num1)
ic(num2)

出力は

ic| num1: 30
ic| num2: 40

非常にシンプルでいいですね！
Icecreamの便利性はこれだけではありません。関数や条件分岐など使えます。
興味ある人はこちらから調べてみてください。

楽天証券の実現損益データをグラフにしてみました。

実現損益とは、保有する株式などを売却したり決済したりして確定する損益のことを意味します。
それまでは含み損益と呼び、現金化されることでその損益を確定することになります。

つまり、実現損益を確認することで、自分が今まで得をしたのか、損をしたのかがわかります。
過去の成績を振り返り、同じミスをしないよう、今後の投資に活かしていきましょう。

※元記事→楽天証券の実現損益データをPythonでグラフ化【学習メモ】

csvファイルの取得

まずは楽天証券から実現損益のデータ（csvファイル）を取得します。
楽天証券にログインし、口座管理→損益・税金履歴→実現損益を開きます。

以下のような画面が表示されたら、画面右下の「CSV形式で保存」をクリックして保存します。

以下のように、楽天証券から実現損益のcsvファイルを取得することができました。（見にくいですが、99%損してます。）

※csvファイルの注意点

楽天証券から持ってくるcsvファイルですが、実は「銘柄コード」以外のすべてのデータ型が「object型」となっています。（dtypesでデータ型を確認）

データ型に関しては、以下の記事をご覧ください。

Pandasでデータ型を確認するdtype/dtypesと型変換を行うastypePandasでデータ型を確認するdtype/dtypesと型変換を行うastype | 侍エンジニア塾ブログ（Samurai Blog） - プログラミング入門者向けサイト

このままだと計算ができないので、「object型」を「int型」に変換する必要があります。
Pythonを使って型を変換しましょう。

と言いたいところでしたが、私にはどうしてもできませんでした。
ということで、今回はスプレッドシートで直接「int型」に変換するという苦肉の策をとりました。

スプレッドシートでファイルを開き、表示形式を「自動」にします。

これで「object型」を「int型」に変換することができました。

今後の課題ですね。

「こうすればできるよ！」という心優しい方がおられましたら、教えていただけたら幸いです。

作成の手順

それではcsvファイルを読み込み、実現損益をグラフ化していきましょう。

完成したコードはこちら。

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

import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_csv('drive/My Drive/csv/rakuten-int.csv', encoding='utf-8')

df.head()

df.dtypes

df['実現損益[円]'].cumsum().plot(figsize=(8,6),fontsize=18)
plt.xlabel("Number of transactions",size="12")
plt.ylabel("Realized profit/loss",size="12")
plt.grid(True)
plt.show()

まず初めに、csvファイルをGoogle Colabで読み込めるようにします。

しかし、Google Colabはインターネット上で実行するため、csvファイルの読み込みには工夫が必要です。

以下のコードを実行

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

「count」ディレクトリ以下の「drive」というフォルダに、Google Driveをマウントするという意味です。
Google Driveのファイルはここで指定したフォルダ「drive」以下の、「drive/My Drive」フォルダに格納されます。
このコードを実行すると認証URLが出てくるので、以下の記事の通り、アップデートまで進めてください。

【Google Colab】ファイルを読み込む方法【3分で出来る】

アップデートまで終わったら、csvファイルを読み込むために、pandasとmatplotlibをインポートします。

pandasとmatplotlibをインポート

import pandas as pd
import matplotlib.pyplot as plt

dfにcsvファイルを代入

df = pd.read_csv('drive/My Drive/csv/rakuten-int.csv', encoding='utf-8
')

read_csv：csvファイルを読み込む際に使用
'drive/My Drive/csv/rakuten-int.csv'：csvファイルまでの相対パス
encoding=’utf-8’：文字コード

楽天証券から持ってきたcsvファイルは文字コードに問題？があるため、エラーが発生します。
相対パスの後に必ず「encoding=’utf-8’」を記入しましょう。
詳しくは以下の記事をご確認ください。

Pythonにおけるunicode decode errorに関する回避方法を現役エンジニアが解説【初心者向け】

データが取得できたか確認

df.head()

データを取得することができました。

実現損益のデータをグラフ化

df['実現損益[円]'].cumsum().plot(figsize=(8,6),fontsize=18)
plt.xlabel("Number of transactions",size="12")
plt.ylabel("Realized profit/loss",size="12")
plt.grid(True)
plt.show()

df['実現損益[円]']：dfの中の['実現損益[円]']
cumsum()：累積和（数列に対して、最初の項からn番目の項まで足したもの）
例えば、日々の売上データからその日までの累計の売上データを出すような計算です。

x軸のラベルをプロット

plt.xlabel("Number of transactions",size="12")

Number of transactons：x軸の名前（日本語で「取引回数」）
size：フォントサイズ

y軸のラベルをプロット

plt.ylabel("Realized profit/loss",size="12")

Realized profit/loss：y軸の名前（日本語で「実現損益」）

軸のラベルが英語表記なのですが、日本語表記にすると文字化けします。
本当は日本語表記にしたかったのですが、長くなりそうなのでまた別の機会にします。

plt.grid(True)：グリッド（格子）をグラフに表示する
plt.show()：作ったグラフを画面に表示する

ここまでの実行結果がこちらです。

実現損益のデータをグラフにすることができました。

見事なWボトムからのジリ貧です。
グラフで可視化すると非常にわかりやすいですね。

一時はマイナス20万円近くまで行きましたが、急激に取り戻してプラス15万円近くまで回復し、調子に乗って結局マイナスで終わりました。

ちなみに、最終の実現損益は−5,091円です。（2020年4月30日でいったん終了）

１年間かけて資産をおよそ5,000円失いました。
勉強代と思えば安いものですね。

今後の課題

• Pythonを使って「object型」を「int型」へ変換する
• 軸のラベルを日本語表記にする
• 実現損益を±０以上にする

もっとPythonのスキルをつけて、詳細な分析ができるようにしたいですね。
そして、恥ずかしくない取引データがお見せできるように、資産運用も引き続き頑張ります。

はじめに

Lightgbmで構築したモデルに対してscikit-learnのplot_partial_dependenceを実行したときにNotFittedErrorが表示され、うまく実行できませんでした。
その対処についての備忘録です。

環境

環境は以下の通り。

$sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G14042

Jupyterlab (Version 0.35.4) 上で作業していたので、python kernelのバージョンも記載しておきます。

Python 3.7.3 (default, Mar 27 2019, 16:54:48) 
IPython 7.4.0 -- An enhanced Interactive Python. Type '?' for help.

やったこと

モデル構築

予測するためのデータとモデルを用意します。
データはscikit-learnで用意されているボストンデータセットを使用しました。

import pandas as pd
import sklearn.datasets as skd

data = skd.load_boston()

df_X = pd.DataFrame(data.data, columns=data.feature_names)
df_y = pd.DataFrame(data.target, columns=['y'])

以下の通り、506行13列のデータで全列がnon-nullのfloat型なのでこのままモデルを作ります。

df_X.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 506 entries, 0 to 505
Data columns (total 13 columns):
CRIM       506 non-null float64
ZN         506 non-null float64
INDUS      506 non-null float64
CHAS       506 non-null float64
NOX        506 non-null float64
RM         506 non-null float64
AGE        506 non-null float64
DIS        506 non-null float64
RAD        506 non-null float64
TAX        506 non-null float64
PTRATIO    506 non-null float64
B          506 non-null float64
LSTAT      506 non-null float64
dtypes: float64(13)
memory usage: 51.5 KB

LightGBMで、ハイパーパラメタはでほぼデフォルトのままでモデルを構築しておきます。

Training APIの場合

Training APIで用意されているtrainの関数でモデル構築した場合は、Boosterが返されます。
結論からいうと、Boosterをそのままscikit-learnのplot_partial_dependenceに渡すとエラーとなります。

どのようなエラーになるのか見ていきたいと思います。
まずはモデルを学習させます。

import lightgbm as lgb
from sklearn.model_selection import train_test_split

df_X_train, df_X_test, df_y_train, df_y_test = train_test_split(df_X, df_y, test_size=0.2, random_state=4)

lgb_train = lgb.Dataset(df_X_train, df_y_train)
lgb_eval = lgb.Dataset(df_X_test, df_y_test)

# Booster

params = {
    'seed':4,
    'objective': 'regression',
    'metric':'rmse'}

lgbm_booster = lgb.train(params,
                lgb_train,
                valid_sets=lgb_eval,
                num_boost_round=200,
                early_stopping_rounds=20,
                verbose_eval=50)

（省略）

Training until validation scores don't improve for 20 rounds
[50]    valid_0's rmse: 3.58803
[100]   valid_0's rmse: 3.39545
[150]   valid_0's rmse: 3.31867
[200]   valid_0's rmse: 3.28222
Did not meet early stopping. Best iteration is:
[192]   valid_0's rmse: 3.27283

これで学習させたモデルlgbm_boosterをplot_partial_dependenceに渡します。

from sklearn.inspection import plot_partial_dependence
plot_partial_dependence(lgbm_booster, df_X_train, ['CRIM', 'ZN']) 

---------------------------------------------------------------------------

ValueError                                Traceback (most recent call last)

<ipython-input-7-9a44e58d2f50> in <module>
      1 from sklearn.inspection import plot_partial_dependence
----> 2 plot_partial_dependence(lgbm, df_X_train, ['CRIM', 'ZN'])

（中略）

/anaconda3/lib/python3.7/site-packages/sklearn/inspection/_partial_dependence.py in partial_dependence(estimator, X, features, response_method, percentiles, grid_resolution, method)
    305     if not (is_classifier(estimator) or is_regressor(estimator)):
    306         raise ValueError(
--> 307             "'estimator' must be a fitted regressor or classifier."
    308         )
    309 


ValueError: 'estimator' must be a fitted regressor or classifier.

というわけでValueErrorが出ることが確認できました。

Scikit-learn APIの場合

それではScikit-learn APIでモデルを学習させます。
回帰問題のため、LGBMRegressorを使います。

lgbm_sk = lgb.LGBMRegressor(objective=params['objective'],
                           random_state=params['seed'],
                           metric=params['metric'])
lgbm_sk.fit(df_X_train, df_y_train)

LGBMRegressor(boosting_type='gbdt', class_weight=None, colsample_bytree=1.0,
              importance_type='split', learning_rate=0.1, max_depth=-1,
              metric='rmse', min_child_samples=20, min_child_weight=0.001,
              min_split_gain=0.0, n_estimators=100, n_jobs=-1, num_leaves=31,
              objective='regression', random_state=4, reg_alpha=0.0,
              reg_lambda=0.0, silent=True, subsample=1.0,
              subsample_for_bin=200000, subsample_freq=0)

これで学習させたモデルlgbm_skをplot_partial_dependenceに渡します。

plot_partial_dependence(lgbm_sk, df_X_train, ['CRIM', 'ZN']) 

---------------------------------------------------------------------------

NotFittedError                            Traceback (most recent call last)

<ipython-input-9-d0724528e406> in <module>
----> 1 plot_partial_dependence(lgbm_sk, df_X_train, ['CRIM', 'ZN'])

（中略）

/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py in check_is_fitted(estimator, attributes, msg, all_or_any)
    965 
    966     if not attrs:
--> 967         raise NotFittedError(msg % {'name': type(estimator).__name__})
    968 
    969 


NotFittedError: This LGBMRegressor instance is not fitted yet. Call 'fit' with appropriate arguments before using this estimator.

このようにNotFittedErrorが発生しました。
ここで使用しているscikit-learnとlightgbmのバージョンはそれぞれ0.22.2.post1と2.3.0です。

from sklearn import __version__ as sk_ver
print(lgb.__version__)
print(sk_ver)

2.3.0
0.22.2.post1

ここで発生したNotFittedErrorについて調べていると、GithubのIssueでこのエラーについてレポートされていました。
どうもscikit-learnのcheck_is_fittedという関数を使ってモデルがフィット済みか確認しているのですが、check_is_fittedで確認するattributeがlightgbmのfitで設定されていないために起こっているようです。
バージョンアップすることで解決できそうなので、lightgbmをバージョンアップします。
2021年1月14日の段階で最新バージョンは3.1.1だったため、このバージョンを指定しました。

conda install -c conda-forge lightgbm=3.1.1

バージョンを確認します。

from sklearn import __version__ as sk_ver
print(lgb.__version__)
print(sk_ver)

3.1.1
0.23.1

改めてモデルを学習させます。

lgbm_sk = lgb.LGBMRegressor(objective=params['objective'],
                           random_state=params['seed'],
                           metric=params['metric'])
lgbm_sk.fit(df_X_train, df_y_train)

LGBMRegressor(metric='rmse', objective='regression', random_state=4)

```python
plot_partial_dependence(lgbm_sk, df_X_train, ['CRIM', 'ZN']) 

今回はエラーなく実行でき、下図のような結果を得られました。

ついでにscikit-learnの他にpdpboxを使ってpartial dependenceを出してみます。

from pdpbox import pdp, get_dataset, info_plots
fig = plt.figure(figsize=(14,5))

pdp_goals = pdp.pdp_isolate(model=lgbm_sk, dataset=df_X_train, model_features=df_X_train.columns, feature='LSTAT')
# plot it
pdp.pdp_plot(pdp_goals,'LSTAT')
plt.show()

下図のような出力を得られました。scikit-learnよりも見た目はいいです。

まとめ

Lightgbmをscikit-learnのplot_partial_dependenceに渡したときに発生したNotFittedErrorはLightgbmをバージョンアップすることで解決しました。

筆者が以前投稿した記事 ( Matplotlib+ObsPyで震源球(beachball)を書いてみた ) の最後で ObsPy と地図ライブラリの Cartopy との連携が未確認であると書きましたが、地図上での描画が可能であるということが判明したので、これについて記述しようと思います。

※表示している地震のデータは防災科学技術研究所が公開している地震のメカニズム ( 防災科研：メカニズム解の検索 ) からランダムに選びダウンロードしたものです。本ページ中の画像は個人がデモンストレーションとして描画したもので、その結果に学術的な意味を求めていません。

※本ページのスクリプトの利用は自己責任でお願いします。
描画例

1. 目次

1. 目次
2. このページの目的
3. 使用環境
4. Cartopy のインストール方法
5. デモンストレーション
6. おわりに
7. 参考文献

2. このページの目的

下半球投影の震源球を、地図の上に描画すること。

3. 使用環境

• Windows 10
• conda 4.9.2
• python 3.8.5
• matplotlib 3.3.2
• obspy 1.2.2
• cartopy 0.18.0

cartopy はpipでインストールしました。

4. Cartopy のインストール方法

公式の記述によりますと Anaconda の環境を前提として、conda コマンドconda install -c conda-forge cartopy で入れるのが手っ取り早いそうです。また、required dependencies が満たされている状態ではpip install cartopy でインストールできるそうです。おそらくですが、Anaconda を使用している場合には pip でインストールしてもこの required dependency は自動で解決してくれるはずです。

ここからは蛇足ですが、公式インストーラーやストア経由で windows 10 に「直に」pythonをインストールした場合は、 Cartopy を導入するのは結構大変かもしれません。Anaconda に切り替える前に私が試したとき、 cartopy を pip で入れようとしたら proj のバージョンが低いと言われ、 proj を pip で入れようとしたら必要なバージョン ( proj4 ) がなくてダメで、proj4 をインストールしようとしたら OSGeo4w を入れる羽目になり、OSGeo4w をインストールしたはいいけどどうやって使うの？という状況に陥り、断念しました。初心者の私にもわかる記事があれば教えてください...

ちなみにUbuntuではAnacondaじゃなくても何の問題もなく入りました。Ubuntu は良い子です。

5. デモンストレーション

というわけで適当なスクリプトを実行してみます。筆者の前作で ObsPy は Matplotlib の Collection としてビーチボールを返してくれることが分かったので、地図上に描けることが分かっていれば意外とあっさりと書けます。

使用しているデータの冒頭部分は以下の通り。

temp.csv

longitude   latitude    depth   strike1 dip1    rake1   strike2 dip2    rake2   mantissa    exponent    type
136.3025    36.2495 11  63  59  94  236 31  84  2.66    21  1
135.5475    35.2028 20  334 85  11  243 79  175 9.92    21  2
136.72      35.1692 14  3   70  105 145 25  55  3.43    21  1
134.8153    34.4015 11  293 87  28  201 62  177 4.23    21  2

これを以下のスクリプトと同じディレクトリに入れ、実行します。Cartopy での地図の描画に関しては、公式のスクリプト例およびこの方のブログ (リンクフリーということでしたので掲載させていただきます...)を参考にさせていただきました。地図の描画を工夫したい人はそちらをご覧ください。

cartopy_obspy.py

import numpy as np
import pandas as pd
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
import matplotlib.ticker as mticker
from obspy.imaging.beachball import beach

# 地図とグリッドの描画
fig = plt.Figure(figsize=(5,4))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree()) # platecarree (正距円筒図法) で表示
ax.set_extent((134, 137, 33.5, 36.5)) # 領域の限定
grids = ax.gridlines(draw_labels=True) # ラベルありで緯線経線を描画
grids.xlocator = mticker.FixedLocator(np.arange(120,150,0.5)) # 0.5度ずつ
grids.ylocator = mticker.FixedLocator(np.arange(20,50,0.5)) 
ax.coastlines(resolution='10m') # 海岸線を分解能 10 m で描画
ax.add_feature(cfeature.LAND, color="lightgrey") # 陸地を灰色に

mfile = pd.read_csv("temp.csv") # データのファイル
mech = mfile[["strike1", "dip1", "rake1"]] # 発震機構解
loc = mfile[["longitude", "latitude"]] # 位置
typ = mfile[["type"]] # 正断層型か、など。0 = 正断層, 1 = 逆断層, 2 = 横ずれ断層
pallet = ["blue", "red", "lime"] # 断層タイプにより色分けする

for i in range(len(loc)) :
    beach1 = beach(mech.iloc[i,:], xy=loc.iloc[i,:], width=0.2,linewidth=0.4,facecolor=pallet[int(typ.iloc[i,0])])
    ax.add_collection(beach1)

fig.savefig("sample.png")
plt.close(fig)

結果は冒頭と同じ。

6. おわりに

やってみたら簡単でした。ただ現状これでは地図上でのプロットのみが可能であり、鉛直断面に震源球を描画する方法は現時点でも模索中です (たぶんやらない)。そういうのもできるという点でいえば GMT の方が圧倒的に自由度が高い気がします。 Python でやらないといけないという縛りがある状況以外では、このスクリプトに利用価値はないかもしれません。

7. 参考文献

• obspy.imaging - Plotting routines for ObsPy — ObsPy Documentation (1.2.0)
  https://docs.obspy.org/packages/obspy.imaging.html
• 19. Beachball Plot — ObsPy Documentation (1.2.0)
  https://docs.obspy.org/tutorial/code_snippets/beachball_plot.html
• obspy.imaging.beachball.beach — ObsPy Documentation (1.2.0)
  https://docs.obspy.org/packages/autoge/obspy.imaging.beachball.beach.html#obspy.imaging.beachball.beach
• 防災科学技術研究所 広帯域地震観測網 (F-net) メカニズム解の検索
  https://www.fnet.bosai.go.jp/event/search.php?LANG=ja
• Introduction — cartopy 0.18.0 documentation
  https://scitools.org.uk/cartopy/docs/latest/index.html
• Using cartopy with matplotlib — cartopy 0.18.0 documentation
  https://scitools.org.uk/cartopy/docs/latest/matplotlib/intro.html
• Python　－　（ベースマップ）　Cartopyの導入 | ふシゼン
  https://www.mitsumatado.com/zen/cartopy/

目次

0.パネル操作を行うサンプルコード
1.Add-onの実行手順
2.サンプルコードの解説

0.パネル操作を行うサンプルコード

↓今回のパネル操作Add-onでできること

↓サンプルコード

#アドオンの概要
#アドオン作成に必須
bl_info = {
    "name" : "Monkey grid",
    "author" : "Yohhei Suzuki",
    "version" : (1,0),
    "blender" : (2,83,0),
    "category" : "Mesh",
    "location" : "Operator Search",
    "description" : "More Monkeys",
    "warning" : "",
    "doc_url" : "",
    "tracker_url" : "",
}


import bpy

#猿を整列させるためのクラス
#前回の記事から変更なし
class MESH_OT_monkey_grid(bpy.types.Operator):
    bl_idname = "mesh.monkey_grid"
    bl_label = "Monkey Grid"
    bl_options = {'REGISTER','UNDO'}

    count_x: bpy.props.IntProperty(
        name = "x",
        description = "Number of Monkeys in the X-direction",
        default = 3,
        min = 1,soft_max = 10,
    )
    count_y: bpy.props.IntProperty(
        name = "y",
        description = "Number of Monkeys in the Y-direction",
        default = 2,
        min = 1,soft_max = 10,
    )
    size: bpy.props.FloatProperty(
        name = "Size",
        description = "Size of each Monkey",
        default = 0.5,
        min = 0.1,soft_max = 1,
    )

    @classmethod
    def poll(cls,context):
        return context.area.type == 'VIEW_3D'

    def execute(self, context):
        for i in range(self.count_x * self.count_y):
            x = i %  self.count_x#only amari
            y = i // self.count_x #only shou
            bpy.ops.mesh.primitive_monkey_add(
                size = self.size,
                location = (x,y,1)
            )

        return {'FINISHED'}


#今回追加したクラス
#パネルでの操作を可能にする
class VIEW3D_PT_monkey_grid(bpy.types.Panel):
    #クラスの情報
    bl_space_type = 'VIEW_3D'
    bl_region_type = 'UI'
    bl_category = "Monkeys"
    bl_label = "Grid"

    #パネルの表示と、パネルが押された時の変数の値を決めている
    def draw(self,context):
        #align = Trueでパネル同士の隙間を埋める
        col = self.layout.column(align = True)
        #パネルの表示と
        props = col.operator('mesh.monkey_grid',
            text = 'Default Grid',
            icon = 'MONKEY')
        #パネルが押された時の変数の値を決定する
        props.count_x = 3
        props.count_y = 2
        props.size = 0.5

        #大きいmonkeyを多く並べる
        props = col.operator('mesh.monkey_grid',
            text = 'Big Grid',
            icon = 'MONKEY')
        props.count_x = 10
        props.count_y = 10
        props.size = 0.8

        #小さいmonkeyを1つ並べる
        props = col.operator('mesh.monkey_grid',
            text = 'Small Grid',
            icon = 'MONKEY')
        props.count_x = 1
        props.count_y = 1
        props.size = 0.2

        #cyclesでのpreview_samplesの値をパネルから操作可能にする
        col = self.layout.column(align = True)
        col.prop(context.scene.cycles, 'preview_samples')

        #アクティブオブジェクトが存在しない場合の表示
        if context.active_object is None:
            col.label(text = '-no active object-')
        #アクティブオブジェクトが存在する場合はそれを隠す
        else:
            col.prop(context.active_object, 'hide_viewport')      

#add→meshからアドオンを選択できるようにする
def mesh_add_menu_draw(self,context):
    col = self.layout.column(align = True)
    props = col.operator('mesh.monkey_grid')
    props.count_x = 3
    props.count_y = 2
    props.size = 0.5    

#Blenderに実行するクラスを登録する
def register():
    bpy.utils.register_class(MESH_OT_monkey_grid)
    bpy.utils.register_class(VIEW3D_PT_monkey_grid)
    bpy.types.VIEW3D_MT_mesh_add.append(mesh_add_menu_draw)

def unregister():
    bpy.utils.unregister_class(MESH_OT_monkey_grid)
    bpy.utils.unregister_class(VIEW3D_PT_monkey_grid)
    bpy.types.VIEW3D_MT_mesh_add.remove(mesh_add_menu_draw)

1.Add-onの実行手順

①コードを保存する

●~~.py(ex:monkey_grid_panel.py)と名前をつける→Text→Save as→任意のフォルダに保存

---

②アドオンを使用可能にする

●Edit→Preferences→Add-ons→Install→保存した~~.py(ex:monkey_grid_panel.py)を選択→Install_Add-on→アドオンにチェックを入れる

---

③アドオンを実行する

●3D viewportで作成したパネルから操作する。もしくは、Add→mesh→アドオン名(ex:Monkey grid)を選択して実行する。

2.サンプルコードの解説

bl_info = {
    #---アドオンの概要---
    #---アドオン作成に必須---
}


import bpy

#猿を整列させるためのクラス
#前回の記事から変更なし
class MESH_OT_monkey_grid(bpy.types.Operator):
    #---monkeyを整列させる処理---



#今回追加したクラス
#パネルでの操作を可能にする
class VIEW3D_PT_monkey_grid(bpy.types.Panel):
    #クラスの情報
    bl_space_type = 'VIEW_3D' #実行する場所
    bl_region_type = 'UI' #クラスの種類
    bl_category = "Monkeys" 
    bl_label = "Grid"


    #押されたパネルに応じた処理を実行する関数
    def draw(self,context):
        #align = Trueでパネル同士の隙間を埋める
        col = self.layout.column(align = True)
        #パネルの表示を決定する
        #パネルが押されたときにMESH_OT_monkey_gridクラスの処理を行う
        #MESH_OT_monkey_gridクラス→bl_idname = "mesh.monkey_grid"
        props = col.operator('mesh.monkey_grid',
            text = 'Default Grid',
            icon = 'MONKEY')
        #パネルが押された時の変数の値を決定する
        props.count_x = 3
        props.count_y = 2
        props.size = 0.5

        #大きいmonkeyを多く並べる
        props = col.operator('mesh.monkey_grid',
            text = 'Big Grid',
            icon = 'MONKEY')
        props.count_x = 10
        props.count_y = 10
        props.size = 0.8

        #小さいmonkeyを1つ並べる
        props = col.operator('mesh.monkey_grid',
            text = 'Small Grid',
            icon = 'MONKEY')
        props.count_x = 1
        props.count_y = 1
        props.size = 0.2

        #cyclesでのpreview_samplesの値をパネルから操作可能にする
        col = self.layout.column(align = True)
        col.prop(context.scene.cycles, 'preview_samples')

        #アクティブオブジェクトが存在しない場合の表示
        if context.active_object is None:
            col.label(text = '-no active object-')
        #アクティブオブジェクトが存在する場合はそれを隠す
        else:
            col.prop(context.active_object, 'hide_viewport')      

#add→meshからアドオンを選択できるようにする
def mesh_add_menu_draw(self,context):
    col = self.layout.column(align = True)
    props = col.operator('mesh.monkey_grid')
    #選択したときの変数の値を指定する
    props.count_x = 3
    props.count_y = 2
    props.size = 0.5    

#Blenderに実行するクラスを登録する
def register():
    #monkeyを整列させるクラス(の登録)
    bpy.utils.register_class(MESH_OT_monkey_grid)
    #パネルの作成とパネルごとの処理を実行するクラス(の登録)
    bpy.utils.register_class(VIEW3D_PT_monkey_grid)
    bpy.types.VIEW3D_MT_mesh_add.append(mesh_add_menu_draw)

def unregister():
    bpy.utils.unregister_class(MESH_OT_monkey_grid)
    bpy.utils.unregister_class(VIEW3D_PT_monkey_grid)
    bpy.types.VIEW3D_MT_mesh_add.remove(mesh_add_menu_draw)

何日で達成できるでしょうか

最適化問題の勉強をしておりまして、思いついてしまった問題。
日本全国を最短距離で巡るのってどれぐらい時間がかかるのだろうか？
という疑問に答えるべく、pythonでちょっと計算してみました。

最適化問題をpythonでハンズオンすることが目的なので、問題設定の甘さはご容赦ください。
また、こちらの記事でも書きましたが、結局本命経路以外で円環経路を作ってしまう問題は解決できませんでした。
最適解ではなく近似解です、この辺りの解決法に詳しい方いらっしゃいましたらコメントお願いします！！

問題設定

算数の問題で良く出る（本当かは知らん）みんな大好きたかしくんを今回は使います。

たかしくんスペック

・たかしくんは強いので海を渡れます。
・たかしくんは強いので山があろうが谷があろうが一定の速度で走れます。
・たかしくんは強いので24時間走れます。

設定

・たかしくんは知的好奇心が旺盛なので日本全国を時速10kmで巡ろうと考えました。
・「日本全国」を今回は「全都道府県庁県庁所在地」と定義します
・たかしくんはせっかちなので移動は全て直線で行うものとします。

やってみた

都道府県庁所在地データの取得

都道府県庁所在地の経度と緯度はこちらの記事から引用しました。
Excelをダウンロード後、所在地名・経度・緯度のデータだけに成形します。

import pandas as pd
import numpy as np
df = pd.read_excel("city_list.xlsx")
display(df)

データはこんな感じ。latiが緯度(latitude)、longが経度(longitude)になります。

各所在地間の距離計算

47か所の所在地のお互いの距離計算をします。
本来であれば、例えば「沖縄-北海道なんて移動しないだろう」と考え組み合わせを絞るのが良い気もします。
今回はその選抜作業が面倒くさかったのと、何かミラクルが起きるかもという淡い期待からすべての組み合わせを網羅して計算しました。

経度・緯度を利用した距離計算はこちらの記事を流用しました。
二点の経度と緯度を引数(それぞれリスト型であることに注意！)として、その間の距離をkmで返す関数です。

def cal_rho(lon_a,lat_a,lon_b,lat_b):
    ra=6378.140  # equatorial radius (km)
    rb=6356.755  # polar radius (km)
    F=(ra-rb)/ra # flattening of the earth
    rad_lat_a=np.radians(lat_a)
    rad_lon_a=np.radians(lon_a)
    rad_lat_b=np.radians(lat_b)
    rad_lon_b=np.radians(lon_b)
    pa=np.arctan(rb/ra*np.tan(rad_lat_a))
    pb=np.arctan(rb/ra*np.tan(rad_lat_b))
    xx=np.arccos(np.sin(pa)*np.sin(pb)+np.cos(pa)*np.cos(pb)*np.cos(rad_lon_a-rad_lon_b))
    c1=(np.sin(xx)-xx)*(np.sin(pa)+np.sin(pb))**2/np.cos(xx/2)**2
    c2=(np.sin(xx)+xx)*(np.sin(pa)-np.sin(pb))**2/np.sin(xx/2)**2
    dr=F/8*(c1-c2)
    rho=ra*(xx+dr)
    return rho

あとはこの関数に全組み合わせをぶち込んでいきます。

#引数となるリストの入れ物を用意
loc_A = []
lon_A = []
lat_A = []

loc_B = []
lon_B = []
lat_B = []

#ひたすらぶち込む
for i in range(47):
    for j in range(47):
        loc_A.append(i)
        lon_A.append(df["long"][i])
        lat_A.append(df["lati"][i])

        loc_B.append(j)
        lon_B.append(df["long"][j])
        lat_B.append(df["lati"][j])

#計算！
rho=cal_rho(lon_A,lat_A,lon_B,lat_B)

#結果をデータフレームに
combi_df = pd.DataFrame([loc_A,loc_B,lon_A,lat_A,lon_B,lat_B,rho]).T
combi_df.columns = ["loc_A","loc_B","long_A","lati_A","long_B","lati_B","Dist/km"]
combi_df = combi_df.fillna(0)
combi_df.head(5)

結果、各地点間の距離を計算することができました。

最短経路の計算

ここについては前回の記事やおおもとの記事とほとんど同じ実装になっていますが、matplotlibによる描画の大きさだけ変えています(参考記事)
図の描画の関数がこちら。

import random, numpy as np, pandas as pd, networkx as nx, matplotlib.pyplot as plt
from itertools import chain, combinations
from pulp import *

def draw(g):
    plt.figure(dpi=150)
    rn = g.nodes() # ノードリスト
    pos = nx.spring_layout(g) # ノード位置
    """描画"""
    nx.draw_networkx_labels(g, pos=pos)
    nx.draw_networkx_nodes(g, node_color='w', pos=pos)
    nx.draw_networkx_edges(g, pos=pos)
    plt.show()

最適化計算のコードはこちら。
Networkxについてはこちらの記事が読みやすいです！
今回は取りうる値を0,1にした離散最適化問題として解きました(参考記事)。

#47地点を定義
n = 47 # ノード数
g = nx.random_graphs.fast_gnp_random_graph(n, 1, 8)

#47地点間の距離をインプット
for a,b in g.edges():
    col_idx = a * 47 + b
    g.edges[a, b]["dist"] = combi_df.loc[col_idx,"Dist/km"]

#いざ計算
source, sink = 0, 46 # 始点, 終点
r = list(enumerate(g.edges()))
m = LpProblem() # 数理モデル
x = [LpVariable('x%d'%k, lowBound=0, upBound=1,cat="Integer" ) for k, (i, j) in r] # 変数(路に入るかどうか)
m += lpSum(x[k] * g.edges[i,j]["dist"] for k, (i, j) in r) # 目的関数
for nd in g.nodes():
    m += lpSum(x[k] for k, (i, j) in r if i == nd) \
        + lpSum(x[k] for k, (i, j) in r if j == nd) == {source:1, sink:1}.get(nd, 2) # 制約
m.solve()

print([(i, j) for k, (i, j) in r if value(x[k]) == 1])
print(LpStatus[m.status])
print(value(m.objective))

一つ目のprintは経路の表示、二つ目は最適化計算がうまくいったか、三つめは最適解の表示です。
出力結果はこちら。

optimalとなっているので最適解が叩き出せたみたい！！
そして最短経路は4304km！！
なのでたかしくんは430時間、つまり約17日22時間で日本全国を巡れるらしい！

やったーーー！計算できた！！！
と思っていたのもつかの間でした。

円環経路問題

どんな経路を巡ったのだろうか？と思い、描画をしてみました。
最適化結果の確認用のコードはこちら。

G = nx.Graph()
for k,(i,j) in r:
    if value(x[k]) ==1:
        nx.add_path(G, [i, j])
draw(G)

描画結果がこちら。
前回の記事でも取り上げた円環生成問題が生じていました。
これは、「ある地点に着目した時に、経由点では2本の道、始点終点では1本の道が伸びている」という制約条件の甘さから生じる問題です。

こちらのサイトから白地図を拝借しプロットしてみました。
東北・関東・中国地方でぐるぐる回ってる場所がありますね。

ここからは経路を少しいじくって、一本の経路にしていこうと思います。
恣意的な操作が加わりますので最適解になるのかが分からなくなりますが、おおよそ近いだろうと踏んで考えてみます。

秋田-新潟間を切る

切りたい経路の長さをとてつもなく長くしてみます。
こんな感じです。

pair_list = [(4,14),(14,4)]
for a,b in pair_list:
    col_idx = a * 47 + b
    combi_df.loc[col_idx,"Dist/km"] = 1000000

これで最適化計算をしてみます。
結果は地図だけ示しますね。
東北地方の統合に成功しました！！やったぜ。
この場合でも4309kmぐらいなので、経路一つだけの変更ではそこまで影響がなさそうですね。

新潟-長野間を切る

関東地方の統合をどうするか悩みました。
そこで、新潟から長野に行かずに群馬へ行き、関東の円環に合流した後に長野へいく、というルートを考えてみました。
なので新潟-長野間の経路を切ってみました。

その結果がこちら。

新潟から栃木へ行き、うまく統合することができました。
距離はやや伸びてしまいましたが4326kmとなりました。

島根-岡山間を切る

最後、中国地方の統合をどうするか。
島根-岡山間の移動が見られますが、岡山-島根を移動するのであれば鳥取か広島を経由するのが良さそう。
なので島根-岡山間を切ってみます。

香川を巻き込みやがった・・・
どれだけそこでこじんまりしたいんだよ。

島根-香川間を切る

ということで島根-香川間を切ってみます。
その結果がこちら！
島根から高知へワープしているものの、ついに一本になりました！！
距離は4361kmに。

ちょっと島根-高知間のワープが気になったのでここもぶち切ってみました。
今度は高知から大分へワープして、山口は九州に行ってから拾いに行くようです。
距離も4361kmと同じぐらいになりました。
中国・四国周りは一筆書きするのが難しそうですね。

結論

走る距離：だいたい4361km
走る時間：だいたい18日4時間

さいごに

最適化問題面白い！
これからもいろいろ実験してみようと思います。
最後まで読んでいただきありがとうございます、LGTMしていただけると励みになります！

AOJ Knapsack Problem

個数制限なしのナップサック問題

価値が $v_i$ 重さが $w_i$ であるような $N$ 種類の品物と、容量が $W$ のナップザックがあります。

次の条件を満たすように、品物を選んでナップザックに入れます：

最適化問題に適用する場合、一般的に、以下の2つが適用する問題に成立していないといけない:

• 選んだ品物の価値の合計をできるだけ高くする。
• 選んだ品物の重さの総和は $W$ を超えない。
• 同じ種類の品物はいくつでも選ぶことができる。

価値の合計の最大値を求めてください。

問題の解法

動的計画法による外部設計は次のようになります。

$dp[i+1][j]$の定義：
$i$ 問目までの品物から重さの総和が$j$以下となるように選んだときの、価値の総和の最大値とする。

dp初期条件：

dp[0][j]=1

dp漸化式の定義：

dp[i+1][j]=max\{dp[i][j-k×w[i]]+k×v[i] \}

これだと三重ループになってしまうため、計算量は $O(nW^2)$ となり、タイムオーバーです。漸化式の式変形を行うことで $k$ のループをなくすことができ、計算量を $O(nW)$ とすることができます。この変形は情報処理における重複計算を省く目的です。理解の仕方としては、$dp[i+1][j]$ を求めるには $dp[i][j]$ と $dp[i+1][j-w[i]]+v[i]$ を見れば十分ということであり、逆にループを回すということでもあります。蟻本の解説が分かり易いです。

dp漸化式の変形：

dp[i+1][j]=max(dp[i][j], dp[i+1][j-w[i]]+v[i])

求める解：$dp[n][W]$

サンプルコード

n = int(input())
w = []
v = []
for i in range(n):
    w_, v_ = map(int, input().split())
    w.append(w_)
    v.append(v_)
W = int(input())

dp = [[0] * (W + 1) for i in range(n + 1)]

for i in range(n):
    for j in range(W + 1):
        if j < w[i]:
            dp[i + 1][j] = dp[i][j]
        else:
            dp[i + 1][j] = max(dp[i][j], dp[i + 1][j - w[i]] + v[i])

print(dp[n][W])

はじめに

最適化アルゴリズムの実装シリーズです。
まずは概要を見てください。

コードはgithubにあります。

くじらさんアルゴリズム

概要

くじらさんアルゴリズム(The Whale Optimization Algorithm; WOA)はクジラの捕食行動に着目して作られたアルゴリズムです。
クジラはバブルネットという群行動によって獲物を捕食します。
バブルネットは図のように円を描きながら獲物を追い詰めて捕食する行動のようです。

参考
・くじらさんアルゴリズム
・Whale Optimization Algorithm

アルゴリズム

獲物に近づく、獲物を探す、円を描くの３つの行動のどれかを選んで獲物を見つけます。

• アルゴリズムのフロー

• 用語の対応

問題	くじらさんアルゴリズム
問題への入力	クジラの位置≒獲物の位置
評価値	クジラの位置における評価値

• ハイパーパラメータに関して

変数名	意味	備考
a_decrease	aの減少値	最初はランダムに移動し、後半は獲物に近づく。aはその移行速度
logarithmic_spiral	対数螺旋の係数	大きいほど回るときに大きく動く

獲物に近づく

獲物に近づくと書いていますが実際に獲物がどこにいるかは分からないので、現状で最も評価が高いクジラを獲物と扱い、そこに近づきます。

まずは $\vec{A}$ と $\vec{C}$ を出します。

$$ \vec{A} = 2a\vec{r_1} - a$$

$$ \vec{C} = 2\vec{r_2}$$

$\vec{r_1}$ と $\vec{r_2}$ は座標と同じ次元の 0～1 の乱数です。
$a$ は初期値2から毎step減少していく変数です。($2 \geqq a \geqq 0$)

ここで $\vec{A}$ のノルム（$|\vec{A}|$）が1以下なら獲物に近づく行動をし、
違う場合は獲物を探す行動をします。

獲物に近づく場合の数式は以下です。

$$ \vec{D} = |\vec{C} \vec{X_{best}}(t) - \vec{X}(t)| $$
$$ \vec{X}(t+1) = \vec{X_{best}}(t) - \vec{A}\vec{D} $$

$\vec{X}$ が自分の座標を表し、$\vec{X_{best}}$ が最も評価の高いクジラです。
コードだと以下になります。

import numpy as np

a = 2  # 毎step減少する

r1 = np.random.rand(problem.size)  # 0-1の乱数を生成
r2 = np.random.rand(problem.size)  # 0-1の乱数を生成

A = (2.0 * np.multiply(a, r1)) - a
C = 2.0 * r2

i = 対象のクジラのindex

if np.linalg.norm(A) < 1:
    pos = np.asarray(whales[i].getArray())  # 対象クジラの座標
    best_pos = np.asarray(best_whale.getArray())  # 最良クジラの座標

    D = np.linalg.norm(np.multiply(C, best_pos) - pos)
    pos = best_pos - np.multiply(A, D)

獲物を探す

獲物に近づくで説明した $|A|$ が1以上の場合は獲物に近づきます。
近づき方は"獲物に近づく"場合と同じで、対象のクジラが最も評価の高いクジラからランダムなクジラに変わっただけです。

$$ \vec{D} = |\vec{C} \vec{X_p}(t) - \vec{X}(t)| $$
$$ \vec{X}(t+1) = \vec{X_p}(t) - \vec{A}\vec{D} $$

$\vec{X_p}$ はランダムなクジラの座標となります。

import numpy as np

i = 対象のクジラのindex

if np.linalg.norm(A) >= 1:
    pos = np.asarray(whales[i].getArray())  # 対象クジラの座標
    p_pos = np.asarray(whales[random.randint(0, len(whales)-1)].getArray())  # ランダムなクジラの座標

    D = np.linalg.norm(np.multiply(C, p_pos) - pos)
    pos = p_pos - np.multiply(A, D)

円を描く

円は獲物（評価値最大のクジラ）を基準に以下の式で描きます。

$$ \vec{D'} = |\vec{X_{best}}(t) - \vec{X}(t)| $$
$$ \vec{X}(t+1) = \vec{D'} e^{bL} \cos(2 \pi L) + \vec{X_{best}}(t) $$

$b$ は対数螺旋の係数(はlogarithmic_spiral)、$L$ は[-1, 1]のランダムな値です。

import numpy as np

i = 対象のクジラのindex

pos = np.asarray(whales[i].getArray())  # 対象クジラの座標
best_pos = np.asarray(best_whale.getArray())  # 最良クジラの座標

# 回る
D = np.linalg.norm(best_pos - pos)
L = np.random.uniform(-1, 1, problem.size)  # [-1,1]の乱数
_b = logarithmic_spiral
pos = np.multiply(np.multiply(D, np.exp(_b*r)), np.cos(2.0*np.pi*r)) + best_pos

スパイラルのグラフ

円を描くがどういう移動をしているかを見てみました。

グラフのx軸が-1から1の乱数($L$）で、y軸が移動距離です。
数式だと以下の部分だけをグラフ化しています。

$$ spiral = \vec{D'} e^{bL} \cos(2 \pi L) $$

• ユークリッド距離による変化

$b$(logarithmic_spiral) を1に固定して、$\vec{D'}$のみを変化させた場合のグラフです。
$\vec{D'}$ は自分と獲物(最良クジラ)とのユークリッド距離を表しているので、獲物との距離による移動の変化を表しているといえます。

離れているほど大きく移動する可能性が高い事がわかりますね。

• 対数螺旋係数による変化

今度は$b$(logarithmic_spiral)を変化させた場合です。
$D'$ は 1 に固定しています。

$D'$ の場合と同じように $b$ が大きくなれば変化が大きいようです。

コード全体

コード全体です。
コードを作成するにあたりここのコードを参考にしています。

import math
import random

import numpy as np

class WOA():
    def __init__(self, 
        whale_max,         # クジラの数
        a_decrease=0.001,  # 変数aの減少値
        logarithmic_spiral=1,  # 対数螺旋の係数
    ):
        self.whale_max = whale_max
        self.a_decrease = a_decrease
        self.logarithmic_spiral = logarithmic_spiral

    def init(self, problem):
        self.problem = problem

        self.best_whale = None
        self.whales = []
        for _ in range(self.whale_max):
            o = problem.create()
            self.whales.append(o)

            if self.best_whale is None or self.best_whale.getScore() < o.getScore():
                self.best_whale = o.copy()
        self._a = 2


    def step(self):

        for whale in self.whales:
            pos = np.asarray(whale.getArray())

            if random.random() < 0.5:
                r1 = np.random.rand(self.problem.size)  # 0-1の乱数
                r2 = np.random.rand(self.problem.size)  # 0-1の乱数

                A = (2.0 * np.multiply(self._a, r1)) - self._a
                C = 2.0 * r2

                if np.linalg.norm(A) < 1:
                    # 獲物に近づく
                    new_pos = self.best_whale.getArray()
                else:
                    # 獲物を探す
                    new_pos = self.whales[random.randint(0, len(self.whales)-1)].getArray()
                new_pos = np.asarray(new_pos)

                D = np.linalg.norm(np.multiply(C, new_pos) - pos)
                pos = new_pos - np.multiply(A, D)

            else:
                # 回る
                best_pos = np.asarray(self.best_whale.getArray())

                D = np.linalg.norm(best_pos - pos)
                L = np.random.uniform(-1, 1, self.problem.size)  # [-1,1]の乱数

                _b = self.logarithmic_spiral
                pos = np.multiply(np.multiply(D, np.exp(_b*L)), np.cos(2.0*np.pi*L)) + best_pos

            whale.setArray(pos)
            if self.best_whale.getScore() < whale.getScore():
                self.best_whale = whale.copy()

        self._a -= self.a_decrease
        if self._a < 0:
            self._a = 0

ハイパーパラメータ例

各問題に対して optuna でハイパーパラメータを最適化した結果です。
最適化の1回の試行は、探索時間を5秒間として結果を出しています。
これを100回実行し、最適なハイパーパラメータを optuna に探してもらいました。

問題	a_decrease	logarithmic_spiral	whale_max
EightQueen	0.020366541568838378	1.5723091675627572	45
function_Ackley	0.012846425355295324	1.817054209171554	49
function_Griewank	0.014055305155620233	1.5507180701447845	48
function_Michalewicz	0.019651164901908023	1.7671279692872293	46
function_Rastrigin	0.0173318428180831	0.7578390758302801	46
function_Schwefel	0.007491644624065206	1.6917050129680944	9
function_StyblinskiTang	0.024426401837372255	1.9474471085818506	50
function_XinSheYang	0.005722874915111857	0.4779624408790928	24
g2048	0.043666294831303784	1.09050483219953	37
LifeGame	0.0054058667884643585	0.06834119701477159	50
OneMax	0.014922301994610046	0.8650190784964947	27
TSP	0.0287871077043457	1.3855183848189636	45

実際の動きの可視化

1次元は6個体、2次元は20個体で50step実行した結果です。
赤い丸がそのstepでの最高スコアを持っている個体となります。

パラメータは以下で実行しました。

WOA(N, a_decrease=2/50, logarithmic_spiral=1)

function_Ackley

• 1次元

• 2次元

function_Rastrigin

• 1次元

• 2次元

function_Schwefel

• 1次元

• 2次元

function_StyblinskiTang

• 1次元

• 2次元

function_XinSheYang

• 1次元

• 2次元

あとがき

ぐるぐる渦巻きながらいい感じに集まっていきますね。
ただ、一度収束してしまうとそこから脱出することはないので a の値をどう制御するかが重要そうです。

np.argsortをnp.sortと同じ勢いで使ったら全く違うよね

下の画像がすごく参考になった

てっきりarrayの順番が返ってくると思ってた(上の写真なら[1,2,0])から、かなり混乱したけど勘違いだった。

参考
https://www.headboost.jp/numpy-argsort/

pythonでシステムトレード：移動平均線を描く

システムトレードの作戦

コロナで世の中大変ですが、株式は２１年１月現在コロナ金融緩和バブルの様相を示しています。今はこのトレンドにのって、イケイケですが「節分天井,彼岸底」という格言通り、イケイケは後少しだと思います。その後は次のトレンドが来るまではしばらくレンジ相場が来るような気がします。レンジ相場、トレンド相場に対応できる機械学習を使ったシステムトレードを考えようと思います。
pythonのライブラリーとしては、TA-Libというライブラリーがありますが、勉強のために、なるべく使わずに検討してみようと思います。

今回の作戦としては、
１．etfの売買を行う
２．ボリンジャーバンドを指標として使う
３．ボリンジャーバンドのバンドの下限で買い、上限で売り

いろいろな指標はありますが、単純なボリンジャーバンド作戦です。
ボリンジャーバンドでpythonの勉強をしたら、他に良い指標（ＭＡＣＤ，ＲＳＩ）がないか、指標を組み合わせたらどうなるかと順番に作っていけたら良いなと思います。

最初にボリンジャーバンドのシステムトレードを作りながら
１．で、どのetfが良いか？topix、日経平均、その他の指標でどれが勝率がよいか。
２．３．で、ボリンジャーバンドの何σが勝率がよいか考える

を確認して、システムトレードを作ろうと思います。

データの入手

まずターゲットを1306の野村アセットのＴＯＰＩＸ連動型上場投信で検討します。
スクレイピング技術は後回しです。クリックして、このサイトから集めました。
https://kabuoji3.com/
少し面倒ですが、クリックするだけで２００１年から毎年ののｃｓｖファイルが入手できます。
1306_2001.csv ～　1306_2021.csv

こういう日々のetfの値動きのファイルが入手できます。

データの結合

まずはデータの結合です。これは簡単にpandasで結合しました。
深く考えずに結合してますが、0行目に「1306 東証ETF TOPIX連動型上場投資信託（ETF）」が入っているので、data = pd.read_csv(file,encoding='Shift-Jis',header = 1)として、header=1で１行目の

をcolumn 名として読み込みました。
column 名が日本語なので、とりあえず適当な英語に直しました。

import pandas as pd
import matplotlib.pyplot as plt
import datetime

import matplotlib.dates as mdates

#read_csvですべてのデータを結合する。
#1306_2001から1306_2021

df = None

for i in range(2001,2021):

    file = '1306_' + str(i) + '.csv'
    print(file)

    data = pd.read_csv(file,encoding='Shift-Jis',header = 1)


    df = pd.concat( [df,data],ignore_index=True)

df.rename(columns={'日付': 'day', '始値': 'open_price',
                   '高値':'high_price', '安値':'low_price',
                   '終値': 'end_price', '出来高':'volume',
                   '終値調整値':'end_price_ad_value'}, inplace=True)

print(df.head())

とりあえずpandasのデータセットができました。

グラフを作ると、こんな感じです。今は１０年間で最高値です。イケイケですね。

次にボリンジャーバンドの描き方です。

ボリンジャーバンドの描き方

ボリンジャーバンドとは、をＭＯＮＥＸ証券のサイトから引用します。
https://info.monex.co.jp/technical-analysis/indicators/003.html
「移動平均線と標準偏差で構成されており、移動平均を表す線とその上下に値動きの幅を示す線を加えた指標で、「価格の大半がこの帯（バンド）の中に収まる」という統計学を応用したテクニカル指標のひとつです。」

つまり、株の価格から移動平均線と標準偏差を計算する必要があります。後、「統計学を応用したテクニカル指標」というところで、機械学習にも相性良いかも知れませんね。
最初に移動平均線を計算します。

単純移動平均線の計算

単純移動平均線は、一定期間の終値の平均価格になります。そこで、pandasのrolling関数を使います。例えばrolling関数で、rolling(5).mean()とすると、５つの数値の平均を返します。そこで、df["end_price"].rolling(5).mean().round(1)で、計算をして、５日移動平均線を描いてみます。

#'day'をdatetimeタイプに変更
df['day'] = pd.to_datetime(df['day'])
#'day'をindexにして置き換える。
df.set_index('day', inplace=True)
#５日間移動平均線を計算
df["5day_ave"]=df["end_price"].rolling(5).mean().round(1)

#2019年の下期７月から１２月を取り出す
df_temp = df['2019-07':'2019-12']

#グラフに描く
plt.plot(df_temp.index, df_temp["end_price"], label="daily")
plt.plot(df_temp.index, df_temp["5day_ave"],  label="sma5")
plt.xlabel("day")
plt.ylabel("price")
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%d-%b-%Y'))
plt.gca().xaxis.set_major_locator(mdates.MonthLocator(interval=2))
plt.gca().xaxis.set_minor_locator(mdates.MonthLocator(interval=1))
plt.gcf().autofmt_xdate()

plt.legend()
plt.show()

移動平均線が描けました。
次は、ボリンジャーバンドを描いていきます。
それは次の「投資quest:pyhtonでシステムトレードを作る(２)」（今下書き中です）で。

はじめに

強化学習の基本となるQ学習について復習も兼ねて書く。

理論

Q学習とは強化学習手法TD学習の一つで、Q値(状態行動価値)をエージェントが行動するたびに更新する手法である。時刻$t$での状態を$s_t$、行動を$a_t$、状態$s_t$のもとで行動$a_t$を起こすことによって得られる報酬を$r_t$とする。また、Q値$Q(s_t,a_t)$とはある状態$s_t$においてある行動$a_t$を取った時の価値のことである。価値の更新は以下のように行われる。
$$
Q(s_t,a_t) \leftarrow Q(s_t,a_t) + \alpha(r_{t+1}+\gamma \max_{a_{t+1}}Q(s_{t+1},a_{t+1})-Q(s_t,a_t))
$$
$\alpha$は学習の大きさを制御するパラメーターで0~1の値を取る。$\gamma$は将来の価値をどれほど考慮するかのパラメータで0~1をとる。$\max_{a_{t+1}}Q(s_{t+1},a_{t+1})$は状態$s_{t+1}$において得られる最大の価値である。この式は$Q(s_t,a_t)$、$\alpha(r_{t+1}-Q(s_t,a_t))$、$\alpha\gamma\max_{a_{t+1}}Q(s_{s_{t+1}}{a_{t+1}})$に分けて考えると理解しやすい。第一式は前回の価値であり、基準と考える。つまり第二式、第三式によって更新の程度が決まる。第二式は前回の価値と報酬の差をなくすように更新している。単に差を加算するだけではその時偶然得た報酬に過剰に依存してしまうため$\alpha$によって制御している。第三式は将来の価値について見積もっている。次に得られる最大の価値を将来の価値として考えており、$\alpha\gamma$によって制御されている。$\alpha$は第二式と同様の制御、$\gamma$は見積もりの不確かさによる割引率である。

実戦(具体例)

具体例として2腕バンディット問題について考える。0~2000の数字がランダムに出てくる山Aと0~1000の数字がランダムに出てくる山Bがあるとする。どちらか一方の山から数字を取得とするとき、より大きな数字(報酬)を得ることができる可能性が高い山を選択できるようなモデルを作成する。更新回数$t$(更新回数$t$は理論での時刻$t$とは異なることに注意)での山A、山Bの価値を$Q_t(A),Q_t(B)$とし、その初期値はどちらも0とする。更新するごとに価値を元とに生成した確率で山を選択し、出た数字によって選択した山に対する価値を更新する。選択した山から出た数字は$r_t$とする。確率は$P_t(A),P_t(B)$とし、$$P_{t+1}(A)=\frac{\exp(\beta Q_t(A))}{\exp(\beta Q_t(A))+\exp(\beta Q_t(B))}$$で計算する。価値の更新は$Q_{t+1}(X)=Q_t(X)+\alpha(r_t-Q_t(X))$によって行う。(時間発展はない系を考えているため理論で扱った第三項はない)
$\alpha=0.05,\beta=0.004$として、以上の設定を100回更新をした。価値は以下のように更新された。

それぞれ山を選ぶ確率は以下のように更新された。

コード

import numpy as np
import matplotlib.pyplot as plt

# seed値固定
np.random.seed(71)

#　引いた総数
N = 100

# 山A:0~2000
pileA = np.array([i for i in range(2001)])
# 山B:0~1000
pileB = np.array([i for i in range(1001)])

# パラメーター
alpha = 0.05
beta = 0.004

"""
変数
Q_A: (時間ごと)山Aの価値
Q_A: (時間ごと)山Bの価値
P_A: (時間ごと)山Aを引く確率
P_B: (時間ごと)山Bを引く確率
select_pile: (時間ごと)選択した山
"""
Q_A = [0]
Q_B = [0]
P_A = []
P_B = []
select_pile = []

for i in range(N):
    P_A.append(np.exp(beta * Q_A[i]) / (np.exp(beta * Q_A[i]) + np.exp(beta * Q_B[i])))
    P_B.append(np.exp(beta * Q_B[i]) / (np.exp(beta * Q_A[i]) + np.exp(beta * Q_B[i])))
    if P_A[i] >= np.random.rand():
        select_pile.append(1)
        Q_A.append(Q_A[i] + alpha * (np.random.choice(pileA, 1)[0] - Q_A[i]))
        Q_B.append(Q_B[i])
    else:
        select_pile.append(0)
        Q_A.append(Q_A[i])
        Q_B.append(Q_B[i] + alpha * (np.random.choice(pileB, 1)[0] - Q_B[i]))

# 時間
t = np.array([i for i in range(N)])

plt.scatter(t, Q_A[1:], c="blue", marker="o", label="A")
plt.scatter(t, 1000*np.ones(N), c="black", marker=".", s=3)
plt.scatter(t, Q_B[1:], c="red", marker="o", label="B")
plt.scatter(t, 500*np.ones(N), c="black", marker=".", s=3)
plt.title("Q")
plt.legend()
plt.savefig("value.png")
plt.show()

plt.scatter(t, select_pile, c="green", marker="s", s=10)
plt.plot(t, P_A, color=(1,0,0), marker="o", label="A")
plt.plot(t, P_B, color=(0,0,1), marker="o", label="B")
plt.title("prob")
plt.legend()
plt.savefig("prob.png")
plt.show()

時間発展のある系ではコードが大量になるのでここでは紹介しなかった。「cartpole qlearning」などで検索すると多くの記事が出てくるので気になる方はそこを見ていただきたい。

参考

理論:Sutton, Richard S. (1998). Reinforcement Learning: An Introduction.
問題引用:データ分析のための数理モデル

この記事はDP初心者向けになります。
そもそもDPって何って人は、けんちょんさんの記事をどうぞ！

最近DPの勉強をはじめて、
DPの次元を減らすテクニック
に非常に感動したので記事にしてみます！！！

以下の2問を実際に解きながら簡単に解説したいと思います。

• TDPC A - コンテスト →DPを2次元から1次元へ
• ABC 015 D - 高橋くんの苦悩 →DPを3次元から2次元へ

TDPC A - コンテスト

典型的な部分和問題です。
今回の問題はNのMAXが100
→bit全探索とかだと計算量的に無理なのでDPで考える。

DP:2次元Ver

普通に、DPを２次元として解くと、こんな感じー
※点数のパターンは0〜100点ではなく、0〜MAX10000(100*100)点であることに注意

DP_2D.py

import itertools,sys
def I(): return int(sys.stdin.readline().rstrip())
def LI(): return list(map(int,sys.stdin.readline().rstrip().split()))
#############################
N = I()
p = LI()

dp = [[0]*10001 for _ in range(N+1)] #0:その点数存在しない、1:その点数存在する
dp[0][0] = 1
for i,j in itertools.product(range(1,N+1),range(10001)):
    pi = p[i-1]
    if j-pi>=0:
        dp[i][j] = max(dp[i][j],dp[i-1][j-pi]) #選んだ場合
    dp[i][j] = max(dp[i][j],dp[i-1][j])        #選ばなかった場合

print(sum(dp[-1])) #1の数字の個数を出力

DP:1次元Ver

これをDP1次元にするにはこうする！

DP_1D.py

import itertools,sys
def I(): return int(sys.stdin.readline().rstrip())
def LI(): return list(map(int,sys.stdin.readline().rstrip().split()))
#############################
N = I()
p = LI()

dp = [0]*10001 #0:その点数存在しない、1:その点数存在する
dp[0] = 1
for i,j in itertools.product(range(1,N+1),range(10001)[::-1]):
    pi = p[i-1]
    if j-pi>=0:
        dp[j] = max(dp[j],dp[j-pi]) #選んだ場合
    #dp[j] = max(dp[j],dp[j])       #選ばなかった場合→不要

print(sum(dp)) #1の数字の個数を出力

ポイントはここ！！！
range(10001)[::-1]
※range(10000,-1,-1)でも可

逆からループを回すことで、dpを使い回すことができる！！！！！
これがテクニック！！！！！

もしdpを使いまわしながら、順番にループを回すと
例えばp=2点の場合に、

dp = [1,0,1,0,1,0,1,0,1 ・・・]

と4点、6点、8点・・・も存在することになってしまう・・・
なぜなら、
dp[j] = max(dp[j],dp[j-pi])
2点が1なので4点も1になる。。。
4点が1なので、6点も1になる。。。
・・・

しかし、逆からループさせると・・・

dp = [1,0,1,0,0,0,0,0,0 ・・・]

dp[j] = max(dp[j],dp[j-pi])

・・・
4点が0なので、6点も0になる！
2点が0なので、4点も0になる！
0点が1なので、2点は1になる！！！
ということでdpを使いまわしながら、正しい結果を得られることができる！！！！！！！！！
→dpの次元を2次元から1次元に落とすことができる！！！！！！！！！！

dpの次元が下がることで、
頭で考える量が減ってめちゃくちゃ嬉しい！

この感動伝わるかな〜
伝わってくれ〜

感動をより実感してもらうためにもう一問！！！

ABC 015 D - 高橋くんの苦悩

Difficulty:1388の水色問題！
ナップサックDPに変数が1つ増えた、少しだけ発展的な問題。

DP:3次元Ver

3次元でDPを考えたらできなくもないけど、頭がパンクしそう・・・・・・・・・・

DP_3D.py

import itertools,sys
def I(): return int(sys.stdin.readline().rstrip())
def LI(): return list(map(int,sys.stdin.readline().rstrip().split()))
#############################
W = I()
N,K = LI()
AB = [LI() for _ in range(N)]

dp = [[[0]*(W+1) for _ in range(K+1)] for _ in range(N)] #k枚使った時の重要度の最大値を保持
for i,k,w in itertools.product(range(N),range(1,K+1),range(W+1)):
    A,B = AB[i]
    if w-A>=0:
        dp[i][k][w] = max(dp[i][k][w],dp[i-1][k-1][w-A]+B)
    dp[i][k][w] = max(dp[i][k][w],dp[i-1][k][w])

print(dp[-1][-1][-1])

しかしこれではPythonだと計算量的にTLE、PyPyだとMLE
オワタ・・・

DP:2次元Ver

ここで！
さきほどの次元削減テクニックをつかって、
なんとかDPを使えまわすことができないか！
と考えてやってみると・・・
どうやら2次元にできそう・・・！

DP_2D.py

import itertools,sys
def I(): return int(sys.stdin.readline().rstrip())
def LI(): return list(map(int,sys.stdin.readline().rstrip().split()))
#############################
W = I()
N,K = LI()
AB = [LI() for _ in range(N)]

dp = [[0]*(W+1) for _ in range(K+1)] #k枚使った時の重要度の最大値を保持
for i,k,w in itertools.product(range(N),range(1,K+1)[::-1],range(W+1)[::-1]):
    A,B = AB[i]
    if w-A>=0:
        dp[k][w] = max(dp[k][w],dp[k-1][w-A]+B)
    #dp[k][w] = max(dp[k][w],dp[k][w])

print(dp[-1][-1])

ポイントはここ！！！
range(1,K+1)[::-1],range(W+1)[::-1]
逆からループで、DPを使いまわせる！！！

気になる結果は・・・

あいかわらず、PythonはTLEでしたがw(次元削減したからといって計算量が削減されるわけではなかったw)
PyPyで無事ACになりました〜！！！

今後DP問題が出た時に、次元が2次元以上になりそうだったら、
この次元削減テクニック使えないか考えてみよう〜

おわり！！！

緯度経度の座標情報をkepler.glで可視化する

kepler.glは地図上で様々なデータを可視化することができるツールである。今回はGPSデータを動画で可視化する。

手元にあるもの

gpxデータから作成したcsvファイル。
内訳はこのようになっている。

ID	時刻	経度	緯度	速度	trackID	segmentID	pointID
int	datetime	float	float	float	int	int	int

このうち、IDとtrackIDで分類すればトリップを一意に特定できるっぽい仕組みになっている。

求められているもの

公式ドキュメントによると、

How to use trip layer to animate path
Data format Currently trip layer support a special geoJSON format where the coordinate linestring has a 4th element denoting timestamp.
In order to animate the path, the geoJSON data needs to contain LineString in its features' geometry, and the coordinates in the LineString need to have 4 elements in the format of [longitude, latitude, altitude, timestamp], with the last element being a timestamp. Valid timestamp formats include unix in seconds such as 1564184363 or in milliseconds such as 1564184363000.

ということらしい。要約すると以下が必要である。

• GeoJsonファイル
• LineString形式
• 座標形式：[経度，緯度，標高，UNIX秒]

ありがたいことにサンプルデータもついており、

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "vendor":  "A"
      },
      "geometry": {
        "type": "LineString",
        "coordinates": [
          [-74.20986, 40.81773, 0, 1564184363],
          [-74.20987, 40.81765, 0, 1564184396],
          [-74.20998, 40.81746, 0, 1564184409]
        ]
      }
    }
  ]
}

この形式が求められている。

GeoJsonについて

だいたいJsonと同じらしい。Googleのロケーション履歴をダウンロードしようとするとJson形式で提供されるし、地理情報を含むJsonファイルを区別の意味を込めてGeoJsonと呼んでいるような雰囲気を感じた。実際、同じ内容で拡張子だけ.jsonにしたファイルも同様に読み込むことができた。

改めて先ほどのサンプルデータに注目すると、

• FeatureCollection
  • type
  • Features
  • type
  • properties
    • type
    • vendor
  • geometry
    • type
    • coordinates

という構造になっている。つまり、これを満たす構造にしてやれば読み込まれるということである。

DataFrameから座標リストの作成

各データの列名は以下の通り。（今回は標高データがないので速度データで代用している）

ID	時刻	経度	緯度	速度	trackID	segmentID	pointID
"ID"	"time"	"longitude"	"latitude"	"speed"	"trackID"	"segmentID"	"pointID"

まず、時刻をdatetime型からtimestamp型(UNIX秒)に変換する。

df["time"] = df["time"].apply(lambda x: x.timestamp())
#DataFrame

次にDataFrameから必要な項目だけを取り出す。リストで指定することで複数項目を一気に取得できる。

df_coordinate = df[["longitude","latitude","speed","time"]]
#DataFrame

次に、IDとtrackIDで分類できるようにする。

g = df_coordinate.groupby(["ID","trackID"])
#DataFrameGroupBy

次に各行(df[["longitude","latitude","speed","time"]]の一行ずつ)をリストにする。
coordinate_listにはdf_coordinateのID，trackID毎のリストのSeriesが作成される。

coordinate_list= g.apply(lambda d : d.values.tolist())
#Series ここ間違えやすい

最後にSeriesをlistに変換する。

coordinate_list = coordinate_list.values.tolist()
#list

これで、手元にリスト[経度，緯度，標高，UNIX秒]の1トリップ分のリストのリストができている。

• すべてのトリップの座リスト
  • 1トリップのリスト
  • 1点のリスト -　経度 -　緯度 -　標高（速度） -　時刻
    

という構造である。

GeoJson形式への変換

前項で作成したリストを用いて、GeoJsonデータを作成する。これにはgeojsonというライブラリを用いた。
素人目の認識としては、json書式に整えてくれるモジュールである。
先ほども書いた通り、求められるデータは以下の形式である。これらのデータを下から組み立てるように作成していくのが着実である。

• FeatureCollection
  • type（自己紹介）
  • Features
  • type（自己紹介）
  • properties（geometoryの詳細）
    • type（自己紹介）
    • vendor
  • geometry（地理データ）
    • type（自己紹介）
    • coordinates（LineString）

このモジュールは単純で、例えばFeatureを作りたければ以下のようにすればよい。
つまり、自己紹介部分や括弧などを勝手に追加してくれるというものである。

> Feature(geometry=my_point)  # doctest: +ELLIPSIS
{"geometry": {"coordinates": [-3.68..., 40.4...], "type": "Point"}, "properties": {}, "type": "Feature"}

この例文を参考に組み立てたものを以下に示す。

#Featureを蓄積するリストを作る
feature_list = list()

#FeatureリストにFeatureを蓄積させる
for i in range(len(coordinate_list)):
    #ptsに1トリップ分のリストを代入
    pts = coordinate_list[i]
    #lsにptsをGeoJson形式に変換したものを代入
    ls = geojson.LineString(pts)
    #lsをgeometryとし、propertiesと合わせてFeatureを作成
    #feature_listの末尾に追加
    feature_list.append(geojson.Feature(properties={'tripID': str(i)},geometry=ls))

#feature_listをFeatureCollectionに変換
Tripdata = geojson.FeatureCollection(feature_list)

これによってTripdataには求められる形式の平文が入ったものとなる。

ちなみにpropertiesにはFeature毎の識別情報を入れることが可能である。工夫次第では人ごとに分類、時間ごとに分類など様々な分類ができる。

書き出し

ここまで出来ればほぼ終わりである。
新たに拡張子に.geojsonをもつファイルを作成し、そこに平文を書き込むことでデータの完成となる。

path = './Tripdata.geojson'
with open(path, 'w') as f:
    f.write(geojson.dumps(Tripdata_dt_ts))
    f.close()

可視化

大げさに可視化と書いたが、先ほど作成したデータをドラッグ＆ドロップするだけで終了である。
propertiesの情報をもとに分類することも可能である。

Tipsとまとめ

すべてのデータについて、年月日だけを同じ日にしてやると時間帯のみの可視化ができるので、移動の多い時間帯が一目でわかりやすい。
変換はdatetimeのうちがやりやすいので、UNIX秒に変換する前にreplace()をするとよい。

df["time"] = df["time"].apply(lambda x: x.replace(year=2021,month=1,day=1))
df["time"] = df["time"].apply(lambda x: x.timestamp())

これら一連の流れでとても見栄えのいい可視化ができるので、ぜひ挑戦してほしい。

A. ある。

関連記事
https://qiita.com/sh1ma/items/3d9280084c340b4cd775

コードみてね

@lambda f: lambda *args: print(f(*args)) # ←動く！！(3.8まではこれ動かなかった)
def f():
    return "hello"

f() # -> hello

ウケるｗ

概要

最新のVideo Instance SegmentationタスクのベースラインとなっているモデルであるMaskTrackRCNNをABCIで動かそうとした際に、度重なるエラーに悩まされたので、ABCIにおける環境構築からモデルの学習までの方法を忘備録として残しておきます。
最後に格闘したエラーとその解決方法をまとめて残しておきました。

前提

• 2021/01/15 時点の情報である。
• ABCIのアカウントを持っており、サーバにリモートアクセスできる。
• anaconda3を導入済である。

環境構築

• ABCIにsshでリモートアクセス

• GPUノードの起動
  以下はGPUを最小構成で最大1時間インタラクティブに利用するコマンド。

  qrsh -g gcc50560 -l rt_G.small=1 -l h_rt=1:00:00 -m abes
  

• 必要なモジュールのロード
  立ち上げたノードで以下のコマンドを実行し、必要なモジュールをロードしておく。

  module load cuda/9.0
  module load cudnn/7.6/7.6.2
  module load nccl/2.3/2.3.7-1
  

• 仮想環境の作成とセットアップ
  以下のbash masktrackrcnn_env.sh を実行してレポジトリのクローン、仮想環境構築を行う。

  masktrackrcnn_env.sh

  # clone MaskTrackRCNN
  git clone https://github.com/youtubevos/MaskTrackRCNN.git
  cd MaskTrackRCNN
  
  # create environment
  conda create -n MaskTrackRCNN python=3.7 -y
  
  # >>> conda init >>>
  __conda_setup="$(CONDA_REPORT_ERRORS=false '$HOME/anaconda3/bin/conda' shell.bash hook 2> /dev/null)"
  if [ $? -eq 0 ]; then
      \eval "$__conda_setup"
  else
      if [ -f "$HOME/anaconda3/etc/profile.d/conda.sh" ]; then
          . "$HOME/anaconda3/etc/profile.d/conda.sh"
          CONDA_CHANGEPS1=false conda activate base
      else
          \export PATH="$PATH:$HOME/anaconda3/bin"
      fi
  fi
  unset __conda_setup
  # <<< conda init <<<
  
  # activate environment
  conda activate MaskTrackRCNN
  
  # setup environment
  conda install -c pytorch pytorch=0.4.1 cudatoolkit=9.0 torchvision -y
  conda install -c conda-forge  opencv -y
  conda install numpy cython -y
  conda install -c psi4 gcc-5 -y
  conda install libgcc -y
  pip install git+https://github.com/youtubevos/cocoapi.git#"egg=pycocotools&subdirectory=PythonAPI"
  bash compile.sh
  pip install . --user
  
  pip uninstall mmcv -y
  pip install mmcv==0.2.0
  

• ライブラリの修正
  ここまでで必要な環境は整ったが、このままでは学習時に不都合があるので最後に一つ修正を加える。
  [python_lib_path]/site-packages/mmcv/runner/checkpoint.py の39行目を以下に変更する。
  ※ [python_lib_path] には、which pythonが返すパスを入れれば良い。

  checkpoint.py

  print('While copying the parameter named {}, '
        'whose dimensions in the model are {} and '
        'whose dimensions in the checkpoint are {}.'
        .format(name, own_state[name].size(),
                param.size()))
  

• 環境構築完了!!

  データセットの準備

• データとラベルのダウンロード
  ここからデータとラベルをダウンロードしておく。

• データセットの配置
  以下の構成になるようにシンボリックリンクを貼る。

  MaskTrackRCNN
  ├── mmdet
  ├── tools
  ├── configs
  ├── data
  │   ├── train
  │   ├── val
  │   ├── annotations
  │   │   ├── instances_train_sub.json
  │   │   ├── instances_val_sub.json
  

  ラベルは元々train.json, valid.jsonというファイル名で与えられるので、（コピーして）改名しておくと良い。configのファイル名を変更しても良い。

  # シンボリックリンクの貼り方の一例 $MaskTrackRCNN は MaskTrackRCNNレポジトリのルートまでのパス
  mkdir $MaskTrackRCNN/data
  ln -s /path/to/original/data_dir/train $MaskTrackRCNN/data/train
  ln -s /path/to/original/data_dir/valid $MaskTrackRCNN/data/val
  ln -s /path/to/original/data_dir/annotations $MaskTrackRCNN/data/annotations
  

モデルの学習

• GCC7.4のロード module load gcc/7.4.0
• 仮想環境のアクティベート conda activate MaskTrackRCNN

• いざ学習
  python tools/train.py configs/masktrack_rcnn_r50_fpn_1x_youtubevos.py
  以下のようなログが出ていれば正しく学習されていると思われます。このままだと学習に24時間くらいかかりそうですが。

  2021-01-15 13:36:50,620 - INFO - Epoch [1][50/7669]    lr: 0.00199, time: 0.774, data_time: 0.045, loss_rpn_cls: 0.0609, loss_rpn_reg: 0.0465, loss_cls: 0.9336, acc: 84.5996, loss_reg: 0.2753, loss_match: 0.4937, match_acc: 89.2641, loss_mask: 0.7734, loss: 2.5835
  2021-01-15 13:37:27,818 - INFO - Epoch [1][100/7669]    lr: 0.00233, time: 0.744, data_time: 0.028, loss_rpn_cls: 0.0469, loss_rpn_reg: 0.0442, loss_cls: 0.7820, acc: 84.7695, loss_reg: 0.3567, loss_match: 0.2895, match_acc: 88.3957, loss_mask: 0.6092, loss: 2.1286
  2021-01-15 13:38:04,878 - INFO - Epoch [1][150/7669]    lr: 0.00266, time: 0.741, data_time: 0.026, loss_rpn_cls: 0.0342, loss_rpn_reg: 0.0342, loss_cls: 0.7171, acc: 85.1309, loss_reg: 0.3467, loss_match: 0.2588, match_acc: 89.5726, loss_mask: 0.5057, loss: 1.8968
  

  configs/masktrack_rcnn_r50_fpn_1x_youtubevos.py をいじればエポック数等学習のハイパーパラメータを変更できそうです。

格闘したエラーまとめ

• unable to execute 'nvcc': No such file or directory

  実行コマンド：bash masktrackrcnn_env.sh
  原因：cuda のロードがうまくできていない
  解決方法：module load cuda/{version} を実行すれば解決するはず。シェルスクリプト(masktrackrcnn_env.sh)中にこのコマンドを書いてもうまくいかなかったが、普通にシェルで実行したらうまくいった。

• AttributeError:module 'torch.nn' has no attribute 'SyncBatchNorm

  実行コマンド： python tools/train.py configs/masktrack_rcnn_r50_fpn_1x_youtubevos.py（学習用スクリプト）
  原因：torchのversionが1.1未満（でも著者はversion0.4.1で実行しているはず。）
  解決方法：mmcvのバージョンを0.2.0に下げる。

  pip uninstall mmcv
  pip install mmcv==0.2.0
  

  Traceback (most recent call last):
    File "tools/train.py", line 4, in <module>
      from mmcv import Config
    File "/home/acb11854zq/.local/lib/python3.7/site-packages/mmcv/__init__.py", line 4, in <module>
      from .fileio import *
    File "/home/acb11854zq/.local/lib/python3.7/site-packages/mmcv/fileio/__init__.py", line 4, in <module>
      from .io import dump, load, register_handler
    File "/home/acb11854zq/.local/lib/python3.7/site-packages/mmcv/fileio/io.py", line 4, in <module>
      from ..utils import is_list_of, is_str
    File "/home/acb11854zq/.local/lib/python3.7/site-packages/mmcv/utils/__init__.py", line 29, in <module>
      from .env import collect_env
    File "/home/acb11854zq/.local/lib/python3.7/site-packages/mmcv/utils/env.py", line 12, in <module>
      from .parrots_wrapper import get_build_config
    File "/home/acb11854zq/.local/lib/python3.7/site-packages/mmcv/utils/parrots_wrapper.py", line 79, in <module>
      _BatchNorm, _InstanceNorm, SyncBatchNorm_ = _get_norm()
    File "/home/acb11854zq/.local/lib/python3.7/site-packages/mmcv/utils/parrots_wrapper.py", line 71, in _get_norm
      SyncBatchNorm_ = torch.nn.SyncBatchNorm
  AttributeError: module 'torch.nn' has no attribute 'SyncBatchNorm'
  

• ImportError: /home/acb11854zq/.local/lib/python3.7/site-packages/mmdet/ops/nms/gpu_nms.cpython-37m-x86_64-linux-gnu.so: undefined symbol: __cudaPopCallConfiguration

  実行コマンド：python tools/train.py configs/masktrack_rcnn_r50_fpn_1x_youtubevos.py
  原因：ビルド時のCUDAのバージョン
  解決方法：ビルド時にCUDA9.2 -> CUDA9.0にすれば直る (source) 。この時、CUDA9.0はGCC<6.0を要求することに注意。

• RuntimeError: While copying the parameter named bbox_head.fc_cls.weight, whose dimensions in the model are torch.Size([41, 1024]) and whose dimensions in the checkpoint are torch.Size([81, 1024]).

  実行コマンド：python tools/train.py configs/masktrack_rcnn_r50_fpn_1x_youtubevos.py
  原因：モデルのbbox_head.fc_cls.weightの次元がチェックポイントのモデルのそれと違う。
  解決方法：[python_lib_path]/site-packages/mmcv/runner/checkpoint.py の39行目を以下に変更する（参考１、参考２）。

  checkpoint.py

  print('While copying the parameter named {}, '
        'whose dimensions in the model are {} and '
        'whose dimensions in the checkpoint are {}.'
        .format(name, own_state[name].size(),
                param.size()))
  

• Segmentation fault

  実行コマンド：python tools/train.py configs/masktrack_rcnn_r50_fpn_1x_youtubevos.py
  原因：issueに依るとbuild時にGCC<4.9.2だったことが問題らしい。
  解決方法：conda install -c psi4 gcc-5conda に gcc を入れる(source)。というのも、これの解決方法と矛盾しないように、4.9.2<=GCC<6.0 を使いたいけどABCIにその選択肢がない 。
  結果：ビルドには成功しているのでおそらくSegmentation faultの問題は解決できているが、新たにopencvのimport errorが出ている。

• ImportError: /home/acb11854zq/anaconda3/envs/test_gcc/bin/../lib/libstdc++.so.6: version `GLIBCXX_3.4.22' not found (required by /home/acb11854zq/.local/lib/python3.7/site-packages/cv2/cv2.cpython-37m-x86_64-linux-gnu.so)

  実行コマンド：python tools/train.py configs/masktrack_rcnn_r50_fpn_1x_youtubevos.py
  原因：condaのlibstdc++のバージョンが低い?
  解決方法： conda install libgcc でとりあえずこの問題は解決するっぽいが、今度はシステムのlibstdc++に関してimport errorが出る。

  Traceback (most recent call last):
    File "tools/train.py", line 4, in <module>
      from mmcv import Config
    File "/home/acb11854zq/anaconda3/envs/test_gcc/lib/python3.7/site-packages/mmcv/__init__.py", line 5, in <module>
      from .opencv_info import *
    File "/home/acb11854zq/anaconda3/envs/test_gcc/lib/python3.7/site-packages/mmcv/opencv_info.py", line 1, in <module>
      import cv2
    File "/home/acb11854zq/.local/lib/python3.7/site-packages/cv2/__init__.py", line 5, in <module>
      from .cv2 import *
  ImportError: /home/acb11854zq/anaconda3/envs/test_gcc/bin/../lib/libstdc++.so.6: version `GLIBCXX_3.4.22' not found (required by /home/acb11854zq/.local/lib/python3.7/site-packages/cv2/cv2.cpython-37m-x86_64-linux-gnu.so)
  

• ImportError: /lib64/libstdc++.so.6: version `CXXABI_1.3.8' not found (required by /home/acb11854zq/.local/lib/python3.7/site-packages/cv2/cv2.cpython-37m-x86_64-linux-gnu.so)

  実行コマンド：python tools/train.py configs/masktrack_rcnn_r50_fpn_1x_youtubevos.py
  原因：おそらくシステムのlibstdc++が古い？ CXXABI_1.3.8はどうやらGCC4.9から導入されている。
  解決方法：モデルの学習時にmodule load gcc/7.4.0としてシステムにGCC7.4をロードすることで無理矢理解消できる。bash compile.shによるビルドが終わってからでないとビルド時にエラーが出るので注意。

   Traceback (most recent call last):
      File "tools/train.py", line 4, in <module>
        from mmcv import Config
      File "/home/acb11854zq/.local/lib/python3.7/site-packages/mmcv/__init__.py", line 5, in <module>
        from .image import *
      File "/home/acb11854zq/.local/lib/python3.7/site-packages/mmcv/image/__init__.py", line 2, in <module>
        from .colorspace import (bgr2gray, bgr2hls, bgr2hsv, bgr2rgb, bgr2ycbcr,
      File "/home/acb11854zq/.local/lib/python3.7/site-packages/mmcv/image/colorspace.py", line 2, in <module>
        import cv2
      File "/home/acb11854zq/.local/lib/python3.7/site-packages/cv2/__init__.py", line 5, in <module>
        from .cv2 import *
    ImportError: /lib64/libstdc++.so.6: version `CXXABI_1.3.8' not found (required by /home/acb11854zq/.local/lib/python3.7/site-packages/cv2/cv2.cpython-37m-x86_64-linux-gnu.so)