theano CNNで抽出した特徴をファイルに出力

タイトルの通り。
CNNで特徴抽出して、今まではそのままMLPで識別していた。だが抽出した特徴を使って別の手法を試す必要が出てきたので、特徴をファイルに出力する方法をメモしておく。

CNNは畳み込みとプーリングの繰り返しにより、画像から自動的に特徴を学習し、抽出する。
オーソドックスなCNNだと、抽出した特徴はそのまま多層パーセプトロン(MLP)に入力されて識別が行われる。なので、MLPに入力される前のデータをそのまま取り出せばそれが特徴となっているわけである。
厳密に言うと別にMLP直前じゃなくても、任意の層の出力を取り出せば特徴が得られるのだが、識別器の違いによる性能の比較をするならば同じ箇所の特徴を使わなければならない。
ここの論文によると、CNNのどの階層から特徴を取るかによって性能がかなり変わってくるらしい。CNNをハナから特徴抽出器として使う場合は注意する必要がある。

参考元は以下
Pythonの基礎 ファイル(CSV)に書き込む編 - Pythonの学習の過程とか

早速やってみたいところだが、とりあえず学習済みネットワークがないと特徴抽出できない。
データセットは適当に用意。以前の記事参照。
で、学習。各layerはここ参照。
DeepLearningTutorialに準拠してるとこ多し

# train_set_x, train_set_yにデータとラベルを格納してある
# conv-pool-conv-pool-conv-pool-MLP-softmax
# image_size: 150x150  color

def gradient_updates_momentum(cost, params, learning_rate, momentum):
    assert momentum < 1 and momentum >= 0
    updates = []
    for param in params:
        param_update = theano.shared(param.get_value()*0., broadcastable=param.broadcastable)
        updates.append((param, param - learning_rate*param_update))
        updates.append((param_update, momentum*param_update + (1. - momentum)*T.grad(cost, param)))
    return updates

def training(learning_rate=0.001, n_epochs=300, momentum=0.9,
             nkerns=[16, 32, 48], batch_size=50, L1_reg=0.00, L2_reg=0.00):
    
    index = T.lscalar()
    
    x = T.matrix("x")
    y = T.ivector("y")
    
    print "building the model..."
    layer0_input = x.reshape((batch_size, 3, 150, 150))
    
    layerC0 = ConvLayer(
        input=layer0_input,
        image_shape=(batch_size, 3, 150, 150),
        filter_shape=(nkerns[0], 3, 11, 11),
        padding=0,
        st=3
        )
    
    layerP0 = PoolLayer(
        input=layerC0.output,
        poolsize=3,
        st=2
        )    
    
    layerC1 = ConvLayer(
        input=layerN_out,
        image_shape=(batch_size, nkerns[0], 24, 24),
        filter_shape=(nkerns[1], nkerns[0], 5, 5),
        padding=2
        )
        
    layerP1 = PoolLayer(
        input=layerC1.output,
        poolsize=3,
        st=2
        )

    
    layerC2 = ConvLayer(
        input=layerP1.output,
        image_shape=(batch_size, nkerns[1], 12, 12),
        filter_shape=(nkerns[2], nkerns[1], 5, 5),
        padding=2
        )
    
    layerP2 = PoolLayer(
        input=layerC2.output,
        poolsize=3,
        st=2
        )
        
        
    # hidden layer 全結合層を定義
    layer3_input = layerP2.output.flatten(2)
    layer3 = HiddenLayer(
        rng=numpy.random.RandomState(23455),
        input=layer3_input,
        n_in=nkerns[2] * 6 * 6,
        n_out=500,
        activation=T.tanh,
        )
        
    # logistic regression 層を定義
    layer4 = LogisticRegression(input=layer3.output, n_in=500, n_out=3)
    
    # 誤差関数NLLの数式を定義
    L1 = ( abs(layerC0.W).sum() + abs(layerC1.W).sum() + abs(layerC2.W).sum() 
            + abs(layer3.W).sum() + abs(layer4.W).sum() )
    L2_sqr = ( abs(layerC0.W).sum() + abs(layerC1.W).sum() + abs(layerC2.W).sum() 
            + abs(layer3.W ** 2).sum()  +  abs(layer4.W ** 2).sum() )

    cost = ( layer4.negative_log_likelihood(y)
        +  L1 * L1_reg
        +  L2_sqr * L2_reg
        )
    params = (layer4.params + layer3.params + layerC2.params 
            + layerC1.params + layerC0.params)

    updates = gradient_updates_momentum(cost, params, learning_rate, momentum)
    train_model = theano.function(
        [index],
        cost,
        updates=updates,
        givens={
            x: train_set_x[index * batch_size: (index + 1) * batch_size],
            y: train_set_y[index * batch_size: (index + 1) * batch_size]
        }
    )

    print 'training ...'
    epoch = 0
    while (epoch < n_epochs):
        epoch = epoch + 1
        for minibatch_index in xrange(n_train_batches):

            iter = (epoch - 1) * n_train_batches + minibatch_index

            if iter % 100 == 0:
                print 'training @ iter = ', iter
            cost_ij = train_model(minibatch_index)
            
    # 学習したネットワークのモデルを保存
    cPickle.dump(layerC0, open("layerC0.pkl","wb"))
    cPickle.dump(layerC1, open("layerC1.pkl", "wb"))
    cPickle.dump(layerC2, open("layerC2.pkl", "wb"))
    cPickle.dump(layer3, open("layer3.pkl", "wb"))
    cPickle.dump(layer4, open("layer4.pkl", "wb"))

パラメータはかなり適当。

その後、モデルを使って特徴を抽出する。
csvのwriterを使うと簡単にできる。

# 要csv, PIL, pylearn2

def feedforward_conv_pool(input, Clayer, pad, conv_st, pool_size, pool_st):
    input_shuffled = input.dimshuffle(1,2,3,0)
    filters_shuffled = Clayer.W.dimshuffle(1,2,3,0)
    conv_op = FilterActs(stride=conv_st, pad=pad, partial_sum=1)   #"full" convolution
    contiguous_input = gpu_contiguous(input_shuffled)
    contiguous_filters = gpu_contiguous(filters_shuffled)
    conv_out_shuffled = conv_op(contiguous_input, contiguous_filters)  
    conv_out = conv_out_shuffled.dimshuffle(3,0,1,2)
    
    pool_input = relu(conv_out + Clayer.b.dimshuffle('x', 0, 'x', 'x'))
    
    input_shuffled = pool_input.dimshuffle(1,2,3,0)
    pool_op = MaxPool(ds=pool_size, stride=pool_st)
    pooled_out_shuffled = pool_op(input_shuffled)
        
    pooled_out = pooled_out_shuffled.dimshuffle(3,0,1,2)
    
    return pooled_out

def feedforward_hidden(input, layer, activation=T.tanh, drop_rate=0):
    hidden_input = input.flatten(2)
    lin_output = T.dot(hidden_input, layer.W) + layer.b
    if drop_rate>0:
        output = lin_output * drop_rate
        return output
    else:
        output = activation(lin_output)
        return output

def feedforward_logistic(input, layer):
    output = T.nnet.softmax(T.dot(input,layer.params[0]) + layer.params[1])
    return output

def feature():
    layerC0 = cPickle.load(open("layerC0.pkl", "rb"))
    layerC1 = cPickle.load(open("layerC1.pkl", "rb"))
    layerC2 = cPickle.load(open("layerC2.pkl", "rb"))
    layer3 = cPickle.load(open("layer3.pkl", "rb"))
    layer4  = cPickle.load(open("layer4.pkl", "rb"))
    
    # 特徴抽出する画像を用意
    img = numpy.array(Image.open(...), dtype=theano.config.floatX)
    img = img.transpose(2,0,1)
    img = img.flatten()

    # 拡張子は.csvとか
    f = open("feature.csv","w")   
    dataWriter = csv.writer(f)
    x = T.matrix("x")
        
    layer0_input = x.reshape((1,3,150,150))
    layer0_out = feedforward_conv_pool(layer0_input, layerC0,
                                       pad=0, conv_st=3, pool_size=3, pool_st=2)
    layer1_out = feedforward_conv_pool(layer0_out, layerC1,
                                       pad=2, conv_st=1, pool_size=3, pool_st=2)
    layer2_out = feedforward_conv_pool(layer1_out, layerC2,
                                       pad=2, conv_st=1, pool_size=3, pool_st=2)
    layer3_out = feedforward_hidden(layer2_out, layer3)
    layer4_out = feedforward_logistic(layer3_out, layer4)
    # theano.functionで途中のlayerの出力を取り出す
    test_model = theano.function(
                    inputs=[layer0_input],
                    outputs=layer2_out,
                    )
    
    ft = test_model(img.reshape((1,3,150,150)))  # typeはnumpy.ndarray
    ft = ft.flatten()
    dataWriter.writerow(ft)
    f.close()

theano.functionでoutputsに途中の層を指定することで、入力データをネットワークの最後まで通さなくても特徴を取り出せる。

ftにpool layerの出力がndarrayで入る。
この時点では[1, channel, row, column]の形になっているので、一次元配列になおしてからwriterow()でファイルに書き込む。
二次元配列の場合はwriterows()を使うとループを使わず一発で書き込める。
実際はこっちのが使うかも。

出力したファイルはcsv形式なので、

[特徴1,特徴2,特徴3,   ....]

のようにカンマ(,)で区切られた数字の羅列になる。
自分はこれをmatlabで使いたかったので、csvread('feature.txt')ですぐ読み込めた。

【追記】
numpyのsavetxt()とloadtxt()でもできた。てかこっちのほうが簡単。

import numpy
a = numpy.zeros((3,3))
numpy.savetxt("test.txt",a)

読み込むときは

import numpy
a = numpy.loadtxt("test.txt")

でおけ
ただしこっちの方法だと全て小数でっぽい？桁数がすごいことになる
多分やり方はあると思うけど必要にかられてから調べよう。。

やっとこさ少しだけtheanoに慣れてきた気がしないでもない。