import pandas as pd
import numpy as np
import pyDOE2
import statsmodels.api as sm
from statsmodels.formula.api import ols
import matplotlib.pyplot as plt


# 定义因子名称
factors = [
    "Glucose",       # 葡萄糖
    "Peptone",       # 蛋白胨
    "KH2PO4",        # 磷酸二氢钾
    "Loading_Vol",   # 装液量/250 mL
    "Speed",         # 转速
    "Inoculum",      # 接种量
    "Temp",          # 温度
    "Time"           # 培养时间
]

# 设定高低水平 (+1 和 -1) 的实际物理值
# 这里定义一个字典来存储每个因子的 [低水平, 高水平]
level_map = {
    "Glucose": [20, 40],          # g/L
    "Peptone": [5, 15],           # g/L
    "KH2PO4": [1, 2],             # g/L
    "Loading_Vol": [50, 100],     # mL
    "Speed": [120, 180],          # rpm
    "Inoculum": [5, 15],          # mL
    "Temp": [22, 28],             # d
    "Time": [5, 9]                # ℃
}

# 结果输入到这里
ergosterol_content = [ 3.05, 3.41, 3.42, 3.24, 3.15, 2.44, 3.02, 3.15, 3.58, 2.32, 3.18, 3.95]






def generate_pb_design(factors, level_map):

    # 1. 生成 Plackett-Burman 设计
    # pbdesign(n) 中的 n 是因子个数
    # pyDOE2 会自动选择大于 n 且为 4 的倍数的最小实验次数 N (这里 N=12)
    pb_matrix = pyDOE2.pbdesign(len(factors))
    
    # 2. 转化为 Pandas DataFrame 方便查看
    df = pd.DataFrame(pb_matrix, columns=factors)
    
    # 3. 生成实际操作表
    real_values_df = df.copy()
    for factor in factors:
        low, high = level_map[factor]
        # 将 -1 替换为低水平，1 替换为高水平
        real_values_df[factor] = real_values_df[factor].map({-1: low, 1: high})
    
    # 添加一个空的响应值列（例如：麦角甾醇含量）
    real_values_df['ergosterol'] = np.nan
    
    return df, real_values_df


# -----------  执行函数，生成PB筛选实验设计  --------------------

coded_df, actual_df = generate_pb_design(factors, level_map)
print("--- PB设计编码矩阵 (-1, 1) ---")
print(coded_df)
print("\n--- 实验记录表 (实际物理值) ---")
print(actual_df)

# 导出为Excel供实验记录使用
coded_df.to_excel("PB_Design_coded_data.xlsx", index=False)
actual_df.to_excel("PB_Design_experimental_data.xlsx", index=False)

# --------------------------  END  ------------------------------



# -----------  进行ANOVA显著性分析，并输出报告  --------------------

# 把结果导入到实验矩阵中
coded_df['ergosterol'] = ergosterol_content

# 构建线性回归模型进行 ANOVA 分析
# 公式：响应值 ~ 因子A + 因子B + ...
formula = 'ergosterol ~ ' + ' + '.join(factors)
model = ols(formula, data=coded_df).fit()

# 输出分析报告
print("\n--- ANOVA 显著性分析报告 ---")
print(model.summary())

# --------------------------  END  ------------------------------




# -----------  提取效应值（回归系数）并绘制帕累托图 (Pareto Chart)  ----------------

# 排除截距项 Intercept
effects = model.params.drop('Intercept').abs().sort_values(ascending=True)

plt.figure(figsize=(10, 6))
effects.plot(kind='barh', color='skyblue')
plt.axvline(x=0.1, color='red', linestyle='--') # 假设的显著性阈值线
plt.title('Pareto Chart of Factor Effects on Ergosterol')
plt.xlabel('Absolute Coefficient (Effect Size)')
plt.ylabel('Factors')
plt.grid(axis='x', linestyle='--', alpha=0.7)
plt.show()

# --------------------------  END  ------------------------------