Python性能优化实践—曼肯德尔趋势检验计算效率提升

曼—肯德尔算法可以进行序列趋势和突变点的检验，实践发现对趋势的检验有较高的准确性，对突变点的检验相对较差。本文参考相关资料通过Python实现了对序列的曼肯德尔检验并进行了运行优化，提升了生产中的应用价值。
本次实践主要使用了算法优化、函数缓存、循环遍历优化、多进程、numba优化。由于整个过程以numpy数组的数值计算为主进一步还可以使用向量化进行加速，或者使用numexpr对Z值计算部分进行优化，但是未实践。
本次优化前将使用第三方库的算法编写函数，并通过apply函数将函数应用到每一列，在我的电脑上每10000行上的计算时间为70秒以上。
方法一：优化循环后为40多秒可以运算完毕，在使用常量法优化然后利用三核多进程计算可以18秒运算完毕。
方法二：利用函数缓存法优化正态分布的分位点部分计算，然后使用numba加速，预热后运算可以达到4.5秒。

主要参考资料：
趋势检验方法（二）MK趋势检验_aaakirito的博客-CSDN博客_mk检验

程序的优化的主要步骤为发现性能瓶颈和对主要瓶颈进行优化，主要方法有算法优化、循环优化、计算优化和编译优化。

注意：由于算法的原因，下面例子中的numpy数组必须是一维行数组，即reshape(1,-1)。

一、程序性能分析

在jupyter中
粗分析：
使用%%timeit对cell进行分析，比较各种实现方法的效率
细分析：
使用line_profiler对函数逐行分析，找到耗时占比最大的行
%load_ext line_profiler 加载line_profiler功能
%lprun -f funcname funcname(argu) 对函数funcname进行分析,获取每行的运行时间

二、程序性能优化

1、算法优化

将使用第三方库的算法，改为手动优化过的算法
进行曼—肯德尔趋势检验可以使用现成的pymannkendall进行计算，这个库中包含了曼—肯德尔趋势检验和他的多种扩展，计算量也较为全面。使用较为方便但是运行较慢。
import pymannkendall as mk
def mk_ku(data):
result=mk.original_test(data)
return (result.trend,result.z,result.slope)
参考CSDN上aaakirito的博客_CSDN博客-ACM算法题,ACM简单题,python领域博主的方法进行改进，提高算法的运行效率，但是功能缩减仅能计算最常使用的trend、z和solope三种结果。

def mk_(y):  
    alpha=0.1 
    n = len(y) 

    # 计算S的值
    v = np.array(y[1:]).reshape((-1, 1))
    v = np.tile(v, len(y) - 1)
    for i in range(len(y) - 1):
        v[:, i] -= y[i]
        v[0:i, i] = 0
    v = np.sign(v)
    s = np.sum(v)

    # 判断x里面是否存在重复的数，输出唯一数队列unique_x,重复数数量队列tp
    unique_x, tp = np.unique(x, return_counts=True)
    g = len(unique_x)
 
    # 计算方差VAR(S)
    if n == g:  # 如果不存在重复点
        var_s = (n * (n - 1) * (2 * n + 5)) / 18
    else:
        var_s = (n * (n - 1) * (2 * n + 5) - np.sum(tp * (tp - 1) * (2 * tp + 5))) / 18
 
    # 计算z_value
    if n <= 10:  # n<=10属于特例
        z = s / (n * (n - 1) / 2)
    else:
        if s > 0:
            z = (s - 1) / np.sqrt(var_s)
        elif s < 0:
            z = (s + 1) / np.sqrt(var_s)
        else:
            z = 0
 
    # 计算p_value，可以选择性先对p_value进行验证
    #p = 2 * (1 - norm.cdf(abs(z)))
 
    # 计算Z(1-alpha/2)
    q=normal_ppf(1 - alpha / 2)
    h = abs(z) >q
    

    # trend趋势判断
    if (z < 0) and h:
        trend ="decreasing"
    elif (z > 0) and h:
        trend ="increasing"
    else:
        trend = "no trend"
        
    # slope     
    n=len(y)   
    t=np.arange(0,n)
    slope, intercept = np.polyfit(y, t, deg=1)
    
    z=round(z,2)
    slope=round(slope,2)
 
    return trend,z,slope

优化算法主要集中在两部分

第一部分：计算S值，趋势检验方法

（二）MK趋势检验_aaakirito的博客-CSDN博客_mk检验的作者aaakirito和pymannkendall库的作者都给出了基于numpy数组的优化计算方法。其中pymannkendall库的计算方法更为优秀

以下是pymannkendall库的方法，利用了全1数组的特性进行了符号化计算和矩阵想减的方法加速了运算。

def mk_score(x):
    s = 0
    n=len(x)
    demo = np.ones(n) 
    for k in range(n-1):
        s = s + np.sum(demo[k+1:n][x[k+1:n] > x[k]]) - np.sum(demo[k+1:n][x[k+1:n] < x[k]])        
    return s

aaakirito的思路为，由双循环改为通过复制np数组错位想减实现，大大提高了效率。作者也提出将循环缩减为log2（len(y)）的方法，但是由于增加了开销，在短序列的时候表现较差，测试时在序列达到10000个样本时依然较差。
# 计算S的值
    v = np.array(y[1:]).reshape((-1, 1))
    v = np.tile(v, len(y) - 1)
    for i in range(len(y) - 1):
        v[:, i] -= y[i]
        v[0:i, i] = 0
    v = np.sign(v)
    s = np.sum(v)

第二部分：优化scipy.stats 模块中标准正态分布的分位点函数norm.ppf()

经过测试该函数效率较低，耗时占总计算时间的40%多，通过在网上找了一个俄罗斯人的写的函数进行替代，效率大大提高。
但是这个函数中exit(‘R8_NORMAL_CDF_INVERSE - Fatal error!’) 语句因为使用了exit()函数导致无法序列化进而不能开启多进程，用return 代替。(jupyter中有此问题，pycharm中没有。)

def r8_huge ( ):
    #返回理论最大值，当p小于等于0或者大于等于1时
    value = 1.79769313486231571E+308
    return value

def r8poly_value ( n, a, x ):
    value = 0.0
    for i in range ( n - 1, -1, -1 ):
        value = value * x + a[i]
    return value

def normal_ppf(p):
    a = np.array([ \
        3.3871328727963666080, 1.3314166789178437745e+2, \
        1.9715909503065514427e+3, 1.3731693765509461125e+4, \
        4.5921953931549871457e+4, 6.7265770927008700853e+4, \
        3.3430575583588128105e+4, 2.5090809287301226727e+3])
    b = np.array([ \
        1.0, 4.2313330701600911252e+1, \
        6.8718700749205790830e+2, 5.3941960214247511077e+3, \
        2.1213794301586595867e+4, 3.9307895800092710610e+4, \
        2.8729085735721942674e+4, 5.2264952788528545610e+3])
    c = np.array([ \
        1.42343711074968357734, 4.63033784615654529590, \
        5.76949722146069140550, 3.64784832476320460504, \
        1.27045825245236838258, 2.41780725177450611770e-1, \
        2.27238449892691845833e-2, 7.74545014278341407640e-4])
    const1 = 0.180625
    const2 = 1.6
    d = np.array([ \
        1.0, 2.05319162663775882187, \
        1.