Compare commits
6 Commits
main
...
fc8fcc2746
| Author | SHA1 | Date | |
|---|---|---|---|
| fc8fcc2746 | |||
| 02e488cd0d | |||
| 484643409d | |||
| dcc023cc04 | |||
|
|
c65e47e6b0 | ||
|
|
05101ac7c7 |
4
.gitignore
vendored
4
.gitignore
vendored
@@ -1,4 +1,4 @@
|
||||
.venv/
|
||||
.venv310/
|
||||
.vscode/
|
||||
src/__pycache__/
|
||||
scripts/__pycache__/
|
||||
old/
|
||||
|
||||
BIN
best_model.pth
BIN
best_model.pth
Binary file not shown.
54
scripts/README.md
Normal file
54
scripts/README.md
Normal file
@@ -0,0 +1,54 @@
|
||||
# 建筑结构简谐激励响应预测 (Harmonic Response Prediction)
|
||||
|
||||
## 项目简介
|
||||
本项目旨在通过数据驱动的方法,预测建筑结构在简谐波底座激振下的顶层加速度响应。工程利用机器学习算法提取核心物理特征,建立激励与响应之间的传递关系,实现对目标结构响应均方根(RMS)的快速、准确预测。
|
||||
|
||||
## 技术路线
|
||||
本项目的核心技术链路包含数据预处理、特征工程、传递率建模与算法回归四个关键环节,逻辑清晰,具备较强的工程可解释性:
|
||||
|
||||
1. **稳态信号截取**
|
||||
由于原始振动数据包含起振与衰减的瞬态过程,系统通过滑动窗口结合变异系数评估,自动截取最平稳的中间段数据,消除非稳态噪声对分析的干扰。
|
||||
2. **物理特征工程**
|
||||
系统从底座激励信号中提取 4 维极具代表性的物理特征作为模型输入:
|
||||
- 频率 ($f$)
|
||||
- 频率的平方 ($f^2$)
|
||||
- 激励均方根 ($x_{rms}$)
|
||||
- 理论加速度幅值 ($(2\pi f)^2 \cdot A$)
|
||||
3. **传递率 (TR) 目标建模**
|
||||
模型将预测目标设定为系统的**传递率 (Transmissibility, TR)**,即输出响应与输入激励的均方根比值 ($TR = y_{rms} / x_{rms}$)。
|
||||
在最终推理阶段,通过公式 $\text{预测响应 } y_{rms} = \text{预测 } TR \times \text{实际输入 } x_{rms}$ 还原最终结果,这种无量纲化的处理极大提升了模型的泛化能力。
|
||||
4. **距离权重 KNN 回归**
|
||||
采用数据标准化 (StandardScaler) 结合距离加权的 K 近邻回归 (KNeighborsRegressor) 算法。通过特征空间中的距离衰减机制,对高维空间中的传递率进行平滑拟合。
|
||||
|
||||
## 核心模块说明
|
||||
|
||||
- **`config.py`**
|
||||
全局配置中心。集中管理数据路径、传感器编号(基座输入与顶层输出)、信号截取比例以及模型超参数。
|
||||
- **`dataset.py`**
|
||||
数据处理引擎。负责读取原始 CSV 振动数据、缺失值插值对齐、稳态窗口搜索、傅里叶主频计算以及核心物理特征的打包提取。
|
||||
- **`train_final.py`**
|
||||
模型训练入口。读取全量简谐波数据,构建特征矩阵与目标向量完成 KNN 模型训练,将最优模型序列化保存,并输出全量数据的拟合评估曲线。
|
||||
- **`evaluate.py`**
|
||||
批量评估工具。加载已保存的模型权重,对全集数据进行预测评估,生成包含理论拟合曲线、散点对比与误差分布的可视化图表。
|
||||
- **`predict_single.py`**
|
||||
单样本预测脚本。支持输入特定的单一 CSV 数据文件,提供端到端的预测功能,并绘制包含时域对比、频域谱线及预测 RMS 对比的综合诊断图。
|
||||
|
||||
## 使用指南
|
||||
|
||||
**1. 训练模型**
|
||||
读取配置文件中指定的数据集,执行全量训练并保存模型权重:
|
||||
```bash
|
||||
python scripts/train_final.py
|
||||
```
|
||||
|
||||
**2. 批量评估**
|
||||
加载已训练的模型,对样本进行批量评估与统计分析:
|
||||
```bash
|
||||
python scripts/evaluate.py
|
||||
```
|
||||
|
||||
**3. 单文件预测**
|
||||
对特定的单一样本进行独立预测,并生成详细的图表报告:
|
||||
```bash
|
||||
python scripts/predict_single.py --file path/to/your/harmonic.csv
|
||||
```
|
||||
BIN
scripts/__pycache__/config.cpython-310.pyc
Normal file
BIN
scripts/__pycache__/config.cpython-310.pyc
Normal file
Binary file not shown.
BIN
scripts/__pycache__/dataset.cpython-310.pyc
Normal file
BIN
scripts/__pycache__/dataset.cpython-310.pyc
Normal file
Binary file not shown.
83
scripts/config.py
Normal file
83
scripts/config.py
Normal file
@@ -0,0 +1,83 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
|
||||
|
||||
CORE_FEATURE_NAMES: tuple[str, ...] = (
|
||||
"frequency_hz",
|
||||
"frequency_squared",
|
||||
"x_rms",
|
||||
"theoretical_accel",
|
||||
)
|
||||
|
||||
|
||||
@dataclass
|
||||
class DataConfig:
|
||||
project_root: Path = field(default_factory=lambda: Path(__file__).resolve().parents[1])
|
||||
data_dir: str = "downloads/Non_TMD"
|
||||
harmonic_pattern: str = "harmonic*.csv"
|
||||
|
||||
code_column: str = "code"
|
||||
time_column: str = "time"
|
||||
base_sensor_code: str = "WSMS00012"
|
||||
base_axis: str = "value1"
|
||||
response_sensor_code: str = "WSMS00007"
|
||||
response_axis: str = "value3"
|
||||
|
||||
middle_segment_start_ratio: float = 0.20
|
||||
middle_segment_end_ratio: float = 0.80
|
||||
min_segment_length: int = 512
|
||||
steady_window_ratio: float = 0.25
|
||||
steady_window_stride_ratio: float = 0.05
|
||||
stability_subwindow_count: int = 4
|
||||
interpolation_method: str = "linear"
|
||||
normalization_eps: float = 1e-6
|
||||
harmonic_amplitude_m: float = 0.005
|
||||
|
||||
data_root: Path = field(init=False)
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
self.project_root = Path(self.project_root).resolve()
|
||||
self.data_root = (self.project_root / self.data_dir).resolve()
|
||||
|
||||
|
||||
@dataclass
|
||||
class ModelConfig:
|
||||
model_name: str = "knn_distance"
|
||||
n_neighbors: int = 4
|
||||
distance_power: int = 2
|
||||
|
||||
|
||||
@dataclass
|
||||
class TrainConfig:
|
||||
checkpoint_dir: str = "checkpoints_final"
|
||||
model_name: str = "task1_final_model.pkl"
|
||||
fit_csv_name: str = "task1_final_fit_all.csv"
|
||||
fit_figure_name: str = "task1_final_fit_curve.png"
|
||||
dense_curve_csv_name: str = "task1_final_dense_curve.csv"
|
||||
evaluation_dense_csv_name: str = "evaluation_dense_curve.csv"
|
||||
dense_curve_points: int = 400
|
||||
|
||||
|
||||
@dataclass
|
||||
class ExperimentConfig:
|
||||
data: DataConfig = field(default_factory=DataConfig)
|
||||
model: ModelConfig = field(default_factory=ModelConfig)
|
||||
train: TrainConfig = field(default_factory=TrainConfig)
|
||||
|
||||
|
||||
def make_experiment_config() -> ExperimentConfig:
|
||||
return ExperimentConfig()
|
||||
|
||||
|
||||
def checkpoint_dir(config: ExperimentConfig) -> Path:
|
||||
path = config.data.project_root / config.train.checkpoint_dir / "task1_final"
|
||||
path.mkdir(parents=True, exist_ok=True)
|
||||
return path
|
||||
|
||||
|
||||
def evaluation_dir(config: ExperimentConfig) -> Path:
|
||||
path = config.data.project_root / "evaluation_outputs" / "task1_final"
|
||||
path.mkdir(parents=True, exist_ok=True)
|
||||
return path
|
||||
294
scripts/dataset.py
Normal file
294
scripts/dataset.py
Normal file
@@ -0,0 +1,294 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
import re
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
try:
|
||||
from .config import CORE_FEATURE_NAMES, DataConfig, ExperimentConfig
|
||||
except ImportError:
|
||||
from config import CORE_FEATURE_NAMES, DataConfig, ExperimentConfig
|
||||
|
||||
|
||||
def calculate_rms(signal: np.ndarray) -> float:
|
||||
signal = np.asarray(signal, dtype=np.float64).reshape(-1)
|
||||
return float(np.sqrt(np.mean(np.square(signal))))
|
||||
|
||||
|
||||
def extract_frequency_hz(file_name: str) -> float | None:
|
||||
match = re.search(r"(\d+(?:\.\d+)?)Hz", file_name, flags=re.IGNORECASE)
|
||||
if match is None:
|
||||
return None
|
||||
return float(match.group(1))
|
||||
|
||||
|
||||
def estimate_sampling_rate(time_values: np.ndarray) -> float:
|
||||
if time_values.size < 2:
|
||||
return 100.0
|
||||
dt = np.diff(time_values)
|
||||
dt = dt[np.isfinite(dt)]
|
||||
dt = dt[dt > 0.0]
|
||||
if dt.size == 0:
|
||||
return 100.0
|
||||
return float(1.0 / np.median(dt))
|
||||
|
||||
|
||||
def get_middle_segment(
|
||||
time_values: np.ndarray,
|
||||
x_values: np.ndarray,
|
||||
y_values: np.ndarray,
|
||||
config: DataConfig,
|
||||
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
|
||||
length = len(x_values)
|
||||
search_start = int(length * config.middle_segment_start_ratio)
|
||||
search_end = int(length * config.middle_segment_end_ratio)
|
||||
search_start = max(0, min(search_start, length - 1))
|
||||
search_end = max(search_start + 1, min(search_end, length))
|
||||
|
||||
if search_end - search_start < config.min_segment_length:
|
||||
center = length // 2
|
||||
half = config.min_segment_length // 2
|
||||
search_start = max(0, center - half)
|
||||
search_end = min(length, search_start + config.min_segment_length)
|
||||
search_start = max(0, search_end - config.min_segment_length)
|
||||
|
||||
search_length = search_end - search_start
|
||||
window_length = max(config.min_segment_length, int(length * config.steady_window_ratio))
|
||||
window_length = min(window_length, search_length)
|
||||
if search_length <= window_length:
|
||||
return (
|
||||
time_values[search_start:search_end],
|
||||
x_values[search_start:search_end],
|
||||
y_values[search_start:search_end],
|
||||
)
|
||||
|
||||
stride = max(1, int(length * config.steady_window_stride_ratio))
|
||||
candidate_y = y_values[search_start:search_end]
|
||||
candidate_rms = calculate_rms(candidate_y)
|
||||
best_score = float("inf")
|
||||
best_slice = slice(search_start, search_start + window_length)
|
||||
|
||||
for window_start in range(search_start, search_end - window_length + 1, stride):
|
||||
window_end = window_start + window_length
|
||||
window_y = y_values[window_start:window_end]
|
||||
split_windows = np.array_split(window_y, config.stability_subwindow_count)
|
||||
split_rms = np.asarray([calculate_rms(chunk) for chunk in split_windows], dtype=np.float64)
|
||||
split_peaks = np.asarray([float(np.max(np.abs(chunk))) for chunk in split_windows], dtype=np.float64)
|
||||
|
||||
rms_cv = float(split_rms.std() / max(split_rms.mean(), config.normalization_eps))
|
||||
peak_cv = float(split_peaks.std() / max(split_peaks.mean(), config.normalization_eps))
|
||||
window_rms = calculate_rms(window_y)
|
||||
|
||||
score = rms_cv + 0.35 * peak_cv - 0.05 * (window_rms / max(candidate_rms, config.normalization_eps))
|
||||
if score < best_score:
|
||||
best_score = score
|
||||
best_slice = slice(window_start, window_end)
|
||||
|
||||
return time_values[best_slice], x_values[best_slice], y_values[best_slice]
|
||||
|
||||
|
||||
def _compute_windowed_spectrum(signal: np.ndarray, sampling_rate: float) -> tuple[np.ndarray, np.ndarray]:
|
||||
signal = np.asarray(signal, dtype=np.float64).reshape(-1)
|
||||
if signal.size < 4:
|
||||
return np.asarray([], dtype=np.float64), np.asarray([], dtype=np.float64)
|
||||
centered = signal - np.mean(signal)
|
||||
window = np.hanning(signal.size)
|
||||
scale = max(np.sum(window), 1e-12)
|
||||
fft_values = np.fft.rfft(centered * window)
|
||||
freqs = np.fft.rfftfreq(signal.size, d=1.0 / sampling_rate)
|
||||
magnitudes = (2.0 / scale) * np.abs(fft_values)
|
||||
return freqs, magnitudes
|
||||
|
||||
|
||||
def _parabolic_peak_frequency(freqs: np.ndarray, magnitudes: np.ndarray, peak_index: int) -> float:
|
||||
if peak_index <= 0 or peak_index >= magnitudes.size - 1:
|
||||
return float(freqs[peak_index])
|
||||
alpha = magnitudes[peak_index - 1]
|
||||
beta = magnitudes[peak_index]
|
||||
gamma = magnitudes[peak_index + 1]
|
||||
denominator = alpha - 2.0 * beta + gamma
|
||||
if abs(denominator) < 1e-12:
|
||||
return float(freqs[peak_index])
|
||||
offset = 0.5 * (alpha - gamma) / denominator
|
||||
bin_width = float(freqs[1] - freqs[0])
|
||||
return float(freqs[peak_index] + offset * bin_width)
|
||||
|
||||
|
||||
def get_dominant_frequency(signal: np.ndarray, sampling_rate: float) -> float:
|
||||
freqs, magnitudes = _compute_windowed_spectrum(signal, sampling_rate)
|
||||
if magnitudes.size == 0:
|
||||
return 0.0
|
||||
magnitudes[0] = 0.0
|
||||
band_mask = (freqs >= 0.1) & (freqs <= 5.0)
|
||||
if not np.any(band_mask):
|
||||
return 0.0
|
||||
band_magnitudes = np.where(band_mask, magnitudes, 0.0)
|
||||
peak_index = int(np.argmax(band_magnitudes))
|
||||
return _parabolic_peak_frequency(freqs, magnitudes, peak_index)
|
||||
|
||||
|
||||
def theoretical_acceleration(frequency_hz: float, amplitude_m: float) -> float:
|
||||
return float((2.0 * np.pi * frequency_hz) ** 2 * amplitude_m)
|
||||
|
||||
|
||||
def extract_core_features(
|
||||
signal: np.ndarray,
|
||||
sampling_rate: float,
|
||||
config: DataConfig,
|
||||
known_frequency_hz: float | None = None,
|
||||
) -> np.ndarray:
|
||||
frequency_hz = known_frequency_hz if known_frequency_hz is not None else get_dominant_frequency(signal, sampling_rate)
|
||||
x_rms = calculate_rms(signal)
|
||||
return np.asarray(
|
||||
[
|
||||
frequency_hz,
|
||||
frequency_hz * frequency_hz,
|
||||
x_rms,
|
||||
theoretical_acceleration(frequency_hz, config.harmonic_amplitude_m),
|
||||
],
|
||||
dtype=np.float32,
|
||||
)
|
||||
|
||||
|
||||
def _value_frame(df: pd.DataFrame, config: DataConfig, sensor_code: str, value_column: str) -> pd.DataFrame:
|
||||
sensor_df = df.loc[df[config.code_column] == sensor_code, [config.time_column, value_column]].copy()
|
||||
sensor_df = sensor_df.sort_values(config.time_column)
|
||||
sensor_df = sensor_df.drop_duplicates(subset=config.time_column, keep="first")
|
||||
sensor_df[config.time_column] = sensor_df[config.time_column].astype("float64")
|
||||
sensor_df[value_column] = sensor_df[value_column].astype("float32")
|
||||
return sensor_df
|
||||
|
||||
|
||||
def load_aligned_signals(file_path: Path, config: DataConfig) -> tuple[np.ndarray, np.ndarray, np.ndarray, int]:
|
||||
df = pd.read_csv(file_path)
|
||||
base_df = _value_frame(df, config, config.base_sensor_code, config.base_axis)
|
||||
response_df = _value_frame(df, config, config.response_sensor_code, config.response_axis)
|
||||
if base_df.empty or response_df.empty:
|
||||
raise ValueError(f"Missing required sensor in {file_path.name}")
|
||||
|
||||
aligned = base_df.rename(columns={config.base_axis: "base_signal"}).merge(
|
||||
response_df.rename(columns={config.response_axis: "response_signal"}),
|
||||
on=config.time_column,
|
||||
how="left",
|
||||
)
|
||||
interpolation_count = int(aligned["response_signal"].isna().sum())
|
||||
aligned["response_signal"] = aligned["response_signal"].interpolate(
|
||||
method=config.interpolation_method,
|
||||
limit_direction="both",
|
||||
).ffill().bfill()
|
||||
if aligned["response_signal"].isna().any():
|
||||
raise ValueError(f"Remaining NaN after interpolation in {file_path.name}")
|
||||
|
||||
time_values = aligned[config.time_column].to_numpy(dtype=np.float64)
|
||||
x_values = aligned["base_signal"].to_numpy(dtype=np.float32)
|
||||
y_values = aligned["response_signal"].to_numpy(dtype=np.float32)
|
||||
return time_values, x_values, y_values, interpolation_count
|
||||
|
||||
|
||||
@dataclass
|
||||
class SampleRecord:
|
||||
file_path: Path
|
||||
frequency_hz: float
|
||||
features: np.ndarray
|
||||
x_rms: float
|
||||
y_rms: float
|
||||
target_tr: float
|
||||
time_middle: np.ndarray
|
||||
x_middle: np.ndarray
|
||||
y_middle: np.ndarray
|
||||
sampling_rate: float
|
||||
interpolation_count: int = 0
|
||||
|
||||
|
||||
@dataclass
|
||||
class LoadReport:
|
||||
loaded_files: list[str] = field(default_factory=list)
|
||||
skipped_files: list[tuple[str, str]] = field(default_factory=list)
|
||||
interpolated_files: dict[str, int] = field(default_factory=dict)
|
||||
|
||||
def to_lines(self) -> list[str]:
|
||||
lines = [f"[all_data] loaded={len(self.loaded_files)} skipped={len(self.skipped_files)}"]
|
||||
for file_name, reason in self.skipped_files:
|
||||
lines.append(f" - skipped {file_name}: {reason}")
|
||||
for file_name, count in self.interpolated_files.items():
|
||||
if count > 0:
|
||||
lines.append(f" - interpolated {file_name}: missing_points={count}")
|
||||
return lines
|
||||
|
||||
|
||||
def list_harmonic_files(config: DataConfig) -> list[Path]:
|
||||
files = sorted(config.data_root.rglob(config.harmonic_pattern))
|
||||
return [file_path for file_path in files if extract_frequency_hz(file_path.name) is not None]
|
||||
|
||||
|
||||
def build_record_from_file(file_path: Path, config: DataConfig) -> SampleRecord:
|
||||
time_values, x_values, y_values, interpolation_count = load_aligned_signals(file_path, config)
|
||||
time_middle, x_middle, y_middle = get_middle_segment(time_values, x_values, y_values, config)
|
||||
sampling_rate = estimate_sampling_rate(time_middle)
|
||||
frequency_hz = extract_frequency_hz(file_path.name)
|
||||
if frequency_hz is None:
|
||||
frequency_hz = get_dominant_frequency(x_middle, sampling_rate)
|
||||
|
||||
x_rms = calculate_rms(x_middle)
|
||||
y_rms = calculate_rms(y_middle)
|
||||
target_tr = float(y_rms / max(x_rms, config.normalization_eps))
|
||||
features = extract_core_features(x_middle, sampling_rate, config, known_frequency_hz=frequency_hz)
|
||||
|
||||
return SampleRecord(
|
||||
file_path=file_path,
|
||||
frequency_hz=float(frequency_hz),
|
||||
features=features,
|
||||
x_rms=float(x_rms),
|
||||
y_rms=float(y_rms),
|
||||
target_tr=target_tr,
|
||||
time_middle=np.asarray(time_middle, dtype=np.float64),
|
||||
x_middle=np.asarray(x_middle, dtype=np.float32).reshape(-1),
|
||||
y_middle=np.asarray(y_middle, dtype=np.float32).reshape(-1),
|
||||
sampling_rate=float(sampling_rate),
|
||||
interpolation_count=interpolation_count,
|
||||
)
|
||||
|
||||
|
||||
def load_all_records(config: ExperimentConfig | DataConfig) -> tuple[list[SampleRecord], LoadReport]:
|
||||
data_config = config.data if isinstance(config, ExperimentConfig) else config
|
||||
records: list[SampleRecord] = []
|
||||
report = LoadReport()
|
||||
|
||||
for file_path in list_harmonic_files(data_config):
|
||||
try:
|
||||
record = build_record_from_file(file_path, data_config)
|
||||
except ValueError as error:
|
||||
report.skipped_files.append((file_path.name, str(error)))
|
||||
continue
|
||||
records.append(record)
|
||||
report.loaded_files.append(file_path.name)
|
||||
if record.interpolation_count > 0:
|
||||
report.interpolated_files[file_path.name] = record.interpolation_count
|
||||
|
||||
if not records:
|
||||
raise RuntimeError("No usable harmonic records found in data_root.")
|
||||
return records, report
|
||||
|
||||
|
||||
def records_to_frame(records: list[SampleRecord]) -> pd.DataFrame:
|
||||
rows = []
|
||||
for record in records:
|
||||
row = {
|
||||
"file_name": record.file_path.name,
|
||||
"frequency_hz": record.frequency_hz,
|
||||
"x_rms": record.x_rms,
|
||||
"y_rms": record.y_rms,
|
||||
"target_tr": record.target_tr,
|
||||
"sampling_rate": record.sampling_rate,
|
||||
}
|
||||
for feature_name, feature_value in zip(CORE_FEATURE_NAMES, record.features.tolist()):
|
||||
row[feature_name] = float(feature_value)
|
||||
rows.append(row)
|
||||
return pd.DataFrame(rows)
|
||||
|
||||
|
||||
def report_to_text(report: LoadReport) -> str:
|
||||
return "\n".join(report.to_lines())
|
||||
141
scripts/evaluate.py
Normal file
141
scripts/evaluate.py
Normal file
@@ -0,0 +1,141 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import pickle
|
||||
from pathlib import Path
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
try:
|
||||
from .config import CORE_FEATURE_NAMES, evaluation_dir, make_experiment_config
|
||||
from .dataset import load_all_records, records_to_frame, report_to_text
|
||||
except ImportError:
|
||||
from config import CORE_FEATURE_NAMES, evaluation_dir, make_experiment_config
|
||||
from dataset import load_all_records, records_to_frame, report_to_text
|
||||
|
||||
|
||||
def parse_args() -> argparse.Namespace:
|
||||
parser = argparse.ArgumentParser(description="Evaluate the final all-data task1 model.")
|
||||
parser.add_argument("--checkpoint", type=str, default=None)
|
||||
return parser.parse_args()
|
||||
|
||||
|
||||
def resolve_checkpoint_path(project_root: Path, checkpoint_arg: str | None) -> Path:
|
||||
if checkpoint_arg:
|
||||
return Path(checkpoint_arg).resolve()
|
||||
return project_root / "checkpoints_final" / "task1_final" / "task1_final_model.pkl"
|
||||
|
||||
|
||||
def relative_error_percent(true_value: float, pred_value: float) -> float:
|
||||
return abs(pred_value - true_value) / max(abs(true_value), 1e-12) * 100.0
|
||||
|
||||
|
||||
def build_dense_curve(model, result_df: pd.DataFrame, config) -> pd.DataFrame:
|
||||
plot_df = result_df.sort_values("frequency_hz").reset_index(drop=True)
|
||||
dense_freq = np.linspace(
|
||||
float(plot_df["frequency_hz"].min()),
|
||||
float(plot_df["frequency_hz"].max()),
|
||||
config.train.dense_curve_points,
|
||||
)
|
||||
dense_x_rms = np.interp(dense_freq, plot_df["frequency_hz"], plot_df["x_rms"])
|
||||
dense_true_rms = np.interp(dense_freq, plot_df["frequency_hz"], plot_df["true_rms"])
|
||||
dense_features = np.column_stack(
|
||||
[
|
||||
dense_freq,
|
||||
dense_freq**2,
|
||||
dense_x_rms,
|
||||
(2.0 * np.pi * dense_freq) ** 2 * config.data.harmonic_amplitude_m,
|
||||
]
|
||||
)
|
||||
dense_pred_tr = model.predict(dense_features).clip(min=config.data.normalization_eps)
|
||||
dense_pred_rms = dense_pred_tr * dense_x_rms
|
||||
return pd.DataFrame(
|
||||
{
|
||||
"frequency_hz": dense_freq,
|
||||
"x_rms_interp": dense_x_rms,
|
||||
"true_rms_interp": dense_true_rms,
|
||||
"pred_tr": dense_pred_tr,
|
||||
"pred_rms": dense_pred_rms,
|
||||
}
|
||||
)
|
||||
|
||||
|
||||
def save_plot(result_df: pd.DataFrame, dense_df: pd.DataFrame, figure_path: Path) -> None:
|
||||
plot_df = result_df.sort_values("frequency_hz").reset_index(drop=True)
|
||||
fig, axes = plt.subplots(2, 1, figsize=(12, 9))
|
||||
fig.suptitle("Task1 Final Model Evaluation | All Harmonic Data")
|
||||
|
||||
axes[0].plot(dense_df["frequency_hz"], dense_df["true_rms_interp"], color="tab:blue", linewidth=2.0, label="True RMS Guide Curve")
|
||||
axes[0].plot(dense_df["frequency_hz"], dense_df["pred_rms"], color="tab:orange", linewidth=2.0, label="Pred RMS Dense Curve")
|
||||
axes[0].scatter(plot_df["frequency_hz"], plot_df["true_rms"], color="tab:blue", s=28, zorder=3, label="True RMS Samples")
|
||||
axes[0].scatter(plot_df["frequency_hz"], plot_df["pred_rms"], color="tab:orange", s=22, zorder=3, label="Pred RMS Samples")
|
||||
axes[0].set_xlabel("Frequency (Hz)")
|
||||
axes[0].set_ylabel("RMS")
|
||||
axes[0].grid(True, alpha=0.3)
|
||||
axes[0].legend()
|
||||
|
||||
axes[1].bar(plot_df["frequency_hz"].astype(str), plot_df["relative_error_percent"], color="tab:orange")
|
||||
axes[1].set_xlabel("Frequency (Hz)")
|
||||
axes[1].set_ylabel("Relative Error (%)")
|
||||
axes[1].grid(True, axis="y", alpha=0.3)
|
||||
axes[1].tick_params(axis="x", labelrotation=45)
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig(figure_path, dpi=180, bbox_inches="tight")
|
||||
plt.close(fig)
|
||||
|
||||
|
||||
def main() -> None:
|
||||
args = parse_args()
|
||||
config = make_experiment_config()
|
||||
ckpt_path = resolve_checkpoint_path(config.data.project_root, args.checkpoint)
|
||||
with ckpt_path.open("rb") as handle:
|
||||
payload = pickle.load(handle)
|
||||
model = payload["model"]
|
||||
|
||||
records, report = load_all_records(config)
|
||||
print(report_to_text(report))
|
||||
data_df = records_to_frame(records)
|
||||
feature_matrix = data_df.loc[:, CORE_FEATURE_NAMES].to_numpy()
|
||||
pred_tr = model.predict(feature_matrix).clip(min=config.data.normalization_eps)
|
||||
pred_rms = pred_tr * data_df["x_rms"].to_numpy()
|
||||
true_rms = data_df["y_rms"].to_numpy()
|
||||
|
||||
result_df = data_df.copy()
|
||||
result_df["pred_tr"] = pred_tr
|
||||
result_df["pred_rms"] = pred_rms
|
||||
result_df["true_rms"] = true_rms
|
||||
result_df["relative_error_percent"] = [
|
||||
relative_error_percent(true, pred) for true, pred in zip(true_rms, pred_rms)
|
||||
]
|
||||
result_df = result_df.sort_values(["frequency_hz", "file_name"]).reset_index(drop=True)
|
||||
|
||||
save_dir = evaluation_dir(config)
|
||||
csv_path = save_dir / "evaluation_all_samples.csv"
|
||||
dense_csv_path = save_dir / config.train.evaluation_dense_csv_name
|
||||
figure_path = save_dir / "evaluation_all_curve.png"
|
||||
dense_df = build_dense_curve(model, result_df, config)
|
||||
result_df.to_csv(csv_path, index=False)
|
||||
dense_df.to_csv(dense_csv_path, index=False)
|
||||
save_plot(result_df, dense_df, figure_path)
|
||||
|
||||
summary = {
|
||||
"count": len(result_df),
|
||||
"mean_error_percent": float(result_df["relative_error_percent"].mean()),
|
||||
"median_error_percent": float(result_df["relative_error_percent"].median()),
|
||||
"max_error_percent": float(result_df["relative_error_percent"].max()),
|
||||
"min_error_percent": float(result_df["relative_error_percent"].min()),
|
||||
}
|
||||
print(f"Checkpoint: {ckpt_path}")
|
||||
print(f"Model: {payload['model_name']}")
|
||||
print(f"Summary: {summary}")
|
||||
print(f"CSV saved to: {csv_path}")
|
||||
print(f"Dense curve CSV saved to: {dense_csv_path}")
|
||||
print(f"Figure saved to: {figure_path}")
|
||||
print(result_df[["file_name", "frequency_hz", "true_rms", "pred_rms", "relative_error_percent"]].to_string(index=False))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
109
scripts/predict_single.py
Normal file
109
scripts/predict_single.py
Normal file
@@ -0,0 +1,109 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import pickle
|
||||
from pathlib import Path
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
try:
|
||||
from .config import CORE_FEATURE_NAMES, evaluation_dir, make_experiment_config
|
||||
from .dataset import build_record_from_file
|
||||
except ImportError:
|
||||
from config import CORE_FEATURE_NAMES, evaluation_dir, make_experiment_config
|
||||
from dataset import build_record_from_file
|
||||
|
||||
|
||||
def parse_args() -> argparse.Namespace:
|
||||
parser = argparse.ArgumentParser(description="Predict task1 RMS from a single harmonic waveform CSV.")
|
||||
parser.add_argument("--file", type=str, required=True)
|
||||
parser.add_argument("--checkpoint", type=str, default=None)
|
||||
return parser.parse_args()
|
||||
|
||||
|
||||
def resolve_checkpoint_path(project_root: Path, checkpoint_arg: str | None) -> Path:
|
||||
if checkpoint_arg:
|
||||
return Path(checkpoint_arg).resolve()
|
||||
return project_root / "checkpoints_final" / "task1_final" / "task1_final_model.pkl"
|
||||
|
||||
|
||||
def save_prediction_figure(file_path: Path, record, pred_rms: float, output_dir: Path) -> Path:
|
||||
time_middle = record.time_middle.reshape(-1)
|
||||
x_middle = record.x_middle.reshape(-1)
|
||||
y_middle = record.y_middle.reshape(-1)
|
||||
sampling_rate = record.sampling_rate
|
||||
fft_values = np.fft.rfft(x_middle - np.mean(x_middle))
|
||||
freqs = np.fft.rfftfreq(x_middle.size, d=1.0 / sampling_rate)
|
||||
amplitudes = np.abs(fft_values)
|
||||
|
||||
fig, axes = plt.subplots(4, 1, figsize=(12, 12))
|
||||
fig.suptitle(f"Task1 Final Single-File Prediction | {file_path.name}")
|
||||
|
||||
axes[0].plot(time_middle, x_middle, color="tab:blue")
|
||||
axes[0].set_title("Input Base Excitation (Middle Segment)")
|
||||
axes[0].set_xlabel("Time")
|
||||
axes[0].set_ylabel("Acceleration")
|
||||
axes[0].grid(True, alpha=0.3)
|
||||
|
||||
axes[1].plot(time_middle, y_middle, color="tab:green")
|
||||
axes[1].set_title("True Top Response (Middle Segment)")
|
||||
axes[1].set_xlabel("Time")
|
||||
axes[1].set_ylabel("Acceleration")
|
||||
axes[1].grid(True, alpha=0.3)
|
||||
|
||||
axes[2].plot(freqs, amplitudes, color="tab:purple")
|
||||
axes[2].axvline(record.frequency_hz, color="tab:red", linestyle="--", label=f"Frequency = {record.frequency_hz:.4f} Hz")
|
||||
axes[2].set_xlim(0.0, 5.0)
|
||||
axes[2].set_title("Input Spectrum")
|
||||
axes[2].set_xlabel("Frequency (Hz)")
|
||||
axes[2].set_ylabel("Amplitude")
|
||||
axes[2].grid(True, alpha=0.3)
|
||||
axes[2].legend()
|
||||
|
||||
axes[3].bar(["True RMS", "Pred RMS"], [record.y_rms, pred_rms], color=["tab:green", "tab:orange"])
|
||||
error_percent = abs(pred_rms - record.y_rms) / max(abs(record.y_rms), 1e-12) * 100.0
|
||||
axes[3].set_title(
|
||||
f"True RMS = {record.y_rms:.6f} | Pred RMS = {pred_rms:.6f} | Error = {error_percent:.2f}%"
|
||||
)
|
||||
axes[3].set_ylabel("RMS")
|
||||
axes[3].grid(True, axis="y", alpha=0.3)
|
||||
|
||||
plt.tight_layout()
|
||||
output_dir.mkdir(parents=True, exist_ok=True)
|
||||
figure_path = output_dir / f"{file_path.stem}_prediction.png"
|
||||
plt.savefig(figure_path, dpi=180, bbox_inches="tight")
|
||||
plt.close(fig)
|
||||
return figure_path
|
||||
|
||||
|
||||
def main() -> None:
|
||||
args = parse_args()
|
||||
config = make_experiment_config()
|
||||
ckpt_path = resolve_checkpoint_path(config.data.project_root, args.checkpoint)
|
||||
with ckpt_path.open("rb") as handle:
|
||||
payload = pickle.load(handle)
|
||||
model = payload["model"]
|
||||
|
||||
file_path = Path(args.file).resolve()
|
||||
record = build_record_from_file(file_path, config.data)
|
||||
pred_tr = max(float(model.predict([record.features.tolist()])[0]), config.data.normalization_eps)
|
||||
pred_rms = pred_tr * record.x_rms
|
||||
|
||||
out_dir = evaluation_dir(config)
|
||||
fig_path = save_prediction_figure(file_path, record, pred_rms, out_dir)
|
||||
|
||||
print(f"Checkpoint: {ckpt_path}")
|
||||
print(f"Model: {payload['model_name']}")
|
||||
print(f"Input file: {file_path}")
|
||||
for feature_name, feature_value in zip(CORE_FEATURE_NAMES, record.features.tolist()):
|
||||
print(f"{feature_name}: {feature_value:.6f}")
|
||||
print(f"Predicted TR: {pred_tr:.6f}")
|
||||
print(f"Predicted RMS: {pred_rms:.6f}")
|
||||
print(f"True RMS: {record.y_rms:.6f}")
|
||||
print(f"Relative RMS Error (%): {abs(pred_rms - record.y_rms) / max(abs(record.y_rms), 1e-12) * 100.0:.4f}")
|
||||
print(f"Figure saved to: {fig_path}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
178
scripts/train_final.py
Normal file
178
scripts/train_final.py
Normal file
@@ -0,0 +1,178 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import pickle
|
||||
from dataclasses import asdict
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
from sklearn.neighbors import KNeighborsRegressor
|
||||
from sklearn.pipeline import Pipeline
|
||||
from sklearn.preprocessing import StandardScaler
|
||||
|
||||
try:
|
||||
from .config import CORE_FEATURE_NAMES, ExperimentConfig, checkpoint_dir, evaluation_dir, make_experiment_config
|
||||
from .dataset import load_all_records, records_to_frame, report_to_text
|
||||
except ImportError:
|
||||
from config import CORE_FEATURE_NAMES, ExperimentConfig, checkpoint_dir, evaluation_dir, make_experiment_config
|
||||
from dataset import load_all_records, records_to_frame, report_to_text
|
||||
|
||||
|
||||
def build_final_model(config: ExperimentConfig) -> Pipeline:
|
||||
return Pipeline(
|
||||
[
|
||||
("scaler", StandardScaler()),
|
||||
(
|
||||
"model",
|
||||
KNeighborsRegressor(
|
||||
n_neighbors=config.model.n_neighbors,
|
||||
weights="distance",
|
||||
p=config.model.distance_power,
|
||||
),
|
||||
),
|
||||
]
|
||||
)
|
||||
|
||||
|
||||
def serialize_for_checkpoint(value: Any) -> Any:
|
||||
if isinstance(value, Path):
|
||||
return str(value)
|
||||
if isinstance(value, dict):
|
||||
return {key: serialize_for_checkpoint(sub_value) for key, sub_value in value.items()}
|
||||
if isinstance(value, tuple):
|
||||
return [serialize_for_checkpoint(item) for item in value]
|
||||
if isinstance(value, list):
|
||||
return [serialize_for_checkpoint(item) for item in value]
|
||||
return value
|
||||
|
||||
|
||||
def relative_error_percent(true_value: float, pred_value: float) -> float:
|
||||
return abs(pred_value - true_value) / max(abs(true_value), 1e-12) * 100.0
|
||||
|
||||
|
||||
def build_dense_curve(model: Pipeline, result_df: pd.DataFrame, config: ExperimentConfig) -> pd.DataFrame:
|
||||
plot_df = result_df.sort_values("frequency_hz").reset_index(drop=True)
|
||||
dense_freq = np.linspace(
|
||||
float(plot_df["frequency_hz"].min()),
|
||||
float(plot_df["frequency_hz"].max()),
|
||||
config.train.dense_curve_points,
|
||||
)
|
||||
dense_x_rms = np.interp(dense_freq, plot_df["frequency_hz"], plot_df["x_rms"])
|
||||
dense_true_rms = np.interp(dense_freq, plot_df["frequency_hz"], plot_df["true_rms"])
|
||||
dense_features = np.column_stack(
|
||||
[
|
||||
dense_freq,
|
||||
dense_freq**2,
|
||||
dense_x_rms,
|
||||
(2.0 * np.pi * dense_freq) ** 2 * config.data.harmonic_amplitude_m,
|
||||
]
|
||||
)
|
||||
dense_pred_tr = model.predict(dense_features).clip(min=config.data.normalization_eps)
|
||||
dense_pred_rms = dense_pred_tr * dense_x_rms
|
||||
return pd.DataFrame(
|
||||
{
|
||||
"frequency_hz": dense_freq,
|
||||
"x_rms_interp": dense_x_rms,
|
||||
"true_rms_interp": dense_true_rms,
|
||||
"pred_tr": dense_pred_tr,
|
||||
"pred_rms": dense_pred_rms,
|
||||
}
|
||||
)
|
||||
|
||||
|
||||
def save_fit_plot(result_df: pd.DataFrame, dense_df: pd.DataFrame, figure_path: Path) -> None:
|
||||
plot_df = result_df.sort_values("frequency_hz").reset_index(drop=True)
|
||||
fig, axes = plt.subplots(2, 1, figsize=(12, 9))
|
||||
fig.suptitle("Task1 Final Model | All Harmonic Data Fit")
|
||||
|
||||
axes[0].plot(dense_df["frequency_hz"], dense_df["true_rms_interp"], color="tab:blue", linewidth=2.0, label="True RMS Guide Curve")
|
||||
axes[0].plot(dense_df["frequency_hz"], dense_df["pred_rms"], color="tab:orange", linewidth=2.0, label="Pred RMS Dense Curve")
|
||||
axes[0].scatter(plot_df["frequency_hz"], plot_df["true_rms"], color="tab:blue", s=28, zorder=3, label="True RMS Samples")
|
||||
axes[0].scatter(plot_df["frequency_hz"], plot_df["pred_rms"], color="tab:orange", s=22, zorder=3, label="Pred RMS Samples")
|
||||
axes[0].set_xlabel("Frequency (Hz)")
|
||||
axes[0].set_ylabel("RMS")
|
||||
axes[0].grid(True, alpha=0.3)
|
||||
axes[0].legend()
|
||||
|
||||
axes[1].bar(plot_df["frequency_hz"].astype(str), plot_df["relative_error_percent"], color="tab:orange")
|
||||
axes[1].set_xlabel("Frequency (Hz)")
|
||||
axes[1].set_ylabel("Relative Error (%)")
|
||||
axes[1].grid(True, axis="y", alpha=0.3)
|
||||
axes[1].tick_params(axis="x", labelrotation=45)
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig(figure_path, dpi=180, bbox_inches="tight")
|
||||
plt.close(fig)
|
||||
|
||||
|
||||
def train_final() -> None:
|
||||
config = make_experiment_config()
|
||||
records, report = load_all_records(config)
|
||||
print(report_to_text(report))
|
||||
|
||||
data_df = records_to_frame(records)
|
||||
feature_matrix = data_df.loc[:, CORE_FEATURE_NAMES].to_numpy()
|
||||
target_tr = data_df["target_tr"].to_numpy()
|
||||
|
||||
model = build_final_model(config)
|
||||
model.fit(feature_matrix, target_tr)
|
||||
|
||||
pred_tr = model.predict(feature_matrix).clip(min=config.data.normalization_eps)
|
||||
pred_rms = pred_tr * data_df["x_rms"].to_numpy()
|
||||
true_rms = data_df["y_rms"].to_numpy()
|
||||
|
||||
result_df = data_df.copy()
|
||||
result_df["pred_tr"] = pred_tr
|
||||
result_df["pred_rms"] = pred_rms
|
||||
result_df["true_rms"] = true_rms
|
||||
result_df["relative_error_percent"] = [
|
||||
relative_error_percent(true, pred) for true, pred in zip(true_rms, pred_rms)
|
||||
]
|
||||
result_df = result_df.sort_values(["frequency_hz", "file_name"]).reset_index(drop=True)
|
||||
|
||||
mean_error = float(result_df["relative_error_percent"].mean())
|
||||
median_error = float(result_df["relative_error_percent"].median())
|
||||
max_error = float(result_df["relative_error_percent"].max())
|
||||
|
||||
ckpt_dir = checkpoint_dir(config)
|
||||
eval_dir = evaluation_dir(config)
|
||||
csv_path = eval_dir / config.train.fit_csv_name
|
||||
fig_path = eval_dir / config.train.fit_figure_name
|
||||
dense_csv_path = eval_dir / config.train.dense_curve_csv_name
|
||||
model_path = ckpt_dir / config.train.model_name
|
||||
|
||||
dense_df = build_dense_curve(model, result_df, config)
|
||||
result_df.to_csv(csv_path, index=False)
|
||||
dense_df.to_csv(dense_csv_path, index=False)
|
||||
save_fit_plot(result_df, dense_df, fig_path)
|
||||
|
||||
payload = {
|
||||
"model_name": config.model.model_name,
|
||||
"model": model,
|
||||
"config": serialize_for_checkpoint(asdict(config)),
|
||||
"feature_names": list(CORE_FEATURE_NAMES),
|
||||
"summary": {
|
||||
"count": len(result_df),
|
||||
"mean_error_percent": mean_error,
|
||||
"median_error_percent": median_error,
|
||||
"max_error_percent": max_error,
|
||||
},
|
||||
}
|
||||
with model_path.open("wb") as handle:
|
||||
pickle.dump(payload, handle)
|
||||
|
||||
print(f"Final model saved to: {model_path}")
|
||||
print(f"Fit CSV saved to: {csv_path}")
|
||||
print(f"Dense curve CSV saved to: {dense_csv_path}")
|
||||
print(f"Fit figure saved to: {fig_path}")
|
||||
print(
|
||||
f"All-data self-fit summary | count={len(result_df)} | "
|
||||
f"mean={mean_error:.4f}% | median={median_error:.4f}% | max={max_error:.4f}%"
|
||||
)
|
||||
print(result_df[["file_name", "frequency_hz", "true_rms", "pred_rms", "relative_error_percent"]].to_string(index=False))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
train_final()
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -1,30 +0,0 @@
|
||||
import os
|
||||
|
||||
# Data Configuration
|
||||
# 使用基于当前文件的绝对路径拼接,以防止你在不同目录下运行报错
|
||||
DATA_DIR = os.path.join(os.path.dirname(os.path.dirname(__file__)), 'downloads')
|
||||
SEQ_LEN = 512 # 滑动窗口的长度 (时间序列步数)
|
||||
STEP_SIZE = 20 # 滑动窗口的步长
|
||||
BATCH_SIZE = 256
|
||||
|
||||
# Features Configuration
|
||||
INPUT_SENSOR = 'WSMS00012'
|
||||
OUTPUT_SENSORS = ['WSMS00007', 'WSMS00008', 'WSMS00009', 'WSMS00010', 'WSMS00011']
|
||||
INPUT_AXIS = 'value1' # 底部传感器输入轴(课程要求:X轴)
|
||||
OUTPUT_AXIS = 'value3' # 目标传感器输出轴
|
||||
|
||||
# Model Configuration
|
||||
CHANNELS = [64, 64, 128, 128, 256, 256] # TCN 各层通道数
|
||||
KERNEL_SIZE = 5
|
||||
DROPOUT = 0.2
|
||||
|
||||
# Training Configuration
|
||||
LEARNING_RATE = 1e-4
|
||||
EPOCHS = 50
|
||||
WEIGHT_DECAY = 1e-3
|
||||
ENABLE_EARLY_STOP = False
|
||||
EARLY_STOP_PATIENCE = 15
|
||||
|
||||
# Device
|
||||
import torch
|
||||
DEVICE = 'cuda' if torch.cuda.is_available() else 'cpu'
|
||||
188
src/dataset.py
188
src/dataset.py
@@ -1,188 +0,0 @@
|
||||
import os
|
||||
import glob
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
import torch
|
||||
from torch.utils.data import Dataset, DataLoader
|
||||
from sklearn.preprocessing import StandardScaler
|
||||
from config import *
|
||||
|
||||
class MultiOutputStandardizer:
|
||||
"""Per-channel standardization that ignores missing labels via masks."""
|
||||
def __init__(self, n_outputs):
|
||||
self.n_outputs = n_outputs
|
||||
self.mean_ = np.zeros(n_outputs, dtype=np.float32)
|
||||
self.scale_ = np.ones(n_outputs, dtype=np.float32)
|
||||
self.fitted = False
|
||||
|
||||
def fit(self, y_sequences, mask_sequences):
|
||||
means = []
|
||||
scales = []
|
||||
for c in range(self.n_outputs):
|
||||
valid_values = []
|
||||
for y_seq, m_seq in zip(y_sequences, mask_sequences):
|
||||
valid = m_seq[:, c] > 0.5
|
||||
if np.any(valid):
|
||||
valid_values.append(y_seq[valid, c])
|
||||
if len(valid_values) == 0:
|
||||
means.append(0.0)
|
||||
scales.append(1.0)
|
||||
continue
|
||||
vals = np.concatenate(valid_values, axis=0)
|
||||
mean = float(np.mean(vals))
|
||||
std = float(np.std(vals))
|
||||
if std < 1e-6:
|
||||
std = 1.0
|
||||
means.append(mean)
|
||||
scales.append(std)
|
||||
|
||||
self.mean_ = np.asarray(means, dtype=np.float32)
|
||||
self.scale_ = np.asarray(scales, dtype=np.float32)
|
||||
self.fitted = True
|
||||
|
||||
def transform(self, y):
|
||||
if not self.fitted:
|
||||
raise RuntimeError("MultiOutputStandardizer must be fitted before transform.")
|
||||
return (y - self.mean_) / self.scale_
|
||||
|
||||
def inverse_transform(self, y):
|
||||
if not self.fitted:
|
||||
raise RuntimeError("MultiOutputStandardizer must be fitted before inverse_transform.")
|
||||
return y * self.scale_ + self.mean_
|
||||
|
||||
|
||||
class BuildingDataset(Dataset):
|
||||
def __init__(self, file_paths, seq_len, step_size, scaler_X=None, scaler_Y=None, fit_scaler=False):
|
||||
self.seq_len = seq_len
|
||||
self.X_data = []
|
||||
self.Y_data = []
|
||||
self.M_data = []
|
||||
self.scaler_X = scaler_X if scaler_X is not None else StandardScaler()
|
||||
self.scaler_Y = scaler_Y if scaler_Y is not None else MultiOutputStandardizer(len(OUTPUT_SENSORS))
|
||||
|
||||
raw_X = []
|
||||
raw_Y = []
|
||||
raw_M = []
|
||||
|
||||
for f in file_paths:
|
||||
# 读取数据
|
||||
df = pd.read_csv(f)
|
||||
# 使用长表格式: code, type, time, value1, value2, value3
|
||||
|
||||
# 提取 012 的输入轴作为基准
|
||||
df_in = df[df['code'] == INPUT_SENSOR][['time', INPUT_AXIS]].rename(columns={INPUT_AXIS: 'input_signal'})
|
||||
if len(df_in) == 0:
|
||||
# 若无输入传感器(如自由衰减数据),则补零
|
||||
df_in = pd.DataFrame({'time': df['time'].unique()})
|
||||
df_in['input_signal'] = 0.0
|
||||
|
||||
# 提取 007~011 的 Z 轴并按时间戳逐步合并 (使用 left join 确保以 df_in 的时间为基准)
|
||||
df_merged = df_in
|
||||
for sens in OUTPUT_SENSORS:
|
||||
df_out_sens = df[df['code'] == sens][['time', OUTPUT_AXIS]].rename(columns={OUTPUT_AXIS: f'out_{sens}'})
|
||||
df_merged = pd.merge(df_merged, df_out_sens, on='time', how='left')
|
||||
|
||||
df_merged = df_merged.sort_values('time').reset_index(drop=True)
|
||||
|
||||
if len(df_merged) == 0:
|
||||
print(f"Warning: Skipping file {f} due to no overlapping timestamps across required sensors.")
|
||||
continue
|
||||
|
||||
x_seq = df_merged['input_signal'].values.reshape(-1, 1).astype(np.float32)
|
||||
|
||||
# 提取所有 target 传感器列与可用性掩码
|
||||
out_cols = [f'out_{sens}' for sens in OUTPUT_SENSORS]
|
||||
y_seq = np.zeros((len(df_merged), len(OUTPUT_SENSORS)), dtype=np.float32)
|
||||
m_seq = np.zeros((len(df_merged), len(OUTPUT_SENSORS)), dtype=np.float32)
|
||||
for c, col in enumerate(out_cols):
|
||||
series = df_merged[col]
|
||||
observed = ~series.isna()
|
||||
m_seq[:, c] = observed.astype(np.float32)
|
||||
if observed.any():
|
||||
filled = series.interpolate(method='linear').bfill().ffill()
|
||||
y_seq[:, c] = filled.fillna(0.0).values.astype(np.float32)
|
||||
else:
|
||||
y_seq[:, c] = 0.0
|
||||
|
||||
raw_X.append(x_seq)
|
||||
raw_Y.append(y_seq)
|
||||
raw_M.append(m_seq)
|
||||
|
||||
if len(raw_X) == 0:
|
||||
raise ValueError("未能从文件中构造出有效序列,请检查数据路径与传感器编码配置。")
|
||||
|
||||
# 拼接所有文件数据进行 fit
|
||||
X_all = np.vstack(raw_X)
|
||||
|
||||
if fit_scaler:
|
||||
self.scaler_X.fit(X_all)
|
||||
self.scaler_Y.fit(raw_Y, raw_M)
|
||||
|
||||
# 切分窗口
|
||||
for x_seq, y_seq, m_seq in zip(raw_X, raw_Y, raw_M):
|
||||
x_seq_scaled = self.scaler_X.transform(x_seq)
|
||||
y_seq_scaled = self.scaler_Y.transform(y_seq)
|
||||
y_seq_scaled = np.where(m_seq > 0.5, y_seq_scaled, 0.0).astype(np.float32)
|
||||
|
||||
for i in range(0, len(x_seq_scaled) - seq_len + 1, step_size):
|
||||
x_win = x_seq_scaled[i:i+seq_len]
|
||||
y_win = y_seq_scaled[i:i+seq_len]
|
||||
m_win = m_seq[i:i+seq_len]
|
||||
if np.sum(m_win) <= 0:
|
||||
continue
|
||||
self.X_data.append(x_win)
|
||||
self.Y_data.append(y_win)
|
||||
self.M_data.append(m_win)
|
||||
|
||||
self.X_data = np.array(self.X_data)
|
||||
self.Y_data = np.array(self.Y_data)
|
||||
self.M_data = np.array(self.M_data)
|
||||
|
||||
def __len__(self):
|
||||
return len(self.X_data)
|
||||
|
||||
def __getitem__(self, idx):
|
||||
return (
|
||||
torch.tensor(self.X_data[idx], dtype=torch.float32),
|
||||
torch.tensor(self.Y_data[idx], dtype=torch.float32),
|
||||
torch.tensor(self.M_data[idx], dtype=torch.float32),
|
||||
)
|
||||
|
||||
def get_dataloaders(condition='Non_TMD', include_free_vib=False):
|
||||
"""
|
||||
condition: 'Non_TMD' 或者是 'TMD'
|
||||
include_free_vib: 是否在训练集中加入自由振动与自由衰减数据
|
||||
"""
|
||||
base_dir = os.path.join(DATA_DIR, condition)
|
||||
|
||||
# 手动区分的子目录
|
||||
train_files = glob.glob(os.path.join(base_dir, 'train', '*.csv'))
|
||||
if include_free_vib:
|
||||
train_files += glob.glob(os.path.join(base_dir, 'free_vib', '*.csv'))
|
||||
|
||||
val_files = glob.glob(os.path.join(base_dir, 'val', '*.csv'))
|
||||
test_files = glob.glob(os.path.join(base_dir, 'test', '*.csv'))
|
||||
|
||||
print(f"[{condition}] Train files: {len(train_files)}, Val files: {len(val_files)}, Test files: {len(test_files)}")
|
||||
|
||||
if len(train_files) == 0:
|
||||
raise ValueError(f"错误: 在 {base_dir}/train 目录下未找到训练文件!请检查路径是否正确。")
|
||||
if len(val_files) == 0:
|
||||
raise ValueError(f"错误: 在 {base_dir}/val 目录下未找到验证文件!请检查路径是否正确。")
|
||||
|
||||
train_dataset = BuildingDataset(train_files, SEQ_LEN, STEP_SIZE, fit_scaler=True)
|
||||
val_dataset = BuildingDataset(val_files, SEQ_LEN, STEP_SIZE,
|
||||
scaler_X=train_dataset.scaler_X,
|
||||
scaler_Y=train_dataset.scaler_Y, fit_scaler=False)
|
||||
# 若某条件(如 TMD)下没有 test 数据,可以处理一下防止报错
|
||||
test_loader = None
|
||||
if len(test_files) > 0:
|
||||
test_dataset = BuildingDataset(test_files, SEQ_LEN, STEP_SIZE,
|
||||
scaler_X=train_dataset.scaler_X,
|
||||
scaler_Y=train_dataset.scaler_Y, fit_scaler=False)
|
||||
test_loader = DataLoader(test_dataset, batch_size=BATCH_SIZE, shuffle=False)
|
||||
|
||||
train_loader = DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True)
|
||||
val_loader = DataLoader(val_dataset, batch_size=BATCH_SIZE, shuffle=False)
|
||||
|
||||
return train_loader, val_loader, test_loader, train_dataset.scaler_X, train_dataset.scaler_Y
|
||||
@@ -1,50 +0,0 @@
|
||||
import torch
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from config import *
|
||||
from dataset import get_dataloaders
|
||||
from model import BuildingTCN
|
||||
|
||||
def evaluate_model():
|
||||
_, _, test_loader, scaler_X, scaler_Y = get_dataloaders()
|
||||
if test_loader is None:
|
||||
raise ValueError("当前数据配置下没有可用的测试集。")
|
||||
|
||||
model = BuildingTCN(input_size=1, output_size=5, num_channels=CHANNELS,
|
||||
kernel_size=KERNEL_SIZE, dropout=DROPOUT).to(DEVICE)
|
||||
model.load_state_dict(torch.load('best_model.pth', map_location=DEVICE))
|
||||
model.eval()
|
||||
|
||||
all_preds = []
|
||||
all_targets = []
|
||||
|
||||
with torch.no_grad():
|
||||
for inputs, targets, masks in test_loader:
|
||||
inputs = inputs.to(DEVICE)
|
||||
outputs = model(inputs)
|
||||
|
||||
all_preds.append(outputs.cpu().numpy())
|
||||
all_targets.append(targets.numpy())
|
||||
|
||||
# Concatenate results
|
||||
all_preds = np.concatenate(all_preds, axis=0)
|
||||
all_targets = np.concatenate(all_targets, axis=0)
|
||||
|
||||
# 取一个 batch 的第一条序列进行可视化
|
||||
sample_pred = scaler_Y.inverse_transform(all_preds[0])
|
||||
sample_target = scaler_Y.inverse_transform(all_targets[0])
|
||||
|
||||
plt.figure(figsize=(15, 10))
|
||||
for i, sens in enumerate(OUTPUT_SENSORS):
|
||||
plt.subplot(5, 1, i+1)
|
||||
plt.plot(sample_target[:, i], label='True', alpha=0.7)
|
||||
plt.plot(sample_pred[:, i], label='Pred', alpha=0.7, linestyle='--')
|
||||
plt.title(f'Sensor {sens} Z-axis Response (Test: Earthquake)')
|
||||
plt.legend()
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig('test_results.png')
|
||||
print("Evaluation done. Result saved to test_results.png")
|
||||
|
||||
if __name__ == '__main__':
|
||||
evaluate_model()
|
||||
76
src/model.py
76
src/model.py
@@ -1,76 +0,0 @@
|
||||
import torch
|
||||
import torch.nn as nn
|
||||
from torch.nn.utils import weight_norm
|
||||
|
||||
class Chomp1d(nn.Module):
|
||||
def __init__(self, chomp_size):
|
||||
super(Chomp1d, self).__init__()
|
||||
self.chomp_size = chomp_size
|
||||
|
||||
def forward(self, x):
|
||||
return x[:, :, :-self.chomp_size].contiguous()
|
||||
|
||||
class TemporalBlock(nn.Module):
|
||||
def __init__(self, n_inputs, n_outputs, kernel_size, stride, dilation, padding, dropout=0.2):
|
||||
super(TemporalBlock, self).__init__()
|
||||
self.conv1 = weight_norm(nn.Conv1d(n_inputs, n_outputs, kernel_size,
|
||||
stride=stride, padding=padding, dilation=dilation))
|
||||
self.chomp1 = Chomp1d(padding)
|
||||
self.relu1 = nn.ReLU()
|
||||
self.dropout1 = nn.Dropout(dropout)
|
||||
|
||||
self.conv2 = weight_norm(nn.Conv1d(n_outputs, n_outputs, kernel_size,
|
||||
stride=stride, padding=padding, dilation=dilation))
|
||||
self.chomp2 = Chomp1d(padding)
|
||||
self.relu2 = nn.ReLU()
|
||||
self.dropout2 = nn.Dropout(dropout)
|
||||
|
||||
self.net = nn.Sequential(self.conv1, self.chomp1, self.relu1, self.dropout1,
|
||||
self.conv2, self.chomp2, self.relu2, self.dropout2)
|
||||
self.downsample = nn.Conv1d(n_inputs, n_outputs, 1) if n_inputs != n_outputs else None
|
||||
self.relu = nn.ReLU()
|
||||
self.init_weights()
|
||||
|
||||
def init_weights(self):
|
||||
self.conv1.weight.data.normal_(0, 0.01)
|
||||
self.conv2.weight.data.normal_(0, 0.01)
|
||||
if self.downsample is not None:
|
||||
self.downsample.weight.data.normal_(0, 0.01)
|
||||
|
||||
def forward(self, x):
|
||||
out = self.net(x)
|
||||
res = x if self.downsample is None else self.downsample(x)
|
||||
return self.relu(out + res)
|
||||
|
||||
class TemporalConvNet(nn.Module):
|
||||
def __init__(self, num_inputs, num_channels, kernel_size=2, dropout=0.2):
|
||||
super(TemporalConvNet, self).__init__()
|
||||
layers = []
|
||||
num_levels = len(num_channels)
|
||||
for i in range(num_levels):
|
||||
dilation_size = 2 ** i
|
||||
in_channels = num_inputs if i == 0 else num_channels[i-1]
|
||||
out_channels = num_channels[i]
|
||||
layers += [TemporalBlock(in_channels, out_channels, kernel_size, stride=1, dilation=dilation_size,
|
||||
padding=(kernel_size-1) * dilation_size, dropout=dropout)]
|
||||
|
||||
self.network = nn.Sequential(*layers)
|
||||
|
||||
def forward(self, x):
|
||||
return self.network(x)
|
||||
|
||||
class BuildingTCN(nn.Module):
|
||||
def __init__(self, input_size, output_size, num_channels, kernel_size=3, dropout=0.2):
|
||||
super(BuildingTCN, self).__init__()
|
||||
self.tcn = TemporalConvNet(input_size, num_channels, kernel_size, dropout=dropout)
|
||||
self.linear = nn.Linear(num_channels[-1], output_size)
|
||||
|
||||
def forward(self, x):
|
||||
# x shape: (batch, seq_len, input_size)
|
||||
# TCN needs shape: (batch, input_size, seq_len)
|
||||
x = x.transpose(1, 2)
|
||||
y = self.tcn(x)
|
||||
# y shape: (batch, num_channels, seq_len)
|
||||
# linear needs shape: (batch, seq_len, num_channels)
|
||||
y = y.transpose(1, 2)
|
||||
return self.linear(y)
|
||||
77
src/train.py
77
src/train.py
@@ -1,77 +0,0 @@
|
||||
import torch
|
||||
import torch.nn as nn
|
||||
import torch.optim as optim
|
||||
from config import *
|
||||
from dataset import get_dataloaders
|
||||
from model import BuildingTCN
|
||||
|
||||
def masked_l1_loss(pred, target, mask):
|
||||
diff = torch.abs(pred - target) * mask
|
||||
denom = mask.sum().clamp(min=1.0)
|
||||
return diff.sum() / denom
|
||||
|
||||
def train_model():
|
||||
train_loader, val_loader, _, _, _ = get_dataloaders()
|
||||
|
||||
# 初始化前向模型 (1 -> 5)
|
||||
model = BuildingTCN(input_size=1, output_size=5, num_channels=CHANNELS,
|
||||
kernel_size=KERNEL_SIZE, dropout=DROPOUT).to(DEVICE)
|
||||
|
||||
optimizer = optim.AdamW(model.parameters(), lr=LEARNING_RATE, weight_decay=WEIGHT_DECAY)
|
||||
scheduler = optim.lr_scheduler.ReduceLROnPlateau(optimizer, 'min', patience=8, factor=0.5)
|
||||
|
||||
best_val_loss = float('inf')
|
||||
no_improve_epochs = 0
|
||||
early_stop_patience = EARLY_STOP_PATIENCE
|
||||
|
||||
for epoch in range(EPOCHS):
|
||||
model.train()
|
||||
train_loss = 0.0
|
||||
|
||||
for inputs, targets, masks in train_loader:
|
||||
inputs = inputs.to(DEVICE)
|
||||
targets = targets.to(DEVICE)
|
||||
masks = masks.to(DEVICE)
|
||||
|
||||
optimizer.zero_grad()
|
||||
outputs = model(inputs)
|
||||
loss = masked_l1_loss(outputs, targets, masks)
|
||||
|
||||
loss.backward()
|
||||
torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)
|
||||
optimizer.step()
|
||||
|
||||
train_loss += loss.item()
|
||||
|
||||
# 验证
|
||||
model.eval()
|
||||
val_loss = 0.0
|
||||
with torch.no_grad():
|
||||
for inputs, targets, masks in val_loader:
|
||||
inputs = inputs.to(DEVICE)
|
||||
targets = targets.to(DEVICE)
|
||||
masks = masks.to(DEVICE)
|
||||
outputs = model(inputs)
|
||||
loss = masked_l1_loss(outputs, targets, masks)
|
||||
val_loss += loss.item()
|
||||
|
||||
train_loss /= len(train_loader)
|
||||
val_loss /= len(val_loader)
|
||||
|
||||
scheduler.step(val_loss)
|
||||
|
||||
print(f"Epoch {epoch+1}/{EPOCHS} | Train Loss: {train_loss:.4f} | Val Loss: {val_loss:.4f}")
|
||||
|
||||
if val_loss < best_val_loss:
|
||||
best_val_loss = val_loss
|
||||
torch.save(model.state_dict(), 'best_model.pth')
|
||||
print(" --> Saved Best Model")
|
||||
no_improve_epochs = 0
|
||||
else:
|
||||
no_improve_epochs += 1
|
||||
if ENABLE_EARLY_STOP and no_improve_epochs >= early_stop_patience:
|
||||
print(f"Early stopping at epoch {epoch+1}")
|
||||
break
|
||||
|
||||
if __name__ == '__main__':
|
||||
train_model()
|
||||
BIN
test_results.png
BIN
test_results.png
Binary file not shown.
|
Before Width: | Height: | Size: 362 KiB |
Reference in New Issue
Block a user