PyTorchでのモデルの中間層へのアクセス方法

  • PyTorchにはいろいろなモデルの書き方があるけど,このメモの内容を考えると学習パラメタが無いレイヤ含めてinitにもたせておくほうが良いのかもしれない。 (と思ったけど,一番最後のレイヤをとりあえずコピーしてforwardを再定義するやり方ならどっちでも良い,と思った)

重みにアクセスしたい場合

  • 重みを特定の値で初期化した場合などが利用としては考えられるか。
  • model.state_dict()が簡単。
  • 辞書で返してくれるのでmodelのinitに付けた名前でアクセスが可能。
model = Net()
model.state_dict() #=> {'conv1.weight':[...], 'conv1.bias':[...],...}

中間層の出力にアクセスしたい場合

  • forwardのhookしたり,色々やり方はあるようだけど,モデルの構築の段階からSequencialを使って分割して作成しておくのがわかりやすい。
  • 隠層の出力をそれぞれ得たい,という場合には後述のようにforwardを再定義するような方法があるようだ。
  • サンプルとしては,torchvisionのモデルが参考になる。例えば,AlexNetは特徴量抽出までのfeaturesと,それを使って分類するclassifeirに分けて実装されている。例えば,AlexNetは下記のようになっている。
class AlexNet(nn.Module):

    def __init__(self, num_classes=1000):
        super(AlexNet, self).__init__()
        self.features = nn.Sequential(
            nn.Conv2d(3, 64, kernel_size=11, stride=4, padding=2),
            nn.ReLU(inplace=True),
            nn.MaxPool2d(kernel_size=3, stride=2),
            nn.Conv2d(64, 192, kernel_size=5, padding=2),
            nn.ReLU(inplace=True),
            nn.MaxPool2d(kernel_size=3, stride=2),
            nn.Conv2d(192, 384, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
            nn.Conv2d(384, 256, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
            nn.Conv2d(256, 256, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
            nn.MaxPool2d(kernel_size=3, stride=2),
        )
        self.classifier = nn.Sequential(
            nn.Dropout(),
            nn.Linear(256 * 6 * 6, 4096),
            nn.ReLU(inplace=True),
            nn.Dropout(),
            nn.Linear(4096, 4096),
            nn.ReLU(inplace=True),
            nn.Linear(4096, num_classes),
        )

    def forward(self, x):
        x = self.features(x)
        x = x.view(x.size(0), 256 * 6 * 6)
        x = self.classifier(x)
    return x
  • これを用いて,学習済みのAlexNetの特徴量抽出までを使用したければ,下記のように使える。
model = models.alexnet(pretrained=True)
y = model.features(x)
  • VGGの各層からの出力を得るには下記のような例がある。ここでは,各層ではなく,特定のReLU層,3,8,15,22を出力している。やっていることは,VGG16のfeaturesのレイヤをそれぞれ通して行きながら,その過程で所望のレイヤの出力を保存している。
  • このやり方は結構汎用的だ。学習済みのモデルをとりあえずもってきて,自分でfowardを書き直しているイメージ。必要なら,適時dropoutを新たに挟む,とかも出来る。
import torch
import torch.nn as nn
from torchvision.models import vgg16
from collections import namedtuple

class Vgg16(torch.nn.Module):
    def __init__(self):
        super(Vgg16, self).__init__()
        features = list(vgg16(pretrained = True).features)[:23]
        self.features = nn.ModuleList(features).eval() 
        
    def forward(self, x):
        results = []
        for ii,model in enumerate(self.features):
            x = model(x)
            if ii in {3,8,15,22}:
                results.append(x)
        return results

CNN in MNIST with PyTorch (PyTorchの基本メモ)

  • PyTorchのチュートリアルとexampleはとても参考になる。0.4から色々変わっているようだし,改めて情報を整理する。
  • まず今回は,MNISTを動かしながら以下の項目についてメモしておく。
    1. MNISTデータのロード方法。可視化方法。
    2. ネットワーク,学習,テストの書き方。
    3. モデルの学習と保存

MNISTデータの利用方法

  • 生のMNISTはバイナリで扱いが面倒。なので,torchvisionのAPIを活用する。(tf.keras, chainerも同様のAPIがあるようだ。)
  • torchvisionはMNISTの他にも,Fasion-MNIST, EMNIST, Cifar-10, STL-10などいくつかのデータセットのダウンロードAPIを用意している。
  • 注意点として,MINST/Fashin-MNIST/EMNISTは同じクラスで実装されているので,ダウンロードで同じディレクトリを指定すると,既存のデータ(/processed/train.ptなど)がすでにあるのでダウンロードをしないことになる。よって,MNIST,EMNIST,Fashion-MNIST毎にrootディレクトリは変更する。transformはダウンロードしてtrain.ptなどを作る段階で適用するわけじゃないので,transformの違いは問題ない。

ダウンロード

  • 下記のコードを実行してMNISTをダウンロードする。rootはデータセットの格納ディレクトリ, downloadをTrueにするとWebから取得する。
  • trainはTrueにするとroot/processed/train.ptからデータを生成し,Falseの場合にはroot/processed/test.pyから生成する。
  • transformはデータに対して共通で適用する前処のの関数を渡す。関数はPILImageを1つ受け取って変換後の画像を返す関数。
from torchvision import datasets, transforms

def main():
    mnist_train = datasets.MNIST(
        root='/home/hoge/PyTorch_MLdata/EMNIST',
        download=True,
        train=True,
        transform=transforms.Compose([
            transforms.ToTensor(),
            transforms.Normalize((0.1307,), (0.3081,))
        ]))
    print(mnist)

    mnist_test = datasets.MNIST(
        root='/home/hoge/PyTorch_MLdata/MNIST',
        download=True,  # If already exist, download is ignored.
        train=False,
        transform=transforms.Compose([
            transforms.ToTensor(),
            transforms.Normalize((0.1307,), (0.3081,))
        ]))        

if __name__ == '__main__':
    main()        
  • 上記のコードを実行するとrootに指定したディレクトリに下記が生成されている。
- root_dir
  - processed
    - test.pt                  # テスト用データがPyTorchデータ向けに生成されている。
    - training.pt              # 画像用データがPyTorchデータ向けに生成されている。
  - raw
    - t10k-images-idx3-ubyte   # テスト用画像
    - t10k-labels-idx1-ubyte   # テスト用ラベル
    - train-images-idx3-ubyte  # 学習用データ画像
    - train-labels-idx1-ubyte  # 学習用データラベル

データローダー

  • torchにはデータローダというバッチ学習に便利なものがある。単純化以外にも読み込みのマルチスレッド化なども処理してくれるので活用するべき。
  • DataSetはデータの集合。集合ではるが,getitemなどが実装されている。DataLoaderはDataSetを内部に持っていて,バッチ単位で分割してくれるイテレータ
mnist_train_loader = torch.utils.data.DataLoader(mnist_train,
                                                 batch_size=32,
                                                 shuffle=True)
mnist_test_loader = torch.utils.data.DataLoader(mnist_test,
                                                batch_size=32,
                                                shuffle=True)
  • MNISTのようなデータ以外に,自分で入力集合Xsと出力集合Ysを生成した場合,それをTensorDatasetにして渡せばデータローダは簡単に使える。渡す際にはXs, Ysは0次元目のサイズが揃ったtorch.tensorであることが条件。0次元目を軸にgetitemされるからね。 (公式Documentを見た限りだと,別に渡すテンソルは2つに限定されず,n個渡せば,そのn個をローダで良きに取り出してくれるみたいだ。)
import torch.utils.data as Data
# x, y are input and label data (MUST be Pytorch Tensors)
torch_dataset = Data.TensorDataset(x, y)
loader = Data.DataLoader(
    dataset=torch_dataset,
    batch_size=32,
    shuffle=True,
    num_workers=2 # subprocess for loading
)
for epoch in range(n_epoch):
    for step, (batch_x, batch_y) in enumerate(loader):
        # training using batch_x, batch_y
        ...
  • ところで,データセットから直接サブセットを取得したい場合に,datasetはスライスオペレータをサポートしていないので少し工夫がいる。PyTorchのForcumに参考になるコードがある。
# データセットmnist_testから1000個取り出す場合
imgs, labels = zip(*[mnist_test[i] for i in range(1000)])

MNISTデータ可視化

  • あまり本論とは関係ないけれど,MNISTの各データ,ラベルの形式を理解するためにも可視化をしてみる。
  • datasetの返り値は(image:torch.Tensor, label:torch.Tensor)になっている。(というよりも,datasetの生成時にtransform.ToTensorを使ってそうしている)。
  • imageは(1,28,28)のテンソルなので,numpyに変換してmatplotlibのimshowに渡せば描画してくれる。ただし,imshowは(H,W), (H,W,3), (H,W,4)の形式で,各要素は0..1のfloatか0..255の整数でなければならないことに注意する。
image, label = mnist_train[0]  # Tensor(1,28,28), Tensor(1)
image_np = image.numpy()       
plt.imshow(image_np.reshape(28,28))
plt.show()

Network,学習,テストの書き方

  • exampleを見ると,ネットワークをクラスで,train/testはそれぞれ関数として実装している。
  • ネットワークの書き方は色々あるが,exampleでは学習パラメタを含むものをinitに書いて,純粋な関数(relu, poolingなど)はforward層で書いている。(nn.Dropout2dもパラメタ無いはずでは?と思うが。)
  • train, testの書き方は参考になる。deviceはGPU/CPUでコードを共通化する書き方。to(device)で必要だったらGPUテンソル化してくる。
  • train/testで重要な項目の1つが不要なところで勾配計算をさせないこと。テンソルの計算は基本的に計算グラフが構成され,計算のたびに値が保存されてメモリを消費する。それを避けるために,trainの中ではloss.itemでlossの値を抽出している。running_lossを計算するとき名などにテンソルのまま計算してしまうと無駄にメモリを消費するので避ける。また,testのようにそもそも勾配計算が不要な場合にはwith no_grad()でrequire_gradをFalseにしている。
  • また,PyTorchのクロスエントロピー周りは少し注意が必要。functional.cross_entropyはnll_lossとsoftmaxを組合せた関数。モデル側にsoftmaxを入れて,損失関数には単純に負のlog尤度関数(negative-log-likelihood)を使う場合にはnll_lossを使う。
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.conv1 = nn.Conv2d(1, 10, kernel_size=5)
        self.conv2 = nn.Conv2d(10, 20, kernel_size=5)
        self.conv2_drop = nn.Dropout2d()
        self.fc1 = nn.Linear(320, 50)
        self.fc2 = nn.Linear(50, 10)

    def forward(self, x):
        x = F.relu(F.max_pool2d(self.conv1(x), 2))
        x = F.relu(F.max_pool2d(self.conv2_drop(self.conv2(x)), 2))
        x = x.view(-1, 320)
        x = F.relu(self.fc1(x))
        x = F.dropout(x, training=self.training)
        x = self.fc2(x)
        return F.log_softmax(x, dim=1)

def train(args, model, device, train_loader, optimizer, epoch):
    model.train()
    for batch_idx, (data, target) in enumerate(train_loader):
        data, target = data.to(device), target.to(device)
        optimizer.zero_grad()
        output = model(data)
        loss = F.nll_loss(output, target)
        loss.backward()
        optimizer.step()
        if batch_idx % args.log_interval == 0:
            print('Train Epoch: {} [{}/{} ({:.0f}%)]\tLoss: {:.6f}'.format(
                epoch, batch_idx * len(data), len(train_loader.dataset),
                100. * batch_idx / len(train_loader), loss.item()))

def test(args, model, device, test_loader):
    model.eval()
    test_loss = 0
    correct = 0
    with torch.no_grad():
        for data, target in test_loader:
            data, target = data.to(device), target.to(device)
            output = model(data)
            # sum up batch loss
            test_loss += F.nll_loss(output, target, reduction='sum').item()
            # get the index of the max log-probability
            pred = output.max(1, keepdim=True)[1]
            correct += pred.eq(target.view_as(pred)).sum().item()

    test_loss /= len(test_loader.dataset)
    print('\nTest set: Average loss: {:.4f}, Accuracy: {}/{} ({:.0f}%)\n'.format(
        test_loss, correct, len(test_loader.dataset),
        100. * correct / len(test_loader.dataset)))

学習とモデルの保存

  • 後は実際に必要なエポック数だけ学習を回して,モデルを保存するだけ。
model = Net().to(device)
optimizer = optim.SGD(model.parameters(), lr=args.lr, momentum=args.momentum)

for epoch in range(1, args.epochs + 1):
    train(args, model, device, train_loader, optimizer, epoch)
    test(args, model, device, test_loader)

# パラメタだけを保存(推奨)
torch.save(model.state_dict(), 'model_param.pkl')
my_model.load_state_dict(torch.load('model_param.pkl'))

# モデル全体をpkl化
# ファイルが大きいこと,GPU/CPU情報を含む依存オブジェクトになってしまうことから非推奨
torch.save(model, 'model.pkl')  # save
my_model = torch.load('model.pkl') # load

クランプ関数(0〜255に値を丸める)のやり方

またまたcodewarsの他人の解答より。 本題はRGBの変換だったけど,その前処理で0以下の数値を0に,255以上の数値を255にクランプする処理がある。

僕は富豪的に下記で。

    lst = map(lambda x: 0 if x<0 else x, lst)
    lst = map(lambda x: 255 if x>255 else x, lst)

でも,一行で書くやり方があって,他の人はみんなこれを使ってた。よくやる処理で当たり前なのかもしれない。

clamp = lambda x: min(255, max(x, 0))

文字列からの数値抜き出しとソート時のインデックス取得

CodeWarsより。 わかりやすいように少し改題。また,CodeWarsはPython2で動いているので下記のコードは通らない。

問題
入力:数値を含む文字列の配列が与えられる。ただし,各文字列には1~9の数字が1つだけ入っており,文字列間に含まれる数値の重複はない。
   例) sentence = ['hoge8hoge', 'aa2bed', '5foobar']
出力:文字列に含まれる数値順にソートせよ。
   例) sentence' = ['aa2bed', '5foobar', 'hoge8hoge']

練習になった点は,(1)文字列から数値部を抜き出すには?,(2)ソートした時に,そのインデックスは取得できる? (2)に関して,numpyだとargsortで簡単にできる。sortした時にそのargも返してくれれば良いけど,それは無いので,両方欲しいときはargsortをもとのnumpy配列にインデックスとして渡す。 matlabとかはmax演算の時にargmaxもくれるけど,その違いに似ている。でも,numpyもCとかに比べれば相当楽だ。Pythonは便利だなぁ。

文字列の抜き出し = re.subを使う。味噌はsubは抽出はできずに代替(substitute)するので,数値以外の部分を削る。

import re

string = 'hoge3'
re.sub(re.compile('\D'), '', string)

# 配列に適用したいならmap
map(lambda x: re.sub(re.compile('\D'), '', x), sentence)
imoprt numpy as np

d = [1, 3, 5, 2]
sorted_index = np.argsort(d)  # [0, 3, 1, 2]
sorted = np.array(d)[sorted_index]  # もしソートしたリストも欲しいなら

配列の中に奇数回出現する1つの整数値

codewarsというので遊んでいたら,他人の良い回答を見つけた。

問題
 整数の配列が与えられる。その中に1つだけ奇数回出現する整数値が存在する。それを求めよ。
 例)[1,3,5,2,6,8,4,3,1,5,2,6,8,3] -> 3

色々とやり方はあると思うけれど,XORを使う方法がなるほどと思った。 対象の整数以外は全て偶数回出現する,というのを利用して,XORとればそれは0になる。 reduceはpython3からfunctoolsに入ったみたい。map/fileterは標準なのに変な感じする。

from operator import xor
from functools import reduce
seq = [1,3,5,2,6,8,4,3,1,5,2,6,8,3]
return reduce(xor, seq)

並列計算の効果検証用プログラム

並列化の試行やJITの効果テストなどのために,Forをたくさん回す計算例が必要になる。 そのような例としてモンテカルロ法を置いておく。

ToDo 並列計算のテスト用としては,モンテカルロ法のような単純に並列化される(データ共有が無い)ものだけでなく,適時データ共有が必要となる例を別途置いておく。

モンテカルロ法による円周率計算

  • モンテカルロ法の例題でよく見る計算。[0,1]の範囲の乱数を生成して,その点が半径1の円の中に入るか否か?(点が円の中か否かは,点の距離dがsqrt(x2+y2)で計算できるから,円周率piを知らなくても判定できる)をランダムに試行して,N回サンプリングすると,円の中に入る比率が1:pi/4になることから推定する(円の面積がpirrなのは知っている)
# monte_carlo.py
import random

def calc_pi(iter):
    cnt = 0
    for _ in range(iter):
        x, y = random.random(), random.random()
        if x*x + y*y < 1:
            cnt += 1
    return (cnt*4.0)/iter
  • 全く同じ書き方でC言語でも。
#include <stdio.h>
#include <stdlib.h>

double calc_pi(int iter)
{
  int cnt=0;
  double x, y;
  
  for (int i=0; i<iter; i++) {
    x = (double)rand()/RAND_MAX;
    y = (double)rand()/RAND_MAX;
    if ((x*x + y*y) < 1)
      cnt++;
  }
  return cnt*4.0/iter;
}

int main(int argc, char** argv) {
  int iter;
  double pi;
  
  iter = atoi(argv[1]);
  pi = calc_pi(iter);
  printf("pi:%f", pi);

  return 0;
}
  • ちなみに,生のPythonは20倍くらい遅い。
# Python
$ time python3 monte_carlo.py 10000000 # 10M
real 0m3.838s                          # 0.4us/iter

# C(Clang, -O3)
$ time ./a.out 10000000                # 10M
real 0m0.214s                          # 0.02us/iter

モンテカルロ法1 (モンテカルロ積分)

Rによるモンテカルロ法入門のメモ,その1。

イントロ

  • 統計的推論では最適化と積分の計算,特に実用上は数値計算が必要になる。
  • 数値積分の手法としてモンテカルロ積分がある。シンプソン法などの区分法は多次元の場合に区分を細かくしようとすると計算量が増加してしまうので適用が困難である。そのような多次元の場合に対してモンテカルロ積分は有用。

モンテカルロ積分

  • モンテカルロ積分の考え方は非常に単純で,下記の期待値計算の積分を対象の確率分布からのサンプリングによる平均値で近似する方法。
  • 例えば,f(x)が確率密度関数で,xがf(x)に従う時のh(x)の期待値を計算したいときに,

 E_f(h(x)) = \int_{x}f(x)dx

 \frac{1}{m}\sum_{i=1}^{m}h(x_j)

で近似する。

  • 同様に考えて,定積分も対象の範囲に対する一様乱数を生成して,その時のh(x)の値を求めて平均取れば,その区間の平均値となるので,積分区間を掛けて上げれば積分値が得られる。(b-a倍するのを忘れないように)
import numpy as np
from scipy import integrate
import math

f = lambda x: math.sin(x)/x  # 被積分関数

def mc_integrate(f, a, b, n=10000):
    def scaling(x):
        '''
        np.random.randで生成する[0,1)の一様乱数を[a,b)に変換する関数
        '''
        return (b-a)*x + a

    xs = np.random.rand(n)
    xs = np.vectorize(scaling)(xs)
    return (b-a) * np.vectorize(f)(xs).sum()/n


def main():
    v = mc_integrate(f, 0, 0.0, 2.0*math.pi)
    
    # Scipyを使って検算
    v_val = integrate.quad(f, 0.0, 2.0*math.pi)[0]
    print('Scipy, v={:.5}, MC, v={:.5}'.format(v_val, v))

if __name__ == '__main__':
    main()