import numpy as np

# 定义矩阵A
A = np.array([[1, 2, 3, 5], [1/2, 1, 1/2, 2], [1/3, 2, 1, 2], [1/5, 1/2, 1/2, 1]])

n = A.shape[0]
eig_val, eig_vec = np.linalg.eig(A)
Max_eig = max(eig_val)

CI = (Max_eig - n) / (n-1)
RI = [0, 0.0001, 0.52, 0.89, 1.12, 1.26, 1.36, 1.41, 1.46, 1.49, 1.52, 1.54, 1.56, 1.58, 1.59]  

CR = CI / RI[n]

print('一致性指标CI=', CI)
print('一致性比例CR=', CR)

if CR < 0.10:
    print('因为CR<0.10，所以该判断矩阵A的一致性可以接受!')
else:
    print('注意：CR >= 0.10，因此该判断矩阵A需要进行修改!')
    
    
# 算术平均
ASum = np.sum(A, axis=0)
Stand_A = A / ASum
weights = np.mean(Stand_A, axis=1)
print(weights)

# 几何平均
n, _ = A.shape
prod_A = np.prod(A, axis=1)
prod_n_A = np.power(prod_A, 1/n)
re_prod_A = prod_n_A / np.sum(prod_n_A)
print(re_prod_A)

# 特征值法
eig_values, eig_vectors = np.linalg.eig(A)
max_index = np.argmax(eig_values)
max_vector = eig_vectors[:, max_index]
weights = max_vector / np.sum(max_vector)
print(weights)

import numpy as np 

def minTomax(maxx, x):
    x = list(x)
    ans = [[(maxx-e)] for e in x]
    return np.array(ans)

def midTomax(bestx, x):
    x = list(x)  
    h = [abs(e-bestx) for e in x] 
    M = max(h)
    if M == 0:
        M = 1  # 防止最大差值为0的情况
    ans = [[(1-e/M)] for e in h]  # 计算每个差值占最大差值的比例，并从1中减去，得到新指标值
    return np.array(ans)

def regTomax(lowx, highx, x):
    x = list(x)
    M = max(lowx-min(x), max(x)-highx)
    if M == 0:
        M = 1  # 防止最大距离为0的情况
    ans = []
    for i in range(len(x)):
        if x[i]<lowx:
            ans.append([(1-(lowx-x[i])/M)]) 
        elif x[i]>highx:
            ans.append([(1-(x[i]-highx)/M)])
        else:
            ans.append([1])
    return np.array(ans)

if __name__ == '__main__':
    
    print("类型矩阵：1:极大型，2：极小型，3：中间型，4：区间型")
    kind = [[1],[2],[3,165],[4,[90,100]]]

    A=[[9,10,175,120],[8,7,164,80],[6,3,157,90]]

    A = np.array(A)
    n, m = A.shape

    X = np.zeros(shape=(n, 1))

    for i in range(m):
        match kind[i][0]:
            case 1:
                v = np.array(A[:, i])
            case 2:
                maxA = max(A[:, i])
                v = minTomax(maxA, A[:, i])
            case 3:
                bestA = kind[i][1]
                v = midTomax(bestA, A[:, i])
            case 4:
                lowA = kind[i][1][0]
                highA = kind[i][1][1]
                v = regTomax(lowA, highA, A[:, i])
        if i==0:
            X = v.reshape(-1, 1)
        else:
            X = np.hstack([X, v.reshape(-1, 1)])
    print("统一指标后矩阵为：\n{}".format(X))  # 打印处理后的矩阵X


    X = X.astype('float')
    for j in range(m):
        X[:, j] = X[:, j]/np.sqrt(sum(X[:, j]**2))  # 欧几里得范数


    x_max = np.max(X, axis=0)
    x_min = np.min(X, axis=0)
    d_z = np.sqrt(np.sum(np.square((X - np.tile(x_max, (n, 1)))), axis=1))
    d_f = np.sqrt(np.sum(np.square((X - np.tile(x_min, (n, 1)))), axis=1))

    s = d_f/(d_z+d_f)
    Score = 100*s/sum(s)
    for i in range(len(Score)):
        print(f"score {i+1}: {Score[i]:.2f}%")

import numpy as np

def entrop(p):
    n = len(p)
    lnp = np.zeros(n)
    for i in range(n):
        if p[i] != 0:
            lnp[i] = p[i]*np.log(p[i])
    return lnp

if __name__ == "__main__":
    X = np.array([[9, 0, 0, 0], [8, 3, 0.9, 0.5], [6, 7, 0.2, 1]])

    Z = X / np.sqrt(np.sum(X*X, axis=0))
    
    n, m = Z.shape
    D = np.zeros(m)

    for i in range(m):
        x = Z[:, i]
        p = x / np.sum(x)
        e = -np.sum(entrop(p)) / np.log(n)
        D[i] = 1 - e

    W = D / np.sum(D)
    print(f"{W=}")

import numpy as np

A = np.array([[55, 24, 10], [65, 38, 22], [75, 40, 18], [100, 50, 20]])
A_norm = A / np.mean(A, axis=0)

Y = A_norm[:, 0].reshape(-1, 1)
X = A_norm[:, 1:]

absX0_Xi = np.abs(X - np.tile(Y, (1, X.shape[1])))


a = np.min(absX0_Xi)
b = np.max(absX0_Xi)

rho = 0.5

gamma = (a + rho * b) / (absX0_Xi + rho * b)

print(f'各指标灰色关联度为：{np.mean(gamma, axis=0)}')

import numpy as np
import matplotlib.pyplot as plt


plt.rcParams['font.sans-serif'] = ['Kaiti']
plt.rcParams['axes.unicode_minus'] = False


def gm11(x0, predict_num):
    """
    传统 GM(1,1) 模型
    参数：
        x0 : (n,) 原始数据（1D 可迭代/ndarray）
        predict_num : 向后预测期数（int）
    返回：
        result : (predict_num,) 未来预测值
        x0_hat : (n,) 对原始数据的拟合值
        relative_residuals : (n-1,) 相对残差（从第2项开始）
        eta : (n-1,) 级比偏差
    """
    x0 = np.asarray(x0, dtype=float).reshape(-1)
    n = x0.size
    if n < 2:
        raise ValueError("x0 的长度至少为 2")


    x1 = np.cumsum(x0)
    z1 = 0.5 * (x1[:-1] + x1[1:])

    #    y = kx + b 且 k = -a
    y = x0[1:]
    x = z1
    nn = n - 1

    denom = (nn * np.sum(x * x) - np.sum(x) ** 2)
    if abs(denom) < 1e-12:
        raise ZeroDivisionError("最小二乘拟合时分母过小，可能数据退化。")

    k = (nn * np.sum(x * y) - np.sum(x) * np.sum(y)) / denom
    b = (np.sum(x * x) * np.sum(y) - np.sum(x) * np.sum(x * y)) / denom
    a = -k

    x0_hat = np.zeros(n)
    x0_hat[0] = x0[0]

    for m in range(1, n):
        x0_hat[m] = (1 - np.exp(a)) * (x0[0] - b / a) * np.exp(-a * (m - 1))


    result = np.zeros(predict_num)
    for i in range(predict_num):
        m = n + i 
        result[i] = (1 - np.exp(a)) * (x0[0] - b / a) * np.exp(-a * (m - 1))


    absolute_residuals = x0[1:] - x0_hat[1:]
    relative_residuals = np.abs(absolute_residuals) / x0[1:]


    class_ratio = x0[1:] / x0[:-1]

    denom_eta = (1 + 0.5 * a)
    if abs(denom_eta) < 1e-12:
        denom_eta = np.sign(denom_eta) * 1e-12
    eta = np.abs(1 - (1 - 0.5 * a) / denom_eta * (1.0 / class_ratio))

    return result, x0_hat, relative_residuals, eta

def check_input(x0):

    SMOOTH_RATIO_TH1 = 0.60 
    SMOOTH_RATIO_TH2 = 0.90 

    print('准指数规律检验')
    x1 = np.cumsum(x0)
    rho = x0[1:] / x1[:-1]  # 光滑度

    plt.figure(2)
    plt.plot(year[1:], rho, 'o-', label='光滑度')
    plt.plot([year[1], year[-1]], [0.5, 0.5], '-', label='临界线 0.5')
    plt.grid(True)
    plt.xticks(year[1:])
    plt.xlabel('年份'); plt.ylabel('原始数据的光滑度')
    plt.legend(loc='best')

    ind1 = 100 * np.mean(rho < 0.5)                     # 指标1：整体占比
    ind2 = 100 * np.mean(rho[2:] < 0.5) if len(x0) - 3 > 0 else np.nan  # 指标2：除前两期
    print(f'指标1：光滑比小于 0.5 的数据占比为 {ind1:.1f}%')
    if not np.isnan(ind2):
        print(f'指标2：除去前两个时期外，光滑比小于 0.5 的数据占比为 {ind2:.1f}%')


    passed = (ind1 >= 100 * SMOOTH_RATIO_TH1) and (np.isnan(ind2) or ind2 >= 100 * SMOOTH_RATIO_TH2)

    return passed


def train_val(x0,year):
    test_num = 3 if len(x0) > 7 else 2
    train_x0 = x0[:-test_num]
    test_x0 = x0[-test_num:]
    test_year = year[-test_num:]

    result1, _, _, _ = gm11(train_x0, test_num)

    # 试验组对比图
    plt.figure(3)
    plt.plot(test_year, test_x0, 'o-', label='试验组真实数据')
    plt.plot(test_year, result1, '*-', label='GM(1,1) 预测结果')
    plt.grid(True)
    plt.xticks(np.arange(test_year[0], test_year[-1] + 1))
    plt.xlabel('x'); plt.ylabel('y')
    plt.legend(loc='best')
    plt.title('试验组：真实 vs 预测')

    plt.show()


def test(x0,year):

    predict_num = 5
    result, x0_hat, relative_residuals, eta = gm11(x0, predict_num)


    # 相对残差与级比偏差图
    plt.figure(4)
    plt.subplot(2, 1, 1)
    plt.plot(year[1:], relative_residuals, '*-')
    plt.grid(True); plt.legend(['相对残差'], loc='best'); plt.xlabel('年份')
    plt.xticks(year[1:])

    plt.subplot(2, 1, 2)
    plt.plot(year[1:], eta, 'o-')
    plt.grid(True); plt.legend(['级比偏差'], loc='best'); plt.xlabel('年份')
    plt.xticks(year[1:])


    # 残差检验
    average_relative_residuals = float(np.mean(relative_residuals))
    print(f'平均相对残差为 {average_relative_residuals:.2f}')
    if average_relative_residuals < 0.1:
        print('模型拟合程度非常不错')
    elif average_relative_residuals < 0.2:
        print('模型拟合程度达到一般要求')
    else:
        print('模型拟合程度不太好，建议使用其他模型')

    # 级比偏差检验
    average_eta = float(np.mean(eta))
    print(f'平均级比偏差为 {average_eta:2f}')
    if average_eta < 0.1:
        print('模型拟合程度非常不错')
    elif average_eta < 0.2:
        print('模型拟合程度达到一般要求')
    else:
        print('模型拟合程度不太好，建议使用其他模型')


    # 最终综合图
    plt.figure(5)
    plt.plot(year, x0, '-o', label='原始数据')
    plt.plot(year, x0_hat, '-*m', label='拟合数据')
    fut_year = np.arange(year[-1] + 1, year[-1] + 1 + predict_num)
    plt.plot(fut_year, result, '-*b', label='预测数据')
    # 连接最后一个已知点与首个预测点
    plt.plot([year[-1], fut_year[0]], [x0[-1], result[0]], '-*b')
    plt.grid(True)
    plt.xticks(np.arange(year[0], fut_year[-1] + 1))
    plt.xlabel('年份'); plt.ylabel('排污总量')
    plt.legend(loc='best')
    plt.title('GM(1,1) 拟合与预测')


    plt.show()
def main(x0,year):
    if np.any(x0 < 0):
        raise ValueError('输入数据不能小于0')
    if len(x0) <= 3:
        raise ValueError('输入数据长度不能小于3')
    if not check_input(x0):
        raise ValueError('输入数据不符合要求')

    train_val(x0,year)
    test(x0,year)

if __name__ == "__main__":
    year = np.arange(1995, 2005, 1, dtype=int)
    x0 = np.array([4,5,6,7,8,9,10,11,12,13], dtype=float)
    main(x0,year)