This repo stores the code during the 2023 Huawei Cup Graduate Mathematical Modeling Competition.
import numpy as np

# 超过33的为发生血肿
def up_relative33(data):
    if data[0]==0:
        print("data is zero, please check")
        return (np.array([], dtype=np.int64), )
    
    return np.where(((data - data[0])/data[0])>=0.33)

process_models/q1/q1a_check_patient.py
TheWangYang/2023_huawei_mathematical_modeling
import pandas as pd
import os
import numpy as np
from process_models.q1.q1a_check_patient import up_relative33
import matplotlib.pyplot as plt

# 设置数据的根路径
root_path = '/data1/wyy/projects/_competition/math_model_v2/data'

# 定义表1变量
excel_name1 = '1.xlsx'

# 定义表2变量
excel_name2 = '2.xlsx'

# 定义附表变量
excel_name3 = 'additional.xlsx'

# 读取原始Excel文件
table1 = pd.read_excel(os.path.join(root_path, excel_name1))
table2 = pd.read_excel(os.path.join(root_path, excel_name2))
table3 = pd.read_excel(os.path.join(root_path, excel_name3))

# 定义保存160个病人的数据的字典
patients_info = []

# 定义保存当前第i个病人信息的字典
info_dict = {}

# 遍历160个病人的数据
for i in range(160):
    # 创建保存病人信息的字典
    info_dict = {}
    
    # 得到三个表格中的每个病人第i行数据
    row1 = table1.loc[i]
    row2 = table2.loc[i]
    row3 = table3.loc[i]
    
    # 首次发病到首次检查时间点
    first_check_timedelta = row1[14]
    print("first_check_timedelta: {}".format(first_check_timedelta))
    
    info_dict['name'] = row2[0]
    
    # 检查号缓存列表
    info_dict['check_seq_buf'] = [row2[1]]
    # 得到首次检查对应的HM体积
    info_dict['HM_volum'] = [row2[2]]
    
    # 得到首次检查对应的时间点
    first_check_time = row3[2]
    
    # 得到首次检查时间点
    info_dict['check_time_buf'] = [first_check_timedelta]

    if not (row2[1] == row3[3] and row1[3]==row2[1]):
        print(f"id not equal, {row2[1]}, {row1[3]}")
    
    # 定义
    j=1
    
    print(i, first_check_timedelta)
    
    # 内层循环实现将病人的所有随访时间点距离首次发病的时间间隔等存储到info_dict中
    while 1+j*23 <len(row2):
        # 判断该病人是否存在随访点i + j * 23
        if not np.isnan(row2[1+j*23]):
            # 得到随访1的流水号
            info_dict['check_seq_buf'].append(int(row2[1+j*23]))
            
            # 得到随访1对应的HM体积
            info_dict['HM_volum'].append(row2[2+j*23])
            
            # 将随访1的时间点-首次检查时间点，换成小时，然后再加上发病到随访1的时间点
            info_dict['check_time_buf'].append((row3[2+2*j]-first_check_time).total_seconds()/3600 + first_check_timedelta)
            
            # 得到后续的随访结点，均按照上述进行处理
            j+=1
        else:
            break
    
    # info_dict['check_time_buf'] = np.array(info_dict['check_time_buf'])
    # 将其中的HM对应的数据转换为numpy数组的形式
    info_dict['HM_volum'] = np.array(info_dict['HM_volum'])
    # 将i编号的病人的信息增加到patients_info字典中
    patients_info.append(info_dict)

# 保存有血肿扩张风险的病人的列表
patient_sign = []
time_buf = []

# 枚举160个病人的信息
for i, info in enumerate(patients_info):
    # 定义
    has_patient = 0
    
    # 找到HM体积相对一首次检查增加超过6000的第一个结果存储到seq_buf中
    seq_buf = np.where((info['HM_volum']-info['HM_volum'][0])>6000)[0]
    
    # 找到HM体积相对一首次检查相对体积增加33%的索引值进行求与（因为是或）
    seq_buf = np.concatenate((seq_buf, up_relative33(info['HM_volum'])[0]))
    
    print("seq_buf: {}".format(seq_buf))
    
    # input()
    
    # 遍历seq_buf数组
    for seq in seq_buf:
        # print(info['check_time_buf'], seq-1)
        if info['check_time_buf'][seq] <= 48:
            has_patient = 1
            # print()
            # print(i+1, seq, info['HM_volum'][seq], info['HM_volum'][0], info['HM_volum'][seq] - info['HM_volum'][0], (info['HM_volum'][seq] - info['HM_volum'][0])/info['HM_volum'][0])
            # print(1)
            break
        
    # 已经小于48小时
    if has_patient:
        # 添加时间
        time_buf.append(info['check_time_buf'][seq])
        print("{:.2f}".format(info['check_time_buf'][seq]))
    else:
        time_buf.append(0)
        print(0)

# 绘图，横轴为病人编号，纵轴为发病时间
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.figure()

plt.plot(time_buf[:100])

# 添加标题和坐标轴标签
plt.title("发病时间")
plt.xlabel("病人编号(0-100)")
plt.ylabel("发病时间(h)")

# plt.show()
plt.savefig("/data1/wyy/projects/_competition/math_model_v2/data_process/output/time_of_patient.jpg")
    # print(has_patient)
    # patient_sign.append(has_patient)
process_models/q1/q1a_read_data.py
TheWangYang/2023_huawei_mathematical_modeling
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from process_models.q1.q1b_dataloader import get_train_valid_q1

from imblearn.over_sampling import SMOTE

import matplotlib.pyplot as plt
import numpy as np
import shap

# 生成示例数据
X_train, X_test, y_train, y_test = get_train_valid_q1()

sm = SMOTE(random_state=0)

X_train, y_train = sm.fit_resample(X_train, y_train)
# X = np.random.random((100, 15))
# Y = np.random.randint(0, 2, (100,))

print(X_train.shape)


# X_train, X_test, y_train, y_test = train_test_split(X, Y, train_size=0.8, test_size=0.2, random_state=0)
# 转换数据为PyTorch的Tensor类型
X_train = torch.Tensor(X_train)
print(X_train.shape)
y_train = torch.LongTensor(y_train)
X_test = torch.Tensor(X_test)
y_test = torch.LongTensor(y_test)

# 创建数据集和数据加载器

X_train = torch.cat([X_train, X_test], dim=0)
y_train = torch.cat([y_train, y_test], dim=0)
train_dataset = TensorDataset(X_train, y_train)
train_loader = DataLoader(train_dataset, batch_size=10, shuffle=False)

# 构建神经网络模型
input_n = X_train.shape[-1]
start_n = 64
middle_n = [64, 64, 32, 16, 2]

print(input_n)
torch.manual_seed(42)

# 构建简单的
class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.fc = nn.Sequential(
            nn.Linear(input_n, start_n),
            # nn.BatchNorm1d(100),
            nn.ReLU(),
        )
        
        for i, num in enumerate(middle_n):
            if i==0:
                self.fc.append(nn.Sequential(nn.Linear(start_n, num),
                                             nn.ReLU()))
            else:
                self.fc.append(nn.Sequential(nn.Linear(middle_n[i-1], num),
                                             nn.ReLU()))

        self.fc.append(nn.Softmax(dim=1))

    def forward(self, x):
        x = self.fc(x)
        return x

model = Net()
model.train()
if True:
    # model.load_state_dict(torch.load('model_weights_base_info.pth'))
    # model.load_state_dict(torch.load('model_weights.pth'))
    print('loaded')

# 加载模型参数
# model.load_state_dict(torch.load('your_model.pth'))

# 初始化SHAP解释器
explainer = shap.DeepExplainer(model, X_train)

# 计算特征重要性
# shap_values = explainer(X_test_tensor)

shap_values = explainer.shap_values(X_test)

print(len(shap_values), X_test.shape, shap_values[1].shape)

# torch.mean(shap_values[0], dim=0)
# 定义柱状图的数据
categories = [f'feature{i}' for i in range(shap_values[0].shape[1])]
values = np.mean(shap_values[0], axis=0)
values2 = np.mean(shap_values[1], axis=0)
# 创建一个图形窗口
plt.figure()

# 绘制横向柱状图
plt.barh(categories[:values.shape[0]], values)
# plt.barh(categories[:values.shape[0]], values)
# 添加标题和坐标轴标签
# plt.title("横向柱状图")
# plt.xlabel("数值")
# plt.ylabel("类别")

# 显示图形
plt.show()

process_models/q1/q1b.py
TheWangYang/2023_huawei_mathematical_modeling
import pandas as pd
import os
import numpy as np
from sklearn.linear_model import LogisticRegression, SGDOneClassSVM
from sklearn.neural_network import MLPClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
import matplotlib.pyplot as plt

from sklearn.preprocessing import MinMaxScaler

from scipy.stats import pearsonr

from sklearn.neighbors import NearestNeighbors

plt.rcParams['font.sans-serif'] = ['SimHei']

# 得到训练和验证集
def get_train_valid_q1():
    root_path = '/data1/wyy/projects/_competition/math_model_v2/Data'
    excel_name1 = '1.xlsx'
    excel_name2 = '2.xlsx'
    excel_name3 = '3.xlsx'
    excel_name_time = 'additional.xlsx'

    excel_answer = '4.xlsx'
    
    # 读取Excel文件
    table1 = pd.read_excel(os.path.join(root_path, excel_name1))
    table2 = pd.read_excel(os.path.join(root_path, excel_name2))
    table3 = pd.read_excel(os.path.join(root_path, excel_name3))
    table_time = pd.read_excel(os.path.join(root_path, excel_name_time))
    table_answer = pd.read_excel(os.path.join(root_path, excel_answer))
    
    all_Base_info = np.array([[]], dtype=np.float64)
    
    # print(table1['年龄'])
    # print(table1['年龄'])
    # series = table1.set_index(['年龄'], drop=True)
    # plt.figure(figsize=(10, 6))
    # series['脑出血前mRS评分'].plot()
    # plt.show()
    
    # 遍历所有160个病人
    for i in range(160):
        info_dict = {}
        x = np.array([], dtype=np.float64)
        row1 = table1.loc[i]
        row2 = table2.loc[i]
        
        # 遍历每个病人的第index=4列到index=19列
        for j in range(4, 4+19):
            if j == 2:
                continue
            elif j == 5: # 对应性别列
                # 男为0，女为1
                sex = 0 if row1[j] == '男' else 1
                x = np.append(x, sex)
            elif j==15:  # 对应血压列
                # 获得血压数值，存储为列表格式[]
                xueya = np.array([eval(k) for k in row1[j].split('/')])
                x = np.append(x, xueya)
            else:
                # 将其他列对应的数值添加到x中
                x = np.append(x, row1[j])
                
        # 将x
        x = x[np.newaxis]

        if all_Base_info.size == 0:
            all_Base_info = x  # 将其中对应的x传递给all_Base_info
        else:
            all_Base_info=np.concatenate((all_Base_info, x), axis=0)
    
    
    # 将table1中的性别列从字符串转换为0或1
    table1['性别'] = table1['性别'].map({'男':1, '女':0})
    
    # 选择的特征列
    info_select = ['年龄', '性别', '高血压病史', '卒中病史', '糖尿病史']
    
    # 得到表格中对应上述选择特征列的信息
    all_Base_info = table1[info_select].values.astype(np.float64)
    
    # 得到表格2中对应的所有的HM体积信息
    all_HM_info = table2['HM_volume'].values.astype(np.float64).reshape(-1, 1)
    
    # 首次检查对应的影像检查结果
    # index from 3 to 23
    for i in range(3, 24):
        tem_data = table2[table2.columns[i]].values.astype(np.float64).reshape(-1, 1)
        # 拼接数据
        all_HM_info = np.concatenate((all_HM_info, tem_data), axis=1)
        # break
    
    print(all_HM_info.shape)

    first_check_seq = table1['入院首次影像检查流水号'].values
    first_check_seq_2 = table2['首次检查流水号'].values
    first_check_seq_3 = table3['流水号'].values
    first_check_seq_time = table_time['入院首次检查流水号'].values
    
    sub_name = table2['ID']
    print("异常数据 in table1, 2\n 流水号不对应:\n", sub_name[first_check_seq!=first_check_seq_2])
    print("异常数据 in table1, time\n 流水号不对应:\n", sub_name[first_check_seq!=first_check_seq_time]) 
    print("异常数据 in table2, time\n 流水号不对应:\n", sub_name[first_check_seq_2!=first_check_seq_time]) 
    print("异常数据: table1与table2=time的流水号存在不对应的情况")

    shape_query_buf = []
    table12_query_buf = []

    # table12_test_query_buf
    for seq in first_check_seq:
        # print(np.where(first_check_seq_3==seq)[0])
        if (np.where(first_check_seq_3==seq)[0]).size == 0:
            shape_query_buf.append(None)
            print(np.where(first_check_seq_3==seq), seq, sub_name[first_check_seq_2==seq])
        else:
            # 得到index索引值在100以内或大于等于132的
            if np.where(first_check_seq==seq)[0][0]<100 or np.where(first_check_seq==seq)[0][0]>=132:
                table12_query_buf.append(np.where(first_check_seq==seq)[0][0])
                shape_query_buf.append(np.where(first_check_seq_3==seq)[0][0])
    
    # X = np.concatenate((X, HM_info), axis=1)

    all_shape_info = table3[table3.columns[2]].values.astype(np.float64).reshape(-1, 1)
    
    for i in range(3, 33):
        # print(table3.columns[i])
        tem_data = table3[table3.columns[i]].values.astype(np.float64).reshape(-1, 1)
        all_shape_info = np.concatenate((all_shape_info, tem_data), axis=1)

    
    shape_info = all_shape_info[shape_query_buf]
    Base_info = all_Base_info[table12_query_buf]
    
    # HM_info = all_HM_info[table12_query_buf]
    print(shape_info.shape)
    
    # X = np.concatenate((X, HM_info), axis=1)
    # X = np.concatenate((X, shape_info), axis=1)
    
    # 构建X和Y作为train和valid数据
    X = Base_info
    # X = np.concatenate((X, HM_info), axis=1)
    # X = np.concatenate((X, shape_info), axis=1)
    # Y = table_answer['是否发生血肿扩张']
    Y = table_answer['问题1：血肿扩张'].values[2:].astype(np.float64)[table12_query_buf]

    seq = []
    p_value_buf = []
    
    for i in range(X.shape[1]):
        pearson_corr, p_value = pearsonr(X[:, i], Y)
        
        if np.isnan(p_value):
            p_value = 0.
        p_value_buf.append(p_value)
        # print(p_value)
        if p_value>0.1:
            seq.append(i)

    p_value_buf = np.array(p_value_buf)
    normalize = plt.Normalize(vmin=p_value_buf.min(), vmax=p_value_buf.max())

    # cmap = plt.cm.get_cmap('viridis')
    # # 将浮点数映射到RGB颜色空间
    # colors = cmap(normalize(p_value_buf))
    # print(colors)
    # colors[seq, 1] = 0
    # colors[seq, 0] = 1.
    # colors[seq, 2] = 0
    # plt.figure()

    # # 添加标题和坐标轴标签
    # plt.title("Pearson 相关系数")
    # plt.xlabel("Features")
    # plt.ylabel("Pearson correlation coefficient")

    # plt.show()
    # plt.bar(np.arange(len(p_value_buf)), p_value_buf, color=colors)
    # plt.savefig("Pearson.jpg")
    # plt.show()


    # X = X[:, seq]
    means = np.mean(X, axis=0)
    # # print(means)
    # stds = np.std(X, axis=0)

    X = (X-means)
    scaler = MinMaxScaler(feature_range=(-1, 1))

    # X = X[:, stds!=0]
    # X = np.random.random((100, 49))
    # Y = np.random.randint(0, 2, (100,))
    scaler = scaler.fit(X)

    # print(stds!=0)
    X_train = X[:100]
    y_train = Y[:100]

    X_valid = X[100:]
    y_valid = Y[100:]
   
    # X_train, X_valid, y_train, y_valid = train_test_split(X, Y, train_size=0.8, test_size=0.2, random_state=0)
    X_train = scaler.transform(X_train)
    X_valid = scaler.transform(X_valid)
    # print(X_train[:, 3])

    # for j in range(72):
    #     print(j)
    #     plt.plot(X_train[:, j])
    #     plt.show()

    # print(X_train[0])
    print("X_train: {}, X_valid: {}, y_train: {}, y_valid: {}".format(X_train, X_valid, y_train, y_valid))
    return X_train, X_valid, y_train, y_valid

if __name__=='__main__':
    get_train_valid_q1()
#     return

# print(X_train.shape)

# model = LogisticRegression(max_iter=1000)
# model.fit(X_train, y_train)

# # Predict on validation set
# y_pred = model.predict(X_valid)

# # Calculate accuracy
# accuracy = accuracy_score(y_valid, y_pred)

# print(accuracy)

# svm = SGDOneClassSVM(max_iter=1000)
# svm.fit(X_train, y_train)
# y_pred_svm = svm.predict(X_valid)

# accuracy_svm = accuracy_score(y_valid, y_pred_svm)
# print(accuracy_svm)


    
process_models/q1/q1b_dataloader.py
TheWangYang/2023_huawei_mathematical_modeling
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from process_models.q1.q1b_dataloader import get_train_valid_q1

from imblearn.over_sampling import SMOTE

import matplotlib.pyplot as plt
import numpy as np
import shap

# 生成示例数据
X_train, X_test, y_train, y_test = get_train_valid_q1()

sm = SMOTE(random_state=0)

X_train, y_train = sm.fit_resample(X_train, y_train)
# X = np.random.random((100, 15))
# Y = np.random.randint(0, 2, (100,))

print(X_train.shape)


# X_train, X_test, y_train, y_test = train_test_split(X, Y, train_size=0.8, test_size=0.2, random_state=0)
# 转换数据为PyTorch的Tensor类型
X_train = torch.Tensor(X_train)
print(X_train.shape)
y_train = torch.LongTensor(y_train)
X_test = torch.Tensor(X_test)
y_test = torch.LongTensor(y_test)

# 创建数据集和数据加载器

# X_train = torch.cat([X_train, X_test], dim=0)
# y_train = torch.cat([y_train, y_test], dim=0)
train_dataset = TensorDataset(X_train, y_train)
train_loader = DataLoader(train_dataset, batch_size=10, shuffle=False)

# 构建神经网络模型
input_n = X_train.shape[-1]
start_n = 64
middle_n = [64, 64, 32, 16, 2]

print(input_n)
torch.manual_seed(42)


class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.fc = nn.Sequential(
            
            nn.Linear(input_n, start_n),
            # nn.BatchNorm1d(100),
            nn.ReLU(),
        )

        for i, num in enumerate(middle_n):
            if i==0:
                self.fc.append(nn.Sequential(nn.Linear(start_n, num),
                                             nn.ReLU()))
            else:
                self.fc.append(nn.Sequential(nn.Linear(middle_n[i-1], num),
                                             nn.ReLU()))

        self.fc.append(nn.Softmax(dim=1))

    def forward(self, x):
        x = self.fc(x)
        return x

model = Net()
model.train()
if True:
    # model.load_state_dict(torch.load('model_weights_base_info.pth'))
    model.load_state_dict(torch.load('model_weights.pth'))
    print('loaded')

# 定义损失函数和优化器
criterion = nn.CrossEntropyLoss()
optimizer = optim.Adam(model.parameters(), lr=1e-4)
# scheduler = optim.lr_scheduler.StepLR(optimizer, step_size=15, gamma=0.1)
# 训练模型
num_epochs = 1000
with torch.no_grad():
    outputs = model(X_test)
    # print(outputs.data)
    _, predicted = torch.max(outputs.data, 1)

    accuracy = (predicted == y_test).sum().item() / len(y_test)

    print("Accuracy:", accuracy)

for epoch in range(num_epochs):
    # scheduler.step()
    for i, (inputs, labels) in enumerate(train_loader):
        optimizer.zero_grad()
        outputs = model(inputs)
        # print(outputs.data)
        # print(outputs.data)
        loss = criterion(outputs, labels)
        loss.backward()
        optimizer.step()
        # print(loss)

    if epoch%10 ==1 :
        with torch.no_grad():
            outputs = model(X_train)
            print(loss)
            _, predicted = torch.max(outputs.data, 1)

            print(outputs.data[-1])
            print((y_train==1).sum()/len(y_train))
            # 输出准确率
            accuracy = (predicted == y_train[:]).sum().item() / len(y_train)
            print("TP:", y_train[predicted == y_train].sum().item()/(y_train==1).sum())
            print("TF:", ((predicted == y_train).sum().item()-y_train[predicted == y_train].sum().item())/(y_train==0).sum())
            
            print("Accuracy:", accuracy)
            # print("Accuracy:", accuracy)

    
    # print(f"Epoch [{epoch+1}/{num_epochs}], Loss: {loss.item()}")

# 在测试集上进行预测
with torch.no_grad():
    outputs = model(X_test)
    
    _, predicted = torch.max(outputs.data, 1)
    print(outputs.data)

    # 输出准确率
    accuracy = (predicted == y_test).sum().item() / len(y_test)

    print("TP:", y_test[predicted == y_test].sum().item()/(y_test==1).sum())
    print("TF:", ((predicted == y_test).sum().item()-y_test[predicted == y_test].sum().item())/(y_test==0).sum())
    
    print("Accuracy:", accuracy)
    print("y_test:", y_test)

# Importance analythis

torch.save(model.state_dict(), 'model_weights.pth')
    # print('X_text', X_test)

# # Importance analythis
# explainer = shap.Explainer(model, X_train)

# shap_values = explainer(X_test)

# # 绘制特征重要性图
# shap.summary_plot(shap_values, X_test_tensor, plot_type='bar')

# # 显示图形
# plt.show()
process_models/q1/q1b_fc_net.py
TheWangYang/2023_huawei_mathematical_modeling
import torch
import torch.nn as nn
import torch.optim as optim


# 构建SVM分类器
class SVM(nn.Module):
    def __init__(self, input_dim):
        super(SVM, self).__init__()
        self.linear = nn.Linear(input_dim, 1)

    def forward(self, x):
        return self.linear(x)
    
# 假设有训练数据和标签
train_data = torch.randn(100, 100)
labels = torch.randint(0, 2, (100,))


svm_model = SVM(100)
optimizer = optim.SGD(svm_model.parameters(), lr=0.01)

num_epochs = 100

for epoch in range(num_epochs):
    # 前向传播
    outputs = svm_model(train_data)

    # 计算损失
    loss = torch.mean(torch.clamp(1 - labels * outputs, min=0))

    # 反向传播和优化
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

    # 打印训练过程中的损失
    if (epoch+1) % 10 == 0:
        print('Epoch [{}/{}], Loss: {:.4f}'.format(epoch+1, num_epochs, loss.item()))
        


process_models/q1/q1b_svm.py
TheWangYang/2023_huawei_mathematical_modeling
import pandas as pd
import os
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.mixture import GaussianMixture


# 定义存放数据的路径
root_path = "/data1/wyy/projects/_competition/mathmodel_code/data"

# 定义表格名称
table1_name = "1.xlsx"
table2_name = "2.xlsx"
table3_name = "3.xlsx"
# 4-答案表
table4_name = "4.xlsx"
table_additional = "additional.xlsx"


# 读取Excel文件
table1 = pd.read_excel(os.path.join(root_path, table1_name))
table2 = pd.read_excel(os.path.join(root_path, table2_name))
table3 = pd.read_excel(os.path.join(root_path, table3_name))
table_addi = pd.read_excel(os.path.join(root_path, table_additional))


'''
得到X轴和y轴需要的数据
root_path为存放数据的路径地址
'''
def get_X_y():
    # 创建需要返回的X和y
    X = []
    y = []
    
    # 根据表2中的每个检查点对应的ED_volume以及（该检查点对应的时间-首次检查点对应的时间+发病到首次检查点对应的时间）
    for i in range(0, 100):
        # 获得表格中每行数据
        table1_row = table1.loc[i]
        table2_row = table2.loc[i]
        # print("ID: {}".format(table2_row[0]))
        
        table3_row = table3.loc[i]
        table_addi_row = table_addi.loc[i]
                
        # 先得到表2中对应的ED_volume中，判断是否为空，不为空得到对应的检查点到发病的时间间隔
        # table2_ED_volume_index = [13, 36, 59, 82, 105, 128, 151, 174, 197]
        # 只记录患者从发病到首次检查时间
        table2_ED_volume_index = [13]
        
        # 遍历其中的每个ED_volume_index得到y并计算对应的检查点时间间隔为x
        for index in table2_ED_volume_index:
            # 判断是否为空，不为空进入逻辑
            if not np.isnan(table2_row[index]):
                # 打印table2_row[index]数据类型
                # print(type(table2_row[index]))
                # 根据此时的检查点计算从发病到此时间点的间隔
                curr_check_number = table2_row[index - 12]
                # 根据流水号在附表中找到对应流水号在附表中的列index
                # 如果没有找到表示出现错误，直接跳过即可
                if (np.where(curr_check_number == table_addi_row)[0]).size != 0:
                    check_number_addi_index = np.where(curr_check_number == table_addi_row)[0][0]
                    # print("check_number_addi_index: {}".format(check_number_addi_index))
                    # 得到检查流水号对应的时间点列索引index
                    time_for_check_number_addi_index = check_number_addi_index - 1
                    # 然后计算当前时间点距离首次检查时间点的小时数
                    interval_to_first_time = (table_addi_row[time_for_check_number_addi_index] - table_addi_row[2]).total_seconds() / 3600
                    # 然后，从表1中获得从发病到首次检查的时间间隔
                    first_check_interval_to_onset = table1_row[14]
                    # 然后得到从当前检查点时间距离发病的时间间隔
                    curr_interval = interval_to_first_time + first_check_interval_to_onset
                    # 将当前检查点对应的时间间隔作为x轴，添加到X列表中
                    # print("curr_interval: {}".format(curr_interval))
                    X.append(curr_interval)
                    # 将该值添加到y列表中
                    y.append(table2_row[index])
        
    # 打印X和y
    print("X size: {}".format(len(X)))
    print("y size: {}".format(len(y)))
    
    return X, y

# 使用曲线拟合
def regression_linear(X, y):
    
    
    print(np.max(X))
    print(np.max(y))
    
    # ------------------------------------使用高斯拟合曲线-------------------------------------
    # 转换为NumPy数组
    X = np.array(X).reshape(-1, 1)
    y = np.array(y)

    # 高斯混合模型拟合
    gmm = GaussianMixture(n_components=2)  # 可根据需要调整混合成分数
    gmm.fit(X, y)

    # # 生成拟合曲线的X值
    X_fit = np.linspace(min(X), max(X), 1000).reshape(-1, 1)
    y_fit = gmm.predict(X_fit)

    # 绘图
    plt.scatter(X, y, label='Original Data')
    plt.plot(X_fit, y_fit, color='red', label='Fitted Curve (Gaussian Mixture)')
    plt.xlabel('X')
    plt.ylabel('y')
    plt.legend()
    plt.savefig("process_models/q2/output/q2a_regression_gaussian.png")
    plt.clf()
    
    # ----------------------------------使用多项式拟合曲线--------------------------------------
    # 多项式回归拟合
    degree = 3  # 多项式阶数，可以根据需要调整
    model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
    model.fit(X, y)

    # 生成拟合曲线的X值
    X_fit = np.linspace(min(X), max(X), 1000).reshape(-1, 1)
    y_fit = model.predict(X_fit)

    # 绘图
    plt.scatter(X, y, label='Original Data')
    plt.plot(X_fit, y_fit, color='red', label=f'Fitted Curve (degree={degree})')
    plt.xlabel('X')
    plt.ylabel('y')
    plt.legend()
    plt.savefig("process_models/q2/output/q2a_regression_polynomial.png")
    
    # ------------------------------------得到预测结果-----------------------------------------
    # 得到从发病到首次检查时间点的时间间隔，得到X轴表示
    # 根据表2中的每个检查点对应的ED_volume以及（该检查点对应的时间-首次检查点对应的时间+发病到首次检查点对应的时间）
    X_new = []
    y_new = []
    
    for i in range(0, 100):
        # 获得表格中每行数据
        table1_row = table1.loc[i]
        table2_row = table2.loc[i]
        
        X_new.append(table1_row[14])
        y_new.append(table2_row[13])
    
    X_new = np.array(X_new).reshape(-1, 1)
    # 使用高斯曲线预测
    # y_pred = gmm.predict(X_new)
    # 使用多项式预测
    y_pred = model.predict(X_new)
    
    # # 计算拟合曲线和真实值之间的残差
    residuals = np.abs(y_new - y_pred)
    
    print(y_new, y_pred, residuals)
    
    # 填充表4答案表格中对应的残差列
    # pd.set_option('display.float_format', lambda x: '%.3f' % x)
    table4_answer = pd.read_excel(os.path.join(root_path, table4_name))
    
    # print(table4_answer)
    
    # 使用 iloc 给特定列赋值
    table4_answer.iloc[2: 102, 5] = residuals
    table4_answer.to_excel(os.path.join(root_path, table4_name), index=False, float_format='%.3f')
    
    print("success")
    return
    
    
# 单元测试代码
if __name__ == "__main__":
    X, y = get_X_y()
    regression_linear(X, y)
    
    
    
process_models/q2/q2a_regression.py
TheWangYang/2023_huawei_mathematical_modeling
import pandas as pd
import os
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.mixture import GaussianMixture


# 定义存放数据的路径
root_path = "/data1/wyy/projects/_competition/mathmodel_code/data"

# 定义表格名称
table1_name = "1.xlsx"
table2_name = "2.xlsx"
table3_name = "3.xlsx"
# 4-答案表
table4_name = "4.xlsx"
table_additional = "additional.xlsx"


# 读取Excel文件
table1 = pd.read_excel(os.path.join(root_path, table1_name))
table2 = pd.read_excel(os.path.join(root_path, table2_name))
table3 = pd.read_excel(os.path.join(root_path, table3_name))
table_addi = pd.read_excel(os.path.join(root_path, table_additional))


# 根据现有的index遍历获得同一个亚组的病人的水肿体积随时间变化情况
def get_X_y(index_keep):
    # 创建需要返回的X和y
    X = []
    y = []
    
    # 根据表2中的每个检查点对应的ED_volume以及（该检查点对应的时间-首次检查点对应的时间+发病到首次检查点对应的时间）
    for i in index_keep:
        # 获得表格中每行数据
        table1_row = table1.loc[i]
        table2_row = table2.loc[i]
        # print("ID: {}".format(table2_row[0]))
        
        table3_row = table3.loc[i]
        table_addi_row = table_addi.loc[i]
                
        # 先得到表2中对应的ED_volume中，判断是否为空，不为空得到对应的检查点到发病的时间间隔
        # table2_ED_volume_index = [13, 36, 59, 82, 105, 128, 151, 174, 197]
        # 只记录患者从发病到首次检查时间
        table2_ED_volume_index = [13]
        
        # 遍历其中的每个ED_volume_index得到y并计算对应的检查点时间间隔为x
        for index in table2_ED_volume_index:
            # 判断是否为空，不为空进入逻辑
            if not np.isnan(table2_row[index]):
                # 打印table2_row[index]数据类型
                # print(type(table2_row[index]))
                # 根据此时的检查点计算从发病到此时间点的间隔
                curr_check_number = table2_row[index - 12]
                # 根据流水号在附表中找到对应流水号在附表中的列index
                # 如果没有找到表示出现错误，直接跳过即可
                if (np.where(curr_check_number == table_addi_row)[0]).size != 0:
                    check_number_addi_index = np.where(curr_check_number == table_addi_row)[0][0]
                    # print("check_number_addi_index: {}".format(check_number_addi_index))
                    # 得到检查流水号对应的时间点列索引index
                    time_for_check_number_addi_index = check_number_addi_index - 1
                    # 然后计算当前时间点距离首次检查时间点的小时数
                    interval_to_first_time = (table_addi_row[time_for_check_number_addi_index] - table_addi_row[2]).total_seconds() / 3600
                    # 然后，从表1中获得从发病到首次检查的时间间隔
                    first_check_interval_to_onset = table1_row[14]
                    # 然后得到从当前检查点时间距离发病的时间间隔
                    curr_interval = interval_to_first_time + first_check_interval_to_onset
                    # 将当前检查点对应的时间间隔作为x轴，添加到X列表中
                    # print("curr_interval: {}".format(curr_interval))
                    X.append(curr_interval)
                    # 将该值添加到y列表中
                    y.append(table2_row[index])
        
    # 打印X和y
    print("X size: {}".format(len(X)))
    print("y size: {}".format(len(y)))
    
    return X, y

# 使用曲线拟合
def regression_linear(label, index_keep, X, y):
    
    print(np.max(X))
    print(np.max(y))
    
    # ------------------------------------使用高斯拟合曲线-------------------------------------
    # 转换为NumPy数组
    X = np.array(X).reshape(-1, 1)
    y = np.array(y)

    # 高斯混合模型拟合
    gmm = GaussianMixture(n_components=2)  # 可根据需要调整混合成分数
    gmm.fit(X, y)

    # # 生成拟合曲线的X值
    X_fit = np.linspace(min(X), max(X), 1000).reshape(-1, 1)
    y_fit = gmm.predict(X_fit)

    # 绘图
    plt.scatter(X, y, label='Original Data')
    plt.plot(X_fit, y_fit, color='red', label='Fitted Curve (Gaussian Mixture)')
    plt.xlabel('X')
    plt.ylabel('y')
    plt.legend()
    # 保存的时候使用label作为区分
    plt.savefig("process_models/q2/output/q2b_regression_gaussian_{}.png".format(label))
    # 清除画布
    plt.clf()
    
    # ----------------------------------使用多项式拟合曲线--------------------------------------
    # 多项式回归拟合
    degree = 3  # 多项式阶数，可以根据需要调整
    model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
    model.fit(X, y)

    # 生成拟合曲线的X值
    X_fit = np.linspace(min(X), max(X), 1000).reshape(-1, 1)
    y_fit = model.predict(X_fit)

    # 绘图
    plt.scatter(X, y, label='Original Data')
    plt.plot(X_fit, y_fit, color='red', label=f'Fitted Curve (degree={degree})')
    plt.xlabel('X')
    plt.ylabel('y')
    plt.legend()
    # 使用label区分是哪个类别
    plt.savefig("process_models/q2/output/q2b_regression_polynomial_{}.png".format(label))
    # 清除画布
    plt.clf()
    
    # ------------------------------------得到预测结果-----------------------------------------
    # 得到从发病到首次检查时间点的时间间隔，得到X轴表示
    # 根据表2中的每个检查点对应的ED_volume以及（该检查点对应的时间-首次检查点对应的时间+发病到首次检查点对应的时间）
    X_new = []
    y_new = []
    
    # 遍历index_keep
    for i in index_keep:
        # 获得表格中每行数据
        table1_row = table1.loc[i]
        table2_row = table2.loc[i]
        
        X_new.append(table1_row[14])
        y_new.append(table2_row[13])
    
    X_new = np.array(X_new).reshape(-1, 1)
    # 使用高斯曲线预测
    # y_pred = gmm.predict(X_new)
    # 使用多项式预测
    y_pred = model.predict(X_new)
    
    # # 计算拟合曲线和真实值之间的残差
    residuals = np.abs(y_new - y_pred)
    
    print(y_new, y_pred, residuals)
    
    # 填充表4答案表格中对应的残差列
    # pd.set_option('display.float_format', lambda x: '%.3f' % x)
    table4_answer = pd.read_excel(os.path.join(root_path, table4_name))
    
    # print(table4_answer)
    
    # 使用 iloc 给特定列赋值
    index_keep_plus_2 = [x + 2 for x in index_keep]
    table4_answer.iloc[index_keep_plus_2, 6] = residuals
    label_list = [label] * len(index_keep)
    table4_answer.iloc[index_keep_plus_2, 7] = label_list
    table4_answer.to_excel(os.path.join(root_path, table4_name), index=False, float_format='%.3f')
    
    print("success")
    return


'''
聚类函数
其中K为聚类数量为存放数据的路径地址
'''

def cluster_K_group():
    # 创建需要返回的X
    X = []

    # 得到表1中病人对应的特征
    for i in range(100):  # You can use range(100) directly
        # 获得表格中每行数据
        table1_row = table1.loc[i]

        # 提取对应的特征列
        feature1 = np.array([table1_row[1]])  # 对应90天mRS
        feature2 = np.array([table1_row[4]])

        # 处理性别特征
        feature_age = np.array([1]) if table1_row[5] == '男' else np.array([0])

        feature3 = table1_row[6:15].values.flatten()

        # 处理血压
        xueya = np.array([eval(k) for k in table1_row[15].split('/')])

        feature5 = table1_row[16:23].values.flatten()

        # 合并特征
        combined_feature = np.concatenate((feature1, feature2, feature_age, feature3, xueya, feature5))
        X.append(combined_feature)

    # 转换为NumPy数组
    X = np.array(X)

    # 使用KMeans进行聚类
    n_clusters = min(5, len(X))  # 不超过5个簇
    kmeans = KMeans(n_clusters=n_clusters)
    kmeans.fit(X)

    # 获取簇标签
    # 对应100前100个病人
    cluster_labels = kmeans.labels_
    print("cluster_labels: {}".format(cluster_labels))
    print("cluster_labels.size: {}".format(cluster_labels.size))
    label_counts = pd.Series(cluster_labels).value_counts()
    print("label_counts: {}".format(label_counts))
    
    
    # 建立保存对应的5个类别[0, 1, 2, 3, 4]对应的病人在表1中的index
    label0_keep = []
    label1_keep = []
    label2_keep = []
    label3_keep = []
    label4_keep = []
    
    # 根据标签，将前100个病人进行分组
    # 共分为5个类别
    for i, label in enumerate(cluster_labels):
        # 得到表格的第i行
        table1_row = table1.loc[i]
        
        if label == 0:
            label0_keep.append(i)
        elif label == 1:
            label1_keep.append(i)
        elif label == 2:
            label2_keep.append(i)
        elif label == 3:
            label3_keep.append(i)
        elif label == 4:
            label4_keep.append(i)
        
    # 对于每个保存的keep列表使用与2a中相同的方法进行回归预测
    X,y = get_X_y(label0_keep)
    regression_linear(0, label0_keep, X, y)
    
    X,y = get_X_y(label1_keep)
    regression_linear(1, label1_keep, X, y)
    
    X,y = get_X_y(label2_keep)
    regression_linear(2, label2_keep, X, y)
    
    X,y = get_X_y(label3_keep)
    regression_linear(3, label3_keep, X, y)
    
    X,y = get_X_y(label4_keep)
    regression_linear(4, label4_keep, X, y)
    
    print("-----------------success---------------")
    

# 单元测试代码
if __name__ == "__main__":
    cluster_K_group()
    
    
    



process_models/q2/q2b_group.py
TheWangYang/2023_huawei_mathematical_modeling
import os
import numpy as np
import pandas as pd
from sklearn.tree import DecisionTreeRegressor
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm



# 定义存放数据的路径
root_path = "/data1/wyy/projects/_competition/mathmodel_code/data"

# 定义表格名称
table1_name = "1.xlsx"
table2_name = "2.xlsx"
table3_name = "3.xlsx"
# 4-答案表
table4_name = "4.xlsx"
table_additional = "additional.xlsx"


# 读取Excel文件
table1 = pd.read_excel(os.path.join(root_path, table1_name))
table2 = pd.read_excel(os.path.join(root_path, table2_name))
table3 = pd.read_excel(os.path.join(root_path, table3_name))
table_addi = pd.read_excel(os.path.join(root_path, table_additional))


def features_impact_ED_volume():
    # 准备存储特征和目标数据的列表
    all_features = []
    all_targets = []

    # 分析前100个病人
    for i in range(100):
        # 获得表格中每行数据
        table1_row = table1.loc[i]
        table2_row = table2.loc[i]
        
        table3_row = table3.loc[i]
        table_addi_row = table_addi.loc[i]
                
        # 得到当前第一次复检点与首次检查得到的ED_volume的差值，获得增量
        curr_to_first_check_interval_ED_volume = table2_row[36] - table2_row[13]
        
        # 提取16-22列特征
        features = table1_row[16:23].values

        # 存储特征和目标数据
        all_features.append(features)
        all_targets.append(curr_to_first_check_interval_ED_volume)

    # 合并特征和目标变量为一个数据集
    features_array = np.array(all_features)
    targets_array = np.array(all_targets)

    # 构建决策树模型
    model = DecisionTreeRegressor()
    model.fit(features_array, targets_array)

    # 获取特征重要性
    feature_importances = model.feature_importances_
    
    print("feature_importances: {}".format(feature_importances))
    
    # 中文显示设置
    fontprop = fm.FontProperties(fname='/data1/wyy/projects/_competition/mathmodel_code/utils/fonts/TIMESBD.TTF')  # 替换成你的字体文件路径
    plt.rcParams['font.family'] = fontprop.get_name()
    
    # 绘制条形图
    plt.figure(figsize=(10, 6))
    # plt.bar(['脑室引流', '止血治疗', '降颅压治疗', '降压治疗', '镇静、镇痛治疗', '止吐护胃', '营养神经'], feature_importances)
    plt.bar(['0', '1', '2', '3', '4', '5', '6'], feature_importances)
    plt.xlabel('Features')
    plt.ylabel('Importance')
    plt.title('Feature Importance from Decision Tree')
    plt.xticks(rotation=45)  # 旋转x轴标签，以避免重叠
    plt.savefig("process_models/q2/output/q2c_features_decision_tree.png")
    plt.clf()
    
    print("--------------------success-------------------")

# 调用函数进行特征影响分析
features_impact_ED_volume()

process_models/q2/q2c_features.py
TheWangYang/2023_huawei_mathematical_modeling
import os
import numpy as np
import pandas as pd
from sklearn.tree import DecisionTreeRegressor
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm


# 题目要求：分析血肿体积，水肿体积，和治疗方法三者之间的关系
# 思路：首先得到（第一复检点的血肿体积-首次检查对应的血肿体积），（第一复检点-首次检查对应的水肿体积），（获得治疗方法特征）
# 然后使用关系分析方法来分析三者之间的关系

# 定义存放数据的路径
root_path = "/data1/wyy/projects/_competition/mathmodel_code/data"

# 定义表格名称
table1_name = "1.xlsx"
table2_name = "2.xlsx"
table3_name = "3.xlsx"
# 4-答案表
table4_name = "4.xlsx"
table_additional = "additional.xlsx"


# 读取Excel文件
table1 = pd.read_excel(os.path.join(root_path, table1_name))
table2 = pd.read_excel(os.path.join(root_path, table2_name))
table3 = pd.read_excel(os.path.join(root_path, table3_name))
table_addi = pd.read_excel(os.path.join(root_path, table_additional))

"""
定义函数：用于分析三者之间的关系
"""
def relationship_between_three():
    # 准备存储特征和目标数据的列表
    all_features = []
    all_HM_targets = []
    all_ED_targets = []
    
    # 分析前100个病人
    for i in range(100):
        # 获得表格中每行数据
        table1_row = table1.loc[i]
        table2_row = table2.loc[i]
        
        table3_row = table3.loc[i]
        table_addi_row = table_addi.loc[i]
                
        # 得到当前第一次复检点与首次检查得到的ED_volume的差值，获得增量
        curr_to_first_check_interval_ED_volume = table2_row[36] - table2_row[13]
        
        # 得到当前第一次复检点与首次检查得到的HM_volume的差值，获得增量
        curr_to_first_check_interval_HM_volume = table2_row[25] - table2_row[2]
        
        # 提取16-22列特征
        features = table1_row[16:23].values

        # 存储特征和目标数据
        all_features.append(features)
        all_HM_targets.append(curr_to_first_check_interval_HM_volume)
        all_ED_targets.append(curr_to_first_check_interval_ED_volume)
        

    # 合并特征和目标变量为一个数据集
    # 修正可能导致错误的代码
    features_array = np.array(all_features).astype('float16')
    HM_targets_array = np.array(all_HM_targets).astype('float16')
    ED_targets_array = np.array(all_ED_targets).astype('float16')
    
    
    # 计算HM目标与特征的相关系数
    hm_correlations = np.corrcoef(features_array.T, HM_targets_array)[0, 1:]

    # 计算ED目标与特征的相关系数
    ed_correlations = np.corrcoef(features_array.T, ED_targets_array)[0, 1:]

    # 打印相关系数
    print("Correlation between features and HM_targets:", hm_correlations)
    print("Correlation between features and ED_targets:", ed_correlations)

    # 绘制散点图
    plt.figure(figsize=(12, 6))
    plt.subplot(1, 2, 1)
    plt.bar([0, 1, 2, 3, 4, 5, 6], hm_correlations, alpha=0.5)
    plt.xlabel('Features')
    plt.ylabel('HM Targets')
    plt.title('Scatter Plot: Features vs HM Targets')

    plt.subplot(1, 2, 2)
    plt.bar([0, 1, 2, 3, 4, 5, 6], ed_correlations, alpha=0.5)
    plt.xlabel('Features')
    plt.ylabel('ED Targets')
    plt.title('Scatter Plot: Features vs ED Targets')
    plt.tight_layout()
    plt.savefig("process_models/q2/output/q2d_relation.png")

    print("--------------------success-------------------")


if __name__ == "__main__":
    relationship_between_three()







process_models/q2/q2d_relation.py
TheWangYang/2023_huawei_mathematical_modeling
import os
import numpy as np
import pandas as pd
from sklearn.tree import DecisionTreeRegressor
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm


# 题目要求：结合多个特征对mRS进行预测
# 思路：保存前100个病人的特征，然后作为X，同时以mRS为目标值


# 定义存放数据的路径
root_path = "/data1/wyy/projects/_competition/mathmodel_code/data"

# 定义表格名称
table1_name = "1.xlsx"
table2_name = "2.xlsx"
table3_name = "3.xlsx"
# 4-答案表
table4_name = "4.xlsx"
table_additional = "additional.xlsx"


# 读取Excel文件
table1 = pd.read_excel(os.path.join(root_path, table1_name))
table2 = pd.read_excel(os.path.join(root_path, table2_name))
table3 = pd.ExcelFile(os.path.join(root_path, table3_name))

table3_ED = None
table3_HM = None

# 遍历表格3中每个工作表
for sheet_name in table3.sheet_names:
    # 在这里可以对每个工作表的 DataFrame 进行相应的操作
    # 例如，打印工作表名称和 DataFrame 的头部
    print(f'工作表名称: {sheet_name}')
    
    if sheet_name == "ED":
        table3_ED = pd.read_excel(table3, "ED")
    elif sheet_name == "Hemo":
        table3_HM = pd.read_excel(table3, "Hemo")

# 读取附加表
table_addi = pd.read_excel(os.path.join(root_path, table_additional))


# 预测mRS
def get_X_y_for_mRS(k):
    # 定义保存前100个病人对应的mRS数值
    mRS_targets = []
    
    # 汇总特征列
    all_features = []
    
    # 分析前100个病人
    for i in range(k):
        # 获得表格中每行数据
        table1_row = table1.loc[i]
        table2_row = table2.loc[i]
        table_addi_row = table_addi.loc[i]
        
        # ---------------------------------处理表格1----------------------------------
        # 年龄列
        age_feature = table1_row[4]
        
        # 性别列
        sex_feature = 1 if table1_row[5] == '男' else 0
        
        # 提取6-14列特征
        feature1 = table1_row[6:15].values
        
        # 血压列
        xueya = np.array([eval(k) for k in table1_row[15].split('/')])
        
        # 提取16-22列特征
        feature2 = table1_row[16:23].values
        
        # ---------------------------------处理表格2----------------------------------
        # 得到table2中对应首次检查结果列
        feature3 = table2_row[2: 24].values
        
        # ---------------------------------处理表格3----------------------------------
        # 得到表格3中对应首次检查结果列
        # 首先根据首次检查流水号找到对应的行
        # 查找表格中第1列中特定值所在的一整行
        # print("table3_ED type: {}".format(type(table3_ED)))
        # print("table3_HM type: {}".format(type(table3_HM)))
        
        # target_table3_ED_row_index = np.where(table3_ED[:, 1] == table1_row[3])[0]
        # target_table3_HM_row_index = np.where(table3_HM[:, 1] == table1_row[3])[0]
        
        target_ED_row_data = table3_ED.loc[table3_ED['流水号'] == table1_row[3]]
        target_HM_row_data = table3_HM.loc[table3_HM['流水号'] == table1_row[3]]
                
        # 获取整行数据
        if len(target_ED_row_data) > 0 and len(target_HM_row_data) > 0:
            # 得到对应ED的特征列
            feature4 = target_ED_row_data[2: 33].values.flatten()
            # 得到对应HM的特征列
            feature5 = target_HM_row_data[2: 33].values.flatten()
            
            # 汇总X特征列
            # 将零维数组扩展为一维数组
            age_feature = np.expand_dims(age_feature, axis=0)
            all_feature = np.concatenate((age_feature, [sex_feature], feature1, xueya, feature2, feature3, feature4, feature5))
            all_features.append(all_feature)
            
            # 添加y目标值列
            # 得到mRS目标值作为y值
            mRS_targets.append(table1_row[1])
            
            # print("Row data for the target value {}: {}".format(table1_row[3], target_row_data))
        else:
            # 如果不存在，那么直接跳过
            print("Target value {} not found in the specified column.".format(table1_row[3]))
    
    return all_features, mRS_targets


# 预测mRS函数
def predict_mRS():
    # ------------------------------得到X和y之后，使用模型预测--------------------------
    from sklearn.tree import DecisionTreeClassifier
    from sklearn.model_selection import train_test_split
    from sklearn.metrics import accuracy_score
    
    # -------------------------首先得到前100个病人的数据，用于训练模型------------------
    X, y = get_X_y_for_mRS(100)
    
    # 将特征数据转换成NumPy数组
    X = np.array(X)
    y = np.array(y)
    
    print(X.shape)
    print(y.shape)

    # 创建决策树模型
    model = DecisionTreeClassifier(random_state=42)
    
    # 训练模型
    model.fit(X, y)

    # 在训练集上进行预测（可以修改为测试集）
    y_pred = model.predict(X)

    # 计算准确率（因为使用的是训练数据，这里计算的是训练准确率，实际应用时应该使用独立的测试集）
    accuracy = accuracy_score(y, y_pred)
    print('模型准确率为:', accuracy)
    
    # # 绘图
    # plt.scatter(X, y, label='Original Data')
    # plt.plot(X, y_pred, color='red', label=f'Fitted Curve (degree={degree})')
    # plt.xlabel('X')
    # plt.ylabel('y')
    # plt.legend()
    # plt.savefig("process_models/q3/output/q3a_predict.png")
    # plt.clf()
    
    # -------------------------使用所有160个病人的数据来进行预测得到mRS-----------------------
    # 得到所有病人的数据
    X_all, y_all = get_X_y_for_mRS(160)
    
    # 调用模型得到预测值
    y_predict_all = model.predict(X_all)
    
    # 填充表4答案表格中对应的残差列
    table4_answer = pd.read_excel(os.path.join(root_path, table4_name))
    
    # print(table4_answer)
    
    # 使用 iloc 给特定列赋值
    table4_answer.iloc[2: 162, 8] = y_predict_all
    table4_answer.to_excel(os.path.join(root_path, table4_name), index=False, float_format='%.3f')
    
    print("success")
    

if __name__ == "__main__":
    predict_mRS()








process_models/q3/q3a_predict.py
TheWangYang/2023_huawei_mathematical_modeling
import os
import numpy as np
import pandas as pd
from sklearn.tree import DecisionTreeRegressor
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm


# 题目要求：结合多个特征对mRS进行预测
# 思路：保存前100个病人的特征，然后作为X，同时以mRS为目标值


# 定义存放数据的路径
root_path = "/data1/wyy/projects/_competition/mathmodel_code/data"

# 定义表格名称
table1_name = "1.xlsx"
table2_name = "2.xlsx"
table3_name = "3.xlsx"
# 4-答案表
table4_name = "4.xlsx"
table_additional = "additional.xlsx"


# Read Excel files and handle empty values
table1 = pd.read_excel(os.path.join(root_path, table1_name))
table2 = pd.read_excel(os.path.join(root_path, table2_name))
table3 = pd.ExcelFile(os.path.join(root_path, table3_name))

table3_ED = None
table3_HM = None

# 遍历表格3中每个工作表
for sheet_name in table3.sheet_names:
    # 在这里可以对每个工作表的 DataFrame 进行相应的操作
    # 例如，打印工作表名称和 DataFrame 的头部
    print(f'工作表名称: {sheet_name}')
    
    if sheet_name == "ED":
        table3_ED = pd.read_excel(table3, "ED")
    elif sheet_name == "Hemo":
        table3_HM = pd.read_excel(table3, "Hemo")

# 读取附加表
table_addi = pd.read_excel(os.path.join(root_path, table_additional))


# 预测mRS
def get_X_y_for_mRS(k):
    
    # 定义保存前100个病人对应的mRS数值
    mRS_targets = []
    
    # 汇总特征列
    all_features = []
    
    # 分析前100个病人
    for i in range(k):
        
        # 获得表格中每行数据
        table1_row = table1.loc[i]
        table2_row = table2.loc[i]
        table_addi_row = table_addi.loc[i]
        
        # ---------------------------------处理表格1----------------------------------
        # 年龄列
        age_feature = table1_row[4]
        
        # 性别列
        sex_feature = 1 if table1_row[5] == '男' else 0
        
        # 提取6-14列特征
        feature1 = table1_row[6:15].values
        
        # 血压列
        xueya = np.array([eval(k) for k in table1_row[15].split('/')])
        
        # 提取表格table1中的16-22列特征
        feature2 = table1_row[16:23].values
        
        # ---------------------------------处理表格2----------------------------------
        # 得到table2中对应首次检查结果列
        feature3 = table2_row[2: 24].values
        
        # ---------------------------------处理表格3----------------------------------
        # 得到ED和HM子表数据
        target_ED_row_data = table3_ED.loc[table3_ED['流水号'] == table1_row[3]]
        target_HM_row_data = table3_HM.loc[table3_HM['流水号'] == table1_row[3]]
                
        # 获取整行数据
        if len(target_ED_row_data) > 0 and len(target_HM_row_data) > 0:
            # 得到对应ED的特征列
            feature4 = target_ED_row_data[2: 33].values.flatten()
            # 得到对应HM的特征列
            feature5 = target_HM_row_data[2: 33].values.flatten()
            
            # print("Row data for the target value {}: {}".format(table1_row[3], target_row_data))
        else:
            # 如果不存在，那么直接跳过
            print("Target value {} not found in the specified column.".format(table1_row[3]))
        
        
        # --------------------------拼接其他随访时间表2和表3对饮的特征列-----------------------------
        # -----------------------------第k次随访对应的表2中特征列-----------------------------------
        # 处理随访特征中的空值，将NaN替换为0
        table2_suifang1_feature = np.nan_to_num(table2_row[25: 47].values)
        table2_suifang2_feature = np.nan_to_num(table2_row[48: 70].values)
        table2_suifang3_feature = np.nan_to_num(table2_row[71: 93].values)
        table2_suifang4_feature = np.nan_to_num(table2_row[94: 116].values)
        table2_suifang5_feature = np.nan_to_num(table2_row[117: 139].values)
        table2_suifang6_feature = np.nan_to_num(table2_row[140: 162].values)
        table2_suifang7_feature = np.nan_to_num(table2_row[163: 185].values)
        table2_suifang8_feature = np.nan_to_num(table2_row[186: 208].values)
        
        
        # ---------------------------------第k次随访对应的表3中ED特征列----------------------------------
        table3_ED_suifang1_row_data = table3_ED.loc[table3_ED['流水号'] == table2_row[24]]
        table3_ED_suifang2_row_data = table3_ED.loc[table3_ED['流水号'] == table2_row[47]]
        table3_ED_suifang3_row_data = table3_ED.loc[table3_ED['流水号'] == table2_row[70]]
        table3_ED_suifang4_row_data = table3_ED.loc[table3_ED['流水号'] == table2_row[93]]
        table3_ED_suifang5_row_data = table3_ED.loc[table3_ED['流水号'] == table2_row[116]]
        table3_ED_suifang6_row_data = table3_ED.loc[table3_ED['流水号'] == table2_row[139]]
        table3_ED_suifang7_row_data = table3_ED.loc[table3_ED['流水号'] == table2_row[162]]
        table3_ED_suifang8_row_data = table3_ED.loc[table3_ED['流水号'] == table2_row[185]]
        
        # 处理表3中ED和HM特征列，将NaN替换为0
        table3_ED_suifang1_feature = np.nan_to_num(table3_ED_suifang1_row_data[2: 33].values.flatten())
        table3_ED_suifang2_feature = np.nan_to_num(table3_ED_suifang2_row_data[2: 33].values.flatten())
        table3_ED_suifang3_feature = np.nan_to_num(table3_ED_suifang3_row_data[2: 33].values.flatten())
        table3_ED_suifang4_feature = np.nan_to_num(table3_ED_suifang4_row_data[2: 33].values.flatten())
        table3_ED_suifang5_feature = np.nan_to_num(table3_ED_suifang5_row_data[2: 33].values.flatten())
        table3_ED_suifang6_feature = np.nan_to_num(table3_ED_suifang6_row_data[2: 33].values.flatten())
        table3_ED_suifang7_feature = np.nan_to_num(table3_ED_suifang7_row_data[2: 33].values.flatten())
        table3_ED_suifang8_feature = np.nan_to_num(table3_ED_suifang8_row_data[2: 33].values.flatten())
        
        
        # -----------------------------第k次随访对应的表3中ED特征列-----------------------------------------
        table3_HM_suifang1_row_data = table3_HM.loc[table3_HM['流水号'] == table2_row[24]]
        table3_HM_suifang2_row_data = table3_HM.loc[table3_HM['流水号'] == table2_row[47]]
        table3_HM_suifang3_row_data = table3_HM.loc[table3_HM['流水号'] == table2_row[70]]
        table3_HM_suifang4_row_data = table3_HM.loc[table3_HM['流水号'] == table2_row[93]]
        table3_HM_suifang5_row_data = table3_HM.loc[table3_HM['流水号'] == table2_row[116]]
        table3_HM_suifang6_row_data = table3_HM.loc[table3_HM['流水号'] == table2_row[139]]
        table3_HM_suifang7_row_data = table3_HM.loc[table3_HM['流水号'] == table2_row[162]]
        table3_HM_suifang8_row_data = table3_HM.loc[table3_HM['流水号'] == table2_row[185]]

        # 处理表3中HM特征列，将NaN替换为0
        table3_HM_suifang1_feature = np.nan_to_num(table3_HM_suifang1_row_data[2: 33].values.flatten())
        table3_HM_suifang2_feature = np.nan_to_num(table3_HM_suifang2_row_data[2: 33].values.flatten())
        table3_HM_suifang3_feature = np.nan_to_num(table3_HM_suifang3_row_data[2: 33].values.flatten())
        table3_HM_suifang4_feature = np.nan_to_num(table3_HM_suifang4_row_data[2: 33].values.flatten())
        table3_HM_suifang5_feature = np.nan_to_num(table3_HM_suifang5_row_data[2: 33].values.flatten())
        table3_HM_suifang6_feature = np.nan_to_num(table3_HM_suifang6_row_data[2: 33].values.flatten())
        table3_HM_suifang7_feature = np.nan_to_num(table3_HM_suifang7_row_data[2: 33].values.flatten())
        table3_HM_suifang8_feature = np.nan_to_num(table3_HM_suifang8_row_data[2: 33].values.flatten())
        
        
        # 汇总X特征列
        # 将零维数组扩展为一维数组
        age_feature = np.expand_dims(age_feature, axis=0)
        # all_feature = np.concatenate((age_feature, [sex_feature], feature1, xueya, feature2, feature3, feature4, feature5))
        
        # 合并所有特征
        all_feature = np.concatenate((age_feature, np.array([sex_feature]), feature1, xueya, feature2, feature3, feature4, feature5, 
                                    table3_ED_suifang1_feature, table3_ED_suifang2_feature,
                                    table3_ED_suifang3_feature, table3_ED_suifang4_feature,
                                    table3_ED_suifang5_feature, table3_ED_suifang6_feature,
                                    table3_ED_suifang7_feature, table3_ED_suifang8_feature,
                                    table3_HM_suifang1_feature, table3_HM_suifang2_feature,
                                    table3_HM_suifang3_feature, table3_HM_suifang4_feature,
                                    table3_HM_suifang5_feature, table3_HM_suifang6_feature,
                                    table3_HM_suifang7_feature, table3_HM_suifang8_feature,
                                    table2_suifang1_feature, table2_suifang2_feature,
                                    table2_suifang3_feature, table2_suifang4_feature,
                                    table2_suifang5_feature, table2_suifang6_feature,
                                    table2_suifang7_feature, table2_suifang8_feature))

        all_features.append(all_feature)
        
        # 添加y目标值列
        # 得到mRS目标值作为y值
        mRS_targets.append(table1_row[1])
    
    
    return all_features, mRS_targets


# 预测mRS函数
def predict_mRS():
    # ------------------------------得到X和y之后，使用模型预测--------------------------
    from sklearn.tree import DecisionTreeClassifier
    from sklearn.model_selection import train_test_split
    from sklearn.metrics import accuracy_score
    
    # -------------------------首先得到前100个病人的数据，用于训练模型------------------
    X, y = get_X_y_for_mRS(100)
    
    X = np.array(X, dtype=np.float64)
    y = np.array(y, dtype=np.float64)
    
    X[np.isnan(X)] = 0
    
    print(X.shape)
    print(y.shape)

    # 创建决策树模型
    model = DecisionTreeClassifier(random_state=42)
    
    # 训练模型
    model.fit(X, y)

    # 在训练集上进行预测（可以修改为测试集）
    y_pred = model.predict(X)

    # 计算准确率（因为使用的是训练数据，这里计算的是训练准确率，实际应用时应该使用独立的测试集）
    accuracy = accuracy_score(y, y_pred)
    print('模型准确率为:', accuracy)
    
    # # 绘图
    # plt.scatter(X, y, label='Original Data')
    # plt.plot(X, y_pred, color='red', label=f'Fitted Curve (degree={degree})')
    # plt.xlabel('X')
    # plt.ylabel('y')
    # plt.legend()
    # plt.savefig("process_models/q3/output/q3a_predict.png")
    # plt.clf()
    
    # -------------------------使用所有160个病人的数据来进行预测得到mRS-----------------------
    # 得到所有病人的数据
    X_all, y_all = get_X_y_for_mRS(160)
    
    X_all = np.array(X_all, dtype=np.float64)
    y_all = np.array(y_all, dtype=np.float64)
    
    X_all[np.isnan(X_all)] = 0
    
    # 调用模型得到预测值
    y_predict_all = model.predict(X_all)
    
    # 填充表4答案表格中对应的残差列
    table4_answer = pd.read_excel(os.path.join(root_path, table4_name))
    
    # print(table4_answer)
    
    # 使用 iloc 给特定列赋值
    table4_answer.iloc[2: 162, 9] = y_predict_all
    table4_answer.to_excel(os.path.join(root_path, table4_name), index=False, float_format='%.3f')
    
    print("success")
    

if __name__ == "__main__":
    predict_mRS()





process_models/q3/q3b_predict.py
TheWangYang/2023_huawei_mathematical_modeling