2017年10月29日日曜日

CRISPR (その1)

クリスピーじゃなくてCRISPR-Cas9(クリスパーキャスナイン)という遺伝子編集技術をつかった様々な結果をこれまで目にしていたのですが、なにか良く分からないが遺伝子を編集できるんだな、という以上の知識はありませんでした。偶然、書店で平積みされていた「CRISPR 究極の遺伝子編集技術の発見(ISBN978-4-16-390738-3)」というタイトルのハードカバーの書籍を見つけ、CRISPR-Cas9を見つけた主要な科学者が著者だったので、購入。翻訳本なので、1年くらい前に元の書籍は発売されていた模様。

全体が物語風味で構成されていますが、専門家ではない人にも技術の内容が分かりやすい構成になっています。第3章までが基礎技術部分の解説で、CRISPR-Cas9がDNAを目的の位置で正確に切断する原理が説明されています。第4章は、CRISPR-Cas9の切断能力をベースにした遺伝子編集技術の解説、それ以降は、すでに実施された様々な応用例と、倫理的な問題についての議論が書かれています。

以下はこの書籍を読んでいる最中にとった専門外の人間のメモなので、そのまま信じないように。ググって見つけたものを真に受けてはいけません。

遺伝子関連の言葉

まず、
ゲノム∋遺伝子∋塩基対
という関係があります。また、
塩基={A, G, C, T}≒ヌクレオチド
塩基対={A=T, C=G}
ヌクレオチドの列=遺伝子=タンパク質=アミノ酸の列
ゲノムの大きさ=塩基対の長さ
とのこと。塩基に色々くっつくとヌクレオチドになるらしい。
遺伝子が翻訳されてタンパク質になるので、1対1に対応している。

遺伝子のサイズ

遺伝子の大体のサイズが記載されていたので、表にしてみました。

名前 遺伝子の数 塩基対の数
ウイルス 数個 数千個
細菌 約4,000個 数百万個
ハエ 1万4千個 数億個
ヒト・ネズミ・蛙 21,000個(*1) 32億個 (*2)
サラマンダー (*1)の10倍 (*2)の10倍
一部の植物 (*1)の100倍 (*2)の100倍
ミトコンドリア ? 16,000個
(*2)このうちの8%(2億5千万個)が人類の祖先がレトロウイルスに感染した名残

おおよそ、ひとつの遺伝子は数千個~1万数千個の塩基対で構成されているようです。アミノ酸が何個の塩基対で作られるかは、この書籍には記載がありませんでした。本題からはずれるのでなくても不思議ではありませんが。


さて、ここからCRISPRに行き着くまでが遠いのですが、遺伝子治療と呼ばれていた遺伝子を挿入する技術とCRISPR-Cas9を使った編集技術は全く異なり、少なくとも細菌の遺伝子を書き換えるときの正確性は極めて高い技術だということです。ただし、人間やサルやマウスの場合は、CRISPR-Cas9を個々の細胞にまで届けないといけないので、色々課題が残っているようです。

MNISTでのmixup法の効果

数日前にarXivに投稿されていたmixup法(https://arxiv.org/abs/1710.09412)をMNISTで試してみました。 手法自体は簡単で、ニューラルネットワークの学習時に使う入力と出力を2つずつランダムに選び、それを比率\(\lambda\)で混ぜるだけです。式で書くと、ランダムに選んだ入力と出力を\((x_i,y_i), (x_j,y_j)\)のように2つ選ぶとき、学習時に渡す入力\(\tilde{x}\)と出力\(\tilde{y}\)はそれぞれ、 \[ \tilde{x} = \lambda x_i + (1-\lambda) x_j \] \[ \tilde{y} = \lambda y_i + (1-\lambda) y_j \] となります。\(\lambda\)はサンプルごとにベータ分布からランダムに決めます。

今回はMNISTで試すので、数字が半透明で重なります。ベータ分布の\(\alpha=\beta=2\)で32個生成した結果が下図です。

\(\alpha=\beta\)の下で、\(\alpha\)を色々変えて学習した結果が下表です。学習は30エポックのみです。

\(\alpha\) Train acc Validation acc
n/a 0.999933 0.991600
0.1 0.999767 0.991800
0.2 0.999483 0.990600
0.3 0.999283 0.990300
0.4 0.998983 0.991600
0.5 0.998633 0.989500
1 0.997483 0.989700
2 0.996383 0.990000

1行目は通常の方法で混ぜずに学習させた場合の結果です。Train accの計算には混ぜていない学習用データを用いています。\(\alpha\)を増やすと半々で混ざるサンプルの確率が増え、Train accが徐々に下がっていきます。一方、Validation accは下がっているようにも見えますが、\(\alpha=0.1\)の場合は若干上がっています。差が小さすぎるので誤差かもしれません。

論文に書かれているような効果は今回のMNIST実験では見られませんでしたが、こういう混ぜ方をしても学習ができることは確認できました。

以下、今回の実験で使ったコードです。122行目のeh.alphaの値を変えることで\(\alpha\)を変更します。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
# -*- coding: utf-8 -*-
import sys,os,math
sys.path.append(os.path.dirname(os.path.abspath(__file__)) + '/../keras-examples/src')
import numpy as np
import keras
from keras.models import Sequential, Model, Input, model_from_json
from keras.layers import Dense, Activation, Reshape, Flatten, Dropout
from keras.layers.normalization import BatchNormalization
from keras.layers.convolutional import Conv2D
from keras.layers.advanced_activations import LeakyReLU
from keras.utils import to_categorical
from util.history import ExperimentHistory
from keras.datasets import mnist
from PIL import Image, ImageDraw, ImageFont

TRAIN_EPOCH = 30
GENERATED_IMAGE_PATH = 'generated_images/'
FONT_DIR = "C:/Windows/Fonts"

# Image size is the same as MNIST
IMAGE_WIDTH = 28
IMAGE_HEIGHT = 28

def make_model():
    model = Sequential()
    model.add(Conv2D(64, (5, 5), strides=(2, 2), padding='same',
                     input_shape=(IMAGE_WIDTH, IMAGE_HEIGHT, 1), 
                     data_format='channels_last'))
    model.add(LeakyReLU(0.2))
    model.add(Conv2D(128, (5, 5), strides=(2, 2), data_format='channels_last'))
    model.add(LeakyReLU(0.2))
    model.add(Flatten())
    model.add(Dense(256))
    model.add(LeakyReLU(0.2))
    model.add(Dropout(0.5))
    model.add(Dense(10))
    model.add(Activation('softmax'))
    print(model.summary())
    return model

def combine_images(generated_images, margin=20):
    total = len(generated_images)
    cols = int(math.sqrt(total))
    rows = math.ceil(float(total)/cols)
    width, height = generated_images[0].shape[0:2]
    print(generated_images[0].shape)
    combined_image = np.ones(((height+margin)*rows, width*cols), dtype=generated_images[0].dtype)*-0.5
    for index, image in enumerate(generated_images):
        i = int(index/cols)
        j = index % cols
        combined_image[(height+margin)*i:(height+margin)*i+height, width*j:width*(j+1)] = image[:, :, 0]
    return combined_image, cols, rows, width, height, margin

# Data generator
def data_gen(x, y, batch_size, alpha, num_classes):
    # Generate data eternally
    assert len(x) == len(y)
    epoch = 0
    index = 0
    font = ImageFont.truetype(FONT_DIR+"/{}".format("Tahoma.ttf"), 10) # For drawing numbers
    while True:
        index += 1
        data = []
        labels = []
        for i in range(batch_size):
            i1 = np.random.randint(0, len(x)) 
            i2 = np.random.randint(0, len(x)) 
            lam = np.random.beta(alpha, alpha)
            # Mixup two samples
            xm = lam * x[i1] + (1.0 - lam) * x[i2]
            ym = lam * y[i1] + (1.0 - lam) * y[i2]
            data.append(xm)
            labels.append(ym)

        ## Draw mixed-up images
        # image, cols, rows, width, height, margin = combine_images(data)
        # image = image*127.5 + 127.5
        # if not os.path.exists(GENERATED_IMAGE_PATH):
        #     os.mkdir(GENERATED_IMAGE_PATH)
        # img = Image.fromarray(image.astype(np.uint8))
        # draw = ImageDraw.Draw(img)
        # for i in range(len(data)):
        #     yi = int(i/cols)
        #     xi = i % cols
        # img.save(GENERATED_IMAGE_PATH+"%04d_%04d.png" % (epoch, index))

        yield np.array(data), np.array(labels).reshape((batch_size, num_classes))

def run(eh):
    # Prepare MNIST data
    (x_train, y_train), (x_val, y_val) = mnist.load_data()
    x_val = (x_val.astype(np.float32) - 127.5)/127.5
    x_val = x_val.reshape(x_val.shape[0], x_val.shape[1], x_val.shape[2], 1)
    y_val = keras.utils.to_categorical(y_val, num_classes=10)

    x_train = (x_train.astype(np.float32) - 127.5)/127.5
    x_train = x_train.reshape(x_train.shape[0], x_train.shape[1], x_train.shape[2], 1)
    y_train = keras.utils.to_categorical(y_train, num_classes=10)

    # Make a model
    opt = keras.optimizers.Adamax(lr=0.002, beta_1=0.9, beta_2=0.999,
                                  epsilon=1e-08, decay=1e-4)
    model = make_model()
    model.compile(loss='categorical_crossentropy', optimizer=opt, metrics=['accuracy'])

    # Train the model
    dg = data_gen(x_train, y_train, eh.batch_size, eh.alpha, 10)
    hist = model.fit_generator(dg, steps_per_epoch=len(x_train)/eh.batch_size, epochs=TRAIN_EPOCH,
                               validation_data=(x_val, y_val)).history

    # Evaluate the model by train and validation data
    score_train = model.evaluate(x_train, y_train, batch_size=eh.batch_size)
    score_val = model.evaluate(x_val, y_val, batch_size=eh.batch_size)

    comments = "#Final model loss_train={0:10.6f} acc_train={1:10.6f} loss_val={2:10.6f} acc_val={3:10.6f}"\
               .format(score_train[0], score_train[1], score_val[0], score_val[1])
    eh.write("history.log", model, opt, hist, comments)

if __name__ == '__main__':
    eh = ExperimentHistory()
    eh.batch_size = 32
    eh.alpha = 2
    run(eh)

2017年10月15日日曜日

機械学習における著作権

投稿時点(2017/10/15)の「日本の法律」では、機械学習に放り込む目的で著作物を使う場合、著作権が制限される条件(著作権法47条の7)があるため、商用・非商用にかかわらず自由に利用することができるようだ。

ITmediaの記事
http://www.itmedia.co.jp/news/articles/1710/10/news040.html

元の記事(上野達弘さんの投稿)
http://rclip.jp/2017/09/09/201708column/

IJCAIでの発表の記事
https://www.businessinsider.jp/post-102878


法律を引用すると
-----------------
http://www.houko.com/00/01/S45/048.HTM
(情報解析のための複製等)
第四七条の七 著作物は、電子計算機による情報解析(多数の著作物その他の大量の情報から、当該情報を構成する言語、音、影像その他の要素に係る情報を抽出し、比較、分類その他の統計的な解析を行うことをいう。以下この条において同じ。)を行うことを目的とする場合には、必要と認められる限度において、記録媒体への記録又は翻案(これにより創作した二次的著作物の記録を含む。)を行うことができる。ただし、情報解析を行う者の用に供するために作成されたデータベースの著作物については、この限りでない。
-----------------
とのことなので、別の法律による制限がないのであれば、確かに機械学習用のデータ利用では著作権の制限を受けないように読める。この規定の例外はデータベースで、おそらく、販売されているデータベース(新聞のコーパスやLDC(https://www.ldc.upenn.edu/language-resources))については、きちんと購入してライセンスを受けてください、ということなのだろう。

2017年10月9日月曜日

加速度センサーで行動分類(CNN編)

前回は加速度センサーデータをLSTMで分類しましたが、今回はCNNで試してみます。

センサーデータは、LSTMのときと同様にランダムに窓で区切ります。ただし、そのまま入力はせず、3次元のヒストグラムを作り、それをCNNの入力とします。センサーデータは3軸あり、振幅は整数でその範囲は[0,64)です。単純に、64x64x64の3次元配列でヒストグラムを作っても良いのですが、大部分が0になるスカスカのヒストグラムになる点と、CNNの規模が大きくなる点を考慮して、振幅を4分の1にしてヒストグラムを作ります。つまり、16x16x16の3次元配列(テンソル)をCNNの入力にします。これでも要素の数は4096個あり、窓幅を64とするとスカスカです。

使用するCNNのモデルは、3Dの畳み込み層3つと全結合層2つです。具体的には、

BIN_SIZE = 4
si = int(64/BIN_SIZE)
model = Sequential()
model.add(Conv3D(32, kernel_size=(5, 5, 5),
                 input_shape=(si, si, si, 1), 
                 data_format='channels_last'))
model.add(LeakyReLU(0.2))
for _ in range(2):
    model.add(Conv3D(32, kernel_size=(5, 5, 5), data_format='channels_last'))
    model.add(LeakyReLU(0.2))
model.add(Flatten())
model.add(Dense(256))
model.add(LeakyReLU(0.2))
model.add(Dropout(0.5))
model.add(Dense(num_classes))
model.add(Activation('softmax'))
の通りです。Conv3Dのフィルタサイズは32としました。Conv3Dの出力は順に、12x12x12x32、8x8x8x32、4x4x4x32です。したがって、1つ目のDenseへの入力サイズは2048となります。Conv3Dで適度なサイズに絞らないと、全結合層のパラメータ数が増えすぎて、GPUで処理できなくなります。

3次元のヒストグラムはデータを増やすためにランダムに回転させます。回転範囲の制限は、回転行列の対角要素の全てが一定値θより大きいという条件を満たす回転行列のみを使用することにより実現します。

前回と同様、バッチサイズ32、1エポックあたりのステップ数は100で、100エポック学習させました。結果は下表の通りです。100エポック終了後のモデルで訓練、検証ともに32000サンプル生成し、評価しています。32000サンプルはランダムに選択しているので毎回基準は揺れます。

Train accuracy





Window width No rotation Random rotation (θ)
0.95 0.9 0.5 0 Full
32 83.7% 57.5% 52.9% 40.2% 34.8% 29.7%
64 93.1% 70.4% 64.4% 49.8% 41.6% 35.9%
128 98.6% 85.5% 80.1% 63.4% 54.6% 47.6%
256 99.9% 95.0% 91.7% 79.5% 70.0% 61.8%
512 100.0% 98.7% 96.9% 89.6% 83.1% 75.6%

Validation accuracy
Window width No rotation Random rotation (θ)
0.95 0.9 0.5 0 Full
32 51.7% 51.2% 50.2% 38.5% 33.5% 25.3%
64 55.9% 59.9% 57.2% 46.1% 36.2% 31.9%
128 65.0% 70.8% 70.3% 58.6% 50.5% 37.3%
256 74.8% 80.6% 82.2% 77.1% 63.1% 51.9%
512 78.4% 80.7% 82.3% 83.1% 71.4% 65.0%

窓幅(Windows width)が32の場合、1秒のデータだけを見て分類していることになります。表から分かるように、入力データのランダム回転の範囲にかかわらず、窓幅が広いほど精度が高くなります。これは、窓幅が広くなるほど入力データのバリエーションが減り、分類に使えるデータが増えるためです。

入力データのランダムな回転の有効性は、窓幅によって変化します。窓幅が狭いと回転範囲が狭い(θが大きい)ほうが検証用データでの精度が良く、窓幅が広くなるにつれて精度が最も良くなる回転範囲が広めになります。ただし、どの窓幅でも回転範囲を広げすぎると精度が下がります。

今回の実験に使ったコードは次の通りです。util.random_rotateとutil.historyはここにあります。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
# This script trains a neural network model using accelerometer data
import sys,os
sys.path.append(os.path.dirname(os.path.abspath(__file__)) + '/../keras-examples/src')

import keras, datetime
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Activation, Conv3D, LeakyReLU, Flatten, Dropout
from keras.layers import LSTM
from util.random_rotate import uniform_random_rotation_matrix_3D
from util.history import ExperimentHistory

TRAIN_EPOCH = 100
BIN_SIZE = 4

def make_model(num_classes):
    si = int(64/BIN_SIZE) # Size of input
    model = Sequential()
    model.add(Conv3D(32, kernel_size=(5, 5, 5),
                     input_shape=(si, si, si, 1), 
                     data_format='channels_last'))
    model.add(LeakyReLU(0.2))
    for _ in range(2):
        model.add(Conv3D(32, kernel_size=(5, 5, 5), data_format='channels_last'))
        model.add(LeakyReLU(0.2))
    model.add(Flatten())
    model.add(Dense(256))
    model.add(LeakyReLU(0.2))
    model.add(Dropout(0.5))
    model.add(Dense(num_classes))
    model.add(Activation('softmax'))
    model.summary()
    return model

def read_accel(filename):
    with open(filename) as f:
        data = []
        for line in f:
            for v in line.rstrip().split(" "):
                if not v.isdigit():
                    print(v, "is not digit")
            data.append([int(v) for v in line.rstrip().split(" ")])
        data = np.array(data)
        return data

# Data generator
def data_gen(filelist, batch_size, window_width, random_rotate=False, rotation_range=0):
    # Make a dict for converting a class ID to files
    num_class = 0
    with open(filelist) as f:
        id2f = dict()
        for line in f:
            id, filename = line.split(" ", maxsplit=2)
            id = int(id)
            id2f.setdefault(id, []).append(read_accel(filename.rstrip()))
            num_class = max(num_class, id)
        num_class += 1
    yield num_class # Return the number of classes at first

    # Generate data eternally
    shift = np.array([32, 32, 32]) # Specify the center of rotation
    while True:
        data = []
        labels = []
        for i in range(batch_size):
            while True:
                c = np.random.randint(0, num_class) # Select a class
                if c in id2f:
                    break
            # Select a file in the class
            f = np.random.randint(0, len(id2f[c]))
            # Select a window position
            p = np.random.randint(0, max(0, len(id2f[c][f])-window_width)+1)
            # Make a label for the selected class
            labels.append(keras.utils.to_categorical(c, num_classes=num_class))
            # Add selected data
            if(random_rotate):
                while True:
                    m = uniform_random_rotation_matrix_3D()
                    if m[0,0] > rotation_range and m[1,1] > rotation_range and m[2,2] > rotation_range:
                        break
                data.append(np.dot(m, (id2f[c][f][p:p+window_width,:]-shift).T).T+shift)
            else:
                data.append(id2f[c][f][p:p+window_width,:])
                
            if(len(data[-1]) < window_width):
                padding = window_width - len(data[-1])
                data[-1] = np.append(np.zeros((padding, data[-1].shape[1]), dtype=np.int),
                                              data[-1], axis=0)
        # Make 3D-histogram (x, y and z range are 0-63)
        hist_data = []
        for i in range(batch_size):
            h = np.zeros(shape=(int(64/BIN_SIZE), int(64/BIN_SIZE), int(64/BIN_SIZE), 1))
            for t in range(data[i].shape[0]):
                if data[i][t,0] > 63 or data[i][t,0] < 0 or \
                   data[i][t,1] > 63 or data[i][t,1] < 0 or \
                   data[i][t,2] > 63 or data[i][t,2] < 0:
                    continue
                x = int(data[i][t,0]/BIN_SIZE)
                y = int(data[i][t,1]/BIN_SIZE)
                z = int(data[i][t,2]/BIN_SIZE)
                h[x, y, z, 0] += 1
            hist_data.append(h)
        yield np.array(hist_data), np.array(labels).reshape((batch_size, num_class))

def get_data(data_generator, steps):
    dgr = [next(data_generator) for _ in range(steps)]
    data = np.array([w for v in dgr for w in v[0]])
    labels = np.array([w for v in dgr for w in v[1]])
    return data, labels

def run(eh):
    dg = data_gen("train.list", eh.batch_size, eh.window_width, eh.random_rotate, eh.rotation_range)
    num_classes = next(dg)
    dgv = data_gen("val.list", eh.batch_size, eh.window_width)
    assert num_classes == next(dgv)

    # Make a model
    opt = keras.optimizers.Adamax(lr=0.002, beta_1=0.9, beta_2=0.999,
                                  epsilon=1e-08, decay=1e-4)
    model = make_model(num_classes)
    model.compile(loss='categorical_crossentropy', optimizer=opt, metrics=['accuracy'])

    # Train the model
    hist = model.fit_generator(dg, steps_per_epoch=100, 
                               validation_data=dgv, validation_steps=100,
                               epochs=TRAIN_EPOCH).history
    # Get data for evaluation
    score_train = model.evaluate(*get_data(dg, 1000))
    score_val = model.evaluate(*get_data(dgv, 1000))
    comments = "#Final model loss_train={0:10.6f} acc_train={1:10.6f} loss_val={2:10.6f} acc_val={3:10.6f}"\
               .format(score_train[0], score_train[1], score_val[0], score_val[1])
    eh.write("history.log", model, opt, hist, comments)

    # Evaluate the last model by train and validation data
    model.save("models/cnn_b{0}_ww{1}_rr{2}.hdf5"\
               .format(BIN_SIZE, eh.window_width, str(eh.rotation_range) if eh.random_rotate else "no"))

if __name__ == '__main__':
    if not os.path.exists("models"):
        os.mkdir("models")

    eh = ExperimentHistory()
    eh.BIN_SIZE = BIN_SIZE
    eh.batch_size = 32
    for (rr, rrange) in [(False, 0), (True, 0.95), (True, 0.9), (True, 0.5), (True, 0), (True, -2)]:
        eh.random_rotate = rr
        eh.rotation_range = rrange
        for ww in [32, 64, 128, 256, 512]:
            eh.window_width = ww
            run(eh)