pythonでjsonの配列にアクセスする時のメモ

jsonファイル

test.json

{
    "section1":{
        "key":"key1",
        "number": [ 
                    ["a", "b"],
                    ["c", "d"],
                    ["e", "f"]
                  ]
    },

    "section2":{
        "key":"key2",
        "number": 2
    }
}

pythonファイル

test.py

import json

json_open = open('test.json', 'r')
json_load = json.load(json_open)

print(json_load)
print(json_load["section1"]["number"][0])
print(json_load["section1"]["number"][1])
print(json_load["section1"]["number"][2])
print(json_load["section1"]["number"][0][0])
print(json_load["section1"]["number"][0][1])

実行

$ python test.py
{'section1': {'key': 'key1', 'number': [['a', 'b'], ['c', 'd'], ['e', 'f']]}, 'section2': {'key': 'key2', 'number': 2}}
['a', 'b']
['c', 'd']
['e', 'f']
a
b

参考

非エンジニアのためのJSON入門

この記事は・・・

pythonで地理空間情報を扱えるGDALをインストールする際に環境設定がうまく行かなかったようで、色々と対処をしたのでその記録です。
結果として、インストール中にしばらく待つだけで解決してたかもしれないので、参考にされる方は自己責任でお願いします。

インストールを行った環境

Anaconda3　ver.1.7.2
python 3.7.9

不具合の内容

まずはcondaコマンドを使ってインストールを試行しました。

conda install -c conda-forge gdal

その結果、「Solving environment： faild ~」となってインストールが完了しません。
確か5分ほどはインストールがすすまずに同じようなメッセージが続きました。

一旦インストールを中断して対処法を探しました。

対処法

対処1

・Anacondaのアップデート
これは当たり前ですかね。以下のコマンドを実行してAnacondaのすべてのパッケージを最新にアップデートしました。

conda update --all

しかし相変わらず同様のエラーが出てインストールが完了しない。

対処2

・環境変数の設定
以下の記事を参考にして環境変数を追加しました。

「anacondaでconda createを実行しようとすると「 Solving environment: failed 」と出て実行できなかった」
https://qiita.com/kizul/items/43a15a21346d121907c6

それでも相変わらずエラーが出ます・・・

対処法3

・バージョンを指定してインストール
以下の記事を参考にしました。

「【超初歩】AnacondaにGDALをインストールしてみた」
https://qiita.com/skperfarming/items/34b18d7cdce20982aa03

どうやらバージョンを指定したほうがいいらしい？
バージョンの確認は以下のコマンドで可能です。

conda search gdal

結果

本記事執筆時点（2020/11）では、3.0.2が最新のようです。
以下を実行する。

conda install gdal==3.0.2

これでも同じようなエラーが出ます。

結局・・・

対処法3を適用してもだめだったので、プロンプトをそのまま放置して更に対処法探しを再開します。
すると10分くらい経ったところで処理が進みました！

え？

インストールするときにある程度のところでプロンプトの表示が更新されるので、ログは拾えていませんが（真っ黒になって遡れなくなる）、インストール自体は完了したらしい。
動作確認のためにpythonコンソールで以下を実行。

>>>from osgeo import gdal
>>>

え？通った・・・？
あまり腑に落ちませんが、時間の問題ってだけで待っていればインストールできたのではないかと思います。
必要だったとしたらバージョンの指定で、環境変数は必要ないかもしれません。
ちょっと疑問は残りましたが、GDALが使えるようになったのでOKとしましょう。

この記事はTakumi Akashiro ひとり Advent Calendar 2020の3日目になります。

前置きの前置き

はい、Advent Calenderの3日目にして、もうネタが無くなりました。
最初の方にはちょっと実用的な記事を書く意識¹がありましたが
仕方ないので、ここから来週の月曜までは好き勝手調べた内容を書いてきます。
TAとかまったく関係ないです。

というわけで、3日目の記事、どうぞ！

前置き

皆さんはウィンドウの透過処理、使ったことありますか？自分は今までの仕事では一回も無いですね！
突然、ウィンドウにオーバーレイさせるアプリケーションとか作りたくなったので、ちょっと覚えていこうと思います(投げやり)

ウィンドウを半透過するには得手不得手があり、使い分けが必要だったので、
今日の記事はその辺を踏まえながら解説していきたいと思います！

本日のサンプルコード

#!python3
# encoding:utf-8

from PySide2 import QtWidgets
from PySide2 import QtCore

class Window(QtWidgets.QWidget):
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

        self.setWindowTitle("Sample")
        self.setStyleSheet("background-color:white;")
        self.setFixedSize(480, 320)

        hbox = QtWidgets.QHBoxLayout()
        vbox = QtWidgets.QVBoxLayout()
        hbox.addLayout(vbox)

        button_first = QtWidgets.QPushButton("Button: 1")
        button_first.setStyleSheet("background-color:gray; color:white;height:100%")
        vbox.addWidget(button_first)

        button_second = QtWidgets.QPushButton("Button: 2")
        button_second.setStyleSheet("background-color:darkgray; color:white;height:100%")
        vbox.addWidget(button_second)

        textedit = QtWidgets.QTextEdit("QtWidgets.QTextEdit")
        textedit.setStyleSheet("background-color:lightgray;")
        hbox.addWidget(textedit)

        self.setLayout(hbox)


if __name__ == "__main__":
    app = QtWidgets.QApplication()
    window = Window()
    window.show()
    exit(app.exec_())

はい、Qt製の何の変哲もない雑なウィンドウですね。
本日はこれを透過させてみたいと思います。

手法1 ウィンドウ背景の透明化オプションをセットする

self.setWindowFlagsに「WA_TranslucentBackground」をセットすることで¹ウィンドウの背景が透過するようになります。

ただしそのままでは、「WA_TranslucentBackground」をセットした時に
追加される背景を描画しない「WA_NoSystemBackground」の影響で背景色が描画されないので、
paintEvent側で背景色を描画する必要があります。

class Window(QtWidgets.QWidget):
    def __init__(self, *args, **kwargs):
        # ~中略~
        self.setAttribute(QtCore.Qt.WA_TranslucentBackground)
        self.setWindowFlags(QtCore.Qt.FramelessWindowHint)
        self.setStyleSheet("background-color: rgba(255, 255, 255, 127);")

    def paintEvent(self, event):
        # NOTE: トップレベルのWidgetはWA_TranslucentBackgroundフラグが立つと、
        #       背景が描画されなくなるので、paintEvent側で描画する。
        painter = QtGui.QPainter(self)
        painter.fillRect(0, 0, self.width(), self.height(), painter.background())

メリット	特定Widgetのみを透過できる
デメリット	ウィンドウのWidgetの透過には WindowFlagsのFramelessWindowHintが必要なため、 タイトルバーが消えてしまう

手法2 ウィンドウに直接透明度を適用する

class Window(QtWidgets.QWidget):
    def __init__(self, *args, **kwargs):
        # ~中略~
        self.setWindowOpacity(0.5)

メリット	めちゃくちゃ簡単に透過を実装できる
デメリット	Widget個別に透過度を設定できない

締め

評価ラベル	ランク(5段階)
おすすめ度	★
難易度	★
ニッチ	★★★★
汎用性	★
もやもや度	★★★★

(上記は加工した画像です)

今回、実はオーバーレイの件とは別件でこういう感じをやりたくて、調べていたんですがダメでした……

一応PySide2+Qtをビルドしなおして、QtWinモジュールを追加すれば、Win7時代によくあったガラスフレームが使えるようで、
上記の画像っぽい感じになるんですが…今回はそこまでの情熱はないです……

まあ手法1に自作のタイトルバーを追加すればいいだけなんで、問題はないんですけどね！！！！(精一杯の負け惜しみ)

というわけで個人的には若干消化不良な結果になりましたが、
今日の記事はここまでで終わりです。

ではでは！

---

1. そこら辺の話は土曜当たりの記事で書く予定です。 ↩

1. はじめに

みなさん、AWSしてますか？
値段も安くて、触りやすく、記事も豊富にあるので、僕はかなり好きです。

ただ、そうは言っても、個人で触りたいとなると料金の面で利用できるものとそうでない物が出てきたり、
一旦作成したサービスを忘れることなく全て削除したりするのは面倒な気がします。
また、WebのコンソールやCLIで作成した物は属人性が高くなり、手順やノウハウが残らないためよくないです。

そこで、便利なのがCDKです！

2. CDKとは

AWSが公式に提供するInfrastructure as Code(=IaC)のライブラリです（GitHub）。
CloudFormationをラップしており、TypescriptやPythonなどを用いて実装することで、サービスを作成することが可能です。

IaCとしては、前述のCloudFormationやTerraformなどもありますが、
これらはJsonで記述する必要があり（VSCodeだとプラグインで入力補助もある。）、
作成するサービスの量が多くなると見通しが悪くなる可能性があります。

2.1 CDKのメリット

• AWSが公式に提供しているので利用できなくなるなどの不安はない（と思う）
• ドキュメントをみながら必要なパラメータを付与していくことになるので、必要な設定などがわかりやすい
• スタックの名前を変えれば同じ内容のサービスを複数個展開できるため、開発と本番の環境分離なども容易に行える
  • ただし、同一の名前のサービスはできないので、環境ごとに名前を変える必要な場合がある（基本的には）
• リソースの作成・削除が容易（IaCのメリット）
• コードでサービスの管理ができるので属人性を排除できる（IaCのメリット）

2.2 CDKのデメリット

• リリースされたばかりのサービス（機能）は対応できないことがある
  • CloudFormationが対応されたのち、CDKの対応がされるため
  • （新サービスのLambda Containerは4時間くらいで更新されていたので、必ずしもそうとは言い切れなくなりましたね・・・笑）
• リリースの速度がかなり早く、過去に作成したコードが利用できなる可能性がある
  • 2020年12月07日時点で、6日に1度のペースでリリース
  • （とはいえ、破壊的な変更はそこまでない印象）
• ローカルからのコマンドでの削除が失敗することがある
  • IAMにCDKで管理していないPolicyを付与などした場合

3. さわってみる

メリット・デメリットは実際に触ってみて感じるのが一番だと思うので、
早速、利用方法を書いていきます。

3.1 環境構築

macの場合はbrewにてCDKをインストールできます。

$ brew install aws-cdk

その後、各プロジェクトに必要なCDKのPythonライブラリを入れていく形になります。

$ pip install aws-cdk.aws-lambda

必要になったら、都度パッケージを入れていくことになるのですが、
pypi aws_cdk athenaなどと必要なサービス名を含めて調べるとどのパッケージを入れればいいかわかります。

3.1 公式サンプルプロジェクト

ゼロから、このライブラリをドキュメントから利用方法を探っていくのはかなり大変だと思うので、その際に有用なのがサンプルプロジェクトです。

AWSの公式が容易したaws-cdk-exampleが本家大元なので、一番信頼できます。
ただ、Pythonのプロジェクトとなると意外と少なかったり、更新日時が古く動かないサンプルなどもあります。

3.2 作成したサンプルプロジェクト

そのため、私も同じようなフォーマットで、
Pythonのみのサンプルプロジェクト一覧aws-cdk-small-examplesを作成してみました。
記事作成時点はCDKのバージョンが1.76.0なのですが、その時点で全てのサンプルプログラムが動くようになっています。
作成してみたサンプルプロジェクトの数は12個です。
（本家公式は25個なので、半分程度はあります。ただ、かぶっているサンプルは2個程度のみです。）

対応しているサービスなどは以下の通りです。

• Lambda
• VPC
• AWS Batch
• Step Functions
• DynamoDB
• API Gateway
• CloudFront
• S3
• Glue

作成したリポジトリにて、公開しているサービスとその関連を表した図です。
（※注意：これらのサービスが一度にデプロイされるわけではありません。）

次の節から公開している、構成図の一部を紹介します。

3.2.1 API Gateway + Lambda + DynamoDB

この例では、3つのLambdaを作成し、DynamoDBにアクセスするapp.pyを公開しています。

3.2.2 Glue + StepFunctions

この例では、GlueのJobを定義して、それをStepFunctionsから実行するapp.pyを公開しています。
意外とGlueの例は少ないので、もし利用するような方がいれば参考になるかなと思います。

3.2.3 Lambda Container

ここでは、reInvent2020で発表された新サービスの1つである、Lambda Containerの例も作成しています。

細かい使い方などは別のQiitaの記事でも紹介しておりますで、合わせてどうぞ！

3.2.4 AWS Batch + StepFunctions

この例では、AWS BatchをStepFunctionsとCloudWatch Eventを利用して定期実行するapp.pyを公開しています。
VPCなども一から作成しており、結構長いコードとなっております。
細かい説明などは、Qiitaの記事にて、紹介しております！

3.2.5 NESTED AWS Batch + StepFunctions

CloudFormationがnestに対応したので、3.2.4の構成をnestを利用した例を作成しております。
nestを利用することのメリットとしては、作成したStackの責任を明確にできることが大きいです。

例では、以下のように分離させています。
* VPCのスタック：ベース
* Bathcのスタック：処理などを実際に行う部分
* StepFunctions + CloudWatch Eventのスタック：定期的に処理を走らせる部分

3.3 おすすめの記事

• API Reference

  • こちら本家のAPI referenceです。Python専用のdocsもありますが、こちらの方が参考になります。

• AWS CDKで快適なクラウド開発

• AWS CDK(Cloud Development Kit)の参考リンク集

4. おわりに

自分の中で、CDKはかなり好きなので、いろんな人にCDKを使っていただきたいなと思います！
IaCを進めるという意味でもかなり良いです！

また、これからも、随時GitHubのサンプルは増やしていく予定なので、ぜひチェックなどをお願いいたします！！

Pythonで音声ファイルのスペクトログラム描画をしてみました。とてもよくあるパターンではありますが、モノラル・ステレオ両対応になっているのが特徴です。

後でも述べますが、スペクトログラム描画により、人間が話しているところと、無音もしくは人間の会話周波数と異なるノイズとを、一目瞭然で区別することが可能です。例えば、話者識別の教師データを作る場合など、スペクログラムを見ることで、やりやすくなるでしょう。

なお、音声分析にはlibrosaを使う場合も多いようです。

細かい説明はコードに書かれていますが、以下概要です。

• pysoundfileでwaveファイルを読み込み
• 冒頭5分のみ切り出し（必須ではありませんが処理速度の関係で時間を短くしました）
• 高速フーリエ変換
• 7kHz未満の周波数のみ切り出し（人間の音声は3.4kHz未満がほとんどで、7kHz未満でも高音質）
• 音圧をdBに変換（相対dB）
• 結果を、縦軸が周波数・横軸が経過時間とし、dB値によるカラーマップを、グラフ描画（ステレオだとCHごとにグラフ描画）

以下がコードです。

import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf
from scipy import signal as sg

# 入力はwavファイルであること（mp3等は入力前にツールで変換する）
x, fs = sf.read('sample.wav')
x = x[:fs*300]   # 最初の5分間だけを分析する
# モノラル音声のndarrayをステレオ音声と同等形式に変換
if len(x.shape) == 1:
    x = x.reshape(-1, 1)

# stftのパラメータ
nperseg = 256
noverlap = nperseg // 2
# ウインドウがぴったりと切り分けられるように入力配列の後ろを少し削る
# stftしたデータを逆stftをするときに、この前処理が効いてくる
x = x[:-((len(x) - nperseg) % (nperseg - noverlap))]
# soundfileの形式をstftするには転置する
f, t, Zxx = sg.stft(x.T, fs=fs, nperseg=nperseg, noverlap=noverlap)

# 人間の声を分析するため7kHz以上の高音はカット
f = f[f < 7000]
Zxx = Zxx[:,:len(f)]
# スペクトログラムを描画するために相対デシベルを求める
dB = 20 * np.log10(np.abs(Zxx))

# スペクトログラムの描画
channels = dB.shape[0]
fig, ax = plt.subplots(channels)
if channels == 1:
    ax = (ax,)  #チャネル1つの時はaxがスカラーなので調整
for i in range(channels):
    ax[i].set_title(f'CH{i+1}-Time-Frequency')
    ax[i].set_xlabel('Time(s)')
    ax[i].set_ylabel('Frequency(Hz)')
    ax[i].pcolormesh(t, f, dB[i], shading='auto', cmap='jet')
plt.subplots_adjust(wspace=0.5, hspace=0.5)
plt.show()

出力結果

赤っぽいところが人が話をしている箇所、青っぽいところが無音です。会話なのか、会話ではないのかの区別は容易になりそうです。しかし、このグラフを人間が見ただけでは、話者の違いまでを判別するのは難しそうです。左右の短い話をしている人と、真ん中の長い話をしている人は、別人なのですが、その違いは、目で見てもわからないでしょう。

はじめに

機械翻訳のTransformerモデルをトレーニングする機会があり，Pytorchベースのfairseqを使ったんですが，アプリ用のコードでモデルのロードにハマってしまいました．備忘録のために記事を書きます．

fairseqとは

Facebookの人工知能研究チームが開発している，機械翻訳用のフレームワークです．Facebookが開発元ということもあり，Pytorchがベースになっています．最近はHuggingfaceのTransformersが人気でTransformerモデルを扱うならPytorchだよね，ということもあり，こちらをフレームワークとして選びました．
その他の機械翻訳フレームワークとしては，MarianNMT，OpenNMT（こちらもPytorchベース）などがあります．基本的な機能はどのフレームワークも大差ない印象ですが，論文実装のコードはfairseqが選ばれている場合が多いです．

モデルのトレーニング方法

fairseqはドキュメントが充実しているため，チュートリアル通りにやれば，自分でも機械翻訳モデルをトレーニングすることができます．チュートリアルは英語のみですが，日本語のドキュメントが不足しているのは，他の機械翻訳フレームワークでも同様です．

トレーニングのためのコード（fairseqの公式ドキュメントから）

> mkdir -p checkpoints/fconv
> CUDA_VISIBLE_DEVICES=0 fairseq-train data-bin/iwslt14.tokenized.de-en \
    --lr 0.25 --clip-norm 0.1 --dropout 0.2 --max-tokens 4000 \
    --arch fconv_iwslt_de_en --save-dir checkpoints/fconv

トレーニングコードのsave-dirは学習済みモデルの格納フォルダですが，このフォルダには，各エポックにおける学習済みモデル（checkpoint(n).pt），精度が一番よかった学習済みモデル（checkpoint_best.pt），最終エポックにおける学習済みモデル（checkpoint_last.pt）が格納されます．

checkpoints_fconv

$ ls -alF
checkpoint1.pt
checkpoint2.pt
...
checkpoint99.pt
checkpoint_bes.pt
checkpoint_last.pt

推論の際は，checkpoint_best.ptを読み込むことになるかと思います．

インタラクティブモードでの推論

fairseqの公式ドキュメントには，シェル上のインタラクティブモードでの推論についてはチュートリアルがありますが，Pythonのコード内でモデルをロードする方法などは紹介されていません．

fairseq-interactiveによ推論（fairseqの公式ドキュメントから）

> MODEL_DIR=wmt14.en-fr.fconv-py
> fairseq-interactive \
    --path $MODEL_DIR/model.pt $MODEL_DIR \
    --beam 5 --source-lang en --target-lang fr \
    --tokenizer moses \
    --bpe subword_nmt --bpe-codes $MODEL_DIR/bpecodes
| loading model(s) from wmt14.en-fr.fconv-py/model.pt
| [en] dictionary: 44206 types
| [fr] dictionary: 44463 types
| Type the input sentence and press return:
Why is it rare to discover new marine mammal species?
S-0     Why is it rare to discover new marine mam@@ mal species ?
H-0     -0.0643349438905716     Pourquoi est-il rare de découvrir de nouvelles espèces de mammifères marins?
P-0     -0.0763 -0.1849 -0.0956 -0.0946 -0.0735 -0.1150 -0.1301 -0.0042 -0.0321 -0.0171 -0.0052 -0.0062 -0.0015

Pythonのコード内で学習済みモデルをロードする方法

modelsクラスを使うのが正解でした．インタラクティブモードの場合と異なり，いくつかの引数は自動でロードされるようです．

>>> from fairseq.models.transformer import TransformerModel
>>> model = TransformerModel.from_pretrained('wmt14.en-fr.fconv-py', 'model.pt', 'wmt14.en-fr.fconv-py')
>>> text = 'Why is it rare to discover new marine mam@@ mal species ?'
>>> model.translate(text, beam=5)
'Pourquoi est @-@ il rare de découvrir de nouvelles espèces de mammifères marins ?'

まとめ

fairseqは日本語ドキュメントが少ないですが，HuggingfaceのTransformersを使っている人であれば使いやすいですし，Pytorchベースということにも将来性を感じます．
論文実装のコードも多いですし，今後も触っていきたいフレームワークですね．

1. はじめに

みなさん、Azureライフ楽しんでいますでしょうか？
私は、Azure FunctionsやDatabricksを主に利用し、とても楽しい毎日を送っております！！

とても、Azureが提供するサービスはどれも素晴らしいのですが、
Azure FunctionsからBlobStorageにアクセスしてファイル操作、
となると既存のライブラリ群ではちょっとラッパー関数が必要になって大変でした。
（もしかしたら、そんな操作をすることは少ないかもなのですが・・・。）

そこで、Azureから提供されているライブラリをラッパーして扱いやすいライブラリを作りました！
その名もAzFS (Azure File System)です。
リポジトリはこちらです。
（ただ、FileSystemとしての機能は2020年11月現在ではほぼないです。）

2. 利用方法

AzFSとは

Azure Storage Account上にあるcsvやparquet、pickleなどのファイルをpandasのDataFrameとして読み込めるライブラリです。
また、jsonの読み書きや、とあるディレクトリ以下のファイルのリストなども対応しています。

2.1 インストール

pipを利用してインストールできます

$ pip install azfs

2.2 インスタンスの作成

動かす環境によって、認証情報の取得方法が異なり、インスタンスに渡す引数が異なります。

現在対応している認証方法は、以下の2つです。

• Azure Active Directory
• Connection String

import azfs

# Azure Active Directory上で動かす場合は、`AzFileClient`への引数は不要です。
azc = azfs.AzFileClient()

# 接続文字列（connection string）を利用する場合
connection_string = "DefaultEndpointsProtocol=https;AccountName=xxxx;AccountKey=xxxx;EndpointSuffix=core.windows.net"
azc = azfs.AzFileClient(connection_string=connection_strin)

また、環境変数に以下の情報を入れてあげると、引数なしでも動くようになります。

2.3 ファイルのリスト

とあるディレクトリ以下のファイルリストを取得する方法は次の通りです。

path = "https://testazfs.blob.core.windows.net/test_caontainer"
file_list = azc.ls(path)

また、特定のファイル名に合致するファイルだけを取得するようなglob関数にも対応しています。

path_pattern = "https://testazfs.blob.core.windows.net/test_caontainer/*.csv"
file_list = azc.glob(path_pattern)

2.4 ファイルの読み込み

csvなどのファイルはpandas-likeとspark-likeの2通りの方法で読み込みことができます。
また、いずれの方法もpandasのDataFrameを返すようになっています。

機能	pandas-like	pyspark-like
ファイル読み込み	単一	複数
pandasの引数を渡す	O	O
読み込み時のフィルタ処理	X	O
並列読み込み	X	O
csv	azc.read_csv()	azc.read().csv()
parquet	azc.read_parquet()	azc.read().parquet()
pickle	azc.read_pickle()	azc.read().pickle()

2.4.1 pandas-likeな読み込み

pandasのメソッドのように記述し、DataFrameとして読み込むことができます。

# Blob Storageのpathを指定
path = "https://testazfs.blob.core.windows.net/test_caontainer/test_file.csv"
df = azc.read_csv(path)

# pandasの引数にも対応
df = azc.read_csv(path, skiprows=4, header=None)

# with句を利用して、pandasのメソッドのように振る舞うことも可能です。
# read_csv_az という別の関数名になってしまうところだけ注意
import pandas as pd
with azc:
    df = pd.read_csv_az(path)

2.4.2 spark-likeな読み込み

複数のデータを一度に読むことがある場合は、こちらの方法を利用した方が効率よく読み込むことが可能です。

path = "https://testazfs.blob.core.windows.net/test_caontainer/*.csv"
df = azc.read().csv(path)

# また、並列処理も引数を指定するだけです。
df = azc.read(use_mp=True).csv(path)

# 並列処理でファイルを読み込んでいる際に、各ファイルに処理を加えたい場合は、
# 関数を定義し、それを渡してあげると処理します。
def some_function(_df: pd.DataFrame) -> pd.DataFrame:
    _df['new_column'] = _df['price'] * 1.1
    return _df

df = azc.read(use_mp=True).apply(function=some_function).csv(path)

2.5 ファイルの書き込み

ファイルの書き込みに関しては、pandas-likeの方法しかないです。

# 書き込み
path = "https://testazfs.blob.core.windows.net/test_caontainer/test_written_file.csv"
azc.write_csv(df, path, index=False)

# readの時と同じようにwith句でpandasのDataFrameから直接関数を実行することも可能です。
# ただし、こちらも関数名がread_csv_azと微妙に異なる点のみ注意してください。
with azc:
    df.write_csv_az(path, index=False)

3. おわりに

冒頭でもお伝えしたように、AzFSは、FileSystemとしての機能はほぼないです。
そのため、Pull RequestやFileSystemにするためのissueを挙げてくれると大変嬉しいです！！

このライブラリを利用して多くの方がAzureライフを謳歌できたらいいいなと思います！

概要

• Python3の都合の良い実行環境が欲しくなったのでdockerで作成しました。
  • いつも使ってるOSがCentOS7だ（ubuntuではない。。）
  • ENVを切るのがめんどくさい。
  • 検証のたびにサーバを準備したくない。
  • 検証が終わったら消したい。

事前準備

• 下記をインストールし実行できること。
  • docker
  • docker-compose

ディレクトリ構成

./---docker-compose.yml
  |--Dockerfile
  |--requirements.txt   # pipインストール用 利用しなければ不要
  |--.env               # 環境変数用 利用しなければ不要
  |--src                # pyファイルを配置するディレクトリ
      |-- sample.py

各種ファイルの配置

docker-compose.yml

docker-compose.yml

version: '3'
services:
  python3:
    restart: always
    build: .
    container_name: 'python3'
    working_dir: '/root/src'
    tty: true
    volumes:
      - ./src:/root/src
    env_file:
      - .env

Dockerfile

FROM python:3
USER root

RUN apt-get update
RUN apt-get -y install locales && \
    localedef -f UTF-8 -i ja_JP ja_JP.UTF-8
ENV LANG ja_JP.UTF-8
ENV LANGUAGE ja_JP:ja
ENV LC_ALL ja_JP.UTF-8
ENV TZ JST-9
ENV TERM xterm

RUN apt-get install -y vim less
RUN pip install --upgrade pip
RUN pip install --upgrade setuptools
ADD requirements.txt /root/
RUN pip install -r /root/requirements.txt

requirements.txt

• pipインストールする必要があれば記載する

requirements.txt

requests

.env

• 環境変数を設定する必要があれば記載する

SAMPLE_USER ='root'

sample.py

sample.py

print("hello world")

実行

コンテナを起動

$ docker-compose up -d

スクリプトの実行

$  docker-compose exec python3 python3 sample.py
hello world

requirements.txtや.envを修正したとき

• 一度停止して再ビルド

$ docker-compose down
$ docker-compose build
$ docker-compose up -d

はじめに

SDFをサクっとCSVに変換するの逆バージョン。

なぜこんなことをしたいかというと、SMILESが含まれたCSVの化合物データに対し、SDF形式でしか表現できない情報（3次元座標や原子プロパティ等）を付加したいため。

仕様

• CSVの中のSMILESプロパティと、プロパティ項目を指定し、SDFとして出力
• 化学構造はSMILESを元に生成
• プロパティ項目は1つだけ指定可能(とりあえず機械学習の目的変数を１つ入れたいだけなので)

ソース

csv2sdf.py

import pandas as pd
from rdkit import Chem
import argparse


def main():

    parser = argparse.ArgumentParser()
    parser.add_argument("-input", type=str, required=True)
    parser.add_argument("-output", type=str, required=True)
    parser.add_argument("-smiles_col", type=str, default="SMILES")
    parser.add_argument("-target_col", type=str, required=True)
    args = parser.parse_args()

    # CSVの読み込み
    print(args.input)
    df = pd.read_csv(args.input)

    # SMILES列を取り出し
    smiles_list = df[args.smiles_col].values

    # target列の取得
    targets = df[args.target_col].values

    mols = []
    for i, (smiles, target) in enumerate(zip(smiles_list, targets)):
        # SMILESからmolを生成
        mol = Chem.MolFromSmiles(smiles)

        # Targetの値を設定
        mol.SetProp(args.target_col, str(target))

        mols.append(mol)

    # SDFへの出力
    sdf_wtr = Chem.SDWriter(args.output)

    for mol in mols:
        sdf_wtr.write(mol)

    sdf_wtr.close()

解説

特になし。
3次元化したいとか、原子プロパティを埋め込みたい等したい場合は、ソース中のRDKitのmolオブジェクトに色々操作を加えてみてください。

0. はじめに

Pythonでべき乗・累乗を計算するには２つの方法があります。
例えば $3^4$ を計算する場合、

• 3 ** 4
• pow(3,4)

のどちらでも計算できます。¹
本記事ではこれらの違いについて説明します。

内容は３章構成で以下の通りです。

第１章     ＊＊とpow()の違い

この章では主にPython初心者を対象に２つの方法の違いについて説明します。

第２章     速度の比較

この章では競技プログラミングなどの理由で計算速度を求める方を対象に２つの方法の実行時間を比較した結果を説明します。

第３章     Python内部のはなし

この章では第２章で見た結果の理由が知りたい方を対象にそれぞれの方法でPython内部で行われていることを説明します。

1. ＊＊とpow()の違い

この章ではPythonの基本事項を含め２つの方法の違いを説明します。

1.1. Python内での立場

１番の大きな違いはPython内での立場です。具体的にいうと

• **は演算子
• pow()は関数

です。

演算子は演算を行う指示で ** は「べき乗という演算をする」という意味があります。
+, %, !=, andなども全て演算子です。

関数は演算子などを用いて行う一連の処理をまとめたものです。pow()は「べき乗計算を行う」という関数なので以下のコードと一緒だと思っておけば良いでしょう。

def pow(a, b):
    result = a ** b

    return result

1.2. 第3の引数

少し発展的な内容となりますが、pow()関数には第3の引数があり、整数を指定することができます。
これを指定した場合、べき乗の結果を指定した整数で割ったあまりが得られます。指定しなかった場合はデフォルトでNoneが設定されているので前節で見た通りです。なお、第3の引数を指定する場合は第１、第２の引数も整数でなければなりません。
なのでpow()関数は実は以下のようなものです。²

def pow(a, b, c=None):
    result = a ** b
    if c != None:
        result = result % c  # %はあまりを求める演算子

    return result

1.3. まとめ

**とpow()は演算子と関数という違いがあります。
また、pow()は第3引数を指定することができ、指定した整数で割ったあまりを求めることができます。
このような違いがありますが、基本的にはどちらを使って書いても良いです。
例えば 「3の４乗を５で割ったあまり」は以下のように様々な書き方ができます。

# 以下のものは全て正しい答え(=1)が得られる 
3 ** 4 % 5
pow(3, 4) % 5
pow(3, 4, 5)

ただし、２番目のようにわざわざ関数と演算子を混ぜる必要はないので１番目か３番目のように書くといいでしょう。

2. 速度の比較

この章では**とpow()の実行速度を比較します。

2.1. 計測方法

正整数 $x, n, m$ について

• x ** n
• pow(x, n)
• x ** n % m
• pow(x, n, m)

の4種類を計算しました。
実行時間の計測にはAtCoderのコードテスト Python(3.8.2)を使用し、各計算法について５回の計測の平均をとりました。

計測に使用したコードは以下のものです。

def main():
    x = 3
    n = 10
    m = 17

    loop = 10000000

    for _ in range(loop):
      result = x ** n
      #result = pow(x, n)
      #result = x ** n % m
      #result = pow(x, n, m)
      #pass


if __name__ == "__main__":
    main()

2.2. 計測結果

結果は以下のようになりました。
$x=3, m = 17$ は固定で n と loop を変え計測しました。

なお、$loop = 10^7$ でfor文内がpassの場合は$196\rm{ms}$ でした。

（i）あまりをとらない場合

n	loop		x ** n	pow(x, n)
$10$	$1$		$21\rm{ms}$	$22\rm{ms}$
$10$	$10^7$		$1885\rm{ms}$	$2163\rm{ms}$
$10^7$	$1$		$2526\rm{ms}$	$2516\rm{ms}$
$10^3$	$2 \times10^6$		$2311\rm{ms}$	$2436\rm{ms}$

（ii）あまりをとる場合

n	loop		x ** n % m	pow(x, n, m)
$10$	$1$		$22\rm{ms}$	$22\rm{ms}$
$10$	$10^7$		$2074\rm{ms}$	$6319\rm{ms}$
$10^7$	$1$		$2553\rm{ms}$	$25\rm{ms}$
$10^3$	$2 \times10^6$		$3652\rm{ms}$	$1430\rm{ms}$

2.3. まとめ

結果をまとめると下図のようになります。
速度を重視したい場合に推奨する方法を示します。
空白部分はどちらでも大きな違いはありません。

当然ですが、べき乗計算にどちらの方法を使うか考えるよりもアルゴリズム自体を洗練させることに時間を使った方が良いです。べき乗計算の違いが重要になってくる場面は特定の場面(nが大きい & あまりをとる)を除いてほとんどないでしょう。
なのでこの結果はなんかちょっと違うらしいな程度にとどめて、実際はコードの統一性などで決めた方がいいと思います。

• n(正確には$x^n$)が大きい & あまりをとる $\rightarrow$ pow()を使う
• その他　$\rightarrow$ コードの統一性などで決める

3. Python内部のはなし

競技プログラミングを始めて間もない頃、べき乗計算について調べてみると次のような文言をよく目にしました。

---

Pythonの組み込み関数pow()は繰り返し二乗法を用いて実装されている

---

当時は「ふーんそうなんだ、便利だな」としか思わなかったのですが、これを確かめたことはないなと思い調べてみたのがこの記事を書いたきっかけです。
この章では実際にPythonのソースコードをみて**とpow()の違いについて調べた結果を説明します。
なお、底、指数は整数のみを想定します。

参考にしたソースコード
オリジナルの、C言語で実装されたいわゆるCPythonで、バージョンは3.8です。以下のリンクからみることができます。

https://github.com/python/cpython/tree/3.8

3.1. 調べる方法

Python標準のdisモジュールを使います。このモジュールを使うことで、コードを実行した際に内部で行われている処理を人間がわかる形でみることができます。
例えば、以下のようなコードを書きます。

test.py

x = 1
y = 2
x + y

そして

terminal

$ python -m dis test.py

と実行すると

terminal

  1           0 LOAD_CONST               0 (1)
              2 STORE_NAME               0 (x)

  2           4 LOAD_CONST               1 (2)
              6 STORE_NAME               1 (y)

  3           8 LOAD_NAME                0 (x)
             10 LOAD_NAME                1 (y)
             12 BINARY_ADD
             14 POP_TOP
             16 LOAD_CONST               2 (None)
             18 RETURN_VALUE

このような結果が表示されると思います。
ここから、BINARY_ADDというのが足し算をしているようだということがわかります。
あとはソースコードからBINARY_ADDの定義を探して追っていけば内部での処理がわかります。

3.2. ＊＊の場合

以下のコードを使って調べました。

x = 3
y = 4
x ** y

すると以下のようになります。

terminal

  1           0 LOAD_CONST               0 (3)
              2 STORE_NAME               0 (x)

  2           4 LOAD_CONST               1 (4)
              6 STORE_NAME               1 (y)

  3           8 LOAD_NAME                0 (x)
             10 LOAD_NAME                1 (y)
             12 BINARY_POWER
             14 POP_TOP
             16 LOAD_CONST               2 (None)
             18 RETURN_VALUE

よってBINARY_POWERというものを探します。
これを追っていくと次のようになります。

BINARY_POWER $\rightarrow$ PyNumber_Power $\rightarrow$ ternary_op $\rightarrow$ long_pow

3.3. pow()の場合

以下のコードを使って調べました。

x = 3
y = 4
pow(x, y)

すると以下のようになります。

terminal

  1           0 LOAD_CONST               0 (3)
              2 STORE_NAME               0 (x)

  2           4 LOAD_CONST               1 (4)
              6 STORE_NAME               1 (y)

  3           8 LOAD_NAME                2 (pow)
             10 LOAD_NAME                0 (x)
             12 LOAD_NAME                1 (y)
             14 CALL_FUNCTION            2
             16 POP_TOP
             18 LOAD_CONST               2 (None)
             20 RETURN_VALUE

pow()の場合はBINARY_POWERではなくCALL_FUNCTIONとなります。
これは、組み込み関数だけでなくコード中に自分で書いた関数など、全ての関数がCALL_FUNCTIONになります。

CALL_FUNCTIONからソースコードを追ってみると次のようになります。

CALL_FUNCTION $\rightarrow$ call_function $\rightarrow$ _PyObject_Vectorcall

そして_PyObject_Vectorcallのなかで以下のことが行われます。

\begin{aligned}
\rm{func} &= \rm{\_PyVectorcall\_Function}\\
&= \rm{builtin\_pow\_impl}\\
&= \rm{PyNumber\_Power}\\
&= \rm{ternary\_op}\\
&= \rm{\mathbf{long\_pow}}\\

\rm{result} &= \rm{func}()\\
&\rm{\_Py\_CheckFunctionResult}
\end{aligned}

3.4. long_pow

実はどちらの方法でべき乗計算をしようとしても結局long_powという関数にたどり着きます。

long_powをPythonで書いてみるとだいたい以下のようになります。(参照カウントや桁が大きい場合の処理などは省略しています。)

def long_pow(base, exp, mod): #base, exp, modはint
    negativeOutput = 0

    # 指数が負でmodがない場合はfloatの累乗へ移る
    if exp < 0 and mod == None:
        return float_pow(base, exp, mod)

    # modが指定されている場合の処理
    if mod != None:
        #mod = 0はエラー
        if mod == 0:
            raise ValueError

        #modが負ならnegativeOutputフラグを立ててmodを正にする
        if mod < 0:
            negativeOutput = 1
            mod = -mod

        # mod = 1なら必ず0
        if mod == 1:
            return 0

        # 指数が負なら指数を正にして底をmodを法とした逆元で置き換える
        if exp < 0:
            exp = -exp
            base = long_invmod(base, mod)

        # 底が負もしくはmodより大きいならmodで割った余りに置き換える
        if base < 0 or base > mod:
            _, base = divmod(base, mod)


    """
    ここまででmodが指定されている場合はbase, exp, modが非負になる。
    modが指定されていない場合はbaseが負の場合もある。
    """

    # modが指定されている場合、xをmodで割った余りに置き換える
    def REDUCE(x):
        if mod != None:
            _, x = divmod(x, mod)
        return x

    # xとyの掛け算。modが指定されている場合はあまりをとる
    def MULT(x, y):
        result = x * y
        result = REDUCE(result)
        return result

    # 右向きバイナリ法（繰り返し二乗法）
    j = 1 << exp.bit_length() - 1
    z = 1
    while j:
        z = MULT(z, z)
        if exp & j:
            z = MULT(z, base)
        j >>= 1

    # negativeOutputの場合は結果を負にする
    if negativeOutput and z != 0:
        z -= mod

    return z

以上より「Pythonの組み込み関数pow()は繰り返し二乗法を用いて実装されている」という文言は正しいですが誤解を含んでいると思います。（私が勝手に誤解していただけかもしれません。）

正しい
確かに組み込み関数pow()は繰り返し二乗法を用いて実装されています。

誤解(？)
pow() ""は"" 繰り返し二乗法を~ということは**は違う?
$\rightarrow$ **も繰り返し二乗法を用いて実装されています。
そもそも、どちらも同じ関数を使っています。

以降では、ここまででわかったことを基に２章の結果を説明します。

3.5. あまりをとらない場合

まず、あまりをとらない場合を見ていきます。
この場合どちらもlong_pow(x, n, None)が実行されるだけなので

• x ** n
• pow(x, n)

はアルゴリズムの面では完全に等しいです。したがって、$n$ がどのような大きさであっても２つに違いはありません。

１章で述べた通り、**は演算子、pow()は関数です。ここではこの違いが重要になります。

ここで3.2, 3.3でわかったPython内部の動きをまとめてみると下図のようになります。

この図からわかるように、関数を呼び出す際にはコスト(関数呼び出しのオーバーヘッド)を支払うことになります。
１回の呼び出しでは体感できる差は生まれませんが、呼び出し回数が増えれば大きな差になります。

このような理由から、

• n の大きさに依存しない
• loopが大きくなるほど演算子である**が関数であるpow()より速度面で優位になる

という2章の結果を説明できます。

3.6. あまりをとる場合

続いて、あまりをとる場合を見ていきます。
この場合も演算子と関数の違いからloopの大きさによって演算子の優位性が見えてきます。

一方、あまりをとらない場合と違うのはアルゴリズムの面でも２つは異なるという部分です。
x ** n % mではあまりを求めるのは最後だけなのに対し、pow(x, n, m)は乗算ごとにあまりを求めます。ここからはこれらの違いが n の大きさによってどのように影響するかを見ていきます。

n が小さい場合
この場合、「割り算のコストは小さくない」という事実が重要になってきます。あまりを求めるというのは割り算をすることに他ならないので何度もあまりを求めるpow(x, n, m)は一度しかあまりを求めないx ** n % mに比べて計算コストが大きくなります。

n が大きい場合
この場合、「大きな数の掛け算のコストが無視できない」という事実が重要になります。$x^n$ は指数関数ですので n が大きくなると $x^n$ は非常に大きくなります。そのため、一時的に $x^n$ そのものを求めるx ** n % mは常に m 未満の数の掛け算をするpow(x, n, m)に比べて計算コストが大きくなります。

これらをまとめると

• n が小さく loop が大きい場合は割り算と関数呼び出しのコストが支配的になりpow()の方が**より遅くなる
• n が大きく loop が小さい場合は大きな数の掛け算のコストが支配的になり**の方がpow()より遅くなる
• n と loop が共に大きい場合はこれらのせめぎ合いになるが、$x^n$ が指数関数だということをふまえると概ね**の方がpow()より遅くなる

というように２章の結果を説明できます。

3.7. まとめ

**とpow()は共にlong_powという関数につながっていて、long_powは繰り返し二乗法を用いて実装されています。
**とpow()の実行速度を比較するポイントは

• あまりをとらない場合は「関数呼び出しのオーバーヘッド」
• あまりをとる場合は加えて「割り算」と「大きな数の掛け算」

でこれらのコストのうちどれが支配的になるかを考えることによって判断できます。

4. おわりに

実用性はあまりなさそうな内容ですが、Pythonについて理解が深まったのでよかったです。

説明の間違いやバグ等ありましたらお知らせください。

---

1. 他にもmath.pow()やnumpy.power()など様々な方法がありますが本記事ではimportなしで扱える２つの方法のみを扱います。 ↩

2. 実際にはpow(a, b, c)と a ** b % cは内部での計算の仕組みが異なりますがa,bが２桁程度の小さな数ならば気にしなくていいです。気になる方は３章を見てください。 ↩

Autodesk Fusion 360にはいつもお世話になっています. (教育機関ライセンス無料なので)

別途生成した座標データを含むCSVがあり, それを点としてインポートする方法がわからなかったのでメモを残します.
自分の用途的に, 同一平面内の点のインポートがスコープです.

課題

Fusionにビルトインであるスクリプト(ImportSplineCSV.py)はCSVの座標をスプラインとしてインポートするもの > 余計なことをせずに点がいい

目的

Pythonスクリプトで, CSVに含まれる座標を点としてxy平面のスケッチにインポートする

入力するCSV

CSVの書式はカンマ区切りで, 1行に1つの点のx,y,z座標をいれる. (zは全て0)

input.csv

1.5,3.6,0
8.4,2.8,0
3.4,6.1,0

スクリプト

フローとしては, ダイアログで指定したCSVファイルに対して, カンマ区切りをばらして座標を抽出し, ポイントとしてスケッチに追加する, という作業を各行に対して行います.

下記コードを.pyファイルに記入し, 適当な場所に配置します.
(自分の場合は C:\Users\USERNAME\AppData\Roaming\Autodesk\Autodesk Fusion 360\API\Scripts あたり)

あとはFusion 360を開き, 画面上部から Tools > ADD-INS > Scripts and Ad-Ins... と進み, My Scripts欄のプラスボタンを押してファイルを保存したディレクトリを指定すればOKです.
ディレクトリ追加直後にはスクリプトが反映されず, ウィンドウを閉じたり開いたりしていたら表示されるようになりました.

ImportCSVPoints.py

import adsk.core, adsk.fusion, traceback
import io

def run(context):
    ui = None
    try:
        app = adsk.core.Application.get()
        ui  = app.userInterface
        product = app.activeProduct
        design = adsk.fusion.Design.cast(product)
        title = 'Import points from csv'
        if not design:
            ui.messageBox('No active Fusion design', title)
            return

        dlg = ui.createFileDialog()
        dlg.title = 'Open CSV File'
        dlg.filter = 'Comma Separated Values (*.csv);;All Files (*.*)'
        if dlg.showOpen() != adsk.core.DialogResults.DialogOK :
            return

        filename = dlg.filename
        with io.open(filename, 'r', encoding='utf-8-sig') as f:
            root = design.rootComponent
            sketch = root.sketches.add(root.xYConstructionPlane)
            line = f.readline()
            data = []
            while line:
                pntStrArr = line.split(',')
                for pntStr in pntStrArr:
                    try:
                        data.append(float(pntStr))
                    except:
                        break

                if len(data) >= 3 :
                    sketchPoints = sketch.sketchPoints
                    point = adsk.core.Point3D.create(data[0], data[1], data[2])
                    sketchPoint = sketchPoints.add(point)
                line = f.readline()
                data.clear()

    except:
        if ui:
            ui.messageBox('Failed:\n{}'.format(traceback.format_exc()))

※未検証なのでわかりませんが… スケッチにxy平面を指定しているので, z座標に何か数字が入ったら動かないかもしれません. あるいは3次元スケッチになるのかな…? あまり使いこなせていないのでよくわかってないです.
検証された方がいればコメントにてご教授ください.

余談

このスクリプトですが, 個数にして700弱の点をインポートするのに使っていますが, 問題なく動いています.

同様の方法でSketchCirclesメソッドを使ってCSVから丸をインポートしようとしましたが, 200ぐらいに数を減らしてもFusionがフリーズしてしまい使い物になりませんでした.
ここを見る限り一般的な問題のようで, 回避するにはスケッチを分けるといいようです.

私の用途としてはインポートした座標を元にバルクに穴を開けたいだけだったので, 円としてインポートしなくても穴径は指定できたのと, スケッチを分けると後工程が面倒になると思ったので思い切って点としてインポートしました.

はじめに

葉っぱや唇のようになだらかな山の曲線を描きたいと思った。ガウス関数（正規分布のアレ）も悪くはないが、両端をまっすぐにしたい。そう考えて思いついたのがベジェ曲線だった。

ベジェ曲線

matplotlibでグラフを描く

ベジェ曲線そのものを勉強したいわけではないので、基本的な部分は先人のプログラムを流用させていただく。
　　Pythonでベジエ曲線を描く（TadaoYamaokaの日記）

bezier_matplotlib.py

#import scipy.misc as scm       # 元サイトのこのインポートは古い
from scipy.special import comb  # 現在はこのcombを使う
import numpy as np
import matplotlib.pyplot as plt

def bernstein(n, i, t):
    return comb(n, i) * t**i * (1 - t)**(n-i)  # 元サイトから修正　scm.comb -> comb

def bezier(n, t, q):
    p = np.zeros(2)
    for i in range(n + 1):
        p += bernstein(n, i, t) * q[i]
    return p

# 当然、制御点は改変している
q = np.array([[0, 0], [3, 0], [4, 1],[5, 0], [8, 0]], dtype=np.float)

list = []
for t in np.linspace(0, 1, 100):
    list.append(bezier(len(q)-1, t, q))  # 元サイトから改良　3を決め打ち -> 制御点の数に従う
P = np.array(list)

plt.plot(P.T[0], P.T[1])
plt.plot(q.T[0], q.T[1], '--o')
plt.show()

両端が水平の山形を作ることができた。線分の両端（0番目と4番目）が決まったら1番目と3番目の制御点を線分上で適当に決めてやれば、自由なのはただ一つ2番目だけとなる。

OpenCVで自由に描く

GUIでxy平面上に自由にこの曲線を描けるようにする。
せっかく細かい計算をしているのにOpenCVにすると整数上の話になってしまうのが残念なところではあるが。
マウスイベントについてはこちらを参照のこと。

bezier_opencv.py

from scipy.special import comb
import numpy as np
import cv2

def bernstein(n, i, t):
    return comb(n, i) * t**i * (1 - t)**(n-i)

def bezier(n, t, pts):
    p = np.zeros(2)
    for i in range(n + 1):
        p += bernstein(n, i, t) * pts[i]
    return p

def MouseEvent(event, x, y, flag, params):
    global pts, cnt, clicking
    # 左クリックの挙動
    if event == cv2.EVENT_LBUTTONDOWN:
        if cnt == 0:  # 最初の点を決める
            cnt += 1
            pts[0] = (x,y)

        elif cnt == 1 and (x,y) is not pts[0]:  # 線分を確定すると同時に制御点も決める
            cnt += 1
            pts[4] = (x,y)
            pts[1] = pts[0]*0.7 + pts[4]*0.3
            pts[2] = pts[0]*0.5 + pts[4]*0.5
            pts[3] = pts[0]*0.3 + pts[4]*0.7
            draw_bezier(pts)

        else:  # 線分が確定したらドラッグモード
            clicking = True

    # 左ボタンから手を離す
    if event == cv2.EVENT_LBUTTONUP:
        clicking = False

    # 左クリック（ドラッグ）中はその座標を中央の制御点とする
    if clicking:
        pts[2] = (x, y)
        draw_bezier(pts)

    # 右クリックでベジェ曲線をリセットする
    if event == cv2.EVENT_RBUTTONDOWN:
        pts[2] = pts[0]*0.5 + pts[1]*0.5
        draw_bezier(pts)

def draw_bezier(pts):
    img = img1.copy()
    curve1 = []
    for t in np.linspace(0, 1, 50):
        curve1.append(bezier(4, t, pts))
    curve = np.array(curve1, np.int32)
    cv2.polylines(img, [curve], False, (255,0,0), thickness = 5)
    cv2.imshow(winname, img)

winname = "kupa"
img1 = np.full((240,320,3), (0,128,128), np.uint8)
pts = np.zeros((5, 2), dtype=np.float32)
cnt = 0
clicking = False
cv2.imshow(winname, img1)
cv2.setMouseCallback(winname, MouseEvent)
cv2.waitKey(0)
cv2.destroyAllWindows()

葉っぱ曲線

幾何の問題

話は変わるが、高校数学で「直線に対して対称な点の座標を求める」という問題がよくあった。だが、そのほとんどは直線や点の座標が具体的に明記されており、汎用化された式を目にすることはなかった。

そこで頑張って解いてみた。
答えは、

X = (2x_1-a) - 2(x_1-x_2)t \\
Y = (2y_1-b) - 2(y_1-y_2)t \\
ここで　t = \frac{(x_1-x_2)(x_1-a) + (y_1-y_2)(y_1-b)}{(x_1-x_2)^2+(y_1-y_2)^2}\\

だ。これから頑張ってtを消すと複雑かつ美しい式が現れるが、数学の問題ではなくプログラムとして実装していく話なのでこの形で十分だろう。それにしてもここまで複雑な式になるとは思わなかった。道理で高校数学の問題としては出てこないわけだ。

葉っぱ曲線を実装する

ソースは略。あとで完成形を紹介します。

応用

閉曲線の内部を別画像にすればいろいろ面白いことができる。
葉っぱとか唇とか書いたが実は私が想定していたのは八雲紫だった。スキマ妖怪に限らずとも空間に生じた裂け目と考えれば、そこから地獄が顔をのぞかせていたり宇宙とつながっていたり巨大な目玉がこちらを睨みつけていたりと様々な妄想をカタチにすることができる。必要とあらば（？）「くぱぁ」も可能だ。

ここでは単一の画像で簡単に実装できるそこそこエモい例として、閉曲線の中は色を反転させるようにした。遊んでいるうちに魔空空間とかポドリアルスペースとかいろいろ懐かしい言葉が思い浮かんできて涙が止まらなくなった。

ソース

kupa.py

from scipy.special import comb
import numpy as np
import cv2

def bernstein(n, i, t):
    # bezier_opencv.pyと同

def bezier(n, t, pts):
    # bezier_opencv.pyと同

def MouseEvent(event, x, y, flag, params):
    # bezier_opencv.pyと同

def get_symmetry_points(pts):
    x1, y1 = pts[0]
    x2, y2 = pts[4]
    a, b = pts[2]
    t = ((x1-x2)*(x1-a)+(y1-y2)*(y1-b))/((x1-x2)**2+(y1-y2)**2)
    X = (2*x1-a)-2*(x1-x2)*t
    Y = (2*y1-b)-2*(y1-y2)*t
    q = pts.copy()
    q[2] = (X, Y)
    return q

def draw_inner_kupa(img1, img2):
    # 元画像に(0,0,0)という色があるとおかしくなる
    img1 = np.where(img1==(0,0,0), img2, img1)
    return img1

def draw_bezier(pts):
    # ptsと両端の線分で対象な点を持つ座標群を作る
    # ptsと対称だからといってqtsという変数名は我ながらどうかと思う    
    qts = get_symmetry_points(pts)

    # ptsおよびqtsを制御点とするベジェ曲線を作る
    curve1, curve2 = [], []
    for t in np.linspace(0, 1, 50):
        curve1.append(bezier(4, t, pts))
        curve2.append(bezier(4, t, qts))

    # curve2の順序を逆にし、curve1とcurve2を連結して閉曲線とする
    curve = np.concatenate((curve1, curve2[::-1]), axis =0)

    curve = np.array(curve, np.int32)
    img = img1.copy()
    cv2.fillPoly(img, [curve], (0,0,0))
    img = draw_inner_kupa(img, img2)
    cv2.imshow(winname, img)

winname = "kupa"
filename = "hoge.jpg"
img1 = cv2.imread(filename)
img2 = 255 - img1
pts = np.zeros((5, 2), dtype=np.float32)
cnt = 0
clicking = False
cv2.imshow(winname, img1)
cv2.setMouseCallback(winname, MouseEvent)
cv2.waitKey(0)
cv2.destroyAllWindows()

実行結果

元画像は新宿の高層ビルの写真素材（ぱくたそ）。いろいろなサイトで使われてるな、この画像。

終わりに

「ギャバンやウィングマンを例に出してるけどそういう演出なら○○のほうが適切ですよ」といったマサカリ、お待ちしています。

この記事は jupyter notebook を 日本語を含む pdf に書き出す方法についてのものです．

環境

	version
OS	macOS Big Sur 11.0.1 (Catalina でも動作確認済)
jupyter notebook	6.1.4
jupyter nbconvert	6.0.7
jupyter nbextensions	6.1.4
texlive	2020 (2017 でも動作確認済)
pandoc	2.11

nbconvert をインストールする

conda や pip を使って nbconvert をインストールする．
詳しくは公式ページを参照．

テンプレートを作る

デフォルトでは，notebook を pdf に書き出す際に日本語が埋め込まれない．
これは，.ipynb -> .tex -> .pdf の順に変換を行っているためである．
そのため， latex に変換する際に日本語用のテンプレートを用意する必要がある．

以下のような jinja2 のテンプレートファイル (.tex.j2 多分拡張子はなんでも良い) を作る．名前はなんでもよいが，以下ではja_template.tex.j2 とする．

ja_template.tex.j2(texlive2020)

((*- extends 'index.tex.j2' -*))

%===============================================================================
% Latex Article
%===============================================================================

((* block docclass *))
\documentclass[xelatex, ja=standard, 11pt]{bxjsarticle}
((* endblock docclass *))

ja_template.tex.j2(texlive2017)

((*- extends 'index.tex.j2' -*))

%===============================================================================
% Latex Article
%===============================================================================

((* block docclass *))
\documentclass[xelatex, 11pt]{bxjsarticle}

\usepackage{zxjatype}
\setjamainfont{ipaexm.ttf}
\setjasansfont{ipaexg.ttf}
\setjamonofont{ipaexg.ttf}
((* endblock docclass *))

texlive 2017 の場合は, ja=standard だと font-not-found のエラーがでた．この場合上記のように対処するとうまくいく．詳しくはtexwikiを参照．

※ 本当は docclass の block の中にフォントの設定を書くのはよろしくないと思いますが，うまく行ってしまったのでちゃんとしたやり方を調べていません...
テンプレートについて詳しく知りたい方は公式ページを参照してください．

ファイルの保存場所カレントディレクトリや ~/.jupyter/ などどこでもよい．

ためしてみる

コマンドラインで以下を実行してうまくいくか確かめる． --no-input はコードを隠すオプションのためお好みで．

jupyter nbconvert --to pdf --no-input --template-file path/to/ja_template.tex.j2 hoge.ipynb

ワンクリックで pdf 化できるようにする

まず nbextensions を入れる．詳しくは公式ページを参照．

jupyter を開き，Nbextensions の設定画面を開く．
Printview を有効にし，nbconvert options に以下を入力する．

--to pdf --no-input --template-file path/to/ja_template.tex.j2

あとは pdf 化したい notebook を開いた状態で，プリントマークをクリックしたら OK.

目的：kaggle学習中に得た新しい知識を記録、再利用
内容：kaggleやプログラミングの知識
備考：はじめのうちは無秩序にメモを保存していく。数が増えてきたら整理したい

【notebookで自作モジュールを扱う】
実行中ipynbファイルがあるフォルダか、pip,condaなどのインストールフォルダにimportしたいパッケージがない場合、sysでパスを通してあげる。
import sys
sys.path.append('/Users/~~~')

【EDAはじめにやること】
・train,testにデータをいれる
・'.shape'で要素数、カラム数を確認
・'.dtypes'で各カラムの型を確認
・'.head()'でdfの内容をざっと確認
・目的変数に対し'.describe()'で統計量確認
・目的変数に対し'.plot.hist()'で分布を確認

【

本日の学習　9:00～12:00 15:00～18：30

昨日に引き続きprogate 有料版　python編をしています。
午前中に最後まで行ったのでもう一周しています。progateを何周もするのは良くないと聞きますが　関数、クラス、継承　などVBAでやらなかったところは、理解が怪しいので2周と言わず3周するかも。
時間がなければ飛ばしたほうが良いのだろうが昨日youtubeでシリコンバレーのエンジニアさんが基礎だけでも３か月はかかるみたいなこと言ってたので安心してじっくり取り掛かりたいと思う。ただ寝てても書けるような事はしないつもりです。
短期的な目標としてはスクレイピングをして作りたいツールがあるのだけど、先日買ったudemyのスクレイピングの講座でHTML、CSSがわからない人はprogateで勉強しろみたいな事書いてたので遠回りになりますがタイミングみてさくっと勉強してみようと思います。
やること多すぎますが楽しんでやってるので良しとしましょう！！

一日の半分をパソコンの前で過ごす生活を一週間続けていたら、両手が腱鞘炎になりました。

はじめに

生化学の小テスト前！
対策問題を見返してみたけど、『DNA』しかわからない！
(あすぱらぎん酸...?めちおにん????おいしい？？)

勉強しなきゃ。。。。

うぅぅ。PC依存症すぎて、手がキーボードから離れない。。。
いや、待てよ。

てきとーにPythonで一問一答作ったら、勉強した気になれるだろ！

???「今回の小テスト、割とまじで終了だねwww」

今回したかったこと

・友達がQiitaデビュー&初投稿していて、なんでもいいから久しぶりにQiitaに投稿したかった。
・嫁(MacBook)に問題を出してもらいたかった。
・今後も使える一問一答ツールを作りたかった。

???「今回の小テストはマジでそんなことしてる場合じゃない」

→毎回の小テストの度に関数を追加したら、最強の一問一答ツールに(?)！

嫁(MacBook)「GitHub使って、クラスのみんなで問題を作れば毎回のテストで勉強できるね！」

そんな目標・モチベを持って、第1回目の関数を記録しておきました。
今回の作成関数：
「順不同の用語を答える問題。すべて正しければ正解とする。」

一問一答ツールプロジェクト第1回目

QandA.py

def gudge(inp, ans):
    inp_list = inp.split("、")
    if len(inp_list) != len(ans):
        print("不正解。（回答数が足りません。）\n答え：")
        print(ans)

    else:
        for i in range(len(ans)-1):
            if ans[i] in inp_list:
                print("正解：", ans[i])
                inp_list.remove(ans[i])
            else:
                print("不正解。\n答え：")
                print(ans)
                break
        print("正解")
    return

def question(que, num, ans):
    print(que)
    num_tsu = str(num) + "つ："
    inp = input(num_tsu)
    gudge(inp, ans)
    return

def q_ex1():
    que = "中性アミド酸。5つの中分類に分けられる。\n( )アミノ酸、( )アミノ酸、( )アミノ酸、( )、( )の5つである。"
    num = 5
    ans = ["イミノ酸","芳香族","含硫","酸アミド","脂肪族"]
    question(que, num, ans)

if __name__ == "__main__":
    q_ex1():

まとめ

???「ふつうにプリントで勉強しろ。」

(『一問一答ツールプロジェクト』第2回目は未定です。だれか作ってください。)

Result

作業	時間
コーディング	5分
デバック	5分
Qiita	50分
合計	60分

俺の一時間返して。（泣）

はじめに

グラフの可視化のツールは様々なものがあるが、graph-toolはその中でも最も強力といっても過言ではないだろう。
いろいろな可視化の方法が実装されているだけでなく、OpenMP化したC++で実装されているのでかなり大きなネットワークでも描画することができる。
例えばこのような可視化ができる。( https://graph-tool.skewed.de/ より )

ここではgraph-toolを用いた可視化について典型的な処理をまとめていく。

reference

https://graph-tool.skewed.de/static/doc/draw.html

動作環境

動作確認環境は以下の通り

• python 3.9.0
• graph-tool 2.35
• numpy, matplotlib

以後のサンプルコードでは、前提条件として以下のモジュールをimportしている。

import graph_tool.all as gt
import numpy as np
import matplotlib.pyplot as plt

またサンプルデータとして以下のものを用いる。（ある小説における名詞と形容詞の隣接関係から構築したネットワーク）

g = gt.collection.data['adjnoun']

Jupyter notebookで実行するかターミナルから実行するか

graph-toolで可視化するコードを実行する場合、jupyter notebookで実行するかターミナルで実行するかに依存して若干挙動が異なる。
可視化するスクリプトをjupyterから実行した場合、ノートブック内にインラインで表示される。

一方でスクリプトをターミナルから実行した場合、interactiveなGUI windowに描画される。ノードをマウスでドラッグ＆ドロップして動かしたり、選択したノードとそのエッジをハイライトしたりできる。用途に応じて使い分けたい。

コード

とりあえず描画する

gt.graph_drawメソッドを呼べば良い

gt.graph_draw( g )

ファイル出力する

outputオプションにファイル名を指定する。拡張子に応じてフォーマットを決めてくれる。
pdf, svg, png, jpgなど様々な形式に対応している。

gt.graph_draw(g, output = 'graph.pdf')

pngなどのラスタ形式の画像で画像の大きさを変えたい場合は、output_size オプションを指定する

gt.graph_draw(g, output = 'graph.png', output_size=(1000,1000))

ノードにラベルを表示する

VertexPropertyで指定する文字列を表示するにはvertex_textオプションを指定する。

gt.graph_draw(g, vertex_text = g.vp.label)

ノードの外側に表示したいときはvertex_text_position = 0を指定する。

gt.graph_draw(g, vertex_text = g.vp.label, vertex_text_position = 0)

その他、オプションをいじることでテキストの色やフォントなど細かく調整できる。

gt.graph_draw(g, vertex_text = g.vp.label, vertex_text_position = 0, vertex_font_size = 16, vertex_text_offset = [-1,1], vertex_text_color=(255,0,255,255))

レイアウトを変更する

graph-toolには様々なlayoutアルゴリズムが用意されている。まずレイアウトを計算し、得られた位置情報をgraph_drawの引数に渡す。

pos = gt.radial_tree_layout(g, 0)
gt.graph_draw(g, pos = pos)

ノードやエッジの見た目を次数などの統計量に応じて変更する

ノードやエッジのサイズや色を次数や他のネットワーク統計量に応じて変更することができる。
例えば、次数に応じてノードの大きさを変更するには以下のようにする。

deg = g.degree_property_map("total")
gt.graph_draw(g, vertex_size = deg)

スケールをうまく調整するには gt.prop_to_size メソッドを使うと良い。指定した最小値、最大値になるように適切にスケーリングしてくれる。logスケールにすることもできる。

gt.graph_draw(g, vertex_size = gt.prop_to_size(deg,mi=5,ma=20,log=False) )

工夫すれば他にも幅や色味など色々と変更できる。

deg = g.degree_property_map("total")
ebet = gt.betweenness(g)[1]   # edge betweeness
gt.graph_draw(g,
              vertex_size=gt.prop_to_size(deg,mi=5,ma=20,log=False),
              vertex_fill_color=gt.prop_to_size(deg,log=False),
              vcmap=(matplotlib.cm.inferno, .9),
              vorder=deg,
              edge_color=gt.prop_to_size(ebet),
              eorder=ebet,
              edge_pen_width=gt.prop_to_size(ebet,mi=1,ma=5,log=False),
              ecmap=(matplotlib.cm.bone, .7)
             )

最適解で何を求めたいか？

みんな大好きBostonの家の価格のデータセットがありますよね。このデータセットを使って、最適解を求めたいと思いました。

カラム	説明
CRIM	犯罪率
ZN	住宅の密集度
INDUS	非小売業の割合
CHAS	チャーリーズ川ダミー変数（川の周辺1、それ以外0）
NOX	酸化窒素濃度
RM	1住戸あたりの平均部屋数
AGE	1940年以前に建設された物件割合
DIS	5つのボストンの雇用センターまでの距離
RAD	大きな道路へのアクセス距離
TAX	10,000ドルあたりの税率
PTRATIO	教師あたりの生徒数
B	町における黒人の割合
LSTAT	人口当たり給料の低い率
MEDV	住宅の中央値

例えば、家の価格（MEDV）が安くて、部屋の数（RM）が多いところに住みたいよね。部屋数が多くて価格が安いのが皆さんうれしいですよね。この最適解モデルができれば、SUUMOとかホームメイトとか、賃貸業者に売れるかもしれません。
とりあえずボストンのデータセットをもとに、遺伝的アルゴリズムでパレート最適解を求めてみよう。

特徴量と目的変数、データの可視化

目的変数が、（MEDV,RM）で、それ以外を特徴量にしました。
その特徴量で、それぞれ目的変数として、サポートベクター回帰を使いMEDVとRMの学習モデルを作りました。

モデル１として、家の価格(MEDV)の予測モデルを作ります。

特徴量	目的変数
CRIM,ZN,INDUS,CHAS,,NOX,AGE,DIS,RAD,TAX,PTRATIO,B,LSTAT	MEDV

モデル２として部屋の数（RM）の予測モデルを作ります。

特徴量	目的変数
CRIM,ZN,INDUS,CHAS,NOX,AGE,DIS,RAD,TAX,PTRATIO,B,LSTAT	,RM

まず家の価格と部屋の数は結構相関があるんですよね。だから、これをお互いのモデルを作るときに、特徴量からはずすと、精度があがらないような気がします。

後、見てわかるように、家の価格の５０のところにグラフの上ではデータが縦に並んでいます。そこで、家の価格５０以上は除きました。

SVRでの予測で価格と部屋の数の予測結果

まずSVRで学習モデルを作りましたが。。。
ここのアルゴリズムは何でもよくて、サポートベクターでもランダムフォレストでもお好きなもので。
サポートベクターのハイパーパラメーターは、optunaで決定しました。

学習モデル１は、optimised_svr_MEDVに、学習モデル２は、optimised_svr_RMに入れました。

予測結果はこういう感じです。
MEDVの学習精度:0.873
MEDVのテスト精度:0.748
RMの精度:0.628
RMの精度:0.439

確かに家の価格は、部屋数がなくても、他の要因で説明できそうです。
でも、部屋の数は、犯罪率や大きな道路へのアクセスとか、関係あるとは思えません。予想通り、精度はあがりませんでした。
まぁ仕方がないので、このまま遺伝的アルゴリズム（GA）にかけようと思います。

NSGAⅡでの価格と部屋のパレート最適解

ここまできて、ふと思った。部屋の数とか整数のはずなのに、dtypesで型を見てみるとfloat型になっている。
とりあえず今回は全部float型で計算をした。
アルゴリズムはNSGAⅡで、遺伝子は１０００個、１００００回繰り返して計算した。

#ライブラリーのインポート
from platypus import NSGAII, Problem, Real

#遺伝子による変数の設定と学習モデル（目的関数）による予測結果

def objective(vars):
    crim = vars[0]
    zn   = vars[1]
    indus = vars[2]
    chas  = vars[3]
    nox  = vars[4]
    age  = vars[5]
    dis  = vars[6]
    rad  = vars[7]
    tax  = vars[8]
    ptratio  = vars[9]
    bb  = vars[10]
    lstat  = vars[11]

    ga_data = [[crim,zn,indus,chas,nox,age,dis,rad,tax,ptratio,bb,lstat]]

    #遺伝子を標準化
    ga_data_std = scaler.transform(ga_data)
    #学習モデルによる予測
    ga_MEDV_np = optimised_svr_MEDV.predict(ga_data_std)
    ga_RM_np = optimised_svr_RM.predict(ga_data_std)

    ga_MEDV = ga_MEDV_np[0]
    ga_RM = ga_RM_np[0]

    return [ga_MEDV,ga_RM] 

#Platypus遺伝的アルゴリズムの設定
#変数が１２個、目的関数が２個
problem = Problem(12, 2)

#変数の設定。すべて実数Realとした。また、最大値最小値は、元データの最大値と最小値にした。
problem.types[:] = [Real(0.006, 88.976),Real(0.00, 100.00),Real(0.740, 27.740),Real(0.00,1.00),Real(0.385, 0.871),Real(2.9000, 100.00),
                    Real(1.137, 12.126),Real(1.00, 24.00),Real(187.00, 711.00),Real(12.60, 22.00),Real(0.32, 396.90),
                    Real(1.980, 37.90)]

problem.function = objective
#価格を最小に、部屋数を最大にする。
problem.directions[:] = [Problem.MINIMIZE, Problem.MAXIMIZE]

algorithm = NSGAII(problem,population_size=1000)
algorithm.run(10000)

計算結果は、後で処理しやすいようにpandasで保存します。空のデータフレームを作って、そこにalgorithm.resultのvariablesの中に変数（遺伝子）が、objectivesに結果が入っていますので、それを空のpandasに取り込んでから、保存しました。

ga_result = pd.DataFrame(columns=['CRIM','ZN','INDUS','CHAS','NOX','AGE','DIS','RAD','TAX','PTRATIO','B','LSTAT','MEDV','RM'])
for i in range(len(algorithm.result)):
    ga_result.loc[i] = algorithm.result[i].variables[:] + algorithm.result[i].objectives[:]
ga_result.to_csv("ga_result.csv")

部屋の数と家の価格のパレート最適解は

青いプロットが元のデータ。赤のプロットが遺伝的アルゴリズムが計算した結果です。
これを見るとほかの条件を気にしなければ、部屋数が６で価格が１０くらいの家が、部屋が多くて家族が多くなっても暮らしやすくて価格が安いということになりますね。

元のデータだと、単純な相関のように見えるけど、遺伝的アルゴリズム的には、変曲点が見えているようです。
草薙素子チックにつぶやくと。。。
「そうささやくのよ、わたしの遺伝子が。。。」

３カ月くらいplatypusいじってます。マニュアルがあっさり書いてあるので、初心者の私にはだいぶ難しいのですが、とりあえずここまでできました。

【注意】本記事は適宜更新します．

概要

　衛星画像には30cmの高分解能のものから，Copernicusで無料で提供されている10mの中分解能のものと，解像度がそれぞれ異なります．高い解像度のものを使いたいのですが，高価なため入手するのは難しいです．でも，目的によってはそれほどの解像度は不要なのではと考え，車を検出することをターゲット として実験しました．

(credit: Cars Overhead With Context (COWC))

　ここで用いたコードは近々整理し，Githubにてシェアします．

1. 衛星画像の分解能と物体検知のレベルについて．

　衛星画像でどこまで見えるのか？例えば，報告されているものに以下があります．

(credit: FAS, Intelligence Resource Program)
　ここでは，10cm分解能であれば一つ一つの車がわかり，25cm分解能で車の区別は可能，50cm分解能では車種まではわからず，1m分解能になると車とはわかるけども，車種までわからない，と報告されています．
　これは20年以上前の評価であり，このときの”見える，見えない”は，人が”見える，見えない”が判断基準であり，人によっては見え方も異なり，統一されていません．また，現在であれば衛星画像から人が判読するよりも，機械に自動で判読させるほうが，品質および処理速度を考えると望ましい（当たり前？）であると思います．
　では，どこまで見えるのだろう？，ということで実験してみます．

1.1　衛星画像の分解能について

　衛星画像の分解能については，以下の宙畑さんの記事がよくまとめられていますので，ご参考にしてください．
人工衛星から人は見える？~衛星別、地上分解能・地方時まとめ~

(クレジット：宙畑）
　この図から，一番高い分解能が30cm（WorldView-3,4)であることがわかります．縦軸は地方時といって，観測時刻と思ってください．（この図のみ横軸のスケールが記載されていたので用いました．）
　また，高い分解能の衛星画像は非常に高価になります．
　
　(クレジット：宙畑）
　高い分解能の画像ほどきれいな画像なのでほしいのですが，目的によってはそこまでの分解能の画像は必要ないかもしれません．価格を低く抑えられるのであれば，より広範囲や異なる時期の多数の画像を入手できますね．

1.2　航空写真による車検知について．

　下記の記事で紹介しましたが，航空写真より車の検出精度を競うプログラムに， Cars Overhead With Context (COWC)があります．
PyTorchによる人工衛星画像から車の推定分布地図を作成してみる．
　ここでは，4カ国・6都市の航空写真（15cm分解能）の画像が提供されており，それぞれ一般自動車（トラックや重機は含まれません）の中心位置と間違うかもしれない物体の位置のアノテーションデータがセットで提供されています．
　今回の評価ではこのデータ・セットを用いました．まずは，15cmの分解能の画像より得られた車の検出率を基準として，画像をResizeすることで擬似的に分解能を落とし，その検出精度がどのようになるのか調べます．
　こちらにこのデータセットを用いた論文が紹介されていますが，極めて高い精度で検出できることが報告されています．興味のある方は，ご参考にしてください．

1.3　物体検知モデル

　物体検知（Object Detection)としては，有名なSSDやその他各種の手法があります．今回は，最近リリースされたYOLO v5 を用いました．
　Yolo(You Only Look Once)は有名なので多くの方がご存知かと思いますが，画像や動画より物体の位置と種類を検出する機械学習アルゴリズムです．そのVersion5が今年の春にリリースされました．（過去のYoloの開発者と異なるため，これをv5と言ってもよいのか，議論されているみたいです．）
　

(credit:YOLO v5)

　このアルゴリズムの精度評価が掲載されています．

　動画での高い精度の適用が考えられており，他のモデル（Efficient Det)と比べて高い精度で速く検知できるかが特徴になります．
　今回は，航空写真や衛星画像は静止画像ですので，一番精度の高いYOLOv5xを用いることとしました．

2.異なる解像度の航空写真・衛星画像から車を検出してみる．

2.1　環境構築（YOLOv5を準備する．）．

　YOLOv5のサイトより以下のコマンドでモデルをダウンロードします．（私はYOLOv5の仮想環境を構築し，その環境で使用しています．）

git clone https://github.com/ultralytics/yolov5.git

次に，yolov5を動作させるために，必要なモジュールを下記のコマンドより準備します．

pip install -r requirements.txt

　こちらのコマンドにより，以下のモジュールがインストールされます．

Cython
matplotlib>=3.2.2
numpy>=1.18.5
opencv-python>=4.1.2
Pillow
PyYAML>=5.3
scipy>=1.4.1
tensorboard>=2.2
torch>=1.6.0
torchvision>=0.7.0
tqdm>=4.41.0

　これでYolov5が利用できる環境が整いました．ここからは，動作を確認するために，jupyter labで実行します．
　次に，今回は画像分類で有名なcocoのパラメータをベースとしたFine tuningを行います．そのため，必要なWeightsをダウンロードします．以下のサイトから必要なWeightsをダウンロードし，YOLOv5に置きます．

YOLOv5のWeightsのダウンロードサイト

　これでYolov5を用いる環境が構築できましたので，次に学習データである航空写真データを取得します．

2.2　COWCより学習データの取得，およびデータの前処理．

　1.1項で述べたように，航空写真およびそのアノテーションデータ（車と非車の位置情報）をCOWCのサイトよりダウロードします．サイトの下部にダウンロード方法がいくつか紹介されていますが，私はhttpよりファイル（cowc-everything.txz）を選びダウンロードしました．
　ダウンロードしたファイルは圧縮されていますので，解凍して以下のファイルを確認します．

　ディレクトリ名は都市名になっており，そこに航空写真(png)とその画像のアノテーションデータ（車(annotated_cars)と非車(annotated_negatives)のpng)があります．
　航空写真は学習時のimage sizeに分割し，アノテーション画像からはYOLOv5に用いるための位置情報のtxtを作成します．その手順を紹介します．
　はじめに使用するモジュールを宣言します．

import os
import shutil
import numpy as np
from PIL import Image

import matplotlib.pyplot as plt
import numpy as np

import cv2

import pandas as pd

Image.MAX_IMAGE_PIXELS = 1000000000

　次に，都市名のディレクトリ構成を確認します．

dir_path = "../cowc/"
dirs =os.listdir(dir_path)
dirs

　出力が以下となります．

['Toronto_ISPRS',
 'Potsdam_ISPRS',
 'Vaihingen_ISPRS',
 'Utah_AGRC',
 'Columbus_CSUAV_AFRL',
 'Selwyn_LINZ',]

　'Vaihingen_ISPRS'と 'Columbus_CSUAV_AFRL'はモノクロ画像でしたので，今回の対象から省きました．
　では，　Toronto_ISPRSのファイルを確認し，航空写真の確認，リサイズによる衛星画像のシミュレーション，アノテーション画像よりYOLOv5を実行するための位置情報のtxtファイルを作成します．

dirs_image = dirs[0]
image_path = dir_path + dirs_image
image_files =os.listdir(image_path)
image_files

こちらを実行することで，ディレクトリ内のファイル情報を取得します．出力結果がこちらとなります．

['03553.png',
 '03553_Annotated.xcf',
 '03553_Annotated_Cars.png',
 '03559_Annotated_Cars.png',
 '03559_Annotated_Negatives.png',
 '03559.png',
 '03747_Annotated_Negatives.png',
 'README.ISPRS.txt',
 '03747.png',
 '03747_Annotated.xcf',
 '03559_Annotated.xcf',
 '03553_Annotated_Negatives.png',
 '03747_Annotated_Cars.png']

　これより，都市の航空写真とそのファイルと同じ名前と車の位置情報画像（Annotated_Cars)と非車の位置情報画像（Annotated_Negatives)で構成されているのがわかります．まず，航空写真のファイルのみを選択します．

#RGB画像のみのファイルリストを作成する．

file_name =[]

for i in image_files:
    if  'Annotated'in i:![Screenshot from 2020-11-28 20-54-02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/37ac50d8-cecd-90a4-eff3-0c8b5f037820.png)

        #print(i)
        pass
    else:
        if 'png' in i:
            file_name.append(i) 
        else:
            pass

file_name

上記を実行すると，以下のファイルであることがわかります．

['03553.png', '03559.png', '03747.png']

では，03553.pngの航空写真を確認します．

file_image_name = file_name[0]

#opencvで画像のよみこみ

im = cv2.imread(dir_path + dirs_image +'/'+ file_image_name)
print(im.shape)
print(file_image_name)
im_rgb = cv2.cvtColor(im, cv2.COLOR_BGR2RGB)

plt.imshow(im_rgb)

カナダの都市トロントの都市域の航空写真です．一部を拡大してみます．

plt.imshow(im_rgb[0:500,100:1500])

　車が多数写っているのがわかります．では，この画像を512x512の画像サイズに分割します．このとき，航空写真から衛星画像の分解能の異なる模擬画像を作成するため分解能を設定します．ここでは，現在の画像サイズの1/2や1/4とすることで，その解像度を低下させます．

#分解能に合わせて，リサイズの分母を決める．
#分解能30cmだと２，60cmで４，90cmで６となる．

rs = 1 #1, 2, 4, 6

im=Image.open(dir_path + dirs_image +'/'+ file_image_name)
img_resize = im.resize((im.size[0]//rs, im.size[1]//rs)) #ここで解像度をリサイズして変更している．
print(file_image_name)

file_name = os.path.splitext(os.path.basename(file_image_name))[0]
print(file_name)

DIR_OUTPUTS = '../cowc/crop_image'

height = 512
width = 512
img_size_v = img_resize.size[1]
img_size_h = img_resize.size[0]

#画像の分割処理関数
def ImgSplit(im):
    # 読み込んだ画像を２５６*２５６のサイズで１４４枚に分割する

    buff = []
    # 縦の分割枚数
    for v1 in range(int(img_size_v/height)):
        # 横の分割枚数
        for h1 in range(int(img_size_h/width)):
            h2 = h1 * width
            v2 = v1 * height
            #print(w2, h2, width + w2, height + h2)
            c = im.crop((h2, v2, width + h2, height + v2))
            buff.append(c)
    return buff


#画像の分割処理の実行
hi=0
for ig in ImgSplit(img_resize):
    hi=hi+1
    #print(hi)
    # 保存先フォルダの指定
    ig_jpeg = ig.convert('RGB')
    ig_jpeg.save(DIR_OUTPUTS + "/" + str(file_name) + '_' + str(hi) + ".jpg") #yolo v5にはpngはNG?

　このとき，2行目のresize関数にて画像の分解能をシミュレートします．後半でResizeした画像を512で分割し，DIR_OUTPUTS に保存します．
　次に，アノテーション画像を同様の処理を行います．

im=Image.open(dir_path + dirs_image +'/'+ file_name + '_Annotated_Cars.png')
img_resize = im.resize((im.size[0]//rs, im.size[1]//rs)) #ここで解像度をリサイズして変更している．
print(file_name + '_Annotated_Cars.png')

print(file_name)

DIR_OUTPUTS_p = '../cowc/crop_annotated'

height = 512
width = 512
img_size_v = img_resize.size[1]
img_size_h = img_resize.size[0]

#画像の分割処理関数
def ImgSplit(im):
    # 読み込んだ画像を２５６*２５６のサイズで１４４枚に分割する

    buff = []
    # 縦の分割枚数
    for v1 in range(int(img_size_v/height)):
        # 横の分割枚数
        for h1 in range(int(img_size_h/width)):
            h2 = h1 * width
            v2 = v1 * height
            c = im.crop((h2, v2, width + h2, height + v2))
            buff.append(c)
    return buff

#画像の分割処理の実行
hi=0
for ig in ImgSplit(img_resize):
    hi=hi+1
    # 保存先フォルダの指定
    ig.save(DIR_OUTPUTS_p + "/" + str(file_name) + '_' + str(hi) + ".png")

　この処理を行うことによって，車のアノテーション画像がリサイズされ512サイズに分割されます．同じく，非車のアノテーション画像も同様の処理を行います．

im=Image.open(dir_path + dirs_image +'/'+ file_name + '_Annotated_Negatives.png')
img_resize = im.resize((im.size[0]//rs, im.size[1]//rs)) #ここで解像度をリサイズして変更している．
print(file_name + '_Annotated_Negatives.png')

print(file_name)

DIR_OUTPUTS_n = '../cowc/crop_annotated_negatives'

height = 512
width = 512
img_size_v = img_resize.size[1]
img_size_h = img_resize.size[0]

#画像の分割処理関数
def ImgSplit(im):
    # 読み込んだ画像を２５６*２５６のサイズで１４４枚に分割する

    buff = []
    # 縦の分割枚数
    for v1 in range(int(img_size_v/height)):
        # 横の分割枚数
        for h1 in range(int(img_size_h/width)):
            h2 = h1 * width
            v2 = v1 * height
            c = im.crop((h2, v2, width + h2, height + v2))
            buff.append(c)
    return buff

#画像の分割処理の実行
hi=0
for ig in ImgSplit(img_resize):
    hi=hi+1
    # 保存先フォルダの指定
    ig.save(DIR_OUTPUTS_n + "/" + str(file_name) + '_' + str(hi) + ".png")

　アノテーション画像は分割後にそれぞれ，'../cowc/crop_annotated'と'../cowc/crop_annotated_negatives'に保存されます．
　次に，分割したアノテーション画像から，車および非車位置のtxtデータを作成します．このtxtデータは，Yoloのフォーマットにあわせて作成します．
　YOLOのアノテーション情報（txt）は以下のフォーマットとなります．

# class x_center y_center width height #

　Classは車であれば0,非車は1に，x-center, y-centerは対象の中心位置，widthとheightは対象物体の周囲枠（Box)のサイズとなります．
　このサイズは，ここでは20ピクセル（300cm）としました．車の縦幅より小さいですが，狭い間隔で駐車している車も識別したいため，まずはこの値としました．解像度を変えた場合，このピクセル値も300cmにあわせてみましたが，それよりも20ピクセルのまま（解像度によってBox値を変えた）場合のほうがPrecisionが高かったため，20ピクセルで固定しました．ここは，解像度によってパラメータを最適化する部分ですね．
　今回は，512サイズで統一しましたので，それぞれを512で割って規格化しました．
　では，アノテーション画像よりYOLOに用いるアノテーションファイルを作成します．

DIR_OUTPUTS_labels = '../cowc/crop_label'

for i in range(len(os.listdir(DIR_OUTPUTS))):
    crop_file_name_jpg = os.path.splitext(os.path.basename(crop_file_name[i]))[0]
    im = cv2.imread(DIR_OUTPUTS_p + '/' + crop_file_name_jpg +'.png')
    im_rgb = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
    row, cow = np.where(im_rgb==76)
    arr = []

    for ii in range(len(row)):
        arr.append([row[ii],cow[ii]])
    arr = np.array(arr)
    if not len(arr) == 0:
        array_pd = pd.DataFrame(arr)
        array_ppd = array_pd.copy()
        array_pd[0] = array_ppd[1]/512 
        array_pd[1] = array_ppd[0]/512
        array_pd[2] = 20/512
        array_pd[3] = 20/512
        array_pd.insert(0,'car',0)
    else:
        array_pd =  pd.DataFrame()

    im = cv2.imread(DIR_OUTPUTS_n + '/'  + crop_file_name_jpg +'.png')
    im_rgb = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
    row, cow = np.where(im_rgb==29)
    arr = []
    for ii in range(len(row)):
        arr.append([row[ii],cow[ii]])
    arr = np.array(arr)

    if not len(arr) == 0:
        array_nd = pd.DataFrame(arr)
        array_nnd = array_nd.copy()
        array_nd[0] = array_nnd[1]/512
        array_nd[1] = array_nnd[0]/512
        array_nd[2] = 20/512
        array_nd[3] = 20/512
        array_nd.insert(0,'car',1)
    else:
        array_nd =  pd.DataFrame()

    array_all = []
    array_all = pd.concat([array_pd, array_nd])

    text_name = crop_file_name_jpg
    array_all.to_csv(DIR_OUTPUTS_labels + '/' + text_name + '.txt', sep=' ', index=False, header=False)

アノテーションファイルは以下に'../cowc/crop_label'保存されます．
例えば，以下となります．

0 0.966796875 0.248046875 0.0390625 0.0390625
0 0.896484375 0.26953125 0.0390625 0.0390625
0 0.82421875 0.291015625 0.0390625 0.0390625
0 0.490234375 0.32421875 0.0390625 0.0390625
0 0.55859375 0.375 0.0390625 0.0390625
0 0.5625 0.375 0.0390625 0.0390625
0 0.470703125 0.408203125 0.0390625 0.0390625
0 0.3828125 0.43359375 0.0390625 0.0390625
0 0.037109375 0.45703125 0.0390625 0.0390625
0 0.29296875 0.4609375 0.0390625 0.0390625
0 0.0625 0.53515625 0.0390625 0.0390625
1 0.83984375 0.541015625 0.0390625 0.0390625
1 0.873046875 0.625 0.0390625 0.0390625
1 0.763671875 0.630859375 0.0390625 0.0390625
1 0.314453125 0.6484375 0.0390625 0.0390625
1 0.65625 0.666015625 0.0390625 0.0390625
1 0.75390625 0.8125 0.0390625 0.0390625
1 0.41796875 0.986328125 0.0390625 0.0390625

　この処理をそれぞれの都市で行い，Train用およびValidation用の画像とアノテーションファイルを準備します．このとき，Train用とValidation用のデータは分け，その割合は約7：3となります．

2.3　航空写真によるモデルの構築と検証．

　YOLOv5のモデルは，githubのチュートリアルを参考に実行します．
はじめに必要なモジュールをインポートします．

import torch
from IPython.display import Image, clear_output  # to display images

clear_output()
print('Setup complete. Using torch %s %s' % (torch.__version__, torch.cuda.get_device_properties(0) if torch.cuda.is_available() else 'CPU'))

　次に，学習結果を可視化するためにTensorboardを実行します．

# Start tensorboard
%load_ext tensorboard
%tensorboard --logdir runs

ここでTensorboardを実行状態すると，学習後にRecallやPrecisionなどの学習結果のグラフ化や異なる学習の比較を行うのに便利です．
　そして，次のコマンドで学習を始めます．

%cd ../yolov5
!python train.py --img 512 --batch 5 --epochs 100 --data ../cowc/data_resize_30cm.yaml --cfg models/yolov5x.yaml --nosave --cache --name cowc_test_15cm_ep100 --weights yolov5x.pt

--img : 学習画像のサイズ
--bath　：バッチサイズ
--epochs : エポック数
--data :train およびvalidationデータの場所とクラスの設定
--cfg :用いるモデルの設定
--name :学習名
--weights : fine tuningを行うときは，ここで設定します．

　詳細は，train.py --helpを実行すると，各パラメータの説明があるのでご参照ください．ここでは，データを示data.yamlについて説明します．

train: ../cowc/data_train
val: ../cowc/data_val

nc: 2
names: ['car', 'not_car']

　上記の通り，trainとvalのデータのアドレス，およびクラス数とその名前を設定しています．
　trainとvalディレクトリにはimagesとlabelsのディレクトリがあり，imagesに画像ファイルが，labelsにアノテーションファイルが入っています．この画像とアノテーションファイルは同じファイル名（拡張子は異なる）である必要があります．ご注意ください．
　Test用の画像およびモデルにて車検出後の画像の合成を行います．

dir_path = "../cowc/"
dirs =os.listdir(dir_path)
dirs_image = dirs[0]　#'Toronto_ISPRS'

image_path = dir_path + dirs_image
image_files =os.listdir(image_path)


#RGB画像のみのファイルリストを作成する．
file_name =[]
for i in image_files:
    if  'Annotated'in i:
        #print(i)
        pass
    else:
        if 'png' in i:
            file_name.append(i) 
        else:
            pass

file_image_name = file_name[2] #'03747.png'

まずは，Train用の画像準備と同じく対象となる画像を選択します．では，画像を確認します．

#opencvで画像のよみこみ

im = cv2.imread(dir_path + dirs_image +'/'+ file_image_name)
print(file_image_name)
im_rgb = cv2.cvtColor(im, cv2.COLOR_BGR2RGB)
#Image.fromarray(im_rgb).save('test.jpg')

plt.imshow(im_rgb)

　それでは，この画像を512サイズに分割してTest用の画像を準備します．方法はTrainと同じです．

DIR_OUTPUTS = '../cowc/test_512/images_resize/'

im=Image.open(dir_path + dirs_image +'/'+ file_image_name)
img_resize = im.resize((im.size[0]//1, im.size[1]//1)) 
print(file_image_name)

file_name = os.path.splitext(os.path.basename(file_image_name))[0]
print(file_name)

height = 512
width = 512
img_size_v = img_resize.size[1]
img_size_h = img_resize.size[0]

#画像の分割処理関数
def ImgSplit(im):
    buff = []
    # 縦の分割枚数
    for v1 in range(int(img_size_v/height)):
        # 横の分割枚数
        for h1 in range(int(img_size_h/width)):
            h2 = h1 * width
            v2 = v1 * height
            c = im.crop((h2, v2, width + h2, height + v2))
            buff.append(c)
    return buff

#画像の分割処理の実行
hi=0
for ig in ImgSplit(img_resize):
    hi=hi+1
    # 保存先フォルダの指定
    ig_jpeg = ig.convert('RGB')
    ig_jpeg.save(DIR_OUTPUTS + "/" + str(file_name) + '_' + str(hi) + ".jpg") 

　その後，分割されたファイル数をx軸(width)とy軸(height）の数を確認します．

img_size_v = int(img_size_v/width)
img_size_h = int(img_size_h/height)

　出力が以下となります．

img_size_v * img_size_h #308
img_size_v #14
img_size_h #22

念の為，ファイル数も確認しておきます．

image_list = os.listdir(DIR_INPUTS)
image_list_sort=natsorted(image_list)
len(image_list_sort)

出力

308

　これでTest用の画像が準備できました．次に，学習したモデルを用いてテスト画像から車を検出させます．ここでは，Trainに用いた画像ではなくValidationで準備した画像を用いました．

!python detect.py --weights runs/train/cowc_test_15cm_ep100/weights/last.pt --img 512 --conf 0.4 --source ../cowc/test_512/images_resize/  --save-txt --class 0 --name cowc_15cm

--weights: train.pyの実行時に作成したWeightsを呼び出します．ファイル名は，train.pyを実行結果の下部に示されますのでご確認ください．
--img : テスト画像のファイルサイズ
--conf : テスト結果で車の有無を判断するためのしきい値
--source : テスト画像があるディレクトリのアドレス
--save-txt : 検出結果をtxtでも保存
-- class : 識別するクラスの指定．ここでは”0”としたため，車のみを検出します．
--name :検出結果の画像を保存するディレクトリ

　解析処理後，データを確認します．このとき，yolov5のdetect.pyをCloneしたまま用いると，検出された車のBoxやその確信度（信頼度，Confidence)が大きく表示されるため，車が密集しているところではその違いが分かり難いです．そのため，utilsディレクトリにあるplot.pyのファイルの解析結果を表示するBoxのライン幅や信頼度のfontサイズを変更しました．コードで提供されているため，自分の利用方法にあわせて変更されることをお薦めします．
　では，処理されたTest画像を確認し，合成します．

height = 512 #512
width = 512 #512

DIR_INPUTS = '../yolov5/runs/detect/cowc_15cmresize/'#detect のディレクトリを指定．
DIR_OUTPUTS = '../cowc/test_512/build_picture/' #合成後の画像ファイルのアドレス
image_list = os.listdir(DIR_INPUTS)
image_list_sort=natsorted(image_list)

max_height = width * img_size_v
total_width = height * img_size_h

new_im = Image.new("RGB", (int(total_width), int(max_height)))

h = 0
w = 0

for i in range(img_size_v*img_size_h):
    new_im.paste(Image.open(DIR_INPUTS + "/"+ image_list_sort[i]), (w*width,h*height))
    w = w +1
    if img_size_h - w  == 0:
        w = 0
        h = h +1

new_im.save(DIR_OUTPUTS + 'test_15cm_resize.jpg')

　画像の合成後にファイルを確認します．

img = Image.open(DIR_OUTPUTS + 'test_90cm_resize.jpg')
#画像をarrayに変換
im_list = np.asarray(img)
#貼り付け
plt.imshow(im_list)
#表示
plt.show()

これで，航空写真（15cm分解能）のモデルの構築，およびそのモデルを用いた解析プロダクトが作成できました．画像を拡大して確認してみます．

　航空写真（15cm分解能）の解析により検出された車

　こちらの航空写真に写っているほぼすべての車が検知されており，その確信度は0.7以上と高い値となっています．
　次に，航空写真をResizeすることで分解能を模擬的に低下させ，同様の処理を行い，結果にどのような違いとなるのか，評価してみます．

2.4　異なる解像度の衛星画像と検知精度

　コードは2.3項で用いたものと同じになります．ただし，模擬する分解能に合わせてrsの値を変えます．例えば，分解能が半分の30cmにするのであればrsを2にします．

#分解能に合わせて，リサイズの分母を決める．
#分解能30cmだと２，60cmで４，90cmで６となる．

rs = 1 #1, 2, 4, 6

　TrainデータおよびValデータが作成したのち，前項とおなじくdata.yamlにて対象とするデータのディレクトリを指定し，train.pyを実行します．そして，実行にて作成されたそれぞれのweightを指定して，Test画像（比較のために，すべて同じカナダ・トロントの画像（Validationで使用））に解析処理を行い車を推定します．
　それでは，すべての処理画像を比較します．

　模擬衛星画像（30cm分解能）の解析により検出された車


　模擬衛星画像（60cm分解能）の解析により検出された車


　模擬衛星画像（90cm分解能）の解析により検出された車

　Test画像のそれぞれの解像度における車の検出結果を比較すると，15cm分解能と比べても，他の解像度でも多くは検出できているように見えます．それぞれの結果を比較してみます．

　
　こちらは，Tesorboardの結果の一部になります．ValidationデータのPrecisionをみると，解像度に変わらずほぼ同じ値になっています．一方，Recallをみると，解像度が15cm（濃い赤）が一番高いですが，ほかは大きく低下しています．これから，解像度に依存せず車と推定されたものは多くは正解となりますが，解像度が低下すると取りこぼしが多いことが予想されます．
　ただ，Trainデータのロス（Class loss)をみるとかなり低い値をとっており，解像度が90cmであっても学習量を100エポックより大きくすることで，より低下することが期待できます．
　Trainの結果は保存されているため，それより比較してみます．

run_dir = "../yolov5/runs/train/"
train_dirs =os.listdir(run_dir )
train_dirs.sort()
train_dirs

　上記を実行すると，それぞれの学習結果が保存されているのがわかります．

['.ipynb_checkpoints',
 'cowc_test_15cm_ep100',
 'cowc_test_30cmresize_ep100',
 'cowc_test_60cmresize_ep100',
 'cowc_test_90cmresize_ep100']

　では，15cmの解像度の画像の学習結果を見てみます．

df_15= pd.read_table(run_dir +train_dir + '/results.txt', header=None, delim_whitespace=True)

columns = ['Epoch','gpu_mem','Box','Objecness','Classsfication','total','targets','img_size','Precison','Recall','mAP@0.5','mAP@0.5:0.95','val Box','val Objectness','val Classification']

RS = '_15cm'
for i in range(len(columns)):
    column = columns[i] + RS
    columns[i] = column

df_15.columns = columns
df_15

出力．

　ここから，PrecisionとRecall,それぞれからf1scoreを求めてグラフにします．

df_15['f1score_15cm'] = 2 * df_15['Precison_15cm'] * df_15['Recall_15cm'] / (df_15['Precison_15cm'] + df_15['Recall_15cm'])

df_15[['Precison_15cm','Recall_15cm', 'f1score_15cm']].plot()

　同様の処理を30cm，60cm，90cmの学習結果よりそれぞれのパラメータを比較します．

#学習結果のグラフ化
fig = plt.figure(figsize=(14, 5))

ax1 = fig.add_subplot(1, 2, 1)
line1, = ax1.plot(df_15.index.values,df_15['Precison_15cm'],label='15cm') 
line2, = ax1.plot(df_30.index.values,df_30['Precison_30cm'],label='30cm')
line3, = ax1.plot(df_60.index.values,df_60['Precison_60cm'],label='60cm')
line4, = ax1.plot(df_90.index.values,df_90['Precison_90cm'],label='90cm')
ax1.set_ylim(0.0, 1.0) 
ax1.set_title("Precision")
ax1.set_xlabel('epoch')
ax1.set_ylabel('Precision')
ax1.legend(loc='lower right')

ax2 = fig.add_subplot(1, 2, 2)
line1, = ax2.plot(df_15.index.values,df_15['Recall_15cm'],label='15cm')
line2, = ax2.plot(df_30.index.values,df_30['Recall_30cm'],label='30cm') 
line2, = ax2.plot(df_30.index.values,df_60['Recall_60cm'],label='60cm') 
line2, = ax2.plot(df_30.index.values,df_90['Recall_90cm'],label='90cm') 
ax2.set_ylim(0.0, 1.0) 
ax2.set_title("Recall")
ax2.set_xlabel('epoch')
ax2.set_ylabel('Recall')
ax2.legend(loc='upper right')

plt.show()

　TensorboardのグラフをMatplotlibでグラフ化しただけですが，出力パラメータの加工や組み合わせによって，より詳細に解像度による違いを考察できそうです．

3．まとめ．

　解像度が15cmの航空写真をベースに，低解像度の画像をResize処理により模擬し，車の検出ができるのかどうか，Yolov5モデルを用いて評価しました．
　1m以下の解像度であれば，学習画像の作り方や，モデルのパラメータを最適化することで，高解像度と変わらない精度で検出できそうです．今回の結果をベースとして，色々試してみます．
　こちらの記事がみなさんの活動のご参考になれば幸いです．間違いやコメントなどありましたら，いただければ励みになりますのでお気軽にご連絡ください．

参考記事:

Deep Learning で航空写真から自動車をカウントする
PyTorchによる人工衛星画像から車の推定分布地図を作成してみる．
Cars Overhead With Context (COWC)
YOLO v5
人工衛星から人は見える？~衛星別、地上分解能・地方時まとめ~
FAS, Intelligence Resource Program