Fri Oct 17 2025 00:00:00 GMT+0800 (中國標準時間)

相当于现在的三个模态的数据：

1. 高光谱：hdr，spe的数据
2. 原始图片：.bmp
3. 荧光值：表格中的item[Fluorescence intensity(A.U.)]

输出：Area/Toxic content（Area公式可以计算得到Toxic content）

想要的结果分析：

除了这个训练得到的准确度之外，

最好还可以分析有关于不同的波长段对应的信息，哪个对于最终的回归结果更有影响（权重更大）

以及以单图像进行输入，锁定目标位置后进行训练，得到的准确度的对比。

甲方要求：模型做一下优化，图像信息采集上能否也做优化

我在 multimodalDataset.py 中添加了 prefilter_valid()（找出 image+hdr+spe 都存在的样本）和 get_splits(val_ratio, test_ratio)（返回 train/val/test 索引）；

• 在本地运行结果：
  • 样本总数（Excel 中）= 1329
  • prefilter_valid 找到可用三模态样本数 = 1329（说明 _find_file 现在能为所有样本定位到文件）
  • 按默认 80/10/10 划分得到 sizes: 1065 / 132 / 132
• 取了一个 DataLoader batch（batch_size=4）并用 collate_fn_mosi_regression 验证：
  • image: torch.Size([4, 3, 960, 960])
  • hyperspectral: torch.Size([4, 960, 960, 600])
  • fluorescence: torch.Size([4, 1])
  • labels: [4], mcs: [4], observeds: [4,3]

我现在需要利用一个具有三个模态的数据，来在I2MoE方法上进行训练测试。第一步是对数据进行编写相关的dataset：你现在先看一下这个表格/Volumes/Extreme Pro/001-实验数据-是这个/001-是这个/20251015-数据统计-是这个.xlsx中的数据，其中有荧光率的数据，以及对应的样品的name，荧光值：表格中的item[Fluorescence intensity(A.U.)]，表格中还有Area还有Toxic content这两个item，这两个item之间存在对应的关系，可以根据Area计算出Toxic content。所以我选定Area作为最终回归的任务的评估的ground truth（label）。此外，在/Volumes/Extreme Pro/001-实验数据-是这个/001-是这个/001-高光谱-是这个/4H/4-1/或事/Volumes/Extreme Pro/001-实验数据-是这个/001-是这个/001-高光谱-是这个/14h-2/14-33这种文件夹具有.bmp文件数据对应原始图片数据，以及两个.hdr和.spe文件对应原始测量 + 参考校准的高光谱和反射率的数据，我现在需要写出一个dataset能够对这三种数据进行处理，然后在getitem返回对应数据以及GT，然后便于后续的训练。

我现在需要你针对这个excel文件中Sample NameFluorescence intensity(A.U.)AreaToxic content的这个item，还有得到的/Volumes/Extreme Pro/001-实验数据-是这个/001-是这个/test/processed_npy/new_metadata.csv中对应的filename，把Fluorescence intensity(A.U.)这个item提取出来写一个csv，这个csv中要有filename，Fluorescence intensity(A.U.)，以及Sample Name。都要对应。

import os
import json
import math
import numpy as np
import pandas as pd
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torch.utils.data import DataLoader, Dataset, Subset
from torchvision import models, transforms
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from sklearn.metrics import mean_absolute_error, r2_score, mean_squared_error
import matplotlib.pyplot as plt
import seaborn as sns
from PIL import Image
from tqdm import tqdm

sns.set_theme(style="whitegrid", font_scale=0.9)
plt.rcParams['font.sans-serif'] = ['DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

def get_best_device():
    """Select the best available device: CUDA > MPS (Apple Silicon) > CPU."""
    try:
        if torch.cuda.is_available():
            return 'cuda'
    except Exception:
        pass
    try:
        if hasattr(torch.backends, 'mps') and torch.backends.mps.is_available():
            return 'mps'
    except Exception:
        pass
    return 'cpu'

class HyperspectralDataset(Dataset):
    def __init__(self, root_folder, metadata_file='new_metadata.csv', transform=None, target_size=(256, 256), scaler=None, feature_scaler=None, target_channels=None):
        self.root_folder = os.path.normpath(root_folder)
        self.transform = transform
        self.target_size = target_size
        self.target_channels = target_channels
        self.file_list = []
        self.labels = []
        self.raw_values = []

        metadata_path = os.path.join(self.root_folder, metadata_file)
        if not os.path.exists(metadata_path):
            raise FileNotFoundError(f"Metadata file {metadata_path} not found")

        metadata = pd.read_csv(
            metadata_path,
            dtype={'toxin_value': float, 'Value_Raw': float},
            encoding='utf-8-sig'
        )
        
        # Filter out abnormal data with toxin_value > 20
        valid_metadata = metadata[metadata['toxin_value'] <= 20]
        print(f"Found {len(valid_metadata)} valid samples (toxin_value ≤ 20)")

        # Count filtered abnormal data
        invalid_count = len(metadata) - len(valid_metadata)
        if invalid_count > 0:
            print(f"⚠️ Warning: Filtered {invalid_count} abnormal samples (toxin_value > 20)")
        
        # Process valid data
        # Prefer using npy_filename (相对于 root_folder 的 .npy 路径)
        has_npy_col = 'npy_filename' in valid_metadata.columns
        for idx, row in valid_metadata.iterrows():
            try:
                # choose path source
                if has_npy_col and isinstance(row['npy_filename'], str) and len(row['npy_filename'].strip()) > 0:
                    rel_path = row['npy_filename'].strip()
                else:
                    # fallback: use 'filename' only if it points to a .npy (relative or absolute)
                    rel_path = str(row['filename']).strip()
                # resolve to absolute path
                if os.path.isabs(rel_path):
                    file_path = os.path.normpath(rel_path)
                else:
                    file_path = os.path.normpath(os.path.join(self.root_folder, rel_path))
                if not os.path.isfile(file_path):
                    print(f"Warning: File {file_path} not found, skipping")
                    continue
                toxin_value = float(row['toxin_value'])
                raw_value = float(row['Value_Raw'])
                
                self.file_list.append(file_path)
                self.labels.append(toxin_value)
                self.raw_values.append(raw_value)
            except Exception as e:
                print(f"Error processing file {file_path}: {e}")
                continue

        if not self.file_list:
            raise RuntimeError("No valid data found")

        print(f"Successfully loaded {len(self.file_list)} samples")
        self.labels = np.array(self.labels).reshape(-1, 1)
        self.raw_values = np.array(self.raw_values).reshape(-1, 1)
        
        # Normalize labels to [0, 1] range using MinMaxScaler
        if scaler:
            self.labels = scaler.transform(self.labels)
        self.labels = self.labels.flatten()
        
        # Normalize features using StandardScaler
        if feature_scaler:
            self.raw_values = feature_scaler.transform(self.raw_values)
        self.raw_values = self.raw_values.flatten()

        # Check sample data shape
        if self.file_list:
            sample_data = np.load(self.file_list[0])
            if sample_data.shape[:2] != self.target_size:
                print(f"⚠️ Warning: Input size {sample_data.shape[:2]} does not match target size {self.target_size}, resizing automatically")

    def __len__(self):
        return len(self.file_list)

    def __getitem__(self, idx):
        file_path = self.file_list[idx]
        data = None
        try:
            # Try to load, handling corrupted files
            try:
                data = np.load(file_path)
            except ValueError as ve:
                if "mmap length is greater than file size" in str(ve) or "cannot reshape" in str(ve):
                    # File has shape mismatch, try to load raw data and reshape correctly
                    print(f"Warning: {file_path} has shape mismatch, attempting manual load")
                    with open(file_path, 'rb') as f:
                        version = np.lib.format.read_magic(f)
                        shape, fortran_order, dtype = np.lib.format.read_array_header_1_0(f)
                        # Read only available data
                        remaining_data = f.read()
                        expected_bytes = np.prod(shape) * np.dtype(dtype).itemsize
                        actual_bytes = len(remaining_data)
                        print(f"Header shape: {shape}, expected {expected_bytes} bytes, got {actual_bytes} bytes")
                        
                        if actual_bytes > 0:
                            # Load what we can
                            elements = actual_bytes // np.dtype(dtype).itemsize
                            flat_data = np.frombuffer(remaining_data[:elements * np.dtype(dtype).itemsize], dtype=dtype)
                            
                            # Try to infer actual dimensions from file size
                            # Common patterns: (957,960,300), (960,960,300), etc.
                            if len(shape) == 3:
                                h_expected, w_expected, c_expected = shape
                                # Calculate what the actual height should be
                                actual_h = elements // (w_expected * c_expected)
                                if actual_h * w_expected * c_expected == elements:
                                    print(f"Reshaping to ({actual_h}, {w_expected}, {c_expected}) instead of {shape}")
                                    data = flat_data.reshape(actual_h, w_expected, c_expected)
                                else:
                                    # Try other common dimensions
                                    for c in [300, 600, 900]:
                                        for w in [960, 957, 950]:
                                            h = elements // (w * c)
                                            if h * w * c == elements and h > 0:
                                                print(f"Reshaping to ({h}, {w}, {c})")
                                                data = flat_data.reshape(h, w, c)
                                                break
                                        if data is not None:
                                            break
                            
                            # If still no success, try fallback reshape
                            if data is None and elements >= 256 * 256:
                                # Find best square dimensions
                                sqrt_elem = int(np.sqrt(elements))
                                while sqrt_elem > 0 and elements % sqrt_elem != 0:
                                    sqrt_elem -= 1
                                if sqrt_elem > 0:
                                    height = sqrt_elem
                                    width = elements // sqrt_elem
                                    data = flat_data[:height*width].reshape(height, width, 1)
                                else:
                                    # Fallback to small square
                                    side = int(np.sqrt(elements))
                                    data = flat_data[:side*side].reshape(side, side, 1)
                            elif data is None:
                                # Very small data, create minimal image
                                data = np.zeros((64, 64, 1), dtype=np.float32)
                        else:
                            raise ValueError("No data to load")
                else:
                    raise ve
            
            # Handle potential shape mismatches with robust reshaping
            if data.ndim == 1:
                # Try to infer shape from file size - assume square spatial dimensions
                total_size = data.size
                # Try common channel counts first
                for channels in [300, 600, 900, 1200]:
                    spatial_size = total_size // channels
                    spatial_dim = int(np.sqrt(spatial_size))
                    if spatial_dim * spatial_dim * channels == total_size:
                        data = data.reshape(spatial_dim, spatial_dim, channels)
                        break
                else:
                    # Fallback: treat as single channel
                    spatial_dim = int(np.sqrt(total_size))
                    if spatial_dim * spatial_dim == total_size:
                        data = data.reshape(spatial_dim, spatial_dim, 1)
                    else:
                        # Force to a reasonable shape
                        target_size = min(256, spatial_dim)
                        data = data[:target_size*target_size].reshape(target_size, target_size, 1)
            elif data.ndim == 2:
                # Add channel dimension
                data = data[:, :, np.newaxis]
            # data should now be (H, W, C)
            
            data = self.resize_data(data)
            data = torch.tensor(data, dtype=torch.float32).permute(2, 0, 1).contiguous()
            # Ensure consistent channel size across samples
            if self.target_channels is not None:
                c, h, w = data.shape
                if c > self.target_channels:
                    data = data[:self.target_channels, :, :].contiguous()
                elif c < self.target_channels:
                    pad = torch.zeros(self.target_channels - c, h, w, dtype=data.dtype)
                    data = torch.cat([data, pad], dim=0).contiguous()
            if self.transform:
                data = self.transform(data)
            label = torch.tensor(self.labels[idx], dtype=torch.float32)
            raw_value = torch.tensor(self.raw_values[idx], dtype=torch.float32)
            return data, label, raw_value
        except Exception as e:
            print(f"Error loading {file_path}: {e}")
            if data is not None:
                print(f"Array shape: {data.shape}, size: {data.size}")
            else:
                print("Failed to load array")
            # Return a dummy tensor to avoid breaking the DataLoader
            dummy_data = torch.zeros(self.target_channels or 300, *self.target_size, dtype=torch.float32)
            dummy_label = torch.tensor(0.0, dtype=torch.float32)
            dummy_raw = torch.tensor(0.0, dtype=torch.float32)
            return dummy_data, dummy_label, dummy_raw

    def resize_data(self, data):
        """Resize images using PyTorch's efficient interpolation"""
        if data.ndim == 2:
            data = data[:, :, np.newaxis]
        data_tensor = torch.tensor(data, dtype=torch.float32).permute(2, 0, 1)
        resized = F.interpolate(
            data_tensor.unsqueeze(0),
            size=self.target_size,
            mode='bilinear',
            align_corners=False
        )
        # Return a fresh, contiguous numpy array to avoid non-resizable storage issues in DataLoader
        return resized.squeeze(0).permute(1, 2, 0).contiguous().cpu().numpy().copy()

def compute_or_load_stats(dataset, cache_path='stats.json', label_scaler=None):
    """Compute or load statistics including normalization parameters for targets"""
    # Check if cache exists
    if os.path.exists(cache_path):
        with open(cache_path, 'r') as f:
            stats = json.load(f)
        print(f"Loaded statistics from {cache_path}")
        
        # Check if cache contains toxin normalization parameters
        has_toxin_params = 'toxin_min' in stats and 'toxin_max' in stats
        
        if has_toxin_params:
            print("Cache contains toxin normalization parameters")
            return (
                stats['mean'], 
                stats['std'], 
                stats['toxin_min'], 
                stats['toxin_max']
            )
        else:
            print("⚠️ Warning: Cache does not contain toxin normalization parameters, will recompute statistics")
    
    # Recompute statistics
    print("Computing dataset statistics (mean and std)... This may take some time")
    # Use single worker for stability on macOS
    loader = DataLoader(dataset, batch_size=16, shuffle=False, num_workers=0, pin_memory=False)
    
    try:
        channels = dataset[0][0].shape[0]
    except IndexError:
        raise RuntimeError("Dataset is empty, cannot compute statistics")

    total = torch.zeros(channels)
    total_sq = torch.zeros(channels)
    num_pixels = 0

    for data, _, _ in tqdm(loader, desc="Computing statistics"):
        data = data.view(data.size(0), channels, -1)
        total += data.sum(dim=[0, 2])
        total_sq += (data ** 2).sum(dim=[0, 2])
        num_pixels += data.shape[0] * data.shape[2]

    mean = (total / num_pixels).tolist()
    std = ((total_sq / num_pixels - torch.tensor(mean) ** 2) ** 0.5).tolist()

    # Save target normalization parameters (min, max)
    toxin_min = float(label_scaler.data_min_[0]) if label_scaler else 0.0
    toxin_max = float(label_scaler.data_max_[0]) if label_scaler else 1.0

    stats = {
        'mean': mean,
        'std': std,
        'toxin_min': toxin_min,
        'toxin_max': toxin_max
    }

    try:
        with open(cache_path, 'w') as f:
            json.dump(stats, f)
        print(f"Statistics computed and saved to {cache_path}")
    except OSError as e:
        if e.errno == 28:  # No space left on device
            print(f"Warning: No space left on device, skipping stats cache. Stats computed successfully.")
        else:
            print(f"Warning: Failed to save stats cache: {e}")
    
    return mean, std, toxin_min, toxin_max

def _infer_common_channels(npy_paths, sample_limit=32):
    """Infer a common channel count across a subset of files (min channels)."""
    chans = []
    for p in npy_paths[:sample_limit]:
        try:
            arr = np.load(p, mmap_mode='r')
            if arr.ndim == 1:
                # Try to infer channels from 1D array
                total_size = arr.size
                for channels in [300, 600, 900, 1200]:
                    spatial_size = total_size // channels
                    spatial_dim = int(np.sqrt(spatial_size))
                    if spatial_dim * spatial_dim * channels == total_size:
                        chans.append(channels)
                        break
                else:
                    # Assume single channel if can't determine
                    spatial_dim = int(np.sqrt(total_size))
                    if spatial_dim * spatial_dim == total_size:
                        chans.append(1)
            elif arr.ndim == 2:
                chans.append(1)
            else:
                chans.append(arr.shape[2])
        except Exception as e:
            print(f"Warning: Failed to process {p}: {e}")
            continue
    if not chans:
        return None
    common_channels = min(chans)
    print(f"Detected channels: {chans[:10]}... (first 10), using min: {common_channels}")
    return int(common_channels)

def preprocess_data(data_folder, batch_size=16, target_size=(256, 256)):
    metadata_path = os.path.join(data_folder, 'new_metadata.csv')
    if not os.path.exists(metadata_path):
        raise FileNotFoundError(f"Metadata file {metadata_path} not found")
    
    metadata = pd.read_csv(metadata_path, encoding='utf-8-sig')
    
    # Check required columns
    required_base = ['toxin_value', 'Value_Raw']
    for col in required_base:
        if col not in metadata.columns:
            raise ValueError(f"Required column '{col}' not found in metadata.csv")
    if ('npy_filename' not in metadata.columns) and ('filename' not in metadata.columns):
        raise ValueError("metadata.csv must contain 'npy_filename' or 'filename'")

    # Create label scaler using MinMaxScaler to [0,1] range
    toxin_values = metadata['toxin_value'].values.reshape(-1, 1)
    label_scaler = MinMaxScaler(feature_range=(0, 1))
    label_scaler.fit(toxin_values)
    
    print(f"Toxin value normalization parameters: min={label_scaler.data_min_[0]}, max={label_scaler.data_max_[0]}")

    # Create feature scaler
    raw_values = metadata['Value_Raw'].values.reshape(-1, 1)
    feature_scaler = StandardScaler().fit(raw_values)
    print(f"Value_Raw normalization: mean={feature_scaler.mean_[0]}, std={feature_scaler.scale_[0]}")

    # Peek metadata to infer channels
    meta = pd.read_csv(metadata_path, encoding='utf-8-sig')
    npy_col = 'npy_filename' if 'npy_filename' in meta.columns else 'filename'
    npy_paths = []
    for _, r in meta.iterrows():
        p = str(r[npy_col]).strip()
        if not os.path.isabs(p):
            p = os.path.join(data_folder, p)
        if os.path.isfile(p):
            npy_paths.append(os.path.normpath(p))
    target_channels = _infer_common_channels(npy_paths)
    if target_channels is None:
        raise RuntimeError("Failed to infer common channel count from npy files")
    print(f"Using common spectral channels: {target_channels}")

    # Create temporary dataset for statistics
    full_dataset_for_stats = HyperspectralDataset(
        data_folder,
        target_size=target_size,
        scaler=None,
        feature_scaler=None,
        target_channels=target_channels
    )
    
    # Compute or load mean and std, and save toxin normalization parameters
    mean, std, toxin_min, toxin_max = compute_or_load_stats(
        full_dataset_for_stats, 
        os.path.join(data_folder, 'stats.json'),
        label_scaler=label_scaler
    )
    
    print(f"Toxin value normalization parameters saved: min={toxin_min}, max={toxin_max}")

    # Data augmentation
    train_transform = transforms.Compose([
        transforms.RandomHorizontalFlip(),
        transforms.RandomVerticalFlip(),
        transforms.RandomRotation(degrees=15),
        transforms.RandomAffine(degrees=0, translate=(0.05, 0.05)),
        transforms.Normalize(mean, std)
    ])
    test_transform = transforms.Compose([
        transforms.Normalize(mean, std)
    ])

    # Create full datasets
    full_dataset_train = HyperspectralDataset(
        data_folder, 
        target_size=target_size, 
        transform=train_transform, 
        scaler=label_scaler,
        feature_scaler=feature_scaler,
        target_channels=target_channels
    )
    full_dataset_test = HyperspectralDataset(
        data_folder, 
        target_size=target_size, 
        transform=test_transform, 
        scaler=label_scaler,
        feature_scaler=feature_scaler,
        target_channels=target_channels
    )

    # Split into train and test sets
    train_idx, test_idx = train_test_split(
        np.arange(len(full_dataset_train)), 
        test_size=0.2, 
        random_state=24
    )
    
    # Create subsets
    train_set = Subset(full_dataset_train, train_idx)
    test_set = Subset(full_dataset_test, test_idx)

    # Create data loaders
    nw = min(4, (os.cpu_count() or 1))
    train_loader = DataLoader(
        train_set,
        batch_size=batch_size,
        shuffle=True,
        num_workers=nw,
        pin_memory=True,
        persistent_workers=(nw > 0)
    )
    test_loader = DataLoader(
        test_set,
        batch_size=batch_size,
        shuffle=False,
        num_workers=nw,
        pin_memory=True,
        persistent_workers=(nw > 0)
    )

    return train_loader, test_loader, label_scaler, feature_scaler

# ========================== CBAM Module Implementation ==========================
class ChannelAttention(nn.Module):
    """Channel Attention Module"""
    def __init__(self, in_channels, reduction_ratio=16):
        super(ChannelAttention, self).__init__()
        self.avg_pool = nn.AdaptiveAvgPool2d(1)
        self.max_pool = nn.AdaptiveMaxPool2d(1)
        
        self.fc = nn.Sequential(
            nn.Linear(in_channels, in_channels // reduction_ratio),
            nn.ReLU(inplace=True),
            nn.Linear(in_channels // reduction_ratio, in_channels),
            nn.Sigmoid()
        )

    def forward(self, x):
        avg_out = self.fc(self.avg_pool(x).view(x.size(0), -1))
        max_out = self.fc(self.max_pool(x).view(x.size(0), -1))
        out = avg_out + max_out
        return out.view(x.size(0), x.size(1), 1, 1)

class SpatialAttention(nn.Module):
    """Spatial Attention Module"""
    def __init__(self, kernel_size=7):
        super(SpatialAttention, self).__init__()
        assert kernel_size in (3, 7), "kernel size must be 3 or 7"
        padding = 3 if kernel_size == 7 else 1
        
        self.conv = nn.Conv2d(2, 1, kernel_size, padding=padding, bias=False)
        self.sigmoid = nn.Sigmoid()

    def forward(self, x):
        avg_out = torch.mean(x, dim=1, keepdim=True)
        max_out, _ = torch.max(x, dim=1, keepdim=True)
        x = torch.cat([avg_out, max_out], dim=1)
        x = self.conv(x)
        return self.sigmoid(x)

class CBAM(nn.Module):
    """CBAM Module: Channel Attention + Spatial Attention"""
    def __init__(self, in_channels, reduction_ratio=16, kernel_size=7):
        super(CBAM, self).__init__()
        self.channel_att = ChannelAttention(in_channels, reduction_ratio)
        self.spatial_att = SpatialAttention(kernel_size)

    def forward(self, x):
        x = x * self.channel_att(x)
        x = x * self.spatial_att(x)
        return x

# ========================== Enhanced Modules ==========================
class PositionalEncoding(nn.Module):
    """Positional Encoding for enhanced spectral information"""
    def __init__(self, d_model, max_len=500):
        super().__init__()
        pe = torch.zeros(max_len, d_model)
        position = torch.arange(0, max_len, dtype=torch.float).unsqueeze(1)
        div_term = torch.exp(torch.arange(0, d_model, 2).float() * (-math.log(10000.0) / d_model))
        pe[:, 0::2] = torch.sin(position * div_term)
        pe[:, 1::2] = torch.cos(position * div_term)
        pe = pe.unsqueeze(0)
        self.register_buffer('pe', pe)
        
    def forward(self, x):
        seq_len = x.size(1)
        pos_emb = self.pe[:, :seq_len, :]
        return x + pos_emb

class GatedLinear(nn.Module):
    """Gated Linear Layer"""
    def __init__(self, in_features, out_features):
        super().__init__()
        self.linear = nn.Linear(in_features, out_features)
        self.gate = nn.Linear(in_features, out_features)
        nn.init.xavier_uniform_(self.linear.weight)
        nn.init.xavier_uniform_(self.gate.weight)
        nn.init.zeros_(self.linear.bias)
        nn.init.zeros_(self.gate.bias)
        
    def forward(self, x):
        return self.linear(x) * torch.sigmoid(self.gate(x))

# ========================== CBAM Model ==========================
class HyperspectralModel(nn.Module):
    def __init__(self, input_channels):
        super(HyperspectralModel, self).__init__()
        
        # Initial feature extraction
        self.conv1 = nn.Sequential(
            nn.Conv2d(input_channels, 64, kernel_size=3, padding=1),
            nn.BatchNorm2d(64),
            nn.ReLU(inplace=True)
        )
        self.cbam1 = CBAM(64)
        
        # Downsample block 1
        self.down1 = nn.Sequential(
            nn.Conv2d(64, 64, kernel_size=3, padding=1),
            nn.BatchNorm2d(64),
            nn.ReLU(inplace=True),
            nn.MaxPool2d(2),
            CBAM(64)
        )
        
        # Downsample block 2
        self.down2 = nn.Sequential(
            nn.Conv2d(64, 128, kernel_size=3, padding=1),
            nn.BatchNorm2d(128),
            nn.ReLU(inplace=True),
            nn.Conv2d(128, 128, kernel_size=3, padding=1),
            nn.BatchNorm2d(128),
            nn.ReLU(inplace=True),
            nn.MaxPool2d(2),
            CBAM(128)
        )
        
        # Final feature extraction
        self.final_conv = nn.Sequential(
            nn.Conv2d(128, 256, kernel_size=3, padding=1),
            nn.BatchNorm2d(256),
            nn.ReLU(inplace=True),
            CBAM(256),
            nn.AdaptiveAvgPool2d(1)
        )
        
        # Fusion FC layers with Value_Raw feature
        self.fc = nn.Sequential(
            nn.Linear(256 + 1, 512),  # 256 image features + 1 Value_Raw
            nn.ReLU(inplace=True),
            nn.Dropout(0.4),
            
            nn.Linear(512, 256),
            nn.ReLU(inplace=True),
            nn.Dropout(0.3),
            
            nn.Linear(256, 128),
            nn.ReLU(inplace=True),
            nn.Dropout(0.2),
            
            nn.Linear(128, 1)
        )
    
    def forward(self, x, feature):
        # Feature extraction + CBAM attention
        x = self.conv1(x)
        x = self.cbam1(x)
        x = self.down1(x)
        x = self.down2(x)
        img_features = self.final_conv(x).squeeze(-1).squeeze(-1)
        
        # Concatenate with Value_Raw feature
        combined = torch.cat([img_features, feature.unsqueeze(1)], dim=1)
        
        # Predict through FC layers
        return self.fc(combined).flatten()
    
    def visualize_attention(self, input_tensor):
        """Visualize CBAM attention effects (for debugging)"""
        with torch.no_grad():
            activations = {}
            
            # First layer
            x1 = self.conv1(input_tensor)
            x1_att = self.cbam1(x1)
            activations['conv1'] = x1
            activations['cbam1'] = x1_att
            
            # Downsample 1
            x2 = self.down1[:-1](x1_att)
            x2_att = self.down1[-1](x2)
            activations['down1'] = x2
            activations['cbam2'] = x2_att
            
            # Downsample 2
            x3 = self.down2[:-1](x2_att)
            x3_att = self.down2[-1](x3)
            activations['down2'] = x3
            activations['cbam3'] = x3_att
            
            # Final layer
            x4 = self.final_conv[:-2](x3_att)
            x4_att = self.final_conv[-2](x4)
            activations['final_conv'] = x4
            activations['cbam4'] = x4_att
            
            return activations

# ========================== Training Functions ==========================
def adjust_learning_rate(optimizer, epoch, warmup_epochs, initial_lr, num_epochs):
    """Custom learning rate scheduler"""
    if epoch < warmup_epochs:
        lr = initial_lr * (epoch + 1) / warmup_epochs
    else:
        progress = (epoch - warmup_epochs) / (num_epochs - warmup_epochs)
        lr = 0.5 * initial_lr * (1 + math.cos(math.pi * progress))
    
    for param_group in optimizer.param_groups:
        param_group['lr'] = lr
    return lr

def evaluate_regression(model, loader, device='cuda', scaler=None):
    """Evaluate regression model performance"""
    model.eval()
    all_labels, all_preds = [], []
    with torch.no_grad():
        for inputs, labels, features in loader:
            inputs = inputs.to(device)
            features = features.to(device)
            outputs = model(inputs, features).flatten()
            all_labels.extend(labels.cpu().numpy())
            all_preds.extend(outputs.cpu().numpy())
    
    all_labels = np.array(all_labels)
    all_preds = np.array(all_preds)

    # Inverse transform if scaler is provided
    if scaler:
        all_labels_orig = scaler.inverse_transform(all_labels.reshape(-1, 1)).flatten()
        all_preds_orig = scaler.inverse_transform(all_preds.reshape(-1, 1)).flatten()
    else:
        all_labels_orig = all_labels
        all_preds_orig = all_preds

    # Calculate metrics
    mae = mean_absolute_error(all_labels_orig, all_preds_orig)
    r2 = r2_score(all_labels_orig, all_preds_orig)
    mse_normalized = np.mean((all_labels - all_preds)**2)

    return all_labels_orig, all_preds_orig, mae, r2, mse_normalized

def train_regression_model(model, train_loader, val_loader, scaler, num_epochs=100, device='cpu'):
    """Train regression model (evaluate every epoch)"""
    # Move model to the selected device
    model.to(device)
    criterion = nn.MSELoss()
    optimizer = optim.AdamW(model.parameters(), lr=1e-4, weight_decay=1e-5)
    
    # Training parameters
    warmup_epochs = 10
    best_val_loss = float('inf')
    best_val_r2 = -float('inf')
    patience = 15
    trigger_times = 0

    # Track training progress
    train_losses = []
    val_losses = []
    val_r2s = []
    best_predictions = None

    print(f"Starting training for {num_epochs} epochs with evaluation every epoch...")
    for epoch in range(num_epochs):
        model.train()
        running_loss = 0.0
        
        # Adjust learning rate
        current_lr = adjust_learning_rate(optimizer, epoch, warmup_epochs, 1e-4, num_epochs)
        
        # Training loop
        for inputs, labels, features in tqdm(train_loader, desc=f'Epoch {epoch + 1}/{num_epochs}'):
            inputs = inputs.to(device)
            labels = labels.to(device)
            features = features.to(device)
            
            optimizer.zero_grad()
            outputs = model(inputs, features)
            loss = criterion(outputs, labels)
            loss.backward()
            
            # Gradient clipping
            torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=2.0)
            
            optimizer.step()
            running_loss += loss.item()
        
        # Calculate average training loss
        avg_train_loss = running_loss / len(train_loader)
        train_losses.append(avg_train_loss)
        
        # Evaluate every epoch
        model.eval()
        with torch.no_grad():
            val_labels, val_preds = [], []
            for val_inputs, val_labels_batch, val_features in val_loader:
                val_inputs = val_inputs.to(device)
                val_features = val_features.to(device)
                val_outputs = model(val_inputs, val_features).flatten()
                val_labels.extend(val_labels_batch.cpu().numpy())
                val_preds.extend(val_outputs.cpu().numpy())
            
            val_labels = np.array(val_labels)
            val_preds = np.array(val_preds)
            
            # Inverse transform if scaler is provided
            if scaler:
                val_labels_orig = scaler.inverse_transform(val_labels.reshape(-1, 1)).flatten()
                val_preds_orig = scaler.inverse_transform(val_preds.reshape(-1, 1)).flatten()
            else:
                val_labels_orig = val_labels
                val_preds_orig = val_preds

            # Calculate metrics
            val_mae = mean_absolute_error(val_labels_orig, val_preds_orig)
            val_r2 = r2_score(val_labels_orig, val_preds_orig)
            val_loss_normalized = np.mean((val_labels - val_preds)**2)
        
        val_losses.append(val_loss_normalized)
        val_r2s.append(val_r2)
        
        # Check for improvement
        r2_improved = (val_r2 - best_val_r2) > 0.001
        loss_improved = (best_val_loss - val_loss_normalized) > 0.01
        improvement = r2_improved or loss_improved
        
        if improvement:
            best_val_r2 = val_r2
            best_val_loss = val_loss_normalized
            trigger_times = 0
            best_predictions = (val_labels_orig, val_preds_orig)
            
            # Save model
            save_dict = {
                'epoch': epoch,
                'model_state_dict': model.state_dict(),
                'optimizer_state_dict': optimizer.state_dict(),
                'loss': best_val_loss,
                'r2': best_val_r2,
            }
            
            # Add toxin normalization parameters
            if hasattr(scaler, 'data_min_') and hasattr(scaler, 'data_max_'):
                save_dict['toxin_min'] = scaler.data_min_[0]
                save_dict['toxin_max'] = scaler.data_max_[0]
            
            try:
                torch.save(save_dict, 'best_model.pth')
                print(f"✔ Validation metrics improved, saving model. R²: {best_val_r2:.4f}, Loss: {val_loss_normalized:.4f} (Epoch {epoch+1})")
            except OSError as e:
                if e.errno == 28:  # No space left on device
                    print(f"✔ Validation metrics improved but cannot save model (no space). R²: {best_val_r2:.4f}, Loss: {val_loss_normalized:.4f} (Epoch {epoch+1})")
                else:
                    print(f"✔ Validation metrics improved but failed to save model: {e}. R²: {best_val_r2:.4f}, Loss: {val_loss_normalized:.4f} (Epoch {epoch+1})")
        else:
            trigger_times += 1
            if trigger_times >= patience:
                print(f"Early stopping triggered! No improvement for {patience} consecutive epochs.")
                break
        
        # Print training status
        print(f"Epoch {epoch + 1}/{num_epochs} | LR: {current_lr:.2e} | "
              f"Train Loss: {avg_train_loss:.4f} | "
              f"Val Loss: {val_loss_normalized:.4f} | "
              f"Val R²: {val_r2:.4f} | "
              f"Trigger: {trigger_times}/{patience}")
    
    # Training complete, load best model
    if os.path.exists('best_model.pth'):
        checkpoint = torch.load('best_model.pth')
        model.load_state_dict(checkpoint['model_state_dict'])
        print(f"\nTraining complete. Best model at epoch {checkpoint['epoch'] + 1}: "
              f"R²: {checkpoint['r2']:.4f}, Loss: {checkpoint['loss']:.4f}")
        if 'toxin_min' in checkpoint and 'toxin_max' in checkpoint:
            print(f"Toxin min={checkpoint['toxin_min']}, Toxin max={checkpoint['toxin_max']}")
    else:
        print("\nTraining complete, but no model saved.")
    
    # Visualize training progress
    visualize_training(train_losses, val_losses, val_r2s, min(num_epochs, len(val_losses)))
    
    return best_predictions

def visualize_training(train_losses, val_losses, val_r2s, num_epochs):
    """Visualize training progress"""
    plt.figure(figsize=(12, 8))
    
    # Loss curves
    epochs = np.arange(1, len(train_losses) + 1)
    plt.subplot(2, 1, 1)
    plt.plot(epochs, train_losses, 'b-', label='Train Loss')
    if val_losses:
        eval_points = np.linspace(0, num_epochs, len(val_losses), endpoint=False).astype(int) + 1
        plt.plot(eval_points, val_losses, 'r-', label='Val Loss')
    plt.xlabel('Epoch')
    plt.ylabel('Loss')
    plt.title('Training and Validation Loss')
    plt.legend()
    plt.grid(True)
    
    # R² curve
    plt.subplot(2, 1, 2)
    if val_r2s:
        plt.plot(eval_points, val_r2s, 'g-', label='Val R²')
    plt.axhline(y=0, color='k', linestyle='--', alpha=0.5)
    plt.xlabel('Epoch')
    plt.ylabel('R² Score')
    plt.title('Validation R² Score')
    plt.ylim(-0.5, 1.0)
    plt.grid(True)
    
    plt.tight_layout()
    plt.savefig('training_metrics.png', dpi=300)
    plt.show()

def enhanced_visualization(true_values, pred_values):
    """Enhanced regression result visualization"""
    plt.figure(figsize=(18, 12))
    residuals = pred_values - true_values
    abs_errors = np.abs(residuals)
    
    # Figure 1: True vs Predicted
    plt.subplot(2, 2, 1)
    sns.regplot(x=true_values, y=pred_values, 
                scatter_kws={'alpha':0.4, 'color':'#4B8BBE'}, 
                line_kws={'color':'#FF0000', 'lw':2},
                ci=95)
    
    # Add metrics
    mae = mean_absolute_error(true_values, pred_values)
    r2 = r2_score(true_values, pred_values)
    rmse = np.sqrt(mean_squared_error(true_values, pred_values))
    textstr = '\n'.join((
        f'MAE = {mae:.2f}',
        f'RMSE = {rmse:.2f}',
        f'R² = {r2:.2f}'))
    
    plt.gca().text(0.05, 0.95, textstr, transform=plt.gca().transAxes,
                   fontsize=12, verticalalignment='top',
                   bbox=dict(facecolor='white', alpha=0.8))
    
    plt.plot([min(true_values), max(true_values)], 
             [min(true_values), max(true_values)], 
             'k--', lw=1, label='Perfect Prediction')
    plt.xlabel('True Values')
    plt.ylabel('Predictions')
    plt.title('True vs Predicted Values')
    plt.legend()
    
    # Figure 2: Residual analysis
    plt.subplot(2, 2, 2)
    sns.residplot(x=true_values, y=residuals, 
                  lowess=True, 
                  scatter_kws={'alpha':0.4, 'color':'#4B8BBE'},
                  line_kws={'color':'#FF0000', 'lw':2})
    
    plt.axhline(y=0, color='k', linestyle='--', lw=1)
    plt.xlabel('True Values')
    plt.ylabel('Residuals')
    plt.title('Residual Analysis')
    
    # Figure 3: Error distribution
    plt.subplot(2, 2, 3)
    sns.histplot(abs_errors, kde=True, 
                 bins=30, color='#4B8BBE',
                 edgecolor='white', linewidth=0.5)
    
    plt.axvline(x=mae, color='r', linestyle='--', lw=2,
                label=f'MAE = {mae:.2f}')
    plt.xlabel('Absolute Error')
    plt.ylabel('Count')
    plt.title('Absolute Error Distribution')
    plt.legend()
    
    # Figure 4: Prediction error
    plt.subplot(2, 2, 4)
    plt.scatter(true_values, residuals, alpha=0.5, color='#4B8BBE')
    plt.axhline(y=0, color='r', linestyle='--', lw=1)
    plt.xlabel('True Values')
    plt.ylabel('Prediction Error')
    plt.title('Prediction Error Distribution')
    plt.grid(True)
    plt.tight_layout()
    plt.savefig('enhanced_regression_analysis.png', dpi=300)
    plt.show()

def main():
    data_folder = "/Volumes/Extreme Pro/001-实验数据-是这个/001-是这个/test/processed_npy" 
    batch_size = 16
    num_epochs = 100

    print("--- Data Preprocessing ---")
    train_loader, test_loader, label_scaler, feature_scaler = preprocess_data(data_folder, batch_size)
    
    # Select device automatically
    device = get_best_device()
    print(f"Using device: {device}")
    
    try:
        sample_data, _, _ = next(iter(train_loader))
        input_channels = sample_data.shape[1]
        print(f"Detected input channels: {input_channels}")
    except StopIteration:
        print("Train loader is empty. Check data folder and metadata file.")
        return

    print("--- Model Creation ---")
    model = HyperspectralModel(input_channels)
    
    # Print model parameters
    total_params = sum(p.numel() for p in model.parameters())
    trainable_params = sum(p.numel() for p in model.parameters() if p.requires_grad)
    print(f"Total parameters: {total_params:,} | Trainable parameters: {trainable_params:,}")

    print("--- Model Training ---")
    best_predictions = train_regression_model(
        model, 
        train_loader, 
        test_loader, 
        label_scaler, 
        num_epochs=num_epochs,
        device=device
    )

    print("\n--- Final Model Evaluation ---")
    # Load best model for final evaluation
    if os.path.exists('best_model.pth'):
        checkpoint = torch.load('best_model.pth')
        model.load_state_dict(checkpoint['model_state_dict'])
        
        final_labels_orig, final_preds_orig, final_mae, final_r2, _ = evaluate_regression(
            model, test_loader, device=device, scaler=label_scaler
        )

        print(f"Final performance on test set:")
        print(f"MAE (original scale): {final_mae:.4f}")
        print(f"R² Score (original scale): {final_r2:.4f}")
        print(f"Best model achieved R²: {checkpoint['r2']:.4f} at epoch {checkpoint['epoch'] + 1}")

        # Visualize predictions
        enhanced_visualization(final_labels_orig, final_preds_orig)
        
        # Visualize best predictions (if available)
        if best_predictions:
            print("\nVisualizing best predictions...")
            best_labels, best_preds = best_predictions
            enhanced_visualization(best_labels, best_preds)
    else:
        print("No saved model found, using current model for evaluation")
        final_labels_orig, final_preds_orig, final_mae, final_r2, _ = evaluate_regression(
            model, test_loader, device=device, scaler=label_scaler
        )
        print(f"Final performance on test set:")
        print(f"MAE (original scale): {final_mae:.4f}")
        print(f"R² Score (original scale): {final_r2:.4f}")
        enhanced_visualization(final_labels_orig, final_preds_orig)

if __name__ == '__main__':
    main()

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
ENVI (.hdr/.spe) 转换为 .npy + 生成元数据 CSV 的预处理脚本

使用示例：
	python scripts/prepare_npy_from_envi.py \
		--excel "/Volumes/Extreme Pro/001-实验数据-是这个/001-是这个/20251015-数据统计-是这个.xlsx" \
		--data_root "/Volumes/Extreme Pro/001-实验数据-是这个/001-是这个/001-高光谱-是这个" \
		--out_dir "/Volumes/Extreme Pro/001-实验数据-是这个/001-是这个/test/processed_npy" \
		--metadata "/Volumes/Extreme Pro/001-实验数据-是这个/001-是这个/test/processed_npy/new_metadata.csv"

功能概述：
	1) 从 Excel 读取三列：
		 - Sample Name：用于定位 ENVI 文件（支持每 40 个样本分卷：如 6H、6H-2，编号在第二卷从 1 重新开始）
		 - Area：写入 Value_Raw
		 - Toxic content：写入 toxin_value
	2) 在 data_root 下寻找对应 .hdr/.spe，读取为 numpy 数组 (H,W,B) float32
	3) 在 out_dir 下以“镜像目录结构”保存为 .npy（例如 14H/14-1.npy 或 46H-2/46-1.npy）
	4) 写出 --metadata 指定的 CSV，字段包括：
		 - filename：样本在原始数据集中的规范路径（绝对路径），如 /.../6H-2/6-1
		 - npy_filename：对应 .npy 在 out_dir 下的相对路径，如 6H-2/6-1.npy
		 - Value_Raw, toxin_value, sample_name
	5) 可选计算数据集通道均值/方差（--compute-stats），写入 stats.json

续跑/容错：
	- --resume：跳过已存在的 .npy，并合并已有的 CSV，支持从中断处继续。
	- --start-index：从指定（Excel 清洗后）行号开始（1-based）。
	- 磁盘满（ENOSPC）时，脚本会先安全落盘当前 CSV 再停止，释放空间后加 --resume 继续。
"""

from __future__ import annotations

import os
import sys
import csv
import json
import math
import glob
import argparse
from typing import Optional, Tuple, List, Dict

import numpy as np
import pandas as pd
import errno


CHUNK_SIZE = 40  # number of samples per "H" folder before rolling to "-2"


def log(msg: str):
	print(msg, flush=True)


def ensure_dir(path: str):
	os.makedirs(path, exist_ok=True)


def parse_excel(excel_path: str) -> pd.DataFrame:
	"""Load the Excel and normalize the required columns.

	We try to support a few possible header variants by case-insensitive matching and
	fallback aliases. Required logical fields:
	  - sample_name
	  - area (Value_Raw)
	  - toxic (toxin_value)
	"""
	df = pd.read_excel(excel_path)
	# normalize column names for robust matching
	col_map = {c: str(c).strip() for c in df.columns}
	df = df.rename(columns=col_map)
	lower_cols = {c.lower(): c for c in df.columns}

	def pick(*candidates) -> Optional[str]:
		for cand in candidates:
			key = cand.lower()
			if key in lower_cols:
				return lower_cols[key]
		# fallback: substring contains
		for c in df.columns:
			lc = c.lower()
			if any(k in lc for k in [cand.lower() for cand in candidates]):
				return c
		return None

	col_sample = pick("Sample Name", "Sample Name(filename)", "filename", "Sample", "样本")
	col_area = pick("Area", "Area(value raw)", "Value_Raw", "value raw", "面积")
	col_toxic = pick("Toxic content", "toxin_value", "Toxic", "毒素", "毒性")

	missing = [
		name for name, val in [
			("Sample Name", col_sample), ("Area", col_area), ("Toxic content", col_toxic)
		] if val is None
	]
	if missing:
		raise ValueError(f"Excel missing required columns: {missing}. Found columns={list(df.columns)}")

	out = pd.DataFrame({
		"sample_name": df[col_sample].astype(str).str.strip(),
		"Value_Raw": pd.to_numeric(df[col_area], errors="coerce"),
		"toxin_value": pd.to_numeric(df[col_toxic], errors="coerce"),
	})
	# drop rows with missing numeric values
	out = out.dropna(subset=["Value_Raw", "toxin_value"]).reset_index(drop=True)
	return out


def parse_hdr(hdr_path: str) -> Dict[str, object]:
	meta: Dict[str, object] = {}
	with open(hdr_path, "r", encoding="utf-8", errors="ignore") as f:
		for raw in f:
			line = raw.strip()
			if not line or "=" not in line:
				continue
			k, v = [x.strip() for x in line.split("=", 1)]
			kl = k.lower()
			if v.startswith("{") and v.endswith("}"):
				v = v[1:-1].strip()
			# try numeric
			try:
				if "." in v:
					vf = float(v)
					if vf.is_integer():
						vf = int(vf)
					meta[kl] = vf
				else:
					meta[kl] = int(v)
				continue
			except Exception:
				pass
			meta[kl] = v

	def get_num(key: str, default=None):
		for cand in (key, key.replace(" ", "")):
			if cand in meta and isinstance(meta[cand], (int, float)):
				return int(meta[cand])
		return default

	samples = get_num("samples")
	lines = get_num("lines")
	bands = get_num("bands")
	interleave = str(meta.get("interleave", meta.get("interleave", "bsq"))).lower()
	data_type = str(meta.get("data type", meta.get("datatype", "4")))
	return {"samples": samples, "lines": lines, "bands": bands, "interleave": interleave, "data type": data_type}


def read_spe(spe_path: str, meta: Dict[str, object]) -> np.ndarray:
	samples = int(meta["samples"]) if meta["samples"] is not None else None
	lines = int(meta["lines"]) if meta["lines"] is not None else None
	bands = int(meta["bands"]) if meta["bands"] is not None else None
	if None in (samples, lines, bands):
		raise ValueError("Missing samples/lines/bands in HDR metadata")
	interleave = str(meta["interleave"]).lower()
	dtype_map = {
		"1": np.uint8,
		"2": np.int16,
		"3": np.int32,
		"4": np.float32,
		"12": np.uint16,
		"13": np.uint32,
	}
	dt = dtype_map.get(str(meta["data type"]).strip(), np.float32)

	arr = np.fromfile(spe_path, dtype=dt)
	expected = lines * samples * bands
	if arr.size < expected:
		raise ValueError(f"SPE size too small: {arr.size} < {expected}")
	if arr.size > expected:
		arr = arr[:expected]

	if interleave == "bsq":
		arr = arr.reshape((bands, lines, samples)).transpose(1, 2, 0)
	elif interleave == "bil":
		arr = arr.reshape((lines, bands, samples)).transpose(0, 2, 1)
	elif interleave == "bip":
		arr = arr.reshape((lines, samples, bands))
	else:
		raise ValueError(f"Unknown interleave: {interleave}")
	return arr.astype(np.float32, copy=False)


def _candidate_dirs(data_root: str, prefix: str, global_idx: Optional[int]) -> List[str]:
	"""Build ordered candidate directories for given prefix (e.g., '14') and global index.
	Handles: 14H, 14h, 14 H, 14H-2, 14h-2, 14H2
	"""
	cands: List[str] = []
	# first include all dirs starting with prefix (robust):
	try:
		for d in os.listdir(data_root):
			if d.lower().startswith(prefix.lower()):
				cands.append(os.path.join(data_root, d))
	except FileNotFoundError:
		return []

	# deterministic ones based on chunk
	if global_idx is not None:
		chunk = (global_idx - 1) // CHUNK_SIZE + 1
		if chunk == 1:
			cands.insert(0, os.path.join(data_root, f"{prefix}H"))
			cands.insert(1, os.path.join(data_root, f"{prefix}h"))
		else:
			cands.insert(0, os.path.join(data_root, f"{prefix}H-{chunk}"))
			cands.insert(1, os.path.join(data_root, f"{prefix}h-{chunk}"))
			cands.insert(2, os.path.join(data_root, f"{prefix}H{chunk}"))
	# remove duplicates, keep order
	uniq = []
	seen = set()
	for p in cands:
		if p not in seen:
			uniq.append(p)
			seen.add(p)
	return uniq


def find_envi_for_sample(data_root: str, sample_name: str) -> Optional[Tuple[str, str, str]]:
	"""Return (folder_rel, hdr_path, spe_path) if found, else None.
	folder_rel is a relative folder like '14H' or '46H-2/46-1'. Final .npy path will be folder_rel + '.npy'
	We compute local index inside chunked folder when sample_name has prefix-suffix pattern like '46-41'.
	"""
	name = sample_name.strip().replace("\\", "/")
	# extract prefix (before '-') and suffix (after '-') if possible
	prefix = None
	suffix_num: Optional[int] = None
	if "-" in name:
		parts = name.split("-", 1)
		prefix = ''.join([ch for ch in parts[0] if ch.isdigit()]) or parts[0]
		try:
			suffix_num = int(''.join([ch for ch in parts[1] if ch.isdigit()]))
		except Exception:
			suffix_num = None
	else:
		# try pattern like '14H/14-1' -> take the last segment
		base = os.path.basename(name)
		if "-" in base:
			return find_envi_for_sample(data_root, base)
		# fallback: take leading digits as prefix
		digits = ''.join([ch for ch in base if ch.isdigit()])
		prefix = digits if digits else base

	candidate_dirs = _candidate_dirs(data_root, prefix, suffix_num)
	local_idx = ((suffix_num - 1) % CHUNK_SIZE + 1) if suffix_num is not None else None

	# search directories
	for d in candidate_dirs:
		if not os.path.isdir(d):
			continue
		# try direct files in dir
		direct_hdr = glob.glob(os.path.join(d, "*.hdr"))
		direct_spe = glob.glob(os.path.join(d, "*.spe"))
		# if there are direct files, try to pick ones containing the local/global index
		if direct_hdr and direct_spe:
			def _pick_match(paths: List[str]) -> Optional[str]:
				# prefer containing '-<local_idx>' or full sample_name
				lc = str(local_idx) if local_idx is not None else None
				sn = name.lower()
				for p in paths:
					b = os.path.basename(p).lower()
					if lc and (f"-{lc}." in b or f"-{lc}_" in b or b.startswith(f"{prefix.lower()}-{lc}")):
						return p
					if prefix and b.startswith(prefix.lower()) and (lc is None or lc in b):
						return p
					if sn in b:
						return p
				return None
			hdr = _pick_match(direct_hdr) or direct_hdr[0]
			spe = _pick_match(direct_spe) or direct_spe[0]
			if hdr and spe:
				# folder_rel is d relative to data_root
				folder_rel = os.path.relpath(d, data_root)
				return folder_rel, hdr, spe

		# try subfolders like '46-1', '46-2', ... and find files inside
		try:
			for sub in os.listdir(d):
				subp = os.path.join(d, sub)
				if not os.path.isdir(subp):
					continue
				# match subfolder by local index
				if local_idx is not None:
					if not (sub.lower() == f"{prefix.lower()}-{local_idx}" or sub.lower() == f"{prefix.lower()}{local_idx}" or sub.lower().startswith(f"{prefix.lower()}-{local_idx}")):
						# not a direct match; still allow scanning but deprioritized
						pass
				hdrs = glob.glob(os.path.join(subp, "*.hdr"))
				spes = glob.glob(os.path.join(subp, "*.spe"))
				if hdrs and spes:
					# pick first or matching
					hdr = hdrs[0]
					spe = spes[0]
					folder_rel = os.path.relpath(subp, data_root)
					return folder_rel, hdr, spe
		except Exception:
			pass

	# Final fallback: recursive glob anywhere containing digits from name
	nums = [tok for tok in name.replace('-', ' ').split() if any(c.isdigit() for c in tok)]
	hdrs = glob.glob(os.path.join(data_root, "**", "*.hdr"), recursive=True)
	spes = glob.glob(os.path.join(data_root, "**", "*.spe"), recursive=True)
	def match_any(paths: List[str]) -> Optional[str]:
		for p in paths:
			b = os.path.basename(p).lower()
			if any(n.lower() in b for n in nums):
				return p
		return None
	hdr = match_any(hdrs)
	spe = match_any(spes)
	if hdr and spe:
		folder_rel = os.path.relpath(os.path.dirname(hdr), data_root)
		return folder_rel, hdr, spe

	return None


def canonical_rel_path(sample_name: str) -> Optional[str]:
	"""Compute canonical relative dataset path like '6H-2/6-1' from '6-41'.
	Returns None if cannot parse numbers.
	"""
	name = sample_name.strip()
	if "-" not in name:
		return None
	p0, p1 = name.split("-", 1)
	prefix_digits = ''.join([ch for ch in p0 if ch.isdigit()]) or p0
	try:
		global_idx = int(''.join([ch for ch in p1 if ch.isdigit()]))
	except Exception:
		return None
	chunk = (global_idx - 1) // CHUNK_SIZE + 1
	local_idx = ((global_idx - 1) % CHUNK_SIZE) + 1
	head = f"{prefix_digits}H" if chunk == 1 else f"{prefix_digits}H-{chunk}"
	leaf = f"{prefix_digits}-{local_idx}"
	return os.path.join(head, leaf)


def compute_stats_over_saved(npy_paths: List[str]) -> Tuple[List[float], List[float]]:
	"""Compute per-channel mean/std across a list of saved .npy files.
	Assumes all arrays share the same last-dimension (bands). If shapes differ in H/W,
	aggregates across all pixels.
	"""
	if not npy_paths:
		return [], []
	# determine bands from first file
	first = np.load(npy_paths[0], mmap_mode='r')
	if first.ndim == 2:
		# upgrade grayscale to one-band
		bands = 1
	else:
		bands = int(first.shape[2])
	total = np.zeros((bands,), dtype=np.float64)
	total_sq = np.zeros((bands,), dtype=np.float64)
	n_pix = 0

	for p in npy_paths:
		arr = np.load(p)
		if arr.ndim == 2:
			arr = arr[:, :, None]
		h, w, b = arr.shape
		if b != bands:
			# skip inconsistent band count
			continue
		flat = arr.reshape(-1, b).astype(np.float64)
		total += flat.sum(axis=0)
		total_sq += (flat ** 2).sum(axis=0)
		n_pix += flat.shape[0]

	if n_pix == 0:
		return [], []
	mean = (total / n_pix)
	var = total_sq / n_pix - mean ** 2
	std = np.sqrt(np.maximum(var, 1e-12))
	return mean.tolist(), std.tolist()


def main():
	ap = argparse.ArgumentParser(description="Prepare .npy and metadata from ENVI dataset")
	ap.add_argument("--excel", required=True, help="Path to the source Excel, e.g., 20251015-数据统计-是这个.xlsx")
	ap.add_argument("--data_root", required=True, help="Root folder of hyperspectral ENVI data (contains *H, *H-2 folders)")
	ap.add_argument("--out_dir", required=True, help="Output folder to write .npy files and (optionally) stats.json")
	ap.add_argument("--metadata", required=True, help="Path to write new_metadata.csv")
	ap.add_argument("--recursive", action="store_true", help="Unused flag (kept for CLI compatibility)")
	ap.add_argument("--limit", type=int, default=0, help="Process only the first N rows for a quick dry run")
	ap.add_argument("--compute-stats", action="store_true", help="Also compute mean/std over saved .npy and write stats.json")
	ap.add_argument("--resume", action="store_true", help="Skip samples whose target .npy already exists; merge with existing metadata if present")
	ap.add_argument("--start-index", type=int, default=1, help="1-based row index in Excel (after cleaning) to start from, e.g., 213")
	args = ap.parse_args()

	excel = os.path.normpath(args.excel)
	data_root = os.path.normpath(args.data_root)
	out_dir = os.path.normpath(args.out_dir)
	metadata_csv = os.path.normpath(args.metadata)
	ensure_dir(out_dir)
	ensure_dir(os.path.dirname(metadata_csv))

	log(f"Excel: {excel}")
	log(f"Data root: {data_root}")
	log(f"Out dir: {out_dir}")
	log(f"Metadata CSV: {metadata_csv}")

	df = parse_excel(excel)
	# Apply start-index (1-based)
	start_idx = max(1, int(args.start_index)) if args.start_index else 1
	if start_idx > 1:
		if start_idx > len(df):
			log(f"start-index {start_idx} beyond dataset length {len(df)}; nothing to do")
			# still flush existing to ensure CSV header exists
			# and exit early
			# create empty CSV if needed
			existing_only = []
			keys = ["filename", "npy_filename", "Value_Raw", "toxin_value", "sample_name"]
			if not os.path.isfile(metadata_csv):
				with open(metadata_csv, "w", newline="", encoding="utf-8-sig") as fw:
					writer = csv.DictWriter(fw, fieldnames=keys)
					writer.writeheader()
			return
		log(f"Starting from row {start_idx}; skipping first {start_idx-1} rows")
		df = df.iloc[start_idx - 1 : ].copy()
	if args.limit and args.limit > 0:
		df = df.iloc[: args.limit].copy()
		log(f"Limiting to first {len(df)} rows for dry run")

	# Load existing metadata (for resume) if any
	existing_records: List[Dict[str, str]] = []
	existing_npy: set = set()
	if os.path.isfile(metadata_csv):
		try:
			with open(metadata_csv, "r", encoding="utf-8-sig") as fr:
				reader = csv.DictReader(fr)
				for r in reader:
					existing_records.append(r)
					if "npy_filename" in r and r["npy_filename"]:
						existing_npy.add(r["npy_filename"].replace("\\", "/"))
			log(f"Loaded existing metadata rows: {len(existing_records)}")
		except Exception as e:
			log(f"Warning: failed to load existing metadata for resume: {e}")

	records = []  # new records in this run
	saved_paths = []  # for stats
	errors = []
	flush_interval = 100

	def flush_metadata():
		# Merge existing + new and write de-duplicated by npy_filename
		if not (existing_records or records):
			return
		merged = []
		seen = set()
		for r in existing_records + records:
			key = (r.get("npy_filename") or r.get("filename") or "").replace("\\", "/")
			if key and key not in seen:
				merged.append(r)
				seen.add(key)
		keys = ["filename", "npy_filename", "Value_Raw", "toxin_value", "sample_name"]
		try:
			with open(metadata_csv, "w", newline="", encoding="utf-8-sig") as fw:
				writer = csv.DictWriter(fw, fieldnames=keys)
				writer.writeheader()
				writer.writerows(merged)
			log(f"Flushed metadata: {len(merged)} rows")
		except Exception as e:
			log(f"Warning: failed to flush metadata: {e}")

	for i, row in df.iterrows():
		sample_name = str(row["sample_name"]).strip()
		value_raw = float(row["Value_Raw"]) if pd.notna(row["Value_Raw"]) else None
		toxin_value = float(row["toxin_value"]) if pd.notna(row["toxin_value"]) else None
		if value_raw is None or toxin_value is None:
			errors.append((sample_name, "missing_value", None))
			continue

		# We first compute the canonical path for consistent metadata & npy placement
		canon_rel = canonical_rel_path(sample_name)
		# If canonical path is computable and npy already exists (resume), rebuild metadata row only
		if canon_rel is not None:
			npy_rel = f"{canon_rel}.npy"
			dataset_rel = canon_rel
			out_path = os.path.join(out_dir, npy_rel)
			if args.resume and os.path.isfile(out_path):
				records.append({
					"filename": os.path.join(data_root, dataset_rel).replace("\\", "/"),
					"npy_filename": npy_rel.replace("\\", "/"),
					"Value_Raw": value_raw,
					"toxin_value": toxin_value,
					"sample_name": sample_name,
				})
				if (i + 1) % 20 == 0:
					log(f"Processed {i + 1} / {len(df)}")
				if (len(records) % flush_interval) == 0 and len(records) > 0:
					flush_metadata()
				continue

		# We still locate actual hdr/spe to read data
		found = find_envi_for_sample(data_root, sample_name)
		if not found:
			errors.append((sample_name, "not_found", None))
			continue
		folder_rel, hdr_p, spe_p = found

		try:
			meta = parse_hdr(hdr_p)
			arr = read_spe(spe_p, meta)  # (H,W,B), float32
		except Exception as e:
			errors.append((sample_name, "read_error", str(e)))
			continue

		# Compose canonical relative filename for npy: e.g., '6H-2/6-1.npy'
		if canon_rel is None:
			# fallback to discovered folder and stem
			leaf = os.path.basename(folder_rel)
			base_name = leaf if ('-' in leaf and any(c.isdigit() for c in leaf)) else os.path.splitext(os.path.basename(spe_p))[0]
			npy_rel = os.path.join(folder_rel, f"{base_name}.npy")
			dataset_rel = os.path.join(folder_rel, base_name)
		else:
			npy_rel = f"{canon_rel}.npy"
			dataset_rel = canon_rel
		out_path = os.path.join(out_dir, npy_rel)

		# Resume: if target npy exists and --resume, skip recompute and just ensure metadata will have the row
		if args.resume and os.path.isfile(out_path):
			pass  # we will still add metadata row if missing below
		else:
			# Make directory and save, with ENOSPC handling
			try:
				ensure_dir(os.path.dirname(out_path))
				np.save(out_path, arr.astype(np.float32, copy=False))
				saved_paths.append(out_path)
			except OSError as e:
				if isinstance(e, OSError) and e.errno == errno.ENOSPC:
					log(f"Disk full when saving {out_path}. Stopping further writes; you can free space and rerun with --resume.")
					# Flush metadata collected so far, then break out of processing loop
					flush_metadata()
					break
				else:
					errors.append((sample_name, "save_error", str(e)))
					continue

		# Store metadata:
		# - filename: absolute canonical dataset path under data_root (folder), e.g., '/.../6H-2/6-1'
		# - npy_filename: path of saved npy relative to out_dir, e.g., '6H-2/6-1.npy'
		records.append({
			"filename": os.path.join(data_root, dataset_rel).replace("\\", "/"),
			"npy_filename": npy_rel.replace("\\", "/"),
			"Value_Raw": value_raw,
			"toxin_value": toxin_value,
			"sample_name": sample_name,
		})

		if (i + 1) % 20 == 0:
			log(f"Processed {i + 1} / {len(df)}")
		# periodic flush to not lose progress
		if (len(records) % flush_interval) == 0 and len(records) > 0:
			flush_metadata()

	# Write metadata CSV
	# Final flush: merge and write all (existing + new)
	flush_metadata()

	# Write errors log if any
	if errors:
		err_path = os.path.join(out_dir, "prepare_errors.csv")
		with open(err_path, "w", newline="", encoding="utf-8-sig") as fw:
			w = csv.writer(fw)
			w.writerow(["sample_name", "error_type", "detail"])
			w.writerows(errors)
		log(f"Logged {len(errors)} issues -> {err_path}")

	# Optional stats
	if args.compute_stats and records:
		mean, std = compute_stats_over_saved([os.path.join(out_dir, r["npy_filename"]) for r in records])
		stats = {"mean": mean, "std": std}
		stats_path = os.path.join(out_dir, "stats.json")
		with open(stats_path, "w", encoding="utf-8") as fw:
			json.dump(stats, fw)
		log(f"Wrote stats -> {stats_path}")

	log("Done.")


if __name__ == "__main__":
	try:
		main()
	except KeyboardInterrupt:
		log("Interrupted.")
		sys.exit(130)
	except Exception as e:
		log(f"Fatal error: {e}")
		sys.exit(1)