# Mann-Kendall change-point detection and plotting utilities.
#--------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False  # Prevent minus sign rendering issues

def Kendall_change_point_detection(index, inputdata, area_simple):
    """Run Mann-Kendall forward/backward tests and plot UF/UB curves."""
    inputdata = np.array(inputdata)
    n = inputdata.shape[0]

    # Forward sequence: cumulative rank statistic Sk and UFk
    Sk = [0]
    UFk = [0]
    s = 0
    Exp_value = [0]
    Var_value = [0]
    # Start at i=1 because UFk(0) is defined as 0 in the MK test
    for i in range(1, n):
        for j in range(i):
            if inputdata[i] > inputdata[j]:
                s = s + 1
            else:
                s = s + 0
        Sk.append(s)
        Exp_value.append((i + 1) * (i + 2) / 4)  # Expected value E[Sk]
        Var_value.append((i + 1) * i * (2 * (i + 1) + 5) / 72)  # Variance of Sk
        UFk.append((Sk[i] - Exp_value[i]) / np.sqrt(Var_value[i]))

    # Backward sequence: reverse the series and compute UBk
    Sk2 = [0]
    UBk = [0]
    UBk2 = [0]
    s2 = 0
    Exp_value2 = [0]
    Var_value2 = [0]
    inputdataT = list(reversed(inputdata))
  # Start at i=1 because UBk(1) is defined as 0 in the reverse test
    for i in range(1, n):
        for j in range(i):
            if inputdataT[i] > inputdataT[j]:
                s2 = s2 + 1
            else:
                s2 = s2 + 0
        Sk2.append(s2)
        Exp_value2.append((i + 1) * (i + 2) / 4)  # Expected value E[Sk2]
        Var_value2.append((i + 1) * i * (2 * (i + 1) + 5) / 72)  # Variance of Sk2
        UBk.append((Sk2[i] - Exp_value2[i]) / np.sqrt(Var_value2[i]))
        UBk2.append(-UBk[i])

    # Re-reverse UBk so it shares the same time axis as UFk
    UBkT = list(reversed(UBk2))
    diff = np.array(UFk) - np.array(UBkT)
    K = list()
    # Identify change points where UF and UB curves cross
    for k in range(1, n):
        if diff[k - 1] * diff[k] < 0:
            K.append(k)

    fig = plt.figure(figsize=(8, 4), dpi=300)
    ax = fig.add_subplot(111)
    bwith = 2  # Spine line width
    ax.spines['bottom'].set_linewidth(bwith)
    ax.spines['left'].set_linewidth(bwith)
    ax.spines['top'].set_linewidth(bwith)
    ax.spines['right'].set_linewidth(bwith)
    ax.spines['bottom'].set_color('black')
    ax.spines['left'].set_color('black')
    ax.spines['top'].set_color('black')
    ax.spines['right'].set_color('black')

    fontdict = {'family': 'Times New Roman', 'size': 16}
    plt.plot(range(n,), UFk, label='UF', color='red', marker='s')
    plt.plot(range(n,), UBkT, label='UB', color='blue', linestyle='--', marker='o')
    plt.ylabel('UFk-UBk')

    # Significance thresholds at +/-1.96 and zero reference line
    x_lim = plt.xlim()
    plt.plot(x_lim, [-1.96, -1.96], ':', color='black')
    plt.plot(x_lim, [0, 0], '--', color='black')
    plt.plot(x_lim, [+1.96, +1.96], ':', color='black')

    plt.xticks(range(n), index, fontdict=fontdict)
    plt.yticks(fontsize=11)
    labels = ax.get_xticklabels() + ax.get_yticklabels()
    [label.set_fontname('Times New Roman') for label in labels]
    plt.ylabel('Mann-Kendall test value', fontdict=fontdict)
    plt.xlabel('Year', fontdict=fontdict)
    plt.text(0.015, 0.97, area_simple, fontdict=fontdict, color='black',
             bbox=dict(facecolor='white', edgecolor='black'),
             transform=plt.gca().transAxes, verticalalignment='top', horizontalalignment='left')

    x_major_locator = MultipleLocator(5)
    ax.xaxis.set_major_locator(x_major_locator)

    plt.show()
    return K

index = [1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020]
input = [18.85047, 18.621725, 18.508102, 18.497377, 18.563808, 18.671936, 18.817886, 19.058601, 18.762302, 18.703215, 18.789518, 18.714691, 19.185438, 18.915833, 18.761465, 18.87459, 18.627039, 18.714598, 18.967249, 18.914642, 18.62189, 18.766188, 19.000858]

Kendall_change_point_detection(index, input, 'AE')
