Compare commits
11 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
3d40e88e35 | ||
|
|
d2dde4f557 | ||
|
|
79f3de05f7 | ||
| f3766c74d6 | |||
| 1b3e836d98 | |||
| fc8fcc2746 | |||
| 02e488cd0d | |||
| 484643409d | |||
| dcc023cc04 | |||
|
|
c65e47e6b0 | ||
|
|
05101ac7c7 |
7
.gitignore
vendored
7
.gitignore
vendored
@@ -1,4 +1,7 @@
|
||||
.venv/
|
||||
.venv310/
|
||||
.vscode/
|
||||
src/__pycache__/
|
||||
scripts/__pycache__/
|
||||
old/
|
||||
evaluation_outputs/
|
||||
.studio_cache/
|
||||
.gitignore
|
||||
265
README.md
265
README.md
@@ -1,193 +1,164 @@
|
||||
# 建筑结构响应预测项目
|
||||
# 建筑结构课程评测代码总览 (Course Evaluation Toolkit)
|
||||
|
||||
## 项目概述
|
||||
## 项目简介
|
||||
本项目面向建筑结构振动课程评测场景,提供三类核心功能:
|
||||
|
||||
本项目使用时间卷积网络(Temporal Convolutional Network, TCN)模型来预测建筑结构在地震波和简谐波激励下的响应。模型基于底部传感器的振动数据预测多个传感器的响应,实现结构健康监测和振动控制。
|
||||
1. 根据基底简谐激励频率,预测模型顶部加速度 Z 向响应的整体均方根 RMS 值。
|
||||
2. 根据顶部自由衰减响应,识别结构一阶自振频率和阻尼比。
|
||||
3. 根据顶部强迫振动响应,反推出输入简谐波的激振频率。
|
||||
|
||||
## 主要特性
|
||||
项目包含三套脚本目录:`scripts`、`scripts_2`、`scripts_3`。三套脚本分别对应 RMS 预测、结构动力特性识别、激振频率逆向预测。
|
||||
|
||||
- **多输出预测**:基于单个输入传感器预测5个输出传感器的响应
|
||||
- **时间序列建模**:使用TCN处理序列数据,支持长距离依赖
|
||||
- **数据增强**:支持有/无TMD(调谐质量阻尼器)条件下的数据训练
|
||||
- **自动化评估**:提供完整的训练和测试流程
|
||||
## 课程评测流程
|
||||
教师评测时可按以下流程使用代码:
|
||||
|
||||
1. **顶部响应 RMS 预测**
|
||||
教师进行一组基底激励输入试验。对每组简谐波频率,使用 `scripts` 目录下的脚本预测顶部加速度 Z 向响应的整体 RMS 值。
|
||||
2. **结构动力参数识别**
|
||||
教师对模型略微做出改动后,采集顶部自由衰减响应,使用 `scripts_2` 目录下的脚本识别一阶自振频率和阻尼比。
|
||||
3. **激振频率逆向预测**
|
||||
教师再次进行一组基底激励输入试验,使用 `scripts_3` 目录下的脚本,根据顶部加速度响应反推输入简谐波的激振频率。
|
||||
|
||||
## 项目结构
|
||||
|
||||
```
|
||||
```text
|
||||
Building-main/
|
||||
├── src/
|
||||
│ ├── __init__.py
|
||||
│ ├── config.py # 配置文件,包含模型和训练参数
|
||||
│ ├── dataset.py # 数据集处理和数据加载器
|
||||
│ ├── model.py # TCN模型定义
|
||||
│ ├── train.py # 训练脚本
|
||||
│ └── evaluate.py # 评估脚本
|
||||
├── downloads/ # 数据目录
|
||||
│ ├── Non_TMD/ # 无TMD数据
|
||||
│ │ ├── train/ # 训练数据
|
||||
│ │ ├── val/ # 验证数据
|
||||
│ │ └── test/ # 测试数据
|
||||
│ └── TMD/ # 有TMD数据
|
||||
├── 数据链接.txt # 数据下载链接配置文件
|
||||
├── download_from_data_links.py # 数据下载脚本
|
||||
├── best_model.pth # 训练好的模型权重
|
||||
└── README.md # 项目说明文档
|
||||
│ └── Non_TMD/
|
||||
│ ├── free_vib/ # 自由衰减数据
|
||||
│ ├── train/ # 简谐波训练数据
|
||||
│ ├── val/ # 简谐波验证数据
|
||||
│ └── test/ # 其他测试数据
|
||||
├── scripts/ # 任务 1:顶部响应 RMS 预测
|
||||
├── scripts_2/ # 任务 2:自振频率与阻尼比识别
|
||||
├── scripts_3/ # 任务 3:激振频率逆向预测
|
||||
├── requirement.txt # Python 依赖清单
|
||||
└── README.md # 总说明文档
|
||||
```
|
||||
|
||||
## 环境要求
|
||||
## 三套脚本说明
|
||||
|
||||
- Python 3.8+
|
||||
- PyTorch 1.9+
|
||||
- CUDA(可选,用于GPU加速)
|
||||
- 其他依赖:numpy, pandas, scikit-learn, matplotlib
|
||||
### `scripts`
|
||||
该目录用于完成顶部响应 RMS 预测任务。
|
||||
|
||||
## 安装依赖
|
||||
- 输入为基底传感器信号与文件频率标签。
|
||||
- 输出为顶部加速度 Z 向响应的整体 RMS 预测值。
|
||||
- 主要脚本包括 `train_final.py`、`evaluate.py`、`predict_single.py`。
|
||||
|
||||
1. 创建虚拟环境:
|
||||
```bash
|
||||
python -m venv .venv310
|
||||
# Windows
|
||||
.venv310\Scripts\activate
|
||||
# Linux/Mac
|
||||
source .venv310/bin/activate
|
||||
```
|
||||
### `scripts_2`
|
||||
该目录用于完成结构动力特性识别任务。
|
||||
|
||||
2. 安装依赖:
|
||||
```bash
|
||||
pip install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu118
|
||||
pip install numpy pandas scikit-learn matplotlib
|
||||
```
|
||||
- 数据来源为 `downloads/Non_TMD/free_vib/` 目录下的自由衰减 CSV 文件。
|
||||
- 使用顶部响应传感器 `WSMS00007` 的 `value3` 列。
|
||||
- 输出一阶自振频率、阻尼比与诊断图像。
|
||||
|
||||
## 数据准备
|
||||
### `scripts_3`
|
||||
该目录用于完成激振频率逆向预测任务。
|
||||
|
||||
### 下载数据
|
||||
- 使用强迫振动数据,如 `downloads/Non_TMD/val/` 或 `downloads/Non_TMD/test/`。
|
||||
- 使用顶部响应传感器 `WSMS00007` 的 `value3` 列识别主频。
|
||||
- 主要脚本为 `task4_predict_freq.py`。
|
||||
|
||||
运行数据下载脚本从配置文件下载数据:
|
||||
## 传感器约束
|
||||
|
||||
- **顶部响应传感器**
|
||||
任务 2 和任务 3 统一使用 `WSMS00007` 的 `value3` 列。
|
||||
- **基底输入传感器**
|
||||
`WSMS00012` 的 `value1` 列用于参考与对照。
|
||||
- **任务 3 频率逆推**
|
||||
激振频率预测直接基于顶部响应频谱结果。
|
||||
|
||||
## 环境准备
|
||||
|
||||
项目建议使用 **Python 3.10** 创建新的虚拟环境,再安装根目录的 `requirement.txt`。
|
||||
|
||||
**1. 创建虚拟环境**
|
||||
|
||||
```bash
|
||||
python download_from_data_links.py
|
||||
py -3.10 -m venv .venv
|
||||
```
|
||||
|
||||
脚本会从 `数据链接.txt` 中读取下载链接,并将数据保存到 `downloads/` 目录。
|
||||
|
||||
### 数据格式
|
||||
|
||||
数据为CSV格式的长表结构,包含以下字段:
|
||||
- `code`: 传感器编码
|
||||
- `type`: 数据类型
|
||||
- `time`: 时间戳
|
||||
- `value1`, `value2`, `value3`: 三个轴的振动值
|
||||
|
||||
### 数据集说明
|
||||
|
||||
- **Non_TMD**: 无调谐质量阻尼器条件下的数据
|
||||
- **TMD**: 安装调谐质量阻尼器条件下的数据
|
||||
- 包含谐波激励(不同频率和振幅)和地震波激励数据
|
||||
|
||||
## 配置说明
|
||||
|
||||
主要配置参数在 `src/config.py` 中:
|
||||
|
||||
- **数据配置**:
|
||||
- `SEQ_LEN`: 序列长度(512)
|
||||
- `STEP_SIZE`: 滑动窗口步长(20)
|
||||
- `BATCH_SIZE`: 批大小(256)
|
||||
|
||||
- **传感器配置**:
|
||||
- `INPUT_SENSOR`: 输入传感器编码('WSMS00012')
|
||||
- `OUTPUT_SENSORS`: 输出传感器编码列表
|
||||
- `INPUT_AXIS`: 输入轴('value1' - X轴)
|
||||
- `OUTPUT_AXIS`: 输出轴('value3' - Z轴)
|
||||
|
||||
- **模型配置**:
|
||||
- `CHANNELS`: TCN各层通道数 [64, 64, 128, 128, 256, 256]
|
||||
- `KERNEL_SIZE`: 卷积核大小(5)
|
||||
- `DROPOUT`: Dropout率(0.2)
|
||||
|
||||
- **训练配置**:
|
||||
- `LEARNING_RATE`: 学习率(1e-4)
|
||||
- `EPOCHS`: 训练轮数(50)
|
||||
- `WEIGHT_DECAY`: 权重衰减(1e-3)
|
||||
|
||||
## 训练模型
|
||||
|
||||
### 基本训练
|
||||
|
||||
运行训练脚本:
|
||||
**2. 激活虚拟环境**
|
||||
|
||||
```bash
|
||||
cd src
|
||||
python train.py
|
||||
.\.venv\Scripts\activate
|
||||
```
|
||||
|
||||
训练过程会:
|
||||
1. 加载Non_TMD条件下的训练和验证数据
|
||||
2. 初始化TCN模型
|
||||
3. 使用AdamW优化器和学习率调度器训练
|
||||
4. 保存最佳模型到 `best_model.pth`
|
||||
|
||||
### 高级训练选项
|
||||
|
||||
修改 `src/config.py` 中的参数来自定义训练:
|
||||
- 调整学习率、批大小等超参数
|
||||
- 启用/禁用早停(`ENABLE_EARLY_STOP`)
|
||||
- 修改模型架构参数
|
||||
|
||||
## 模型评估
|
||||
|
||||
运行评估脚本:
|
||||
**3. 安装依赖**
|
||||
|
||||
```bash
|
||||
cd src
|
||||
python evaluate.py
|
||||
python -m pip install -r requirement.txt
|
||||
```
|
||||
|
||||
评估过程会:
|
||||
1. 加载测试数据
|
||||
2. 加载训练好的模型权重
|
||||
3. 对测试数据进行预测
|
||||
4. 生成可视化结果保存为 `test_results.png`
|
||||
`requirement.txt` 已包含三套脚本使用到的核心依赖:
|
||||
|
||||
## 模型架构
|
||||
- `numpy`
|
||||
- `pandas`
|
||||
- `scipy`
|
||||
- `matplotlib`
|
||||
- `scikit-learn`
|
||||
|
||||
### BuildingTCN
|
||||
## 使用方法
|
||||
|
||||
- **输入**: (batch_size, seq_len, 1) - 底部传感器X轴数据
|
||||
- **输出**: (batch_size, seq_len, 5) - 5个传感器的Z轴预测
|
||||
- **架构**: 多层TCN + 线性输出层
|
||||
### 1. 顶部响应 RMS 预测
|
||||
对单个简谐波文件预测顶部响应 RMS:
|
||||
|
||||
### TemporalConvNet
|
||||
```bash
|
||||
python scripts\predict_single.py --file downloads\Non_TMD\val\harmonic_5mm_0.75Hz.csv
|
||||
```
|
||||
|
||||
- 使用膨胀卷积(dilation)实现长距离依赖
|
||||
- 残差连接保证梯度传播
|
||||
- 权重归一化提高训练稳定性
|
||||
对全量样本进行评估:
|
||||
|
||||
## 使用说明
|
||||
```bash
|
||||
python scripts\evaluate.py
|
||||
```
|
||||
|
||||
### 自定义数据集
|
||||
如需重新训练模型:
|
||||
|
||||
1. 将新数据放入 `downloads/Non_TMD/` 或 `downloads/TMD/` 相应子目录
|
||||
2. 确保数据格式符合要求
|
||||
3. 修改 `src/config.py` 中的传感器配置(如需要)
|
||||
4. 重新运行训练脚本
|
||||
```bash
|
||||
python scripts\train_final.py
|
||||
```
|
||||
|
||||
### 预测新数据
|
||||
### 2. 自由衰减法识别自振频率和阻尼比
|
||||
对自由衰减文件执行识别:
|
||||
|
||||
1. 准备输入数据(CSV格式)
|
||||
2. 修改 `src/dataset.py` 中的数据加载逻辑
|
||||
3. 使用加载的模型进行推理
|
||||
```bash
|
||||
python scripts_2\task3_identify.py --file downloads\Non_TMD\free_vib\free_decay.csv
|
||||
```
|
||||
|
||||
## 注意事项
|
||||
脚本输出内容包括:
|
||||
|
||||
- 确保数据目录结构正确
|
||||
- GPU内存不足时可适当减小 `BATCH_SIZE`
|
||||
- 模型收敛可能需要调整学习率
|
||||
- 测试数据可视化仅显示第一个batch的第一条序列
|
||||
- 一阶自振频率 `f_n`
|
||||
- 对数衰减率 `delta`
|
||||
- 阻尼比 `zeta`
|
||||
- 原始波形、滤波结果与包络线图像
|
||||
|
||||
## 贡献
|
||||
### 3. 根据顶部响应反推激振频率
|
||||
对单个强迫振动文件执行频率逆推:
|
||||
|
||||
欢迎提交Issue和Pull Request来改进项目。
|
||||
```bash
|
||||
python scripts_3\task4_predict_freq.py --file downloads\Non_TMD\val\harmonic_5mm_0.75Hz.csv
|
||||
```
|
||||
|
||||
## 许可证
|
||||
脚本输出内容包括:
|
||||
|
||||
本项目仅供学习和研究使用。
|
||||
- 稳态段采样范围
|
||||
- 预测激振频率
|
||||
- 文件名中的频率真值
|
||||
- 绝对误差
|
||||
- 时域与频域诊断图像
|
||||
|
||||
## 结果输出
|
||||
|
||||
运行脚本后,图像与评测结果默认保存到 `evaluation_outputs/` 目录下,对应子目录如下:
|
||||
|
||||
- `evaluation_outputs/task1_final/`
|
||||
- `evaluation_outputs/task3_identify/`
|
||||
- `evaluation_outputs/task4_predict_freq/`
|
||||
|
||||
## 说明
|
||||
|
||||
- `scripts` 目录负责 RMS 预测。
|
||||
- `scripts_2` 目录负责自由衰减识别。
|
||||
- `scripts_3` 目录负责激振频率逆推。
|
||||
- 数据文件命名中包含频率标签时,脚本会自动提取真值用于误差统计。
|
||||
|
||||
BIN
__pycache__/convert_teacher_csv.cpython-314.pyc
Normal file
BIN
__pycache__/convert_teacher_csv.cpython-314.pyc
Normal file
Binary file not shown.
BIN
__pycache__/evaluation_studio.cpython-314.pyc
Normal file
BIN
__pycache__/evaluation_studio.cpython-314.pyc
Normal file
Binary file not shown.
BIN
best_model.pth
BIN
best_model.pth
Binary file not shown.
BIN
checkpoints_final/task1_final/task1_final_model.pkl
Normal file
BIN
checkpoints_final/task1_final/task1_final_model.pkl
Normal file
Binary file not shown.
BIN
checkpoints_final/task1_final/task1_tmd_adapter.pkl
Normal file
BIN
checkpoints_final/task1_final/task1_tmd_adapter.pkl
Normal file
Binary file not shown.
200
convert_teacher_csv.py
Normal file
200
convert_teacher_csv.py
Normal file
@@ -0,0 +1,200 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
from pathlib import Path
|
||||
import re
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
|
||||
VIB_COLUMN_RE = re.compile(r"^(WSMS\d+)\.Acc([XYZ])$", flags=re.IGNORECASE)
|
||||
SGD_COLUMN_RE = re.compile(r"^(WSGD\d+)\.chdata$", flags=re.IGNORECASE)
|
||||
OUTPUT_COLUMNS = ["code", "type", "time", "value1", "value2", "value3"]
|
||||
|
||||
|
||||
def parse_args() -> argparse.Namespace:
|
||||
parser = argparse.ArgumentParser(
|
||||
description="Convert teacher wide-format CSV into downloads-compatible long-format CSV."
|
||||
)
|
||||
parser.add_argument("--input", type=str, required=True, help="Input wide-format CSV path.")
|
||||
parser.add_argument("--output", type=str, required=True, help="Output long-format CSV path.")
|
||||
parser.add_argument(
|
||||
"--time-column",
|
||||
type=str,
|
||||
default="time",
|
||||
help="Time column name in the input CSV (default: time).",
|
||||
)
|
||||
parser.add_argument(
|
||||
"--timezone-offset-hours",
|
||||
type=float,
|
||||
default=8.0,
|
||||
help="Offset used when parsing naive datetime strings to epoch seconds (default: 8).",
|
||||
)
|
||||
return parser.parse_args()
|
||||
|
||||
|
||||
def discover_sensor_columns(columns: list[str]) -> tuple[dict[str, dict[str, str]], dict[str, str]]:
|
||||
vibr_map: dict[str, dict[str, str]] = {}
|
||||
sgd_map: dict[str, str] = {}
|
||||
for column in columns:
|
||||
vibr_match = VIB_COLUMN_RE.match(column)
|
||||
if vibr_match is not None:
|
||||
sensor_code = vibr_match.group(1).upper()
|
||||
axis = vibr_match.group(2).upper()
|
||||
vibr_map.setdefault(sensor_code, {})[axis] = column
|
||||
continue
|
||||
sgd_match = SGD_COLUMN_RE.match(column)
|
||||
if sgd_match is not None:
|
||||
sensor_code = sgd_match.group(1).upper()
|
||||
sgd_map[sensor_code] = column
|
||||
return vibr_map, sgd_map
|
||||
|
||||
|
||||
def _strip_excel_datetime_wrapper(values: pd.Series) -> pd.Series:
|
||||
text = values.astype(str).str.strip()
|
||||
text = text.str.replace(r'^="(.*)"$', r"\1", regex=True)
|
||||
text = text.str.replace(r'^"(.*)"$', r"\1", regex=True)
|
||||
return text
|
||||
|
||||
|
||||
def to_epoch_seconds(time_values: pd.Series, timezone_offset_hours: float) -> pd.Series:
|
||||
numeric_time = pd.to_numeric(time_values, errors="coerce")
|
||||
if numeric_time.notna().mean() > 0.8:
|
||||
return numeric_time.astype("float64")
|
||||
|
||||
text_time = _strip_excel_datetime_wrapper(time_values)
|
||||
dt_values = pd.to_datetime(text_time, errors="coerce")
|
||||
epoch = (dt_values - pd.Timestamp("1970-01-01")) / pd.Timedelta(seconds=1)
|
||||
epoch = epoch - timezone_offset_hours * 3600.0
|
||||
epoch = pd.to_numeric(epoch, errors="coerce")
|
||||
|
||||
# Fall back to numeric parsing for rows that were not parsed as datetime.
|
||||
fallback_numeric = pd.to_numeric(text_time, errors="coerce")
|
||||
epoch = epoch.where(epoch.notna(), fallback_numeric)
|
||||
return epoch.astype("float64")
|
||||
|
||||
|
||||
def _empty_output_frame() -> pd.DataFrame:
|
||||
return pd.DataFrame(columns=OUTPUT_COLUMNS)
|
||||
|
||||
|
||||
def build_vibr_rows(df: pd.DataFrame, time_column: str, vibr_map: dict[str, dict[str, str]]) -> list[pd.DataFrame]:
|
||||
rows: list[pd.DataFrame] = []
|
||||
for sensor_code in sorted(vibr_map.keys()):
|
||||
axis_columns = vibr_map[sensor_code]
|
||||
use_columns = [time_column]
|
||||
if "X" in axis_columns:
|
||||
use_columns.append(axis_columns["X"])
|
||||
if "Y" in axis_columns:
|
||||
use_columns.append(axis_columns["Y"])
|
||||
if "Z" in axis_columns:
|
||||
use_columns.append(axis_columns["Z"])
|
||||
if len(use_columns) == 1:
|
||||
continue
|
||||
|
||||
sensor_df = df.loc[:, use_columns].copy()
|
||||
sensor_df = sensor_df.rename(
|
||||
columns={
|
||||
axis_columns.get("X", ""): "value1",
|
||||
axis_columns.get("Y", ""): "value2",
|
||||
axis_columns.get("Z", ""): "value3",
|
||||
}
|
||||
)
|
||||
for column in ["value1", "value2", "value3"]:
|
||||
if column not in sensor_df.columns:
|
||||
sensor_df[column] = np.nan
|
||||
sensor_df["value1"] = pd.to_numeric(sensor_df["value1"], errors="coerce")
|
||||
sensor_df["value2"] = pd.to_numeric(sensor_df["value2"], errors="coerce")
|
||||
sensor_df["value3"] = pd.to_numeric(sensor_df["value3"], errors="coerce")
|
||||
sensor_df = sensor_df[sensor_df[["value1", "value2", "value3"]].notna().any(axis=1)]
|
||||
if sensor_df.empty:
|
||||
continue
|
||||
|
||||
sensor_df.insert(0, "type", "vibr")
|
||||
sensor_df.insert(0, "code", sensor_code)
|
||||
rows.append(sensor_df[OUTPUT_COLUMNS])
|
||||
return rows
|
||||
|
||||
|
||||
def build_sgd_rows(df: pd.DataFrame, time_column: str, sgd_map: dict[str, str]) -> list[pd.DataFrame]:
|
||||
rows: list[pd.DataFrame] = []
|
||||
for sensor_code in sorted(sgd_map.keys()):
|
||||
sensor_column = sgd_map[sensor_code]
|
||||
sensor_df = df.loc[:, [time_column, sensor_column]].copy()
|
||||
sensor_df = sensor_df.rename(columns={sensor_column: "value1"})
|
||||
sensor_df["value1"] = pd.to_numeric(sensor_df["value1"], errors="coerce")
|
||||
sensor_df = sensor_df[sensor_df["value1"].notna()]
|
||||
if sensor_df.empty:
|
||||
continue
|
||||
|
||||
sensor_df["value2"] = np.nan
|
||||
sensor_df["value3"] = np.nan
|
||||
sensor_df.insert(0, "type", "sgd")
|
||||
sensor_df.insert(0, "code", sensor_code)
|
||||
rows.append(sensor_df[OUTPUT_COLUMNS])
|
||||
return rows
|
||||
|
||||
|
||||
def convert_dataframe(input_df: pd.DataFrame, time_column: str, timezone_offset_hours: float) -> pd.DataFrame:
|
||||
if time_column not in input_df.columns:
|
||||
raise ValueError(f"Time column '{time_column}' not found in input CSV.")
|
||||
|
||||
vibr_map, sgd_map = discover_sensor_columns(input_df.columns.tolist())
|
||||
if not vibr_map and not sgd_map:
|
||||
raise ValueError("No sensor columns detected. Expected WSMSxxxx.AccX/Y/Z or WSGDxxxx.chdata columns.")
|
||||
|
||||
df = input_df.copy()
|
||||
df[time_column] = to_epoch_seconds(df[time_column], timezone_offset_hours=timezone_offset_hours)
|
||||
df = df[df[time_column].notna()]
|
||||
if df.empty:
|
||||
return _empty_output_frame()
|
||||
|
||||
rows: list[pd.DataFrame] = []
|
||||
rows.extend(build_vibr_rows(df, time_column=time_column, vibr_map=vibr_map))
|
||||
rows.extend(build_sgd_rows(df, time_column=time_column, sgd_map=sgd_map))
|
||||
if not rows:
|
||||
return _empty_output_frame()
|
||||
|
||||
output_df = pd.concat(rows, ignore_index=True)
|
||||
output_df["time"] = pd.to_numeric(output_df["time"], errors="coerce").astype("float64")
|
||||
output_df = output_df[output_df["time"].notna()]
|
||||
output_df = output_df.sort_values(["time", "code"], kind="mergesort").reset_index(drop=True)
|
||||
return output_df.loc[:, OUTPUT_COLUMNS]
|
||||
|
||||
|
||||
def print_summary(df: pd.DataFrame, output_path: Path) -> None:
|
||||
print(f"Output saved to: {output_path}")
|
||||
print(f"Rows: {len(df)}")
|
||||
if df.empty:
|
||||
print("No data rows were generated.")
|
||||
return
|
||||
print(f"Sensors: {df['code'].nunique()}")
|
||||
print("Rows by sensor:")
|
||||
counts = df.groupby(["code", "type"]).size().sort_index()
|
||||
for (code, sensor_type), count in counts.items():
|
||||
print(f" - {code} ({sensor_type}): {int(count)}")
|
||||
|
||||
|
||||
def main() -> None:
|
||||
args = parse_args()
|
||||
input_path = Path(args.input).resolve()
|
||||
output_path = Path(args.output).resolve()
|
||||
|
||||
if not input_path.exists():
|
||||
raise FileNotFoundError(f"Input file not found: {input_path}")
|
||||
|
||||
input_df = pd.read_csv(input_path)
|
||||
output_df = convert_dataframe(
|
||||
input_df=input_df,
|
||||
time_column=args.time_column,
|
||||
timezone_offset_hours=float(args.timezone_offset_hours),
|
||||
)
|
||||
|
||||
output_path.parent.mkdir(parents=True, exist_ok=True)
|
||||
output_df.to_csv(output_path, index=False)
|
||||
print_summary(output_df, output_path=output_path)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
63151
downloads/Non_TMD/free_vib/data_converted_free.csv
Normal file
63151
downloads/Non_TMD/free_vib/data_converted_free.csv
Normal file
File diff suppressed because it is too large
Load Diff
63151
downloads/TMD/val/data_converted_tmd.csv
Normal file
63151
downloads/TMD/val/data_converted_tmd.csv
Normal file
File diff suppressed because it is too large
Load Diff
735
evaluation_studio.py
Normal file
735
evaluation_studio.py
Normal file
@@ -0,0 +1,735 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
from datetime import datetime
|
||||
import math
|
||||
import pickle
|
||||
from pathlib import Path
|
||||
import shutil
|
||||
import threading
|
||||
import traceback
|
||||
import tkinter as tk
|
||||
from tkinter import filedialog, messagebox, ttk
|
||||
from tkinter.scrolledtext import ScrolledText
|
||||
|
||||
import matplotlib
|
||||
|
||||
matplotlib.use("Agg")
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
|
||||
PROJECT_ROOT = Path(__file__).resolve().parent
|
||||
CHECKPOINT_DEFAULT = PROJECT_ROOT / "checkpoints_final" / "task1_final" / "task1_final_model.pkl"
|
||||
ADAPTER_DEFAULT = PROJECT_ROOT / "checkpoints_final" / "task1_final" / "task1_tmd_adapter.pkl"
|
||||
APP_OUTPUT_DIR = PROJECT_ROOT / "evaluation_outputs" / "studio_outputs"
|
||||
APP_CACHE_DIR = PROJECT_ROOT / ".studio_cache"
|
||||
|
||||
TASK1_LABEL = "任务1:顶部RMS预测(Non_TMD)"
|
||||
TASK2_LABEL = "任务2:自由衰减识别自振频率和阻尼比"
|
||||
TASK3_LABEL = "任务3:顶部响应反推激振频率"
|
||||
TASK4_LABEL = "挑战任务:TMD RMS预测(基础模型+适配器)"
|
||||
|
||||
TASK_CHOICES = [TASK1_LABEL, TASK2_LABEL, TASK3_LABEL, TASK4_LABEL]
|
||||
|
||||
HAS_DND = False
|
||||
TkRootClass: type[tk.Tk] = tk.Tk
|
||||
try:
|
||||
from tkinterdnd2 import DND_FILES, TkinterDnD # type: ignore
|
||||
|
||||
HAS_DND = True
|
||||
TkRootClass = TkinterDnD.Tk
|
||||
except Exception:
|
||||
DND_FILES = "DND_Files" # type: ignore
|
||||
|
||||
|
||||
def pretty_json_like_lines(payload: dict) -> str:
|
||||
lines: list[str] = []
|
||||
for key in sorted(payload.keys()):
|
||||
value = payload[key]
|
||||
if isinstance(value, float):
|
||||
lines.append(f"- {key}: {value:.6f}")
|
||||
else:
|
||||
lines.append(f"- {key}: {value}")
|
||||
return "\n".join(lines)
|
||||
|
||||
|
||||
def ensure_long_format_csv(input_path: Path, force_convert: bool, tz_offset_hours: float) -> tuple[Path, bool]:
|
||||
from convert_teacher_csv import convert_dataframe
|
||||
|
||||
df_head = pd.read_csv(input_path, nrows=5)
|
||||
lower_columns = {column.lower() for column in df_head.columns}
|
||||
long_required = {"code", "type", "time", "value1", "value2", "value3"}
|
||||
is_long = long_required.issubset(lower_columns)
|
||||
if is_long and not force_convert:
|
||||
return input_path, False
|
||||
|
||||
APP_CACHE_DIR.mkdir(parents=True, exist_ok=True)
|
||||
full_df = pd.read_csv(input_path)
|
||||
converted_df = convert_dataframe(
|
||||
input_df=full_df,
|
||||
time_column="time",
|
||||
timezone_offset_hours=tz_offset_hours,
|
||||
)
|
||||
stamp = datetime.now().strftime("%Y%m%d_%H%M%S")
|
||||
output_path = APP_CACHE_DIR / f"{input_path.stem}_converted_{stamp}.csv"
|
||||
converted_df.to_csv(output_path, index=False)
|
||||
return output_path, True
|
||||
|
||||
|
||||
def _load_sensor_frame(df: pd.DataFrame, sensor_code: str, value_col: str, time_col: str = "time") -> pd.DataFrame:
|
||||
sensor_df = df.loc[df["code"] == sensor_code, [time_col, value_col]].copy()
|
||||
if sensor_df.empty:
|
||||
return sensor_df
|
||||
sensor_df = sensor_df.sort_values(time_col)
|
||||
sensor_df = sensor_df.drop_duplicates(subset=time_col, keep="first")
|
||||
sensor_df[time_col] = pd.to_numeric(sensor_df[time_col], errors="coerce")
|
||||
sensor_df[value_col] = pd.to_numeric(sensor_df[value_col], errors="coerce")
|
||||
sensor_df = sensor_df[sensor_df[time_col].notna() & sensor_df[value_col].notna()]
|
||||
return sensor_df
|
||||
|
||||
|
||||
def build_prediction_input(csv_path: Path, for_tmd: bool) -> dict:
|
||||
from scripts.config import make_experiment_config
|
||||
from scripts.dataset import (
|
||||
calculate_rms,
|
||||
estimate_sampling_rate,
|
||||
extract_core_features,
|
||||
get_dominant_frequency,
|
||||
get_middle_segment,
|
||||
)
|
||||
|
||||
config = make_experiment_config().data
|
||||
if for_tmd:
|
||||
config.data_dir = "downloads/TMD"
|
||||
config.__post_init__()
|
||||
|
||||
df = pd.read_csv(csv_path)
|
||||
base_df = _load_sensor_frame(df, sensor_code=config.base_sensor_code, value_col=config.base_axis, time_col=config.time_column)
|
||||
if base_df.empty:
|
||||
raise ValueError(
|
||||
f"缺少基底传感器 {config.base_sensor_code}/{config.base_axis},"
|
||||
"无法执行RMS预测任务。"
|
||||
)
|
||||
|
||||
response_df = _load_sensor_frame(
|
||||
df,
|
||||
sensor_code=config.response_sensor_code,
|
||||
value_col=config.response_axis,
|
||||
time_col=config.time_column,
|
||||
)
|
||||
|
||||
if response_df.empty:
|
||||
time_values = base_df[config.time_column].to_numpy(dtype=np.float64)
|
||||
x_values = base_df[config.base_axis].to_numpy(dtype=np.float32)
|
||||
time_middle, x_middle, _ = get_middle_segment(time_values, x_values, x_values, config)
|
||||
y_middle = None
|
||||
has_truth = False
|
||||
else:
|
||||
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",
|
||||
)
|
||||
aligned["response_signal"] = aligned["response_signal"].interpolate(method=config.interpolation_method, limit_direction="both")
|
||||
aligned["response_signal"] = aligned["response_signal"].ffill().bfill()
|
||||
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)
|
||||
time_middle, x_middle, y_middle = get_middle_segment(time_values, x_values, y_values, config)
|
||||
has_truth = True
|
||||
|
||||
sampling_rate = estimate_sampling_rate(np.asarray(time_middle, dtype=np.float64))
|
||||
frequency_hz = get_dominant_frequency(np.asarray(x_middle, dtype=np.float32), sampling_rate)
|
||||
features = extract_core_features(np.asarray(x_middle, dtype=np.float32), sampling_rate, config, known_frequency_hz=frequency_hz)
|
||||
x_rms = float(calculate_rms(np.asarray(x_middle, dtype=np.float32)))
|
||||
|
||||
true_y_rms = None
|
||||
if has_truth and y_middle is not None:
|
||||
true_y_rms = float(calculate_rms(np.asarray(y_middle, dtype=np.float32)))
|
||||
|
||||
return {
|
||||
"config": config,
|
||||
"time_middle": np.asarray(time_middle, dtype=np.float64),
|
||||
"x_middle": np.asarray(x_middle, dtype=np.float32),
|
||||
"y_middle": None if y_middle is None else np.asarray(y_middle, dtype=np.float32),
|
||||
"sampling_rate": float(sampling_rate),
|
||||
"frequency_hz": float(frequency_hz),
|
||||
"features": np.asarray(features, dtype=np.float32),
|
||||
"x_rms": x_rms,
|
||||
"true_y_rms": true_y_rms,
|
||||
"has_truth": has_truth,
|
||||
}
|
||||
|
||||
|
||||
def save_prediction_overview_figure(
|
||||
*,
|
||||
file_path: Path,
|
||||
output_dir: Path,
|
||||
title_prefix: str,
|
||||
time_middle: np.ndarray,
|
||||
x_middle: np.ndarray,
|
||||
sampling_rate: float,
|
||||
frequency_hz: float,
|
||||
pred_rms: float,
|
||||
true_rms: float | None,
|
||||
) -> Path:
|
||||
output_dir.mkdir(parents=True, exist_ok=True)
|
||||
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(3, 1, figsize=(10, 8))
|
||||
fig.suptitle(f"{title_prefix} | {file_path.name}")
|
||||
|
||||
axes[0].plot(time_middle, x_middle, color="tab:blue")
|
||||
axes[0].set_title("Input Base Signal (Middle Segment)")
|
||||
axes[0].set_xlabel("Time")
|
||||
axes[0].set_ylabel("Acceleration")
|
||||
axes[0].grid(True, alpha=0.25)
|
||||
|
||||
axes[1].plot(freqs, amplitudes, color="tab:purple")
|
||||
axes[1].axvline(frequency_hz, color="tab:red", linestyle="--", linewidth=1.2, label=f"f={frequency_hz:.4f}Hz")
|
||||
axes[1].set_xlim(0.0, 5.0)
|
||||
axes[1].set_title("Input Spectrum")
|
||||
axes[1].set_xlabel("Frequency (Hz)")
|
||||
axes[1].set_ylabel("Amplitude")
|
||||
axes[1].grid(True, alpha=0.25)
|
||||
axes[1].legend()
|
||||
|
||||
if true_rms is None:
|
||||
axes[2].bar(["Pred RMS"], [pred_rms], color=["tab:orange"])
|
||||
axes[2].set_title(f"Pred RMS = {pred_rms:.6f} (truth sensor not provided)")
|
||||
else:
|
||||
axes[2].bar(["True RMS", "Pred RMS"], [true_rms, pred_rms], color=["tab:green", "tab:orange"])
|
||||
err_pct = abs(pred_rms - true_rms) / max(abs(true_rms), 1e-12) * 100.0
|
||||
axes[2].set_title(f"True={true_rms:.6f} | Pred={pred_rms:.6f} | Error={err_pct:.2f}%")
|
||||
axes[2].set_ylabel("RMS")
|
||||
axes[2].grid(True, axis="y", alpha=0.25)
|
||||
|
||||
plt.tight_layout()
|
||||
figure_path = output_dir / f"{file_path.stem}_{title_prefix.replace(' ', '_').lower()}.png"
|
||||
plt.savefig(figure_path, dpi=170, bbox_inches="tight")
|
||||
plt.close(fig)
|
||||
return figure_path
|
||||
|
||||
|
||||
def run_task1(csv_path: Path) -> tuple[dict, Path]:
|
||||
from scripts.config import CORE_FEATURE_NAMES
|
||||
from scripts.inference_utils import predict_task1_tr
|
||||
|
||||
with CHECKPOINT_DEFAULT.open("rb") as handle:
|
||||
payload = pickle.load(handle)
|
||||
model = payload["model"]
|
||||
pred_input = build_prediction_input(csv_path=csv_path, for_tmd=False)
|
||||
pred_tr, tr_source = predict_task1_tr(
|
||||
model,
|
||||
feature_vector=pred_input["features"],
|
||||
frequency_hz=float(pred_input["frequency_hz"]),
|
||||
normalization_eps=float(pred_input["config"].normalization_eps),
|
||||
project_root=PROJECT_ROOT,
|
||||
prefer_curve=True,
|
||||
)
|
||||
pred_rms = pred_tr * float(pred_input["x_rms"])
|
||||
true_rms = pred_input["true_y_rms"]
|
||||
rel_err_pct = None if true_rms is None else abs(pred_rms - true_rms) / max(abs(true_rms), 1e-12) * 100.0
|
||||
|
||||
out_dir = APP_OUTPUT_DIR / "task1"
|
||||
fig_path = save_prediction_overview_figure(
|
||||
file_path=csv_path,
|
||||
output_dir=out_dir,
|
||||
title_prefix="Task1 Prediction",
|
||||
time_middle=pred_input["time_middle"],
|
||||
x_middle=pred_input["x_middle"],
|
||||
sampling_rate=float(pred_input["sampling_rate"]),
|
||||
frequency_hz=float(pred_input["frequency_hz"]),
|
||||
pred_rms=float(pred_rms),
|
||||
true_rms=true_rms,
|
||||
)
|
||||
result = {
|
||||
"task": "task1_rms_prediction",
|
||||
"file_path": str(csv_path),
|
||||
"frequency_hz": float(pred_input["frequency_hz"]),
|
||||
"x_rms": float(pred_input["x_rms"]),
|
||||
"true_y_rms": None if true_rms is None else float(true_rms),
|
||||
"pred_tr": float(pred_tr),
|
||||
"pred_tr_source": str(tr_source),
|
||||
"pred_y_rms": float(pred_rms),
|
||||
"relative_error_percent": None if rel_err_pct is None else float(rel_err_pct),
|
||||
"truth_available": bool(pred_input["has_truth"]),
|
||||
"model_name": str(payload.get("model_name", "unknown")),
|
||||
"feature_names": ",".join(CORE_FEATURE_NAMES),
|
||||
}
|
||||
return result, Path(fig_path)
|
||||
|
||||
|
||||
def run_task2(csv_path: Path) -> tuple[dict, Path]:
|
||||
from scripts_2.task3_identify import analyze_free_vibration
|
||||
|
||||
real_csv = csv_path
|
||||
if "free_vib" not in {part.lower() for part in csv_path.parts}:
|
||||
temp_dir = APP_CACHE_DIR / "free_vib"
|
||||
temp_dir.mkdir(parents=True, exist_ok=True)
|
||||
real_csv = temp_dir / csv_path.name
|
||||
shutil.copyfile(csv_path, real_csv)
|
||||
|
||||
output_dir = APP_OUTPUT_DIR / "task2"
|
||||
result = analyze_free_vibration(
|
||||
file_path=real_csv,
|
||||
output_dir=output_dir,
|
||||
bandwidth_hz=0.3,
|
||||
show=False,
|
||||
)
|
||||
fig_path = Path(str(result["figure_path"]))
|
||||
return result, fig_path
|
||||
|
||||
|
||||
def run_task3(csv_path: Path) -> tuple[dict, Path]:
|
||||
from scripts.config import make_experiment_config
|
||||
from scripts_2.signal_utils import (
|
||||
TOP_RESPONSE_AXIS,
|
||||
TOP_RESPONSE_SENSOR,
|
||||
calculate_rms,
|
||||
get_middle_segment,
|
||||
read_sensor_signal,
|
||||
)
|
||||
from scripts_3.task4_predict_freq import analyze_forced_vibration
|
||||
|
||||
output_dir = APP_OUTPUT_DIR / "task3"
|
||||
result = analyze_forced_vibration(
|
||||
file_path=csv_path,
|
||||
output_dir=output_dir,
|
||||
show=False,
|
||||
)
|
||||
|
||||
# Optional inverse estimation:
|
||||
# infer base RMS from top-response RMS with fixed-point iteration:
|
||||
# x = y_rms / TR_model(f, x)
|
||||
try:
|
||||
config = make_experiment_config()
|
||||
with CHECKPOINT_DEFAULT.open("rb") as handle:
|
||||
payload = pickle.load(handle)
|
||||
model = payload["model"]
|
||||
|
||||
time_values, raw_signal = read_sensor_signal(
|
||||
csv_path,
|
||||
sensor_code=TOP_RESPONSE_SENSOR,
|
||||
value_column=TOP_RESPONSE_AXIS,
|
||||
)
|
||||
_, signal_segment = get_middle_segment(time_values, raw_signal)
|
||||
top_rms = float(calculate_rms(signal_segment))
|
||||
freq_hz = float(result["predicted_frequency_hz"])
|
||||
|
||||
min_x = config.data.normalization_eps
|
||||
x_est = max(top_rms, min_x)
|
||||
converged = False
|
||||
rel_change = float("inf")
|
||||
iter_count = 0
|
||||
|
||||
for i in range(1, 41):
|
||||
feat = np.asarray(
|
||||
[
|
||||
[
|
||||
freq_hz,
|
||||
freq_hz * freq_hz,
|
||||
x_est,
|
||||
(2.0 * np.pi * freq_hz) ** 2 * config.data.harmonic_amplitude_m,
|
||||
]
|
||||
],
|
||||
dtype=np.float64,
|
||||
)
|
||||
tr_est = max(float(model.predict(feat)[0]), config.data.normalization_eps)
|
||||
x_new = max(top_rms / max(tr_est, config.data.normalization_eps), min_x)
|
||||
if not np.isfinite(x_new):
|
||||
break
|
||||
rel_change = abs(x_new - x_est) / max(abs(x_est), min_x)
|
||||
# Damping update to improve stability in fixed-point iteration.
|
||||
x_est = 0.6 * x_est + 0.4 * x_new
|
||||
iter_count = i
|
||||
if rel_change < 1e-4:
|
||||
converged = True
|
||||
break
|
||||
|
||||
feat_final = np.asarray(
|
||||
[
|
||||
[
|
||||
freq_hz,
|
||||
freq_hz * freq_hz,
|
||||
x_est,
|
||||
(2.0 * np.pi * freq_hz) ** 2 * config.data.harmonic_amplitude_m,
|
||||
]
|
||||
],
|
||||
dtype=np.float64,
|
||||
)
|
||||
tr_final = max(float(model.predict(feat_final)[0]), config.data.normalization_eps)
|
||||
y_reconstructed = tr_final * x_est
|
||||
residual_pct = abs(y_reconstructed - top_rms) / max(abs(top_rms), min_x) * 100.0
|
||||
|
||||
result["top_response_rms"] = float(top_rms)
|
||||
result["estimated_base_rms"] = float(x_est)
|
||||
result["estimated_tr_from_model"] = float(tr_final)
|
||||
result["inversion_converged"] = bool(converged)
|
||||
result["inversion_iterations"] = int(iter_count)
|
||||
result["inversion_relative_change"] = float(rel_change)
|
||||
result["inversion_residual_percent"] = float(residual_pct)
|
||||
result["inversion_note"] = "Estimated by fixed-point inversion using task1 TR model; use as an estimate only."
|
||||
except Exception as error:
|
||||
result["inversion_note"] = f"Base RMS inversion unavailable: {error}"
|
||||
|
||||
fig_path = Path(str(result["figure_path"]))
|
||||
return result, fig_path
|
||||
|
||||
|
||||
def run_task4_tmd(csv_path: Path) -> tuple[dict, Path]:
|
||||
from scripts_4.adapter import load_adapter
|
||||
from scripts.inference_utils import predict_task1_tr
|
||||
|
||||
with CHECKPOINT_DEFAULT.open("rb") as handle:
|
||||
payload = pickle.load(handle)
|
||||
model = payload["model"]
|
||||
adapter, adapter_extra = load_adapter(ADAPTER_DEFAULT)
|
||||
pred_input = build_prediction_input(csv_path=csv_path, for_tmd=True)
|
||||
pred_tr, tr_source = predict_task1_tr(
|
||||
model,
|
||||
feature_vector=pred_input["features"],
|
||||
frequency_hz=float(pred_input["frequency_hz"]),
|
||||
normalization_eps=float(pred_input["config"].normalization_eps),
|
||||
project_root=PROJECT_ROOT,
|
||||
prefer_curve=True,
|
||||
)
|
||||
pred_base = pred_tr * float(pred_input["x_rms"])
|
||||
pred_tmd = adapter.predict(pred_base, float(pred_input["frequency_hz"]))
|
||||
true_rms = pred_input["true_y_rms"]
|
||||
rel_err_pct = None if true_rms is None else abs(pred_tmd - true_rms) / max(abs(true_rms), 1e-12) * 100.0
|
||||
|
||||
out_dir = APP_OUTPUT_DIR / "task4"
|
||||
fig_path = save_prediction_overview_figure(
|
||||
file_path=csv_path,
|
||||
output_dir=out_dir,
|
||||
title_prefix="Task4 TMD Prediction",
|
||||
time_middle=pred_input["time_middle"],
|
||||
x_middle=pred_input["x_middle"],
|
||||
sampling_rate=float(pred_input["sampling_rate"]),
|
||||
frequency_hz=float(pred_input["frequency_hz"]),
|
||||
pred_rms=float(pred_tmd),
|
||||
true_rms=true_rms,
|
||||
)
|
||||
result = {
|
||||
"task": "task4_tmd_rms_prediction",
|
||||
"file_path": str(csv_path),
|
||||
"frequency_hz": float(pred_input["frequency_hz"]),
|
||||
"x_rms": float(pred_input["x_rms"]),
|
||||
"true_y_rms": None if true_rms is None else float(true_rms),
|
||||
"pred_base_y_rms": float(pred_base),
|
||||
"pred_tr_source": str(tr_source),
|
||||
"pred_tmd_y_rms": float(pred_tmd),
|
||||
"adapter_scale_k": float(adapter.scale_at(float(pred_input["frequency_hz"]))),
|
||||
"relative_error_percent": None if rel_err_pct is None else float(rel_err_pct),
|
||||
"truth_available": bool(pred_input["has_truth"]),
|
||||
"model_name": str(payload.get("model_name", "unknown")),
|
||||
"adapter_extra": str(adapter_extra),
|
||||
}
|
||||
return result, Path(fig_path)
|
||||
|
||||
|
||||
def build_key_metrics(task_label: str, result: dict) -> str:
|
||||
if task_label == TASK1_LABEL:
|
||||
lines = [
|
||||
f"预测RMS: {result.get('pred_y_rms')}",
|
||||
f"预测TR: {result.get('pred_tr')}",
|
||||
f"TR来源: {result.get('pred_tr_source')}",
|
||||
f"激振频率(Hz): {result.get('frequency_hz')}",
|
||||
f"输入RMS(x): {result.get('x_rms')}",
|
||||
f"是否有真值: {'是' if result.get('truth_available') else '否'}",
|
||||
]
|
||||
if result.get("truth_available"):
|
||||
lines.append(f"真值RMS(y): {result.get('true_y_rms')}")
|
||||
lines.append(f"误差(%): {result.get('relative_error_percent')}")
|
||||
return "\n".join(lines)
|
||||
if task_label == TASK2_LABEL:
|
||||
lines = [
|
||||
f"自振频率f_n(Hz): {result.get('natural_frequency_hz')}",
|
||||
f"局部阻尼比zeta: {result.get('damping_ratio')}",
|
||||
f"全局阻尼比zeta: {result.get('global_damping_ratio')}",
|
||||
f"最高峰后衰减阻尼比zeta: {result.get('max_tail_damping_ratio')}",
|
||||
f"局部精确阻尼比: {result.get('damping_ratio_exact')}",
|
||||
f"全局精确阻尼比: {result.get('global_damping_ratio_exact')}",
|
||||
f"最高峰后衰减精确阻尼比: {result.get('max_tail_damping_ratio_exact')}",
|
||||
f"局部包络R²: {result.get('envelope_r2')}",
|
||||
f"全局包络R²: {result.get('global_envelope_r2')}",
|
||||
f"最高峰后衰减包络R²: {result.get('max_tail_envelope_r2')}",
|
||||
f"采样率(Hz): {result.get('sampling_rate_hz')}",
|
||||
]
|
||||
return "\n".join(lines)
|
||||
if task_label == TASK3_LABEL:
|
||||
lines = [
|
||||
f"预测激振频率(Hz): {result.get('predicted_frequency_hz')}",
|
||||
f"稳态采样点数: {result.get('segment_sample_count')}",
|
||||
f"采样率(Hz): {result.get('sampling_rate_hz')}",
|
||||
]
|
||||
if result.get("top_response_rms") is not None:
|
||||
lines.append(f"顶部RMS(y): {result.get('top_response_rms')}")
|
||||
if result.get("estimated_base_rms") is not None:
|
||||
lines.append(f"反演底部RMS(x): {result.get('estimated_base_rms')}")
|
||||
lines.append(f"反演TR: {result.get('estimated_tr_from_model')}")
|
||||
lines.append(f"反演残差(%): {result.get('inversion_residual_percent')}")
|
||||
lines.append(f"反演收敛: {result.get('inversion_converged')}")
|
||||
if result.get("true_frequency_hz") is not None:
|
||||
lines.append(f"文件名真值(Hz): {result.get('true_frequency_hz')}")
|
||||
lines.append(f"绝对误差(Hz): {result.get('absolute_error_hz')}")
|
||||
return "\n".join(lines)
|
||||
if task_label == TASK4_LABEL:
|
||||
lines = [
|
||||
f"挑战任务预测RMS: {result.get('pred_tmd_y_rms')}",
|
||||
f"基础模型RMS: {result.get('pred_base_y_rms')}",
|
||||
f"TR来源: {result.get('pred_tr_source')}",
|
||||
f"适配系数k(f): {result.get('adapter_scale_k')}",
|
||||
f"激振频率(Hz): {result.get('frequency_hz')}",
|
||||
f"是否有真值: {'是' if result.get('truth_available') else '否'}",
|
||||
]
|
||||
if result.get("truth_available"):
|
||||
lines.append(f"真值RMS(y): {result.get('true_y_rms')}")
|
||||
lines.append(f"误差(%): {result.get('relative_error_percent')}")
|
||||
return "\n".join(lines)
|
||||
return pretty_json_like_lines(result)
|
||||
|
||||
|
||||
class EvaluationStudio:
|
||||
def __init__(self) -> None:
|
||||
self.root = TkRootClass()
|
||||
self.root.title("建筑评测一体式测试系统")
|
||||
self.root.geometry("1180x760")
|
||||
self.root.minsize(980, 620)
|
||||
|
||||
self.task_var = tk.StringVar(value=TASK1_LABEL)
|
||||
self.file_var = tk.StringVar(value="")
|
||||
self.force_convert_var = tk.BooleanVar(value=False)
|
||||
self.timezone_var = tk.StringVar(value="8")
|
||||
self.status_var = tk.StringVar(value="就绪")
|
||||
self.preview_image: tk.PhotoImage | None = None
|
||||
self.original_image: tk.PhotoImage | None = None
|
||||
self.current_image_path: Path | None = None
|
||||
|
||||
self._build_ui()
|
||||
|
||||
def _build_ui(self) -> None:
|
||||
main = ttk.Frame(self.root, padding=10)
|
||||
main.pack(fill="both", expand=True)
|
||||
main.columnconfigure(1, weight=1)
|
||||
main.rowconfigure(4, weight=1)
|
||||
|
||||
ttk.Label(main, text="评测任务").grid(row=0, column=0, sticky="w", padx=(0, 8), pady=5)
|
||||
task_box = ttk.Combobox(main, textvariable=self.task_var, values=TASK_CHOICES, state="readonly")
|
||||
task_box.grid(row=0, column=1, sticky="ew", pady=5)
|
||||
|
||||
ttk.Label(main, text="CSV 文件").grid(row=1, column=0, sticky="w", padx=(0, 8), pady=5)
|
||||
file_entry = ttk.Entry(main, textvariable=self.file_var)
|
||||
file_entry.grid(row=1, column=1, sticky="ew", pady=5)
|
||||
ttk.Button(main, text="选择文件", command=self._choose_file).grid(row=1, column=2, padx=(8, 0), pady=5)
|
||||
|
||||
self.drop_label = ttk.Label(
|
||||
main,
|
||||
text="将 CSV 拖到这里(若无拖拽支持,可用“选择文件”)",
|
||||
relief="groove",
|
||||
anchor="center",
|
||||
)
|
||||
self.drop_label.grid(row=2, column=0, columnspan=3, sticky="ew", pady=(4, 8), ipady=12)
|
||||
if HAS_DND:
|
||||
self.drop_label.drop_target_register(DND_FILES)
|
||||
self.drop_label.dnd_bind("<<Drop>>", self._on_drop)
|
||||
|
||||
config_line = ttk.Frame(main)
|
||||
config_line.grid(row=3, column=0, columnspan=3, sticky="ew", pady=5)
|
||||
config_line.columnconfigure(6, weight=1)
|
||||
ttk.Checkbutton(
|
||||
config_line,
|
||||
text="强制按老师宽表格式转换",
|
||||
variable=self.force_convert_var,
|
||||
).grid(row=0, column=0, sticky="w")
|
||||
ttk.Label(config_line, text="字符串时间时区偏移(小时)").grid(row=0, column=1, sticky="w", padx=(12, 4))
|
||||
ttk.Entry(config_line, textvariable=self.timezone_var, width=7).grid(row=0, column=2, sticky="w")
|
||||
ttk.Button(config_line, text="开始评测", command=self._run_clicked).grid(row=0, column=3, padx=(12, 0))
|
||||
ttk.Label(config_line, textvariable=self.status_var).grid(row=0, column=4, sticky="w", padx=(12, 0))
|
||||
if not HAS_DND:
|
||||
ttk.Label(
|
||||
config_line,
|
||||
text="提示:未检测到 tkinterdnd2,拖拽不可用;可先 pip install tkinterdnd2",
|
||||
foreground="#aa5500",
|
||||
).grid(row=0, column=5, sticky="w", padx=(12, 0))
|
||||
|
||||
paned = ttk.Panedwindow(main, orient="horizontal")
|
||||
paned.grid(row=4, column=0, columnspan=3, sticky="nsew", pady=(6, 0))
|
||||
|
||||
left = ttk.Frame(paned, padding=8)
|
||||
right = ttk.Frame(paned, padding=8)
|
||||
paned.add(left, weight=3)
|
||||
paned.add(right, weight=2)
|
||||
|
||||
ttk.Label(left, text="关键结果").pack(anchor="w")
|
||||
self.key_text = ScrolledText(left, wrap="word", height=8)
|
||||
self.key_text.pack(fill="x", expand=False, pady=(0, 6))
|
||||
|
||||
ttk.Label(left, text="详细结果").pack(anchor="w")
|
||||
self.result_text = ScrolledText(left, wrap="word", height=16)
|
||||
self.result_text.pack(fill="both", expand=True)
|
||||
|
||||
ttk.Label(right, text="图像预览(最新结果)").pack(anchor="w")
|
||||
self.image_label = ttk.Label(right, relief="sunken", anchor="center")
|
||||
self.image_label.pack(fill="both", expand=True)
|
||||
self.image_label.bind("<Configure>", self._on_image_area_resize)
|
||||
|
||||
def _choose_file(self) -> None:
|
||||
file_path = filedialog.askopenfilename(
|
||||
title="选择CSV文件",
|
||||
filetypes=[("CSV Files", "*.csv"), ("All Files", "*.*")],
|
||||
initialdir=str(PROJECT_ROOT),
|
||||
)
|
||||
if file_path:
|
||||
self.file_var.set(file_path)
|
||||
|
||||
def _on_drop(self, event) -> None: # noqa: ANN001
|
||||
raw = event.data
|
||||
if not raw:
|
||||
return
|
||||
paths = self.root.tk.splitlist(raw)
|
||||
if not paths:
|
||||
return
|
||||
self.file_var.set(str(paths[0]).strip("{}"))
|
||||
|
||||
def _append_result(self, text: str) -> None:
|
||||
self.result_text.insert("end", text + "\n")
|
||||
self.result_text.see("end")
|
||||
|
||||
def _on_image_area_resize(self, _event) -> None: # noqa: ANN001
|
||||
self._refresh_image_fit()
|
||||
|
||||
def _refresh_image_fit(self) -> None:
|
||||
if self.original_image is None:
|
||||
return
|
||||
width = max(self.image_label.winfo_width(), 1)
|
||||
height = max(self.image_label.winfo_height(), 1)
|
||||
orig_w = max(self.original_image.width(), 1)
|
||||
orig_h = max(self.original_image.height(), 1)
|
||||
ratio_w = math.ceil(orig_w / width)
|
||||
ratio_h = math.ceil(orig_h / height)
|
||||
ratio = max(1, ratio_w, ratio_h)
|
||||
scaled = self.original_image.subsample(ratio, ratio)
|
||||
self.preview_image = scaled
|
||||
self.image_label.configure(image=scaled, text="")
|
||||
|
||||
def _set_image(self, image_path: Path | None) -> None:
|
||||
if image_path is None or not image_path.exists():
|
||||
self.image_label.configure(image="", text="暂无图像")
|
||||
self.preview_image = None
|
||||
self.original_image = None
|
||||
self.current_image_path = None
|
||||
return
|
||||
try:
|
||||
self.original_image = tk.PhotoImage(file=str(image_path))
|
||||
self.current_image_path = image_path
|
||||
self._refresh_image_fit()
|
||||
except Exception:
|
||||
self.image_label.configure(text=f"图像已生成:\n{image_path}", image="")
|
||||
self.preview_image = None
|
||||
self.original_image = None
|
||||
self.current_image_path = image_path
|
||||
|
||||
def _run_clicked(self) -> None:
|
||||
csv_text = self.file_var.get().strip()
|
||||
if not csv_text:
|
||||
messagebox.showerror("缺少输入", "请先选择或拖入CSV文件。")
|
||||
return
|
||||
csv_path = Path(csv_text).expanduser().resolve()
|
||||
if not csv_path.exists():
|
||||
messagebox.showerror("文件不存在", f"找不到文件:\n{csv_path}")
|
||||
return
|
||||
|
||||
try:
|
||||
tz_offset = float(self.timezone_var.get().strip())
|
||||
except ValueError:
|
||||
messagebox.showerror("参数错误", "时区偏移必须是数字。")
|
||||
return
|
||||
|
||||
self.status_var.set("运行中...")
|
||||
self.key_text.delete("1.0", "end")
|
||||
self.result_text.delete("1.0", "end")
|
||||
self._set_image(None)
|
||||
|
||||
thread = threading.Thread(
|
||||
target=self._run_pipeline,
|
||||
args=(csv_path, self.task_var.get(), self.force_convert_var.get(), tz_offset),
|
||||
daemon=True,
|
||||
)
|
||||
thread.start()
|
||||
|
||||
def _run_pipeline(self, csv_path: Path, task_label: str, force_convert: bool, tz_offset: float) -> None:
|
||||
try:
|
||||
APP_OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
|
||||
|
||||
use_csv, converted = ensure_long_format_csv(csv_path, force_convert=force_convert, tz_offset_hours=tz_offset)
|
||||
notes = [
|
||||
f"输入文件: {csv_path}",
|
||||
f"实际评测文件: {use_csv}",
|
||||
f"是否执行宽转长: {'是' if converted else '否'}",
|
||||
f"任务: {task_label}",
|
||||
"",
|
||||
]
|
||||
|
||||
if task_label == TASK1_LABEL:
|
||||
result, fig_path = run_task1(use_csv)
|
||||
elif task_label == TASK2_LABEL:
|
||||
result, fig_path = run_task2(use_csv)
|
||||
elif task_label == TASK3_LABEL:
|
||||
result, fig_path = run_task3(use_csv)
|
||||
elif task_label == TASK4_LABEL:
|
||||
result, fig_path = run_task4_tmd(use_csv)
|
||||
else:
|
||||
raise ValueError(f"不支持的任务:{task_label}")
|
||||
|
||||
output_text = "\n".join(notes) + pretty_json_like_lines(result)
|
||||
key_text = build_key_metrics(task_label, result)
|
||||
self.root.after(0, lambda: self.key_text.insert("end", key_text + "\n"))
|
||||
self.root.after(0, lambda: self.result_text.insert("end", output_text + "\n"))
|
||||
self.root.after(0, lambda: self._set_image(fig_path))
|
||||
self.root.after(0, lambda: self.status_var.set("完成"))
|
||||
except Exception as error:
|
||||
stack = traceback.format_exc()
|
||||
message = f"运行失败:{error}\n\n{stack}"
|
||||
self.root.after(0, lambda: self.key_text.insert("end", "运行失败,请查看详细错误。\n"))
|
||||
self.root.after(0, lambda: self.result_text.insert("end", message + "\n"))
|
||||
self.root.after(0, lambda: self.status_var.set("失败"))
|
||||
|
||||
def run(self) -> None:
|
||||
self.root.mainloop()
|
||||
|
||||
|
||||
def parse_args() -> argparse.Namespace:
|
||||
parser = argparse.ArgumentParser(description="Integrated visual evaluation studio for tasks 1~4.")
|
||||
parser.add_argument("--task", type=str, default=None, choices=TASK_CHOICES)
|
||||
parser.add_argument("--file", type=str, default=None, help="Optional CSV file path to preload.")
|
||||
parser.add_argument(
|
||||
"--force-convert",
|
||||
action="store_true",
|
||||
help="Force teacher wide-format conversion before running tasks.",
|
||||
)
|
||||
parser.add_argument("--timezone-offset-hours", type=float, default=8.0)
|
||||
return parser.parse_args()
|
||||
|
||||
|
||||
def main() -> None:
|
||||
args = parse_args()
|
||||
app = EvaluationStudio()
|
||||
if args.task:
|
||||
app.task_var.set(args.task)
|
||||
if args.file:
|
||||
app.file_var.set(str(Path(args.file).expanduser().resolve()))
|
||||
if args.force_convert:
|
||||
app.force_convert_var.set(True)
|
||||
app.timezone_var.set(str(args.timezone_offset_hours))
|
||||
app.run()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
6
requirement.txt
Normal file
6
requirement.txt
Normal file
@@ -0,0 +1,6 @@
|
||||
# Python 3.10
|
||||
numpy>=1.24,<2.0
|
||||
pandas>=2.0,<3.0
|
||||
scipy>=1.10,<2.0
|
||||
matplotlib>=3.7,<4.0
|
||||
scikit-learn>=1.3,<2.0
|
||||
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
|
||||
292
scripts/dataset.py
Normal file
292
scripts/dataset.py
Normal file
@@ -0,0 +1,292 @@
|
||||
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 files
|
||||
|
||||
|
||||
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 = 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()
|
||||
72
scripts/inference_utils.py
Normal file
72
scripts/inference_utils.py
Normal file
@@ -0,0 +1,72 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from functools import lru_cache
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
|
||||
def _task1_curve_candidates(project_root: Path) -> tuple[Path, ...]:
|
||||
base_dir = project_root / "evaluation_outputs" / "task1_final"
|
||||
return (
|
||||
base_dir / "task1_final_dense_curve.csv",
|
||||
base_dir / "evaluation_dense_curve.csv",
|
||||
)
|
||||
|
||||
|
||||
@lru_cache(maxsize=8)
|
||||
def load_task1_dense_curve(project_root_str: str) -> tuple[np.ndarray, np.ndarray, np.ndarray] | None:
|
||||
project_root = Path(project_root_str)
|
||||
for path in _task1_curve_candidates(project_root):
|
||||
if not path.exists():
|
||||
continue
|
||||
try:
|
||||
curve_df = pd.read_csv(path)
|
||||
except Exception:
|
||||
continue
|
||||
required = {"frequency_hz", "pred_tr", "x_rms_interp"}
|
||||
if not required.issubset(set(curve_df.columns)):
|
||||
continue
|
||||
cleaned = curve_df[["frequency_hz", "pred_tr", "x_rms_interp"]].copy()
|
||||
cleaned["frequency_hz"] = pd.to_numeric(cleaned["frequency_hz"], errors="coerce")
|
||||
cleaned["pred_tr"] = pd.to_numeric(cleaned["pred_tr"], errors="coerce")
|
||||
cleaned["x_rms_interp"] = pd.to_numeric(cleaned["x_rms_interp"], errors="coerce")
|
||||
cleaned = cleaned.dropna()
|
||||
if cleaned.empty:
|
||||
continue
|
||||
cleaned = cleaned.sort_values("frequency_hz").drop_duplicates(subset="frequency_hz", keep="first")
|
||||
return (
|
||||
cleaned["frequency_hz"].to_numpy(dtype=np.float64),
|
||||
cleaned["pred_tr"].to_numpy(dtype=np.float64),
|
||||
cleaned["x_rms_interp"].to_numpy(dtype=np.float64),
|
||||
)
|
||||
return None
|
||||
|
||||
|
||||
def predict_task1_tr(
|
||||
model,
|
||||
*,
|
||||
feature_vector: np.ndarray,
|
||||
frequency_hz: float,
|
||||
normalization_eps: float,
|
||||
project_root: Path,
|
||||
prefer_curve: bool = True,
|
||||
) -> tuple[float, str]:
|
||||
features = np.asarray(feature_vector, dtype=np.float64).reshape(1, -1)
|
||||
tr_value = max(float(model.predict(features)[0]), normalization_eps)
|
||||
source = "model_direct"
|
||||
|
||||
if prefer_curve:
|
||||
curve = load_task1_dense_curve(str(project_root.resolve()))
|
||||
if curve is not None and features.shape[1] >= 3:
|
||||
freq_grid, _tr_grid, x_ref_grid = curve
|
||||
query_freq = float(np.clip(frequency_hz, freq_grid.min(), freq_grid.max()))
|
||||
x_ref = float(np.interp(query_freq, freq_grid, x_ref_grid))
|
||||
x_obs = float(features[0, 2])
|
||||
ratio = x_obs / max(x_ref, normalization_eps)
|
||||
if ratio > 1.15 or ratio < 0.85:
|
||||
tr_value = max(tr_value / ratio, normalization_eps)
|
||||
source = "model_amp_adapt"
|
||||
|
||||
return tr_value, source
|
||||
119
scripts/predict_single.py
Normal file
119
scripts/predict_single.py
Normal file
@@ -0,0 +1,119 @@
|
||||
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
|
||||
from .inference_utils import predict_task1_tr
|
||||
except ImportError:
|
||||
from config import CORE_FEATURE_NAMES, evaluation_dir, make_experiment_config
|
||||
from dataset import build_record_from_file
|
||||
from inference_utils import predict_task1_tr
|
||||
|
||||
|
||||
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, tr_source = predict_task1_tr(
|
||||
model,
|
||||
feature_vector=record.features,
|
||||
frequency_hz=float(record.frequency_hz),
|
||||
normalization_eps=float(config.data.normalization_eps),
|
||||
project_root=config.data.project_root,
|
||||
prefer_curve=True,
|
||||
)
|
||||
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"TR Source: {tr_source}")
|
||||
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()
|
||||
43
scripts_2/README.md
Normal file
43
scripts_2/README.md
Normal file
@@ -0,0 +1,43 @@
|
||||
# 结构动力特性识别脚本 (Free-Vibration Identification)
|
||||
|
||||
## 项目简介
|
||||
本目录用于完成建筑结构动力特性识别任务。脚本基于专门的自由衰减数据文件,提取顶部响应传感器 `WSMS00007` 的 Z 向信号,识别结构一阶自振频率与阻尼比,并输出可视化诊断结果。
|
||||
|
||||
## 技术路线
|
||||
本目录的分析流程围绕自由衰减信号展开,包含信号读取、主频识别、滤波降噪与阻尼估计四个步骤:
|
||||
|
||||
1. **顶部响应读取**
|
||||
从 `free_vib` 目录下的 CSV 文件中读取 `WSMS00007` 的 `value3` 列,作为结构顶部响应信号。
|
||||
2. **频域主频识别**
|
||||
对原始自由衰减信号执行 FFT,在 `0.1Hz ~ 5.0Hz` 范围内寻找主峰,识别一阶自振频率。
|
||||
3. **带通滤波**
|
||||
以识别出的主频为中心,构造巴特沃斯带通滤波器,对原始信号进行滤波,得到平滑的衰减振动波形。
|
||||
4. **阻尼比计算**
|
||||
从滤波信号中提取波峰,自动选择包络线质量较高的一段,通过对数衰减率计算平均衰减量与阻尼比。
|
||||
|
||||
## 核心模块说明
|
||||
|
||||
- **`signal_utils.py`**
|
||||
共享信号处理工具。负责传感器常量定义、采样率估计、稳态段截取、FFT 主频识别、抛物线插值、带通滤波与波峰筛选。
|
||||
- **`task3_identify.py`**
|
||||
任务主入口。负责校验自由衰减数据来源、执行一阶频率与阻尼比识别,并生成包含原始波形、滤波波形和包络线的诊断图。
|
||||
|
||||
## 使用指南
|
||||
|
||||
**1. 运行任务脚本**
|
||||
对单个自由衰减 CSV 文件执行结构动力特性识别:
|
||||
```bash
|
||||
python scripts_2\task3_identify.py --file downloads\Non_TMD\free_vib\free_decay.csv
|
||||
```
|
||||
|
||||
**2. 指定输出目录**
|
||||
将图像与诊断结果输出到自定义目录:
|
||||
```bash
|
||||
python scripts_2\task3_identify.py --file downloads\Non_TMD\free_vib\free_decay.csv --output-dir evaluation_outputs\task3_custom
|
||||
```
|
||||
|
||||
**3. 调整滤波带宽**
|
||||
通过命令行参数修改主频两侧的带宽范围:
|
||||
```bash
|
||||
python scripts_2\task3_identify.py --file downloads\Non_TMD\free_vib\free_decay.csv --bandwidth 0.25
|
||||
```
|
||||
1
scripts_2/__init__.py
Normal file
1
scripts_2/__init__.py
Normal file
@@ -0,0 +1 @@
|
||||
"""Signal-processing scripts for free-vibration identification."""
|
||||
BIN
scripts_2/__pycache__/__init__.cpython-310.pyc
Normal file
BIN
scripts_2/__pycache__/__init__.cpython-310.pyc
Normal file
Binary file not shown.
BIN
scripts_2/__pycache__/__init__.cpython-314.pyc
Normal file
BIN
scripts_2/__pycache__/__init__.cpython-314.pyc
Normal file
Binary file not shown.
BIN
scripts_2/__pycache__/signal_utils.cpython-310.pyc
Normal file
BIN
scripts_2/__pycache__/signal_utils.cpython-310.pyc
Normal file
Binary file not shown.
BIN
scripts_2/__pycache__/signal_utils.cpython-314.pyc
Normal file
BIN
scripts_2/__pycache__/signal_utils.cpython-314.pyc
Normal file
Binary file not shown.
BIN
scripts_2/__pycache__/task3_identify.cpython-310.pyc
Normal file
BIN
scripts_2/__pycache__/task3_identify.cpython-310.pyc
Normal file
Binary file not shown.
BIN
scripts_2/__pycache__/task3_identify.cpython-314.pyc
Normal file
BIN
scripts_2/__pycache__/task3_identify.cpython-314.pyc
Normal file
Binary file not shown.
273
scripts_2/signal_utils.py
Normal file
273
scripts_2/signal_utils.py
Normal file
@@ -0,0 +1,273 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from pathlib import Path
|
||||
import re
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
from scipy.signal import butter, filtfilt, find_peaks
|
||||
|
||||
|
||||
TOP_RESPONSE_SENSOR = "WSMS00007"
|
||||
TOP_RESPONSE_AXIS = "value3"
|
||||
BASE_REFERENCE_SENSOR = "WSMS00012"
|
||||
BASE_REFERENCE_AXIS = "value1"
|
||||
TIME_COLUMN = "time"
|
||||
CODE_COLUMN = "code"
|
||||
|
||||
|
||||
def ensure_parent_dir(path: Path) -> Path:
|
||||
path = Path(path).resolve()
|
||||
path.mkdir(parents=True, exist_ok=True)
|
||||
return path
|
||||
|
||||
|
||||
def parse_frequency_from_filename(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:
|
||||
time_values = np.asarray(time_values, dtype=np.float64).reshape(-1)
|
||||
if time_values.size < 2:
|
||||
raise ValueError("At least two time samples are required to estimate sampling rate.")
|
||||
dt = np.diff(time_values)
|
||||
dt = dt[np.isfinite(dt)]
|
||||
dt = dt[dt > 0.0]
|
||||
if dt.size == 0:
|
||||
raise ValueError("Failed to infer a positive sampling interval from the CSV time column.")
|
||||
return float(1.0 / np.median(dt))
|
||||
|
||||
|
||||
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 read_sensor_signal(
|
||||
file_path: Path,
|
||||
sensor_code: str = TOP_RESPONSE_SENSOR,
|
||||
value_column: str = TOP_RESPONSE_AXIS,
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
df = pd.read_csv(file_path)
|
||||
sensor_df = df.loc[df[CODE_COLUMN] == sensor_code, [TIME_COLUMN, value_column]].copy()
|
||||
if sensor_df.empty:
|
||||
raise ValueError(f"Sensor {sensor_code}/{value_column} was not found in {file_path.name}.")
|
||||
|
||||
sensor_df = sensor_df.sort_values(TIME_COLUMN)
|
||||
sensor_df = sensor_df.drop_duplicates(subset=TIME_COLUMN, keep="first")
|
||||
sensor_df[TIME_COLUMN] = sensor_df[TIME_COLUMN].astype("float64")
|
||||
sensor_df[value_column] = sensor_df[value_column].astype("float64")
|
||||
|
||||
time_values = sensor_df[TIME_COLUMN].to_numpy(dtype=np.float64)
|
||||
signal = sensor_df[value_column].to_numpy(dtype=np.float64)
|
||||
if time_values.size < 8:
|
||||
raise ValueError(f"Too few usable samples in {file_path.name}.")
|
||||
return time_values, signal
|
||||
|
||||
|
||||
def get_middle_segment(
|
||||
time_values: np.ndarray,
|
||||
signal: np.ndarray,
|
||||
start_ratio: float = 0.20,
|
||||
end_ratio: float = 0.80,
|
||||
steady_window_ratio: float = 0.25,
|
||||
stride_ratio: float = 0.05,
|
||||
min_segment_length: int = 512,
|
||||
subwindow_count: int = 4,
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
signal = np.asarray(signal, dtype=np.float64).reshape(-1)
|
||||
time_values = np.asarray(time_values, dtype=np.float64).reshape(-1)
|
||||
if signal.size != time_values.size:
|
||||
raise ValueError("Time and signal arrays must have the same length.")
|
||||
|
||||
length = signal.size
|
||||
search_start = int(length * start_ratio)
|
||||
search_end = int(length * 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 < min_segment_length:
|
||||
center = length // 2
|
||||
half = min_segment_length // 2
|
||||
search_start = max(0, center - half)
|
||||
search_end = min(length, search_start + min_segment_length)
|
||||
search_start = max(0, search_end - min_segment_length)
|
||||
|
||||
search_length = search_end - search_start
|
||||
window_length = max(min_segment_length, int(length * steady_window_ratio))
|
||||
window_length = min(window_length, search_length)
|
||||
if search_length <= window_length:
|
||||
return time_values[search_start:search_end], signal[search_start:search_end]
|
||||
|
||||
stride = max(1, int(length * stride_ratio))
|
||||
candidate_signal = signal[search_start:search_end]
|
||||
candidate_rms = calculate_rms(candidate_signal)
|
||||
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_signal = signal[window_start:window_end]
|
||||
split_windows = np.array_split(window_signal, 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(), 1e-8))
|
||||
peak_cv = float(split_peaks.std() / max(split_peaks.mean(), 1e-8))
|
||||
window_rms = calculate_rms(window_signal)
|
||||
score = rms_cv + 0.35 * peak_cv - 0.05 * (window_rms / max(candidate_rms, 1e-8))
|
||||
if score < best_score:
|
||||
best_score = score
|
||||
best_slice = slice(window_start, window_end)
|
||||
|
||||
return time_values[best_slice], signal[best_slice]
|
||||
|
||||
|
||||
def compute_windowed_spectrum(
|
||||
signal: np.ndarray,
|
||||
sampling_rate: float,
|
||||
zero_padding_factor: int = 8,
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
signal = np.asarray(signal, dtype=np.float64).reshape(-1)
|
||||
if signal.size < 4:
|
||||
raise ValueError("Signal is too short for FFT analysis.")
|
||||
|
||||
centered = signal - np.mean(signal)
|
||||
window = np.hanning(centered.size)
|
||||
fft_size = int(2 ** np.ceil(np.log2(max(centered.size * zero_padding_factor, centered.size))))
|
||||
fft_values = np.fft.rfft(centered * window, n=fft_size)
|
||||
freqs = np.fft.rfftfreq(fft_size, d=1.0 / sampling_rate)
|
||||
magnitudes = 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 = float(magnitudes[peak_index - 1])
|
||||
beta = float(magnitudes[peak_index])
|
||||
gamma = float(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 dominant_frequency_in_band(
|
||||
signal: np.ndarray,
|
||||
sampling_rate: float,
|
||||
min_hz: float = 0.1,
|
||||
max_hz: float = 5.0,
|
||||
zero_padding_factor: int = 8,
|
||||
) -> tuple[float, np.ndarray, np.ndarray, int]:
|
||||
freqs, magnitudes = compute_windowed_spectrum(
|
||||
signal=signal,
|
||||
sampling_rate=sampling_rate,
|
||||
zero_padding_factor=zero_padding_factor,
|
||||
)
|
||||
band_mask = (freqs >= min_hz) & (freqs <= max_hz)
|
||||
if not np.any(band_mask):
|
||||
raise ValueError(f"No FFT bins fall inside {min_hz:.2f} Hz to {max_hz:.2f} Hz.")
|
||||
|
||||
band_indices = np.flatnonzero(band_mask)
|
||||
peak_index = int(band_indices[np.argmax(magnitudes[band_mask])])
|
||||
peak_frequency = parabolic_peak_frequency(freqs, magnitudes, peak_index)
|
||||
return peak_frequency, freqs, magnitudes, peak_index
|
||||
|
||||
|
||||
def bandpass_filter(
|
||||
signal: np.ndarray,
|
||||
sampling_rate: float,
|
||||
center_hz: float,
|
||||
bandwidth_hz: float = 0.3,
|
||||
order: int = 4,
|
||||
) -> np.ndarray:
|
||||
nyquist = 0.5 * sampling_rate
|
||||
low_hz = max(0.05, center_hz - bandwidth_hz)
|
||||
high_hz = min(center_hz + bandwidth_hz, nyquist * 0.98)
|
||||
if not low_hz < high_hz:
|
||||
raise ValueError("Invalid band-pass range. Check the estimated natural frequency and sampling rate.")
|
||||
|
||||
b, a = butter(order, [low_hz / nyquist, high_hz / nyquist], btype="bandpass")
|
||||
return filtfilt(b, a, np.asarray(signal, dtype=np.float64))
|
||||
|
||||
|
||||
def detect_decay_peaks(filtered_signal: np.ndarray, sampling_rate: float, frequency_hz: float) -> np.ndarray:
|
||||
filtered_signal = np.asarray(filtered_signal, dtype=np.float64).reshape(-1)
|
||||
samples_per_cycle = max(int(round(sampling_rate / max(frequency_hz, 1e-6))), 3)
|
||||
min_peak_distance = max(1, int(0.65 * samples_per_cycle))
|
||||
peaks, _ = find_peaks(filtered_signal, distance=min_peak_distance)
|
||||
if peaks.size < 4:
|
||||
peaks, _ = find_peaks(np.abs(filtered_signal), distance=max(1, min_peak_distance // 2))
|
||||
return peaks.astype(int)
|
||||
|
||||
|
||||
def select_decay_peak_window(
|
||||
filtered_signal: np.ndarray,
|
||||
peak_indices: np.ndarray,
|
||||
sampling_rate: float,
|
||||
min_peaks: int = 5,
|
||||
max_peaks: int = 12,
|
||||
) -> dict[str, np.ndarray | float | int]:
|
||||
filtered_signal = np.asarray(filtered_signal, dtype=np.float64).reshape(-1)
|
||||
peak_indices = np.asarray(peak_indices, dtype=int).reshape(-1)
|
||||
if peak_indices.size < min_peaks:
|
||||
raise ValueError("Too few peaks were detected for log-decrement damping estimation.")
|
||||
|
||||
amplitudes = np.abs(filtered_signal[peak_indices])
|
||||
tail_start = int(filtered_signal.size * 0.85)
|
||||
noise_slice = filtered_signal[tail_start:] if tail_start < filtered_signal.size else filtered_signal
|
||||
noise_floor = max(float(np.median(np.abs(noise_slice))), 1e-8)
|
||||
useful_mask = amplitudes > max(3.0 * noise_floor, 0.05 * float(amplitudes.max()))
|
||||
useful_indices = peak_indices[useful_mask]
|
||||
useful_amplitudes = amplitudes[useful_mask]
|
||||
if useful_indices.size < min_peaks:
|
||||
useful_indices = peak_indices
|
||||
useful_amplitudes = amplitudes
|
||||
|
||||
best_window: dict[str, np.ndarray | float | int] | None = None
|
||||
best_score = -float("inf")
|
||||
max_window_size = min(max_peaks, useful_indices.size)
|
||||
|
||||
for window_size in range(min_peaks, max_window_size + 1):
|
||||
for start in range(1, useful_indices.size - window_size + 1):
|
||||
stop = start + window_size
|
||||
idx_window = useful_indices[start:stop]
|
||||
amp_window = useful_amplitudes[start:stop]
|
||||
if np.any(amp_window <= 0.0):
|
||||
continue
|
||||
|
||||
time_window = idx_window.astype(np.float64) / sampling_rate
|
||||
log_amp = np.log(amp_window)
|
||||
slope, intercept = np.polyfit(time_window, log_amp, 1)
|
||||
if slope >= 0.0:
|
||||
continue
|
||||
|
||||
fit_values = slope * time_window + intercept
|
||||
ss_res = float(np.sum((log_amp - fit_values) ** 2))
|
||||
ss_tot = float(np.sum((log_amp - np.mean(log_amp)) ** 2))
|
||||
r_squared = 1.0 - ss_res / max(ss_tot, 1e-12)
|
||||
amplitude_ratio = float(amp_window[0] / max(amp_window[-1], 1e-12))
|
||||
score = r_squared + 0.02 * np.log(max(amplitude_ratio, 1.0))
|
||||
if score > best_score:
|
||||
best_score = score
|
||||
best_window = {
|
||||
"peak_indices": idx_window,
|
||||
"peak_amplitudes": amp_window,
|
||||
"time_window": time_window,
|
||||
"fit_log_amplitudes": fit_values,
|
||||
"fit_amplitudes": np.exp(fit_values),
|
||||
"r_squared": r_squared,
|
||||
"noise_floor": noise_floor,
|
||||
}
|
||||
|
||||
if best_window is None:
|
||||
raise ValueError("Failed to isolate a clean decay envelope segment for damping estimation.")
|
||||
return best_window
|
||||
381
scripts_2/task3_identify.py
Normal file
381
scripts_2/task3_identify.py
Normal file
@@ -0,0 +1,381 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
from pathlib import Path
|
||||
import sys
|
||||
|
||||
import matplotlib
|
||||
|
||||
matplotlib.use("Agg")
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
if __package__ is None or __package__ == "":
|
||||
sys.path.append(str(Path(__file__).resolve().parents[1]))
|
||||
|
||||
from scripts_2.signal_utils import (
|
||||
TOP_RESPONSE_AXIS,
|
||||
TOP_RESPONSE_SENSOR,
|
||||
bandpass_filter,
|
||||
detect_decay_peaks,
|
||||
dominant_frequency_in_band,
|
||||
ensure_parent_dir,
|
||||
estimate_sampling_rate,
|
||||
read_sensor_signal,
|
||||
select_decay_peak_window,
|
||||
)
|
||||
|
||||
|
||||
def parse_args() -> argparse.Namespace:
|
||||
project_root = Path(__file__).resolve().parents[1]
|
||||
parser = argparse.ArgumentParser(description="Identify the first natural frequency and damping ratio from free-decay data.")
|
||||
parser.add_argument(
|
||||
"--file",
|
||||
type=str,
|
||||
default=str(project_root / "downloads" / "Non_TMD" / "free_vib" / "free_decay.csv"),
|
||||
help="Path to a free-vibration CSV inside the dedicated free_vib folder.",
|
||||
)
|
||||
parser.add_argument(
|
||||
"--output-dir",
|
||||
type=str,
|
||||
default=str(project_root / "evaluation_outputs" / "task3_identify"),
|
||||
help="Directory for plots and optional exported diagnostics.",
|
||||
)
|
||||
parser.add_argument(
|
||||
"--bandwidth",
|
||||
type=float,
|
||||
default=0.3,
|
||||
help="Half-bandwidth (Hz) of the Butterworth band-pass filter around the identified natural frequency.",
|
||||
)
|
||||
parser.add_argument(
|
||||
"--show",
|
||||
action="store_true",
|
||||
help="Display the generated figure in addition to saving it.",
|
||||
)
|
||||
return parser.parse_args()
|
||||
|
||||
|
||||
def validate_free_vibration_source(file_path: Path) -> None:
|
||||
if "free_vib" not in {part.lower() for part in file_path.parts}:
|
||||
raise ValueError(
|
||||
f"{file_path.name} is not inside a dedicated free_vib folder. Forced-vibration tails are forbidden for task 3."
|
||||
)
|
||||
|
||||
|
||||
def build_envelope_curve(
|
||||
time_values: np.ndarray,
|
||||
peak_times: np.ndarray,
|
||||
fit_amplitudes: np.ndarray,
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
if peak_times.size < 2:
|
||||
return time_values, np.zeros_like(time_values)
|
||||
envelope = np.interp(time_values, peak_times, fit_amplitudes, left=fit_amplitudes[0], right=fit_amplitudes[-1])
|
||||
return time_values, envelope
|
||||
|
||||
|
||||
def log_decrement_metrics(peak_amplitudes: np.ndarray) -> tuple[float, float, float]:
|
||||
amplitudes = np.asarray(peak_amplitudes, dtype=np.float64).reshape(-1)
|
||||
if amplitudes.size < 2 or np.any(amplitudes <= 0.0):
|
||||
raise ValueError("At least two positive peak amplitudes are required for damping estimation.")
|
||||
log_decrements = np.log(amplitudes[:-1] / amplitudes[1:])
|
||||
mean_log_decrement = float(np.mean(log_decrements))
|
||||
damping_ratio = float(mean_log_decrement / (2.0 * np.pi))
|
||||
damping_ratio_exact = float(
|
||||
mean_log_decrement / np.sqrt((2.0 * np.pi) ** 2 + mean_log_decrement**2)
|
||||
)
|
||||
return mean_log_decrement, damping_ratio, damping_ratio_exact
|
||||
|
||||
|
||||
def estimate_global_decay_metrics(
|
||||
filtered_signal: np.ndarray,
|
||||
peak_indices: np.ndarray,
|
||||
sampling_rate: float,
|
||||
min_peaks: int = 5,
|
||||
) -> dict[str, np.ndarray | float]:
|
||||
filtered_signal = np.asarray(filtered_signal, dtype=np.float64).reshape(-1)
|
||||
peak_indices = np.asarray(peak_indices, dtype=int).reshape(-1)
|
||||
if peak_indices.size < min_peaks:
|
||||
raise ValueError("Too few peaks were detected for global damping estimation.")
|
||||
|
||||
amplitudes = np.abs(filtered_signal[peak_indices])
|
||||
tail_start = int(filtered_signal.size * 0.85)
|
||||
noise_slice = filtered_signal[tail_start:] if tail_start < filtered_signal.size else filtered_signal
|
||||
noise_floor = max(float(np.median(np.abs(noise_slice))), 1e-8)
|
||||
useful_mask = amplitudes > max(3.0 * noise_floor, 0.05 * float(amplitudes.max()))
|
||||
useful_indices = peak_indices[useful_mask]
|
||||
useful_amplitudes = amplitudes[useful_mask]
|
||||
if useful_indices.size < min_peaks:
|
||||
useful_indices = peak_indices
|
||||
useful_amplitudes = amplitudes
|
||||
|
||||
time_window = useful_indices.astype(np.float64) / sampling_rate
|
||||
log_amp = np.log(np.maximum(useful_amplitudes, 1e-12))
|
||||
slope, intercept = np.polyfit(time_window, log_amp, 1)
|
||||
fit_log = slope * time_window + intercept
|
||||
ss_res = float(np.sum((log_amp - fit_log) ** 2))
|
||||
ss_tot = float(np.sum((log_amp - np.mean(log_amp)) ** 2))
|
||||
r_squared = 1.0 - ss_res / max(ss_tot, 1e-12)
|
||||
mean_log_decrement, damping_ratio, damping_ratio_exact = log_decrement_metrics(useful_amplitudes)
|
||||
|
||||
return {
|
||||
"peak_indices": useful_indices,
|
||||
"peak_amplitudes": useful_amplitudes,
|
||||
"fit_amplitudes": np.exp(fit_log),
|
||||
"r_squared": float(r_squared),
|
||||
"mean_log_decrement": float(mean_log_decrement),
|
||||
"damping_ratio": float(damping_ratio),
|
||||
"damping_ratio_exact": float(damping_ratio_exact),
|
||||
}
|
||||
|
||||
|
||||
def estimate_tail_from_max_peak_metrics(
|
||||
filtered_signal: np.ndarray,
|
||||
peak_indices: np.ndarray,
|
||||
sampling_rate: float,
|
||||
min_peaks: int = 5,
|
||||
) -> dict[str, np.ndarray | float]:
|
||||
filtered_signal = np.asarray(filtered_signal, dtype=np.float64).reshape(-1)
|
||||
peak_indices = np.asarray(peak_indices, dtype=int).reshape(-1)
|
||||
if peak_indices.size < min_peaks:
|
||||
raise ValueError("Too few peaks were detected for max-peak-tail damping estimation.")
|
||||
|
||||
amplitudes = np.abs(filtered_signal[peak_indices])
|
||||
max_peak_pos = int(np.argmax(amplitudes))
|
||||
tail_indices = peak_indices[max_peak_pos:]
|
||||
tail_amplitudes = amplitudes[max_peak_pos:]
|
||||
if tail_indices.size < min_peaks:
|
||||
raise ValueError("Too few peaks after the maximum peak for tail damping estimation.")
|
||||
|
||||
tail_start = int(filtered_signal.size * 0.85)
|
||||
noise_slice = filtered_signal[tail_start:] if tail_start < filtered_signal.size else filtered_signal
|
||||
noise_floor = max(float(np.median(np.abs(noise_slice))), 1e-8)
|
||||
useful_mask = tail_amplitudes > max(3.0 * noise_floor, 0.05 * float(tail_amplitudes.max()))
|
||||
useful_indices = tail_indices[useful_mask]
|
||||
useful_amplitudes = tail_amplitudes[useful_mask]
|
||||
if useful_indices.size < min_peaks:
|
||||
useful_indices = tail_indices
|
||||
useful_amplitudes = tail_amplitudes
|
||||
|
||||
time_window = useful_indices.astype(np.float64) / sampling_rate
|
||||
log_amp = np.log(np.maximum(useful_amplitudes, 1e-12))
|
||||
slope, intercept = np.polyfit(time_window, log_amp, 1)
|
||||
fit_log = slope * time_window + intercept
|
||||
ss_res = float(np.sum((log_amp - fit_log) ** 2))
|
||||
ss_tot = float(np.sum((log_amp - np.mean(log_amp)) ** 2))
|
||||
r_squared = 1.0 - ss_res / max(ss_tot, 1e-12)
|
||||
mean_log_decrement, damping_ratio, damping_ratio_exact = log_decrement_metrics(useful_amplitudes)
|
||||
|
||||
return {
|
||||
"peak_indices": useful_indices,
|
||||
"peak_amplitudes": useful_amplitudes,
|
||||
"fit_amplitudes": np.exp(fit_log),
|
||||
"r_squared": float(r_squared),
|
||||
"mean_log_decrement": float(mean_log_decrement),
|
||||
"damping_ratio": float(damping_ratio),
|
||||
"damping_ratio_exact": float(damping_ratio_exact),
|
||||
}
|
||||
|
||||
|
||||
def save_task3_figure(
|
||||
file_path: Path,
|
||||
output_dir: Path,
|
||||
time_values: np.ndarray,
|
||||
raw_signal: np.ndarray,
|
||||
filtered_signal: np.ndarray,
|
||||
peak_indices: np.ndarray,
|
||||
selected_peak_indices: np.ndarray,
|
||||
selected_peak_amplitudes: np.ndarray,
|
||||
fit_amplitudes: np.ndarray,
|
||||
natural_frequency_hz: float,
|
||||
local_damping_ratio: float,
|
||||
global_peak_indices: np.ndarray,
|
||||
global_peak_amplitudes: np.ndarray,
|
||||
global_fit_amplitudes: np.ndarray,
|
||||
global_damping_ratio: float,
|
||||
tail_peak_indices: np.ndarray,
|
||||
tail_peak_amplitudes: np.ndarray,
|
||||
tail_fit_amplitudes: np.ndarray,
|
||||
tail_damping_ratio: float,
|
||||
) -> Path:
|
||||
output_dir = ensure_parent_dir(output_dir)
|
||||
peak_times = time_values[selected_peak_indices]
|
||||
envelope_time, envelope = build_envelope_curve(time_values, peak_times, fit_amplitudes)
|
||||
|
||||
fig, axes = plt.subplots(3, 1, figsize=(12, 10), sharex=False)
|
||||
fig.suptitle(f"Task 3 | Free-Vibration Identification | {file_path.name}")
|
||||
|
||||
axes[0].plot(time_values, raw_signal, color="tab:gray", linewidth=1.0, label="Raw top response")
|
||||
axes[0].set_title(f"Top Response ({TOP_RESPONSE_SENSOR} / {TOP_RESPONSE_AXIS})")
|
||||
axes[0].set_ylabel("Acceleration")
|
||||
axes[0].grid(True, alpha=0.25)
|
||||
axes[0].legend()
|
||||
|
||||
axes[1].plot(time_values, filtered_signal, color="tab:blue", linewidth=1.0, label="Band-pass filtered")
|
||||
axes[1].scatter(time_values[peak_indices], filtered_signal[peak_indices], color="tab:orange", s=14, label="Detected peaks")
|
||||
axes[1].set_title(f"Filtered around f_n = {natural_frequency_hz:.4f} Hz")
|
||||
axes[1].set_ylabel("Acceleration")
|
||||
axes[1].grid(True, alpha=0.25)
|
||||
axes[1].legend()
|
||||
|
||||
axes[2].plot(time_values, filtered_signal, color="tab:blue", linewidth=1.0, alpha=0.65, label="Filtered signal")
|
||||
axes[2].scatter(peak_times, selected_peak_amplitudes, color="tab:red", s=28, label="Selected envelope peaks")
|
||||
axes[2].plot(envelope_time, envelope, color="tab:green", linewidth=1.8, label="Fitted positive envelope")
|
||||
axes[2].plot(envelope_time, -envelope, color="tab:green", linewidth=1.2, linestyle="--", label="Fitted negative envelope")
|
||||
global_peak_times = time_values[global_peak_indices]
|
||||
axes[2].scatter(global_peak_times, global_peak_amplitudes, color="tab:purple", s=14, alpha=0.65, label="Global peaks")
|
||||
axes[2].plot(global_peak_times, global_fit_amplitudes, color="tab:purple", linewidth=1.2, linestyle="-.", label="Global envelope fit")
|
||||
tail_peak_times = time_values[tail_peak_indices]
|
||||
axes[2].scatter(tail_peak_times, tail_peak_amplitudes, color="tab:brown", s=14, alpha=0.65, label="Max-peak tail peaks")
|
||||
axes[2].plot(tail_peak_times, tail_fit_amplitudes, color="tab:brown", linewidth=1.2, linestyle=":", label="Max-peak tail fit")
|
||||
axes[2].set_title(
|
||||
"Damping ratio | "
|
||||
f"local zeta={local_damping_ratio:.5f} | "
|
||||
f"global zeta={global_damping_ratio:.5f} | "
|
||||
f"max-tail zeta={tail_damping_ratio:.5f}"
|
||||
)
|
||||
axes[2].set_xlabel("Time (s)")
|
||||
axes[2].set_ylabel("Acceleration")
|
||||
axes[2].grid(True, alpha=0.25)
|
||||
axes[2].legend()
|
||||
|
||||
plt.tight_layout()
|
||||
figure_path = output_dir / f"{file_path.stem}_task3_identify.png"
|
||||
plt.savefig(figure_path, dpi=180, bbox_inches="tight")
|
||||
return figure_path
|
||||
|
||||
|
||||
def analyze_free_vibration(
|
||||
file_path: Path,
|
||||
output_dir: Path,
|
||||
bandwidth_hz: float = 0.3,
|
||||
show: bool = False,
|
||||
) -> dict[str, float | int | str | Path]:
|
||||
file_path = Path(file_path).resolve()
|
||||
validate_free_vibration_source(file_path)
|
||||
|
||||
time_values, raw_signal = read_sensor_signal(file_path, sensor_code=TOP_RESPONSE_SENSOR, value_column=TOP_RESPONSE_AXIS)
|
||||
sampling_rate = estimate_sampling_rate(time_values)
|
||||
natural_frequency_hz, freqs, magnitudes, peak_index = dominant_frequency_in_band(
|
||||
raw_signal,
|
||||
sampling_rate=sampling_rate,
|
||||
min_hz=0.1,
|
||||
max_hz=5.0,
|
||||
zero_padding_factor=16,
|
||||
)
|
||||
filtered_signal = bandpass_filter(
|
||||
raw_signal,
|
||||
sampling_rate=sampling_rate,
|
||||
center_hz=natural_frequency_hz,
|
||||
bandwidth_hz=bandwidth_hz,
|
||||
order=4,
|
||||
)
|
||||
peak_indices = detect_decay_peaks(filtered_signal, sampling_rate=sampling_rate, frequency_hz=natural_frequency_hz)
|
||||
peak_window = select_decay_peak_window(filtered_signal, peak_indices, sampling_rate=sampling_rate)
|
||||
|
||||
selected_peak_indices = np.asarray(peak_window["peak_indices"], dtype=int)
|
||||
selected_peak_amplitudes = np.asarray(peak_window["peak_amplitudes"], dtype=np.float64)
|
||||
fit_amplitudes = np.asarray(peak_window["fit_amplitudes"], dtype=np.float64)
|
||||
mean_log_decrement, damping_ratio, damping_ratio_exact = log_decrement_metrics(selected_peak_amplitudes)
|
||||
global_metrics = estimate_global_decay_metrics(
|
||||
filtered_signal=filtered_signal,
|
||||
peak_indices=peak_indices,
|
||||
sampling_rate=sampling_rate,
|
||||
)
|
||||
tail_metrics = estimate_tail_from_max_peak_metrics(
|
||||
filtered_signal=filtered_signal,
|
||||
peak_indices=peak_indices,
|
||||
sampling_rate=sampling_rate,
|
||||
)
|
||||
|
||||
figure_path = save_task3_figure(
|
||||
file_path=file_path,
|
||||
output_dir=output_dir,
|
||||
time_values=time_values,
|
||||
raw_signal=raw_signal,
|
||||
filtered_signal=filtered_signal,
|
||||
peak_indices=peak_indices,
|
||||
selected_peak_indices=selected_peak_indices,
|
||||
selected_peak_amplitudes=selected_peak_amplitudes,
|
||||
fit_amplitudes=fit_amplitudes,
|
||||
natural_frequency_hz=natural_frequency_hz,
|
||||
local_damping_ratio=damping_ratio,
|
||||
global_peak_indices=np.asarray(global_metrics["peak_indices"], dtype=int),
|
||||
global_peak_amplitudes=np.asarray(global_metrics["peak_amplitudes"], dtype=np.float64),
|
||||
global_fit_amplitudes=np.asarray(global_metrics["fit_amplitudes"], dtype=np.float64),
|
||||
global_damping_ratio=float(global_metrics["damping_ratio"]),
|
||||
tail_peak_indices=np.asarray(tail_metrics["peak_indices"], dtype=int),
|
||||
tail_peak_amplitudes=np.asarray(tail_metrics["peak_amplitudes"], dtype=np.float64),
|
||||
tail_fit_amplitudes=np.asarray(tail_metrics["fit_amplitudes"], dtype=np.float64),
|
||||
tail_damping_ratio=float(tail_metrics["damping_ratio"]),
|
||||
)
|
||||
if show:
|
||||
plt.show()
|
||||
plt.close("all")
|
||||
|
||||
return {
|
||||
"file_path": str(file_path),
|
||||
"sensor_code": TOP_RESPONSE_SENSOR,
|
||||
"axis": TOP_RESPONSE_AXIS,
|
||||
"sampling_rate_hz": float(sampling_rate),
|
||||
"natural_frequency_hz": float(natural_frequency_hz),
|
||||
"fft_peak_bin_hz": float(freqs[peak_index]),
|
||||
"fft_peak_magnitude": float(magnitudes[peak_index]),
|
||||
"detected_peak_count": int(len(peak_indices)),
|
||||
"selected_peak_count": int(len(selected_peak_indices)),
|
||||
"mean_log_decrement": mean_log_decrement,
|
||||
"damping_ratio": damping_ratio,
|
||||
"damping_ratio_exact": damping_ratio_exact,
|
||||
"global_selected_peak_count": int(len(np.asarray(global_metrics["peak_indices"], dtype=int))),
|
||||
"global_mean_log_decrement": float(global_metrics["mean_log_decrement"]),
|
||||
"global_damping_ratio": float(global_metrics["damping_ratio"]),
|
||||
"global_damping_ratio_exact": float(global_metrics["damping_ratio_exact"]),
|
||||
"global_envelope_r2": float(global_metrics["r_squared"]),
|
||||
"max_tail_selected_peak_count": int(len(np.asarray(tail_metrics["peak_indices"], dtype=int))),
|
||||
"max_tail_mean_log_decrement": float(tail_metrics["mean_log_decrement"]),
|
||||
"max_tail_damping_ratio": float(tail_metrics["damping_ratio"]),
|
||||
"max_tail_damping_ratio_exact": float(tail_metrics["damping_ratio_exact"]),
|
||||
"max_tail_envelope_r2": float(tail_metrics["r_squared"]),
|
||||
"envelope_r2": float(peak_window["r_squared"]),
|
||||
"figure_path": str(figure_path),
|
||||
}
|
||||
|
||||
|
||||
def print_report(result: dict[str, float | int | str | Path]) -> None:
|
||||
print("=== Task 3: Structural Dynamic Identification ===")
|
||||
print(f"Input file: {result['file_path']}")
|
||||
print(f"Top response sensor: {result['sensor_code']} / {result['axis']}")
|
||||
print(f"Sampling rate: {result['sampling_rate_hz']:.4f} Hz")
|
||||
print(f"Natural frequency f_n: {result['natural_frequency_hz']:.6f} Hz")
|
||||
print(f"Local mean log decrement delta: {result['mean_log_decrement']:.6f}")
|
||||
print(f"Local damping ratio zeta (delta / 2pi): {result['damping_ratio']:.6f}")
|
||||
print(f"Local damping ratio exact: {result['damping_ratio_exact']:.6f}")
|
||||
print(f"Global mean log decrement delta: {result['global_mean_log_decrement']:.6f}")
|
||||
print(f"Global damping ratio zeta (delta / 2pi): {result['global_damping_ratio']:.6f}")
|
||||
print(f"Global damping ratio exact: {result['global_damping_ratio_exact']:.6f}")
|
||||
print(f"Max-tail mean log decrement delta: {result['max_tail_mean_log_decrement']:.6f}")
|
||||
print(f"Max-tail damping ratio zeta (delta / 2pi): {result['max_tail_damping_ratio']:.6f}")
|
||||
print(f"Max-tail damping ratio exact: {result['max_tail_damping_ratio_exact']:.6f}")
|
||||
print(f"Detected peaks: {result['detected_peak_count']}")
|
||||
print(f"Selected peaks for local envelope: {result['selected_peak_count']}")
|
||||
print(f"Selected peaks for global envelope: {result['global_selected_peak_count']}")
|
||||
print(f"Selected peaks for max-tail envelope: {result['max_tail_selected_peak_count']}")
|
||||
print(f"Local envelope linearity R^2: {result['envelope_r2']:.6f}")
|
||||
print(f"Global envelope linearity R^2: {result['global_envelope_r2']:.6f}")
|
||||
print(f"Max-tail envelope linearity R^2: {result['max_tail_envelope_r2']:.6f}")
|
||||
print(f"Figure saved to: {result['figure_path']}")
|
||||
|
||||
|
||||
def main() -> None:
|
||||
args = parse_args()
|
||||
result = analyze_free_vibration(
|
||||
file_path=Path(args.file),
|
||||
output_dir=Path(args.output_dir),
|
||||
bandwidth_hz=float(args.bandwidth),
|
||||
show=bool(args.show),
|
||||
)
|
||||
print_report(result)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
33
scripts_3/README.md
Normal file
33
scripts_3/README.md
Normal file
@@ -0,0 +1,33 @@
|
||||
# 激振频率逆向预测 (Forced-Vibration Frequency Prediction)
|
||||
|
||||
## 项目简介
|
||||
本目录用于完成强迫振动频率逆向预测任务。脚本基于建筑顶部响应传感器 `WSMS00007` 的 Z 向数据,从稳态振动段中识别主频,反推出地震台输入的简谐波频率,并输出频率误差与可视化结果。
|
||||
|
||||
## 技术路线
|
||||
本目录的分析流程覆盖稳态段提取、主频精细识别与结果校验三个部分:
|
||||
|
||||
1. **稳态中段截取**
|
||||
对强迫振动时程执行滑动窗口搜索,从中间区域选出波动最稳定的响应片段,用于频谱分析。
|
||||
2. **顶部响应主频识别**
|
||||
读取 `WSMS00007` 的 `value3` 列,对稳态段进行加窗 FFT 与补零计算,并结合抛物线插值输出高精度激振频率。
|
||||
3. **文件名真值校验**
|
||||
从样本文件名中提取频率标签,计算预测值与真值之间的绝对误差。
|
||||
|
||||
## 核心模块说明
|
||||
|
||||
- **`task4_predict_freq.py`**
|
||||
任务主脚本。负责稳态段截取、频谱分析、抛物线插值、误差计算与频域图像生成。
|
||||
|
||||
## 使用指南
|
||||
|
||||
**1. 运行任务脚本**
|
||||
对单个强迫振动样本预测激振频率:
|
||||
```bash
|
||||
python scripts_3\task4_predict_freq.py --file downloads\Non_TMD\val\harmonic_5mm_0.75Hz.csv
|
||||
```
|
||||
|
||||
**2. 指定输出目录**
|
||||
将图像保存到自定义目录:
|
||||
```bash
|
||||
python scripts_3\task4_predict_freq.py --file downloads\Non_TMD\val\harmonic_5mm_0.75Hz.csv --output-dir evaluation_outputs\task4_custom
|
||||
```
|
||||
1
scripts_3/__init__.py
Normal file
1
scripts_3/__init__.py
Normal file
@@ -0,0 +1 @@
|
||||
"""Signal-processing scripts for forced-vibration frequency prediction."""
|
||||
BIN
scripts_3/__pycache__/__init__.cpython-310.pyc
Normal file
BIN
scripts_3/__pycache__/__init__.cpython-310.pyc
Normal file
Binary file not shown.
BIN
scripts_3/__pycache__/__init__.cpython-314.pyc
Normal file
BIN
scripts_3/__pycache__/__init__.cpython-314.pyc
Normal file
Binary file not shown.
BIN
scripts_3/__pycache__/task4_predict_freq.cpython-310.pyc
Normal file
BIN
scripts_3/__pycache__/task4_predict_freq.cpython-310.pyc
Normal file
Binary file not shown.
BIN
scripts_3/__pycache__/task4_predict_freq.cpython-314.pyc
Normal file
BIN
scripts_3/__pycache__/task4_predict_freq.cpython-314.pyc
Normal file
Binary file not shown.
221
scripts_3/task4_predict_freq.py
Normal file
221
scripts_3/task4_predict_freq.py
Normal file
@@ -0,0 +1,221 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
from pathlib import Path
|
||||
import sys
|
||||
|
||||
import matplotlib
|
||||
|
||||
matplotlib.use("Agg")
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
from scipy.signal import find_peaks
|
||||
|
||||
if __package__ is None or __package__ == "":
|
||||
sys.path.append(str(Path(__file__).resolve().parents[1]))
|
||||
|
||||
from scripts_2.signal_utils import (
|
||||
TOP_RESPONSE_AXIS,
|
||||
TOP_RESPONSE_SENSOR,
|
||||
dominant_frequency_in_band,
|
||||
ensure_parent_dir,
|
||||
estimate_sampling_rate,
|
||||
get_middle_segment,
|
||||
parse_frequency_from_filename,
|
||||
read_sensor_signal,
|
||||
)
|
||||
|
||||
|
||||
def parse_args() -> argparse.Namespace:
|
||||
project_root = Path(__file__).resolve().parents[1]
|
||||
parser = argparse.ArgumentParser(description="Predict the harmonic input frequency from top-response steady-state data only.")
|
||||
parser.add_argument(
|
||||
"--file",
|
||||
type=str,
|
||||
default=str(project_root / "downloads" / "Non_TMD" / "val" / "harmonic_5mm_0.75Hz.csv"),
|
||||
help="Path to a forced-vibration CSV such as data from Non_TMD val/test.",
|
||||
)
|
||||
parser.add_argument(
|
||||
"--output-dir",
|
||||
type=str,
|
||||
default=str(project_root / "evaluation_outputs" / "task4_predict_freq"),
|
||||
help="Directory for plots and optional exported diagnostics.",
|
||||
)
|
||||
parser.add_argument(
|
||||
"--show",
|
||||
action="store_true",
|
||||
help="Display the generated figure in addition to saving it.",
|
||||
)
|
||||
return parser.parse_args()
|
||||
|
||||
|
||||
def save_task4_figure(
|
||||
file_path: Path,
|
||||
output_dir: Path,
|
||||
time_segment: np.ndarray,
|
||||
signal_segment: np.ndarray,
|
||||
freqs: np.ndarray,
|
||||
magnitudes: np.ndarray,
|
||||
predicted_frequency_hz: float,
|
||||
true_frequency_hz: float | None,
|
||||
) -> Path:
|
||||
output_dir = ensure_parent_dir(output_dir)
|
||||
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
|
||||
fig.suptitle(f"Task 4 | Forced-Vibration Frequency Prediction | {file_path.name}")
|
||||
|
||||
axes[0].plot(time_segment, signal_segment, color="tab:blue", linewidth=1.0)
|
||||
axes[0].set_title(f"Steady middle segment from {TOP_RESPONSE_SENSOR} / {TOP_RESPONSE_AXIS}")
|
||||
axes[0].set_xlabel("Time (s)")
|
||||
axes[0].set_ylabel("Acceleration")
|
||||
axes[0].grid(True, alpha=0.25)
|
||||
|
||||
band_mask = (freqs >= 0.1) & (freqs <= 5.0)
|
||||
axes[1].plot(freqs[band_mask], magnitudes[band_mask], color="tab:purple", linewidth=1.1, label="Windowed FFT")
|
||||
axes[1].axvline(
|
||||
predicted_frequency_hz,
|
||||
color="tab:red",
|
||||
linestyle="--",
|
||||
linewidth=1.6,
|
||||
label=f"Predicted = {predicted_frequency_hz:.5f} Hz",
|
||||
)
|
||||
if true_frequency_hz is not None:
|
||||
axes[1].axvline(
|
||||
true_frequency_hz,
|
||||
color="tab:green",
|
||||
linestyle=":",
|
||||
linewidth=1.6,
|
||||
label=f"Filename truth = {true_frequency_hz:.5f} Hz",
|
||||
)
|
||||
axes[1].set_title("High-resolution spectrum from top steady-state response")
|
||||
axes[1].set_xlabel("Frequency (Hz)")
|
||||
axes[1].set_ylabel("Magnitude")
|
||||
axes[1].grid(True, alpha=0.25)
|
||||
axes[1].legend()
|
||||
|
||||
plt.tight_layout()
|
||||
figure_path = output_dir / f"{file_path.stem}_task4_predict_freq.png"
|
||||
plt.savefig(figure_path, dpi=180, bbox_inches="tight")
|
||||
return figure_path
|
||||
|
||||
|
||||
def pick_fundamental_peak_index(
|
||||
freqs: np.ndarray,
|
||||
magnitudes: np.ndarray,
|
||||
min_hz: float = 0.1,
|
||||
max_hz: float = 5.0,
|
||||
relative_threshold: float = 0.20,
|
||||
relative_prominence: float = 0.08,
|
||||
) -> int:
|
||||
band_mask = (freqs >= min_hz) & (freqs <= max_hz)
|
||||
band_indices = np.where(band_mask)[0]
|
||||
if band_indices.size == 0:
|
||||
raise ValueError("No FFT bins fall inside the requested frequency band.")
|
||||
|
||||
band_magnitudes = magnitudes[band_mask]
|
||||
max_magnitude = float(np.max(band_magnitudes))
|
||||
height_threshold = relative_threshold * max_magnitude
|
||||
prominence_threshold = relative_prominence * max_magnitude
|
||||
local_peaks, _ = find_peaks(
|
||||
band_magnitudes,
|
||||
height=height_threshold,
|
||||
prominence=prominence_threshold,
|
||||
)
|
||||
|
||||
if local_peaks.size == 0:
|
||||
# Fallback for very smooth spectra: choose the strongest bin in band.
|
||||
return int(band_indices[int(np.argmax(band_magnitudes))])
|
||||
|
||||
# Fundamental-first rule: choose the left-most valid peak.
|
||||
return int(band_indices[int(local_peaks[0])])
|
||||
|
||||
|
||||
def analyze_forced_vibration(
|
||||
file_path: Path,
|
||||
output_dir: Path,
|
||||
show: bool = False,
|
||||
) -> dict[str, float | int | str]:
|
||||
file_path = Path(file_path).resolve()
|
||||
time_values, raw_signal = read_sensor_signal(file_path, sensor_code=TOP_RESPONSE_SENSOR, value_column=TOP_RESPONSE_AXIS)
|
||||
time_segment, signal_segment = get_middle_segment(time_values, raw_signal)
|
||||
sampling_rate = estimate_sampling_rate(time_segment)
|
||||
_, freqs, magnitudes, _ = dominant_frequency_in_band(
|
||||
signal_segment,
|
||||
sampling_rate=sampling_rate,
|
||||
min_hz=0.1,
|
||||
max_hz=5.0,
|
||||
zero_padding_factor=32,
|
||||
)
|
||||
peak_index = pick_fundamental_peak_index(
|
||||
freqs=freqs,
|
||||
magnitudes=magnitudes,
|
||||
min_hz=0.1,
|
||||
max_hz=5.0,
|
||||
relative_threshold=0.20,
|
||||
relative_prominence=0.08,
|
||||
)
|
||||
predicted_frequency_hz = float(freqs[peak_index])
|
||||
true_frequency_hz = parse_frequency_from_filename(file_path.name)
|
||||
absolute_error_hz = None if true_frequency_hz is None else abs(predicted_frequency_hz - true_frequency_hz)
|
||||
|
||||
figure_path = save_task4_figure(
|
||||
file_path=file_path,
|
||||
output_dir=output_dir,
|
||||
time_segment=time_segment,
|
||||
signal_segment=signal_segment,
|
||||
freqs=freqs,
|
||||
magnitudes=magnitudes,
|
||||
predicted_frequency_hz=predicted_frequency_hz,
|
||||
true_frequency_hz=true_frequency_hz,
|
||||
)
|
||||
if show:
|
||||
plt.show()
|
||||
plt.close("all")
|
||||
|
||||
return {
|
||||
"file_path": str(file_path),
|
||||
"sensor_code": TOP_RESPONSE_SENSOR,
|
||||
"axis": TOP_RESPONSE_AXIS,
|
||||
"sampling_rate_hz": float(sampling_rate),
|
||||
"segment_sample_count": int(signal_segment.size),
|
||||
"segment_start_time": float(time_segment[0]),
|
||||
"segment_end_time": float(time_segment[-1]),
|
||||
"predicted_frequency_hz": float(predicted_frequency_hz),
|
||||
"fft_peak_bin_hz": float(freqs[peak_index]),
|
||||
"fft_peak_magnitude": float(magnitudes[peak_index]),
|
||||
"true_frequency_hz": None if true_frequency_hz is None else float(true_frequency_hz),
|
||||
"absolute_error_hz": None if absolute_error_hz is None else float(absolute_error_hz),
|
||||
"figure_path": str(figure_path),
|
||||
"note": "Prediction uses only the top response sensor WSMS00007/value3. Base sensor WSMS00012/value1 is not used to infer the answer.",
|
||||
}
|
||||
|
||||
|
||||
def print_report(result: dict[str, float | int | str]) -> None:
|
||||
print("=== Task 4: Excitation Frequency Reverse Prediction ===")
|
||||
print(f"Input file: {result['file_path']}")
|
||||
print(f"Top response sensor: {result['sensor_code']} / {result['axis']}")
|
||||
print(f"Sampling rate: {result['sampling_rate_hz']:.4f} Hz")
|
||||
print(f"Steady segment samples: {result['segment_sample_count']}")
|
||||
print(f"Steady segment time range: {result['segment_start_time']:.4f} s to {result['segment_end_time']:.4f} s")
|
||||
print(f"Predicted excitation frequency: {result['predicted_frequency_hz']:.6f} Hz")
|
||||
if result["true_frequency_hz"] is not None:
|
||||
print(f"Filename frequency truth: {result['true_frequency_hz']:.6f} Hz")
|
||||
print(f"Absolute error: {result['absolute_error_hz']:.6f} Hz")
|
||||
else:
|
||||
print("Filename frequency truth: unavailable")
|
||||
print(f"Figure saved to: {result['figure_path']}")
|
||||
print(result["note"])
|
||||
|
||||
|
||||
def main() -> None:
|
||||
args = parse_args()
|
||||
result = analyze_forced_vibration(
|
||||
file_path=Path(args.file),
|
||||
output_dir=Path(args.output_dir),
|
||||
show=bool(args.show),
|
||||
)
|
||||
print_report(result)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
45
scripts_4/README.md
Normal file
45
scripts_4/README.md
Normal file
@@ -0,0 +1,45 @@
|
||||
# scripts_4: TMD 小样本适配
|
||||
|
||||
本目录实现二阶段方案:
|
||||
|
||||
1. 基础模型(`Non_TMD` 训练得到的 `task1_final_model.pkl`)先给出基线预测;
|
||||
2. 用极少量 `TMD` 样本拟合频率相关比例系数 `k(f)`;
|
||||
3. 输出 `TMD` 预测:`y_tmd = k(f) * y_base`。
|
||||
|
||||
## 为什么这样做
|
||||
|
||||
- 仅 5 个 TMD 点直接重训主模型容易过拟合;
|
||||
- 保留基础模型泛化能力,只学习域偏移;
|
||||
- 通过 `shrinkage` 把校正系数向 `1.0` 收缩,降低小样本波动影响;
|
||||
- 频率外推时会自动向 `1.0` 回归,避免越界夸张。
|
||||
|
||||
## 脚本说明
|
||||
|
||||
- `train_tmd_adapter.py`:拟合并保存适配器,同时输出 LOOCV 报告。
|
||||
- `predict_tmd_single.py`:加载基础模型 + 适配器,对单个 TMD 文件推理。
|
||||
- `adapter.py`:频率校正器定义与序列化工具。
|
||||
|
||||
## 使用方法(请使用虚拟环境 Python)
|
||||
|
||||
```bash
|
||||
.venv/bin/python scripts_4/train_tmd_adapter.py
|
||||
```
|
||||
|
||||
```bash
|
||||
.venv/bin/python scripts_4/predict_tmd_single.py \
|
||||
--file downloads/TMD/val/harmonic_5mm_1.25Hz_TMD.csv
|
||||
```
|
||||
|
||||
可选参数示例:
|
||||
|
||||
```bash
|
||||
.venv/bin/python scripts_4/train_tmd_adapter.py \
|
||||
--shrinkage 0.25 \
|
||||
--extrapolation-decay-hz 0.3
|
||||
```
|
||||
|
||||
## 输出文件
|
||||
|
||||
- 适配器:`checkpoints_final/task1_final/task1_tmd_adapter.pkl`
|
||||
- 报告:`evaluation_outputs/task1_final/tmd_adapter_report.json`
|
||||
|
||||
2
scripts_4/__init__.py
Normal file
2
scripts_4/__init__.py
Normal file
@@ -0,0 +1,2 @@
|
||||
"""Task4: Tiny-sample TMD adaptation on top of baseline RMS model."""
|
||||
|
||||
BIN
scripts_4/__pycache__/__init__.cpython-314.pyc
Normal file
BIN
scripts_4/__pycache__/__init__.cpython-314.pyc
Normal file
Binary file not shown.
BIN
scripts_4/__pycache__/adapter.cpython-314.pyc
Normal file
BIN
scripts_4/__pycache__/adapter.cpython-314.pyc
Normal file
Binary file not shown.
101
scripts_4/adapter.py
Normal file
101
scripts_4/adapter.py
Normal file
@@ -0,0 +1,101 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from pathlib import Path
|
||||
import pickle
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
def _safe_float(value: float) -> float:
|
||||
return float(np.asarray(value, dtype=np.float64).reshape(-1)[0])
|
||||
|
||||
|
||||
@dataclass
|
||||
class FrequencyScaleAdapter:
|
||||
frequencies_hz: np.ndarray
|
||||
scale_values: np.ndarray
|
||||
shrinkage: float = 0.2
|
||||
extrapolation_decay_hz: float = 0.25
|
||||
min_scale: float = 0.4
|
||||
max_scale: float = 2.5
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
freqs = np.asarray(self.frequencies_hz, dtype=np.float64).reshape(-1)
|
||||
scales = np.asarray(self.scale_values, dtype=np.float64).reshape(-1)
|
||||
if freqs.size == 0:
|
||||
raise ValueError("frequencies_hz cannot be empty.")
|
||||
if freqs.size != scales.size:
|
||||
raise ValueError("frequencies_hz and scale_values must have the same length.")
|
||||
order = np.argsort(freqs)
|
||||
freqs = freqs[order]
|
||||
scales = np.clip(scales[order], self.min_scale, self.max_scale)
|
||||
# Pull scales toward 1.0 to reduce overfitting with tiny TMD samples.
|
||||
scales = 1.0 + (scales - 1.0) * (1.0 - float(np.clip(self.shrinkage, 0.0, 0.95)))
|
||||
self.frequencies_hz = freqs
|
||||
self.scale_values = scales
|
||||
|
||||
def _inside_range_weight(self, frequency_hz: float) -> float:
|
||||
left = self.frequencies_hz[0]
|
||||
right = self.frequencies_hz[-1]
|
||||
f = _safe_float(frequency_hz)
|
||||
if left <= f <= right:
|
||||
return 1.0
|
||||
distance = min(abs(f - left), abs(f - right))
|
||||
decay = max(_safe_float(self.extrapolation_decay_hz), 1e-6)
|
||||
return float(np.exp(-distance / decay))
|
||||
|
||||
def scale_at(self, frequency_hz: float) -> float:
|
||||
f = _safe_float(frequency_hz)
|
||||
boundary_interp = float(np.interp(f, self.frequencies_hz, self.scale_values))
|
||||
in_range_weight = self._inside_range_weight(f)
|
||||
scale = 1.0 + in_range_weight * (boundary_interp - 1.0)
|
||||
return float(np.clip(scale, self.min_scale, self.max_scale))
|
||||
|
||||
def predict(self, base_rms: float, frequency_hz: float) -> float:
|
||||
return max(_safe_float(base_rms), 0.0) * self.scale_at(frequency_hz)
|
||||
|
||||
def to_payload(self) -> dict:
|
||||
return {
|
||||
"frequencies_hz": self.frequencies_hz.tolist(),
|
||||
"scale_values": self.scale_values.tolist(),
|
||||
# scale_values are already post-shrinkage values after __post_init__.
|
||||
"shrinkage": 0.0,
|
||||
"extrapolation_decay_hz": float(self.extrapolation_decay_hz),
|
||||
"min_scale": float(self.min_scale),
|
||||
"max_scale": float(self.max_scale),
|
||||
"already_shrunk": True,
|
||||
}
|
||||
|
||||
@classmethod
|
||||
def from_payload(cls, payload: dict) -> "FrequencyScaleAdapter":
|
||||
shrinkage = float(payload.get("shrinkage", 0.2))
|
||||
if bool(payload.get("already_shrunk", False)):
|
||||
shrinkage = 0.0
|
||||
return cls(
|
||||
frequencies_hz=np.asarray(payload["frequencies_hz"], dtype=np.float64),
|
||||
scale_values=np.asarray(payload["scale_values"], dtype=np.float64),
|
||||
shrinkage=shrinkage,
|
||||
extrapolation_decay_hz=float(payload.get("extrapolation_decay_hz", 0.25)),
|
||||
min_scale=float(payload.get("min_scale", 0.4)),
|
||||
max_scale=float(payload.get("max_scale", 2.5)),
|
||||
)
|
||||
|
||||
|
||||
def save_adapter(adapter: FrequencyScaleAdapter, output_path: Path, extra: dict | None = None) -> None:
|
||||
output_path.parent.mkdir(parents=True, exist_ok=True)
|
||||
payload = {
|
||||
"adapter_type": "frequency_scale_interp_v1",
|
||||
"adapter": adapter.to_payload(),
|
||||
"extra": extra or {},
|
||||
}
|
||||
with output_path.open("wb") as handle:
|
||||
pickle.dump(payload, handle)
|
||||
|
||||
|
||||
def load_adapter(adapter_path: Path) -> tuple[FrequencyScaleAdapter, dict]:
|
||||
with adapter_path.open("rb") as handle:
|
||||
payload = pickle.load(handle)
|
||||
adapter = FrequencyScaleAdapter.from_payload(payload["adapter"])
|
||||
extra = payload.get("extra", {})
|
||||
return adapter, extra
|
||||
78
scripts_4/predict_tmd_single.py
Normal file
78
scripts_4/predict_tmd_single.py
Normal file
@@ -0,0 +1,78 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import importlib.util
|
||||
import pickle
|
||||
from pathlib import Path
|
||||
import sys
|
||||
|
||||
PROJECT_ROOT = Path(__file__).resolve().parents[1]
|
||||
|
||||
try:
|
||||
from .adapter import load_adapter
|
||||
except ImportError:
|
||||
from adapter import load_adapter
|
||||
|
||||
|
||||
def _load_module(module_name: str, file_path: Path):
|
||||
spec = importlib.util.spec_from_file_location(module_name, file_path)
|
||||
if spec is None or spec.loader is None:
|
||||
raise ImportError(f"Cannot load module from {file_path}")
|
||||
module = importlib.util.module_from_spec(spec)
|
||||
sys.modules[module_name] = module
|
||||
spec.loader.exec_module(module)
|
||||
return module
|
||||
|
||||
|
||||
def parse_args() -> argparse.Namespace:
|
||||
parser = argparse.ArgumentParser(description="Predict TMD RMS using baseline model + frequency adapter.")
|
||||
parser.add_argument("--file", type=str, required=True)
|
||||
parser.add_argument("--checkpoint", type=str, default="checkpoints_final/task1_final/task1_final_model.pkl")
|
||||
parser.add_argument("--adapter", type=str, default="checkpoints_final/task1_final/task1_tmd_adapter.pkl")
|
||||
parser.add_argument("--data-dir", type=str, default="downloads/TMD")
|
||||
return parser.parse_args()
|
||||
|
||||
|
||||
def main() -> None:
|
||||
args = parse_args()
|
||||
config_module = _load_module("config", PROJECT_ROOT / "scripts" / "config.py")
|
||||
dataset_module = _load_module("task1_dataset", PROJECT_ROOT / "scripts" / "dataset.py")
|
||||
make_experiment_config = getattr(config_module, "make_experiment_config")
|
||||
build_record_from_file = getattr(dataset_module, "build_record_from_file")
|
||||
|
||||
config = make_experiment_config()
|
||||
config.data.data_dir = args.data_dir
|
||||
config.data.__post_init__()
|
||||
|
||||
ckpt_path = (PROJECT_ROOT / args.checkpoint).resolve()
|
||||
adapter_path = (PROJECT_ROOT / args.adapter).resolve()
|
||||
with ckpt_path.open("rb") as handle:
|
||||
payload = pickle.load(handle)
|
||||
model = payload["model"]
|
||||
adapter, extra = load_adapter(adapter_path)
|
||||
|
||||
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_base = pred_tr * record.x_rms
|
||||
pred_tmd = adapter.predict(pred_base, record.frequency_hz)
|
||||
|
||||
true_rms = float(record.y_rms)
|
||||
base_err = abs(pred_base - true_rms) / max(abs(true_rms), 1e-12) * 100.0
|
||||
tmd_err = abs(pred_tmd - true_rms) / max(abs(true_rms), 1e-12) * 100.0
|
||||
|
||||
print(f"Checkpoint: {ckpt_path}")
|
||||
print(f"Adapter: {adapter_path}")
|
||||
print(f"Adapter extra: {extra}")
|
||||
print(f"Input file: {file_path}")
|
||||
print(f"Frequency (Hz): {record.frequency_hz:.6f}")
|
||||
print(f"x_rms: {record.x_rms:.6f}")
|
||||
print(f"True y_rms: {true_rms:.6f}")
|
||||
print(f"Base pred y_rms: {pred_base:.6f} | error={base_err:.4f}%")
|
||||
print(f"Adapted pred y_rms: {pred_tmd:.6f} | error={tmd_err:.4f}%")
|
||||
print(f"Applied scale k(f): {adapter.scale_at(record.frequency_hz):.6f}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
164
scripts_4/train_tmd_adapter.py
Normal file
164
scripts_4/train_tmd_adapter.py
Normal file
@@ -0,0 +1,164 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import json
|
||||
import pickle
|
||||
from pathlib import Path
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
PROJECT_ROOT = Path(__file__).resolve().parents[1]
|
||||
sys.path.insert(0, str(PROJECT_ROOT / "scripts"))
|
||||
|
||||
from config import make_experiment_config
|
||||
from dataset import build_record_from_file
|
||||
|
||||
try:
|
||||
from .adapter import FrequencyScaleAdapter, save_adapter
|
||||
except ImportError:
|
||||
from adapter import FrequencyScaleAdapter, save_adapter
|
||||
|
||||
|
||||
def parse_args() -> argparse.Namespace:
|
||||
parser = argparse.ArgumentParser(description="Fit a tiny-sample TMD adapter on top of the Non_TMD baseline .pkl model.")
|
||||
parser.add_argument("--checkpoint", type=str, default=None, help="Path to baseline model .pkl.")
|
||||
parser.add_argument("--tmd-dir", type=str, default="downloads/TMD", help="Directory that contains TMD harmonic csv files.")
|
||||
parser.add_argument("--output", type=str, default="checkpoints_final/task1_final/task1_tmd_adapter.pkl")
|
||||
parser.add_argument("--report-json", type=str, default="evaluation_outputs/task1_final/tmd_adapter_report.json")
|
||||
parser.add_argument("--shrinkage", type=float, default=0.20, help="Scale shrinkage toward 1.0.")
|
||||
parser.add_argument("--extrapolation-decay-hz", type=float, default=0.25)
|
||||
parser.add_argument("--min-scale", type=float, default=0.40)
|
||||
parser.add_argument("--max-scale", type=float, default=2.50)
|
||||
return parser.parse_args()
|
||||
|
||||
|
||||
def resolve_checkpoint_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 load_baseline_model(model_path: Path):
|
||||
with model_path.open("rb") as handle:
|
||||
payload = pickle.load(handle)
|
||||
return payload["model"], payload.get("summary", {})
|
||||
|
||||
|
||||
def base_predict_rms(model, features: np.ndarray, x_rms: float, eps: float) -> float:
|
||||
pred_tr = max(float(model.predict([features.tolist()])[0]), eps)
|
||||
return pred_tr * float(x_rms)
|
||||
|
||||
|
||||
def collect_tmd_rows(model, tmd_dir: Path, eps: float) -> list[dict]:
|
||||
config = make_experiment_config().data
|
||||
config.data_dir = str(Path(tmd_dir).relative_to(PROJECT_ROOT))
|
||||
config.__post_init__()
|
||||
rows: list[dict] = []
|
||||
for file_path in sorted(config.data_root.rglob("harmonic_5mm_*Hz_TMD.csv")):
|
||||
record = build_record_from_file(file_path, config)
|
||||
base_pred = base_predict_rms(model, record.features, record.x_rms, eps)
|
||||
scale = float(record.y_rms / max(base_pred, eps))
|
||||
rows.append(
|
||||
{
|
||||
"file_name": file_path.name,
|
||||
"frequency_hz": float(record.frequency_hz),
|
||||
"true_y_rms": float(record.y_rms),
|
||||
"base_pred_y_rms": float(base_pred),
|
||||
"scale_true_over_base": float(scale),
|
||||
}
|
||||
)
|
||||
if not rows:
|
||||
raise RuntimeError("No TMD harmonic files found.")
|
||||
return rows
|
||||
|
||||
|
||||
def loocv_error(rows: list[dict], shrinkage: float, extrapolation_decay_hz: float, min_scale: float, max_scale: float) -> list[float]:
|
||||
errors: list[float] = []
|
||||
eps = 1e-12
|
||||
for holdout in range(len(rows)):
|
||||
train_rows = [row for i, row in enumerate(rows) if i != holdout]
|
||||
adapter = FrequencyScaleAdapter(
|
||||
frequencies_hz=np.asarray([row["frequency_hz"] for row in train_rows], dtype=np.float64),
|
||||
scale_values=np.asarray([row["scale_true_over_base"] for row in train_rows], dtype=np.float64),
|
||||
shrinkage=shrinkage,
|
||||
extrapolation_decay_hz=extrapolation_decay_hz,
|
||||
min_scale=min_scale,
|
||||
max_scale=max_scale,
|
||||
)
|
||||
sample = rows[holdout]
|
||||
pred = adapter.predict(sample["base_pred_y_rms"], sample["frequency_hz"])
|
||||
err = abs(pred - sample["true_y_rms"]) / max(abs(sample["true_y_rms"]), eps) * 100.0
|
||||
errors.append(float(err))
|
||||
return errors
|
||||
|
||||
|
||||
def main() -> None:
|
||||
args = parse_args()
|
||||
ckpt_path = resolve_checkpoint_path(args.checkpoint)
|
||||
model, baseline_summary = load_baseline_model(ckpt_path)
|
||||
rows = collect_tmd_rows(model, (PROJECT_ROOT / args.tmd_dir).resolve(), eps=1e-6)
|
||||
rows = sorted(rows, key=lambda x: x["frequency_hz"])
|
||||
|
||||
adapter = FrequencyScaleAdapter(
|
||||
frequencies_hz=np.asarray([row["frequency_hz"] for row in rows], dtype=np.float64),
|
||||
scale_values=np.asarray([row["scale_true_over_base"] for row in rows], dtype=np.float64),
|
||||
shrinkage=args.shrinkage,
|
||||
extrapolation_decay_hz=args.extrapolation_decay_hz,
|
||||
min_scale=args.min_scale,
|
||||
max_scale=args.max_scale,
|
||||
)
|
||||
|
||||
adapted_rows = []
|
||||
for row in rows:
|
||||
adapted_pred = adapter.predict(row["base_pred_y_rms"], row["frequency_hz"])
|
||||
adapted_err = abs(adapted_pred - row["true_y_rms"]) / max(abs(row["true_y_rms"]), 1e-12) * 100.0
|
||||
adapted_rows.append({**row, "adapted_pred_y_rms": adapted_pred, "adapted_err_pct": adapted_err})
|
||||
|
||||
loocv_errors = loocv_error(
|
||||
rows=rows,
|
||||
shrinkage=args.shrinkage,
|
||||
extrapolation_decay_hz=args.extrapolation_decay_hz,
|
||||
min_scale=args.min_scale,
|
||||
max_scale=args.max_scale,
|
||||
)
|
||||
report = {
|
||||
"checkpoint": str(ckpt_path),
|
||||
"baseline_summary": baseline_summary,
|
||||
"samples": adapted_rows,
|
||||
"adapter_scales_after_shrinkage": [
|
||||
{"frequency_hz": float(f), "scale": float(s)}
|
||||
for f, s in zip(adapter.frequencies_hz.tolist(), adapter.scale_values.tolist())
|
||||
],
|
||||
"metrics": {
|
||||
"base_mean_err_pct": float(
|
||||
np.mean(
|
||||
[
|
||||
abs(row["base_pred_y_rms"] - row["true_y_rms"]) / max(abs(row["true_y_rms"]), 1e-12) * 100.0
|
||||
for row in rows
|
||||
]
|
||||
)
|
||||
),
|
||||
"adapted_mean_err_pct": float(np.mean([row["adapted_err_pct"] for row in adapted_rows])),
|
||||
"loocv_mean_err_pct": float(np.mean(loocv_errors)),
|
||||
"loocv_max_err_pct": float(np.max(loocv_errors)),
|
||||
},
|
||||
}
|
||||
|
||||
output_path = (PROJECT_ROOT / args.output).resolve()
|
||||
save_adapter(adapter, output_path, extra={"report_metrics": report["metrics"], "checkpoint": str(ckpt_path)})
|
||||
|
||||
report_path = (PROJECT_ROOT / args.report_json).resolve()
|
||||
report_path.parent.mkdir(parents=True, exist_ok=True)
|
||||
report_path.write_text(json.dumps(report, ensure_ascii=False, indent=2), encoding="utf-8")
|
||||
|
||||
print(f"Adapter saved to: {output_path}")
|
||||
print(f"Report saved to: {report_path}")
|
||||
print(pd.DataFrame(adapted_rows).to_string(index=False))
|
||||
print(f"LOOCV mean error (%): {report['metrics']['loocv_mean_err_pct']:.4f}")
|
||||
print(f"LOOCV max error (%): {report['metrics']['loocv_max_err_pct']:.4f}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
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