[solafune] 衛星画像から空港利用者数予測ベースモデル

先日新しいデータ分析コンペサイトsolafuneがリリースされました。衛星データに特化したサイトということで、今後の展開が非常に楽しみです。

第1回目の題材が衛星画像から空港利用者数予測ということでベースモデルを作成してみました。GPU環境下でのColabにて上から順に実行することを想定しています。

また想定読者層としては、コンペを始めたいけれども、どのように手をつけたらいいかわからない方に向けて書いています。もし上級者がご覧になられた場合には有用な知見は皆無なので、気が向いたら修正箇所などのアドバイスをいただけますと幸いです。

(1) ファイルのデータ一式をinputディレクトリに入れています。Colabの場合には、Google Driveに入れて、毎回の実行時にローカルにコピーをすることをお勧めします。

!ls input
---
(出力)
evaluatemodel.zip	  testimage.zip		     trainimage.zip
testdataset_anotated.csv  traindataset_anotated.csv  uploadfile.csv

(2) 今回使用するライブラリをインストールします。（Colabでない場合には、torchは別途インストールする必要があります。PyTorch公式を参照ください。）

!pip install pytorch_pfn_extras
!pip install torch_optimizer

(3) ファイルを解凍します。

!cd input; unzip -q trainimage.zip
!cd input; unzip -q testimage.zip
!cd input; unzip -q evaluatemodel.zip

(4) ライブラリをインポートします。

import numpy as np
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt

from tqdm.auto import tqdm
from PIL import Image

import torch
import torch.nn as nn
import torch.nn.functional as F
from torchvision import datasets, transforms, models

import torch_optimizer as optim
import pytorch_pfn_extras as ppe
import pytorch_pfn_extras.training.extensions as extensions

(5) ラベルデータをロードします。

df_train = pd.read_csv("input/traindataset_anotated.csv", header=None)
df_test = pd.read_csv("input/testdataset_anotated.csv", header=None)

(6) データサイズが統一されていないので適当なサイズでクロップし、ラベルデータをロードして対数に変換しておきます。（再学習がしやすいように保存してあります。再学習する際には書いていないですがnp.loadで読み込んでください）

def center_crop(img, Lx, Ly):
    x, y, _ = img.shape
    startx = x//2 - (Lx//2)
    starty = y//2 - (Ly//2)    
    return img[startx:startx+Lx, starty:starty+Ly]
size = 1700

# 訓練データ
x_train = []
y_train = []
for i in tqdm(range(60)):
    img = np.array(Image.open(f'input/image/train_{i}.jpg'))
    x_train.append(center_crop(img, size, size))
    iy = np.log10(df_train[1][i])
    y_train.append(iy)
x_train = np.stack(x_train)
y_train = np.stack(y_train)
np.save("x_train.npy", x_train)
np.save("y_train.npy", y_train)

# 評価用データ
x_valid = []
y_valid = []
for i in tqdm(range(20)):
    img = np.array(Image.open(f'input/image/test_{i}.jpg'))
    x_valid.append(center_crop(img, size, size))
    iy = np.log10(df_test[1][i])
    y_valid.append(iy)
x_valid = np.stack(x_valid)
y_valid = np.stack(y_valid)
np.save("x_valid.npy", x_valid)
np.save("y_valid.npy", y_valid)

# テストデータ
x_test = []
for i in tqdm(range(58)):
    img = np.array(Image.open(f'input/image/comp_{i}.jpg'))
    x_test.append(center_crop(img, size, size))
x_test = np.stack(x_test)
np.save("x_test.npy", x_test)

(7) ベースモデルを定義します。ここではResNet18で試してみました。

class Net(nn.Module):
    def __init__(self):
        super().__init__()
        self.resnet = models.resnet18(pretrained=True)
        self.resnet.fc = nn.Sequential(
            nn.Linear(512, 1024),
            nn.ReLU(inplace=True),
            nn.Linear(1024, 1)
        )

    def forward(self, x):
        return self.resnet(x).squeeze(1)

(8) Augmentationを定義します。train時にはランダムに切り出し、反転、回転、明るさなどを変換し、validation時には変動しないように中心でクロップしています。albumentationsを使ってもっと多様に変動することもできます。albumentationsでできることに関してはこの記事がおすすめです、

transform_train = transforms.Compose([
                  transforms.ToPILImage(),
                  transforms.Resize(size=(720, 720)),
                  transforms.RandomCrop(size=(512, 512)),
                  transforms.RandomHorizontalFlip(p=0.5),
                  transforms.RandomVerticalFlip(p=0.5),
                  transforms.RandomRotation(15),
                  transforms.ColorJitter(brightness=0.5, contrast=0.5, saturation=0.5, hue=0.1),
                  transforms.ToTensor(),
              ])
transform_valid = transforms.Compose([
                  transforms.ToPILImage(),
                  transforms.Resize(size=(720, 720)),
                  transforms.CenterCrop(size=(512,512)),
                  transforms.ToTensor(),
              ])

(9) データセットの定義、生成します。

class Dataset(torch.utils.data.Dataset):
    def __init__(self, x_data, y_data=None, transform=None):
        self.x_data = x_data
        self.y_data = y_data
        self.transform = transform

    def __len__(self):
        return self.x_data.shape[0]

    def __getitem__(self, idx):
        ix = self.x_data[idx]

        if self.transform:
            ix = self.transform(ix)

        if self.y_data is not None:
            iy = self.y_data[idx]
            return ix, iy
        else:
            return ix

train_dataset = Dataset(x_train, y_train, transform_train)
valid_dataset = Dataset(x_valid, y_valid, transform_valid)

(10) データローダーの定義します。num_workersは実行環境のスレッド数、batch_sizeはメモリ数に応じて適当に変えてください。

train_loader = torch.utils.data.DataLoader(
    train_dataset,
    batch_size=30,
    num_workers=4,
    pin_memory=True,
    shuffle=True)
valid_loader = torch.utils.data.DataLoader(
    valid_dataset,
    batch_size=20,
    num_workers=4,
    pin_memory=True,
    shuffle=False)

(11) 訓練・評価ループを定義します。pytorch_pfn_extrasを用いて出力部分を簡略化しています。

def train(manager, model, loss_fn, scheduler, device, train_loader):
    while not manager.stop_trigger:
        model.train()
        for batch_idx, (data, target) in enumerate(train_loader):
            with manager.run_iteration(step_optimizers=['main']):
                data, target = data.to(device), target.to(device)
                output = model(data)
                loss = loss_fn(output, target.float())
                ppe.reporting.report({'train/loss': loss.item()})
                loss.backward()
        scheduler.step()

def validation(model, device, data, target):
    model.eval()
    with torch.no_grad():
        test_loss = 0
        data, target = data.to(device), target.to(device)
        output = model(data)
        test_loss += F.mse_loss(output, target.float()).item()
        ppe.reporting.report({'val/loss': test_loss})

(12) モデル、最適化、スケジューラを生成します。

n_epochs = 120
device = torch.device("cuda")
model = Net()
model.to(device)
loss_fn = torch.nn.MSELoss()
optimizer = optim.RAdam(
    model.parameters(),
    lr= 1e-3,
    betas=(0.9, 0.999),
    eps=1e-8,
    weight_decay=1.0e-3,
)
scheduler = torch.optim.lr_scheduler.OneCycleLR(
    optimizer,
    max_lr=5.0e-3,
    steps_per_epoch=1,
    epochs=n_epochs
)

(13) pytorch_pfn_extrasのextensionsを記載します。これが実行時のレポートやグラフの描画、snapshotの作成などを行ってくれます。

my_extensions = [
    extensions.LogReport(),
    #extensions.ProgressBar(),
    extensions.observe_lr(optimizer=optimizer),
    extensions.Evaluator(
        valid_loader, model,
        eval_func=lambda data, target:
            validation(model, device, data, target),
        progress_bar=False),
    extensions.PlotReport(
        ['train/loss', 'val/loss'], 'epoch', filename='loss.png'),
    extensions.PrintReport(['epoch', 'iteration', 'train/loss', 'val/loss', 'lr', 'elapsed_time']),
    extensions.snapshot(),
]

manager = ppe.training.ExtensionsManager(
    model, optimizer, n_epochs,
    extensions=my_extensions,
    iters_per_epoch=len(train_loader),
    stop_trigger=trigger)

(14) 学習を実行します。私の環境では120epochが663秒で完了しました。

train(manager, model, loss_fn, scheduler, device, train_loader)

(15) 最終モデルを保存します。

torch.save(model.state_dict(), "last.pth")

(16) 評価用データでの予測をグラフで確認し、予測が外れていないかを念のため確認します（任意）。

model.eval()
with torch.no_grad():
    for ix, iy in valid_loader:
        ix = ix.to(device)
        y_valid_pred = model(ix).cpu().detach().numpy()
        y_valid = iy.cpu().detach().numpy()

print(np.mean(np.square(y_valid - y_valid_pred)))
plt.plot(y_valid_pred, y_valid, "o")
plt.xlabel('Prediction')
plt.ylabel('True')
x = np.arange(5, 8)
plt.plot(x, x, "--", color="grey")
plt.grid()
plt.show()
---
（出力）
0.5196750537803555

予測が６付近に固まっていていまいちですね。

(17) 提出ファイルを作成します。ターゲットを常用対数としていたので10の肩に載せることを忘れないでください。

# 保存したモデルをロードする場合
# model.state_dict(torch.load("last.pth"))
x_test = np.load("x_test.npy")
test_dataset = Dataset(x_test, transform=transform_valid)
test_loader = torch.utils.data.DataLoader(
    test_dataset,
    batch_size=58,
    shuffle=False)
model.eval()
with torch.no_grad():
    for ix in test_loader:
        ix = ix.to(device)
        y_pred = model(ix)
    y_pred_np = y_pred.cpu().detach().numpy()
with open("submit_v1.csv", "w") as fp:
    for i in range(58):
        fp.write(f"comp_{i}.jpg,{np.power(10, y_pred_np[i])}\n")

(18) Colabでは再起動のたびにデータが消えてしまうので、再現性を確保するためのファイル一覧をGoogleDriveに保存します。（ColabでGoogle Driveに保存することを想定して書いています。パスなどは個々の環境に変えて書いてださい）

!mkdir "drive/My Drive/solafune/v1"
!cp last.pth "drive/My Drive/solafune/v1/"
!cp "submit_v1.csv" "drive/My Drive/solafune/v1/"
!cp "result/log" "drive/My Drive/solafune/v1/"
# snapshotファイルを保存したい場合
!cp "result/snapshot_iter_240" "drive/My Drive/solafune/v1/"

(19) サブミットしてスコアを確認します。

2020/10/24,14:00時点で7位でした。validationでのスコアが0.51だったので、結果にズレがありますね。testデータはvalidationデータよりも来場者の幅が広いのかもしれません。画像のクロップが適当だったので、取り込み位置などを変更してどう影響するかを分析してみたいですね。

以上です。

他に気力の関係で書ききれなかったのですが、neptune.aiなどで実験管理することもお勧めです。ここでは深層学習を使いましたがもちろんこれが唯一の回答ではないことは自明なので思いついた方法をいろいろ試していくのが良いでしょう。特に今回はデータ数が少ないので、画像から特徴量を取得するなどし深層学習以外の手法ですすめるのが良いかもしれません。コンテストが終わった際には手法などを公開・共有し、データ分析コンテスト界隈も盛り上げていただければ幸いです！

今後のコンテストに期待すること

最後にsolafuneの運営さんに個人的な期待が２つあります。

１つ目は、本当にやる意義のある題材を選んでほしいという点です。今回の題材でいうと統計情報として得られる空港利用者を衛星画像から求める必要があるのかが疑問です。これは正解を調べてそれにoverfitさせる不正にも繋がってしまいます。私がコンペに参加するのは普段扱えないデータを扱うことで得られる知見やそのモデルを作る意義を感じたときで、競うために必要のないモデルを作ることにモチベーションを持っていません。もちろんこれは人によって参加する理由は様々なので競技としてテーマを選ぶのも一つだとは思います。ただせっかく参加者が多くの時間や労力をかけて作るモデルを争って終わりにするのではなく、今後の衛星データの分析に活用できる題材にして欲しいと期待しています。

２つ目は、ルールはSIGNATEをベースに作られていますが、以下のルールが考えた末入れたのでなければ今後は外してほしいというものです。

学習データのラベルを書き換えてモデルを学習することは禁止します。ただし、学習用画像データやラベルデータを画像処理手法により自動で水増しすることや、加工して学習用データとして利用することは可能です。

例えば今回のコンテストでは空港をランク付けしてそのランクを予測することや、空港の位置座標を用いることや、山の中かどうか、人工島かどうか、都市があるかないか、滑走路の数などをラベルとして予測して汎化性能を上げるなどのアプローチが取れなくなります。確かにコンテストとして限られたアプローチの中でモデリングするというのも1つの競技性ではあるのですが、解法のバリエーションを減らしてまで導入したいかどうかを考えてからルールを入れてほしいと考えています。
（ただ途中からルールが変わるのはこれまでモデリングを行っていた人の作業が無駄になりかねないので次回から考えていただければ）

今後、面白そうなコンペをどんどん開いていっていただけることを期待しています！！

cf. @solafune (https://solafune.com)