import pandas as pd
from sklearn.model_selection import train_test_split
from matplotlib import rcParams
from sklearn.preprocessing import MinMaxScaler
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score

config = {
            "font.family": 'serif',
            "font.size": 10,               # 相当于六号大小
            "mathtext.fontset": 'stix',    # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大
            "font.serif": ['SimSun'],      # 宋体
            'axes.unicode_minus': False    # 处理负号，即-号
         }
rcParams.update(config)

# 数据处理
pathFile = "F:\郑州市1991-2022气象数据(3降雨1积温).xlsx"
df = pd.read_excel(pathFile, sheet_name='Sheet1')
m1 = df[['年份1']]
m2 = df[['年份2']]
x = df[['播种—越冬期', '拔节—抽穗期', '开花—成熟期', '全生育期的积温']]
y = df[['小麦产量']]

# 截取训练数据
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=23)

# 归一化
x_scaler = MinMaxScaler(feature_range=(-1, 1))
y_scaler = MinMaxScaler(feature_range=(-1, 1))
x_train = x_scaler.fit_transform(x_train)
y_train = y_scaler.fit_transform(y_train)
x_test = x_scaler.fit_transform(x_test)
y_test = y_scaler.fit_transform(y_test)

# 对训练集、测试集数据进行转置，方便后续的矩阵乘法的操作
sample_in = x_train.T
sample_out = y_train.T
ceshi_sample_in = x_test.T
ceshi_sample_out = y_test.T

# 网络参数配置
max_epochs = 100000     # 最大迭代次数
learn_rate = 0.01      # 学习率
mse_final = 6.5e-5      # 迭代结束条件
sample_numble = x_train.shape[0]
input_number = 4        # 输入特征属性(输入层)
output_number = 1       # 目标属性(输出层)
hidden_unit_number = 9  # 隐含层神经元个数
l2_lambda = 0.01        # L2正则化参数

# 超参
# 8*4
w1 = 0.5 * np.random.rand(hidden_unit_number, input_number) - 0.1

# 8*1
b1 = 0.5 * np.random.rand(hidden_unit_number, 1) - 0.1

# 8*1
w2 = 0.5 * np.random.rand(output_number, hidden_unit_number) - 0.1

# 2*1
b2 = 0.5 * np.random.rand(output_number, 1) - 0.1

def sigmoid(z):
    return 1.0 / (1 + np.exp(-z))

mse_history = []
for i in range(max_epochs):
    # 前向传播
    hidden_out = sigmoid(np.dot(w1, sample_in) + b1)
    network_out = np.dot(w2, hidden_out) + b2
    # 误差
    err = sample_out - network_out
    mse = np.average(np.square(err)) + l2_lambda * (np.sum(np.square(w1)) + np.sum(np.square(w2))) / (2 * sample_numble)
    mse_history.append(mse)
    if mse < mse_final:
        break
    # 后向传播
    delta2 = -err
    delta1 = np.dot(w2.transpose(), delta2) * hidden_out * (1 - hidden_out)

    delta_w2 = np.dot(delta2, hidden_out.transpose()) + l2_lambda * w2 / sample_numble
    delta_b2 = np.dot(delta2, np.ones((sample_numble, 1)))

    delta_w1 = np.dot(delta1, sample_in.transpose()) + l2_lambda * w1 / sample_numble
    delta_b1 = np.dot(delta1, np.ones((sample_numble, 1)))

    w2 -= learn_rate * delta_w2
    b2 -= learn_rate * delta_b2
    w1 -= learn_rate * delta_w1
    b1 -= learn_rate * delta_b1

# 隐藏层输出
hidden_out = sigmoid((np.dot(w1, sample_in) + b1))

# 输出层输出
network_out = np.dot(w2, hidden_out) + b2

# 翻转获取实际值
network_out = y_scaler.inverse_transform(network_out.T)
sample_out = y_scaler.inverse_transform(y_train)

print(network_out)
print("R^2=", r2_score(sample_out, network_out))
print("MRE=", np.sum(abs(sample_out - network_out) / sample_out) / len(network_out))

plt.figure(figsize=(8, 7))

plt.plot(np.arange(1991, 2016, 1), network_out[:, 0], color='red', marker='o', linestyle='-', lw=2, label="预测值",
         markerfacecolor='white')
plt.plot(np.arange(1991, 2016, 1), sample_out[:, 0], color='blue', marker='x', linestyle='-', lw=2, label="实际值")
# plt.plot(np.arange(1991, 2013, 1), network_out[:, 0], color='tomato', alpha=0, marker='o', label="R²=0.99794")
# plt.plot(np.arange(1991, 2013, 1), sample_out[:, 0], color='darkviolet', alpha=0, marker='*', markersize=9,
#          label="MRE=0.31%")

plt.xlabel("年份\nYear")
plt.ylabel("单位面积产量kg/h㎡\nYield per unit area kg/h㎡")
plt.title("郑州麦区\nZhengzhou Wheat District", fontdict={'fontsize': 10, 'color': 'black'})
m_ticks = np.arange(1991, 2025, 1)
plt.xticks(m_ticks, rotation=90)

plt.xlim(1991, 2016)
plt.ylim(2500, 5500)

plt.legend()

# *********************测试*********************

# 隐藏层输出
ceshi_hidden_out = sigmoid((np.dot(w1, ceshi_sample_in) + b1))

# 输出层输出
ceshi_network_out = np.dot(w2, ceshi_hidden_out) + b2

# 翻转获取实际值
ceshi_network_out = y_scaler.inverse_transform(ceshi_network_out.T)
ceshi_sample_out = y_scaler.inverse_transform(y_test)
print(ceshi_network_out)
print("R^2=", r2_score(ceshi_sample_out, ceshi_network_out))
print("MRE=", np.sum(abs(ceshi_sample_out - ceshi_network_out) / ceshi_sample_out) / len(ceshi_network_out))

plt.figure(figsize=(8, 7))

plt.plot(np.arange(2016, 2023, 1), ceshi_network_out[:, 0], color='red', marker='o', linestyle='-', lw=2, label="预测值", markerfacecolor='white')
plt.plot(np.arange(2016, 2023, 1), ceshi_sample_out[:, 0], color='blue', marker='x', linestyle='-', lw=2, label="实际值")
# plt.plot(np.arange(2016, 2023, 1), ceshi_network_out[:, 0], color='tomato', alpha=0, marker='o', label="R²=0.99794")
# plt.plot(np.arange(2016, 2023, 1), ceshi_sample_out[:, 0], color='darkviolet', alpha=0, marker='*', markersize=9, label="MRE=0.31%")

plt.xlabel("年份\nYear")
plt.ylabel("单位面积产量kg/h㎡\nYield per unit area kg/h㎡")
plt.title("郑州麦区\nZhengzhou Wheat District", fontdict={'fontsize': 10, 'color': 'black'})
m_ticks = np.arange(1991, 2025, 1)
plt.xticks(m_ticks, rotation=90)

plt.xlim(2016, 2022)
plt.ylim(1500, 5500)

plt.legend()


# ********************************************************************
# 截取训练数据  x1  1991-,   y1  -2022
x1 = np.array(x[:25])
y1 = np.array(y[:25])
x2 = np.array(x[25:])
y2 = np.array(y[25:])

# 归一化
x_scaler = MinMaxScaler(feature_range=(-1, 1))
y_scaler = MinMaxScaler(feature_range=(-1, 1))
x1 = x_scaler.fit_transform(x1)
x2 = x_scaler.fit_transform(x2)
y1 = y_scaler.fit_transform(y1)
y2 = y_scaler.fit_transform(y2)

# 对训练集、测试集数据进行转置，方便后续的矩阵乘法的操作
sample_in1 = x1.T
sample_out1 = y1.T
sample_in2 = x2.T
sample_out2 = y2.T

# 隐藏层输出
hidden_out2 = sigmoid((np.dot(w1, sample_in2) + b1))
# 输出层输出
network_out2 = np.dot(w2, hidden_out2) + b2
# 翻转获取实际值
network_out2 = y_scaler.inverse_transform(network_out2.T)
sample_out2 = y_scaler.inverse_transform(y2)

print(network_out2)
print()
print(sample_out2)
print("R^2=", r2_score(sample_out2, network_out2))
print("MRE=", np.sum(abs(sample_out2 - network_out2) / sample_out2) / len(network_out2))

plt.figure(figsize=(8, 7))

plt.plot(np.arange(2016, 2023, 1), network_out2[:, 0], color='red', marker='o', linestyle='-', lw=2, label="预测值", markerfacecolor='white')
plt.plot(np.arange(2016, 2023, 1), sample_out2[:, 0], color='blue', marker='x', linestyle='-', lw=2, label="实际值")
# plt.plot(np.arange(2016, 2023, 1), ceshi_network_out[:, 0], color='tomato', alpha=0, marker='o', label="R²=0.99794")
# plt.plot(np.arange(2016, 2023, 1), ceshi_sample_out[:, 0], color='darkviolet', alpha=0, marker='*', markersize=9, label="MRE=0.31%")

plt.xlabel("年份\nYear")
plt.ylabel("单位面积产量kg/h㎡\nYield per unit area kg/h㎡")
plt.title("郑州麦区\nZhengzhou Wheat District", fontdict={'fontsize': 10, 'color': 'black'})
m_ticks = np.arange(1991, 2025, 1)
plt.xticks(m_ticks, rotation=90)

plt.xlim(2016, 2022)
plt.ylim(2500, 5500)

plt.legend()
plt.show()




将该代码加入遗传算法优化，输出最终结果