python
# Create the single all-in-one script requested by the user
from pathlib import Path

script_path = Path("/mnt/data/fecam_tcn_itransformer_single.py")
code = r'''# -*- coding: utf-8 -*-
"""
单文件版：FECAM-based TCN + iTransformer（仅新模型）
包含：训练(train)、预测(predict)、误差评估(eval)
---------------------------------------------------
你当前的数据格式沿用你原来的脚本：
- 训练：Input_GNSS.csv + Input_INS.csv
- 测试：出城INS.csv（包含 GPS周内秒(s), 纬度(°), 经度(°) 以及12维INS特征）
- 预测目标：GNSS相邻点的北向/东向位移增量（delta_north, delta_east），单位 m

关键变化（按你要求）：
- 输入从“固定20点窗口”改为：每个输出时刻（GNSS增量时刻）前 3 秒内的所有12维INS数据
  为了能喂给TCN/Transformer，仍需固定长度：脚本会根据 INS 的采样间隔自动计算 3s 对应的样本数 L
  训练和预测都使用同一个 L（训练阶段计算并保存到 meta.json）

用法（推荐）：
1) 训练 + 预测 + 误差评估（默认）:
   python fecam_tcn_itransformer_single.py --mode all

2) 只训练:
   python fecam_tcn_itransformer_single.py --mode train

3) 只预测并输出误差（需要已有模型与scaler）:
   python fecam_tcn_itransformer_single.py --mode eval

注意：
- BASE_DIR 请改成你电脑上的工程目录
"""

import os
import json
import math
import argparse
from pathlib import Path

import numpy as np
import pandas as pd
import joblib
import matplotlib.pyplot as plt

import torch
import torch.nn as nn
import torch.optim as optim
from torch.nn.utils import weight_norm
from torch.utils.data import DataLoader, TensorDataset, random_split

from sklearn.preprocessing import MinMaxScaler


# ============================================================
# 1)  模型：FECAM-based TCN + iTransformer
# ============================================================

class CausalConv1d(nn.Module):
    """因果卷积：只在左侧padding，避免看见未来。"""
    def __init__(self, in_ch: int, out_ch: int, kernel_size: int, dilation: int = 1):
        super().__init__()
        self.kernel_size = kernel_size
        self.dilation = dilation
        self.pad = (kernel_size - 1) * dilation
        self.conv = weight_norm(nn.Conv1d(in_ch, out_ch, kernel_size=kernel_size, dilation=dilation))

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (B, C, T)
        x = nn.functional.pad(x, (self.pad, 0))  # 左pad
        return self.conv(x)


class TCNResidualUnit(nn.Module):
    """标准TCN残差单元：两层Dilated Causal Conv + ReLU + Dropout + 残差。"""
    def __init__(self, in_ch: int, out_ch: int, kernel_size: int, dilation: int, dropout: float):
        super().__init__()
        self.conv1 = CausalConv1d(in_ch, out_ch, kernel_size, dilation)
        self.relu1 = nn.ReLU()
        self.drop1 = nn.Dropout(dropout)

        self.conv2 = CausalConv1d(out_ch, out_ch, kernel_size, dilation)
        self.relu2 = nn.ReLU()
        self.drop2 = nn.Dropout(dropout)

        self.downsample = None
        if in_ch != out_ch:
            self.downsample = nn.Conv1d(in_ch, out_ch, kernel_size=1)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (B, C, T)
        y = self.conv1(x)
        y = self.relu1(y)
        y = self.drop1(y)

        y = self.conv2(y)
        y = self.relu2(y)
        y = self.drop2(y)

        res = x if self.downsample is None else self.downsample(x)
        return y + res


def dct_2(x: torch.Tensor) -> torch.Tensor:
    """
    DCT-II（近似实现）沿最后一维做DCT，用FFT的偶延拓实现。
    x: (..., T)
    return: (..., T)
    """
    N = x.size(-1)
    v = torch.cat([x, x.flip(dims=[-1])], dim=-1)  # (..., 2N)
    Vc = torch.fft.rfft(v, dim=-1)  # (..., N+1)
    k = torch.arange(N, device=x.device, dtype=torch.float32)
    factor = torch.exp(-1j * math.pi * k / (2.0 * N))
    # 取前N项并取实部
    out = (Vc[..., :N] * factor).real
    return out * 2.0


class FECAM(nn.Module):
    """
    Frequency-Enhanced Channel Attention Mechanism（简化实现）：
    - 对每个通道做DCT，取前freq_bins个频率分量作为描述
    - MLP 生成通道权重（类似SE，但用频域描述）
    - 小FCNN（1x1 conv）生成增强张量，并与输入残差相加
    """
    def __init__(self, channels: int, freq_bins: int = 16, mlp_ratio: float = 0.5, dropout: float = 0.1):
        super().__init__()
        self.channels = channels
        self.freq_bins = freq_bins
        hidden = max(8, int(channels * mlp_ratio))

        self.mlp = nn.Sequential(
            nn.Linear(channels * freq_bins, hidden),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(hidden, channels),
            nn.Sigmoid()
        )

        self.fcnn = nn.Sequential(
            nn.Conv1d(channels, channels, kernel_size=1),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Conv1d(channels, channels, kernel_size=1)
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (B, C, T)
        B, C, T = x.shape
        F = min(self.freq_bins, T)
        # DCT沿时间维
        dct = dct_2(x)  # (B,C,T)
        desc = dct[..., :F]  # (B,C,F)
        desc = desc.reshape(B, C * F)  # (B, C*F)
        w = self.mlp(desc).unsqueeze(-1)  # (B,C,1)
        x_scaled = x * w
        enh = self.fcnn(x_scaled)  # (B,C,T)
        return x + enh


class FECAM_TCN(nn.Module):
    """
    FECAM-based TCN：两个TCN堆栈，中间插FECAM（对应论文的“两个TCN residual blocks + FECAM串联”）
    """
    def __init__(
        self,
        in_dim: int,
        channels: int,
        kernel_size: int,
        n_layers_per_block: int,
        dropout_tcn: float,
        dropout_fecam: float,
        freq_bins: int = 16
    ):
        super().__init__()
        self.in_proj = nn.Conv1d(in_dim, channels, kernel_size=1)

        def make_stack():
            layers = []
            for i in range(n_layers_per_block):
                dilation = 2 ** i
                layers.append(TCNResidualUnit(
                    in_ch=channels,
                    out_ch=channels,
                    kernel_size=kernel_size,
                    dilation=dilation,
                    dropout=dropout_tcn
                ))
            return nn.Sequential(*layers)

        self.tcn_block1 = make_stack()
        self.fecam = FECAM(channels=channels, freq_bins=freq_bins, dropout=dropout_fecam)
        self.tcn_block2 = make_stack()

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (B, T, D) -> (B, C, T)
        x = x.transpose(1, 2)  # (B,D,T)
        x = self.in_proj(x)     # (B,C,T)
        x = self.tcn_block1(x)
        x = self.fecam(x)
        x = self.tcn_block2(x)
        return x  # (B,C,T)


class iTransformerEncoderLayer(nn.Module):
    def __init__(self, d_model: int, n_heads: int, dropout: float = 0.1, ff_ratio: float = 4.0):
        super().__init__()
        self.attn = nn.MultiheadAttention(d_model, n_heads, dropout=dropout, batch_first=True)
        self.norm1 = nn.LayerNorm(d_model)
        self.drop1 = nn.Dropout(dropout)

        ff_hidden = int(d_model * ff_ratio)
        self.ffn = nn.Sequential(
            nn.Linear(d_model, ff_hidden),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(ff_hidden, d_model),
        )
        self.norm2 = nn.LayerNorm(d_model)
        self.drop2 = nn.Dropout(dropout)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (B, V, d_model)  V=变量token数（这里就是TCN通道数）
        attn_out, _ = self.attn(x, x, x, need_weights=False)
        x = self.norm1(x + self.drop1(attn_out))
        ff_out = self.ffn(x)
        x = self.norm2(x + self.drop2(ff_out))
        return x


class iTransformer(nn.Module):
    """
    iTransformer（简化版）：
    - 输入 (B,C,T) ：每个通道的整段时间序列 -> 变量token
    - Embedding：Linear(T -> d_model) 逐变量
    - Encoder-only：在变量token维度做自注意力
    - Projection：聚合变量token -> 输出(2)
    """
    def __init__(self, seq_len: int, n_vars: int, d_model: int = 128, n_heads: int = 8, n_layers: int = 2, dropout: float = 0.1):
        super().__init__()
        self.seq_len = seq_len
        self.n_vars = n_vars
        self.embed = nn.Linear(seq_len, d_model)
        self.enc_layers = nn.ModuleList([
            iTransformerEncoderLayer(d_model=d_model, n_heads=n_heads, dropout=dropout)
            for _ in range(n_layers)
        ])
        self.proj = nn.Sequential(
            nn.LayerNorm(d_model),
            nn.Linear(d_model, d_model),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(d_model, 2)
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (B, C, T)
        B, C, T = x.shape
        if T != self.seq_len:
            raise ValueError(f"iTransformer expects seq_len={self.seq_len}, got T={T}")
        # 变量token：每个变量一条长度T的序列
        tokens = self.embed(x)  # (B,C,d_model)
        for layer in self.enc_layers:
            tokens = layer(tokens)  # (B,C,d_model)
        pooled = tokens.mean(dim=1)  # (B,d_model)
        return self.proj(pooled)     # (B,2)


class FECAM_TCN_iTransformer(nn.Module):
    """
    总模型：FECAM-based TCN -> iTransformer -> (delta_north, delta_east)
    """
    def __init__(
        self,
        input_dim: int,
        seq_len: int,
        tcn_channels: int = 64,
        tcn_kernel: int = 3,
        tcn_layers_per_block: int = 5,
        dropout_tcn: float = 0.2,
        dropout_fecam: float = 0.1,
        itr_d_model: int = 128,
        itr_heads: int = 8,
        itr_layers: int = 2,
        dropout_itr: float = 0.1,
        fecam_freq_bins: int = 16
    ):
        super().__init__()
        self.backbone = FECAM_TCN(
            in_dim=input_dim,
            channels=tcn_channels,
            kernel_size=tcn_kernel,
            n_layers_per_block=tcn_layers_per_block,
            dropout_tcn=dropout_tcn,
            dropout_fecam=dropout_fecam,
            freq_bins=fecam_freq_bins
        )
        self.head = iTransformer(
            seq_len=seq_len,
            n_vars=tcn_channels,
            d_model=itr_d_model,
            n_heads=itr_heads,
            n_layers=itr_layers,
            dropout=dropout_itr
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (B,T,D)
        feat = self.backbone(x)  # (B,C,T)
        out = self.head(feat)    # (B,2)
        return out


# ============================================================
# 2)  数据处理：按“输出时刻前3秒”构造样本
# ============================================================

FEATURE_COLS = [
    '北向速度(m/s)', '东向速度(m/s)', '地向速度(m/s)',
    '航向角(°)', '俯仰角(°)', '横滚角(°)',
    '加表x（g）', '加表y（g）', '加表z（g）',
    '陀螺x（°/s）', '陀螺y（°/s）', '陀螺z（°/s）'
]

EARTH_RADIUS = 6371000.0


def median_dt(times: np.ndarray) -> float:
    d = np.diff(times.astype(np.float64))
    d = d[np.isfinite(d) & (d > 0)]
    if len(d) == 0:
        raise ValueError("时间列无法计算采样间隔（diff为空/非正）。")
    return float(np.median(d))


def compute_gnss_deltas(gnss_df: pd.DataFrame) -> pd.DataFrame:
    """从GNSS纬经度算北向/东向位移增量（米），并保留对应的“增量时刻”=后一时刻。"""
    df = gnss_df.sort_values('GNSS周内秒(s)').copy()
    lat = df['GNSS纬度(°)'].to_numpy(np.float64)
    lon = df['GNSS经度(°)'].to_numpy(np.float64)
    t = df['GNSS周内秒(s)'].to_numpy(np.float64)

    dlat = np.diff(lat) * (np.pi / 180.0)
    dlon = np.diff(lon) * (np.pi / 180.0)
    lat_mid = lat[1:]  # 用后一时刻的纬度近似
    dn = dlat * EARTH_RADIUS
    de = dlon * EARTH_RADIUS * np.cos(lat_mid * np.pi / 180.0)

    out = pd.DataFrame({
        't_gnss': t[1:],
        'delta_north': dn,
        'delta_east': de
    })
    return out


def build_windows_from_ins(
    ins_times: np.ndarray,
    ins_features: np.ndarray,
    target_times: np.ndarray,
    window_len: int
) -> np.ndarray:
    """
    对每个target_time，取“<= target_time”的最后window_len个INS样本作为输入窗口。
    若不足则左侧padding（用首样本重复）。
    返回: (N, window_len, D)
    """
    N = len(target_times)
    D = ins_features.shape[1]
    X = np.zeros((N, window_len, D), dtype=np.float32)

    for i, t in enumerate(target_times):
        end_idx = np.searchsorted(ins_times, t, side='right') - 1
        if end_idx < 0:
            # 全padding为0
            continue
        start_idx = end_idx - window_len + 1
        if start_idx >= 0:
            X[i] = ins_features[start_idx:end_idx + 1]
        else:
            pad = -start_idx
            # 左侧padding：用第0行重复（更稳一点，也可以改成0）
            X[i, :pad] = ins_features[0]
            X[i, pad:] = ins_features[0:end_idx + 1]
    return X


# ============================================================
# 3)  训练 / 预测 / 误差
# ============================================================

def set_seed(seed: int = 42):
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)


def train_only(args):
    base_dir = Path(args.base_dir)
    base_dir.mkdir(parents=True, exist_ok=True)

    device = torch.device("cuda" if torch.cuda.is_available() and not args.force_cpu else "cpu")
    print(f"[Train] Using device: {device}")

    # 读取数据
    gnss_data = pd.read_csv(base_dir / args.gnss_file)
    ins_data = pd.read_csv(base_dir / args.ins_file)

    # 训练时间范围（你要求：输出仍是 199900-200000 的GNSS增量）
    gnss_filtered = gnss_data[(gnss_data['GNSS周内秒(s)'] >= args.train_gnss_start) &
                              (gnss_data['GNSS周内秒(s)'] <= args.train_gnss_end)].copy()

    # INS 需要额外往前扩 3 秒（否则前几条GNSS增量没法取满窗口）
    ins_start = args.train_gnss_start - args.window_seconds - 0.5
    ins_end = args.train_gnss_end + 0.5
    ins_filtered = ins_data[(ins_data['GPS周内秒(s)'] >= ins_start) &
                            (ins_data['GPS周内秒(s)'] <= ins_end)].copy()
    ins_filtered = ins_filtered.sort_values('GPS周内秒(s)')

    # 计算GNSS增量标签
    deltas_df = compute_gnss_deltas(gnss_filtered)
    y = deltas_df[['delta_north', 'delta_east']].to_numpy(np.float32)
    t_targets = deltas_df['t_gnss'].to_numpy(np.float64)

    # INS特征
    ins_times = ins_filtered['GPS周内秒(s)'].to_numpy(np.float64)
    ins_feat = ins_filtered[FEATURE_COLS].to_numpy(np.float32)

    # 自动算 3 秒窗口对应长度
    dt_ins = median_dt(ins_times)
    window_len = int(round(args.window_seconds / dt_ins))
    window_len = max(window_len, 8)  # 给个下限
    # 为了保证“包含整3秒”，+1 更保险
    window_len = window_len + 1

    # GNSS采样间隔（用于预测时的步长估计）
    dt_gnss = median_dt(gnss_filtered.sort_values('GNSS周内秒(s)')['GNSS周内秒(s)'].to_numpy(np.float64))
    step_samples = int(max(1, round(dt_gnss / dt_ins)))

    print(f"[Train] dt_ins≈{dt_ins:.6f}s, window_seconds={args.window_seconds}s => window_len={window_len} samples")
    print(f"[Train] dt_gnss≈{dt_gnss:.6f}s, step_samples≈{step_samples}")

    # 构造输入窗口
    X = build_windows_from_ins(ins_times, ins_feat, t_targets, window_len=window_len)

    # 去掉“几乎全padding”的样本（可选）：end_idx<0 或有效长度太短
    # 这里用一个简单策略：窗口中是否存在非零（如果你INS数据本身可能为0，可关闭）
    valid_mask = np.isfinite(X).all(axis=(1, 2))
    X = X[valid_mask]
    y = y[valid_mask]
    t_targets = t_targets[valid_mask]
    print(f"[Train] 样本数: {len(X)}")

    # 归一化
    feature_scaler = MinMaxScaler()
    X_reshape = X.reshape(-1, X.shape[-1])
    X_norm = feature_scaler.fit_transform(X_reshape).reshape(X.shape)

    label_scaler = MinMaxScaler()
    y_norm = label_scaler.fit_transform(y)

    # 保存scaler + meta
    joblib.dump(feature_scaler, base_dir / args.feature_scaler_out)
    joblib.dump(label_scaler, base_dir / args.label_scaler_out)
    meta = {
        "window_seconds": args.window_seconds,
        "dt_ins": dt_ins,
        "dt_gnss": dt_gnss,
        "window_len": window_len,
        "step_samples": step_samples,
        "feature_cols": FEATURE_COLS,
        "input_dim": len(FEATURE_COLS)
    }
    (base_dir / args.meta_out).write_text(json.dumps(meta, ensure_ascii=False, indent=2), encoding="utf-8")
    print(f"[Train] Saved scalers + meta: {args.feature_scaler_out}, {args.label_scaler_out}, {args.meta_out}")

    # Dataset
    X_tensor = torch.tensor(X_norm, dtype=torch.float32)
    y_tensor = torch.tensor(y_norm, dtype=torch.float32)
    dataset = TensorDataset(X_tensor, y_tensor)

    train_size = int((1.0 - args.val_ratio) * len(dataset))
    val_size = len(dataset) - train_size
    train_ds, val_ds = random_split(dataset, [train_size, val_size], generator=torch.Generator().manual_seed(args.seed))

    train_loader = DataLoader(train_ds, batch_size=args.batch_size, shuffle=True, drop_last=False)
    val_loader = DataLoader(val_ds, batch_size=args.batch_size, shuffle=False, drop_last=False)

    # Model
    model = FECAM_TCN_iTransformer(
        input_dim=len(FEATURE_COLS),
        seq_len=window_len,
        tcn_channels=args.tcn_channels,
        tcn_kernel=args.tcn_kernel,
        tcn_layers_per_block=args.tcn_layers_per_block,
        dropout_tcn=args.dropout_tcn,
        dropout_fecam=args.dropout_fecam,
        itr_d_model=args.itr_d_model,
        itr_heads=args.itr_heads,
        itr_layers=args.itr_layers,
        dropout_itr=args.dropout_itr,
        fecam_freq_bins=args.fecam_freq_bins
    ).to(device)

    criterion = nn.MSELoss()
    optimizer = optim.Adam(model.parameters(), lr=args.lr, weight_decay=args.weight_decay)

    train_losses, val_losses = [], []
    best_val = float("inf")

    for epoch in range(args.epochs):
        model.train()
        total = 0.0
        for xb, yb in train_loader:
            xb = xb.to(device)
            yb = yb.to(device)

            optimizer.zero_grad(set_to_none=True)
            pred = model(xb)
            loss = criterion(pred, yb)
            loss.backward()

            if args.grad_clip > 0:
                nn.utils.clip_grad_norm_(model.parameters(), args.grad_clip)

            optimizer.step()
            total += float(loss.item())
        train_loss = total / max(1, len(train_loader))

        model.eval()
        vtotal = 0.0
        with torch.no_grad():
            for xb, yb in val_loader:
                xb = xb.to(device)
                yb = yb.to(device)
                pred = model(xb)
                vtotal += float(criterion(pred, yb).item())
        val_loss = vtotal / max(1, len(val_loader))

        train_losses.append(train_loss)
        val_losses.append(val_loss)

        if val_loss < best_val:
            best_val = val_loss
            torch.save(model.state_dict(), base_dir / args.model_out)

        if (epoch + 1) % args.log_every == 0 or epoch == 0:
            print(f"[Train] Epoch {epoch+1:03d}/{args.epochs} | train={train_loss:.6f} | val={val_loss:.6f} | best={best_val:.6f}")

    # 训练曲线
    if args.save_plots:
        plt.figure(figsize=(8, 4))
        plt.plot(train_losses, label="Train Loss")
        plt.plot(val_losses, label="Val Loss")
        plt.xlabel("Epoch")
        plt.ylabel("MSE Loss")
        plt.title("FECAM-TCN-iTransformer Training Curve")
        plt.grid(True, alpha=0.3)
        plt.legend()
        plt.tight_layout()
        plt.savefig(base_dir / args.train_curve_out, dpi=150)
        plt.close()
        print(f"[Train] Saved training curve: {args.train_curve_out}")

    print(f"[Train] Done. Best model saved to: {args.model_out}")
    return meta


def deltas_to_coords(lat, lon, dn, de):
    dlat = dn / EARTH_RADIUS * (180.0 / np.pi)
    dlon = de / (EARTH_RADIUS * np.cos(lat * np.pi / 180.0)) * (180.0 / np.pi)
    return lat + dlat, lon + dlon


def vincenty_offset(lat1, lon1, lat2, lon2):
    """
    Vincenty（简化迭代）计算从点1到点2的north/east offset（米）
    用法与你 compare_errors.py 一致。
    """
    from math import radians, sin, cos, sqrt, atan2

    a = 6378137.0
    f = 1/298.257223563
    b = a * (1 - f)

    lat1, lon1, lat2, lon2 = map(radians, [lat1, lon1, lat2, lon2])

    U1 = atan2((1 - f) * sin(lat1), cos(lat1))
    U2 = atan2((1 - f) * sin(lat2), cos(lat2))
    L = lon2 - lon1

    sinU1, cosU1 = sin(U1), cos(U1)
    sinU2, cosU2 = sin(U2), cos(U2)

    lam = L
    for _ in range(100):
        sinlam = sin(lam)
        coslam = cos(lam)
        sinsig = sqrt((cosU2 * sinlam) ** 2 +
                      (cosU1 * sinU2 - sinU1 * cosU2 * coslam) ** 2)
        if sinsig == 0:
            return 0.0, 0.0
        cossig = sinU1 * sinU2 + cosU1 * cosU2 * coslam
        sig = atan2(sinsig, cossig)
        sinalpha = cosU1 * cosU2 * sinlam / sinsig
        cos2alpha = 1 - sinalpha ** 2
        if cos2alpha == 0:
            cos2sigm = 0
        else:
            cos2sigm = cossig - 2 * sinU1 * sinU2 / cos2alpha
        C = f / 16 * cos2alpha * (4 + f * (4 - 3 * cos2alpha))
        lam_new = L + (1 - C) * f * sinalpha * (
            sig + C * sinsig * (cos2sigm + C * cossig * (-1 + 2 * cos2sigm ** 2)))
        if abs(lam_new - lam) < 1e-12:
            lam = lam_new
            break
        lam = lam_new

    az = atan2(cosU2 * sin(lam), cosU1 * sinU2 - sinU1 * cosU2 * cos(lam))
    dist = b * sig
    north = dist * cos(az)
    east = dist * sin(az)
    return north, east


def eval_on_test(args):
    base_dir = Path(args.base_dir)
    device = torch.device("cuda" if torch.cuda.is_available() and not args.force_cpu else "cpu")
    print(f"[Eval] Using device: {device}")

    # load meta + scaler
    meta = json.loads((base_dir / args.meta_out).read_text(encoding="utf-8"))
    window_len = int(meta["window_len"])
    step_samples = int(meta["step_samples"])

    feature_scaler = joblib.load(base_dir / args.feature_scaler_out)
    label_scaler = joblib.load(base_dir / args.label_scaler_out)

    # load test data
    test_df = pd.read_csv(base_dir / args.test_file)
    test_df = test_df[(test_df['GPS周内秒(s)'] >= args.test_start) &
                      (test_df['GPS周内秒(s)'] <= args.test_end)].copy()
    test_df = test_df.sort_values('GPS周内秒(s)').reset_index(drop=True)

    if len(test_df) < window_len + step_samples + 1:
        raise ValueError(f"测试数据太短：len={len(test_df)}, 需要至少 window_len({window_len})+step_samples({step_samples})")

    times = test_df['GPS周内秒(s)'].to_numpy(np.float64)
    feat = test_df[FEATURE_COLS].to_numpy(np.float32)

    # 输出索引：每 step_samples 输出一次（对应GNSS增量周期），且每次输入是前window_len样本
    end_indices = list(range(window_len - 1, len(test_df), step_samples))

    X_list = []
    out_times = []
    base_coords = []
    true_coords = []

    for end_idx in end_indices:
        prev_idx = end_idx - step_samples
        if prev_idx < 0:
            continue
        start_idx = end_idx - window_len + 1
        if start_idx < 0:
            continue
        X_list.append(feat[start_idx:end_idx + 1])
        out_times.append(times[end_idx])
        base_coords.append([test_df.loc[prev_idx, '纬度(°)'], test_df.loc[prev_idx, '经度(°)']])
        true_coords.append([test_df.loc[end_idx, '纬度(°)'], test_df.loc[end_idx, '经度(°)']])

    X = np.array(X_list, dtype=np.float32)  # (N,L,12)

    # normalize
    Xn = feature_scaler.transform(X.reshape(-1, X.shape[-1])).reshape(X.shape)
    X_tensor = torch.tensor(Xn, dtype=torch.float32).to(device)

    # model
    model = FECAM_TCN_iTransformer(
        input_dim=len(FEATURE_COLS),
        seq_len=window_len,
        tcn_channels=args.tcn_channels,
        tcn_kernel=args.tcn_kernel,
        tcn_layers_per_block=args.tcn_layers_per_block,
        dropout_tcn=args.dropout_tcn,
        dropout_fecam=args.dropout_fecam,
        itr_d_model=args.itr_d_model,
        itr_heads=args.itr_heads,
        itr_layers=args.itr_layers,
        dropout_itr=args.dropout_itr,
        fecam_freq_bins=args.fecam_freq_bins
    ).to(device)

    model.load_state_dict(torch.load(base_dir / args.model_out, map_location=device))
    model.eval()

    with torch.no_grad():
        pred_norm = model(X_tensor).cpu().numpy()
    pred = label_scaler.inverse_transform(pred_norm)  # (N,2)

    # coords update + errors
    updated_coords = []
    north_err = []
    east_err = []
    total_err = []

    for (lat0, lon0), (lat_true, lon_true), (dn, de) in zip(base_coords, true_coords, pred):
        lat_upd, lon_upd = deltas_to_coords(lat0, lon0, float(dn), float(de))
        updated_coords.append([lat_upd, lon_upd])

        n_off, e_off = vincenty_offset(lat_true, lon_true, lat_upd, lon_upd)
        north_err.append(n_off)
        east_err.append(e_off)
        total_err.append(math.sqrt(n_off ** 2 + e_off ** 2))

    results = pd.DataFrame({
        '时间戳': out_times,
        '基准纬度': [c[0] for c in base_coords],
        '基准经度': [c[1] for c in base_coords],
        '真实纬度': [c[0] for c in true_coords],
        '真实经度': [c[1] for c in true_coords],
        '预测北向增量(m)': pred[:, 0],
        '预测东向增量(m)': pred[:, 1],
        '更新后纬度': [c[0] for c in updated_coords],
        '更新后经度': [c[1] for c in updated_coords],
        '北向误差(m)': north_err,
        '东向误差(m)': east_err,
        '总误差(m)': total_err
    })

    out_csv = base_dir / args.pred_out
    results.to_csv(out_csv, index=False)
    print(f"[Eval] Saved predictions+errors: {args.pred_out}")

    # stats
    te = np.array(total_err, dtype=np.float64)
    stats = {
        "模型": "FECAM_TCN_iTransformer",
        "最小误差(m)": float(np.min(te)),
        "最大误差(m)": float(np.max(te)),
        "平均误差(m)": float(np.mean(te)),
        "标准差(m)": float(np.std(te, ddof=1)) if len(te) > 1 else 0.0,
        "RMSE(m)": float(np.sqrt(np.mean(te ** 2)))
    }
    pd.DataFrame([stats]).to_csv(base_dir / args.stats_out, index=False, encoding="utf-8-sig")
    print("[Eval] Error stats:")
    for k, v in stats.items():
        if k == "模型":
            print(f"  {k}: {v}")
        else:
            print(f"  {k}: {v:.6f}")

    # plot
    if args.save_plots:
        plt.figure(figsize=(10, 4))
        plt.plot(results.index, results['总误差(m)'].values)
        plt.xlabel("时间步")
        plt.ylabel("总误差 (m)")
        plt.title("FECAM-TCN-iTransformer Total Error")
        plt.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.savefig(base_dir / args.err_curve_out, dpi=150)
        plt.close()
        print(f"[Eval] Saved error curve: {args.err_curve_out}")

    return results, stats


# ============================================================
# 4)  CLI
# ============================================================

def build_argparser():
    p = argparse.ArgumentParser(description="Single-file FECAM-TCN-iTransformer: train + predict + eval")

    # paths
    p.add_argument("--base_dir", type=str,
                   default=r"D:\09-毕业论文\2025.07.08 小论文代码\神经网络缝模块\多模型一键运行",
                   help="你的数据目录（包含Input_GNSS.csv, Input_INS.csv, 出城INS.csv）")
    p.add_argument("--gnss_file", type=str, default="Input_GNSS.csv")
    p.add_argument("--ins_file", type=str, default="Input_INS.csv")
    p.add_argument("--test_file", type=str, default="出城INS.csv")

    # train range (GNSS)
    p.add_argument("--train_gnss_start", type=float, default=199900.0)
    p.add_argument("--train_gnss_end", type=float, default=200000.0)

    # test range
    p.add_argument("--test_start", type=float, default=129171.003)
    p.add_argument("--test_end", type=float, default=129181.003)

    # window
    p.add_argument("--window_seconds", type=float, default=3.0)

    # training hyperparams
    p.add_argument("--epochs", type=int, default=80)
    p.add_argument("--batch_size", type=int, default=64)
    p.add_argument("--lr", type=float, default=1e-3)
    p.add_argument("--weight_decay", type=float, default=1e-4)
    p.add_argument("--grad_clip", type=float, default=1.0)
    p.add_argument("--val_ratio", type=float, default=0.2)
    p.add_argument("--seed", type=int, default=42)
    p.add_argument("--log_every", type=int, default=10)

    # model hyperparams
    p.add_argument("--tcn_channels", type=int, default=64)
    p.add_argument("--tcn_kernel", type=int, default=3)
    p.add_argument("--tcn_layers_per_block", type=int, default=5)  # 5层 dilations=1,2,4,8,16 感受野足够覆盖3s@20Hz
    p.add_argument("--dropout_tcn", type=float, default=0.2)
    p.add_argument("--dropout_fecam", type=float, default=0.1)
    p.add_argument("--fecam_freq_bins", type=int, default=16)

    p.add_argument("--itr_d_model", type=int, default=128)
    p.add_argument("--itr_heads", type=int, default=8)
    p.add_argument("--itr_layers", type=int, default=2)
    p.add_argument("--dropout_itr", type=float, default=0.1)

    # outputs
    p.add_argument("--model_out", type=str, default="fecam_tcn_itransformer_model.pth")
    p.add_argument("--feature_scaler_out", type=str, default="fecam_feature_scaler.pkl")
    p.add_argument("--label_scaler_out", type=str, default="fecam_label_scaler.pkl")
    p.add_argument("--meta_out", type=str, default="fecam_meta.json")

    p.add_argument("--pred_out", type=str, default="predict_data_fecam_tcn_itransformer.csv")
    p.add_argument("--stats_out", type=str, default="error_statistics_fecam_tcn_itransformer.csv")

    p.add_argument("--train_curve_out", type=str, default="fecam_training_curve.png")
    p.add_argument("--err_curve_out", type=str, default="fecam_error_curve.png")

    # run control
    p.add_argument("--mode", type=str, default="all", choices=["all", "train", "eval"],
                   help="all=训练+评估；train=只训练；eval=只评估（需要已有模型+scaler+meta）")
    p.add_argument("--save_plots", action="store_true", help="保存训练曲线/误差曲线png")
    p.add_argument("--force_cpu", action="store_true", help="强制用CPU")
    return p


def main():
    args = build_argparser().parse_args()
    set_seed(args.seed)

    if args.mode in ("all", "train"):
        train_only(args)

    if args.mode in ("all", "eval"):
        eval_on_test(args)


if __name__ == "__main__":
    main()
'''
script_path.write_text(code, encoding="utf-8")
str(script_path)

'/mnt/data/fecam_tcn_itransformer_single.py'

python
import textwrap, os, json, re, math, random, pathlib
code = r"""
# -*- coding: utf-8 -*-
"""

python
from pathlib import Path
script_path = Path("/mnt/data/fecam_tcn_transformer_single.py")
script = r'''# -*- coding: utf-8 -*-
"""
FECAM-based TCN + Transformer (single file)
- Train: uses Input_GNSS.csv + Input_INS.csv (time range 199900.0-200000.0 by default)
- Input for each GNSS delta time: all 12 INS channels within previous 3 seconds (fixed window length inferred from INS dt)
- Output: GNSS position increments (delta_north, delta_east) in meters
- Eval: uses 出城INS.csv (or your test INS file) with lat/lon columns to compute coordinate error via Vincenty

Usage:
  python fecam_tcn_transformer_single.py --mode all --save_plots
  python fecam_tcn_transformer_single.py --mode train
  python fecam_tcn_transformer_single.py --mode eval

Notes:
- Column names are the same as your existing scripts.
- The script saves:
    fecam_tcn_transformer_model.pth
    fecam_feature_scaler.pkl
    fecam_label_scaler.pkl
    fecam_meta.json
    predict_data_fecam_tcn_transformer.csv
    error_statistics_fecam_tcn_transformer.csv
    (optional) plots
"""

import argparse
import json
import math
import os
from pathlib import Path
from typing import Dict, List, Optional, Tuple

import joblib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import torch
import torch.nn as nn
import torch.optim as optim
from sklearn.preprocessing import MinMaxScaler
from torch.nn.utils import weight_norm
from torch.utils.data import DataLoader, TensorDataset, random_split


# =========================
# Config / Defaults
# =========================
DEFAULT_BASE_DIR = r"D:\09-毕业论文\2025.07.08 小论文代码\神经网络缝模块\多模型一键运行"

TRAIN_GNSS_FILE = "Input_GNSS.csv"
TRAIN_INS_FILE = "Input_INS.csv"
TEST_INS_FILE = "出城INS.csv"

# Train range (same as your multi-model scripts)
DEFAULT_GNSS_START = 199900.0
DEFAULT_GNSS_END = 200000.0
DEFAULT_INS_START = 199900.008
DEFAULT_INS_END = 200000.008

# Evaluation range (same as your predict_all_models.py defaults)
DEFAULT_TEST_START = 129171.003
DEFAULT_TEST_END = 129181.003

FEATURE_COLS = [
    "北向速度(m/s)", "东向速度(m/s)", "地向速度(m/s)",
    "航向角(°)", "俯仰角(°)", "横滚角(°)",
    "加表x（g）", "加表y（g）", "加表z（g）",
    "陀螺x（°/s）", "陀螺y（°/s）", "陀螺z（°/s）",
]
INS_TIME_COL = "GPS周内秒(s)"
GNSS_TIME_COL = "GNSS周内秒(s)"
GNSS_LAT_COL = "GNSS纬度(°)"
GNSS_LON_COL = "GNSS经度(°)"

# Earth
EARTH_RADIUS = 6371000.0  # for approximate delta computation
WGS84_A = 6378137.0
WGS84_F = 1 / 298.257223563
WGS84_B = WGS84_A * (1 - WGS84_F)


# =========================
# Utilities
# =========================
def set_seed(seed: int = 42) -> None:
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)


def infer_dt(times: np.ndarray) -> float:
    """Infer sampling interval from time series (seconds). Robust to occasional jitter."""
    times = np.asarray(times, dtype=float)
    diffs = np.diff(times)
    diffs = diffs[np.isfinite(diffs)]
    diffs = diffs[diffs > 0]
    if len(diffs) == 0:
        raise ValueError("Cannot infer dt: time diffs are empty/non-positive.")
    return float(np.median(diffs))


def nearest_row_by_time(df: pd.DataFrame, time_col: str, t: float) -> pd.Series:
    """Return row with time closest to t (df must be sorted by time_col)."""
    arr = df[time_col].to_numpy(dtype=float)
    idx = np.searchsorted(arr, t)
    if idx <= 0:
        return df.iloc[0]
    if idx >= len(arr):
        return df.iloc[-1]
    # choose closer of idx-1 and idx
    if abs(arr[idx] - t) < abs(arr[idx - 1] - t):
        return df.iloc[idx]
    return df.iloc[idx - 1]


def calculate_gnss_position_deltas(gnss_df: pd.DataFrame) -> pd.DataFrame:
    """Compute north/east displacement between consecutive GNSS points (meters)."""
    df = gnss_df.sort_values(GNSS_TIME_COL).copy()
    # radians
    lat_rad = np.deg2rad(df[GNSS_LAT_COL].to_numpy(dtype=float))
    lon_rad = np.deg2rad(df[GNSS_LON_COL].to_numpy(dtype=float))
    dlat = np.diff(lat_rad)
    dlon = np.diff(lon_rad)
    # north approx
    dn = dlat * EARTH_RADIUS
    # east approx with cos(lat)
    lat_mid = lat_rad[1:]  # close enough for small step
    de = dlon * EARTH_RADIUS * np.cos(lat_mid)
    out = pd.DataFrame({
        GNSS_TIME_COL: df[GNSS_TIME_COL].iloc[1:].to_numpy(dtype=float),
        "delta_north": dn,
        "delta_east": de,
    })
    return out


def deltas_to_coords(lat_deg: float, lon_deg: float, delta_north_m: float, delta_east_m: float) -> Tuple[float, float]:
    """Approximate lat/lon update from north/east displacements."""
    lat = float(lat_deg)
    lon = float(lon_deg)
    dlat = (delta_north_m / EARTH_RADIUS) * (180.0 / math.pi)
    dlon = (delta_east_m / (EARTH_RADIUS * math.cos(lat * math.pi / 180.0))) * (180.0 / math.pi)
    return lat + dlat, lon + dlon


def vincenty_offset_north_east(lat1_deg: float, lon1_deg: float, lat2_deg: float, lon2_deg: float) -> Tuple[float, float]:
    """
    Vincenty inverse: returns (north_offset_m, east_offset_m) from point1 -> point2 (WGS84).
    """
    a = WGS84_A
    f = WGS84_F
    b = WGS84_B

    phi1, lam1, phi2, lam2 = map(math.radians, [lat1_deg, lon1_deg, lat2_deg, lon2_deg])
    U1 = math.atan((1 - f) * math.tan(phi1))
    U2 = math.atan((1 - f) * math.tan(phi2))
    L = lam2 - lam1
    lam = L

    sinU1, cosU1 = math.sin(U1), math.cos(U1)
    sinU2, cosU2 = math.sin(U2), math.cos(U2)

    for _ in range(200):
        sin_lam = math.sin(lam)
        cos_lam = math.cos(lam)
        sin_sigma = math.sqrt((cosU2 * sin_lam) ** 2 + (cosU1 * sinU2 - sinU1 * cosU2 * cos_lam) ** 2)
        if sin_sigma == 0:
            return 0.0, 0.0
        cos_sigma = sinU1 * sinU2 + cosU1 * cosU2 * cos_lam
        sigma = math.atan2(sin_sigma, cos_sigma)
        sin_alpha = cosU1 * cosU2 * sin_lam / sin_sigma
        cos2_alpha = 1 - sin_alpha ** 2
        if cos2_alpha == 0:
            cos2_sigma_m = 0.0
        else:
            cos2_sigma_m = cos_sigma - 2 * sinU1 * sinU2 / cos2_alpha
        C = f / 16 * cos2_alpha * (4 + f * (4 - 3 * cos2_alpha))
        lam_new = L + (1 - C) * f * sin_alpha * (
            sigma + C * sin_sigma * (cos2_sigma_m + C * cos_sigma * (-1 + 2 * cos2_sigma_m ** 2))
        )
        if abs(lam_new - lam) < 1e-12:
            lam = lam_new
            break
        lam = lam_new

    # azimuth from point1 to point2
    azimuth = math.atan2(
        cosU2 * math.sin(lam),
        cosU1 * sinU2 - sinU1 * cosU2 * math.cos(lam),
    )
    distance = b * sigma  # adequate for our usage here

    north = distance * math.cos(azimuth)
    east = distance * math.sin(azimuth)
    return float(north), float(east)


def find_lat_lon_cols(df: pd.DataFrame) -> Tuple[str, str]:
    """Find lat/lon columns in test INS file."""
    candidates = [
        ("纬度(°)", "经度(°)"),
        ("GNSS纬度(°)", "GNSS经度(°)"),
        ("Latitude", "Longitude"),
        ("lat", "lon"),
    ]
    for lat_c, lon_c in candidates:
        if lat_c in df.columns and lon_c in df.columns:
            return lat_c, lon_c
    raise KeyError(f"Cannot find lat/lon columns in test file. Tried: {candidates}")


# =========================
# Model: FECAM + TCN + Transformer
# =========================
class Chomp1d(nn.Module):
    """Chomp off padding on the right to keep causality."""
    def __init__(self, chomp_size: int):
        super().__init__()
        self.chomp_size = chomp_size

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        if self.chomp_size == 0:
            return x
        return x[:, :, :-self.chomp_size].contiguous()


class TemporalBlock(nn.Module):
    """TCN residual block with dilated causal conv, weight norm, dropout."""
    def __init__(self, in_ch: int, out_ch: int, kernel_size: int, dilation: int, dropout: float):
        super().__init__()
        padding = (kernel_size - 1) * dilation

        self.conv1 = weight_norm(nn.Conv1d(in_ch, out_ch, kernel_size, padding=padding, dilation=dilation))
        self.chomp1 = Chomp1d(padding)
        self.relu1 = nn.ReLU()
        self.drop1 = nn.Dropout(dropout)

        self.conv2 = weight_norm(nn.Conv1d(out_ch, out_ch, kernel_size, padding=padding, dilation=dilation))
        self.chomp2 = Chomp1d(padding)
        self.relu2 = nn.ReLU()
        self.drop2 = nn.Dropout(dropout)

        self.downsample = nn.Conv1d(in_ch, out_ch, 1) if in_ch != out_ch else nn.Identity()
        self.final_relu = nn.ReLU()

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        out = self.conv1(x)
        out = self.chomp1(out)
        out = self.relu1(out)
        out = self.drop1(out)

        out = self.conv2(out)
        out = self.chomp2(out)
        out = self.relu2(out)
        out = self.drop2(out)

        res = self.downsample(x)
        return self.final_relu(out + res)


class DCTCache:
    """Cache cosine basis for DCT-II for a given (T, K)."""
    def __init__(self):
        self._cache: Dict[Tuple[int, int, torch.device, torch.dtype], torch.Tensor] = {}

    def get_basis(self, T: int, K: int, device: torch.device, dtype: torch.dtype) -> torch.Tensor:
        key = (T, K, device, dtype)
        if key in self._cache:
            return self._cache[key]
        n = torch.arange(T, device=device, dtype=dtype).view(1, 1, T)  # (1,1,T)
        k = torch.arange(K, device=device, dtype=dtype).view(1, K, 1)  # (1,K,1)
        # DCT-II basis: cos(pi/T*(n+0.5)*k)
        basis = torch.cos(math.pi / T * (n + 0.5) * k)  # (1,K,T)
        self._cache[key] = basis
        return basis


_DCT_CACHE = DCTCache()


class FECAM(nn.Module):
    """
    Frequency-Enhanced Channel Attention (practical implementation):
    - DCT along time dimension to get low-frequency descriptors
    - channel attention weights from frequency descriptors
    - small 1x1 Conv "FCNN" to generate enhancement tensor
    """
    def __init__(self, channels: int, dct_k: int = 16, reduction: int = 8, dropout: float = 0.1):
        super().__init__()
        self.channels = channels
        self.dct_k = dct_k
        hidden = max(4, channels // reduction)

        # Channel attention from frequency features (B,C) -> (B,C)
        self.ca = nn.Sequential(
            nn.Linear(channels, hidden),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(hidden, channels),
            nn.Sigmoid(),
        )

        # FCNN-like tensor generator (1x1 convs)
        self.fcnn = nn.Sequential(
            nn.Conv1d(channels, channels, kernel_size=1),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Conv1d(channels, channels, kernel_size=1),
            nn.Sigmoid(),
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        """
        x: (B,C,T)
        """
        B, C, T = x.shape
        K = min(self.dct_k, T)

        basis = _DCT_CACHE.get_basis(T, K, x.device, x.dtype)  # (1,K,T)
        # DCT per channel: (B,C,T) -> (B,C,K) via matmul with basis
        # Expand basis for broadcasting
        # x_unsq: (B,C,1,T); basis: (1,1,K,T) -> multiply & sum over T
        dct = (x.unsqueeze(2) * basis.unsqueeze(1)).sum(dim=-1)  # (B,C,K)
        freq_desc = dct.abs().mean(dim=-1)  # (B,C)

        w = self.ca(freq_desc).unsqueeze(-1)  # (B,C,1)
        x_scaled = x * w

        enh = self.fcnn(x_scaled)  # (B,C,T) in [0,1]
        return x_scaled * (1.0 + enh)  # residual enhancement


class PositionalEncoding(nn.Module):
    """Standard sinusoidal positional encoding for Transformer."""
    def __init__(self, d_model: int, max_len: int = 5000):
        super().__init__()
        pe = torch.zeros(max_len, d_model, dtype=torch.float32)
        position = torch.arange(0, max_len, dtype=torch.float32).unsqueeze(1)
        div_term = torch.exp(torch.arange(0, d_model, 2, dtype=torch.float32) * (-math.log(10000.0) / d_model))
        pe[:, 0::2] = torch.sin(position * div_term)
        pe[:, 1::2] = torch.cos(position * div_term)
        self.register_buffer("pe", pe.unsqueeze(0), persistent=False)  # (1,max_len,d_model)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        """
        x: (B,T,D)
        """
        T = x.size(1)
        if T > self.pe.size(1):
            raise ValueError(f"Sequence length {T} exceeds max_len {self.pe.size(1)} for positional encoding.")
        return x + self.pe[:, :T, :].to(x.dtype)


class FECAMTCNTransformerModel(nn.Module):
    """
    FECAM-based TCN front-end + standard Transformer encoder (temporal tokens).
    Input: (B,T,input_dim)
    Output: (B,2)
    """
    def __init__(
        self,
        input_dim: int,
        seq_length: int,
        tcn_channels: int = 64,
        kernel_size: int = 3,
        dilation1: int = 1,
        dilation2: int = 2,
        fecam_dct_k: int = 16,
        fecam_dropout: float = 0.1,
        tcn_dropout: float = 0.2,
        d_model: int = 64,
        nhead: int = 8,
        num_layers: int = 2,
        ff_dim: int = 256,
        tr_dropout: float = 0.1,
    ):
        super().__init__()
        self.seq_length = seq_length

        # Input projection for TCN
        self.in_proj = nn.Conv1d(input_dim, tcn_channels, kernel_size=1)

        # Two TCN residual blocks with inserted FECAM
        self.tcn1 = TemporalBlock(tcn_channels, tcn_channels, kernel_size, dilation1, tcn_dropout)
        self.fecam = FECAM(tcn_channels, dct_k=fecam_dct_k, dropout=fecam_dropout)
        self.tcn2 = TemporalBlock(tcn_channels, tcn_channels, kernel_size, dilation2, tcn_dropout)

        # Transformer encoder (temporal tokens)
        self.to_d_model = nn.Linear(tcn_channels, d_model)
        self.pos_enc = PositionalEncoding(d_model=d_model, max_len=max(5000, seq_length + 10))
        encoder_layer = nn.TransformerEncoderLayer(
            d_model=d_model,
            nhead=nhead,
            dim_feedforward=ff_dim,
            dropout=tr_dropout,
            batch_first=True,
            activation="gelu",
            norm_first=True,
        )
        self.transformer = nn.TransformerEncoder(encoder_layer, num_layers=num_layers)

        # Output head (use last token)
        self.head = nn.Sequential(
            nn.LayerNorm(d_model),
            nn.Linear(d_model, 128),
            nn.GELU(),
            nn.Dropout(tr_dropout),
            nn.Linear(128, 2),
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        """
        x: (B,T,input_dim)
        """
        # TCN expects (B,C,T)
        x = x.transpose(1, 2)  # (B,input_dim,T)
        x = self.in_proj(x)    # (B,tcn_channels,T)
        x = self.tcn1(x)
        x = self.fecam(x)
        x = self.tcn2(x)

        # Transformer expects (B,T,D)
        x = x.transpose(1, 2)          # (B,T,tcn_channels)
        x = self.to_d_model(x)         # (B,T,d_model)
        x = self.pos_enc(x)
        x = self.transformer(x)        # (B,T,d_model)
        last = x[:, -1, :]             # (B,d_model)
        return self.head(last)


# =========================
# Data building
# =========================
def build_train_samples(
    gnss_df: pd.DataFrame,
    ins_df: pd.DataFrame,
    window_seconds: float,
) -> Tuple[np.ndarray, np.ndarray, Dict]:
    """
    Create (X, y) for training:
      - y: GNSS displacement between consecutive GNSS points at time t_k
      - X: INS window within [t_k - window_seconds, t_k], sampled at INS rate (fixed length)
    """
    gnss_df = gnss_df.sort_values(GNSS_TIME_COL).reset_index(drop=True)
    ins_df = ins_df.sort_values(INS_TIME_COL).reset_index(drop=True)

    # infer dt for INS
    ins_dt = infer_dt(ins_df[INS_TIME_COL].to_numpy())
    window_len = int(round(window_seconds / ins_dt)) + 1

    # infer GNSS step (for meta only)
    gnss_times = gnss_df[GNSS_TIME_COL].to_numpy(dtype=float)
    gnss_diffs = np.diff(gnss_times)
    gnss_diffs = gnss_diffs[np.isfinite(gnss_diffs) & (gnss_diffs > 0)]
    gnss_step = float(np.median(gnss_diffs)) if len(gnss_diffs) else 1.0

    # estimate time offset between INS and GNSS within the filtered window
    # (matches your current filtering: INS range is GNSS range + ~0.008)
    offset = float(ins_df[INS_TIME_COL].iloc[0] - gnss_df[GNSS_TIME_COL].iloc[0])

    deltas_df = calculate_gnss_position_deltas(gnss_df)  # has GNSS_TIME_COL + delta_north/east at t_k

    ins_times = ins_df[INS_TIME_COL].to_numpy(dtype=float)
    ins_feat = ins_df[FEATURE_COLS].to_numpy(dtype=float)

    X_list: List[np.ndarray] = []
    y_list: List[np.ndarray] = []
    used_times: List[float] = []

    for _, row in deltas_df.iterrows():
        t_gnss = float(row[GNSS_TIME_COL])
        t_end = t_gnss + offset
        # take all INS samples up to t_end, then last window_len
        end_idx = np.searchsorted(ins_times, t_end, side="right") - 1
        if end_idx < 0:
            continue
        start_idx = end_idx - window_len + 1
        if start_idx < 0:
            continue

        window = ins_feat[start_idx:end_idx + 1]  # (window_len, 12)
        if window.shape[0] != window_len:
            continue

        X_list.append(window.astype(np.float32))
        y_list.append(np.array([row["delta_north"], row["delta_east"]], dtype=np.float32))
        used_times.append(t_gnss)

    if len(X_list) < 10:
        raise RuntimeError(f"Too few training samples built: {len(X_list)}. Check time alignment / ranges.")

    meta = {
        "window_seconds": window_seconds,
        "ins_dt": ins_dt,
        "window_len": window_len,
        "gnss_step_seconds": gnss_step,
        "ins_gnss_offset_seconds": offset,
        "num_samples": len(X_list),
        "label_time_definition": "delta between consecutive GNSS points, timestamp at later GNSS time",
        "features": FEATURE_COLS,
    }
    return np.stack(X_list, axis=0), np.stack(y_list, axis=0), meta


def build_eval_samples(
    ins_df: pd.DataFrame,
    window_len: int,
    window_seconds: float,
    step_seconds: float,
    time_col: str,
    feature_cols: List[str],
    lat_col: str,
    lon_col: str,
) -> Tuple[np.ndarray, List[float], List[Tuple[float, float]], List[Tuple[float, float]]]:
    """
    Build eval windows from a test INS file (with lat/lon ground truth):
      - For each step end time t_end, predict delta from t0=t_end-step_seconds to t_end
      - window ends at t_end, length=window_len
      - base_coords: (lat/lon at t0), true_coords: (lat/lon at t_end)
    Returns:
      X (N, window_len, 12), end_times, base_coords, true_coords
    """
    ins_df = ins_df.sort_values(time_col).reset_index(drop=True)
    times = ins_df[time_col].to_numpy(dtype=float)
    feats = ins_df[feature_cols].to_numpy(dtype=float)

    # Choose evaluation times: start after enough history and after 1 step
    t_min = float(times[0])
    t_max = float(times[-1])

    # Ensure we have at least window history before first prediction
    first_end = t_min + window_seconds
    first_end = max(first_end, t_min + step_seconds)

    end_times = np.arange(first_end, t_max + 1e-9, step_seconds, dtype=float)
    if len(end_times) < 3:
        raise RuntimeError("Too few evaluation points. Check your TEST_START/END or step_seconds.")

    X_list: List[np.ndarray] = []
    used_end_times: List[float] = []
    base_coords: List[Tuple[float, float]] = []
    true_coords: List[Tuple[float, float]] = []

    for t_end in end_times:
        # window: last window_len samples <= t_end
        end_idx = np.searchsorted(times, t_end, side="right") - 1
        if end_idx < 0:
            continue
        start_idx = end_idx - window_len + 1
        if start_idx < 0:
            continue

        window = feats[start_idx:end_idx + 1]
        if window.shape[0] != window_len:
            continue

        t0 = t_end - step_seconds
        base_row = nearest_row_by_time(ins_df, time_col, t0)
        true_row = nearest_row_by_time(ins_df, time_col, t_end)

        base_coords.append((float(base_row[lat_col]), float(base_row[lon_col])))
        true_coords.append((float(true_row[lat_col]), float(true_row[lon_col])))

        X_list.append(window.astype(np.float32))
        used_end_times.append(float(t_end))

    if len(X_list) < 3:
        raise RuntimeError("Could not build evaluation windows. Check window_len/step_seconds and time coverage.")

    return np.stack(X_list, axis=0), used_end_times, base_coords, true_coords


# =========================
# Train / Eval
# =========================
def train(
    base_dir: Path,
    window_seconds: float,
    epochs: int,
    batch_size: int,
    lr: float,
    weight_decay: float,
    device: torch.device,
    save_plots: bool,
) -> Dict:
    gnss_path = base_dir / TRAIN_GNSS_FILE
    ins_path = base_dir / TRAIN_INS_FILE
    if not gnss_path.exists() or not ins_path.exists():
        raise FileNotFoundError(f"Training files not found: {gnss_path} / {ins_path}")

    gnss = pd.read_csv(gnss_path)
    ins = pd.read_csv(ins_path)

    # Filter time range (same as your old scripts)
    gnss = gnss[(gnss[GNSS_TIME_COL] >= DEFAULT_GNSS_START) & (gnss[GNSS_TIME_COL] <= DEFAULT_GNSS_END)].copy()
    ins = ins[(ins[INS_TIME_COL] >= DEFAULT_INS_START) & (ins[INS_TIME_COL] <= DEFAULT_INS_END)].copy()

    X, y, meta = build_train_samples(gnss, ins, window_seconds=window_seconds)
    seq_len = int(meta["window_len"])

    # Normalize
    feat_scaler = MinMaxScaler()
    X_2d = X.reshape(-1, X.shape[-1])
    Xn = feat_scaler.fit_transform(X_2d).reshape(X.shape)

    label_scaler = MinMaxScaler()
    yn = label_scaler.fit_transform(y)

    joblib.dump(feat_scaler, base_dir / "fecam_feature_scaler.pkl")
    joblib.dump(label_scaler, base_dir / "fecam_label_scaler.pkl")

    # Save meta
    meta_path = base_dir / "fecam_meta.json"
    with open(meta_path, "w", encoding="utf-8") as f:
        json.dump(meta, f, ensure_ascii=False, indent=2)

    # Dataset
    X_tensor = torch.tensor(Xn, dtype=torch.float32)
    y_tensor = torch.tensor(yn, dtype=torch.float32)

    dataset = TensorDataset(X_tensor, y_tensor)
    train_size = int(0.8 * len(dataset))
    val_size = len(dataset) - train_size
    train_ds, val_ds = random_split(dataset, [train_size, val_size], generator=torch.Generator().manual_seed(42))

    train_loader = DataLoader(train_ds, batch_size=batch_size, shuffle=True, drop_last=False)
    val_loader = DataLoader(val_ds, batch_size=batch_size, shuffle=False, drop_last=False)

    # Model
    model = FECAMTCNTransformerModel(
        input_dim=len(FEATURE_COLS),
        seq_length=seq_len,
        tcn_channels=64,
        kernel_size=3,
        dilation1=1,
        dilation2=2,
        fecam_dct_k=16,
        fecam_dropout=0.1,
        tcn_dropout=0.2,
        d_model=64,
        nhead=8,
        num_layers=2,
        ff_dim=256,
        tr_dropout=0.1,
    ).to(device)

    criterion = nn.MSELoss()
    optimizer = optim.AdamW(model.parameters(), lr=lr, weight_decay=weight_decay)

    train_losses: List[float] = []
    val_losses: List[float] = []

    best_val = float("inf")
    best_path = base_dir / "fecam_tcn_transformer_model.pth"

    for ep in range(1, epochs + 1):
        model.train()
        total = 0.0
        for xb, yb in train_loader:
            xb = xb.to(device)
            yb = yb.to(device)

            optimizer.zero_grad(set_to_none=True)
            pred = model(xb)
            loss = criterion(pred, yb)
            loss.backward()
            torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)
            optimizer.step()
            total += float(loss.item())
        train_loss = total / max(1, len(train_loader))
        train_losses.append(train_loss)

        model.eval()
        total = 0.0
        with torch.no_grad():
            for xb, yb in val_loader:
                xb = xb.to(device)
                yb = yb.to(device)
                pred = model(xb)
                loss = criterion(pred, yb)
                total += float(loss.item())
        val_loss = total / max(1, len(val_loader))
        val_losses.append(val_loss)

        if val_loss < best_val:
            best_val = val_loss
            torch.save(model.state_dict(), best_path)

        if ep % 10 == 0 or ep == 1 or ep == epochs:
            print(f"[Train] Epoch {ep}/{epochs}  train={train_loss:.6f}  val={val_loss:.6f}  best_val={best_val:.6f}")

    if save_plots:
        plt.figure(figsize=(8, 5))
        plt.plot(range(1, epochs + 1), train_losses, label="Train Loss")
        plt.plot(range(1, epochs + 1), val_losses, label="Val Loss", linestyle="--")
        plt.xlabel("Epoch")
        plt.ylabel("MSE Loss")
        plt.title("FECAM-TCN-Transformer Training Curve")
        plt.grid(True, alpha=0.3)
        plt.legend()
        plt.tight_layout()
        plt.savefig(base_dir / "fecam_tcn_transformer_training_curve.png", dpi=200)
        plt.close()

    out = {
        "meta": meta,
        "best_model_path": str(best_path),
        "train_losses": train_losses,
        "val_losses": val_losses,
    }
    return out


def eval_and_report(
    base_dir: Path,
    test_file: str,
    test_start: float,
    test_end: float,
    device: torch.device,
    save_plots: bool,
) -> Dict:
    # Load meta & scalers & model
    meta_path = base_dir / "fecam_meta.json"
    if not meta_path.exists():
        raise FileNotFoundError(f"Missing meta file: {meta_path} (train first).")
    meta = json.loads(meta_path.read_text(encoding="utf-8"))

    feat_scaler = joblib.load(base_dir / "fecam_feature_scaler.pkl")
    label_scaler = joblib.load(base_dir / "fecam_label_scaler.pkl")

    model_path = base_dir / "fecam_tcn_transformer_model.pth"
    if not model_path.exists():
        raise FileNotFoundError(f"Missing model: {model_path} (train first).")

    window_len = int(meta["window_len"])
    window_seconds = float(meta["window_seconds"])
    step_seconds = float(meta.get("gnss_step_seconds", 1.0))

    # Load and filter test file
    test_path = base_dir / test_file
    if not test_path.exists():
        raise FileNotFoundError(f"Test file not found: {test_path}")
    test_df = pd.read_csv(test_path)
    if INS_TIME_COL not in test_df.columns:
        raise KeyError(f"Test file missing time col '{INS_TIME_COL}'")

    test_df = test_df[(test_df[INS_TIME_COL] >= test_start) & (test_df[INS_TIME_COL] <= test_end)].copy()
    test_df = test_df.sort_values(INS_TIME_COL).reset_index(drop=True)
    if len(test_df) < window_len + 5:
        raise RuntimeError(f"Test data too short after filtering: {len(test_df)} rows.")

    lat_col, lon_col = find_lat_lon_cols(test_df)

    # Build eval windows
    X, end_times, base_coords, true_coords = build_eval_samples(
        ins_df=test_df,
        window_len=window_len,
        window_seconds=window_seconds,
        step_seconds=step_seconds,
        time_col=INS_TIME_COL,
        feature_cols=FEATURE_COLS,
        lat_col=lat_col,
        lon_col=lon_col,
    )

    # Normalize features
    Xn = feat_scaler.transform(X.reshape(-1, X.shape[-1])).reshape(X.shape)
    X_tensor = torch.tensor(Xn, dtype=torch.float32).to(device)

    # Load model
    model = FECAMTCNTransformerModel(
        input_dim=len(FEATURE_COLS),
        seq_length=window_len,
        tcn_channels=64,
        kernel_size=3,
        dilation1=1,
        dilation2=2,
        fecam_dct_k=16,
        fecam_dropout=0.1,
        tcn_dropout=0.2,
        d_model=64,
        nhead=8,
        num_layers=2,
        ff_dim=256,
        tr_dropout=0.1,
    ).to(device)
    model.load_state_dict(torch.load(model_path, map_location=device))
    model.eval()

    with torch.no_grad():
        pred_norm = model(X_tensor).cpu().numpy()
    pred = label_scaler.inverse_transform(pred_norm)  # (N,2) in meters

    # Generate predicted coords + errors
    pred_coords: List[Tuple[float, float]] = []
    err_n: List[float] = []
    err_e: List[float] = []
    err_total: List[float] = []

    true_dn: List[float] = []
    true_de: List[float] = []

    for (lat0, lon0), (lat_true, lon_true), (dn_pred, de_pred) in zip(base_coords, true_coords, pred):
        lat_pred, lon_pred = deltas_to_coords(lat0, lon0, float(dn_pred), float(de_pred))
        pred_coords.append((lat_pred, lon_pred))

        # error: true -> pred
        n, e = vincenty_offset_north_east(lat_true, lon_true, lat_pred, lon_pred)
        err_n.append(n)
        err_e.append(e)
        err_total.append(math.sqrt(n * n + e * e))

        # (optional) true delta between base and true (for reference)
        dn_t, de_t = vincenty_offset_north_east(lat0, lon0, lat_true, lon_true)
        true_dn.append(dn_t)
        true_de.append(de_t)

    results = pd.DataFrame({
        "时间戳": end_times,
        "基准纬度": [c[0] for c in base_coords],
        "基准经度": [c[1] for c in base_coords],
        "真值纬度": [c[0] for c in true_coords],
        "真值经度": [c[1] for c in true_coords],
        "北向增量真值(m)": true_dn,
        "东向增量真值(m)": true_de,
        "北向增量预测(m)": pred[:, 0],
        "东向增量预测(m)": pred[:, 1],
        "预测纬度": [c[0] for c in pred_coords],
        "预测经度": [c[1] for c in pred_coords],
        "北向误差(m)": err_n,
        "东向误差(m)": err_e,
        "总误差(m)": err_total,
    })

    pred_csv = base_dir / "predict_data_fecam_tcn_transformer.csv"
    results.to_csv(pred_csv, index=False, encoding="utf-8-sig")

    # Stats
    min_err = float(np.min(err_total))
    max_err = float(np.max(err_total))
    mean_err = float(np.mean(err_total))
    std_err = float(np.std(err_total, ddof=1)) if len(err_total) > 1 else 0.0
    rmse = float(np.sqrt(np.mean(np.square(err_total))))

    stats = pd.DataFrame([{
        "模型": "FECAM_TCN_Transformer",
        "最小误差(m)": min_err,
        "最大误差(m)": max_err,
        "平均误差(m)": mean_err,
        "标准差(m)": std_err,
        "RMSE(m)": rmse,
        "step_seconds": step_seconds,
        "window_seconds": window_seconds,
        "window_len": window_len,
    }])
    stats_csv = base_dir / "error_statistics_fecam_tcn_transformer.csv"
    stats.to_csv(stats_csv, index=False, encoding="utf-8-sig")

    print("=" * 60)
    print("FECAM_TCN_Transformer 误差统计")
    print("=" * 60)
    print(f"最小误差: {min_err:.4f} m")
    print(f"最大误差: {max_err:.4f} m")
    print(f"平均误差: {mean_err:.4f} m")
    print(f"标准差:   {std_err:.4f} m")
    print(f"RMSE:     {rmse:.4f} m")
    print(f"\n预测结果已保存: {pred_csv}")
    print(f"误差统计已保存: {stats_csv}")

    if save_plots:
        plt.figure(figsize=(10, 4))
        plt.plot(results.index, results["总误差(m)"].to_numpy(), linewidth=1.2)
        plt.xlabel("时间步")
        plt.ylabel("总误差 (m)")
        plt.title("FECAM-TCN-Transformer 总误差曲线")
        plt.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.savefig(base_dir / "fecam_tcn_transformer_error_curve.png", dpi=200)
        plt.close()

    return {
        "predict_csv": str(pred_csv),
        "stats_csv": str(stats_csv),
        "min_err": min_err,
        "max_err": max_err,
        "mean_err": mean_err,
        "rmse": rmse,
        "window_len": window_len,
        "step_seconds": step_seconds,
    }


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--base_dir", type=str, default=DEFAULT_BASE_DIR, help="Directory containing csv files")
    parser.add_argument("--mode", type=str, default="all", choices=["train", "eval", "all"])
    parser.add_argument("--window_seconds", type=float, default=3.0, help="INS history length (seconds)")
    parser.add_argument("--epochs", type=int, default=80)
    parser.add_argument("--batch_size", type=int, default=64)
    parser.add_argument("--lr", type=float, default=1e-3)
    parser.add_argument("--weight_decay", type=float, default=1e-4)
    parser.add_argument("--test_file", type=str, default=TEST_INS_FILE)
    parser.add_argument("--test_start", type=float, default=DEFAULT_TEST_START)
    parser.add_argument("--test_end", type=float, default=DEFAULT_TEST_END)
    parser.add_argument("--save_plots", action="store_true")
    parser.add_argument("--seed", type=int, default=42)
    args = parser.parse_args()

    set_seed(args.seed)
    base_dir = Path(args.base_dir)
    base_dir.mkdir(parents=True, exist_ok=True)

    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
    print(f"Using device: {device}")

    if args.mode in ("train", "all"):
        train(
            base_dir=base_dir,
            window_seconds=args.window_seconds,
            epochs=args.epochs,
            batch_size=args.batch_size,
            lr=args.lr,
            weight_decay=args.weight_decay,
            device=device,
            save_plots=args.save_plots,
        )

    if args.mode in ("eval", "all"):
        eval_and_report(
            base_dir=base_dir,
            test_file=args.test_file,
            test_start=args.test_start,
            test_end=args.test_end,
            device=device,
            save_plots=args.save_plots,
        )


if __name__ == "__main__":
    main()
'''
script_path.write_text(script, encoding="utf-8")
str(script_path), script_path.stat().st_size


('/mnt/data/fecam_tcn_transformer_single.py', 33097)