【论文实现】一篇Sigkdd的弹幕分析论文的python实现【LDA 实践者】

Author : Jasper Yang
School : Bupt

warning : 此篇文章基于较为熟悉GibbsLDA++的源码的前提下阅读。另外，这篇文章是我的一个很不成熟的笔记，里面的很多东西和我实现的最终版本的TPTM（这篇论文的模型）已经大相径庭了，所以这篇文章就当成一篇简单的记录吧，我还是放在我的blog里，不太想丢弃。-_-

论文link

在开始之前，我想提前说明一下，这篇文章是我实现了基于用户喜好和评论的前后影响的改进LDA，因为实验室在做弹幕分析这方面的研究，就找到了这篇A类会议的文章，本来找他们要源码的，但是他们一作并不理我。。。然后我只好自己实现了。

由于是第一次写这类文章，不知道在逻辑上表达的清不清楚，~.~，望有看的各位给点建议，哈。

怎么实现

first of first

    from scipy.special import digamma
    import eulerlib
    import numpy as np
    import math
    from scipy.special import gamma
    import jieba
    import pandas as pd
    from scipy import stats

以上是本文所需库。

let's begin(前期准备工作)

首先是处理数据

数据来源：bilibili（爬虫爬取）

数据格式如下。

                                                                                                    <?xml version="1.0" encoding="UTF-8"?><i>   <chatserver>chat.bilibili.com</chatserver><chatid>1007373</chatid><mission>0</mission><maxlimit>8000</maxlimit><source>e-r</source><ds>274694331</ds><de>3034701550</de><max_count>8000</max_count>
                                                                                                    <d p="72.409,1,25,16777215,1375542175,0,7526c714,274694331">我来组成弹幕..................</d>
                                                                                                    
<d p="33.551,1,25,16777215,1375543533,0,925384b2,274711742">大佬系边 </d>

<d p="117.977,1,25,16777215,1375543631,0,925384b2,274712904">甘嗨假噶 </d>

<d p="134.849,1,25,16777215,1375547487,0,D3dfe4a5,274755463">呢个日文..一个字都听唔明</d>
...

下面写了段代码解析获得上面的用户id以及每个用户的评论。
最后形成一个所有评论的列表以及每个用户对应其评论的字典。

    user = {}
    comments = []
    split_num = 10 # 分割一部电影为几个片段
    
    # 数据和停用词
    danmu = open('danmu/1007373.xml')
    stopwords = {}.fromkeys([ line.rstrip().decode('gbk') for line in open('stopwords.txt') ])
    # 读取文件，分析后存储到 user 和 comments
    for line in danmu.readlines()[:-1]:
        start = line.find('p=')
        stop = line.find('">')
        sub1 = line[start+3:stop]
        time = sub1.split(',')[0]
        sub1 = sub1.split(',')[6]
        start = line.find('">')
        stop = line.find('</d>')
        sub2 = line[start+2:stop].decode('utf-8')
        comments.append((float(time),sub2))
        temp = []
        if not user.has_key(sub1) :
            temp.append(sub2)
            user[str(sub1)] = temp
        else:
            user[str(sub1)].append(sub2)

经过以上处理后，我们还需要分片，我这里默认分割成十份。

    # 统计user的个数 , 现在统计的是这个文档里的user，后期要做成对所有文档的统计量，还要能支持增量
    user_num = len(user)
    
    # comments的数量
    comments_num = len(comments)
            
    # 排序，分割comments ---> shots
    comments = sorted(comments)
    spli = (comments[-1][0]-comments[0][0])/split_num
    shots = []
    for i in range(10):
        shots.append([x[1] for x in comments if x[0] > i*spli and x[0] <= (i+1)*spli ])

分割之后就是分词，我这里用的是python版的jieba，效果很一般，以后再更新。

注意：这里的切词我分成了两个部分。因为在文中认为片段(shot)和片段之间有时间上的影响，然后每个片段里的每一个comment之间有时间上的影响。但是不同的片段之间的comment是没有关系的。

def cut_word(x):
    words = jieba.cut(x, cut_all=False)
    final = []
    for word in words:
        if word not in stopwords:
            final.append(word)
    return final

def null(l):
    return len(l) > 0 or l[0] == ' '

for i in range(split_num):
    shots[i] = map(cut_word,shots[i])
    shots[i] = filter(null,shots[i])

上面的代码对每一个shot里面的所有comment进行切词。

准备工作到此为止~

real work

由于这是我第一次还原论文的实现内容，一开始也是一头雾水，无从下手。慢慢多看了几遍公式，我发现了一些规律。

这里写图片描述

了解LDA的同学对上面的一系列定义应该不陌生。
这是一个我们要实现的所有的矩阵的描述。

从头开始，对每一个用户生成一个用户对应主题的矩阵（$x_u$）。每一个用户都有一个自己的$sigma_u^2$。但是，由于我们是在初步实验阶段，我就写了个文件，里面存储所有的$sigma^2$值并且都为0.1。

    # 制造假的sigma文件
    user_num = len(user)
    f = open('sigma_u_t.csv','w')
    f.write(',user')
    for i in range(split_num*10):
        f.write(',topic'+str(i))
    f.write('\n')
    for key in user.keys():
        f.write(','+key)
        for j in range(split_num*10):
            f.write(',0.1')
        f.write('\n')


    # 每一个用户的user-topic分布
    # sigma_u_t 是每个用户对于每一个topic的sigma值
    # 从文件里面读取每一个用户的每一个topic的sigma值
    # 每一行一个用户 (顺序就是下面生成的 user_ 中的顺序)
    user_sigma = pd.read_csv('sigma_u_t.csv')
    user_sigma = user_sigma.drop(['Unnamed: 0'],1)
    user_sigma.fillna(0.1)
    
            
    # 利用上面的用户对应评论的字典 make 一个 dataframe
    user_ = pd.DataFrame()
    temp1 = []
    temp2 = []
    for key in user.keys():
        for i in range(len(user[key])):
            temp1.append(key)
            temp2.append(user[key][i])
    user_['user'] = temp1
    user_['comment'] = temp2

【后期的做法可以是去通过标签统计分析每一个用户对于主题的概率再去确定每一个用户的$sigma^2$值】

下面我们需要实现$lambda_s = N(m_{pre_s},sigma_sI_K)$。这里就会发现我们要先去实现$m_{pre_s}$

这里写图片描述

可以从文中读到这个$m_{pre_s}$就是当前的shot的之前所有的shot对其的影响值相加，具体怎么影响的公式用到了 exponential decay 。（其中$Delta(s,s^{'})$是两个shot之间的绝对差，我这里的实现是用cos，也就是余弦值）

因为comment的$Delta$和shot的$Delta$计算方式相近，我就在下面的实现中一起实现了。

实现 $Delta$

先讲下余弦值计算的原理：

一个简单的例子（后面第三部分是相对复杂的例子）：

　　句子A：我喜欢看电视，不喜欢看电影。

　　句子B：我不喜欢看电视，也不喜欢看电影。

请问怎样才能计算上面两句话的相似程度？
基本思路是：如果这两句话的用词越相似，它们的内容就应该越相似。因此，可以从词频入手，计算它们的相似程度。

第一步，分词。

　　句子A：我/喜欢/看/电视，不/喜欢/看/电影。

　　句子B：我/不/喜欢/看/电视，也/不/喜欢/看/电影。

第二步，列出所有的词。

　　我，喜欢，看，电视，电影，不，也。

第三步，计算词频。

　　句子A：我 1，喜欢 2，看 2，电视 1，电影 1，不 1，也 0。

　　句子B：我 1，喜欢 2，看 2，电视 1，电影 1，不 2，也 1。

第四步，写出词频向量。

　　句子A：[1, 2, 2, 1, 1, 1, 0]

　　句子B：[1, 2, 2, 1, 1, 2, 1]

到这里，问题就变成了如何计算这两个向量的相似程度。

使用余弦这个公式，我们就可以得到，句子A与句子B的夹角的余弦。

这里写图片描述


余弦值越接近1，就表明夹角越接近0度，也就是两个向量越相似，这就叫余弦相似性。所以，上面的句子A和句子B是很相似的，事实上它们的夹角大约为20.3度。

从上面可以看出如果套用到我们的模型中，分成两类：

每个shot有一个词向量矩阵
每个comment有一个词向量矩阵

如下图。
左边是单词，右边是单词出现的次数。
map

因为我们是对一部视频也就是一篇文档处理，所以词向量可以是这篇文档里出现过的词。因此用下面的代码找出所有的词并存储成一个pandas的dataframe（使用这个的原因是处理起来很快而且方便）

# 统计关键词及个数 （根据文件）
def CountKey(fileName, resultName):
    try:
        #计算文件行数
        lineNums = len(open(fileName,'rU').readlines())
        # print u'文件行数: ' + str(lineNums)

        #统计格式 格式<Key:Value> <属性:出现个数>
        i = 0
        table = {}
        source = open(fileName,"r")
        result = open(resultName,"w")
        
        while i < lineNums:
            line = source.readline()
            line = line.rstrip()
            # print line
            
            words = line.split(" ")  #空格分隔
            # print str(words).decode('string_escape') #list显示中文
            
            #字典插入与赋值
            for word in words:
                if word!="" and table.has_key(word):      #如果存在次数加1
                    num = table[word]
                    table[word] = num + 1
                elif word!="":                            #否则初值为1
                    table[word] = 1
            i = i + 1

        #键值从大到小排序 函数原型：sorted(dic,value,reverse)
        dic = sorted(table.iteritems(), key = lambda asd:asd[1], reverse = True)
        word_fre = pd.DataFrame(dic)
        for i in range(len(dic)):
            #print 'key=%s, value=%s' % (dic[i][0],dic[i][1])
            result.write("<"+dic[i][0]+":"+str(dic[i][1])+">\n")
        return word_fre
        
    except Exception,e:    
        print 'Error:',e
    finally:
        source.close()
        result.close()
        # print 'END\n\n'


f = open('comments.txt','w')
for i in range(split_num):
    for x in shots[i]:
        for word in x:
            f.write(word.encode('utf-8') + ' ')
    f.write('\n')

word_fre = CountKey('comments.txt','comments_map')

最后得到的 word_fre 就是一个词频向量（全集），其实并不需要计算全集的词频。

    0    1
0    好    120
1    哈哈哈    116
2    哈哈哈哈    72
3    吴妈    50
4    卧槽    48
5    神父    48
6    人    41
7    黑社会    37
8    靓坤    35
9    真的    34
10    死    33
11    叻    31
12    说    30
13    君    28
14    一个    25
15    太    23
16    想    22
17    大佬    20
18    卖    20
19    吴    20
20    坤    20
21    香港    19
22    樽    19
23    爆    19
24    古惑仔    18
25    2333333    17
26    233333    17
27    笑    16
28    可爱    16
29    李丽珍    16
...    ...    ...
1986    额滴    1
1987    痛    1
1988    死于    1
1989    递纸    1
1990    hahahahhahahah8    1
1991    扭    1
1992    扑    1
1993    却    1
1994    扛    1
1995    阿公    1
1996    头子    1
1997    交个    1
1998    对手    1
1999    解构    1
2000    改一改    1
2001    惹不起    1
2002    湖地    1
2003    把持    1
2004    布吉岛    1
2005    傻仔    1
2006    莫名    1
2007    ′    1
2008    ‵    1
2009    陸仔    1
2010    兴趣    1
2011    祛湿    1
2012    君比靓    1
2013    培养    1
2014    不卡    1
2015    留学    1

我的构思是这样的。
构建每一个shot的词向量，就去统计每个shot里面的每个词的词频，没在该shot里出现过的但是在全集有的为0，词向量的顺序就和上面的 word_fre 一样，这样后面的计算直接就是处理两个dataframe就可以了。

同理，对每一个comment也是一样的。每一个comment都有一个词向量dataframe（这里会造成对内存的大量使用，但是计算起来快）

# 计算每一个shot里面的所有的单词的词频 ------->   缺点：执行速度实在太慢了，后期需要修改
result_s = []
for i in range(split_num):
    shot_word_fre = word_fre.copy()
    shot_word_fre['time'] = 0
    for x in shots[i]:
        for word in x:
            index = shot_word_fre[word_fre[0] == word.encode('utf-8')].index
            shot_word_fre['time'][index] = shot_word_fre['time'][index] + 1
    shot_word_fre = shot_word_fre.drop(1,1)
    result_s.append(shot_word_fre)
    
# 计算每一个comment的词频向量  -----------> 现在的办法是每个 comment 都有一个完整的词向量，便于后面的计算，问题是这样很占内存资源
# 按照每一个shot分片后内部的comment之间的delta计算
# result_c = []
# for i in range(split_num):
#     temp = []
#     for j in range(len(shots[i])):
#         shot_word_fre = word_fre.copy()
#         shot_word_fre['time'] = 0
#         for x in shots[i][j]:
#             for word in x:
#                 index = shot_word_fre[word_fre[0] == word.encode('utf-8')].index
#                 shot_word_fre['time'][index] = shot_word_fre['time'][index] + 1
#         shot_word_fre = shot_word_fre.drop(1,1)
#         temp.append(shot_word_fre)
#     result_c.append(temp)


# 计算每一个comment的词频向量  -----------> 现在的办法是每个 comment 都有一个完整的词向量，便于后面的计算，问题是这样很占内存资源
# 不按照每一个shot分片后内部的comment之间的delta计算，所有的comment进行计算
result_c = []
for i in range(split_num):
    for j in range(len(shots[i])):
        shot_word_fre = word_fre.copy()
        shot_word_fre['time'] = 0
        for x in shots[i][j]:
            for word in x:
                index = shot_word_fre[word_fre[0] == word.encode('utf-8')].index
                shot_word_fre['time'][index] = shot_word_fre['time'][index] + 1
        shot_word_fre = shot_word_fre.drop(1,1)
        result_c.append(shot_word_fre)

有了词向量之后就可以计算$Delta$了
我这里是将所有的$Delta$都计算出来了，做成了一个下三角矩阵，具体为什么是下三角矩阵仔细想想就知道了。这样做的好处是我后面使用$Delta(s,s^{'})$直接就可以变成delta_s(i,j)或是delta_c(i,j)。

p.s. 我做了修改，之前我理解的是每个shot里面的所有的comment之间计算$Delta$值，但是后来我想想不是这样的，对于comment应该还是所有的comment前后进行计算。因此，上面的result_c要改，这里的delta_c也要改，我把原先的代码注释掉了。

# 计算delta<s,_s> : 这里用的是词频向量 余弦值    -----> 下三角矩阵，后面方便计算
# 从后面的shot往前计算
delta_s = np.zeros((split_num,split_num))
seq = range(split_num)
# 修改 time 的数据类型 to float64
for shot in result_s:
    shot.time = shot.time.astype('float64')

seq.reverse()
for i in seq:
    for j in range(i):
        numerator = np.sum(result_s[i].time*result_s[j].time)
        denominator = pow(np.sum(pow(result_s[i].time,2)),0.5)*pow(np.sum(pow(result_s[j].time,2)),0.5)
        if denominator != 0:
            cos = numerator/denominator
        else:
            cos = 0
        delta_s[i][j] = cos
        

# 计算delta<c,_c> : 这里用的是词频向量 余弦值    -----> 下三角矩阵，后面方便计算
# 从后往前
# 这里是按照每个shot分开然后计算里面的comment
# seq = range(len(result_c))
# # 修改 time 的数据类型 to float64
# for i in seq:
#     for comment in result_c[i]:
#         comment.time = comment.time.astype('float64')

# # 创建每个shot的一个矩阵，用list存储
# delta_c = []
# for i in seq:
#     length = len(result_c[i])
#     delta_c_temp = np.zeros((length,length))
#     delta_c.append(delta_c_temp)

# for i in seq:
#     seq2 = range(len(result_c[i]))
#     seq2.reverse()
#     for j in seq2:
#         for k in range(j):
#             numerator = np.sum(result_c[i][j].time*result_c[i][k].time)
#             denominator = pow(np.sum(pow(result_c[i][j].time,2)),0.5)*pow(np.sum(pow(result_c[i][i].time,2)),0.5)
#             if denominator != 0:
#                 cos = numerator/denominator
#             else:
#                 cos = 0
#             delta_c[i][j][k] = cos
            
    
# 计算delta<c,_c> : 这里用的是词频向量 余弦值    -----> 下三角矩阵，后面方便计算
# 从后往前
# 这里是不按照每个shot分开然后计算里面的comment
seq = range(len(result_c))
# 修改 time 的数据类型 to float64
for i in seq:
    result_c[i].time = result_c[i].time.astype('float64')

# list存储
delta_c = np.zeros((len(result_c),len(result_c)))

for i in seq:
    for k in range(i):
        numerator = np.sum(result_c[i].time*result_c[k].time)
        denominator = pow(np.sum(pow(result_c[i].time,2)),0.5)*pow(np.sum(pow(result_c[j].time,2)),0.5)
        if denominator != 0:
            cos = numerator/denominator
        else:
            cos = 0
        delta_c[i][k] = cos

由于第一个shot没有在它之前的shot，所以第一个shot的$m_{pre_s}$等于零。
接下来的每个$m_{pre_s}$都与之前的有关，并且是针对每一个topic的，所以$m_{pre_s}$应该是一个矩阵（这点想清楚才能编程）

# 有了上面的矩阵后，计算论文中提到的 M_pre_s 以及 M_pre_c
# 需要两个衰减参数 gamma_s 以及 gamma_c
# M_pre_s 比较好计算，M_pre_c 比较复杂一点，因为涉及到了每一个shot
gamma_s = 0.5 # 我自己设的
gamma_c = 0.3 # 论文中做实验得到的最好的值
M_pre_s = np.zeros((split_num,total_topic))  # 行：shot个数    列：topic个数
lambda_s = np.zeros((split_num,total_topic))
sigma_s = 0.1 # 应该是每个片段的都不一样，但是这里我认为其实每个片段的topic分布没有统计可能性，不合理，都设成一样的了

# 先初始化 M_pre_s[0] 以及 lambda_s[0]
mu = 0 # 初始的 M_pre_s[0] 都是0
s = np.random.normal(mu,sigma_s,total_topic) # 不知道这个做法对不对，用正态生成x坐标，再用pdf去生成y值
lambda_s[0] = st.norm(mu, sigma_s).pdf(s)

# 从 第1的开始
for i in range(1,split_num):
    for topic in range(total_topic): # 先循环topic
        for j in range(i):
            numerator = np.exp(-gamma_s*delta_s[i][j])*lambda_s[j][topic]
            denominator = np.exp(-gamma_s*delta_s[i][j])
        M_pre_s[i][topic] = numerator/denominator
        s = np.random.normal(M_pre_s[i][topic],sigma_s,1)
        lambda_s[i][topic] = st.norm(M_pre_s[i][topic], sigma_s).pdf(s)

需要提一句，我里面可能会有些变量没定义就使用了，毕竟这是我的一个心路历程的总结，不是完整的源代码，如果需要看源代码可以去我的 Github 上看。

接下来就是计算 $m_{pre_c}$和$pi_c$了，处理起来会比较复杂一点，因为里面涉及了评论的用户以及用户对应的topic分布。这时候如果只是匹配的话程序会慢到死的，我的做法就是先处理出一张大表（dataframe）之后，每条评论以及对应的user以及对应的topic分布就可以很轻松快速地查到了。

# 总的topic个数，我在这里才填了total_topic这个参数，是有点晚了，不过，我之前在这里还遇到了一些问题，我以为是每个shot里面有固定的topic数，然后总的topic数是相乘的结果，后来经过一番认真思考，我才悔悟到原LDA中的topic数是固定的，然后不管你输入了多少文档，这个也应该一样，只不过文档变成了shot。
total_topic = 10


# 每一个用户的user-topic分布
# sigma_u_t 是每个用户对于每一个topic的sigma值
# 从文件里面读取每一个用户的每一个topic的sigma值
# 每一行一个用户 (顺序就是下面生成的 user_ 中的顺序)
user_sigma = pd.read_csv('sigma_u_t.csv')
user_sigma = user_sigma.drop(['Unnamed: 0'],1)
user_sigma.fillna(0.1)

        
# 利用上面的用户对应评论的字典 make 一个 dataframe
user_ = pd.DataFrame()
temp1 = []
temp2 = []
for key in user.keys():
    for i in range(len(user[key])):
        temp1.append(key)
        temp2.append(user[key][i])
user_['user'] = temp1
user_['comment'] = temp2


# 处理得到一个大表，里面包括所有评论以及评论的人，和每个人对应的所有的topic的sigma值
# 这里处理之后好像有点问题，有些用户没有，下面我直接就都填充0.1了
comment_per_shot = []
for i in range(split_num):
    temp = pd.DataFrame(com[i])
    u = []
    tem = pd.DataFrame()
    for j in range(len(temp)):
        user_id = user_[user_.comment == temp[0][j]].iloc[0][0]
        u.append(user_id)
        a = user_sigma[user_sigma.user == user_id].iloc[:,1:]
        tem = [tem,a]
        tem = pd.concat(tem)
    tem = tem.reset_index().drop(['index'],1)
    temp['user'] = pd.DataFrame(u)
    temp = temp.join(tem)
    comment_per_shot.append(temp)

# 所有的 comment 的一个 dataframe ,comment-user_id-topic0,1,2...99 ，后面的topic分布是user_id的
comment_all = pd.concat(comment_per_shot).reset_index().drop('index',1)
# 给那些没有topic分布的用户填充0.1 ----> 缺失值（就是生成用户的topic分布表没有生成全）
comment_all = comment_all.fillna(0.1) # 没有topic分布的都填充为0.1
comment_all = comment_all.rename(columns={0:'comment'})

上面的 comment_all 的结构基本上就是

index - comment - user - user's topic distribution(列数是总的topic个数)

这结构下面我还会更新用作更方便的计算。
然后有个这个 dataframe 之后我们就可以计算 $m_{pre_c}$和$pi_c$

# 生成 pi_c 和 M_pre_c 不同于上面，因为这里是对每个shot的面的comment进行操作
# 先初始化 M_pre_c[0] 和 第0个 shot 里的第一个 comment 对应的 pi_c[0]
M_pre_c = np.zeros((len(comment_all),total_topic))  # 行：shot个数    列：topic个数
pi_c = np.zeros((len(comment_all),total_topic))
for i in range(total_topic):
    pi_c[0][i] = lambda_s[0][i]*comment_all.iloc[0][i+2] + M_pre_c[0][i]

start = 0 # shot 之间的位移
for q in range(split_num):
    if q == 0:
        for i in range(1,len(com[q])):
            for topic in range(total_topic): # 先循环topic
                numerator = 0
                denominator = 0
                for j in range(i):
                    numerator += np.exp(-gamma_c*delta_c[i][j])*pi_c[j][topic]
                    denominator += np.exp(-gamma_c*delta_c[i][j])
                M_pre_c[i][topic] = numerator/denominator
                pi_c[i][topic] = lambda_s[q][topic]*comment_all.iloc[i][topic+2] + M_pre_c[i][topic]
        start += len(com[q])
    else:
        for i in range(start,start+len(com[q])):
            for topic in range(total_topic): # 先循环topic
                numerator = 0
                denominator = 0
                for j in range(i):
                    numerator += np.exp(-gamma_c*delta_c[i][j])*pi_c[j][topic]
                    denominator += np.exp(-gamma_c*delta_c[i][j])
                M_pre_c[i][topic] = numerator/denominator
                pi_c[i][topic] = lambda_s[q][topic]*comment_all.iloc[i][topic+2] + M_pre_c[i][topic]
        start += len(com[q])

基本任务上面就完成了，下面的是两个更新公式的实现（下面两幅图），看上去只有两个，但是实现起来特别困难，而且还要考虑时间复杂的问题。

这里写图片描述

在我第一次实现出来的版本需要两个多小时的执行时间（-_-|），后来进行了 dataframe的更新以及采用多个线程池的方式提高了运行的速度。

最后的最后，我终于知道为什么了，我把topic设置的太多了，我认为一部电影分片10片，每一片里面有10个topic，这样一部电影里面就有了100个topic，计算起来时间很久，这段心路历程在上面的total_topic这个变量定义的地方我有详述了，所以在优化之后，我修改了topic的总数

下面我直接贴出我优化后的版本，再去讲怎么优化的。

下面代码的前半段是生成一些需要的矩阵，当然，我这里做成了 dataframe。

注意：下面的代码很长很长，但是，我做了很多的注释，我相信能够解释的清楚。

里面涉及了 GibbsLDApy （我模仿 GibbsLDA++ 实现的python版本）的内容。大家也可以去看看我实现的这个版本 GibbsLDApy，顺便点点赞：）。后面我会整合所有代码形成新的 danmuLDA 做成分支。

# 生成 trndocs.dat 文件
# 该文件就是视频的剪切 -----> 分成了 split_num 份数，每一份代表一篇文档
f = open('test_data/trndocs.dat','w')
f.write(str(split_num)+'\n')
for i in range(split_num):
    for j in range(len(shots[i])):
        for k in range(len(shots[i][j])):
            f.write(shots[i][j][k].encode('utf-8')+' ')
    f.write('\n')


import time # 用来记录代码执行时间
# 欧拉函数的定义
eur = eulerlib.numtheory.Divisors(10000) # maxnum


#　执行 model 初始化
#  因为现在还是实验阶段，我没有和原LDA整合成一个完整的LDA，所以我这里用了 GibbsLDApy的初始化model的功能
argv = ['-est', '-alpha', '0.5', '-beta', '0.1', '-ntopics', '100', '-niters',
        '1000', '-savestep', '100', '-twords', '20', '-dfile', 'trndocs.dat', '-dir', 'test_data/',
        '-model', 'trnmodel']
pmodel = Model()
pmodel.init(len(argv),argv)


# 将 comment_all 升级成一个新的大表 comment_all_sort 结构为 {comment,user_id,user_id的topic,该comment属于的shot的topic分布},有了这个表，后面的处理会很方便
a1 = pd.concat([comment_all,pd.DataFrame(M_pre_c)],axis=1)
temp = []
for i in range(split_num):
    for j in range(len(shots[i])):
        t = pd.DataFrame(lambda_s)[i:i+1]
        t['shot'] = i
        t['com'] = j
        temp.append(t)
a2 = pd.concat(temp)
a2 = a2.reset_index().drop('index',1)
comment_all_sort = pd.concat([a1,a2],axis=1)
comment_all_sort = comment_all.sort_values('user') # 按照 user 排序


# 生成 user-topic 分布的 dataframe
x_u_c_t = np.zeros((len(comment_all_sort),total_topic))
for i in range(len(comment_all_sort)):
    for topic in range(total_topic):
        s = np.random.normal(mu,comment_all_sort.icol(topic+2)[i],1)
        x_u_c_t[i][topic] = st.norm(mu, comment_all_sort.icol(topic+2)[i]).pdf(s)
user_id = comment_all_sort.drop_duplicates('user')['user'].reset_index().drop('index',1)
x_u_c_t = user_id.join(pd.DataFrame(x_u_c_t))

        
def lgt(y):
    return math.log(1+math.exp(y))
def dlgt(y):
    return 1/((1+math.exp(y))*np.log(10))

word2id = pd.read_csv('test_data/wordmap.txt',sep=' ') # 读取单词对应id的表
column = list(word2id)[0]   # 这个是因为第一行是单词的个数，会变成index，下面转换成字典后出现二级索引，所以做了处理
word2id = word2id.to_dict()
yita_lambda_s = lambda_s.copy()

# 线程函数 --> 计算 yita_lambda_s
def calculate_lambda_s(shot,start):
    for topic in range(total_topic):
        result = 0
        lam_s = lambda_s[shot][topic]
        for comment in range(len(shots[shot])):
            u = x_u_c_t[x_u_c_t.user == comment_all.iloc[comment+start][topic+1]]
            x_u = u.iloc[0][topic+1]
            m_pre_c = M_pre_c[comment+start][topic]
            t1 = x_u*dlgt(x_u*lam_s+m_pre_c)
            t2 = 0
            for t in range(total_topic):
                t2 += lgt(comment_all.iloc[comment+start][t+2]*lam_s+M_pre_c[comment+start][t])
            t3 =t2
            t2 = eur.phi(t2)
            t3 = eur.phi(t3+len(shots[shot][comment]))
            n_tc = 0
            for word in shots[shot][comment]:
                word = word.encode('utf-8')
                if word != ' ' :
                    try:
                        num = word2id[column][word]
                        n_tc += pmodel.nw[num][topic]
                    except Exception,e:  
                        print Exception,":",e
            t4 = eur.phi(lgt(x_u*lam_s+m_pre_c) + n_tc)
            t5 = eur.phi(lgt(x_u*lam_s+m_pre_c))
        result += t1 * (t2 - t3 + t4 - t5)
        yita_lambda_s[shot][topic] = -(lam_s+M_pre_s[shot][topic])/(lam_s*lam_s) + result
        
        


# 定义容量比视频片段一样多一些的线程池
pool = threadpool.ThreadPool(split_num+2)     

start_time = time.time() # 下面的多线程开始执行的时间
start = 0 # 初始化，用于控制在哪一个shot里面
for shot in range(len(shots)):
    lst_vars = [shot,start]
    func_var = [(lst_vars, None)]
    start += len(shots[shot]) # start 增加位移，移动一个shot
    requests = threadpool.makeRequests(calculate_lambda_s, func_var)
    [pool.putRequest(req) for req in requests]
pool.wait()
print 'updating lambda_s %d second'% (time.time()-start_time)



# 定义容量为 total_topic 的一半
pool_cal = threadpool.ThreadPool(total_topic/2)  

# x_u_c_t 的更新代码
# 注意 ：这里的 comment_all 已经排过序了，和上面的不一样
def calculate_x_u_c_t(i,start):
    for topic in range(total_topic):
        result = 0
        for j in range(start,start+user_ct.iloc[i]):
            lambda_s_t = comment_all_sort.iloc[j,topic+total_topic+total_topic+2]
            m_pre_c_t = comment_all_sort.iloc[j,topic+total_topic+2]
            x_u = x_u_c_t.iloc[j,topic+1]
            print(lambda_s_t)
            print(m_pre_c_t)
            print(x_u)
            t1 = lambda_s_t*dlgt(x_u*lambda_s_t + m_pre_c_t)
            t2 = []
            for t in range(total_topic):
                lst_vars = [comment_all_sort.iloc[j,t+2]*comment_all_sort.iloc[j,t+total_topic+total_topic+2]+comment_all_sort.iloc[j,t+total_topic+2],t2]
                func_var = [(lst_vars, None)]
                requests = threadpool.makeRequests(add_t2, func_var)
                [pool_cal.putRequest(req) for req in requests]
            pool_cal.wait()
            t2 = sum(t2)
            print(t2)
            t3 = eur.phi(t2+len(shots[int(comment_all_sort.ix[j,['shot']])][int(comment_all_sort.ix[j,['com']])]))
            t2 = eur.phi(t2)
            n_tc = 0
            for word in shots[int(comment_all_sort.ix[j,['shot']])][int(comment_all_sort.ix[j,['com']])]:
                word = word.encode('utf-8')
                if word != ' ' :
                    try:
                        num = word2id[column][word]
                        n_tc += pmodel.nw[num][topic]
                    except Exception,e:  
                        print Exception,":",e
            t4 = eur.phi(lgt(x_u*lambda_s_t + m_pre_c_t)+ n_tc)
            t5 = eur.phi(lgt(x_u*lambda_s_t + m_pre_c_t))
            result += t1 * (t2 - t3 + t4 - t5)
        x_u_c_t.iloc[j,topic+1] = x_u_c_t.iloc[j,topic+1] - yita*(-x_u/(comment_all_sort.iloc[j,topic+2]*comment_all_sort.iloc[j,topic+2]) + result)
        print(x_u_c_t.iloc[j,topic+1])

# 定义容量比用户数量十分之一多一些的线程池
pool = threadpool.ThreadPool(len(x_u_c_t)/10+2)     

user_ct = comment_all_sort.groupby('user').count()['topic0']
yita_x_u_c_t = x_u_c_t.copy()
yita = 0.3
start_time = time.time() # 下面的多线程开始执行的时间
start = 0 # 初始化，用于控制在哪一个shot里面
for i in range(len(user_ct)):
    lst_vars = [i,start]
    func_var = [(lst_vars, None)]
    start += user_ct.iloc[i] # start 增加位移，移动一个shot
    requests = threadpool.makeRequests(calculate_x_u_c_t, func_var)
    [pool.putRequest(req) for req in requests]
pool.wait()
print 'updating x_u_c_t %d second'% (time.time()-start_time)

到现在为止我们所需要的新功能就都实现啦，后面需要的就是遵循之前的伪代码图去阶段更新参数~

beta 版本先写到这，后面我还会补充的，代码整合过几天再做，服务器被老师关掉了 -_-||

paper done 2017/05/12

【论文实现】一篇Sigkdd的弹幕分析论文的python实现【LDA 实践者】