import os
import json
import np
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
from scipy import signal
from scipy.optimize import curve_fit
import pandas as pd
import warnings
import math
warnings.filterwarnings("ignore")
np.random.seed(42)
output_dirs = ['out', 'out/data', 'out/figs', 'out/tables', 'out/models']
for directory in output_dirs:
os.makedirs(directory, exist_ok=True)
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
def load_or_simulate_data(np=None):
"""
读取附件1和附件2的数据,或者在文件不存在时模拟生成数据
返回:两个DataFrame,分别包含波数和反射率
"""
try:
# 尝试读取实际数据
data_10deg = pd.read_excel(r"C:\Users\R7000P\Downloads\CUMCM2025Problems\B题\附件\附件1.xlsx", header=None)
data_15deg = pd.read_excel(r"C:\Users\R7000P\Downloads\CUMCM2025Problems\B题\附件\附件2.xlsx", header=None)
# 重命名列
data_10deg.columns = ['wavenumber', 'reflectivity']
data_15deg.columns = ['wavenumber', 'reflectivity']
print("成功读取实际数据文件")
# 保存一份到输出目录
data_10deg.to_csv(f"out/data/data_10deg_{timestamp}.csv", index=False)
data_15deg.to_csv(f"out/data/data_15deg_{timestamp}.csv", index=False)
return data_10deg, data_15deg
except FileNotFoundError:
print("未找到实际数据文件,生成模拟数据...")
# 这里可以补充模拟生成数据的代码,比如用随机数生成等
# 示例:生成模拟数据
import numpy as np
wavenumber = np.linspace(1000, 2000, 100)
reflectivity_10 = np.random.rand(100)
reflectivity_15 = np.random.rand(100)
data_10deg = pd.DataFrame({'wavenumber': wavenumber, 'reflectivity': reflectivity_10})
data_15deg = pd.DataFrame({'wavenumber': wavenumber, 'reflectivity': reflectivity_15})
# 保存模拟数据
data_10deg.to_csv(f"out/data/sim_data_10deg_{timestamp}.csv", index=False)
data_15deg.to_csv(f"out/data/sim_data_15deg_{timestamp}.csv", index=False)
return data_10deg, data_15deg
except FileNotFoundError:
print ("未找到实际数据文件,生成模拟数据...")
# 模拟碳化硅样本的波数和反射率数据
# 设定参数:厚度约35μm,折射率约2.6
d = 35.0 # 微米
n_eff = 2.6
theta_1_rad = np.deg2rad(10)
theta_2_rad = np.deg2rad(15)
# 生成波数范围 (400-4000 cm^-1)
wavenumbers = np.linspace(400, 4000, 1000)
# 计算光程差
opd_10deg = 2 * d * np.sqrt(n_eff**2 - np.sin(theta_1_rad)**2)
opd_15deg = 2 * d * np.sqrt(n_eff**2 - np.sin(theta_2_rad)**2)
# 生成反射率,模拟干涉条纹
# 使用余弦函数模拟干涉条纹,并加入随机噪声
reflectivity_10deg = 30 + 20 * np.cos(2 * np.pi * wavenumbers * opd_10deg) + np.random.normal(0, 2, len(wavenumbers))
reflectivity_15deg = 30 + 20 * np.cos(2 * np.pi * wavenumbers * opd_15deg) + np.random.normal(0, 2, len(wavenumbers))
# 创建DataFrame
data_10deg = (pd.DataFrame({
'wavenumber': wavenumbers,
'reflectivity': reflectivity_10deg
}))
data_15deg = pd.DataFrame({
'wavenumber': wavenumbers,
'reflectivity': reflectivity_15deg
})
# 保存模拟数据
data_10deg.to_csv(f"out/data/simulated_data_10deg_{timestamp}.csv", index=False)
data_15deg.to_csv(f"out/data/simulated_data_15deg_{timestamp}.csv", index=False)
return data_10deg, data_15deg
def explore_data(data_10deg, data_15deg):
"""
对数据进行基本探索和审计
"""
print("\n数据审计:")
# 打印数据维度
data_10deg = np.array([[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12]])
data_15deg = np.array([[10, 20, 30],
[40, 50, 60]])
print(f"10度入射角数据维度: {data_10deg.shape}")
print(f"15度入射角数据维度: {data_15deg.shape}")
# 打印样例数据
print("\n10度入射角数据样例:")
print(data_10deg.head())
print("\n15度入射角数据样例:")
print(data_15deg.head())
# 基本统计信息
print("\n10度入射角数据统计信息:")
print(data_10deg.describe())
print("\n15度入射角数据统计信息:")
print(data_15deg.describe())
# 检查缺失值
print(f"\n10度入射角数据缺失值: {data_10deg.isnull().sum().sum()}")
print(f"15度入射角数据缺失值: {data_15deg.isnull().sum().sum()}")
# 检查波数范围
print(f"\n10度入射角波数范围: {data_10deg['wavenumber'].min()} - {data_10deg['wavenumber'].max()} cm^-1")
print(f"15度入射角波数范围: {data_15deg['wavenumber'].min()} - {data_15deg['wavenumber'].max()} cm^-1")
# 保存数据统计信息
stats_10deg = data_10deg.describe().to_dict()
stats_15deg = data_15deg.describe().to_dict()
with open(f"out/tables/data_stats_{timestamp}.json", 'w') as f:
json.dump({
"10deg": stats_10deg,
"15deg": stats_15deg
}, f, indent=4)
# 绘制反射率随波数变化的图像
plt.figure(figsize=(10, 6))
plt.plot(data_10deg['wavenumber'], data_10deg['reflectivity'], label='10°入射角')
plt.xlabel('波数 (cm^-1)')
plt.ylabel('反射率 (%)')
plt.title('10°入射角下的反射率光谱')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(f"out/figs/reflectivity_10deg_{timestamp}.png", dpi=300)
plt.figure(figsize=(10, 6))
plt.plot(data_15deg['wavenumber'], data_15deg['reflectivity'], label='15°入射角')
plt.xlabel('波数 (cm^-1)')
plt.ylabel('反射率 (%)')
plt.title('15°入射角下的反射率光谱')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(f"out/figs/reflectivity_15deg_{timestamp}.png", dpi=300)
# 绘制波数分布直方图
plt.figure(figsize=(8, 5))
plt.hist(data_10deg['wavenumber'], bins=50, alpha=0.7)
plt.xlabel('波数 (cm^-1)')
plt.ylabel('频数')
plt.title('10°入射角波数分布')
plt.grid(True)
plt.tight_layout()
plt.savefig(f"out/figs/wavenumber_hist_10deg_{timestamp}.png", dpi=300)
plt.figure(figsize=(8, 5))
plt.hist(data_15deg['wavenumber'], bins=50, alpha=0.7)
plt.xlabel('波数 (cm^-1)')
plt.ylabel('频数')
plt.title('15°入射角波数分布')
plt.grid(True)
plt.tight_layout()
plt.savefig(f"out/figs/wavenumber_hist_15deg_{timestamp}.png", dpi=300)
def build_features(data_10deg, data_15deg):
"""
对数据进行处理,提取干涉条纹的极值点
"""
# 使用Savitzky-Golay滤波器平滑数据
window_length = 21 # 窗口长度
polyorder = 3 # 多项式阶数
# 对10°数据进行平滑
reflectivity_smooth_10deg = signal.savgol_filter(
data_10deg['reflectivity'],
window_length=11,
polyorder=2
)
# 对15°数据进行平滑
reflectivity_smooth_15deg = signal.savgol_filter(
data_15deg['reflectivity'],
window_length=11,
polyorder=2
)
# 创建包含平滑数据的新DataFrame
data_10deg_smooth = data_10deg.copy()
data_10deg_smooth['reflectivity_smooth'] = reflectivity_smooth_10deg
data_15deg_smooth = data_15deg.copy()
data_15deg_smooth['reflectivity_smooth'] = reflectivity_smooth_15deg
# 绘制原始数据和平滑后的数据对比图
plt.figure(figsize=(12, 6))
plt.plot(data_10deg['wavenumber'], data_10deg['reflectivity'], 'b-', alpha=0.5, label='原始数据')
plt.plot(data_10deg_smooth['wavenumber'], data_10deg_smooth['reflectivity_smooth'], 'r-', label='平滑后数据')
plt.xlabel('波数 (cm^-1)')
plt.ylabel('反射率 (%)')
plt.title('10°入射角数据平滑对比')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(f"out/figs/smoothed_data_10deg_{timestamp}.png", dpi=300)
plt.figure(figsize=(12, 6))
plt.plot(data_15deg['wavenumber'], data_15deg['reflectivity'], 'b-', alpha=0.5, label='原始数据')
plt.plot(data_15deg_smooth['wavenumber'], data_15deg_smooth['reflectivity_smooth'], 'r-', label='平滑后数据')
plt.xlabel('波数 (cm^-1)')
plt.ylabel('反射率 (%)')
plt.title('15°入射角数据平滑对比')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(f"out/figs/smoothed_data_15deg_{timestamp}.png", dpi=300)
# 寻找极大值点(反射率峰值)
# 使用SciPy的find_peaks函数
# 设置prominence参数以控制峰的显著性
peaks_10deg, _ = signal.find_peaks(reflectivity_smooth_10deg, prominence=2)
peaks_15deg, _ = signal.find_peaks(reflectivity_smooth_15deg, prominence=2)
# 提取峰值对应的波数和反射率
peak_wavenumbers_10deg = data_10deg['wavenumber'].iloc[peaks_10deg].values
peak_reflectivity_10deg = reflectivity_smooth_10deg[peaks_10deg]
peak_wavenumbers_15deg = data_15deg['wavenumber'].iloc[peaks_15deg].values
peak_reflectivity_15deg = reflectivity_smooth_15deg[peaks_15deg]
# 绘制找到的峰值点
plt.figure(figsize=(12, 6))
plt.plot(data_10deg_smooth['wavenumber'], data_10deg_smooth['reflectivity_smooth'], 'b-')
plt.plot(peak_wavenumbers_10deg, peak_reflectivity_10deg, 'ro', label='峰值点')
plt.xlabel('波数 (cm^-1)')
plt.ylabel('反射率 (%)')
plt.title('10°入射角数据中的峰值点')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(f"out/figs/peak_points_10deg_{timestamp}.png", dpi=300)
plt.figure(figsize=(12, 6))
plt.plot(data_15deg_smooth['wavenumber'], data_15deg_smooth['reflectivity_smooth'], 'b-')
plt.plot(peak_wavenumbers_15deg, peak_reflectivity_15deg, 'ro', label='峰值点')
plt.xlabel('波数 (cm^-1)')
plt.ylabel('反射率 (%)')
plt.title('15°入射角数据中的峰值点')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(f"out/figs/peak_points_15deg_{timestamp}.png", dpi=300)
# 创建包含峰值点信息的DataFrame
peaks_data_10deg = pd.DataFrame({
'peak_index': np.arange(len(peak_wavenumbers_10deg)),
'wavenumber': peak_wavenumbers_10deg,
'reflectivity': peak_reflectivity_10deg
})
peaks_data_15deg = pd.DataFrame({
'peak_index': np.arange(len(peak_wavenumbers_15deg)),
'wavenumber': peak_wavenumbers_15deg,
'reflectivity': peak_reflectivity_15deg
})
# 保存峰值点数据
peaks_data_10deg.to_csv(f"out/data/peak_points_10deg_{timestamp}.csv", index=False)
peaks_data_15deg.to_csv(f"out/data/peak_points_15deg_{timestamp}.csv", index=False)
print("\n峰值点提取完成:")
print(f"10°入射角数据中找到 {len(peaks_data_10deg)} 个峰值点")
print(f"15°入射角数据中找到 {len(peaks_data_15deg)} 个峰值点")
def build_and_evaluate_model(peaks_data_10deg, peaks_data_15deg):
"""
基于双角度测量的极值点线性回归法,计算外延层厚度
"""
print("\n开始建模与求解...")
# 线性回归函数
def linear_fit(x, a, b):
return a * x + b
# 对10°数据进行线性回归
popt_10deg, pcov_10deg = curve_fit(
linear_fit,
peaks_data_10deg['wavenumber'],
peaks_data_10deg['peak_index']
)
# 对15°数据进行线性回归
popt_15deg, pcov_15deg = curve_fit(
linear_fit,
peaks_data_15deg['wavenumber'],
peaks_data_15deg['peak_index']
)
# 提取斜率
S_1 = popt_10deg[0]
S_2 = popt_15deg[0]
print(f"10°入射角数据的斜率 S_1 = {S_1}")
print(f"15°入射角数据的斜率 S_2 = {S_2}")
# 计算厚度d
theta_1_rad = np.deg2rad(10)
theta_2_rad = np.deg2rad(15)
# 使用公式: d = (1/2) * sqrt((S_1^2 - S_2^2) / (sin^2(θ_2) - sin^2(θ_1)))
numerator = S_1**2 - S_2**2
denominator = np.sin(theta_2_rad)**2 - np.sin(theta_1_rad)**2
d = 0.5 * np.sqrt(np.abs(numerator) / denominator) # 使用abs确保开方的值为正
# 计算有效折射率
n_eff_squared = (1/(4*d**2)) * (S_1**2 / (1 - np.sin(theta_1_rad)**2))
n_eff = np.sqrt(n_eff_squared)
print(f"\n计算结果:")
print(f"外延层厚度 d = {d:.6f} cm = {d*10000:.2f} μm")
print(f"有效折射率 n_eff = {n_eff:.4f}")
# 绘制线性回归结果
# 10°数据
plt.figure(figsize=(10, 6))
x_fit = np.linspace(
peaks_data_10deg['wavenumber'].min(),
peaks_data_10deg['wavenumber'].max(),
100
)
y_fit = linear_fit(x_fit, *popt_10deg)
plt.scatter(
peaks_data_10deg['wavenumber'],
peaks_data_10deg['peak_index'],
color='blue',
label='峰值点'
)
plt.plot(
x_fit,
y_fit,
'r-',
label=f'线性拟合: y = {popt_10deg[0]:.8f}x + {popt_10deg[1]:.4f}'
)
plt.xlabel('波数 (cm^-1)')
plt.ylabel('峰值序号')
plt.title('10°入射角数据的线性回归')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(f"out/figs/linear_regression_10deg_{timestamp}.png", dpi=300)
# 15°数据
plt.figure(figsize=(10, 6))
x_fit = np.linspace(
peaks_data_15deg['wavenumber'].min(),
peaks_data_15deg['wavenumber'].max(),
100
)
y_fit = linear_fit(x_fit, *popt_15deg)
plt.scatter(
peaks_data_15deg['wavenumber'],
peaks_data_15deg['peak_index'],
color='blue',
label='峰值点'
)
plt.plot(
x_fit,
y_fit,
'r-',
label=f'线性拟合: y = {popt_15deg[0]:.8f}x + {popt_15deg[1]:.4f}'
)
plt.xlabel('波数 (cm^-1)')
plt.ylabel('峰值序号')
plt.title('15°入射角数据的线性回归')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(f"out/figs/linear_regression_15deg_{timestamp}.png", dpi=300)
# 计算拟合优度 R^2
def r_squared(y_true, y_pred):
ss_tot = np.sum((y_true - np.mean(y_true))**2)
ss_res = np.sum((y_true - y_pred)**2)
return 1 - (ss_res / ss_tot)
y_pred_10deg = linear_fit(peaks_data_10deg['wavenumber'], *popt_10deg)
r2_10deg = r_squared(peaks_data_10deg['peak_index'], y_pred_10deg)
y_pred_15deg = linear_fit(peaks_data_15deg['wavenumber'], *popt_15deg)
r2_15deg = r_squared(peaks_data_15deg['peak_index'], y_pred_15deg)
print(f"\n拟合优度:")
print(f"10°数据 R^2 = {r2_10deg:.6f}")
print(f"15°数据 R^2 = {r2_15deg:.6f}")
# 计算斜率的标准误差
slope_std_err_10deg = np.sqrt(pcov_10deg[0, 0])
slope_std_err_15deg = np.sqrt(pcov_15deg[0, 0])
print(f"\n斜率标准误差:")
print(f"10°数据斜率标准误差 = {slope_std_err_10deg:.8f}")
print(f"15°数据斜率标准误差 = {slope_std_err_15deg:.8f}")
# 计算厚度的不确定度(误差传播)
# 厚度计算公式: d = (1/2) * sqrt((S_1^2 - S_2^2) / (sin^2(θ_2) - sin^2(θ_1)))
# 对斜率的误差进行传播
d_uncertainty = d * np.sqrt(
(2 * S_1 * slope_std_err_10deg)**2 + (2 * S_2 * slope_std_err_15deg)**2
) / (2 * np.abs(S_1**2 - S_2**2))
print(f"\n厚度不确定度:")
print(f"厚度不确定度 = {d_uncertainty:.8f} cm = {d_uncertainty*10000:.2f} μm")
print(f"相对不确定度 = {(d_uncertainty/d)*100:.2f}%")
# 返回结果字典
def calculate_results(d, n_eff, S_1, S_2, r2_10deg, r2_15deg, d_uncertainty):
results = {
"thickness_cm": float(d),
"thickness_um": float(d * 10000),
"effective_refractive_index": float(n_eff),
"slope_10deg": float(S_1),
"slope_15deg": float(S_2),
"r_squared_10deg": float(r2_10deg),
"r_squared_15deg": float(r2_15deg),
"thickness_uncertainty_cm": float(d_uncertainty),
"thickness_uncertainty_um": float(d_uncertainty * 10000),
"relative_uncertainty_percent": float((d_uncertainty / d) * 100)
}
return results
def robustness_checks(data_10deg_smooth, data_15deg_smooth, results):
"""
进行稳健性检验,分析结果对参数变化的敏感性
"""
print("\n开始稳健性检验...")
# 测试不同峰值检测参数对结果的影响
prominence_values = [1.0, 2.0, 3.0, 4.0, 5.0]
thickness_results = []
n_eff_results = []
for prominence in prominence_values:
# 寻找10°数据的峰值点
peaks_10deg, _ = signal.find_peaks(
data_10deg_smooth['reflectivity_smooth'],
prominence=prominence
)
peak_wavenumbers_10deg = data_10deg_smooth['wavenumber'].iloc[peaks_10deg].values
# 寻找15°数据的峰值点
peaks_15deg, _ = signal.find_peaks(
data_15deg_smooth['reflectivity_smooth'],
prominence=prominence
)
peak_wavenumbers_15deg = data_15deg_smooth['wavenumber'].iloc[peaks_15deg].values
# 创建峰值数据DataFrame
peaks_data_10deg = pd.DataFrame({
'peak_index': np.arange(len(peak_wavenumbers_10deg)),
'wavenumber': peak_wavenumbers_10deg
})
peaks_data_15deg = pd.DataFrame({
'peak_index': np.arange(len(peak_wavenumbers_15deg)),
'wavenumber': peak_wavenumbers_15deg
})
# 如果找到的峰值点太少,跳过此次计算
if len(peaks_data_10deg) < 3 or len(peaks_data_15deg) < 3:
thickness_results.append(None)
n_eff_results.append(None)
continue
# 线性回归
def linear_fit(x, a, b):
return a * x + b
try:
# 对10°数据进行线性回归
popt_10deg, _ = curve_fit(
linear_fit,
peaks_data_10deg['wavenumber'],
peaks_data_10deg['peak_index']
)
# 对15°数据进行线性回归
popt_15deg, _ = curve_fit(
linear_fit,
peaks_data_15deg['wavenumber'],
peaks_data_15deg['peak_index']
)
# 提取斜率
S_1 = popt_10deg[0]
S_2 = popt_15deg[0]
# 计算厚度d
theta_1_rad = np.deg2rad(10)
theta_2_rad = np.deg2rad(15)
numerator = S_1**2 - S_2**2
denominator = np.sin(theta_2_rad)**2 - np.sin(theta_1_rad)**2
d = 0.5 * np.sqrt(np.abs(numerator) / denominator)
# 计算有效折射率
n_eff_squared = (1/(4*d**2)) * (S_1**2 / (1 - np.sin(theta_1_rad)**2))
n_eff = np.sqrt(n_eff_squared)
thickness_results.append(d*10000) # 转换为微米
n_eff_results.append(n_eff)
except Exception as e:
print(f"prominence={prominence}时拟合失败: {e}")
thickness_results.append(None)
n_eff_results.append(None)
# 绘制稳健性检验结果
valid_indices = [i for i, val in enumerate(thickness_results) if val is not None]
valid_prominence = [prominence_values[i] for i in valid_indices]
valid_thickness = [thickness_results[i] for i in valid_indices]
valid_n_eff = [n_eff_results[i] for i in valid_indices]
if len(valid_indices) > 1:
plt.figure(figsize=(10, 6))
plt.plot(valid_prominence, valid_thickness, 'bo-')
plt.xlabel('峰值检测参数 (prominence)')
plt.ylabel('计算的厚度 (μm)')
plt.title('峰值检测参数对厚度计算结果的影响')
plt.grid(True)
plt.tight_layout()
plt.savefig(f"out/figs/robustness_thickness_{timestamp}.png", dpi=300)
plt.figure(figsize=(10, 6))
plt.plot(valid_prominence, valid_n_eff, 'ro-')
plt.xlabel('峰值检测参数 (prominence)')
plt.ylabel('计算的有效折射率')
plt.title('峰值检测参数对有效折射率计算结果的影响')
plt.grid(True)
plt.tight_layout()
plt.savefig(f"out/figs/robustness_n_eff_{timestamp}.png", dpi=300)
# 计算结果的稳定性指标
thickness_std = np.std(valid_thickness)
thickness_mean = np.mean(valid_thickness)
thickness_cv = thickness_std / thickness_mean # 变异系数
n_eff_std = np.std(valid_n_eff)
n_eff_mean = np.mean(valid_n_eff)
n_eff_cv = n_eff_std / n_eff_mean # 变异系数
print("\n稳健性检验结果:")
print(f"厚度平均值 = {thickness_mean:.2f} μm")
print(f"厚度标准差 = {thickness_std:.2f} μm")
print(f"厚度变异系数 = {thickness_cv*100:.2f}%")
print(f"有效折射率平均值 = {n_eff_mean:.4f}")
print(f"有效折射率标准差 = {n_eff_std:.4f}")
print(f"有效折射率变异系数 = {n_eff_cv*100:.2f}%")
# 添加稳健性结果到结果字典
some_threshold = 10
def process_robustness_data(thickness_mean, thickness_std, thickness_cv, n_eff_mean, n_eff_std, n_eff_cv,
data_count):
if data_count >= some_threshold: # some_threshold 是你设定的判断数据是否充足的阈值,需要自行定义
robustness_results = {
"thickness_mean_um": float(thickness_mean),
"thickness_std_um": float(thickness_std),
"thickness_cv_percent": float(thickness_cv * 100),
"n_eff_mean": float(n_eff_mean),
"n_eff_std": float(n_eff_std),
"n_eff_cv_percent": float(n_eff_cv * 100)
}
else:
print("稳健性检验数据不足,无法生成图表和统计信息")
robustness_results = {
"note": "稳健性检验数据不足"
}
return robustness_results
def export_results(results, robustness_results):
print("\n导出结果...")
all_results = {
"thickness_results": results,
"robustness_results": robustness_results
}
metrics_path = f"out/metrics_{timestamp}.json"
with open(metrics_path, 'w') as f:
json.dump(all_results, f, indent=4)
results_table = pd.DataFrame([
{
"参数": "外延层厚度",
"数值": f"{results['thickness_um']:.2f} μm",
"不确定度": f"{results['thickness_uncertainty_um']:.2f} μm",
"相对不确定度": f"{results['relative_uncertainty_percent']:.2f}%"
},
{
"参数": "有效折射率",
"数值": f"{results['effective_refractive_index']:.4f}",
"不确定度": "N/A",
"相对不确定度": "N/A"
}
])
table_path = f"out/tables/results_table_{timestamp}.csv"
results_table.to_csv(table_path, index=False)
return metrics_path, table_path
def main():
print("开始执行碳化硅外延层厚度确定的计算流程...")
# 1. 读取或模拟数据
data_10deg, data_15deg = load_or_simulate_data()
# 2. 数据审计
explore_data(data_10deg, data_15deg)
# 3. 特征工程(提取峰值点)
peaks_data_10deg, peaks_data_15deg, data_10deg_smooth, data_15deg_smooth = build_features(data_10deg,
data_15deg)
# 4. 建模与求解
results = build_and_evaluate_model(peaks_data_10deg, peaks_data_15deg)
# 5. 稳健性检验
robustness_results = robustness_checks(data_10deg_smooth, data_15deg_smooth, results)
# 6. 导出结果
metrics_path, table_path = export_results(results, robustness_results)
print("\n计算完成!关键文件的绝对路径:")
print(f"指标文件: {os.path.abspath(metrics_path)}")
print(f"结果表格: {os.path.abspath(table_path)}")
return results
解释代码
最新发布