随机森林回归算法（CART决策树）【numpy版】

1.实现思路

（1）选择CART回归树的训练集：从有N个样本的样本集中，随机有放回地抽取N个样本作为当前回归树的训练样本
（2）执行CART回归树生成算法：
（a）选择特征：从数据集的M个特征中选取m个特征作为树节点的划分特征
（b）寻找最佳切分特征m和切分点s，求解：

遍历m个特征，对于每个特征寻找一个切分点s，求得能使上述式子达到最小值的值对（m,s），其中R1、R2为切分后的样本数据。
（c）用（b）中求得的值对（m,s）划分该节点，并满足以下条件：

并求得当前节点的输出：

（d）重复步骤（a）-（c），直到满足停止条件。
（e）最终决策树将输入样本划分成L个区域，分别对应决策树上的L个节点。决策树的非叶子节点保存当前节点的最佳切分特征和切分点，叶子节点则保存节点内所有样本的标签的均值。
（3）随机森林将多棵决策树集成起来，在预测时，随机森林中的每一棵决策树都将测试样本作为输入，每一棵决策树都会返回一个预测结果。随机森林将这些结果的平均值作为这个测试样本的预测输出。

2.并行化思路

本实验中采用了三种并行化思路：

（1）对决策树的训练使用多进程并行

这个思路比较简单。由于每棵决策树所需要的训练集都是不同的，那么可以认为每棵决策树在训练的过程中所拥有的系统资源都是独立的，互不干扰，因此可以充分使用多进程的优势而不必担心出现死锁问题。
示例代码：

（2）对决策树的预测使用多进程并行

这里和训练决策树思路相似，因为每棵决策树进行预测输出时的运算过程都是独立的，数据处理互不干扰。因此我们可以在计算每棵决策树的预测输出时为每棵决策树使用一个进程进行处理，可以加快预测速度。

（3）对决策树节点的最佳切分特征和切分点的寻找使用numpy并行计算

通过实验发现，在训练单棵决策树时，最耗时的地方就在于寻找每个节点的最佳分割特征和分割阈值。
在一般的算法中，通常使用两个for循环进行遍历：遍历所有特征，对于每个特征遍历节点内的数据集子集，找到当前特征的最佳分割阈值，最后选择损失函数最小的特征及其对应的分割阈值作为最佳分割特征与分割阈值。这种方法在小型数据集上效率还能接受，但在大数据集下运行效率极低。

于是，针对这一情况，对上述流程做出改进。通过观察可以发现，对于同一特征的每个分割阈值来说，计算切分后的方差误差这一运算是相互独立的。基于这个结论，我们可以把待遍历的所有阈值组合成一个大矩阵，即一个的矩阵，为切分阈值的数量。而对每个分割阈值进行分割操作后的输出MSE可以看作是这个矩阵按列经过一系列运算后的结果。因此我们可以将这层循环等价为一个N维矩阵的按列运算。于是，我们可以使用numpy库的矩阵运算对这部分运算进行优化，充分发挥numpy库矩阵并行化运算的优势。改进后的运算方法如下：

对于待遍历的特征列表，依次选择一个特征，重复以下操作：
① 对数据集中的训练样本根据当前特征进行排序，得到递增的分割阈值序列。
② 将排序后的数据集的y标签进行一个累加操作，得到一个数列。数列中的第n个元素代表数据集中前n个数据样本的y标签的和。
③ 将①中得到的分割阈值序列去重，得到去重后的分割阈值序列。
④ 找到中的元素在中第一次出现的位置p（从0开始算），记录数列中对应位置的元素，得到数列。中的第p-1个元素代表当前数据集以为切割阈值时，特征值小于的样本的y标签的总和，样本数为。
⑤ 对于中的每个元素，进行如下操作：用数据集中所有样本的y标签的总和减去即可得到特征值大于等于的样本的y标签的总和，样本数为。
⑥ 将数据集中所有的样本的y标签进行平方运算，得到数列，重复②③④⑤操作得到和。
⑦ 对于s’中每个元素，利用方差公式即可分别求得以当前元素（切割阈值）切割数据集后得到的两个数据集的方差。其中：

⑧ 根据⑦中的运算结果，很容易求得每个元素（切割阈值）的MSE：

⑨ 最终，我们要找的最佳分割阈值即为MSE最小值所对应的中的元素。

3.代码实现

import random, os, time, json
from numpy import *
import numpy as np
import pandas as pd
import multiprocessing, threading

from sklearn.externals.joblib import Parallel, delayed

'''
Author: ZTao-z
Time: 2019/06/30
'''

'''
dataset为matrix，且最后一列为label

算法流程：
# 随机筛选N个样本
# 随机抽取M个特征
# 遍历N个样本上的特征m，找到最优分割阈值，分割成两堆
    记录此特征及阈值
    若满足分裂条件：
        生成此结点的左右子节点，左节点小于该阈值，右节点大于该阈值
        将数据分成两堆，分别保存到左右节点
    否则：
        终止分裂
'''

class CART_tree:
    def __init__(self):
        # data为样本序号
        self.min_leaf_size = 100
        self.varThres = 0

    # 划分数据集
    def datasetSplit(self, dataset, feaNum, thres):
        dataL = dataset[nonzero(dataset[:,feaNum] < thres)[0],:]
        dataR = dataset[nonzero(dataset[:,feaNum] >= thres)[0],:]
        return dataL, dataR

    # 计算总的方差
    def getAllVar(self, dataset):
        return var(dataset[:,-1]) * shape(dataset)[0]

    def findFeatureAndThresParallel(self, feature, dataset):
        # 样本子集根据特征进行排序
        m = shape(dataset)[0]
        dataset_t = dataset[np.lexsort(dataset[:,feature].T)]
        # 分割阈值列表
        thresList = np.array(dataset_t[0][:,feature].T.tolist()[0])
        # 累加y标签和y^2标签
        sum_List = np.array(np.cumsum(dataset_t[0][:,-1]).tolist()[0])
        sq_sum_List = np.array(np.cumsum(np.square(dataset_t[0][:,-1])).tolist()[0])
        # 获取标签总和
        sum = sum_List[-1]
        sq_sum = sq_sum_List[-1]
        # 分割阈值去重
        new_thresList, index = np.unique(thresList, return_index=True)
        # 计算分割后数据集的大小N1，N2
        left_size = index
        right_size = m - left_size
        # 计算N1的y总和以及y^2总和
        left_sum = sum_List[left_size-1]
        left_sq_sum = sq_sum_List[left_size-1]
        # 计算N2的y总和以及y^2总和
        right_sum = sum - left_sum
        right_sq_sum = sq_sum - left_sq_sum
        # 防止0做除数
        left_size[0] = 1
        # 计算N1和N2的方差
        var_left = left_sq_sum / left_size - np.square(left_sum / left_size)
        var_right = right_sq_sum / right_size - np.square(right_sum / right_size)
        # 计算MSE
        total_lost = var_left * left_size + var_right * right_size
        # 计算最小MSE及其对应分割阈值
        if shape(thresList)[0] <= 2 * self.min_leaf_size:
            return
        l = index >= self.min_leaf_size
        r = index < m - self.min_leaf_size
        listRange = nonzero(l&r)[0]
        if shape(listRange)[0] == 0:
            return
        index = np.argmin(total_lost[listRange], axis=0)
        if total_lost[listRange[0]+index] < self.lowError:
            self.lowError = total_lost[listRange[0]+index]
            self.bestFeature = feature
            self.bestThres = new_thresList[listRange[0]+index]

    def chooseBestFeature(self, dataset, featureList, max_depth):
        # 停止条件 1：标签相同
        if len(set(dataset[:,-1].T.tolist()[0])) == 1 or max_depth == 0:
            regLeaf = mean(dataset[:,-1])
            return None, regLeaf 

        # 停止条件 2：已完成所有标签分类
        if len(featureList) == 1:
            regLeaf = mean(dataset[:,-1])
            return None, regLeaf 

        m,n = shape(dataset)
        totalVar = self.getAllVar(dataset)
        self.bestFeature = -1
        self.bestThres = float('-inf')
        self.lowError = inf

        # 遍历剩余划分特征 i
        for feature in featureList:
            self.findFeatureAndThresParallel(feature, dataset)

        # 停止条件3：划分后方差更大，则取消划分
        if totalVar - self.lowError < self.varThres:
            return None, mean(dataset[:,-1])

        # 停止条件4：划分后数据集太小
        dataL, dataR = self.datasetSplit(dataset, self.bestFeature, self.bestThres)
        if shape(dataL)[0] < self.min_leaf_size or shape(dataR)[0] < self.min_leaf_size:
            return None, mean(dataset[:,-1])

        # 成功则返回最佳特征和最小方差
        return self.bestFeature, self.bestThres

    # dataset: 数据集, featureList: 随机特征
    def createTree(self, dataset, max_depth=100):
        m, n = shape(dataset)
        featureList = list(range(n-1))
        bestFeat, bestThres = self.chooseBestFeature(dataset, featureList, max_depth) #最耗时

        if bestFeat == None:
            return bestThres
        regTree = {}
        # 记录此特征及阈值
        regTree['spliteIndex'] = bestFeat
        regTree['spliteValue'] = bestThres
        # 划分数据集
        dataL,dataR = self.datasetSplit(dataset, bestFeat, bestThres)
        regTree['left'] = self.createTree(dataL, max_depth-1)
        regTree['right'] = self.createTree(dataR, max_depth-1)
        return regTree

    def isTree(self, tree):
        return type(tree).__name__=='dict'

    def predict(self, tree, testData):
        if not self.isTree(tree):
            return float(tree)
        if testData[tree['spliteIndex']] < tree['spliteValue']:
            if not self.isTree(tree['left']):
                return float(tree['left'])
            else:
                return self.predict(tree['left'], testData)
        else:
            if not self.isTree(tree['right']):
                return float(tree['right'])
            else:
                return self.predict(tree['right'], testData)

class RandomForest:
    def __init__(self, n):
        self.treeNum = n
        self.treeList = []
        self.ct = CART_tree()
        self.n_jobs = 1

    def fit(self, dataset, jobs=1):
        m, n = shape(dataset)
        self.n_jobs = jobs
        pool = multiprocessing.Pool(processes=self.n_jobs)
        for i in range(self.treeNum):
            # 有放回采样
            #data_t = np.random.choice(range(m), 2000000).tolist()
            data_t = list(range(2000000))
            random_dataset = dataset[data_t,:]
            tt = createTreeThread(random_dataset, i)
            self.treeList.append(pool.apply_async(tt.run))
        pool.close()
        pool.join()
        for treeNum in range(len(self.treeList)):
            self.treeList[treeNum] = self.treeList[treeNum].get()

    def writeToFile(self, tree):
        if not os.path.exists("./model"):
            os.mkdir("./model")
        with open("./model/tree1.json", "w") as f:
            json.dump(tree,f)

    def loadFromFile(self):
        for root, dirs, files in os.walk("./model"):
            for file in files:
                with open(file, "r") as f:
                    for line in f.readlines():
                        self.treeList.append(json.loads(line))

    def predict(self, testData):
        result = []
        for i in range(len(testData)):
            res = []
            for tree in self.treeList:
                res.append(self.ct.predict(tree, testData[i]))
            result.append(res)
            
        return np.matrix(result).mean(1).T.tolist()
    
    def predictParell(self, testData):
        result = []
        pool = multiprocessing.Pool(processes=self.n_jobs)
        for tree in self.treeList:
            tt = predictTreeThread(tree, testData)
            result.append(pool.apply_async(tt.run))
        pool.close()
        pool.join()
        for i in range(len(result)):
            result[i] = result[i].get()
        return np.matrix(result).mean(0).tolist()

    def saveModel(self):
        self.writeToFile(self.treeList)

class createTreeThread:
    def __init__(self, dataset, number=0):
        self.data = dataset
        self.ct = CART_tree()
        self.n = number

    def run(self):
        begin = time.time()
        self.tree = self.ct.createTree(self.data)
        end = time.time()
        print("Tree", self.n," Finish in :", end-begin)
        return self.tree

class predictTreeThread:
    def __init__(self, tree, testData):
        self.ct = CART_tree()
        self.tree = tree
        self.datas = testData
    
    def run(self):
        self.datas = np.matrix(self.datas)
        index = np.array(range(shape(self.datas)[0]))
        self.result = np.array([0.0] * shape(self.datas)[0])
        self.predict_v2(self.tree, index)
        return self.result.tolist()
        

    def predict_v2(self, tree, index):
        if not self.ct.isTree(tree):
            self.result[index] = float(tree)
            return
        temp_index_smaller = nonzero(self.datas[index, tree['spliteIndex']] < tree['spliteValue'])[0]
        temp_index_bigger = nonzero(self.datas[index, tree['spliteIndex']] >= tree['spliteValue'])[0]
        if shape(temp_index_smaller)[0] > 0:
            self.predict_v2(tree['left'], index[temp_index_smaller])
        if shape(temp_index_bigger)[0] > 0:
            self.predict_v2(tree['right'], index[temp_index_bigger])

def readDataset():
    trainSet = []
    labelSet = []
    for i in range(1, 6):
        trainData = pd.read_csv(os.path.join('data/', 'train{}.csv'.format(i)), header=None, \
            delimiter="\t", quoting=3)
        labelData = pd.read_csv(os.path.join('data/', 'label{}.csv'.format(i)), header=None, \
            delimiter="\t", quoting=3)
        for example in list(trainData[0]):
            cur_example = example.strip().split(',')
            fin_example = map(float, cur_example)
            trainSet.append(list(fin_example))
        for label in list(labelData[0]):
            labelSet.append(float(label))

    tS_matrix = np.matrix(trainSet)
    tL_matrix = np.matrix(labelSet)
    final_trainSet = np.insert(tS_matrix, 13, values=tL_matrix, axis=1)

    return final_trainSet

def readTestData():
    testSet = []
    for i in range(1, 7):
        testData = pd.read_csv(os.path.join('data/', 'test{}.csv'.format(i)), header=None, \
            delimiter="\t", quoting=3)
        for example in list(testData[0]):
            cur_example = example.strip().split(',')
            fin_example = map(float, cur_example)
            testSet.append(list(fin_example))

    return testSet

def R2Loss(y_test, y_true):
    return 1 - (np.sum(np.square(y_test-y_true))) / (np.sum(np.square(y_true-np.mean(y_true))))

if __name__=="__main__":
    print("read dataset")
    trainData = readDataset()
    
    print("begin generate forest")
    rf = RandomForest(1)
    rf.fit(trainData, jobs=1)
    
    
    print("begin predict")
    begin = time.time()
    result = rf.predictParell(trainData.tolist())
    end = time.time()
    print("Predict use:", end-begin,"second")
    
    print("R2 Loss:",R2Loss(np.array(result[0]), np.array(trainData[:,-1].T.tolist()[0])))
    '''
    begin = time.time()
    test = readTestData()
    result = rf.predictParell(test)
    end = time.time()
    print("Predict use:", end-begin,"second")
    '''
    '''
    r = list(range(1,len(test)+1))
    output = pd.DataFrame( data={"id": r, "Predicted":result[0]} )
    # Use pandas to write the comma-separated output file
    output.to_csv(os.path.join(os.path.dirname(__file__), 'data', 'result_my_1.csv'), index=False, quoting=3)

    rf.saveModel()
    '''