# -*- coding: utf-8 -*-

'''
描述：基于RFM及K-means聚类的客户价值度分析
时间：2019-2-19
程序开发环境：JointQuant
程序输入：sales.csv(包含用户ID，订单日期，订单金额等客户属性)
程序输出：RFM得分数据写本地文件rfm_score.csv
参考文献：
1 宋天龙. 《Python数据分析与数据化运营》
2 HarveyLau.《RFM-Clustering》 https://github.com/HarveyLau/RFM-Clustering
3 watermelon12138.《K-means聚类算法、Pandas绘制概率密度图和TSNE展示聚类结果》
https://blog.csdn.net/watermelon12138/article/details/86549474
4 Wilame Lima Vallantin.《Apply RFM principles to cluster customers with K-Means》
'''

# 导入库
import time  # 导入时间库
import numpy as np  # 导入numpy库
import pandas as pd  # 导入pandas库
from scipy import stats
from sklearn.cluster import KMeans
from mpl_toolkits.mplot3d import Axes3D

#-- 参数初始化 --#
# input,output文件指定
inputfile = 'sales.csv'  # 销量及其他属性数据
outputfile = 'rfm_score.csv' # 保存结果的文件名

# 权重
W = [0.30, 0.54, 0.16]

# 原始数据参数
dtypes = {'ORDERDATE': object, 'ORDERID': object, 'AMOUNTINFO': np.float32}  # 设置每列数据类型_
index_column = 'USERID'

### 数据读取
data = pd.read_csv(inputfile, dtype=dtypes, index_col=index_column)  # 读取数据

# 原始数据审查和校验
# 原始数据概览
print ('Data Overview:')
print (data.head(4))  # 打印原始数据前4条
print ('-' * 60)
print ('Data before initial:')
print (data.describe())  # 打印原始数据基本描述性信息
print ('-' * 60)
#print ('Unique data:')
#print (data.nunique())  # 
#print ('-' * 60)

### 数据初始化
# 1，缺失值处理
# 2，无效业务数据处理，比如订单金额<=1的订单
# 3，日期格式转换
#
def initial_data(raw_data):
    # 缺失值审查
    na_cols = raw_data.isnull().any(axis=0)  # 查看每一列是否具有缺失值
    print ('NA Cols:')
    print (na_cols)  # 查看具有缺失值的列
    print ('-' * 60)
    na_lines = raw_data.isnull().any(axis=1)  # 查看每一行是否具有缺失值
    print ('NA Recors:')
    print ('Total number of NA lines is: {0}'.format(na_lines.sum()))  # 查看具有缺失值的行总记录数
    print (raw_data[na_lines])  # 只查看具有缺失值的行信息
    print ('-' * 60)    
    # 缺失值处理
    data = raw_data.dropna()  # 丢弃带有缺失值的行记录
    data = data[data['AMOUNTINFO'] > 1]  # 丢弃订单金额<=1的记录
    data = data[data['AMOUNTINFO'] < 9000]  # 丢弃订单金额>=9000的记录    
    
    # 日期格式转换
    data['ORDERDATE'] = pd.to_datetime(data['ORDERDATE'], format='%Y-%m-%d')  # 将字符串转换为日期格
    
    return data

# 原始数据初始化
data = initial_data(data)
print ('Data after initial:')
print (data.describe())  # 打印原始数据基本描述性信息
print ('-' * 60)

print('Most recent invoice is from:')
print(data['ORDERDATE'].max())
print ('-' * 60)

# To calculate recency, we’ll look into the ‘ORDERDATE’. 
# Since the date of our last invoice was 2018–12–06, 
# we’ll consider it as the most recent one. 
# Then, we’ll subtract each day from the day after to calculate the other ‘recencies’.
deadline_date = pd.datetime(2018, 12, 5)  # 指定一个时间节点，用于计算其他时间与该时间的距离

def tran2rfm_pd(sales_data):
    # 数据转换
    recency_value = sales_data['ORDERDATE'].groupby(sales_data.index).max()  # 计算原始最近一次订单时间
    frequency_value = sales_data['ORDERDATE'].groupby(sales_data.index).count()  # 计算原始订单频率
    monetary_value = sales_data['AMOUNTINFO'].groupby(sales_data.index).sum()  # 计算原始订单总金额
    
    # 计算RFM得分
    # 分别计算R、F、M得分
    r_interval = (deadline_date - recency_value).dt.days  # 计算R间隔
    
    r_score = pd.cut(r_interval, 5, labels=[5, 4, 3, 2, 1])  # 计算R得分
    f_score = pd.cut(frequency_value, 5, labels=[1, 2, 3, 4, 5])  # 计算F得分
    m_score = pd.cut(monetary_value, 5, labels=[1, 2, 3, 4, 5])  # 计算M得分

    # R、F、M标签数据
    rfm_list = [r_score, f_score, m_score]  # 将r、f、m三个维度组成列表
    rfm_cols = ['R', 'F', 'M']  # 设置r、f、m三个维度列名
    rfm_pd = pd.DataFrame(np.array(rfm_list).transpose(), dtype=np.int32, columns=rfm_cols,
                          index=frequency_value.index)  # 建立r、f、m数据框
    # R、F、M实际数据
    rfm_pd['recency'] = r_interval
    rfm_pd['frequency'] = frequency_value
    rfm_pd['monetary'] = monetary_value

    return rfm_pd

# RFM DataFrame做成
rfm_pd = tran2rfm_pd(data)   

# 打印输出和保存结果
# 打印结果
print ('RFM Scores Overview:')
print (rfm_pd.head())  
print ('-' * 60)
print ('RFM Scores Describe:')
print (rfm_pd.describe())    
print ('-' * 60)

# 散点图展示
rfm_pd.plot.scatter(x='recency',y='frequency',c='monetary',cmap='viridis_r', s=50)

# 3D图展示
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(rfm_pd['recency'], rfm_pd['frequency'], rfm_pd['monetary'], s=50)
ax.set_xlabel('Recency')
ax.set_ylabel('Frequency')
ax.set_zlabel('Monetary')

# RFM总分计算
def compute(r,f,m):
    return r*W[0] + f*W[1] + m*W[2]

#------------------ 方法一：基于RFM区间标志 + 分类规则的分类  ------------------#

# 客户分类
H = [u'重要价值用户', #分类1
     u'重要发展用户',   #分类2
     u'重要保持用户',  #分类3
     u'重要挽留用户',  #分类4
     u'一般价值用户',  #分类5
     u'一般发展用户', #分类6
     u'一般保持用户', #分类7
     u'一般挽留用户'] #分类8

def set_rfm_level(rs):
    '''
    r_mean = rs['R'][100:-100].mean()        #剔除最高，前100 降低平均值
    f_mean = rs['F'][100:-1].mean()
    m_mean = rs['M'][100:-350].mean() #剔除最高，前350 降低平均值
    '''
    r_mean = rs['R'].mean() 
    f_mean = rs['F'].mean()
    m_mean = rs['M'].mean() 
   
    for index in rs.index:
        rs.loc[index, u'用户价值分数'] = compute(rs.loc[index]['R'],rs.loc[index]['F'],rs.loc[index]['M'])    
        if rs.loc[index]['R'] < r_mean:
            if rs.loc[index]['F'] > f_mean:
                if rs.loc[index]['M'] > m_mean:
                    rs.loc[index, u'用户分类'] = H[0]
                    rs.loc[index, u'用户分类id'] = 0
                else:
                    rs.loc[index, u'用户分类'] = H[1]
                    rs.loc[index, u'用户分类id'] = 1

            else:
                if rs.loc[index]['M'] > m_mean:
                    rs.loc[index, u'用户分类'] = H[2]
                    rs.loc[index, u'用户分类id'] = 2
                else:
                    rs.loc[index, u'用户分类'] = H[3]
                    rs.loc[index, u'用户分类id'] = 3
        else:
            if rs.loc[index]['F'] > f_mean:
                if rs.loc[index]['M'] > m_mean:
                    rs.loc[index, u'用户分类'] = H[4]
                    rs.loc[index, u'用户分类id'] = 4
                else:
                    rs.loc[index, u'用户分类'] = H[5]
                    rs.loc[index, u'用户分类id'] = 5
            else:
                if rs.loc[index]['M'] > m_mean:
                    rs.loc[index, u'用户分类'] = H[6]
                    rs.loc[index, u'用户分类id'] = 6
                else:
                    rs.loc[index, u'用户分类'] = H[7]
                    rs.loc[index, u'用户分类id'] = 7
    return rs

#
rs = set_rfm_level(rfm_pd)
print ('Fianl RFM Scores Overview:')    
print(rs.head(10))
print ('-' * 60)

# 3D图表示
fig = plt.figure(figsize=(15,10))
dx = fig.add_subplot(111, projection='3d')
# 根据分类种类，定义区分颜色
colors = ['blue', 'white', 'red', 'black' , 'cyan', 'magenta' , 'yellow', 'green']
# c+: cyan + + 表示形式
# y-: yellow + - 表现形式

for i in range(0,8):
    '''
    dx.scatter(rs[rs[u'用户分类id'] == i].recency, 
               rs[rs[u'用户分类id'] == i].frequency, 
               rs[rs[u'用户分类id'] == i].monetary, 
               c = colors[i], 
               label = '分类 ' + str(i+1), 
               s=50)
    '''
    dx.scatter(rs[rs[u'用户分类id'] == i].recency, 
               rs[rs[u'用户分类id'] == i].frequency, 
               rs[rs[u'用户分类id'] == i].monetary, 
               c = colors[i], 
               label = H[i],
               s=50)

#dx.set_title('Clusters of clients')
dx.set_xlabel('Recency')
dx.set_ylabel('Frequency')
dx.set_zlabel('Monetary')
dx.legend()

#------------------ 方法二：基于RFM实际数值 + K-means的分类  ------------------#

def winsorize(rfm_pd):
    # 离散值处理
    rfm_pd['recency'] = np.log10(rfm_pd['recency'])
    rfm_pd['frequency'] = np.log10(rfm_pd['frequency'])
    rfm_pd['monetary'] = np.log10(rfm_pd['monetary'])
    return rfm_pd

# 离散值处理
rfm_pd = winsorize(rfm_pd)
rfm_pd.plot.scatter(x='recency',y='frequency',c='monetary',cmap='viridis_r', s=50)

def standardlize(rfm_pd):
    # 标准化处理
    rfm_pd['recency'] = stats.zscore(rfm_pd['recency'])
    rfm_pd['frequency'] = stats.zscore(rfm_pd['frequency'])
    rfm_pd['monetary'] = stats.zscore(rfm_pd['monetary'])
    return rfm_pd

# 标准化处理
rfm_pd = standardlize(rfm_pd)
rfm_pd.plot.scatter(x='recency',y='frequency',c='monetary',cmap='viridis_r', s=50)

def get_wcss_kmeans(X_train):
    wcss = []
    for i in range(1,8):
        kmeans = KMeans(n_clusters=i, init='k-means++', random_state=0)
        kmeans.fit(X_train)
        wcss.append(kmeans.inertia_)
    return wcss

X = rfm_pd.iloc[:,3:] # 取RFM实际数据作为训练集

plt.plot(range(1,8), get_wcss_kmeans(X))
#plt.title('Elbow graph')
plt.xlabel('Cluster number')
plt.ylabel('WCSS')
plt.show()

# K-means聚类
# Elbow graph(wcss) shows that an optimal number of clusters would be ?. 
# This is the point where line starts to do a soft descent. 
k = 4
iterations = 500
rank = pd.Series([4, 3, 2, 1])    #根据k值，设定排名，从到小

# 客户分类
H = [ #分类1
     u'4星会员',   #分类2
     u'3星会员',  #分类3
     u'2星会员',  #分类4
     u'1星会员'] 

def get_kmeans(X_train):
    model = KMeans(n_clusters=k, n_jobs=-1, max_iter=iterations)  # 分为k类
    #model = KMeans(n_clusters=k)  # 分为k类
    model.fit(X_train)  # 开始聚类
    return model

# 
def set_rfm_level(data,model):
    rfm = pd.concat([data, pd.Series(model.labels_, index=data.index)], axis=1)  # 详细输出每个样本对应的类别
    rfm.columns = list(data.columns) + [u'聚类类别']  # 重命名表头
    return rfm

kmeans = get_kmeans(X)
rfm_level_pd = set_rfm_level(rfm_pd,kmeans)

# 散点图表示
#rfm_level_pd.plot.scatter(x='recency',y='frequency',c=kmeans.labels_,cmap='viridis_r', s=50)
#rfm_level_pd.plot.scatter(x='recency',y='frequency',c='聚类类别',cmap='viridis_r', s=50)

# 3D图表示
fig = plt.figure(figsize=(15,10))
dx = fig.add_subplot(111, projection='3d')
# 根据分类种类，定义区分颜色
colors = ['blue', 'yellow', 'green', 'red', 'black']

for i in range(0,k):
    dx.scatter(rfm_level_pd[rfm_level_pd[u'聚类类别'] == i].recency, 
               rfm_level_pd[rfm_level_pd[u'聚类类别'] == i].frequency, 
               rfm_level_pd[rfm_level_pd[u'聚类类别'] == i].monetary, 
               c = colors[i], 
               label = '类别 ' + str(i+1), 
               s=50)

dx.set_title('Clusters of clients')
dx.set_xlabel('Recency')
dx.set_ylabel('Frequency')
dx.set_zlabel('Monetary')
dx.legend()

def clear_k_means(model):
    r1 = pd.Series(model.labels_).value_counts()  # 统计各个类别的数目
    r2 = pd.DataFrame(model.cluster_centers_)  # 找出聚类中心
    #r3 = pd.Series(r2[0] * W[0] + r2[1] * W[1] + r2[2] * W[2])
    r3 = pd.Series(compute(r2[0],r2[1],r2[2]))    
    r = pd.concat([r2, r1, r3], axis=1)  # axis=1 横向连接（0是纵向），得到聚类中心对应的类别下的数目
    r.columns = [u'R_质心'] + [u'F_质心'] + [u'M_质心'] + [u'类别数目'] + [u'分数']  # 重命名表头
    r = r.sort_values(by=[u'分数'])
    s = pd.Series(pd.DataFrame(r).index)
    s1 = pd.concat([s,rank],axis=1)
    s1.columns = [u'聚类类别'] + [u'排名']
    # rs = pd.concat([r, r4],axis=1)
    # rs.columns = [u'R_质心'] + [u'F_质心'] + [u'M_质心'] + [u'类别数目'] + [u'分数'] + [u'排名']
    return r,s1

r,s = clear_k_means(kmeans)
print(r)
print ('-' * 60)
print(s)
print ('-' * 60)

# 设置排名
def init_clear_data(rs,s):
    rs['USERID'] = rs.index
    # 按照rs,s相同名字的列合并
    rs = pd.merge(rs,s,how='left')
    rs.sort_values(by='R')
    rs.index = rs['USERID']
    del rs['USERID']

    return rs

#rs = init_clear_data(data=rfm_level_pd,model=kmeans,s=s)
rs2 = init_clear_data(rs=rfm_level_pd,s=s)
print(rs2.head(10))
print ('-' * 60)

# 3D图表示
fig = plt.figure(figsize=(15,10))
dx = fig.add_subplot(111, projection='3d')
# 根据分类种类，定义区分颜色
colors = ['blue', 'yellow', 'green', 'red']

for i in range(0,k):
    dx.scatter(rfm_pd[rs2[u'排名'] == i+1].recency, 
           rfm_pd[rs2[u'排名'] == i+1].frequency, 
           rfm_pd[rs2[u'排名'] == i+1].monetary, 
           c = colors[i], 
           label = H[i],
           s=50)

#dx.set_title('Clusters of clients')
dx.set_xlabel('Recency')
dx.set_ylabel('Frequency')
dx.set_zlabel('Monetary')
dx.legend()

def density_plot(data): 
    import matplotlib.pyplot as plt 
    plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 
    plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 
    # subplots=True表示dataframe格式的数据中每一列绘制一幅子图 
    p = data.plot(kind='kde', linewidth=2, subplots=True, sharex=False) 
    # p[1]代表第1个子图 
    [p[i].set_ylabel(u'密度') for i in range(3)] 
    plt.legend() 
    return plt

# 取RFM实际数据
X2 = rfm_pd.iloc[:,3:]

pic_output = 'pd_' #概率密度图文件名前缀

for i in range(k):
    # rs[u'聚类类别'] == i返回的是index+boolen（注意：rs的index事先不能丢失）
    # X2[]返回的是只有true的行
    #density_plot(X2[rs[u'聚类类别'] == i]).show()
    density_plot(X2[rs2[u'排名']==i+1]).savefig(u'%s%s.png' %(pic_output, i))

#write_csv(rs_=rs)