Face Generation Using DCGAN in PyTorch based on CelebA image dataset 使用PyTorch打造基于CelebA图片集的DCGAN生成人脸

千呼万唤始出来的iPhone X有没有惊艳到你呢?还记得苹果的软件总监Craig Federighi发布会现场第一次尝试Face ID解锁失败的窘状吗?好在机智的Craig再次尝试的时候果断拿起桌上的另一只iPhone X, 并完美成功刷脸。当然,他还象征性地擦了下脸颊(冷汗!?)。😂😂😂

哈,言归正传。早在今年7月,苹果就发表了一篇名为“改进合成图像的仿真性”的技术博客。其中就提到利用GAN来生成更多的合成人眼图像。虽然未经证实,我相信这也是成就Face ID 的关键技术之一。

GAN全称是Generative Adversarial Nets,中文叫生成对抗网络,其本质是一种生成模型。而生成模型可以理解成非监督学习,即只有数据,没有标签。表征数据的一种有效方式是其概率密度分度。而为求其概率分布函数的参数,最大似然估计是一种有效的方式。如此,我们可以得到如下所示的公式。

Screen Shot 2017-09-23 at 15.25.11

然而,并非任意数据的概率分布都可以显示表征,因而生成模型可以进一步划分为一下类型。下面,让我们追随Stanford的CS231n课程来深入理解下3个典型的生成模型吧。

Screen Shot 2017-09-23 at 15.36.29

  • PixelRNN/CNN

首先,PixelRNN/CNN都属于Fully Visible Belief Nets。它利用链式法则将全概率展开成条件概率的乘积。显然,这种表征是显示而且复杂的。再有,这个表征不是跟上次讨论的语言模型如出一辙吗?只不过语言模型做了简化,只考虑与相邻前N个的依赖关系,即所谓N-gram模型。那么,顺理成章的上节讨论的RNN也适用与这里的Fully Visible Belief Nets生成模型,故而有了基于LSTM的PixelRNN。

Screen Shot 2017-09-23 at 15.52.20

PixelRNN选取图片的左上角元素作为起始点, 然后像训练序列一样训练LSTM网络,只不过相较于语言序列,图片序列是2维的。

Screen Shot 2017-09-23 at 15.56.26

既然是图片,为何不尝试用卷积呢?是啊,于是又有了PixelCNN。相较于PixelRNN,PixelCNN的训练速度更快,这和我们的一贯认知也一致,即卷积网络快过循环网络。不过,二者生成都是序列的,因而都不够快。

Screen Shot 2017-09-23 at 16.01.32

  • Variational AutoEncoder

首先不得不提标准的AutoEncoder。AutoEncoder的思想是训练一个Encoder使得输入x可以转化为维度更小的隐藏变量(latent variable)z,并且同时训练一个Decoder使得这个z确实可以还原原始的x。实际使用的时候,我们一般只保留Encoder,这相当于对数据x做了降维处理,并且使得得到的z可以最大程度保留x的variance。

显然,标准的AutoEncoder不足以成为一个生成模型,因为Decoder的输入z已经被限制成最大程度代表x。那是不是可以通过对z做某种处理来达到生成模型的效果呢?事实上,由此衍生出来的Variational AutoEncoder就是通过限定隐藏变量z服从高斯分布来实现一个生成模型的。具体的细节和公式推导就不在展开了。

Screen Shot 2017-09-23 at 17.28.29

有一点值得一提,Variational AutoEncoder是对目标函数的一个下限(lower bound)做优化,因此其效果不见得比上面提到的PixelRNN/CNN好。

  • GAN

不论是PixelRNN/CNN,还是VAE,它们都试图定义一个显示的概率密度函数。那么,如果我们换个思路呢?GAN就是DL界青年翘楚Ian Goodfellow于2014年在一场酒吧party后构思诞生的。

Screen Shot 2017-09-23 at 18.00.53

GAN的核心思想是采取博弈论(game theory)的手法,给定2个角色彼此对抗,最终达到均衡(equilibrium)状态。具体来说,角色1是生成器(Generator),角色2是判决器(Discriminaotr)。生成器的目标是生成以假乱真足以骗过判决器的图片,而判决器的目标则是慧眼识珠鉴定真伪。

Screen Shot 2017-09-23 at 18.58.48.png

GAN的minmax目标函数乍看起来可谓惊悚。不过稍加分解也不难理解。首先,判决器的输出是[0,1]区间的“概率”值,并且约定1表示真,0表示假,这可以通过sigmoid函数实现。其次,我们采用对数概率。如此,我们可以对判决器的输出采样,E(logD(x))是对真实图像的判决值的期望,E(log(D(G(z))))是对生成器生成的假图像的判决值得期望,显然,判决器的优化目标是使得前者值越大越好,后者越小越好。为了数学上的便利,我们将后者改写成E(log(1-D(G(z)))),所以,判决器的优化目标可写成max[E(logD(x)) + E(log(1-D(G(z))))],也就是下面总的minmax目标函数的内层。显然,外层的min是针对生成器的,这时,E(logD(x))由于与G无关,可以直接忽略,那么min[E(log(1-D(G(z))))]的含义如何理解呢?这也正是GAN的博弈所在,同样是E(log(1-D(G(z)))),判决器希望它越达越好,而生成器希望它越小越好,因为它的值越小,即E(log(D(G(z))))值越大,表示判决器误把假图当真了。

Screen Shot 2017-09-23 at 19.07.33.png

不过,log(1-D(G(z)))的数学性质不佳,表现在当1-D(G(z))的值较小时,即判决器误判时,log(1-D(G(z)))的梯度过于平坦,这对于仰仗梯度反向传播训练的神经网络太多宛如晴天霹雳,这样会直接导致训练初始由于梯度变化不大,参数的更新也会变化迟缓,整体的Loss也会久久维持在高位。好在已经有聪明的人发现改用-log(D(G(z)))就可以有效规避这个问题。

Screen Shot 2017-09-23 at 19.32.28

  • DCGAN

有了以上目标函数,接下来考虑下生成器和判决器的网络各自怎么搭建。DCGAN就是将卷积层引入生成器和判决器之中。这里,判决器输入的可以使原始图像,故其卷积层更类似普通的CNN。但是,生成器的输入则是高斯噪声,而输出则要与原始图像等大小。如下如所示,DCGAN的生成器像是一个逆向卷积网络,相对应的运算有的地方称之DeConv,有的为ConvTranspose(Tensorflow和PyTorch)。第一层将输入的100行噪音vector重塑成4x4x1024。然后紧接着我们使用batch normalization和leaky ReLU 激活函数。下一层即为卷积的转置,输入数据的深度会减半,而单张“图片”的长和宽则会双倍。紧接着依旧是batch normalization和leaky ReLU 激活函数。如此堆叠到生成“图片”的尺寸变成32x32x128。最后一层则转化为原始图片的大小64x64x3。当然,如果原始图片就是32x32x3,那么可以舍弃DCGAN的最后一层,让倒数第2曾直接输出32x32x3即可。

Screen Shot 2017-09-23 at 19.50.20

下面我重点看下DCGAN的mini-batch训练代码。与下图所示的GAN的mini_batch训练算法基本一致,不同的是DCGAN的目标函数改成Binary Cross Entropy, 因此自带负号,因此可以通过梯度下降寻找最小值。

Screen Shot 2017-09-23 at 21.41.10.png

criterion = nn.BCELoss() #Binary Cross Entropy做目标函数

input = torch.FloatTensor(batchSize, 1, imageSize, imageSize)
noise = torch.FloatTensor(batchSize, nz, 1, 1)
fixed_noise = torch.FloatTensor(batchSize, nz, 1, 1).normal_(0, 1)
label = torch.FloatTensor(batchSize)
real_label = 1
fake_label = 0

fixed_noise = Variable(fixed_noise)

#setup optimizer

optimizerD = optim.Adam(netD.parameters(), lr=lr, betas=(beta1, 0.999))
optimizerG = optim.Adam(netG.parameters(), lr=lr, betas=(beta1, 0.999))

for epoch in range(niter):
for i, data in enumerate(dataloader, 0):
############################

#(1) 更新判决器: 最小化- log(D(x)) – log(1 – D(G(z)))

###########################

#真实图像

netD.zero_grad()
real_cpu, _ = data
batch_size = real_cpu.size(0)

input.resize_as_(real_cpu).copy_(real_cpu)
label.resize_(batch_size).fill_(real_label)
inputv = Variable(input)
labelv = Variable(label)

output = netD(inputv)

#errD_real is binary cross entropy

errD_real = criterion(output, labelv)
errD_real.backward()
D_x = output.data.mean()

#生成器合成图像

noise.resize_(batch_size, nz, 1, 1).normal_(0, 1)
noisev = Variable(noise)
fake = netG(noisev)
labelv = Variable(label.fill_(fake_label))
output = netD(fake.detach())
errD_fake = criterion(output, labelv)
errD_fake.backward()
D_G_z1 = output.data.mean()

#判决器总的损失函数 = errD_real + errD_fake

errD = errD_real + errD_fake
optimizerD.step()

############################

#(2) 更新生成器: 最小化- log(D(G(z)))

###########################
netG.zero_grad()
labelv = Variable(label.fill_(real_label)) # fake labels are real for generator cost
output = netD(fake)
errG = criterion(output, labelv)
errG.backward()
D_G_z2 = output.data.mean()
optimizerG.step()

下图是利用GCGAN训练MNIST 60000个数据5个epoch的结果。

output.gif

下图是利用GCGAN训练CelebA20000个数据21个epoch的结果。

output.gif

DCGAN的一个缺点就是没有一个有效的性能指标。比如,上面CelebA的训练结果中,Epoch0和Epoch21下的生成器与判决器的loss均看不出明显下降,然而却是此时生成的人脸确实可以做到以假乱真了。性能的豪华评估基本上就是靠人肉眼看了。😶😶😶

[0/50][0/313] Loss_D: 1.9310 Loss_G: 5.5808 D(x): 0.6005 D(G(z)): 0.6739 / 0.0052
[0/50][1/313] Loss_D: 0.6302 Loss_G: 6.7364 D(x): 0.9122 D(G(z)): 0.3692 / 0.0016
[0/50][2/313] Loss_D: 0.3457 Loss_G: 5.9983 D(x): 0.8600 D(G(z)): 0.1234 / 0.0038
[0/50][3/313] Loss_D: 0.4428 Loss_G: 6.2101 D(x): 0.8989 D(G(z)): 0.2130 / 0.0028
[0/50][4/313] Loss_D: 0.3170 Loss_G: 6.2056 D(x): 0.8961 D(G(z)): 0.1283 / 0.0027
[0/50][5/313] Loss_D: 0.2949 Loss_G: 7.0226 D(x): 0.9300 D(G(z)): 0.1543 / 0.0012
[0/50][6/313] Loss_D: 0.3280 Loss_G: 6.4508 D(x): 0.9103 D(G(z)): 0.0819 / 0.0022
[0/50][7/313] Loss_D: 0.3757 Loss_G: 8.7778 D(x): 0.9562 D(G(z)): 0.2443 / 0.0002
[0/50][8/313] Loss_D: 0.1574 Loss_G: 7.7744 D(x): 0.9128 D(G(z)): 0.0337 / 0.0008
[0/50][9/313] Loss_D: 0.2521 Loss_G: 6.3284 D(x): 0.8888 D(G(z)): 0.0744 / 0.0026
[20/50][306/313] Loss_D: 0.4354 Loss_G: 3.0039 D(x): 0.8895 D(G(z)): 0.2443 / 0.0642
[20/50][307/313] Loss_D: 0.5098 Loss_G: 2.6394 D(x): 0.7712 D(G(z)): 0.1861 / 0.0986
[20/50][308/313] Loss_D: 0.3870 Loss_G: 3.6130 D(x): 0.9053 D(G(z)): 0.2309 / 0.0356
[20/50][309/313] Loss_D: 1.1355 Loss_G: 0.9685 D(x): 0.4427 D(G(z)): 0.1269 / 0.4505
[20/50][310/313] Loss_D: 1.4342 Loss_G: 5.7484 D(x): 0.9524 D(G(z)): 0.6657 / 0.0056
[20/50][311/313] Loss_D: 0.8923 Loss_G: 1.7463 D(x): 0.4965 D(G(z)): 0.0415 / 0.2323
[20/50][312/313] Loss_D: 0.3199 Loss_G: 2.1348 D(x): 0.8800 D(G(z)): 0.1577 / 0.1307
[21/50][0/313] Loss_D: 0.6148 Loss_G: 4.3992 D(x): 0.9183 D(G(z)): 0.3477 / 0.0207
  • WGAN

WGAN基于DCGAN对GAN的训练过程做了如下调整。

  1. 判决器的最后输出不再经过sigmoid,因此其值也不再受限于0到1之间;
  2. 判决器和生成器的损失函数相应也发生变化,不再取对数;
  3. 优化器建议选择RMSprop而不是Adam;
  4. 判决器每次更新时参数截断在[-c, c]范围内;
  5. 不同于DCGAN训练时,判决器和生成器每个迭代里面都只做一次,WGAN在训练时,每个迭代里判决器训练5次,生成器训练1次。

Paper给出的WGAN训练算法如下。

WGAN.jpg

input = torch.FloatTensor(batchSize, 1, imageSize, imageSize)
noise = torch.FloatTensor(batchSize, nz, 1, 1)
fixed_noise = torch.FloatTensor(batchSize, nz, 1, 1).normal_(0, 1)
label = torch.FloatTensor(batchSize)
real_label = 1
fake_label = 0

 

fixed_noise = Variable(fixed_noise)

#setup optimizer

#WGAN Change — 不要用基于动量的优化算法(包括momentum和Adam),推荐RMSProp

optimizerD = optim.RMSprop(netD_wgan.parameters(), lr = lr)
optimizerG = optim.RMSprop(netG_wgan.parameters(), lr = lr)

for epoch in range(niter):
for i, data in enumerate(dataloader, 0):
############################

#(1) Update D network: maximize D(x) + (1 – D(G(z)))

###########################

make_trainable(netD_wgan, True)

#WGAN change — 判别器训练n_critic次

n_critic = (100 if (i < 25) or i % 500 == 0 else 5)

for j in range(n_critic):

#WGAN change — 每次更新判别器的参数之后把它们的绝对值截断到不超过一个固定常数c

for p in netD_wgan.parameters():
p.data.clamp_(-0.01, 0.01)

#train with real

netD_wgan.zero_grad()
real_img, _ = data
batch_size = real_img.size(0)

input.resize_as_(real_img).copy_(real_img)
label.resize_(batch_size).fill_(real_label)
inputv = Variable(input)
labelv = Variable(label)

output = netD_wgan(inputv)

#errD_real is binary cross entropy

#errD_real = criterion(output, labelv)

#WGAN change — 生成器和判别器的loss不取log

errD_real = output.mean()
errD_real.backward(labelv)
D_x = output.data.mean()

#train with fake

noise.resize_(batch_size, nz, 1, 1).normal_(0, 1)
noisev = Variable(noise)
fake = netG_wgan(noisev)

labelv = Variable(label.fill_(fake_label))
output = netD_wgan(fake.detach())

#errD_fake = criterion(output, labelv)

#WGAN change — 生成器和判别器的loss不取log

errD_fake = output.mean()
errD_fake.backward(labelv)
D_G_z1 = output.data.mean()

#errD = errD_real + errD_fake

#WGAN change — 判决器的loss新方法

#to minimize errD

errD = errD_fake – errD_real
Wasserstein_D = (errD_real – errD_fake).data.mean()
optimizerD.step()

############################

(2) Update G network: maximize D(G(z))

###########################
make_trainable(netD_wgan, False)

netG_wgan.zero_grad()

labelv = Variable(label.fill_(real_label)) # fake labels are real for generator cost

noise.resize_(batch_size, nz, 1, 1).normal_(0, 1)
noisev = Variable(noise)
fake = netG_wgan(noisev)
output = netD_wgan(fake)

#errG = criterion(output, labelv)

#WGAN change — 生成器和判别器的loss不取log

#to minimize errG

errG = -(output.mean())
errG.backward(labelv)
D_G_z2 = output.data.mean()
optimizerG.step()

 

然而,我实际运行的结果是WGAN没法收敛。不知道是代码哪里存在bug。知乎上有一篇中山大学郑华滨的博文总结的WGAN写得非常透彻,推荐一读。

Advertisements

Chinese WuYan Poetry Writing using LSTM 用LSTM写五言绝句

在介入LSTM之前,让我们回顾下Language Modeling的历程。先不论采取何种模型,Language Modeling的实质是建立单词序列的联合概率分布, 而实际中的做法则是利用链式法则展开成条件概率的乘积,表示如下。

Screen Shot 2017-09-04 at 00.50.33

  • Count based N-gram

该模型假定单词序列可以用马尔科夫链表示,进而将上面的联合概率分布简化为如下(2-gram)。其缺点很显然就是舍弃了单词序列可能存在的长依赖关系。

Screen Shot 2017-09-04 at 00.52.47

  • N-gram Neural Network

神经网络(3-gram)的输入是前2个单词,通过隐藏层再经由softmax得到向量P,其维度和词汇量相等,所以输出的单词就是P中概率最大的那个所对应的单词。不难看出,该模型依然难以掌握单词序列可能存在的长依赖关系(除非加长输入单词序列量,但这样也会显著增大神经网络的参数量)。Screen Shot 2017-09-04 at 01.05.40

  • Recurrent Neural Network

RNN设计之初就是为了处理序列,那应用到Language Model上也是自然而然的事。对比上下两张图,RNN不仅输入当前单词,而且会把前一个隐藏层输出h当成输入,这样的话,RNN就自然胜任不同的N-gram model,而不需要像前两种模型那样受限。注意上下两张图都是timestep的展开图,换句话说,h1到h4对应的参数矩阵V和向量c都是共享的。Screen Shot 2017-09-04 at 01.17.34.png

然而,RNN想要学习单词序列可能存在的长依赖关系也不太容易。这里就要谈谈梯度消失和梯度爆炸的问题。通过链式法则,可以推导出costN对h1的偏导包含参数矩阵Vh的N-1次累乘。所以,当Vh的最大特征值大于1时,梯度也会随之呈指数增长(爆炸),当Vh的最大特征值小于1时,梯度也会随之呈指数减小(消失)。

Screen Shot 2017-09-04 at 02.20.47Screen Shot 2017-09-04 at 02.21.21

解决梯度爆炸,我们尚且可以通过添加一个门阈值来控制。但是,梯度消失则需要调整从神经网络本身结构。Long Short-Term Memory应运而生。事实上,正如其名,LSTN就是为了解决RNN存在的梯度消失问题,从而更好的学习单词序列可能存在的长依赖关系。

下面我们来剖析下LSTM的内在结构。LSTM的核心是所谓的“细胞状态”,如下图的Ct。它与Xt一同决定ht,一个LSTM“细胞”单元将以tuple的形式存储(Ct,ht)。

Screen Shot 2017-09-04 at 02.43.10

首先,遗忘门 (forget gate)控制多大程度保留Ct-1。ft取值范围是sigmoid函数的输出范围,即[0,1]。其意义体现在LSTM多大程度上需要保留之前的单词序列的信息。

Screen Shot 2017-09-04 at 02.51.54

接着,输入门(input gate)控制多大程度上保留当前输入信息。而Ct(加上标)则代表当前的输入信息。

Screen Shot 2017-09-04 at 02.56.40.png

现在,我们可以计算得到Ct了。Ct由过去的ftCt-1和现在的itCt(上标)共同组成。

Screen Shot 2017-09-04 at 03.01.54.png

最后,输出门(output gate)控制多大程度上输出Ct成为ht。

Screen Shot 2017-09-04 at 03.04.37.png

所以,说了这么多,LSTM究竟是如何解决梯度消失的问题的呢?不难看出,如果遗忘门的值尽量接近1的话,下图的偏导数是可以成立的。(顺便也符合Andrew Karpathy吐槽下这个遗忘们门的叫法,明明叫记忆门才更加准确嘛!)

Screen Shot 2017-09-04 at 03.14.03

好了,铺垫这么多,无非是说明现在LSTM基本是RNN的标配。那具体RNN的应用场景也是多样的。下图所示的后两种many to many也经常叫做seq2seq。如果输入输出的序列长度是相同的,那么可以直接采取many to many的策略;如果输入输出的序列长度不相同,比如机器翻译,那么需要先采取many to one的策略,将变长的序列“编码”为固定长度,再采取many to many的策略“解码”为输出变长序列。

Screen Shot 2017-09-04 at 03.33.12.png

本文的应用即简单的固定长度的seq2seq模型。要说固定序列,第一反应就是我们的五言/七言绝句。下面就介绍下我是如何利用LSTM学习写诗的。

  • 首先,读入诗句。

# \u3002是中文句号的unicode编码,\uff0c是中文逗号的unicode编码。这里,检查当前诗句是不是五言绝句。
for p in content.split(u’\u3002′)[:-1]:
for s in p.split(u’\uff0c’):
if len(s) != 5:
wuyan = False
break

  • 然后,将文本编码成数字。topNWords是指当前诗文经常使用的前N的单字。

def vocabularyDict(content, topNWords):
cnt = Counter()
for p in content:
for w in p:
cnt[w] +=1

content_vocabulary = [i for (i,j) in sorted(dict(cnt).items(), key=lambda x:x[1], reverse=True)[:topNWords]]
#Note that the index here starts from 1
word_to_int_dict = dict((w,i) for i,w in enumerate(content_vocabulary,1))
int_to_word_dict = dict((i,w) for i,w in enumerate(content_vocabulary,1))

#Insert UNK for unknown word at index 0
word_to_int_dict[‘UNK’]=0
int_to_word_dict[0] = ‘UNK’

return word_to_int_dict, int_to_word_dict

def encode2sequence(content, encoder,sequence_length, timestep_gap):
content_ful = ”.join(content)
content_ful_ex = [w for w in content_ful if w != u’\u3002′ and w != u’\uff0c’]
content_int = [encoder[w] if w in encoder else encoder[‘UNK’] for w in content_ful_ex]

x_rnn, y_rnn = divide_sequences(content_int,sequence_length,timestep_gap)
return content_int, x_rnn, y_rnn

  • 接着,生成LSTM批训练时需要的输入序列矩阵和输出序列矩阵。

def divide_sequences(content_int, sequence_length, timestep_gap):
#Both c_in_dat and c_out_dat are list of sequence_length elements and each element will have n_samples values
c_in_dat = [[content_int[i+n] for i in xrange(0, len(content_int)-1-sequence_length, sequence_length)] for n in range(sequence_length)]
c_out_dat = [[content_int[i+n] for i in xrange(timestep_gap, len(content_int)-sequence_length, sequence_length)] for n in range(sequence_length)]

xs = [np.stack(c) for c in c_in_dat]
ys = [np.stack(c) for c in c_out_dat]

#Magic happens with this stack operation
#x_rnn will turn into an 2-d array of shape [n_samples,sequence_length]
x_rnn=np.stack(xs, axis=1)
y_rnn=np.expand_dims(np.stack(ys, axis=1), -1)
return x_rnn, y_rnn

  • 再就是定义LSTM的RNN。

这里要提下,Keras虽然简洁,但是确实牺牲了一些变通性。比如Keras必须在定义model的时候定义batch_input_shape的batch_size,那么在做预测的时候也必须给足batch_size的输入才行。解决方法我参考了stackoverflow,这也是为什么我会定义2个model,一个用来训练,一个用来预测,二者结构和参数完全一样,只有batch_input_shape不同而而已。不知道Keras的Functional API是否可以解决这一问题。

  • 开始批训练。

除了慢,我真的没什么好说的。相较之下,CNN的卷积计算真是快的多了。要知道一个LSTM里面的参数矩阵可多了去了。

Epoch 1/5
57856/57856 [==============================] - 563s - loss: 6.4563   
Epoch 2/5
57856/57856 [==============================] - 561s - loss: 5.8773   
Epoch 3/5
57856/57856 [==============================] - 556s - loss: 5.5700   
Epoch 4/5
57856/57856 [==============================] - 556s - loss: 5.3272   
Epoch 5/5
57856/57856 [==============================] - 547s - loss: 5.1238

download

  • 大功告成,小试牛刀。

来欣赏下LSTM的大作吧。有句话怎说来着,一本正经胡说八道。哈哈。

床前明月光 疑是地上霜
尘十重君天 愁年人风道 
************
白日依山尽 黄河入海流
去离水庐来 杯歌原风看 
************
危楼高百尺 手可摘星辰
回照照稽山 容兴自旧客 
************
春眠不觉晓 处处闻啼鸟
人相知得然 风草无手去 
************
夕阳无限好 只是近黄昏
项金犹勤山 云为成目秋 
************
故国三千里 深宫二十年
峰听春晚君 叶自山上何 
************
人闲桂花落 夜静春山空
必邓城中柴 林洲心茅遥

具体的代码可以在我的GitHub找到。

Image Style Transfer Using Keras and Tensorflow 使用Keras和Tensorflow生成风格转移图片

fast.ai有关Image Style Transfer的介绍是我见过最好的。Image Style Transfer实质上就是生成一副新的Image,使得它与给定Content Image的内容可以尽量接近,同时又兼具Stype Image的风格。为简单考虑,我们分开考虑如何分别生成Content Image和Style Image。不过,这里所讲的“生成”并不是train一个神经网络Model的参数,相反,我们会从一副随机生成的图片开始,通过L-BFGS的迭代优化算法来最小化这幅图片和给定Content Image之间的基于像素的MSE,此时,每次迭代都会利用梯度下降更新我们生成的图片。

最初提到Image Style Transfer的论文设定我们利用卷积神经网络,那么我们这里也选取熟悉的VGG16。为了更好的保留图片信息,我们需要将原始VGG16的max pooling层替换成average pooling层。

#Max pooling after Convolution
x = MaxPooling2D((2, 2), strides=(2, 2), name=’block1_pool’)(x)

#Average pooling after Convolution
x = AveragePooling2D((2, 2), strides=(2, 2), name=’block1_pool’)(x)

下面我们考察下如何生成Content Image。
#Step 0 — Create a pre-trained VGG16 model with last 3 FC layers removed
#and pooling layer replaced with average pooling 首先当然是创建一个VGG16网络,include_top=False表示舍弃后3层全连层
model = VGG16_Avg(include_top=False)

#Step 1 — Load Content Image and Image Preprocessing 接着加载Content图片,crop剪裁并无他意,由于我们并不是用VGG16来做图片分类,因此没有必要一定遵循VGG16规定的224乘224的输入图片尺寸
img = Image.open(‘JasonRidge.jpeg’)
img = img.crop(box=(30,60,430,460))
img

Oricontent

#Image Preprocessing — Subtract mean value and convert from RGB to BGR to feed VGG16 VGG16要求输入图片都减去ImageNet的均值并转化为BGR格式
rn_mean = np.array([123.68, 116.779, 103.939], dtype=np.float32)
preproc = lambda x: (x – rn_mean)[:, :, :, ::-1]
deproc = lambda x,s: np.clip(x.reshape(s)[:, :, :, ::-1] + rn_mean, 0, 255)

img_arr = preproc(np.expand_dims(np.array(img), 0))
shp = img_arr.shape

%matplotlib inline
import matplotlib.pyplot as plt
plt.imshow(np.squeeze(img_arr))

Preprocontent.png

#Step 2 — Build Neural Network Graph (Model)
from scipy.optimize import fmin_l_bfgs_b
from scipy.misc import imsave
from keras import metrics
#We take the first Convolution output of the 5th also the last Convolution Block
#We can latet verify with earlier layer 我们选取最后一个卷积块的第一层输出作为衡量标准
layer = model.get_layer(‘block5_conv1’).output
layer_model = Model(model.input, layer)
targ = K.variable(layer_model.predict(img_arr))

#Step 3 — Build Neural Network Loss and Gradient (Optimization)
#Loss is a scalar value — MSE 注意loss是一个标量,而梯度是一个ndarrdy
loss = metrics.mse(layer, targ)
grads = K.gradients(loss, model.input)
loss_fn = K.function([model.input], [loss])
grad_fn = K.function([model.input], grads)

#Create a container class for both Loss and Gradients Evaluator类可以在优化阶段方便计算Loss和梯度
class Evaluator(object):
def init(self, loss_f, grad_f, shp):
self.loss_f, self.grad_f, self.shp = loss_f, grad_f, shp

def loss(self, x):
loss_ = self.loss_f([x.reshape(self.shp)])
return np.array(loss_).astype(np.float64)

def grads(self, x):
grad_ = self.grad_f([x.reshape(self.shp)])
return np.array(grad_).flatten().astype(np.float64)
evaluator = Evaluator(loss_fn, grad_fn, shp)

#Step 4 — Time to Optimize
#Define a method which will optimize (minimize) the Loss using L_BFGS algorithm solve_image方法会将每次迭代更新的图片保存下来
def solve_image(eval_obj, niter, x):
for i in range(niter):
x, min_val, info = fmin_l_bfgs_b(eval_obj.loss, x.flatten(),
fprime=eval_obj.grads, maxfun=20)

x = np.clip(x, -127,127)
print(‘Current loss value:’, min_val)
imsave(‘results/Content_gen_conv5_1_{i}.png’.format(i=i+1), deproc(x.copy(), shp)[0])
return x

#Generate a random ‘noise’ image which will be the starting point of the optimization iteration 这里生成我们最初的图片,其实就是随机数组成的噪音图片
rand_img = lambda shape: np.random.uniform(-2.5, 2.5, shape)/100
x = rand_img(shp)
imsave(‘results/Content_gen_conv5_1_0.png’,x[0])

#Start the optimization with 10 iterations
iterations = 10
x = solve_image(evaluator, iterations, x)

('Current loss value:', array([ 92.49624634]))
('Current loss value:', array([ 38.59172058]))
('Current loss value:', array([ 25.30461121]))
('Current loss value:', array([ 19.56064987]))
('Current loss value:', array([ 16.33964348]))
('Current loss value:', array([ 14.33158779]))
('Current loss value:', array([ 12.92133713]))
('Current loss value:', array([ 11.83940792]))
('Current loss value:', array([ 10.99672127]))
('Current loss value:', array([ 10.28565693]))

Step 5 — Visulization

from IPython.display import HTML
from matplotlib import animation, rc

fig, ax = plt.subplots()
def animate(i):ax.imshow(Image.open(‘results/Content_gen_conv5_1_{i}.png’.format(i=i)))

anim = animation.FuncAnimation(fig, animate, frames=iterations+1, interval=500)
HTML(anim.to_html5_video())

可惜wordpress没法上传视频,不然那种从无到有做转变的一定会让你惊艳到。

Image.open(‘results/Content_gen_conv5_1_10.png’)

上面左图就是最初的噪音图片,右图则是10次优化后生成的图片。

有了生成Content Image的经验,我们就来考察下如何做Style Transfer。其实核心的改动就在Step2和Step3。下面就着重介绍下要点。

#Step 2 — Build Neural Network Graph (Model)

outputs = {l.name: l.output for l in model.layers}
#Content Model和上面完全一致,只是我们选了更早的卷积层输出,以期可以保留尽量多的图片内容
content_name = ‘block3_conv2’
content_layer = outputs[content_name]
content_model = Model(model.input, content_layer)
content_targ = K.variable(content_model.predict(content_arr))
#Style Model则不同,它需要5个卷积层的输出
style_layers = [outputs[‘block{}_conv2’.format(o)] for o in range(1,6)]
style_model = Model(model.input, style_layers)
style_targs = [K.variable(o) for o in style_model.predict(style_arr)]

#Step 3 — Build Neural Network Loss and Gradient (Optimization)
#First create a fucntion to calculate the Gram Matrix for Style Loss 下面新定义的方法专门计算格拉姆矩阵
def gram_matrix(x):
#We want each row to be a channel, and the columns to be flattened x,y locations
features = K.batch_flatten(K.permute_dimensions(x, (2, 0, 1)))
#The dot product of this with its transpose shows the correlation
#between each pair of channels
return K.dot(features, K.transpose(features)) / x.get_shape().num_elements()

def style_loss(x, targ): return metrics.mse(gram_matrix(x), gram_matrix(targ))

#Style Loss is a scalar
style_wgts = [0.05,0.2,0.2,0.25,0.3]
loss_sty = sum(style_loss(l1[0], l2[0])*w for l1,l2,w in zip(style_layers, style_targs, style_wgts))

#Content Loss is a scalar
loss_con = metrics.mse(content_layer, content_targ)

#Total Loss is the linear sum of Style Loss and Content Loss
#But here the denominator 10 is just used to balance whether
#Style or Content is more emphasized 这里,我们兼顾Content Loss和Style Loss的方法就是取2者的和,Content Loss下面的分母10是用来控制2者的权重,分母越大,则Content Loss的权重就越小,相应生成的图片也会更加侧重保留风格而非内容。
loss = loss_sty + loss_con/10

#K.gradients = tf.contrib.keras.backend.gradients 有了总的Loss,梯度的处理和上面Image处理也完全一致。
grads = K.gradients(loss, model.input)

#K.function = tf.contrib.keras.backend.function
#Note that loss_fn will return a scalar Loss while grad_fn a numpy array Gradients
loss_fn = K.function([model.input], [loss])
grad_fn = K.function([model.input], grads)

好了,那我们来看下最终的成果吧。

看下帅大叔和百花图以及小萌妹的搭配是什么效果吧。

 

 

Image Classification based on VGG16 Transfer Learning using Tensorflow 使用Tensorflow打造基于VGG16的图片分类器

时下最流行的机器学习非深度学习莫属。而想要事半功倍入门深度学习就非Transfer Learning莫属。fast.ai和Udacity的深度学习入门课程不约而同都以VGG16为基础,介绍了如何使用卷积神经网络实现图片分类器。我也参照着打造了以下专门识别男女的图片分类器。

作为基础的VGG16,其模型结构如下图所示。它先后由2个2层卷积、3个3层卷积和3个全连层组成,共16层。其中第1层和第2层卷积都是64个卷积核,所以输出的维度也相应变成(224,224,64),再经由(2,2)max pooling变成(112,112,64)。以后可如此类推。

vgg16

Transfer Learning 的核心思想就是,最大限度利用已经训练好的模型,为了实现特定分类功能,只需要舍弃原模型的最后一层或若干层,然后拼上新定义层并做fine tuning以获得新定义层的weights。当然,保留下来的原模型的所有层在fine tuning过程中均保持weights不变。这样显然可以显著减少模型的训练开销,并且也可以取得不错的结果,从而达到事半功倍的效果。

在fast.ai的Keras示例中,其做法如下代码所示。

model.pop() #舍弃当前模型的最后一层
for layer in model.layers: layer.trainable=False #当前模型的所有层的weights保持不变
model.add(Dense(num, activation=’softmax’)) #拼上新定义的全连层
self.compile() #使模型的新网络结构生效

而在Udacity的Tensorflow实例中,其做法如下代码所示。

feed_dict = {input_: images} #输入图片到VGG16模型
codes = sess.run(vgg.relu6, feed_dict=feed_dict) #截取VGG16模型的第一层全连层输出

以下代码重新定义一个2层全连层网络,输入即上面截取的VGG16模型的第一层全连层输出,而输出则为我们需要的特定分类器输出标签。

inputs_ = tf.placeholder(tf.float32, shape=[None, codes.shape[1]])
labels_ = tf.placeholder(tf.int64, shape=[None, labels_vecs.shape[1]])

fc = tf.contrib.layers.fully_connected(inputs_, 256)
logits = tf.contrib.layers.fully_connected(fc, labels_vecs.shape[1], activation_fn=None)

两者相比较,示例1中fine tuning是在包含VGG16前15层和新加1层全连层的模型上开展,当然真正可调的也只是最后一层新加的全连层的参数;示例2则需要先单独截取VGG16的前14层输出,再作为新定义的2层全连层的输入做模型训练。看似做法略有不同,其实本质上思想是一致的。

下面,看下图片经过VGG16不同卷积层的输出大概是什么样子的。原始图片是图像处理领域的“明星”, Lena。

第1层max pooling输出,共64个卷积核。可以看出不同的卷积核关注的焦点也不同,得到的输出也是各有侧重。

conl1

第2层max pooling输出,共128个卷积核。

convl2

第3层max pooling输出,共256个卷积核。

convl3

第5层max pooling输出,共512个卷积核。好吧,现在的输出每张image只有7乘7的像素。肉眼已经完全“沦陷”了。

convl5

那我们看下打造的性别分类器的输出吧。

Screen Shot 2017-08-07 at 23.35.36

哈,话说金星老师你即便还是他的时候,分类器还是可以“慧眼识珠”的。

Screen Shot 2017-08-07 at 23.35.52.png

额,下面这位。。。哈哈,我已无言。

Screen Shot 2017-08-07 at 23.37.41.png

至此,简单几十行代码便可以实现图片分类,并且还可以任意推广到其他图片分类应用。看来深度学习的魅力和潜力真不是盖的。那么,就继续加油吧。

fast.ai Deep Learning Courses – Making neural nets uncool again fast.ai 深度学习课程学习自省和小结

今天很开心发现fast.ai更新了深度学习的第2期课程。看来沉寂了几个月的深度学习又可以鼓噪起来了。接下来,我会定期记录深度学习的各种重点和启示。

真心佩服和感谢课程的讲师JEREMY HOWARD和RACHEL THOMAS,堪称当代AI福音传道者。

1期课程:http://course.fast.ai/lessons/lessons.html

2期课程:http://course.fast.ai/lessons/lessons2.html

Topic Extraction on Medium.com Documents using NMF and LDA in Python使用Python NMF和LDA对Medium文章主题提取

最近发现一个满是资源和干货的博客网站Medium.com,正好也在倒腾自然语言处理,决定写个网络爬虫做点主题提取分析。Medium.com可以按topic添加自己感兴趣的文章,比如下图的人工智能topic。Medium.com不会一次性加载所有的文章列表,而是通过用户下拉经由JavaScript动态加载。简单起见,我写的爬虫也就刚好爬取首次加载的文章。
Screen Shot 2017-07-31 at 14.57.42

# Define a Web Scraper which will get documents from Medium.com for a given topic
def mediumScraper(topic):
    links = []
    content = []
    
    response = requests.get('https://medium.com/topic/'+topic)
    soup = BeautifulSoup(response.text,'lxml')
    
    for post in soup.findAll("a", { "data-action" : "open-post" }):
        links.append(post['href'].split('?')[0])
    
    for link in links:
        response2 = requests.get(link)
        soup2 = BeautifulSoup(response2.text,'lxml')
        paragraphs = ''
        for paragraph in soup2.findAll("p", { "class" : "graf graf--p graf-after--p" }):
            paragraphs += paragraph.text

        content.append(paragraphs)
        print('Found one document on topic %s' %(topic))
        
    return pd.DataFrame({'topic':topic,'links':links,'content':content},columns=['topic','links','content'])

# Scrape 8 topics and return as a data frame
topics = ['data-science','politics','technology','sexuality','artificial-intelligence','culture','software-engineering','equality']
df = mediumScraper(topics[0])
for topic in topics[1:]:
    df = pd.concat([df, mediumScraper(topic)],ignore_index=True)

由于有些文章只针对会员开放,所以content内容为空。删除2篇空文章后,最终得到的data frame共200行,既每个topic有25篇文章。

接下来就是NLP的tokenization标准流程了。得到的corpus只包含去除非字符和stopwords并取字根的转小写单词。

corpus = []
for i in range(0,200):
    # Tokenization
    # Remove any character except letters
    # Convert into lower case
    content = re.sub('[^a-zA-Z]',' ',df['content'][i]).lower()
    # Convert into stem word
    ps = PorterStemmer()
    content = ' '.join([ps.stem(word) for word in content.split(' ') if not word in set(stopwords.words('english'))])
    corpus.append(content)

紧接着就可以对corpus做数值化矩阵转化。简单来说就是词频转化。

# 简单的词频转化;适用LDA
cv = CountVectorizer(max_features=2000)
X = cv.fit_transform(corpus)
# TF-IDF转化;适用NMF
tf = TfidfVectorizer(max_features=2000)
X_2 = tf.fit_transform(corpus)

说实话,Latent Dirichlet Allocation要透彻理解真不是件容易的事。简单来说,LDA是一个三层贝叶斯概率模型,包含词、主题和文档三层结构。

Non-Negative Matrix Factorisation非负矩阵分解则要直观许多。如下图,不过这里需要做下调整,V就是X或者X_2,一行是一篇文档,而一行包含的列就是文档包含的词;W就是一会作用LDA或者NMF(transform方法)得到的矩阵,它把列从词转变成主题,显然,文档可以被划归为所在行主题权值最大的那个;H就是作用LDA或者NMF得到的矩阵componets_,我们可以取每行词权值最大的几个作为该主题的一个概要。
Screen Shot 2017-07-31 at 14.14.30

下面我们可以比较下LDA和NMF的结果。

  • 上面8组是LDA得到的主题包含的top 50词的wordcloud。下面8组是NMF得到的主题包含的top 50词的wordcloud。看起来,NMF的分组更能贴近Medium.com本身的topic [‘data-science’,‘politics’,‘technology’,‘sexuality’,‘artificial-intelligence’,‘culture’,‘software-engineering’,‘equality’]。

 

  • 下面左图是LDA,右图是NMF。图片是8行25列,既每个像素代表一个文档,且每行的25列来源于df的同一个topic。像素的颜色则代表LDA或者NMF给出的topic。若能同行每列颜色一致,两行之间颜色异同,则越能接近Medium.com自身的topic给定。显然,NMF的表现优于LDA。

     

  • 另外,NMF的计算复杂度也明显小于LDA。
    LDA done in 1.403s.
    NMF done in 0.022s.

最后,小结一下,以SVD,NMF为代表的矩阵讲解运算在机器学习,图像识别,文字挖掘等诸多领域都有广泛应用。单就文本主题提取而言,NMF在这里由于LDA。不过实际使用中,不知道应该怎么取舍。哈,我会持续关注和提高的。

 

Predict Shipment Time in Supply Chain Using Random Forest Regressor 用随机森林回归预测供应链货运时长

最近参加了公司组织的Hackthon活动。其中一项的任务就是数字化供应量的货运执行与管理。作为切入点,我选择了考察货运的运输时间及其影响因素,并利用机器学习的回归算法train了一个模型,专门用来预测货运的运输时间。下面总结下过程中遇到的问题以及解决方法。

其实,最初我们是打算做一个分类器来预测货运的及时交付与否。然而,在于供应链的一线同事交流之后才发现,实际中的运营是约20%的货运可以保证on time delivery,剩余80%的只是on time delivery。所以他们并不是特别在意货运是否一定要及时交付。既然如此,我们也就改成做货运时间的回归预测。

遇到的第一个问题,就是数据集里面并没有直接给出货运所花的时间。不过,倒是可以找出特定carton的进出站记录。例如,下面截取carton ID 为#5##74@52的记录。1月8号准备运输到1月16号送达总共花费8天。
Screen Shot 2017-07-17 at 09.24.24
#简单起见,我们先对数据集按carton ID分组再取时间戳的最大值与最小值之差作为每个carton ID的运输时间。
def getTimeNMode(group):
try:
result = pd.Series([(max(group[‘DATE’])-min(group[‘DATE’])).days,group[‘MODE’].unique()[0]],index=[‘ShippingDays’,’ShippingMode’])
return result
except:
return pd.Series([0,”],index=[‘ShippingDays’,’ShippingMode’])
#group_keys=False可以将’CONTAINER ID’,’STATUS’,’DATE’,’MODE’任然保留为columns而不是index
df1 = df[[‘CONTAINER ID’,’STATUS’,’DATE’,’MODE’]].groupby(‘CONTAINER ID’,group_keys=False).apply(getTimeNMode)

下图是ShippingDays的分布图。显然,不同时间段的ShippingDays也会略有不同。
ShippingDays

有了lable值,接下来开始选取feature。这里自然少不了对供应链物流环节的业务知识。目前,公司依赖一套OTM系统来完成货运的运输安排,对应数据集里的OTM_ROUTE_CODE。比如,
CC-SJZ,从工厂直接运到客户
COC-SLCHONGKONG,从工厂先运到就近的物流中心,再运到客户
CDC-SLCROERMOND,从工厂先运到客户就近的物流中心,再运到客户
CODC-SLCHONGKONG-SLCROERMOND,从工厂先运到就近的物流中心,再运到客户就近的物流中心,最后运到客户
#提取OTM_ROUTE_CODE的3个字段,作为模型的3个feature
def getInTransitStopOvers(row):
try:
code = row[‘OTM_ROUTE_CODE’]
#code = row
codes = []
if len(code) > 0 and code.count(‘-‘) >0:
codes = code.split(‘-‘)
if codes[0] in [‘COC’,’CDC’]:
codes.append(‘NA’)
elif codes[0] == ‘CC’:
codes[1] = ‘NA’
codes.append(‘NA’)
else:
codes = [”,”,”]
return pd.Series(codes,index=[‘RouteType’,’OriSLC’,’DesSLC’])
except Exception, e:
print code, ‘ failed in parsing In-Transit StopOvers with the error ‘+str(e)
return pd.Series([”,”,”],index=[‘RouteType’,’OriSLC’,’DesSLC’])

df[‘OTM_ROUTE_CODE’].fillna(”, inplace=True)
df[[‘RouteType’,’OriSLC’,’DesSLC’]] = df.apply(getInTransitStopOvers,axis=1)

然后考察运输承运商集运输方式。为此,考察SHIPPING_METHOD_CODE的字段。例如,
000001_DHLC_P_1D中,DHLC表示承运商为DHL,运输方式是P(普通包裹),优先级为1D(1天到达)。
#提取SHIPPING_METHOD_CODE的3个字段,作为模型的3个feature
def getDeliveryData(row):
try:
method = row[‘SHIPPING_METHOD_CODE’]
#method = row
methods = []
if method.count(‘‘) >2:
methods = method.split(‘
‘,3)
else:
methods = [”,”,”,”]
return pd.Series(methods[1:],index=[‘Carrier’,’DeliveryType’,’DeliveryPriority’])
except Exception, e:
print method + ‘ failed in parsing Delivery Data with the error ‘+str(e)
return pd.Series([”,”,”],index=[‘Carrier’,’DeliveryType’,’DeliveryPriority’])

df[‘SHIPPING_METHOD_CODE’].fillna(”, inplace=True)
df[[‘Carrier’,’DeliveryType’,’DeliveryPriority’]] = df.apply(getDeliveryData, axis=1)

接着将carton的尺寸拆分为长宽高3个地段,作为另外3个feature。
def getCarbonDimensions(row):
try:
dimension = str(row[‘CART_DIM’])
#dimension = row
dimensions = []
if dimension.count(‘X’) >1:
dimensions = [float(d.strip()) for d in dimension.split(‘X’,3)]
else:
dimensions = [0.0,0.0,0.0]
return pd.Series(dimensions,index=[‘Length’,’Width’,’Height’])
except Exception, e:
print dimension, ‘ failed in parsing dimensions with the error ‘+str(e)
return pd.Series([0.0,0.0,0.0],index=[‘Length’,’Width’,’Height’])

df[‘CART_DIM’].fillna(0.0, inplace=True)
df[[‘Length’,’Width’,’Height’]] = df.apply(getCarbonDimensions, axis=1)

最终我们选取下面12个feature。
SHIP_FROM_ORG_ID
TO_CITY
WEIGHT
ORDERED_ITEM
Width
Length
Height
RouteType
OriSLC
DesSLC
Carrier
DeliveryType

我们选取的RandomForestRegressor模型要求所有的feature必须是数值型的。为此,需要先对Categorical类型的feature进行数值编码。这里,我也遇到本次Hackthon最大的难点。我发现sklearn自带的labelEncoder虽然好用,但是,接下来的难题着实耗费了很多时间。

  • 由于我们做的是offline training,需要把模型序列化成一个可导入的文件。这样的话,每个做了转化的faeture的labelEncoder都要单独序列化,甚是麻烦。
  • 在做预测的时候,我的理解是输入的值需要经由同样的labelEncoder做转化。然而,此时labelEncoder永远输出0。一阵吐槽纠结后才发现,labelEncoder的fit或者fit_transform好像只能针对一个array起作用。换句话说,labelEncoder把预测数据当成一个新的一维数组进行重现转化,所以输出值永远为0。最终,我写了一个小方法,可能存在问题,但是能达到我需要的目的。
    from sklearn.externals import joblib
    def categoryEncoding(col,ColName):
    colDict = {}
    colSet = set([i for i in col])
    i = 1
    for j in colSet:
    colDict[j] = i
    i += 1
    joblib.dump(colDict, ‘/home/felix/jupyter/data/dict_’+ColName)
    return colDictdef encodingLookup(ColName, val):
    colDict = joblib.load(‘/home/felix/jupyter/data/dict_’+ColName)
    if val in colDict:
    return colDict[val]
    else:
    return colDict[max(colDict)] +1

接下来就是水到渠成的training和tuning。
from sklearn.cross_validation import train_test_split
from sklearn.model_selection import cross_val_score
from sklearn.ensemble import RandomForestRegressor
from sklearn.grid_search import GridSearchCV
cols = [‘TO_CITY’,’DeliveryType’,’RouteType’,’OriSLC’,’DesSLC’,’Carrier’,’ORDERED_ITEM’]
for colName in cols:
#le = LabelEncoder()
#df_model[colName] = le.fit_transform(df_model[colName])
colDict = categoryEncoding(df_model[colName],colName)
df_model[colName] = df_model[colName].apply(lambda x:colDict[x])

Split into Training and Testing set

train_set, test_set = train_test_split(df_model, test_size=0.2)
X_train = train_set.iloc[:,1:]
y_train = train_set.iloc[:,0]

X_test = test_set.iloc[:,1:]
y_test = test_set.iloc[:,0]

param_grid = [{‘n_estimators’: [400,600,800], ‘max_depth’:[3,5,7]}]
forest_reg = RandomForestRegressor()
grid_search = GridSearchCV(forest_reg, param_grid, cv=5)
grid_search.fit(X_train, y_train)

最终,我选取n_estimators=500,性能如下。
R2 = 0.62546409754312426
相应的Feature Importances排序如下。
TO_CITY
0.26507
WEIGHT
0.191827
ORDERED_ITEM
0.148944
Width
0.075103
Length
0.065954
Height
0.065687
OriSLC
0.050146
Carrier
0.041274
DesSLC
0.040447
SHIP_FROM_ORG_ID
0.038849
DeliveryType
0.015408
RouteType
0.001292

完整的代码参加我的GitHub。经历这次Hackthon,以下几点思考:

  • 你以为的不一定就是业界关心的问题。比如,我们从最初的分类转成回归分析。
  • 字符转数值的labelEncoder在真正的线上系统是如何应用的?
  • 随机森林回归需要考虑regularization吗?
  • 从Feature Importances来看,长宽高貌似对货运时间的影响还比较大,值得再做深入挖掘。

Cisco Encrypted Traffic Analytics 思科机密流量分析

思科近日对外宣布下一代基于可编程、自适应和高安全的革新网络(DNA+ETA)。

其中,最引人注目就是思科在业界率先利用机器学习算法实现了对网络加密流量的监督学习及分类,从而可以在无需解密的前提下判定网络攻击(malware)。具体的实现细节尚未完全披露,不过可以从公开发布的白皮书里略窥一二。

nb-09-encrytd-traf-anlytcs-wp-cte-en

Data Insights into Provisioning Orders Using R and Python 使用R和Python数说Provisioning Orders

数据源来自于Oracle DB导出的2份数据文件。其中,orderList是order流水记录,二orderRequests是每个order的详细XML格式request payload。
首先,借助Python将2份文件读入Pandas Data Frame。然后遍历orderRequest利用xml.etree.ElementTree解析order包含的具体服务类型和license,并添加到orderList作为新的列。最终汇总得到19633行119维的数据集。数据清洗的部分此处先忽略。

  • pandas.DataFrame.diff
    在分析change order中license的变化时,需要比较同一客户的service license前后变化。如果用SQL,得在做完group by service之后再重新join查询得到license,然后再写DB function实现一个group里面的前后license值变化。简直如同噩梦,果断放弃。这一功能恰好可以直接在分组上使用Pandas的diff方法,只需对差值稍作转换即可得到想要的license变化类型,可谓简洁明了。
    大体的Python代码如下。
    #根据diff差值解释license是Upsell,Downsell还是Unchanged
    def licenseUpsellDownsellUnit(x):

if x > 0:
return ‘Upsell’
elif x 1:
group_df = group.sort_values(by=’ServiceOrderID’).reset_index()
group_df[serviceName+’.License Change’] = group_df[serviceName+’.TOTAL LICENSE VOLUME’].dropna().diff().apply(licenseUpsellDownsellUnit)
group = group_df.set_index(keys=’index’)
else:
group[serviceName+’.License Change’] = ‘Only One Order Found’
return group
#按照产品(这里是SiteURL)分类
def licenseUpsellDownsellAll(df_service, serviceName):
groups = df_service[[‘ServiceOrderID’,’TransactionType’,serviceTotalLicense,’OrderRecieveDate’]].groupby(df_service[‘SiteURL’])
licenseChangeGroups = groups.apply(licenseUpsellDownsellByGroup)
return licenseChangeGroups

  • 等到图形化license change数据的时候,我犯了一个小错误。
    如果我告诉你下面2张图我的本意实际想要传递同样的信息,你是不是会会对图1表示抓狂?!
    其实图1看起来凌乱,实则包含一个很重要的信息,就是同一产品的licens change频率,这点可以从曲线的长度“直观”体现。折线每次上升或者下降即表明一次Upsell与Downsell,所以折线的一段越长,代表保持的时间越长,反之,越短则表示license很快变化。而放弃曲线改为散点的图2就丢失了这部分信息。当然,图2的纵轴可以直接看到license change的百分比。所以嘛,并不是信息量越大的可视化才叫好,相反,可视化的目标不就正是希望同图表代替数据来说话吗,越少烧脑,才能更好传达信息,不是吗?
    图1
    Plot9CMRUpsellDownsell-2
    图2
    Plot9CMRUpsellDownsell-1
  • 往R dyplyr方法内传入列名变量
    这点上Python可以碾压R了。在Python里面,把DataFrame的列名作为函数入参传入是很正常的做法,但是R里面却没那么直观了。
    比如我想定义一个对order里面的服务按周分组的function,当我把serviceName传入的时候,正常使用dyplyr里面的n()和sum()来计数是不行的。这时要改用SE function, 即summarise变成summarise_。同时,由于summarise_里面既有特殊变量serviceName也有一般的totalCount和serviceBookingRatio,参考vignette(“nse”)文档使用interp。
    #Use as.name if you have a character string that gives a variable name
    interp(~ mean(var), var = as.name(“mpg”))
    #> ~mean(mpg)
    #or supply the quoted name directly
    interp(~ mean(var), var = quote(mpg))

#将Data.Frame的列名serviceName作为入参传入
plotServiceBookingTrendByWeek<- function(instanceName,serviceName) { aggregatedServices % filter(InstanceName==instanceName, TransactionType==’Provide’|TransactionType==’Trransfer’|TransactionType==’Renewal’,ProductModule!=”) %>% group_by(year,week,ProductModule) %>% summarise_(totalCount = ~n(), serviceBookingRatio = interp(~ sum(service,na.rm=TRUE)/totalCount, service = as.name(serviceName)))
… ..

  • 绘制order在世界地图上的分布图
    最简单应该是调用GoogleVis封装好的gvisGeoChart,可惜它自身的iso2c识别针对小国家不太理想,比如新加坡和香港在生成的图里面竟然没有显示。
    最终,通过引入rworldmap和countrycode这2个库成功绘制。

#Count of unique Webex Sites per each Country
webexAnnuityOrdersAggrbyCountry % filter(InstanceName==’SITE’, CountryCode != ”) %>% group_by(CountryCode) %>% summarise(WebExSitesCount = n_distinct(SiteURL)) %>% arrange(desc(WebExSitesCount))

#使用GoogleVis封装好的gvisGeoChart
webexAnnuityOrdersAggrbyCountry$CountryName <- countrycode(webexAnnuityOrdersAggrbyCountry$CountryCode,’iso2c’, ‘country.name’)
Geo=gvisGeoChart(webexAnnuityOrdersAggrbyCountry, locationvar=”CountryName”,
colorvar=”WebExSitesCount”,
options=list(projection=”kavrayskiy-vii”))
plot(Geo)

#使用rworldmap和countrycode绘制
annuityOrders_map_poly has both world map and annuityOrdersAggrbyCountry data merged together
webeAnnuityOrders_map <- joinCountryData2Map(webexAnnuityOrdersAggrbyCountry, joinCode = “ISO2”, nameJoinColumn = “CountryCode”)
webexAnnuityOrders_map_poly <- ggplot2::fortify(webeAnnuityOrders_map)
webexAnnuityOrders_map_poly <- merge(webexAnnuityOrders_map_poly, webeAnnuityOrders_map@data, by.x=”id”, by.y=”ADMIN”, all.x=TRUE)

#Ashmore and Cartier Islands作为澳大利亚的附属岛屿可以被忽略不计
webexAnnuityOrders_map_poly$WebExSitesCount[webexAnnuityOrders_map_poly$id==’Ashmore and Cartier Islands’] <- NA
webeAnnuityOrders_map@data$WebExSitesCount[webeAnnuityOrders_map@data$ADMIN==’Ashmore and Cartier Islands’] <- NA
webeAnnuityOrders_map@data[‘Text’] <- paste(webeAnnuityOrders_map@data$ADMIN,webeAnnuityOrders_map@data$WebExSitesCount,sep = ‘:’) ggplot() + coord_map(xlim = c(-180, 180), ylim = c(-60, 75)) + geom_polygon(data = webexAnnuityOrders_map_poly, aes(long, lat, group=group,fill=WebExSitesCount), color=’white’, size = 0.3) + scale_fill_gradient2(trans= ‘log’, breaks=c(0,10,100,200,4000), labels=c(‘4000′)) + xlab(NULL) + ylab(NULL) + guides(fill = guide_legend( title=’Number of Annunity WebEx Sites Around the World \n(June 1,2016 – June 1,2017)’, keywidth = 0.7, keyheight = 0.7, reverse=F, title.position=”top”)) + theme( plot.background = element_blank() ,panel.grid.major = element_blank() ,panel.grid.minor = element_blank() ,panel.border = element_blank() ,axis.ticks = element_blank() ,axis.text.x = element_blank() ,axis.text.y = element_blank() ,legend.position = “bottom” ,legend.direction = “horizontal” ) + geom_text(data=subset(webeAnnuityOrders_map@data,WebExSitesCount>0), aes(x=LON, y=LAT,label=Text), size=2.0)
最终图如下。
Plot1WorldMapofWebexSites-1

  •  Sankey绘制
    Sankey又叫RiverFlow,是一种反应不同类型转化的有效可视化形式。考查change order里面同一service的不同type的变化情况可以使用。

MCSessionSwappedCount % filter(TransactionType==’Amend’|TransactionType==’Transfer’,WEBEX.CONFERENCING.Session.Swap==’Swapped’) %>% group_by(WEBEX.CONFERENCING.FORMER.SESSION.ID,WEBEX.CONFERENCING.SESSION.ID) %>% summarise(count=n())
#这里的一个trick就是前后的名字必须不同,否则绘制会报错。简单的处理方法就是认为添加空格制造不同。
MCSessionSwappedCount$WEBEX.CONFERENCING.FORMER.SESSION.ID[MCSessionSwappedCount$WEBEX.CONFERENCING.FORMER.SESSION.ID==3] <- ‘Pro Meeting ‘
… …

MCSessionSwappedCount$WEBEX.CONFERENCING.SESSION.ID[MCSessionSwappedCount$WEBEX.CONFERENCING.SESSION.ID==3] <- ‘Pro Meeting’
… …

sk_MC <- gvisSankey(MCSessionSwappedCount, from=”WEBEX.CONFERENCING.FORMER.SESSION.ID”, to=”WEBEX.CONFERENCING.SESSION.ID”, weight=”count”)

print(sk_MC, tag = ‘chart’)

Sankey

  • Treemap绘制
    Treemap尤其适合展示多层级关系。这里我用它来显示提交order的合作商及其对应的系统账号的分布。
    treemap(subset(webexAnnuityOrdersAggrbyPartnerAccount,partnerAccountName!=’Direct’),
    index=c(“partnerName”,”partnerAccountName”),

vSize = “CCWOrderCount”,
title=”Partner Orders”,
type=”index”,
title.legend=”number of Partner orders”,
fontsize.labels=10,
align.labels = list(c(“center”, “center”), c(“right”, “bottom”)),
overlap.labels=1,
lowerbound.cex.labels=0,
inflate.labels=FALSE,
aspRatio=1.4)

Plot2PartnerTreeMap-1

当然,完整的分析还涉及到很多metric的考察,由于涉及到商业信息,就不在一一赘述了。完整的Rmd代码上传至我的Github账号。