用过pytorch之后，现在使用prototxt构建网络，用shell脚本对网络训练总是感觉别扭。
幸好caffe提供了python接口，这篇博客就是使用python脚本实现Mnist分类，可以和上一篇文章对比学习Caffe:实现LeNet网络对Mnist做分类

#-*-coding:utf-8-*-
import os
# os.chdir('..')

import sys
# sys.path.insert(0, './python')
import caffe

from pylab import *
# %matplotlib inline

from caffe import layers as L
from caffe import params as P
import cv2

'''
我们将以一种简洁自然的方式将网络编写为python代码，
它将序列化为caffe的protobuf模型格式。
该网络期望从预生成的lmdbs中读取数据，
但使用memoryDatalayer也可以直接从ndarrays中读取数据。
'''
def lenet(lmdb, batch_size):
    # our version of LeNet: a series of linear and simple nonlinear transformations

    # ctor
    n = caffe.NetSpec()
    # data layer, output shape N,C,H,W=64,1,28,28
    n.data, n.label = L.Data(batch_size=batch_size, backend=P.Data.LMDB, source=lmdb,
                             transform_param=dict(scale=1. / 255), ntop=2)
    # conv1 layer, output shape N,C,H,W=64,20,24,24
    n.conv1 = L.Convolution(n.data, kernel_size=5, num_output=20, weight_filler=dict(type='xavier'))
    # pool1 layer, output shape N,C,H,W=64,20,13,13
    n.pool1 = L.Pooling(n.conv1, kernel_size=2, stride=2, pool=P.Pooling.MAX)
    # conv2 layer, output shape N,C,H,W=64,50,9,9
    n.conv2 = L.Convolution(n.pool1, kernel_size=5, num_output=50, weight_filler=dict(type='xavier'))
    # pool2 layer, output shape N,C,H,W=64,50,4,4
    n.pool2 = L.Pooling(n.conv2, kernel_size=2, stride=2, pool=P.Pooling.MAX)
    # full connect 1 layer(inner product 1 layer), output shape N,C*H*W=64,500
    n.ip1 = L.InnerProduct(n.pool2, num_output=500, weight_filler=dict(type='xavier'))
    # relu1 layer, output shape N,C*H*W=64,500
    n.relu1 = L.ReLU(n.ip1, in_place=True)
    # full connect 2 layer(inner product 2 layer), output shape N,C*H*W=64,10
    n.ip2 = L.InnerProduct(n.relu1, num_output=10, weight_filler=dict(type='xavier'))
    # loss layer(softmax & loss layer), output shape N,C*H*W=64,1
    n.loss = L.SoftmaxWithLoss(n.ip2, n.label)

    return n.to_proto()

'''
使用谷歌的Protobuf库将网络以序列化格式写入磁盘。
可以直接读取、写入和修改此描述文件。
'''
with open('../caffe/examples/mnist/lenet_auto_train.prototxt', 'w') as f:
    f.write(str(lenet('../caffe/examples/mnist/mnist_train_lmdb', 64)))  # batch_size = 64

with open('../caffe/examples/mnist/lenet_auto_test.prototxt', 'w') as f:
    f.write(str(lenet('../caffe/examples/mnist/mnist_test_lmdb', 100)))

'''
让我们看看训练网络。
'''
# os.system("cat ../caffe/examples/mnist/lenet_auto_train.prototxt")
'''
现在让我们看看学习参数，它也被写为prototxt文件。
我们使用的是带有动量、weights衰减和特定学习速率的SGD。
'''
# os.system("cat ../caffe/examples/mnist/lenet_auto_solver.prototxt")

'''
让我们选择一个设备并加载solver。
我们将使用SGD（带动量），但也可以使用Adagrad和Nesterov的加速梯度。
'''
caffe.set_device(0)
caffe.set_mode_gpu()

solver = caffe.SGDSolver('../caffe/examples/mnist/lenet_auto_solver.prototxt')

'''
为了了解我们网络的架构，
我们可以检查中间层的feature map（blob）和参数的维度（这些在以后操作数据时也将很有用）。

网络结构都存储在solver.net.blobs中,
solver.net.blobs是一个OrderedDict(有序字典).
'''

# each output is (batch size, feature dim, spatial dim)
print [(k, v.data.shape) for k, v in solver.net.blobs.items()]

# just print the weight sizes (not biases)
print [(k, v[0].data.shape) for k, v in solver.net.params.items()]

'''
开始前，让我们检查一下所有东西是否如我们所期望的那样装好。

我们将在训练和测试网络上进行一次forward，并检查它们是否包含了我们的数据。
'''
print solver.net.forward()  # train net
print solver.test_nets[0].forward()  # test net (there can be more than one)

'''
显示训练集前八张图片和对应的标签

solver.net.blobs['data'].data能够取到一个minibatch，即64张图片.

验证:
In [13]: solver.net.blobs['data'].data.shape
Out[13]: (64, 1, 28, 28)  # (N C H W)

In [16]: solver.net.blobs['data'].data[:8,0].shape
Out[16]: (8, 28, 28)

In [20]: solver.net.blobs['data'].data[:8,0].transpose(1, 0, 2).shape
Out[20]: (28, 8, 28)

In [21]: solver.net.blobs['data'].data[:8, 0].transpose(1, 0, 2).reshape(28, 8*28).shape
Out[21]: (28, 224) # (h,n*w)

其他层同理，比如conv1层：
In [4]: solver.net.blobs['conv1'].data.shape
Out[4]: (64, 20, 24, 24)
'''
# we use a little trick to tile the first eight images
imshow(solver.net.blobs['data'].data[:8, 0].transpose(1, 0, 2).reshape(28, 8*28), cmap='gray')
show()
# cv2.imshow("imshow1",solver.net.blobs['data'].data[:8, 0].transpose(1, 0, 2).reshape(28, 8*28))
# cv2.waitKey(0)
print solver.net.blobs['label'].data[:8]

'''
显示测试集前八张图片和对应的标签
'''

imshow(solver.test_nets[0].blobs['data'].data[:8, 0].transpose(1, 0, 2).reshape(28, 8*28), cmap='gray')
show()
# cv2.imshow("imshow2",solver.test_nets[0].blobs['data'].data[:8, 0].transpose(1, 0, 2).reshape(28, 8*28))
# cv2.waitKey(0)
print solver.test_nets[0].blobs['label'].data[:8]
'''
训练网络和测试网络似乎都在加载数据，并且都有正确的标签。
让我们从（minibatch）SGD中进行一个step，看看会发生什么。
'''
solver.step(1)
'''
我们是否有梯度通过filter传播？
让我们看看第一层(conv1)的更新.

conv1有20个filter,显然也有20个filter对应的梯度
所有的filter都有5*5的窗口大小,对应的filter的梯度显然也有5*5的窗口大小。
我们在4*5网格中显示所有的filter梯度。

In [24]: solver.net.params['conv1'][0].diff.shape
Out[24]: (20, 1, 5, 5)


'''

# reshape 4*5 grid,transpose gy h gx w,reshape gy*h gx*w
imshow(solver.net.params['conv1'][0].diff[:, 0].reshape(4, 5, 5, 5)
       .transpose(0, 2, 1, 3).reshape(4*5, 5*5), cmap='gray')
show()
#cv2.imshow("imshow3",solver.net.params['conv1'][0].diff[:, 0].reshape(4, 5, 5, 5).transpose(0, 2, 1, 3).reshape(4*5, 5*5))
#cv2.waitKey(0)

'''
现在让网络运行一段时间，我们来跟踪一些事情。请注意，此过程将与通过CAFFE二进制可执行文件进行训练的过程相同。特别地：
* 日志记录将继续正常进行
* 快照将在solver prototxt中指定的时间间隔内拍摄（此处，每5000次迭代一次）
* 测试将在指定的时间间隔内进行（此处，每500次迭代一次）

因为我们在python中控制了循环，所以我们可以自由地计算额外的东西，如下所示。我们也可以做很多其他事情，例如：
* 编写自定义停止条件
* 通过更新循环中的网络更改求解过程
'''


# % % time
niter = 200
test_interval = 25
# losses will also be stored in the log
train_loss = zeros(niter)
test_acc = zeros(int(np.ceil(niter / test_interval)))
output = zeros((niter, 8, 10))

# the main solver loop
for it in range(niter):
    solver.step(1)  # SGD by Caffe

    # store the train loss
    train_loss[it] = solver.net.blobs['loss'].data

    # store the output on the first test batch
    # (start the forward pass at conv1 to avoid loading new data)
    solver.test_nets[0].forward(start='conv1')
    output[it] = solver.test_nets[0].blobs['ip2'].data[:8]

    # run a full test every so often
    # (Caffe can also do this for us and write to a log, but we show here
    #  how to do it directly in Python, where more complicated things are easier.)
    if it % test_interval == 0:
        print 'Iteration', it, 'testing...'
        correct = 0
        for test_it in range(100):
            solver.test_nets[0].forward()
            correct += sum(solver.test_nets[0].blobs['ip2'].data.argmax(1)
                           == solver.test_nets[0].blobs['label'].data)
        test_acc[it // test_interval] = correct / 1e4

'''
让我们绘制训练损失和测试精度。
'''
_, ax1 = subplots()  # 左侧的y轴
ax2 = ax1.twinx()  # 右侧的y轴
ax1.plot(arange(niter), train_loss)
ax2.plot(test_interval * arange(len(test_acc)), test_acc, 'r')
ax1.set_xlabel('iteration')
ax1.set_ylabel('train loss')
ax2.set_ylabel('test accuracy')
show()
'''
损失似乎很快就消失了（除了随机性），而准确度也相应上升。万岁！

因为我们保存了第一批测试的结果，所以我们可以观察预测分数是如何演变的。
我们将在轴上绘制时间，并在上绘制每个可能的标签，其中的亮度表示信心。
'''
for i in range(8):
    figure(figsize=(2, 2))
    imshow(solver.test_nets[0].blobs['data'].data[i, 0], cmap='gray')
    figure(figsize=(10, 2))
    imshow(output[:50, i].T, interpolation='nearest', cmap='gray')
    xlabel('iteration')
    ylabel('label')
show()

'''
我们一开始对这些数字一无所知，最后对每个数字都进行了正确的分类。
如果你一直在跟踪，你会发现最后一个数字是最困难的，倾斜的“9”和“4”最混淆（可以理解）。
请注意，这些是“原始”输出分数，而不是softmax计算得到的概率向量。

softmax计算得到的概率向量如下图所示，使我们更容易看到我们网络的置信度
（但更难看到不太可能数字的分数）。
'''
for i in range(8):
    figure(figsize=(2, 2))
    imshow(solver.test_nets[0].blobs['data'].data[i, 0], cmap='gray')
    figure(figsize=(10, 2))
    imshow(exp(output[:50, i].T) / exp(output[:50, i].T).sum(0), interpolation='nearest', cmap='gray')
    xlabel('iteration')
    ylabel('label')
show()