{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"702acf32-bab0-46c0-bc2c-d3fe074ebd1a","cell_type":"markdown","source":"# Inversion de Forme d'Onde Sismique par Deep Learning\n## Full Waveform Inversion (FWI) — ResUNet + CBAM Attention\n\n---\n\n### Contexte du problème\n\nL'**inversion de forme d'onde sismique** (FWI) est l'un des problèmes les plus complexes de la géophysique appliquée. L'objectif est de reconstruire une **carte de vitesse de propagation des ondes** dans le sous-sol à partir de signaux sismiques enregistrés en surface.\n\nConcrètement :\n- **Input** : 5 sismogrammes (récepteurs) enregistrant les réflexions d'ondes sismiques — shape `(5, 1000, 70)`\n- **Output** : carte 2D de vitesse en m/s représentant la géologie — shape `(1, 70, 70)`\n\nLes vitesses varient de **1500 m/s** (sédiments meubles) à **4500 m/s** (roches dures), et chaque interface visible correspond à un changement lithologique réel.\n\n### Notre approche\n\nNous utilisons un **ResUNet** (U-Net avec blocs résiduels) couplé à des modules d'attention **CBAM** (Convolutional Block Attention Module), entraîné en précision float32 complète pour une stabilité maximale.\n\n| Composante | Choix | Justification |\n|---|---|---|\n| Architecture | ResUNet + CBAM | Skip connections + attention spatiale/canal |\n| Paramètres | 33.5M | Capacité suffisante pour les structures latérales |\n| Loss | MSE + L1 + SSIM + Gradient | Précision globale + robustesse + structure |\n| Optimiseur | AdamW + Warmup + Cosine | Convergence stable sur 99 epochs |\n| Précision | Float32 (pas d'AMP) | Évite les instabilités numériques de float16 |","metadata":{}},{"id":"682fe1f8-1e07-450b-8085-2dcf43e2b091","cell_type":"code","source":"# ============================================================\n#  INVERSION DE FORME D'ONDE SISMIQUE\n#  Approche deep learning : ResUNet + CBAM\n#  Float32 pur — stabilite maximale\n# ============================================================\nimport os, time, warnings\nimport numpy as np\nimport pandas as pd\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader\nimport torch.optim.lr_scheduler as lr_scheduler\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nimport matplotlib.patches as mpatches\nfrom matplotlib.colors import TwoSlopeNorm\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error\nfrom scipy import stats\nwarnings.filterwarnings('ignore')\n\nSEED = 42\ntorch.manual_seed(SEED)\nnp.random.seed(SEED)\nif torch.cuda.is_available():\n    torch.cuda.manual_seed_all(SEED)\n\nprint(f'PyTorch : {torch.__version__}')\nprint(f'CUDA    : {torch.cuda.is_available()}')\nif torch.cuda.is_available():\n    print(f'GPU     : {torch.cuda.get_device_name(0)}')\n    print(f'VRAM    : {torch.cuda.get_device_properties(0).total_memory/1e9:.1f} GB')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:20:01.406320Z","iopub.execute_input":"2026-06-02T16:20:01.406585Z","iopub.status.idle":"2026-06-02T16:20:07.684761Z","shell.execute_reply.started":"2026-06-02T16:20:01.406549Z","shell.execute_reply":"2026-06-02T16:20:07.683993Z"}},"outputs":[],"execution_count":null},{"id":"5b4d4e93-7924-48d9-b891-7a8ca9712b3a","cell_type":"code","source":"# ============================================================\n#  CONFIGURATION CENTRALE\n# ============================================================\nconfig = {\n    'batch_size': 8,\n    'validation_split': 0.2,\n    'num_workers': 4,\n    'pin_memory': True,\n    'learning_rate': 5e-5,\n    'epochs': 100,\n    'weight_decay': 1e-4,\n    'gradient_clip': 1.0,\n    'patience': 15,\n    'n_channels': 5,\n    'n_classes': 1,\n    'base_filters': 64,\n    'use_amp': False,\n    'device': 'cuda' if torch.cuda.is_available() else 'cpu',\n    'model_path': 'best_resunet_seismic.pth',\n    'clip_percentile': 1,\n}\nprint(f'Device      : {config[\"device\"]}')\nprint(f'LR          : {config[\"learning_rate\"]}')\nprint(f'Base filters: {config[\"base_filters\"]}')\nprint(f'AMP         : {config[\"use_amp\"]}')\nprint(f'Batch size  : {config[\"batch_size\"]}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:20:07.686165Z","iopub.execute_input":"2026-06-02T16:20:07.686629Z","iopub.status.idle":"2026-06-02T16:20:07.693041Z","shell.execute_reply.started":"2026-06-02T16:20:07.686603Z","shell.execute_reply":"2026-06-02T16:20:07.692204Z"}},"outputs":[],"execution_count":null},{"id":"4ab98de1-5f41-4fd9-b4cc-2984d10246ca","cell_type":"markdown","source":"## Section 1 — Chargement des données\n\nLe dataset provient de la compétition Kaggle **Waveform Inversion**.  \nLes données sont organisées en paires fichiers `.npy` :\n- `data_*.npy` : sismogrammes (input X) — shape `(500, 5, 1000, 70)` par fichier\n- `model_*.npy` : cartes de vitesse (target Y) — shape `(500, 1, 70, 70)` par fichier\n\nLe chargement concatène tous les fichiers pour obtenir le dataset complet.","metadata":{}},{"id":"6d3040cf-68f9-4dd1-9d26-d8e84de23694","cell_type":"code","source":"# ============================================================\n#  SECTION 1 — CHARGEMENT DES DONNEES\n# ============================================================\nbase_dir = '/kaggle/input/competitions/waveform-inversion/train_samples'\nseismic_files, velocity_files = [], []\n\nfor root, dirs, files in os.walk(base_dir):\n    for f in sorted(files):\n        if f.endswith('.npy'):\n            full_path = os.path.join(root, f)\n            if 'data' in full_path:\n                seismic_files.append(full_path)\n            elif 'model' in full_path:\n                velocity_files.append(full_path)\n\nseismic_files  = sorted(seismic_files)\nvelocity_files = sorted(velocity_files)\nassert len(seismic_files) == len(velocity_files)\n\nprint(f'Fichiers sismiques : {len(seismic_files)}')\nprint(f'Fichiers vitesse   : {len(velocity_files)}')\n\nt0 = time.time()\nx_raw = np.concatenate([np.load(f) for f in seismic_files], axis=0)\ny_raw = np.concatenate([np.load(f) for f in velocity_files], axis=0)\nprint(f'Chargement : {time.time()-t0:.1f}s')\nprint(f'Shape X : {x_raw.shape}')\nprint(f'Shape Y : {y_raw.shape}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:20:07.693889Z","iopub.execute_input":"2026-06-02T16:20:07.694173Z","iopub.status.idle":"2026-06-02T16:21:10.304843Z","shell.execute_reply.started":"2026-06-02T16:20:07.694151Z","shell.execute_reply":"2026-06-02T16:21:10.303999Z"}},"outputs":[],"execution_count":null},{"id":"39c126ea-228b-4483-b71b-7e42577d1c7c","cell_type":"markdown","source":"## Section 2 — Exploration des Données (EDA)\n\nAvant tout entraînement, nous analysons visuellement et statistiquement les données pour comprendre leur structure, leur distribution et leurs caractéristiques physiques. Cette étape est essentielle pour justifier les choix de preprocessing et d'architecture.","metadata":{}},{"id":"a234ef38-72af-48ba-bba6-431129c98d06","cell_type":"code","source":"# ============================================================\n#  SECTION 2 — EXPLORATION DES DONNEES (EDA)\n# ============================================================\nprint('Statistiques X :')\nprint(f'  Min={x_raw.min():.4f} Max={x_raw.max():.4f} Mean={x_raw.mean():.4f} Std={x_raw.std():.4f}')\nprint('Statistiques Y (m/s) :')\nprint(f'  Min={y_raw.min():.2f} Max={y_raw.max():.2f} Mean={y_raw.mean():.2f} Std={y_raw.std():.2f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:10.305935Z","iopub.execute_input":"2026-06-02T16:21:10.306375Z","iopub.status.idle":"2026-06-02T16:21:17.407715Z","shell.execute_reply.started":"2026-06-02T16:21:10.306338Z","shell.execute_reply":"2026-06-02T16:21:17.406922Z"}},"outputs":[],"execution_count":null},{"id":"fc13aba8-45f2-4513-b4ea-bba47da2abce","cell_type":"code","source":"# EDA 2.1 — Sismogrammes bruts\nfig, axes = plt.subplots(1, 5, figsize=(20, 5))\nfig.suptitle('Sismogrammes bruts — 5 recepteurs (echantillon #0)', fontsize=13, fontweight='bold')\nfor ch in range(5):\n    im = axes[ch].imshow(x_raw[0, ch, :, :], aspect='auto', cmap='RdBu_r')\n    axes[ch].set_title(f'Recepteur {ch+1}', fontsize=11)\n    axes[ch].set_xlabel('Position (x)', fontsize=9)\n    if ch == 0: axes[ch].set_ylabel('Temps (samples)', fontsize=9)\n    plt.colorbar(im, ax=axes[ch], fraction=0.046, pad=0.04)\nplt.tight_layout()\nplt.savefig('eda_sismogrammes.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:17.409780Z","iopub.execute_input":"2026-06-02T16:21:17.410133Z","iopub.status.idle":"2026-06-02T16:21:19.383690Z","shell.execute_reply.started":"2026-06-02T16:21:17.410106Z","shell.execute_reply":"2026-06-02T16:21:19.382874Z"}},"outputs":[],"execution_count":null},{"id":"44cef883-e0a8-4469-90a3-ece0653f0d08","cell_type":"markdown","source":"**Lecture des sismogrammes** : chaque panneau représente un récepteur différent. On observe clairement les **arrivées hyperboliques** caractéristiques de la propagation d'ondes dans un milieu stratifié. Le récepteur central (R3) présente le signal le plus symétrique — les récepteurs latéraux enregistrent les ondes avec un décalage temporel croissant selon leur distance à la source.","metadata":{}},{"id":"82208400-cc28-485e-a73c-2b96a62a8d6c","cell_type":"code","source":"# EDA 2.2 — Cartes de vitesse\nfig, axes = plt.subplots(1, 3, figsize=(15, 4))\nfig.suptitle('Cartes de vitesse (target Y) — 3 echantillons\\nLes couches = interfaces geologiques', fontsize=13, fontweight='bold')\nfor i, idx in enumerate([0, 5, 10]):\n    im = axes[i].imshow(y_raw[idx, 0], cmap='jet', aspect='auto')\n    axes[i].set_title(f'Echantillon {idx+1}', fontsize=11)\n    axes[i].set_xlabel('Position X', fontsize=9)\n    if i == 0: axes[i].set_ylabel('Profondeur Z', fontsize=9)\n    cbar = plt.colorbar(im, ax=axes[i], fraction=0.046, pad=0.04)\n    cbar.set_label('Vitesse (m/s)', fontsize=8)\nplt.tight_layout()\nplt.savefig('eda_cartes_vitesse.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:19.384598Z","iopub.execute_input":"2026-06-02T16:21:19.384889Z","iopub.status.idle":"2026-06-02T16:21:20.571255Z","shell.execute_reply.started":"2026-06-02T16:21:19.384865Z","shell.execute_reply":"2026-06-02T16:21:20.570342Z"}},"outputs":[],"execution_count":null},{"id":"813ae02c-eb60-4bde-9f95-8d5d82796887","cell_type":"markdown","source":"**Cartes de vitesse** : les interfaces géologiques apparaissent comme des **transitions de couleur nettes**. La vitesse augmente avec la profondeur (bleu → rouge), ce qui reflète la compaction naturelle des sédiments. Les variations latérales (ondulations des couches) représentent la complexité structurale du sous-sol — c'est précisément ce que le modèle doit apprendre à reconstruire.","metadata":{}},{"id":"ecec8da4-fa22-47c4-8383-5d79a0913907","cell_type":"code","source":"# EDA 2.3 — Distribution des vitesses\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\nfig.suptitle('Distribution des vitesses sismiques', fontsize=13, fontweight='bold')\naxes[0].hist(y_raw.flatten(), bins=80, color='steelblue', edgecolor='white', linewidth=0.3, alpha=0.85)\naxes[0].axvline(y_raw.mean(), color='red', linestyle='--', linewidth=2, label=f'Moyenne : {y_raw.mean():.0f} m/s')\naxes[0].axvline(np.percentile(y_raw, 5),  color='orange', linestyle=':', linewidth=1.5, label=f'P5  : {np.percentile(y_raw,5):.0f} m/s')\naxes[0].axvline(np.percentile(y_raw, 95), color='orange', linestyle=':', linewidth=1.5, label=f'P95 : {np.percentile(y_raw,95):.0f} m/s')\naxes[0].set_xlabel('Vitesse (m/s)', fontsize=11)\naxes[0].set_ylabel('Frequence', fontsize=11)\naxes[0].set_title('Distribution globale', fontsize=11)\naxes[0].legend(fontsize=9)\naxes[0].grid(True, alpha=0.3)\nn_depth = y_raw.shape[2]\ndepth_levels = np.linspace(0, n_depth-1, 10, dtype=int)\ndepth_data = [y_raw[:, 0, d, :].flatten() for d in depth_levels]\nbp = axes[1].boxplot(depth_data, patch_artist=True, labels=[f'z={d}' for d in depth_levels])\nfor patch in bp['boxes']:\n    patch.set_facecolor('lightsteelblue')\n    patch.set_alpha(0.8)\naxes[1].set_xlabel('Niveau de profondeur', fontsize=11)\naxes[1].set_ylabel('Vitesse (m/s)', fontsize=11)\naxes[1].set_title('Vitesses par profondeur', fontsize=10)\naxes[1].grid(True, alpha=0.3)\naxes[1].tick_params(axis='x', rotation=45)\nplt.tight_layout()\nplt.savefig('eda_distribution_vitesses.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:20.572347Z","iopub.execute_input":"2026-06-02T16:21:20.572733Z","iopub.status.idle":"2026-06-02T16:21:23.177759Z","shell.execute_reply.started":"2026-06-02T16:21:20.572708Z","shell.execute_reply":"2026-06-02T16:21:23.176967Z"}},"outputs":[],"execution_count":null},{"id":"c9594f30-5214-444d-8b8a-a06ddb301a18","cell_type":"markdown","source":"**Distribution des vitesses** : la distribution est quasi-uniforme entre 1500 et 4500 m/s, avec une moyenne à **2850 m/s**. Le boxplot par profondeur confirme que les vitesses augmentent systématiquement avec la profondeur (médiane de ~1800 m/s en surface contre ~3500 m/s en profondeur), ce qui est cohérent avec la physique de la compaction sédimentaire.","metadata":{}},{"id":"9cf70a97-1ce4-4e94-b1bd-cbf985f1b19e","cell_type":"code","source":"# EDA 2.4 — Correlation inter-canaux\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\nfig.suptitle('Analyse de correlation entre recepteurs', fontsize=13, fontweight='bold')\nsample_flat = x_raw[0].reshape(5, -1)\ncorr_matrix = np.corrcoef(sample_flat)\nim = axes[0].imshow(corr_matrix, cmap='coolwarm', vmin=-1, vmax=1)\naxes[0].set_xticks(range(5))\naxes[0].set_yticks(range(5))\naxes[0].set_xticklabels([f'R{i+1}' for i in range(5)])\naxes[0].set_yticklabels([f'R{i+1}' for i in range(5)])\naxes[0].set_title('Correlation inter-recepteurs', fontsize=11)\nplt.colorbar(im, ax=axes[0], fraction=0.046)\nfor i in range(5):\n    for j in range(5):\n        axes[0].text(j, i, f'{corr_matrix[i,j]:.2f}', ha='center', va='center',\n                     fontsize=9, color='white' if abs(corr_matrix[i,j]) > 0.5 else 'black')\ncolors = plt.cm.viridis(np.linspace(0, 1, 5))\nfor ch in range(5):\n    signal = x_raw[0, ch, :, 35]\n    axes[1].plot(np.arange(x_raw.shape[2])[::10], signal[::10], color=colors[ch], alpha=0.8, label=f'R{ch+1}')\naxes[1].set_xlabel('Temps (samples)', fontsize=11)\naxes[1].set_ylabel('Amplitude', fontsize=11)\naxes[1].set_title('Signaux temporels — colonne centrale', fontsize=10)\naxes[1].legend(fontsize=9)\naxes[1].grid(True, alpha=0.3)\nplt.tight_layout()\nplt.savefig('eda_correlation.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:23.178658Z","iopub.execute_input":"2026-06-02T16:21:23.178971Z","iopub.status.idle":"2026-06-02T16:21:23.975206Z","shell.execute_reply.started":"2026-06-02T16:21:23.178935Z","shell.execute_reply":"2026-06-02T16:21:23.974571Z"}},"outputs":[],"execution_count":null},{"id":"a5d528de-cdea-414c-9559-4f9920ef7740","cell_type":"markdown","source":"**Corrélation inter-récepteurs** : la matrice de corrélation montre que les 5 récepteurs sont **quasi-indépendants** (corrélations ≈ 0). Chaque récepteur apporte donc une information complémentaire sur la structure du sous-sol, ce qui justifie l'utilisation des 5 canaux comme entrée du modèle. Les signaux temporels montrent des impulsions brèves — la majorité de l'énergie est concentrée dans les premiers 300 samples.","metadata":{}},{"id":"304cea48-d849-4f8c-a6f4-ddf5ef5ee130","cell_type":"code","source":"# EDA 2.5 — Paires Input->Output\nfig = plt.figure(figsize=(20, 8))\nfig.suptitle('La tache : Sismogrammes -> Carte de vitesse\\nLe modele apprend a voir a travers les ondes',\n             fontsize=13, fontweight='bold')\ngs = gridspec.GridSpec(2, 4, figure=fig, hspace=0.4, wspace=0.3)\nfor i in range(4):\n    ax_x = fig.add_subplot(gs[0, i])\n    im = ax_x.imshow(x_raw[i, 0, ::10, :], aspect='auto', cmap='RdBu_r')\n    ax_x.set_title(f'Input #{i+1} (Sismo R1)', fontsize=9)\n    ax_x.set_xlabel('Position', fontsize=8)\n    if i == 0: ax_x.set_ylabel('Temps', fontsize=8)\n    plt.colorbar(im, ax=ax_x, fraction=0.046, pad=0.04)\n    ax_y = fig.add_subplot(gs[1, i])\n    im2 = ax_y.imshow(y_raw[i, 0], cmap='jet', aspect='auto')\n    ax_y.set_title(f'Target #{i+1} (Vitesse)', fontsize=9)\n    ax_y.set_xlabel('Position X', fontsize=8)\n    if i == 0: ax_y.set_ylabel('Profondeur', fontsize=8)\n    plt.colorbar(im2, ax=ax_y, fraction=0.046, pad=0.04)\nplt.savefig('eda_paires.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:23.976269Z","iopub.execute_input":"2026-06-02T16:21:23.976585Z","iopub.status.idle":"2026-06-02T16:21:26.178671Z","shell.execute_reply.started":"2026-06-02T16:21:23.976562Z","shell.execute_reply":"2026-06-02T16:21:26.177839Z"}},"outputs":[],"execution_count":null},{"id":"a0882485-cded-435e-a8e7-ea050676b259","cell_type":"markdown","source":"**La tâche visualisée** : la difficulté du problème est immédiatement visible. Le sismogramme (input) est un signal 1D projeté en 2D — il ne ressemble en rien visuellement à la carte de vitesse (output). Le modèle doit apprendre une transformation hautement non-linéaire, en \"inversant\" la physique de propagation des ondes pour remonter à la géologie source.","metadata":{}},{"id":"a6066c07-31b2-49e2-8475-b606381dfc59","cell_type":"markdown","source":"## Section 3 — Preprocessing\n\nDeux opérations sont appliquées séquentiellement :\n\n1. **Downsampling** : les sismogrammes `(N, 5, 1000, 70)` sont réduits à `(N, 5, 70, 70)` par interpolation bilinéaire. La dimension temporelle (1000 → 70) compresse l'information temporelle tout en préservant les structures principales.\n\n2. **Normalisation robuste** : les statistiques (mean, std) sont calculées **uniquement sur le train set** pour éviter toute fuite de données vers la validation. La normalisation par percentile (P1/P99) élimine les outliers d'acquisition avant de centrer-réduire.","metadata":{}},{"id":"258e61ef-ee46-48fd-a76a-70f997c0f680","cell_type":"code","source":"# ============================================================\n#  SECTION 3 — PREPROCESSING\n# ============================================================\ndef downsample_seismic(arr, target_size=(70, 70)):\n    tensor = torch.tensor(arr, dtype=torch.float32)\n    tensor = tensor.permute(0, 1, 3, 2)\n    tensor = F.interpolate(tensor, size=target_size, mode='bilinear', align_corners=False)\n    tensor = tensor.permute(0, 1, 3, 2)\n    return tensor.numpy()\n\ndef robust_normalize(data, clip_percentile=1):\n    lower = np.percentile(data, clip_percentile)\n    upper = np.percentile(data, 100 - clip_percentile)\n    data_clipped = np.clip(data, lower, upper)\n    mean = np.mean(data_clipped)\n    std  = np.std(data_clipped) + 1e-8\n    return (data_clipped - mean) / std, mean, std\n\ndef apply_normalize(data, mean, std, clip_percentile=1):\n    lower = np.percentile(data, clip_percentile)\n    upper = np.percentile(data, 100 - clip_percentile)\n    return (np.clip(data, lower, upper) - mean) / std\n\nprint('Downsampling...')\nt0 = time.time()\nx_data = downsample_seismic(x_raw)\nprint(f'Termine en {time.time()-t0:.1f}s | Shape : {x_data.shape}')\n\nx_data_raw_copy = x_data.copy()\n\nprint('Split train/val avant normalisation...')\nn_total    = len(x_data)\nval_size   = int(config['validation_split'] * n_total)\ntrain_size = n_total - val_size\ngenerator  = torch.Generator().manual_seed(SEED)\nindices    = torch.randperm(n_total, generator=generator).tolist()\ntrain_idx  = indices[:train_size]\nval_idx    = indices[train_size:]\n\nx_train_raw = x_data[train_idx].copy()\nx_val_raw   = x_data[val_idx].copy()\n\nx_train_norm = x_train_raw.copy()\nx_val_norm   = x_val_raw.copy()\nx_norm_stats = {}\nfor i in range(x_data.shape[1]):\n    x_train_norm[:, i], m, s = robust_normalize(x_train_raw[:, i])\n    x_val_norm[:, i]          = apply_normalize(x_val_raw[:, i], m, s)\n    x_norm_stats[i]           = (m, s)\n\ny_train_raw = y_raw[train_idx]\ny_mean = float(np.mean(y_train_raw))\ny_std  = float(np.std(y_train_raw))\n\ndef normalize_y(y, mean, std):\n    return (np.clip(y, np.percentile(y, 1), np.percentile(y, 99)) - mean) / (std + 1e-8)\n\ny_train_norm = normalize_y(y_raw[train_idx], y_mean, y_std)\ny_val_norm   = normalize_y(y_raw[val_idx],   y_mean, y_std)\n\nprint(f'X train norm : mean={x_train_norm.mean():.4f} std={x_train_norm.std():.4f}')\nprint(f'X val norm   : mean={x_val_norm.mean():.4f}   std={x_val_norm.std():.4f}')\nprint(f'Y params (train) : mean={y_mean:.2f} m/s  std={y_std:.2f} m/s')\nprint(f'Train : {train_size} | Val : {val_size}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:26.179652Z","iopub.execute_input":"2026-06-02T16:21:26.179963Z","iopub.status.idle":"2026-06-02T16:21:35.901909Z","shell.execute_reply.started":"2026-06-02T16:21:26.179939Z","shell.execute_reply":"2026-06-02T16:21:35.901116Z"}},"outputs":[],"execution_count":null},{"id":"8144e3b4-fc79-4916-b4df-5d53d7118b0c","cell_type":"code","source":"# Preprocessing 3.1 — Avant / Apres normalisation\nfig, axes = plt.subplots(2, 3, figsize=(16, 8))\nfig.suptitle('Effet de la normalisation robuste sur les donnees sismiques', fontsize=13, fontweight='bold')\nfor ch_idx, ax_row in enumerate([0, 2]):\n    axes[ch_idx, 0].hist(x_data_raw_copy[:, ax_row].flatten(), bins=60, color='salmon', edgecolor='white', linewidth=0.3, alpha=0.85)\n    axes[ch_idx, 0].set_title(f'Canal {ax_row+1} — AVANT', fontsize=10)\n    axes[ch_idx, 0].set_xlabel('Amplitude brute', fontsize=9)\n    axes[ch_idx, 0].grid(True, alpha=0.3)\n    axes[ch_idx, 1].hist(x_train_norm[:, ax_row].flatten(), bins=60, color='steelblue', edgecolor='white', linewidth=0.3, alpha=0.85)\n    axes[ch_idx, 1].set_title(f'Canal {ax_row+1} — APRES', fontsize=10)\n    axes[ch_idx, 1].set_xlabel('Amplitude normalisee', fontsize=9)\n    axes[ch_idx, 1].grid(True, alpha=0.3)\n    axes[ch_idx, 2].imshow(np.hstack([x_data_raw_copy[0, ax_row], x_train_norm[0, ax_row]]), cmap='RdBu_r', aspect='auto')\n    axes[ch_idx, 2].axvline(70, color='yellow', linewidth=2, linestyle='--')\n    axes[ch_idx, 2].set_title(f'Canal {ax_row+1} : Brut | Normalise', fontsize=10)\n    axes[ch_idx, 2].set_xlabel('<- Brut  |  Normalise ->', fontsize=9)\nplt.tight_layout()\nplt.savefig('preprocessing_normalisation.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:35.902870Z","iopub.execute_input":"2026-06-02T16:21:35.903146Z","iopub.status.idle":"2026-06-02T16:21:39.403401Z","shell.execute_reply.started":"2026-06-02T16:21:35.903114Z","shell.execute_reply":"2026-06-02T16:21:39.402708Z"}},"outputs":[],"execution_count":null},{"id":"5cc109cb-dc78-4207-91a6-b4aa7bc85dba","cell_type":"markdown","source":"**Effet de la normalisation** : avant normalisation, les amplitudes sont fortement concentrées autour de zéro avec des queues épaisses (outliers d'acquisition). Après normalisation robuste, la distribution est gaussienne centrée, ce qui facilite l'optimisation du réseau. La ligne jaune sépare visuellement les deux versions du signal.","metadata":{}},{"id":"79052502-57ea-4f7f-92d4-2367d9adfd7b","cell_type":"markdown","source":"## Section 4 — Architecture : ResUNet + CBAM\n\n### Pourquoi ResUNet ?\n\nL'architecture **U-Net** est le standard de référence pour les tâches de régression pixel-à-pixel. Nous l'enrichissons avec des **blocs résiduels** (ResBlock) qui permettent :\n- Un meilleur flux de gradient lors de la rétropropagation\n- Une convergence plus rapide et plus stable\n- Une capacité à apprendre des transformations d'identité\n\n### Pourquoi CBAM ?\n\nLe module **CBAM** (Convolutional Block Attention Module) combine deux mécanismes d'attention :\n- **Channel Attention** : \"quels canaux de features sont importants ?\"\n- **Spatial Attention** : \"où regarder dans l'image ?\"\n\nCes deux modules permettent au décodeur de se concentrer sur les zones d'interfaces géologiques — là où l'information est la plus critique pour la reconstruction de vitesse.","metadata":{}},{"id":"f0e0a212-5c04-43af-9f83-2ce057a217cc","cell_type":"code","source":"# ============================================================\n#  SECTION 4 — ARCHITECTURE : ResUNet + CBAM\n# ============================================================\nclass ResBlock(nn.Module):\n    def __init__(self, in_ch, out_ch):\n        super().__init__()\n        self.conv_block = nn.Sequential(\n            nn.Conv2d(in_ch, out_ch, 3, padding=1, bias=False),\n            nn.BatchNorm2d(out_ch), nn.ReLU(inplace=True),\n            nn.Conv2d(out_ch, out_ch, 3, padding=1, bias=False),\n            nn.BatchNorm2d(out_ch),\n        )\n        self.skip = nn.Sequential(\n            nn.Conv2d(in_ch, out_ch, 1, bias=False), nn.BatchNorm2d(out_ch)\n        ) if in_ch != out_ch else nn.Identity()\n        self.relu = nn.ReLU(inplace=True)\n    def forward(self, x):\n        return self.relu(self.conv_block(x) + self.skip(x))\n\nclass ChannelAttention(nn.Module):\n    def __init__(self, channels, ratio=8):\n        super().__init__()\n        self.avg_pool = nn.AdaptiveAvgPool2d(1)\n        self.max_pool = nn.AdaptiveMaxPool2d(1)\n        mid = max(1, channels // ratio)\n        self.mlp = nn.Sequential(\n            nn.Conv2d(channels, mid, 1, bias=False), nn.ReLU(inplace=True),\n            nn.Conv2d(mid, channels, 1, bias=False)\n        )\n        self.sigmoid = nn.Sigmoid()\n    def forward(self, x):\n        return x * self.sigmoid(self.mlp(self.avg_pool(x)) + self.mlp(self.max_pool(x)))\n\nclass SpatialAttention(nn.Module):\n    def __init__(self):\n        super().__init__()\n        self.conv = nn.Conv2d(2, 1, kernel_size=7, padding=3, bias=False)\n        self.sigmoid = nn.Sigmoid()\n    def forward(self, x):\n        avg = torch.mean(x, dim=1, keepdim=True)\n        mx, _ = torch.max(x, dim=1, keepdim=True)\n        return x * self.sigmoid(self.conv(torch.cat([avg, mx], dim=1)))\n\nclass CBAM(nn.Module):\n    def __init__(self, channels):\n        super().__init__()\n        self.ca = ChannelAttention(channels)\n        self.sa = SpatialAttention()\n    def forward(self, x):\n        return self.sa(self.ca(x))\n\nclass EncoderBlock(nn.Module):\n    def __init__(self, in_ch, out_ch):\n        super().__init__()\n        self.res  = ResBlock(in_ch, out_ch)\n        self.pool = nn.MaxPool2d(2)\n    def forward(self, x):\n        skip = self.res(x)\n        return self.pool(skip), skip\n\nclass DecoderBlock(nn.Module):\n    def __init__(self, in_ch, skip_ch, out_ch):\n        super().__init__()\n        self.up   = nn.Upsample(scale_factor=2, mode='bilinear', align_corners=True)\n        self.res  = ResBlock(in_ch + skip_ch, out_ch)\n        self.cbam = CBAM(out_ch)\n    def forward(self, x, skip):\n        x = self.up(x)\n        if x.shape != skip.shape:\n            dY = skip.shape[2] - x.shape[2]\n            dX = skip.shape[3] - x.shape[3]\n            x  = F.pad(x, [dX//2, dX-dX//2, dY//2, dY-dY//2])\n        return self.cbam(self.res(torch.cat([skip, x], dim=1)))\n\nclass ResUNet(nn.Module):\n    def __init__(self, n_channels=5, n_classes=1, base=64):\n        super().__init__()\n        f = base\n        self.enc1       = EncoderBlock(n_channels, f)\n        self.enc2       = EncoderBlock(f,   f*2)\n        self.enc3       = EncoderBlock(f*2, f*4)\n        self.enc4       = EncoderBlock(f*4, f*8)\n        self.bottleneck = nn.Sequential(ResBlock(f*8, f*16), CBAM(f*16))\n        self.dec4       = DecoderBlock(f*16, f*8,  f*8)\n        self.dec3       = DecoderBlock(f*8,  f*4,  f*4)\n        self.dec2       = DecoderBlock(f*4,  f*2,  f*2)\n        self.dec1       = DecoderBlock(f*2,  f,    f)\n        self.head       = nn.Sequential(\n            nn.Conv2d(f, f//2, 3, padding=1),\n            nn.BatchNorm2d(f//2), nn.ReLU(inplace=True),\n            nn.Dropout2d(0.05),\n            nn.Conv2d(f//2, n_classes, 1)\n        )\n    def forward(self, x):\n        x, s1 = self.enc1(x)\n        x, s2 = self.enc2(x)\n        x, s3 = self.enc3(x)\n        x, s4 = self.enc4(x)\n        x     = self.bottleneck(x)\n        x     = self.dec4(x, s4)\n        x     = self.dec3(x, s3)\n        x     = self.dec2(x, s2)\n        x     = self.dec1(x, s1)\n        return self.head(x)\n\ndevice = config['device']\n_t = ResUNet(config['n_channels'], config['n_classes'], config['base_filters']).to(device)\nwith torch.no_grad():\n    _d = torch.randn(2, 5, 70, 70).to(device)\n    _o = _t(_d)\n    assert _o.shape == (2, 1, 70, 70)\n    print(f'Input : {_d.shape} -> Output : {_o.shape}  OK')\nn_params = sum(p.numel() for p in _t.parameters())\nprint(f'Parametres : {n_params:,}')\ndel _t, _d, _o","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:39.404391Z","iopub.execute_input":"2026-06-02T16:21:39.404758Z","iopub.status.idle":"2026-06-02T16:21:40.871129Z","shell.execute_reply.started":"2026-06-02T16:21:39.404733Z","shell.execute_reply":"2026-06-02T16:21:40.870482Z"}},"outputs":[],"execution_count":null},{"id":"e5553662-81e3-4601-8fa6-55f31f951bba","cell_type":"code","source":"# Architecture 4.1 — Schema ResUNet\nfig, ax = plt.subplots(figsize=(14, 9))\nax.set_xlim(0, 14); ax.set_ylim(0, 9); ax.axis('off')\nax.set_facecolor('#f8f9fa'); fig.patch.set_facecolor('#f8f9fa')\nax.set_title('Architecture ResUNet pour lInversion Sismique\\nEncodeur residuel -> Bottleneck CBAM -> Decodeur avec attention',\n             fontsize=13, fontweight='bold', pad=15)\nenc_color='#4C72B0'; dec_color='#DD8452'; bot_color='#55A868'; skip_color='#C44E52'\ndef draw_block(ax, x, y, w, h, label, sublabel, color, fs=8):\n    rect = mpatches.FancyBboxPatch((x-w/2, y-h/2), w, h,\n        boxstyle='round,pad=0.05', facecolor=color, edgecolor='white', linewidth=1.5, alpha=0.85)\n    ax.add_patch(rect)\n    ax.text(x, y+0.05, label, ha='center', va='center', fontsize=fs, fontweight='bold', color='white')\n    ax.text(x, y-0.25, sublabel, ha='center', va='center', fontsize=6.5, color='white', alpha=0.9)\nenc_data = [(2.0,7.5,'Input','5x70x70','#2d6a9f'),(2.0,6.2,'EncBlock 1','64x35x35',enc_color),\n            (2.0,5.0,'EncBlock 2','128x17x17',enc_color),(2.0,3.8,'EncBlock 3','256x8x8',enc_color),\n            (2.0,2.6,'EncBlock 4','512x4x4',enc_color)]\ndec_data = [(12.0,2.6,'DecBlock 4','512x8x8',dec_color),(12.0,3.8,'DecBlock 3','256x17x17',dec_color),\n            (12.0,5.0,'DecBlock 2','128x35x35',dec_color),(12.0,6.2,'DecBlock 1','64x70x70',dec_color),\n            (12.0,7.5,'Output','1x70x70','#b35900')]\nfor x,y,l,s,c in enc_data: draw_block(ax,x,y,2.2,0.75,l,s,c)\ndraw_block(ax,7.0,1.6,2.6,0.85,'Bottleneck','1024x4x4 + CBAM',bot_color,fs=9)\nfor x,y,l,s,c in dec_data: draw_block(ax,x,y,2.2,0.75,l,s,c)\nfor i in range(len(enc_data)-1):\n    ax.annotate('',xy=(2.0,enc_data[i+1][1]+0.38),xytext=(2.0,enc_data[i][1]-0.38),\n                arrowprops=dict(arrowstyle='->',color='#333',lw=1.5))\nax.annotate('',xy=(5.7,1.6),xytext=(3.1,2.6),arrowprops=dict(arrowstyle='->',color='#333',lw=1.5))\nax.annotate('',xy=(10.9,2.6),xytext=(8.3,1.6),arrowprops=dict(arrowstyle='->',color='#333',lw=1.5))\nfor i in range(len(dec_data)-1):\n    ax.annotate('',xy=(12.0,dec_data[i+1][1]-0.38),xytext=(12.0,dec_data[i][1]+0.38),\n                arrowprops=dict(arrowstyle='->',color='#333',lw=1.5))\nskip_pairs = [(enc_data[1][1],dec_data[3][1]),(enc_data[2][1],dec_data[2][1]),\n              (enc_data[3][1],dec_data[1][1]),(enc_data[4][1],dec_data[0][1])]\nfor ye,yd in skip_pairs:\n    ax.annotate('',xy=(10.9,yd),xytext=(3.1,ye),\n                arrowprops=dict(arrowstyle='->',color=skip_color,lw=1.5,linestyle='dashed',\n                                connectionstyle='arc3,rad=0.2'))\nlegend_elements = [\n    mpatches.Patch(facecolor=enc_color, label='Encodeur (ResBlock+Pool)'),\n    mpatches.Patch(facecolor=bot_color, label='Bottleneck (ResBlock+CBAM)'),\n    mpatches.Patch(facecolor=dec_color, label='Decodeur (Up+ResBlock+CBAM)'),\n    mpatches.Patch(facecolor=skip_color,label='Skip connections'),\n]\nax.legend(handles=legend_elements, loc='lower center', ncol=2, fontsize=9,\n          framealpha=0.8, bbox_to_anchor=(0.5,-0.02))\nplt.tight_layout()\nplt.savefig('architecture_resunet.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:40.872202Z","iopub.execute_input":"2026-06-02T16:21:40.872957Z","iopub.status.idle":"2026-06-02T16:21:41.540433Z","shell.execute_reply.started":"2026-06-02T16:21:40.872921Z","shell.execute_reply":"2026-06-02T16:21:41.539699Z"}},"outputs":[],"execution_count":null},{"id":"dec9f748-36ac-443f-ab4c-bcba9c386ba6","cell_type":"markdown","source":"**Flux de l'architecture** : l'encodeur compresse progressivement l'information spatiale (70×70 → 4×4) en extrayant des features de plus en plus abstraites. Le bottleneck (1024×4×4 + CBAM) est le point de compression maximale. Le décodeur reconstruit la résolution spatiale en intégrant les skip connections — qui transmettent directement les détails fins de l'encodeur pour éviter leur perte lors de la compression.","metadata":{}},{"id":"22e02f13-f660-466a-aea8-34ff941b7fc6","cell_type":"code","source":"# Architecture 4.2 — Repartition des parametres\nmodel_viz     = ResUNet(config['n_channels'], config['n_classes'], config['base_filters'])\nmodule_params = {name: sum(p.numel() for p in m.parameters()) for name, m in model_viz.named_children()}\nfig, ax = plt.subplots(figsize=(10, 5))\nnames      = list(module_params.keys())\nvals       = [v/1e6 for v in module_params.values()]\ncolors_bar = plt.cm.RdYlGn(np.linspace(0.2, 0.9, len(names)))\nbars = ax.bar(names, vals, color=colors_bar, edgecolor='white', linewidth=0.8)\nax.set_xlabel('Module', fontsize=11)\nax.set_ylabel('Parametres (millions)', fontsize=11)\nax.set_title(f'Repartition des parametres par module\\nTotal : {sum(module_params.values())/1e6:.2f}M',\n             fontsize=12, fontweight='bold')\nax.tick_params(axis='x', rotation=30)\nax.grid(True, alpha=0.3, axis='y')\nfor bar, val in zip(bars, vals):\n    ax.text(bar.get_x()+bar.get_width()/2, bar.get_height()+0.02, f'{val:.2f}M',\n            ha='center', va='bottom', fontsize=9)\nplt.tight_layout()\nplt.savefig('architecture_params.png', dpi=120, bbox_inches='tight')\nplt.show()\ndel model_viz","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:41.543611Z","iopub.execute_input":"2026-06-02T16:21:41.543852Z","iopub.status.idle":"2026-06-02T16:21:42.183417Z","shell.execute_reply.started":"2026-06-02T16:21:41.543830Z","shell.execute_reply":"2026-06-02T16:21:42.182744Z"}},"outputs":[],"execution_count":null},{"id":"59e7411d-59f5-4961-a90e-56672c0b5fd4","cell_type":"markdown","source":"**Répartition des paramètres** : le bottleneck (14.95M) et le premier bloc du décodeur dec4 (10.29M) concentrent ~75% des paramètres totaux. C'est cohérent — ces niveaux opèrent sur les représentations les plus abstraites et les plus informationnellement riches. Les niveaux superficiels (enc1, dec1) ont peu de paramètres car ils traitent des features bas-niveau (détecteurs de bords, gradients).","metadata":{}},{"id":"12aae55d-ebac-40b9-afe9-8ab49886acdb","cell_type":"markdown","source":"## Section 5 — Fonction de Perte Composite\n\nUne loss unique (ex: MSE pure) est insuffisante pour la FWI car elle ne pénalise pas de la même façon :\n- Les erreurs d'amplitude globale\n- Les erreurs sur les discontinuités (interfaces)\n- Les erreurs de structure spatiale\n\nNotre loss composite combine **4 termes complémentaires** :\n\n| Terme | Poids | Rôle |\n|---|---|---|\n| **MSE** | 0.40 | Précision globale des amplitudes de vitesse |\n| **L1** | 0.35 | Robustesse aux valeurs extrêmes (outliers géologiques) |\n| **SSIM** | 0.15 | Cohérence structurelle des couches (perceptual loss) |\n| **Gradient** | 0.10 | Préservation des interfaces nettes entre couches |\n\n> **Note** : la Focal Frequency Loss (FFT) initialement envisagée a été abandonnée car cuFFT en float16 ne supporte pas les dimensions non-puissances-de-2 (70×70). La Gradient Loss offre une alternative stable et géophysiquement motivée.","metadata":{}},{"id":"09bad770-cddf-420f-9633-cceb123fc9d6","cell_type":"code","source":"# ============================================================\n#  SECTION 5 — FONCTION DE PERTE COMPOSITE\n# ============================================================\nclass GradientLoss(nn.Module):\n    def __init__(self):\n        super().__init__()\n    def forward(self, pred, target):\n        pred_dx   = pred[:, :, :, 1:]   - pred[:, :, :, :-1]\n        pred_dy   = pred[:, :, 1:, :]   - pred[:, :, :-1, :]\n        target_dx = target[:, :, :, 1:] - target[:, :, :, :-1]\n        target_dy = target[:, :, 1:, :] - target[:, :, :-1, :]\n        return F.l1_loss(pred_dx, target_dx) + F.l1_loss(pred_dy, target_dy)\n\nclass SeismicLoss(nn.Module):\n    def __init__(self, w_mse=0.40, w_l1=0.35, w_ssim=0.15, w_grad=0.10):\n        super().__init__()\n        self.mse  = nn.MSELoss()\n        self.l1   = nn.L1Loss()\n        self.grad = GradientLoss()\n        self.w_mse=w_mse; self.w_l1=w_l1\n        self.w_ssim=w_ssim; self.w_grad=w_grad\n    def ssim_loss(self, pred, target, w=3):\n        mu1 = F.avg_pool2d(pred,   w, 1, w//2)\n        mu2 = F.avg_pool2d(target, w, 1, w//2)\n        s1  = F.avg_pool2d(pred**2,     w, 1, w//2) - mu1**2\n        s2  = F.avg_pool2d(target**2,   w, 1, w//2) - mu2**2\n        s12 = F.avg_pool2d(pred*target, w, 1, w//2) - mu1*mu2\n        C1, C2 = 0.01**2, 0.03**2\n        ssim = ((2*mu1*mu2+C1)*(2*s12+C2)) / ((mu1**2+mu2**2+C1)*(s1+s2+C2))\n        return 1 - ssim.mean()\n    def forward(self, pred, target):\n        l_mse  = self.mse(pred, target)\n        l_l1   = self.l1(pred, target)\n        l_ssim = self.ssim_loss(pred, target)\n        l_grad = self.grad(pred, target)\n        total  = self.w_mse*l_mse + self.w_l1*l_l1 + self.w_ssim*l_ssim + self.w_grad*l_grad\n        return total, {'mse':l_mse.item(),'l1':l_l1.item(),'ssim':l_ssim.item(),'grad':l_grad.item()}\n\nprint('SeismicLoss (MSE+L1+SSIM+Grad) definie.')\n\n# Section 5.1 — Visualisation LR schedule\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\nfig.suptitle('Strategie doptimisation : AdamW + Warmup 5 epochs + CosineAnnealing', fontsize=13, fontweight='bold')\n_sm = nn.Linear(10, 1)\n_so = torch.optim.AdamW(_sm.parameters(), lr=config['learning_rate'])\nwarmup_ep = 5\ncosine_ep = config['epochs'] - warmup_ep\n_sw = lr_scheduler.LinearLR(_so, start_factor=0.1, end_factor=1.0, total_iters=warmup_ep)\n_sc = lr_scheduler.CosineAnnealingLR(_so, T_max=cosine_ep, eta_min=1e-7)\n_ss = lr_scheduler.SequentialLR(_so, schedulers=[_sw, _sc], milestones=[warmup_ep])\nlrs_sim = []\nfor ep in range(config['epochs']):\n    lrs_sim.append(_so.param_groups[0]['lr'])\n    _ss.step()\nep_ax = np.arange(1, len(lrs_sim)+1)\naxes[0].plot(ep_ax, lrs_sim, color='steelblue', linewidth=2)\naxes[0].fill_between(ep_ax, lrs_sim, alpha=0.2, color='steelblue')\naxes[0].axvline(warmup_ep, color='red', linestyle='--', alpha=0.7, label=f'Fin warmup ep{warmup_ep}')\naxes[0].set_xlabel('Epoque'); axes[0].set_ylabel('Learning Rate')\naxes[0].set_title('Warmup 5 epochs puis CosineAnnealing')\naxes[0].legend(); axes[0].grid(True, alpha=0.3)\naxes[1].bar(['ResUNet 33M\\n(notre choix)', 'ResUNet 8M\\n(trop petit)', 'U-Net original\\n(reference)'],\n            [33, 8, 33], color=['#55A868', '#C44E52', '#4C72B0'], alpha=0.8, edgecolor='white')\naxes[1].set_ylabel('Millions de parametres'); axes[1].grid(True, alpha=0.3, axis='y')\naxes[1].set_title('Comparaison capacite des modeles')\nfor i,(v,l) in enumerate(zip([33,8,33],['optimal','insuffisant','reference'])):\n    axes[1].text(i, v+0.3, l, ha='center', fontsize=10, fontweight='bold')\nplt.tight_layout()\nplt.savefig('lr_schedule.png', dpi=120, bbox_inches='tight')\nplt.show()\ndel _sm, _so, _ss, _sw, _sc","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:42.184215Z","iopub.execute_input":"2026-06-02T16:21:42.184541Z","iopub.status.idle":"2026-06-02T16:21:45.740916Z","shell.execute_reply.started":"2026-06-02T16:21:42.184503Z","shell.execute_reply":"2026-06-02T16:21:45.740233Z"}},"outputs":[],"execution_count":null},{"id":"4cba81fd-5355-4b81-ae2c-c208ce64870f","cell_type":"markdown","source":"**Stratégie d'optimisation** : le warmup linéaire sur 5 epochs permet d'initialiser les poids progressivement, évitant les mises à jour trop agressives en début d'entraînement. Le CosineAnnealing réduit ensuite le learning rate de façon douce et monotone — sans aucun \"choc\" de réduction brutale qui pourrait déstabiliser la convergence.\n\nLe choix du **float32 pur** (sans AMP) est délibéré : nos tests ont montré que la SSIM loss combinée à l'optimisation float16 génère des instabilités numériques sur des tenseurs 70×70. La légère perte de vitesse (~2×) est entièrement compensée par la stabilité de la convergence.","metadata":{}},{"id":"bf0dfc08-50f1-441d-ad50-6a3b20389908","cell_type":"code","source":"# ============================================================\n#  SECTION 6 — DATASET & DATALOADERS\n# ============================================================\nclass SeismicDataset(Dataset):\n    def __init__(self, x, y, augment=False, noise_std_max=0.03):\n        self.x = torch.tensor(x, dtype=torch.float32)\n        self.y = torch.tensor(y, dtype=torch.float32)\n        self.augment      = augment\n        self.noise_std_max = noise_std_max\n    def __len__(self):\n        return len(self.x)\n    def __getitem__(self, idx):\n        x, y = self.x[idx].clone(), self.y[idx].clone()\n        if self.augment:\n            if torch.rand(1).item() > 0.5:\n                std = torch.rand(1).item() * self.noise_std_max\n                x = x + torch.randn_like(x) * std\n            if torch.rand(1).item() > 0.7:\n                x = torch.flip(x, dims=[2])\n                y = torch.flip(y, dims=[2])\n        return x, y\n\ntrain_ds = SeismicDataset(x_train_norm, y_train_norm, augment=True)\nval_ds   = SeismicDataset(x_val_norm,   y_val_norm,   augment=False)\n\ntrain_loader = DataLoader(train_ds, batch_size=config['batch_size'], shuffle=True,\n                          num_workers=config['num_workers'], pin_memory=config['pin_memory'], drop_last=True)\nval_loader   = DataLoader(val_ds,   batch_size=config['batch_size'], shuffle=False,\n                          num_workers=config['num_workers'], pin_memory=config['pin_memory'])\n\nprint(f'Train    : {len(train_ds)}')\nprint(f'Val      : {len(val_ds)}')\nprint(f'Batches  : {len(train_loader)} train / {len(val_loader)} val')\nprint(f'Augmentation train : {train_ds.augment}')\nprint(f'Augmentation val   : {val_ds.augment}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:45.741769Z","iopub.execute_input":"2026-06-02T16:21:45.742260Z","iopub.status.idle":"2026-06-02T16:21:46.008718Z","shell.execute_reply.started":"2026-06-02T16:21:45.742233Z","shell.execute_reply":"2026-06-02T16:21:46.007952Z"}},"outputs":[],"execution_count":null},{"id":"2598c90a-796d-4689-983f-0ad17b10aeb4","cell_type":"markdown","source":"## Section 7 — Entraînement\n\nLe modèle est entraîné sur **4800 échantillons** (80% du dataset) avec une validation sur **1200 échantillons**. Le split est effectué **avant** la normalisation pour éviter toute fuite de données statistiques du train vers la validation.\n\nPoints clés de la boucle d'entraînement :\n- `optimizer.zero_grad(set_to_none=True)` : libération mémoire optimisée\n- `clip_grad_norm_(max_norm=1.0)` : contrôle des gradients explosifs\n- Sauvegarde du meilleur modèle à chaque amélioration de la val loss\n- Early stopping avec patience=15 epochs","metadata":{}},{"id":"cefcb0f1-e7ad-43c0-b629-f650b5696bb1","cell_type":"code","source":"# ============================================================\n#  SECTION 7 — ENTRAINEMENT SANS AMP\n#  Float32 pur — stable avec loss composite SSIM+Gradient\n#  Warmup 5 epochs + CosineAnnealing\n# ============================================================\ndef train_model(model, train_loader, val_loader, config):\n    device    = config['device']\n    criterion = SeismicLoss()\n    optimizer = torch.optim.AdamW(\n        model.parameters(),\n        lr=config['learning_rate'],\n        weight_decay=config['weight_decay'],\n        betas=(0.9, 0.999),\n        eps=1e-8\n    )\n\n    warmup_epochs = 5\n    cosine_epochs = config['epochs'] - warmup_epochs\n\n    warmup_scheduler = lr_scheduler.LinearLR(\n        optimizer, start_factor=0.1, end_factor=1.0, total_iters=warmup_epochs\n    )\n    cosine_scheduler = lr_scheduler.CosineAnnealingLR(\n        optimizer, T_max=cosine_epochs, eta_min=1e-7\n    )\n    scheduler = lr_scheduler.SequentialLR(\n        optimizer,\n        schedulers=[warmup_scheduler, cosine_scheduler],\n        milestones=[warmup_epochs]\n    )\n\n    best_val, patience_counter = float('inf'), 0\n    train_losses, val_losses, lr_history, epoch_times = [], [], [], []\n    loss_components = {'mse':[], 'l1':[], 'ssim':[], 'grad':[]}\n\n    print(f'Debut entrainement sur {device}')\n    print(f'AMP desactive — float32 pur pour stabilite')\n    print(f'Warmup {warmup_epochs} ep → CosineAnnealing {cosine_epochs} ep')\n    print(f'LR max={config[\"learning_rate\"]} | clip={config[\"gradient_clip\"]}')\n    print('='*65)\n\n    for epoch in range(config['epochs']):\n        t_start = time.time()\n        model.train()\n        train_loss, epoch_comps = 0.0, {k:0.0 for k in loss_components}\n\n        for xb, yb in train_loader:\n            xb, yb = xb.to(device), yb.to(device)\n            optimizer.zero_grad(set_to_none=True)\n            pred = model(xb)\n            loss, comps = criterion(pred, yb)\n            loss.backward()\n            torch.nn.utils.clip_grad_norm_(model.parameters(), config['gradient_clip'])\n            optimizer.step()\n            train_loss += loss.item()\n            for k in comps: epoch_comps[k] += comps[k]\n\n        model.eval()\n        val_loss = 0.0\n        with torch.no_grad():\n            for xb, yb in val_loader:\n                xb, yb = xb.to(device), yb.to(device)\n                pred = model(xb)\n                loss, _ = criterion(pred, yb)\n                val_loss += loss.item()\n\n        avg_train  = train_loss / len(train_loader)\n        avg_val    = val_loss   / len(val_loader)\n        elapsed    = time.time() - t_start\n        current_lr = optimizer.param_groups[0]['lr']\n\n        train_losses.append(avg_train); val_losses.append(avg_val)\n        lr_history.append(current_lr);  epoch_times.append(elapsed)\n        for k in epoch_comps:\n            loss_components[k].append(epoch_comps[k] / len(train_loader))\n\n        scheduler.step()\n\n        print(f'Ep {epoch+1:3d}/{config[\"epochs\"]} | '\n              f'Train:{avg_train:.5f} Val:{avg_val:.5f} '\n              f'LR:{current_lr:.2e} {elapsed:.0f}s')\n\n        if avg_val < best_val:\n            best_val = avg_val; patience_counter = 0\n            torch.save({\n                'model_state_dict':     model.state_dict(),\n                'optimizer_state_dict': optimizer.state_dict(),\n                'epoch':  epoch,\n                'loss':   best_val,\n                'y_mean': y_mean,\n                'y_std':  y_std,\n                'config': config\n            }, config['model_path'])\n            print(f'  Meilleur modele sauvegarde (val={best_val:.5f})')\n        else:\n            patience_counter += 1\n            if patience_counter >= config['patience']:\n                print(f'Early stopping ep {epoch+1}')\n                break\n\n    print('='*65)\n    print(f'Termine. Meilleure val loss : {best_val:.5f}')\n    return train_losses, val_losses, loss_components, lr_history, epoch_times\n\n\nmodel = ResUNet(\n    config['n_channels'],\n    config['n_classes'],\n    config['base_filters']\n).to(device)\n\nn_params = sum(p.numel() for p in model.parameters())\nprint(f'ResUNet base={config[\"base_filters\"]} : {n_params:,} parametres')\n\ntrain_losses, val_losses, loss_components, lr_history, epoch_times = \\\n    train_model(model, train_loader, val_loader, config)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T16:21:46.009767Z","iopub.execute_input":"2026-06-02T16:21:46.010146Z","iopub.status.idle":"2026-06-02T18:01:49.065628Z","shell.execute_reply.started":"2026-06-02T16:21:46.010104Z","shell.execute_reply":"2026-06-02T18:01:49.064866Z"}},"outputs":[],"execution_count":null},{"id":"f14a485a-3e19-4b08-9434-0a0d818b6017","cell_type":"code","source":"# Section 7.1 — Courbes d entrainement\ndef moving_average(x, w=3):\n    return np.convolve(x, np.ones(w)/w, mode='valid')\nn_ep    = len(train_losses)\nep_axis = np.arange(1, n_ep+1)\nfig, axes = plt.subplots(2, 3, figsize=(18, 10))\nfig.suptitle('Monitoring complet de lentrainement ResUNet', fontsize=14, fontweight='bold')\naxes[0,0].plot(ep_axis, train_losses, color='steelblue', alpha=0.4, lw=1, label='Train brut')\naxes[0,0].plot(ep_axis, val_losses,   color='salmon',   alpha=0.4, lw=1, label='Val brut')\nif n_ep > 3:\n    axes[0,0].plot(ep_axis[2:], moving_average(train_losses,3), color='steelblue', lw=2.5, label='Train lisse')\n    axes[0,0].plot(ep_axis[2:], moving_average(val_losses,3),   color='red',       lw=2.5, label='Val lisse')\nbest_ep = np.argmin(val_losses)+1\naxes[0,0].axvline(best_ep, color='green', linestyle='--', alpha=0.8, label=f'Best ep{best_ep}')\naxes[0,0].set_xlabel('Epoque'); axes[0,0].set_ylabel('Loss')\naxes[0,0].set_title('Loss globale Train vs Val'); axes[0,0].legend(fontsize=8); axes[0,0].grid(True,alpha=0.3)\naxes[0,1].plot(ep_axis, lr_history, color='purple', lw=2)\naxes[0,1].fill_between(ep_axis, lr_history, alpha=0.2, color='purple')\naxes[0,1].set_xlabel('Epoque'); axes[0,1].set_ylabel('Learning Rate')\naxes[0,1].set_title('Evolution du Learning Rate'); axes[0,1].grid(True,alpha=0.3)\naxes[0,2].bar(ep_axis, epoch_times, color='teal', alpha=0.7, edgecolor='white')\naxes[0,2].axhline(np.mean(epoch_times), color='red', linestyle='--', lw=2, label=f'Moy:{np.mean(epoch_times):.0f}s')\naxes[0,2].set_xlabel('Epoque'); axes[0,2].set_ylabel('Temps (s)')\naxes[0,2].set_title(f'Duree par epoque | Total:{sum(epoch_times)/60:.1f} min')\naxes[0,2].legend(); axes[0,2].grid(True,alpha=0.3,axis='y')\ncomp_colors = {'mse':'#E74C3C','l1':'#3498DB','ssim':'#2ECC71','grad':'#F39C12'}\nfor i,(key,vals) in enumerate(loss_components.items()):\n    row,col = 1, i%3\n    c = comp_colors[key]\n    axes[row,col].plot(ep_axis, vals, color=c, lw=2, label=key.upper())\n    if n_ep>3: axes[row,col].plot(ep_axis[2:], moving_average(vals,3), color=c, lw=3, alpha=0.6, ls='--')\n    axes[row,col].set_xlabel('Epoque'); axes[row,col].set_ylabel('Loss')\n    axes[row,col].set_title(f'Composante {key.upper()}'); axes[row,col].legend(fontsize=9)\n    axes[row,col].grid(True,alpha=0.3)\nplt.tight_layout()\nplt.savefig('training_monitoring.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:01:49.067159Z","iopub.execute_input":"2026-06-02T18:01:49.067589Z","iopub.status.idle":"2026-06-02T18:01:51.196639Z","shell.execute_reply.started":"2026-06-02T18:01:49.067559Z","shell.execute_reply":"2026-06-02T18:01:51.195954Z"}},"outputs":[],"execution_count":null},{"id":"56d5ab52-bfef-4d18-a2d2-089d270461f3","cell_type":"markdown","source":"**Analyse des courbes** :\n\n- **Loss globale** : la convergence est remarquablement stable sur 99 epochs sans aucun spike. Le gap train/val se creuse légèrement après l'epoch 60 — signe d'un début d'overfitting contrôlé, le meilleur modèle (epoch 84) se situe précisément avant ce point.\n\n- **Learning Rate** : le warmup 5 epochs est clairement visible (montée linéaire), suivi d'une descente cosinus propre jusqu'à ~1e-7. Cette stratégie explique la convergence stable.\n\n- **Durée** : 60s/epoch × 99 epochs = **99.5 minutes** total. Raisonnable pour un modèle de 33M paramètres en float32 sur GPU T4.\n\n- **Composantes MSE et L1** : descente régulière, montrant que le modèle améliore constamment sa précision d'amplitude.\n\n- **Composante SSIM** : descend de 0.87 à 0.27 — le modèle apprend progressivement à reproduire la structure spatiale des couches géologiques.\n\n- **Composante Gradient** : stabilisation rapide, le modèle préserve bien les interfaces dès les premières epochs.","metadata":{}},{"id":"d23408a3-9033-42e0-b0d2-e5e7866c556c","cell_type":"code","source":"# ============================================================\n#  SECTION 8 — CHARGEMENT DU MEILLEUR MODELE\n# ============================================================\ncheckpoint   = torch.load(config['model_path'], weights_only=False)\nmodel.load_state_dict(checkpoint['model_state_dict'])\nsaved_y_mean = checkpoint['y_mean']\nsaved_y_std  = checkpoint['y_std']\nsaved_epoch  = checkpoint['epoch']\nprint(f'Modele recharge depuis epoque {saved_epoch+1}')\nprint(f'Denorm params : mean={saved_y_mean:.2f}  std={saved_y_std:.2f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:01:51.197560Z","iopub.execute_input":"2026-06-02T18:01:51.198265Z","iopub.status.idle":"2026-06-02T18:01:51.559628Z","shell.execute_reply.started":"2026-06-02T18:01:51.198238Z","shell.execute_reply":"2026-06-02T18:01:51.558689Z"}},"outputs":[],"execution_count":null},{"id":"3434c7bb-7f88-4810-b7e3-b802d9a2a6b6","cell_type":"markdown","source":"## Section 8 — Chargement du Meilleur Modèle\n\nLe modèle sauvegardé correspond à l'**epoch 84** — point optimal entre convergence maximale et début d'overfitting. Les paramètres de dénormalisation (mean=2849.47 m/s, std=776.54 m/s) permettent de reconvertir les prédictions normalisées vers les vitesses réelles en m/s.","metadata":{}},{"id":"f335e85c-123c-4f4b-a18e-491870c29d0b","cell_type":"code","source":"# ============================================================\n#  SECTION 9 — EVALUATION ENRICHIE\n# ============================================================\nmodel.eval()\nall_preds, all_targets = [], []\nwith torch.no_grad():\n    for xb, yb in val_loader:\n        xb = xb.to(device)\n        pred = model(xb)\n        all_preds.append(pred.cpu().squeeze(1).numpy())\n        all_targets.append(yb.squeeze(1).numpy())\nall_preds      = np.concatenate(all_preds,   axis=0)\nall_targets    = np.concatenate(all_targets, axis=0)\npreds_denorm   = all_preds   * saved_y_std + saved_y_mean\ntargets_denorm = all_targets * saved_y_std + saved_y_mean\nmse_all  = mean_squared_error(targets_denorm.flatten(), preds_denorm.flatten())\nmae_all  = mean_absolute_error(targets_denorm.flatten(), preds_denorm.flatten())\nrmse_all = np.sqrt(mse_all)\ncorr_all = np.corrcoef(targets_denorm.flatten(), preds_denorm.flatten())[0,1]\nprint('='*60)\nprint(f'METRIQUES VALIDATION ({len(all_preds)} echantillons)')\nprint('='*60)\nprint(f'  RMSE        : {rmse_all:.2f} m/s')\nprint(f'  MAE         : {mae_all:.2f} m/s')\nprint(f'  Correlation : {corr_all:.4f}')\nprint('='*60)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:01:51.561358Z","iopub.execute_input":"2026-06-02T18:01:51.561874Z","iopub.status.idle":"2026-06-02T18:01:56.728385Z","shell.execute_reply.started":"2026-06-02T18:01:51.561847Z","shell.execute_reply":"2026-06-02T18:01:56.727509Z"}},"outputs":[],"execution_count":null},{"id":"6692a86a-93fd-4490-a9ce-630e9dc26534","cell_type":"markdown","source":"## Section 9 — Évaluation Enrichie\n\n### Métriques globales sur 1200 échantillons de validation\n\n| Métrique | Valeur | Interprétation |\n|---|---|---|\n| **RMSE** | 292 m/s | Erreur quadratique moyenne — ≈10% de la plage 1500-4500 m/s |\n| **MAE** | 160 m/s | Erreur absolue moyenne — erreur typique par pixel |\n| **Corrélation** | 0.9262 | Très forte corrélation linéaire entre prédit et réel |\n\nUne corrélation de **0.926** est un résultat solide pour de la FWI par deep learning. Le MAE de 160 m/s représente environ **5.6%** de la plage de vitesses totale.","metadata":{}},{"id":"7cacee5f-52b8-4316-afb8-6cc0f2107eb9","cell_type":"code","source":"# Evaluation 9.1 — Grille visuelle\ntest_indices = [i for i in [0,1,2,5,10,15] if i < len(all_preds)]\nfig, axes = plt.subplots(4, len(test_indices), figsize=(22, 14))\nfig.suptitle('Evaluation visuelle — Reel vs Predit vs Erreurs\\nL1:Reel | L2:Predit | L3:Erreur abs | L4:Erreur rel(%)',\n             fontsize=12, fontweight='bold')\nrow_labels = ['Reel (m/s)','Predit (m/s)','Erreur |delta| (m/s)','Erreur rel (%)']\nfor i,idx in enumerate(test_indices):\n    t=targets_denorm[idx]; p=preds_denorm[idx]\n    err_abs=np.abs(t-p); err_rel=err_abs/(np.abs(t)+1e-6)*100\n    mse_i=mean_squared_error(t.flatten(),p.flatten())\n    mae_i=mean_absolute_error(t.flatten(),p.flatten())\n    corr_i=np.corrcoef(t.flatten(),p.flatten())[0,1]\n    vmin,vmax=min(t.min(),p.min()),max(t.max(),p.max())\n    im0=axes[0,i].imshow(t, cmap='jet', vmin=vmin, vmax=vmax, aspect='auto')\n    axes[0,i].set_title(f'#{idx+1} MSE={mse_i:.0f}', fontsize=8)\n    plt.colorbar(im0, ax=axes[0,i], fraction=0.046, pad=0.04)\n    im1=axes[1,i].imshow(p, cmap='jet', vmin=vmin, vmax=vmax, aspect='auto')\n    axes[1,i].set_title(f'Corr={corr_i:.3f}', fontsize=8)\n    plt.colorbar(im1, ax=axes[1,i], fraction=0.046, pad=0.04)\n    im2=axes[2,i].imshow(err_abs, cmap='Reds', aspect='auto')\n    axes[2,i].set_title(f'MAE={mae_i:.0f}', fontsize=8)\n    plt.colorbar(im2, ax=axes[2,i], fraction=0.046, pad=0.04)\n    im3=axes[3,i].imshow(err_rel, cmap='plasma', vmax=30, aspect='auto')\n    axes[3,i].set_title(f'Max={err_rel.max():.1f}%', fontsize=8)\n    plt.colorbar(im3, ax=axes[3,i], fraction=0.046, pad=0.04)\n    if i==0:\n        for row in range(4): axes[row,i].set_ylabel(row_labels[row], fontsize=8)\nplt.tight_layout()\nplt.savefig('evaluation_grille.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:01:56.729787Z","iopub.execute_input":"2026-06-02T18:01:56.730090Z","iopub.status.idle":"2026-06-02T18:02:04.259829Z","shell.execute_reply.started":"2026-06-02T18:01:56.730062Z","shell.execute_reply":"2026-06-02T18:02:04.259050Z"}},"outputs":[],"execution_count":null},{"id":"53c8c242-4235-46e2-a9f3-b4bb60adf3d0","cell_type":"markdown","source":"**Évaluation visuelle** : les 6 échantillons montrent une qualité variable selon la complexité géologique :\n- **Échantillons 3, 4, 6, 11** (structures en couches horizontales) : prédictions quasi-parfaites, corrélations > 0.93, erreurs relatives < 15%\n- **Échantillons 1, 2** (structures ondulées complexes) : le modèle capture la tendance générale mais lisse les détails fins des ondulations latérales\n\nCette différence de performance selon la complexité géologique est attendue et documentée dans la littérature FWI — les structures non-horizontales restent un défi ouvert.","metadata":{}},{"id":"bebfc455-a4e8-43c7-8480-3a4cd165d9f7","cell_type":"code","source":"# Evaluation 9.2 — Scatter Predit vs Reel\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\nfig.suptitle('Scatter Plot : Vitesse Predite vs Vitesse Reelle\\nDiagonale parfaite = prediction parfaite', fontsize=13, fontweight='bold')\nfor i,idx in enumerate([0,5,10]):\n    if idx >= len(all_preds): continue\n    t_f=targets_denorm[idx].flatten(); p_f=preds_denorm[idx].flatten()\n    axes[i].scatter(t_f, p_f, s=2, alpha=0.3, color='steelblue')\n    lims=[min(t_f.min(),p_f.min()),max(t_f.max(),p_f.max())]\n    axes[i].plot(lims, lims, 'r--', lw=2, label='y=x (parfait)')\n    axes[i].set_xlabel('Vitesse reelle (m/s)', fontsize=10)\n    axes[i].set_ylabel('Vitesse predite (m/s)', fontsize=10)\n    corr=np.corrcoef(t_f,p_f)[0,1]; mae=mean_absolute_error(t_f,p_f)\n    axes[i].set_title(f'Echantillon #{idx+1}\\nCorr={corr:.4f} MAE={mae:.0f} m/s', fontsize=10)\n    axes[i].legend(fontsize=9); axes[i].grid(True,alpha=0.3)\nplt.tight_layout()\nplt.savefig('scatter_pred_vs_real.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:02:04.260796Z","iopub.execute_input":"2026-06-02T18:02:04.261142Z","iopub.status.idle":"2026-06-02T18:02:05.471379Z","shell.execute_reply.started":"2026-06-02T18:02:04.261116Z","shell.execute_reply":"2026-06-02T18:02:05.470526Z"}},"outputs":[],"execution_count":null},{"id":"7f01e0cb-44fd-4c2e-8e0c-e413844937cd","cell_type":"markdown","source":"**Scatter plot** : la diagonale rouge (prédiction parfaite) sert de référence. L'échantillon 1 présente une dispersion notable (Corr=0.827) due à sa structure géologique complexe. Les échantillons 6 et 11 montrent des corrélations exceptionnelles (0.987 et 0.973) — les points s'alignent quasi-parfaitement sur la diagonale. Les colonnes verticales visibles correspondent aux interfaces géologiques discrètes (vitesses constantes par couche).","metadata":{}},{"id":"4593471a-347c-4888-a71f-fd26c06c7f22","cell_type":"code","source":"# Evaluation 9.3 — Distribution des erreurs\nerrors_flat = (preds_denorm - targets_denorm).flatten()\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\nfig.suptitle('Analyse de la distribution des erreurs', fontsize=13, fontweight='bold')\naxes[0].hist(errors_flat, bins=80, color='steelblue', edgecolor='white', linewidth=0.3, alpha=0.85)\naxes[0].axvline(0, color='red', lw=2, label='Erreur zero')\naxes[0].axvline(errors_flat.mean(), color='orange', lw=2, linestyle='--', label=f'Biais:{errors_flat.mean():.1f} m/s')\naxes[0].set_xlabel('Erreur (m/s)'); axes[0].set_ylabel('Frequence')\naxes[0].set_title('Distribution erreurs signees'); axes[0].legend(fontsize=9); axes[0].grid(True,alpha=0.3)\ndeciles=np.percentile(targets_denorm.flatten(), np.linspace(0,100,11))\nmae_per_decile=[]\nfor j in range(10):\n    mask=(targets_denorm.flatten()>=deciles[j])&(targets_denorm.flatten()<deciles[j+1])\n    mae_per_decile.append(mean_absolute_error(targets_denorm.flatten()[mask],preds_denorm.flatten()[mask]) if mask.sum()>0 else 0)\ndecile_centers=[(deciles[j]+deciles[j+1])/2 for j in range(10)]\naxes[1].bar(range(10), mae_per_decile, color='salmon', edgecolor='white', alpha=0.85)\naxes[1].set_xticks(range(10))\naxes[1].set_xticklabels([f'{v:.0f}' for v in decile_centers], rotation=45, fontsize=8)\naxes[1].set_xlabel('Vitesse reelle (m/s) - deciles'); axes[1].set_ylabel('MAE (m/s)')\naxes[1].set_title('MAE par plage de vitesse'); axes[1].grid(True,alpha=0.3,axis='y')\n(osm,osr),(slope,intercept,r)=stats.probplot(errors_flat, dist='norm')\naxes[2].scatter(osm[::100], osr[::100], s=3, alpha=0.4, color='steelblue')\nx_line=np.array([osm.min(),osm.max()])\naxes[2].plot(x_line, slope*x_line+intercept, 'r-', lw=2)\naxes[2].set_xlabel('Quantiles theoriques'); axes[2].set_ylabel('Quantiles observes')\naxes[2].set_title(f'Q-Q Plot R2={r**2:.4f}'); axes[2].grid(True,alpha=0.3)\nplt.tight_layout()\nplt.savefig('distribution_erreurs.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:02:05.472479Z","iopub.execute_input":"2026-06-02T18:02:05.472827Z","iopub.status.idle":"2026-06-02T18:02:07.883711Z","shell.execute_reply.started":"2026-06-02T18:02:05.472801Z","shell.execute_reply":"2026-06-02T18:02:07.882961Z"}},"outputs":[],"execution_count":null},{"id":"6aa6e4b0-151d-4a4a-9f39-6bdadca28f45","cell_type":"markdown","source":"**Distribution des erreurs** :\n- **Histogramme** : distribution bien centrée sur zéro avec un biais de seulement -18 m/s — le modèle ne sur-estime ni ne sous-estime systématiquement\n- **MAE par décile** : l'erreur augmente avec la vitesse (100 m/s pour les basses vitesses, ~280 m/s pour les hautes vitesses >4000 m/s) — le modèle a plus de difficulté sur les roches rapides en profondeur\n- **Q-Q plot** : R²=0.774, indiquant des queues légèrement plus lourdes qu'une gaussienne parfaite — cohérent avec la présence de cas géologiquement atypiques","metadata":{}},{"id":"304669d4-b4a7-4b50-92db-71b8ab706453","cell_type":"code","source":"# Evaluation 9.4 — Profils de vitesse verticaux\nfig, axes = plt.subplots(2, 3, figsize=(16, 10))\nfig.suptitle('Profils de vitesse verticaux — Predit vs Reel\\n(validation style geophysicien)', fontsize=13, fontweight='bold')\nfor i,(ax,idx) in enumerate(zip(axes.flatten(),[0,1,2,5,8,10])):\n    if idx >= len(all_preds): ax.axis('off'); continue\n    t=targets_denorm[idx]; p=preds_denorm[idx]\n    depth=np.arange(t.shape[0])\n    for x_pos,color,ls in [(17,'blue','-'),(35,'red','--'),(52,'green',':')]:\n        ax.plot(t[:,x_pos], depth, color=color, lw=2, linestyle=ls, label=f'Reel x={x_pos}')\n        ax.plot(p[:,x_pos], depth, color=color, lw=1.5, linestyle=ls, alpha=0.5, marker='o', markersize=2, label=f'Predit x={x_pos}')\n    ax.invert_yaxis()\n    ax.set_xlabel('Vitesse (m/s)', fontsize=9); ax.set_ylabel('Profondeur', fontsize=9)\n    ax.set_title(f'Echantillon #{idx+1}', fontsize=10); ax.grid(True,alpha=0.3)\n    if i==0: ax.legend(fontsize=7, loc='lower right')\nplt.tight_layout()\nplt.savefig('profils_vitesse.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:02:07.884677Z","iopub.execute_input":"2026-06-02T18:02:07.884960Z","iopub.status.idle":"2026-06-02T18:02:09.716500Z","shell.execute_reply.started":"2026-06-02T18:02:07.884929Z","shell.execute_reply":"2026-06-02T18:02:09.715555Z"}},"outputs":[],"execution_count":null},{"id":"e36d89fa-be8c-4bdd-b10a-4068e3c2b8a0","cell_type":"markdown","source":"**Profils de vitesse** — validation style géophysicien :\n\nC'est la visualisation la plus pertinente pour un géophysicien. Les profils 1D à 3 positions latérales (x=17, 35, 52) montrent :\n- **Échantillons 6, 9, 11** : les échelons de vitesse (interfaces) sont parfaitement localisés en profondeur. Les pointillés (prédits) suivent exactement les courbes pleines (réels).\n- **Échantillon 1** : les oscillations latérales complexes sont partiellement capturées mais lissées — limite connue des architectures convolutionnelles sur les géologies fortement hétérogènes.\n- **Échantillon 2** : les interfaces sont bien détectées mais légèrement décalées en profondeur sur certaines positions.","metadata":{}},{"id":"ac6106e0-eb36-4ed8-be08-74ab16abca27","cell_type":"code","source":"# Evaluation 9.5 — Meilleures et Pires predictions\nsample_maes   = [mean_absolute_error(targets_denorm[i].flatten(), preds_denorm[i].flatten()) for i in range(len(all_preds))]\nworst_indices = np.argsort(sample_maes)[-2:][::-1]\nbest_indices  = np.argsort(sample_maes)[:2]\nall_show      = list(best_indices) + list(worst_indices)\nfig, axes = plt.subplots(3, 4, figsize=(20, 12))\nfig.suptitle('Analyse des extremes : Meilleures et Pires predictions', fontsize=13, fontweight='bold')\nfor col,idx in enumerate(all_show):\n    t=targets_denorm[idx]; p=preds_denorm[idx]\n    err=np.abs(t-p); mae_i=sample_maes[idx]\n    quality='MEILLEURE' if idx in best_indices else 'PIRE'\n    color='green' if idx in best_indices else 'red'\n    vmin,vmax=min(t.min(),p.min()),max(t.max(),p.max())\n    im0=axes[0,col].imshow(t, cmap='jet', vmin=vmin, vmax=vmax, aspect='auto')\n    axes[0,col].set_title(f'{quality} #{col+1}\\nMAE={mae_i:.0f}', fontsize=9, color=color, fontweight='bold')\n    plt.colorbar(im0, ax=axes[0,col], fraction=0.046, pad=0.04)\n    im1=axes[1,col].imshow(p, cmap='jet', vmin=vmin, vmax=vmax, aspect='auto')\n    axes[1,col].set_title('Predit', fontsize=9)\n    plt.colorbar(im1, ax=axes[1,col], fraction=0.046, pad=0.04)\n    diff=p-t; norm=TwoSlopeNorm(vmin=-err.max(), vcenter=0, vmax=err.max())\n    im2=axes[2,col].imshow(diff, cmap='RdBu_r', norm=norm, aspect='auto')\n    axes[2,col].set_title('Erreur signee', fontsize=9)\n    plt.colorbar(im2, ax=axes[2,col], fraction=0.046, pad=0.04)\nplt.tight_layout()\nplt.savefig('analyse_extremes.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:02:09.717894Z","iopub.execute_input":"2026-06-02T18:02:09.718436Z","iopub.status.idle":"2026-06-02T18:02:14.843922Z","shell.execute_reply.started":"2026-06-02T18:02:09.718408Z","shell.execute_reply":"2026-06-02T18:02:14.843026Z"}},"outputs":[],"execution_count":null},{"id":"0427f2ea-b620-4183-ac13-d46eddabbecc","cell_type":"markdown","source":"**Analyse des cas extrêmes** :\n\n- **Meilleures prédictions** (MAE=17 et 18 m/s) : structures à couches horizontales simples. L'erreur signée est quasi-nulle sur toute la carte — performance remarquable.\n- **Pires prédictions** (MAE=734 et 692 m/s) : structures géologiques à ondulations latérales profondes et multiples discontinuités. Le modèle capture la tendance verticale mais rate les détails latéraux fins.\n\nCette analyse des cas extrêmes est utile pour identifier les types de géologies où davantage de données d'entraînement ou une architecture spécialisée serait nécessaire.","metadata":{}},{"id":"3da01fbe-4ee9-4b62-b6c4-6da68c49fd2d","cell_type":"code","source":"# Evaluation 9.6 — Monte Carlo Dropout\ndef mc_dropout_predict(model, x_tensor, n_samples=20):\n    model.eval()\n    for m in model.modules():\n        if isinstance(m, nn.Dropout2d): m.train()\n    preds=[]\n    with torch.no_grad():\n        for _ in range(n_samples):\n            out = model(x_tensor)\n            preds.append(out.cpu().squeeze().numpy())\n    preds=np.stack(preds,axis=0)\n    return preds.mean(axis=0), preds.std(axis=0)\n\nfig, axes = plt.subplots(3, 4, figsize=(20, 12))\nfig.suptitle('Incertitude par Monte Carlo Dropout (20 passes)\\nEcart-type pixel = confiance du modele', fontsize=13, fontweight='bold')\nfor col,idx in enumerate(range(min(4,len(all_preds)))):\n    x_t,y_t = val_ds[idx]\n    x_tensor = x_t.unsqueeze(0).to(device)\n    mean_pred,std_pred = mc_dropout_predict(model, x_tensor)\n    mean_d=mean_pred*saved_y_std+saved_y_mean; std_d=std_pred*saved_y_std\n    tgt_d=y_t.squeeze().numpy()*saved_y_std+saved_y_mean\n    vmin,vmax=min(tgt_d.min(),mean_d.min()),max(tgt_d.max(),mean_d.max())\n    im0=axes[0,col].imshow(tgt_d, cmap='jet', vmin=vmin, vmax=vmax, aspect='auto')\n    axes[0,col].set_title(f'Reel #{idx+1}', fontsize=9)\n    plt.colorbar(im0, ax=axes[0,col], fraction=0.046, pad=0.04)\n    im1=axes[1,col].imshow(mean_d, cmap='jet', vmin=vmin, vmax=vmax, aspect='auto')\n    axes[1,col].set_title('Prediction (moyenne MC)', fontsize=9)\n    plt.colorbar(im1, ax=axes[1,col], fraction=0.046, pad=0.04)\n    im2=axes[2,col].imshow(std_d, cmap='hot', aspect='auto')\n    axes[2,col].set_title(f'Incertitude sigma\\nmax={std_d.max():.0f} m/s', fontsize=9)\n    plt.colorbar(im2, ax=axes[2,col], fraction=0.046, pad=0.04)\nplt.tight_layout()\nplt.savefig('mc_dropout_uncertainty.png', dpi=120, bbox_inches='tight')\nplt.show()\nmodel.eval()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:02:14.845012Z","iopub.execute_input":"2026-06-02T18:02:14.845342Z","iopub.status.idle":"2026-06-02T18:02:19.641942Z","shell.execute_reply.started":"2026-06-02T18:02:14.845319Z","shell.execute_reply":"2026-06-02T18:02:19.641344Z"}},"outputs":[],"execution_count":null},{"id":"72ef5c23-ce6a-4622-9fe0-2a2dedf5bc58","cell_type":"markdown","source":"**Incertitude par Monte Carlo Dropout** :\n\nEn activant le dropout en mode inférence et en faisant 20 passes forward, on obtient une **carte d'incertitude pixel par pixel**. L'écart-type σ représente la confiance du modèle :\n\n- σ max de **69-139 m/s** — bien inférieur au MAE global de 160 m/s, cohérent\n- L'incertitude est **concentrée aux interfaces** entre couches géologiques — exactement là où la prédiction est la plus difficile physiquement\n- Les zones à vitesse uniforme (intérieur des couches) ont une incertitude très faible (zones noires sur les cartes hot)\n\nCette carte d'incertitude pourrait être utilisée opérationnellement pour guider l'acquisition sismique supplémentaire vers les zones incertaines.","metadata":{}},{"id":"0c7b1665-a88c-4826-81d6-5c49e9a812b5","cell_type":"code","source":"# Evaluation 9.7 — Cartes d attention CBAM\ndef get_attention_maps(model, x_tensor):\n    attention_maps=[]\n    def hook_fn(module, input, output): attention_maps.append(output.detach().cpu())\n    hooks=[module.register_forward_hook(hook_fn) for name,module in model.named_modules() if isinstance(module, SpatialAttention)]\n    model.eval()\n    with torch.no_grad(): _ = model(x_tensor)\n    for h in hooks: h.remove()\n    return attention_maps\n\nx_sample,_ = val_ds[0]\nx_tensor   = x_sample.unsqueeze(0).to(device)\natt_maps   = get_attention_maps(model, x_tensor)\nif att_maps:\n    n_maps = min(len(att_maps),4)\n    fig, axes = plt.subplots(1, n_maps+1, figsize=(5*(n_maps+1),4))\n    fig.suptitle('Cartes dattention spatiale CBAM\\nBlanc=forte attention | Noir=ignore', fontsize=12, fontweight='bold')\n    axes[0].imshow(x_sample[0].numpy(), cmap='RdBu_r', aspect='auto')\n    axes[0].set_title('Sismogramme input\\n(canal 1)', fontsize=10)\n    for i in range(n_maps):\n        att=att_maps[i][0,0].numpy()\n        att_r=F.interpolate(torch.tensor(att).unsqueeze(0).unsqueeze(0),\n                            size=(70,70),mode='bilinear',align_corners=False)[0,0].numpy()\n        im=axes[i+1].imshow(att_r, cmap='hot', aspect='auto', vmin=0, vmax=1)\n        axes[i+1].set_title(f'Attention niveau {i+1}\\n(decodeur D{4-i})', fontsize=10)\n        plt.colorbar(im, ax=axes[i+1], fraction=0.046, pad=0.04)\n    plt.tight_layout()\n    plt.savefig('attention_maps.png', dpi=120, bbox_inches='tight')\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:02:19.642832Z","iopub.execute_input":"2026-06-02T18:02:19.643252Z","iopub.status.idle":"2026-06-02T18:02:21.345247Z","shell.execute_reply.started":"2026-06-02T18:02:19.643225Z","shell.execute_reply":"2026-06-02T18:02:21.344275Z"}},"outputs":[],"execution_count":null},{"id":"c36353d2-9087-4755-9f8e-1e0db3236b63","cell_type":"markdown","source":"**Cartes d'attention CBAM** :\n\nLes 4 niveaux d'attention du décodeur révèlent comment le modèle \"regarde\" l'image à différentes résolutions :\n- **Niveau 1 (D4, basse résolution)** : attention concentrée sur la zone centrale — détection des structures profondes\n- **Niveau 2 (D3)** : attention sur les zones d'interface principales\n- **Niveaux 3-4 (D2, D1, haute résolution)** : attention diffuse sur les discontinuités latérales fines\n\nCe comportement est géophysiquement cohérent : le modèle focalise son attention précisément sur les zones où l'information est la plus utile pour la reconstruction.","metadata":{}},{"id":"5cb780d6-2344-4283-ac85-0305e97d06fe","cell_type":"code","source":"# Evaluation 9.8 — Feature maps encodeur\ndef get_encoder_features(model, x_tensor):\n    features=[]\n    def hook_fn(module,input,output): features.append(output.detach().cpu())\n    hook=model.enc1.res.register_forward_hook(hook_fn)\n    model.eval()\n    with torch.no_grad(): _ = model(x_tensor)\n    hook.remove()\n    return features[0]\n\nfeats    = get_encoder_features(model, x_tensor)\nfeats_np = feats[0].numpy()\nn_show   = 16\nfig, axes = plt.subplots(2, n_show//2, figsize=(20,6))\nfig.suptitle(f'Feature maps 1er ResBlock encodeur ({feats_np.shape[0]} filtres)\\nChaque filtre detecte un pattern different',\n             fontsize=12, fontweight='bold')\nfor i,ax in enumerate(axes.flatten()):\n    ax.imshow(feats_np[i], cmap='viridis', aspect='auto')\n    ax.set_title(f'Filtre {i+1}', fontsize=8)\n    ax.axis('off')\nplt.tight_layout()\nplt.savefig('feature_maps.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:02:21.346482Z","iopub.execute_input":"2026-06-02T18:02:21.346862Z","iopub.status.idle":"2026-06-02T18:02:22.774861Z","shell.execute_reply.started":"2026-06-02T18:02:21.346828Z","shell.execute_reply":"2026-06-02T18:02:22.774042Z"}},"outputs":[],"execution_count":null},{"id":"a43176a3-38a7-4324-b7b9-b65a40f108c5","cell_type":"markdown","source":"**Feature maps du premier bloc encodeur** :\n\nLes 16 premiers filtres (sur 64) du ResBlock initial montrent que le réseau a appris des **détecteurs spécialisés** :\n- Certains filtres détectent les **arrivées directes** (ligne diagonale principale)\n- D'autres détectent les **réflexions** (arrivées hyperboliques secondaires)\n- Certains sont sensibles aux **gradients d'amplitude** en x ou en temps\n\nCette diversité de features dès le premier niveau confirme que le réseau exploite pleinement les 5 canaux d'entrée et extrait des informations géophysiquement pertinentes.","metadata":{}},{"id":"bc69ce8d-9cda-44cb-bc3c-901268ad0223","cell_type":"markdown","source":"## Section 10 — Inférence sur le Test Set et Soumission\n\nLe test set contient **65 818 fichiers** individuels (1 échantillon par fichier, shape `(5, 1000, 70)`). Pour éviter l'erreur mémoire (OOM), le chargement et l'inférence sont effectués par **chunks de 500 fichiers** avec libération explicite de la mémoire entre chaque chunk.\n\nLe même pipeline de preprocessing (downsampling + normalisation avec les stats du train) est appliqué identiquement sur le test set pour garantir la cohérence train/test.","metadata":{}},{"id":"cb210f00-5eb6-420c-b2cf-17537ce65102","cell_type":"code","source":"# ============================================================\n#  SECTION 10 — INFERENCE TEST SET & SOUMISSION\n#  Chargement par chunks — evite l OOM sur 65818 fichiers\n# ============================================================\ntest_dir   = '/kaggle/input/competitions/waveform-inversion/test'\ntest_files = sorted([\n    os.path.join(test_dir, f)\n    for f in os.listdir(test_dir)\n    if f.endswith('.npy')\n])\nprint(f'Fichiers test : {len(test_files)}')\n\n_sample = np.load(test_files[0])\nprint(f'Shape fichier test (1 fichier) : {_sample.shape}')\n_single_file = (_sample.ndim == 3)\nprint(f'Un echantillon par fichier : {_single_file}')\ndel _sample\n\nCHUNK_SIZE      = 500\nmodel.eval()\nall_predictions = []\nall_test_names  = []\n\nfor chunk_start in range(0, len(test_files), CHUNK_SIZE):\n    chunk_files = test_files[chunk_start : chunk_start + CHUNK_SIZE]\n    chunks = []\n    for f in chunk_files:\n        arr = np.load(f)\n        if arr.ndim == 3:\n            arr = np.expand_dims(arr, axis=0)\n        chunks.append(arr)\n    chunk_data = np.concatenate(chunks, axis=0)\n    del chunks\n    chunk_ds   = downsample_seismic(chunk_data)\n    del chunk_data\n    chunk_norm = chunk_ds.copy()\n    del chunk_ds\n    for i in range(chunk_norm.shape[1]):\n        m, s = x_norm_stats[i]\n        chunk_norm[:, i] = apply_normalize(chunk_norm[:, i], m, s)\n    chunk_tensor = torch.tensor(chunk_norm, dtype=torch.float32)\n    del chunk_norm\n    chunk_loader = DataLoader(\n        torch.utils.data.TensorDataset(chunk_tensor),\n        batch_size=config['batch_size'] * 2,\n        shuffle=False, num_workers=2, pin_memory=config['pin_memory']\n    )\n    chunk_preds = []\n    with torch.no_grad():\n        for (xb,) in chunk_loader:\n            xb = xb.to(device)\n            out = model(xb)\n            chunk_preds.append(out.squeeze(1).cpu().numpy())\n    chunk_preds = np.concatenate(chunk_preds, axis=0)\n    chunk_preds = chunk_preds * saved_y_std + saved_y_mean\n    all_predictions.append(chunk_preds)\n    all_test_names.extend(chunk_files)\n    del chunk_tensor, chunk_preds\n    torch.cuda.empty_cache()\n    n_done = min(chunk_start + CHUNK_SIZE, len(test_files))\n    print(f'Chunk {chunk_start}->{n_done} ({n_done/len(test_files)*100:.1f}%)')\n\npredictions = np.concatenate(all_predictions, axis=0)\nprint(f'Predictions totales : {predictions.shape}')\nprint(f'Plage : [{predictions.min():.0f}, {predictions.max():.0f}] m/s')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:02:22.775940Z","iopub.execute_input":"2026-06-02T18:02:22.776262Z","iopub.status.idle":"2026-06-02T18:30:17.808349Z","shell.execute_reply.started":"2026-06-02T18:02:22.776237Z","shell.execute_reply":"2026-06-02T18:30:17.806377Z"}},"outputs":[],"execution_count":null},{"id":"f8e20af2-0f1d-4a6f-b52b-4e2b4446d998","cell_type":"code","source":"# Section 10.1 — Distribution test vs train\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\nfig.suptitle('Verification coherence : Test vs Train/Val', fontsize=12, fontweight='bold')\naxes[0].hist(targets_denorm.flatten(), bins=60, alpha=0.6, color='steelblue', label='Train/Val reel', density=True)\naxes[0].hist(predictions.flatten(),   bins=60, alpha=0.6, color='salmon',    label='Test predit',   density=True)\naxes[0].set_xlabel('Vitesse (m/s)'); axes[0].set_ylabel('Densite')\naxes[0].set_title('Distribution Train vs Test'); axes[0].legend(fontsize=9); axes[0].grid(True,alpha=0.3)\nim1=axes[1].imshow(predictions.mean(axis=0), cmap='jet', aspect='auto')\naxes[1].set_title('Carte moyenne test set', fontsize=10)\naxes[1].set_xlabel('Position X'); axes[1].set_ylabel('Profondeur')\nplt.colorbar(im1, ax=axes[1], fraction=0.046, label='Vitesse (m/s)')\nim2=axes[2].imshow(predictions.std(axis=0), cmap='hot', aspect='auto')\naxes[2].set_title('Ecart-type predictions test', fontsize=10)\naxes[2].set_xlabel('Position X')\nplt.colorbar(im2, ax=axes[2], fraction=0.046, label='Std (m/s)')\nplt.tight_layout()\nplt.savefig('test_distribution_check.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:30:17.814720Z","iopub.execute_input":"2026-06-02T18:30:17.815686Z","iopub.status.idle":"2026-06-02T18:30:27.516627Z","shell.execute_reply.started":"2026-06-02T18:30:17.815651Z","shell.execute_reply":"2026-06-02T18:30:27.515798Z"}},"outputs":[],"execution_count":null},{"id":"801477b6-85cc-4cb4-9230-9cfca7614f79","cell_type":"markdown","source":"**Vérification de cohérence train/test** :\n- **Distributions** : les distributions train/val (bleu) et test prédit (rouge) se chevauchent bien sur la plage 1500-4500 m/s — pas de distribution shift majeur\n- **Carte moyenne** : le gradient vertical (bleu en surface, rouge en profondeur) confirme que le modèle prédit physiquement des vitesses croissantes avec la profondeur sur l'ensemble du test set\n- **Écart-type élevé** (480-600 m/s) : le modèle prédit des structures géologiques variées selon les échantillons — pas de prédiction uniforme, signe d'une bonne généralisation","metadata":{}},{"id":"25a9ae2b-5fec-42b1-8f7f-c7a0a2038b61","cell_type":"code","source":"# Section 10.2 — Exemples predictions test\n_viz_chunks = []\nfor f in test_files[:4]:\n    arr = np.load(f)\n    if arr.ndim == 3:\n        arr = np.expand_dims(arr, axis=0)\n    _viz_chunks.append(arr)\n_viz_data = np.concatenate(_viz_chunks, axis=0)\n_viz_ds   = downsample_seismic(_viz_data)\n_viz_norm = _viz_ds.copy()\nfor i in range(_viz_norm.shape[1]):\n    m, s = x_norm_stats[i]\n    _viz_norm[:, i] = apply_normalize(_viz_norm[:, i], m, s)\ndel _viz_data, _viz_ds, _viz_chunks\nfig, axes = plt.subplots(2, 4, figsize=(20, 8))\nfig.suptitle('Exemples de predictions sur le test set\\n(pas de ground truth — evaluation visuelle)',\n             fontsize=12, fontweight='bold')\nfor i in range(4):\n    axes[0,i].imshow(_viz_norm[i,0], cmap='RdBu_r', aspect='auto')\n    axes[0,i].set_title(f'Input test #{i+1}', fontsize=9)\n    axes[0,i].set_xlabel('Position', fontsize=8)\n    if i==0: axes[0,i].set_ylabel('Temps', fontsize=8)\n    im=axes[1,i].imshow(predictions[i], cmap='jet', aspect='auto')\n    axes[1,i].set_title(f'Predit #{i+1}\\n[{predictions[i].min():.0f}-{predictions[i].max():.0f}] m/s', fontsize=9)\n    axes[1,i].set_xlabel('Position X', fontsize=8)\n    if i==0: axes[1,i].set_ylabel('Profondeur', fontsize=8)\n    plt.colorbar(im, ax=axes[1,i], fraction=0.046, pad=0.04)\ndel _viz_norm\nplt.tight_layout()\nplt.savefig('test_predictions_examples.png', dpi=120, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:30:27.517749Z","iopub.execute_input":"2026-06-02T18:30:27.518190Z","iopub.status.idle":"2026-06-02T18:30:30.221492Z","shell.execute_reply.started":"2026-06-02T18:30:27.518163Z","shell.execute_reply":"2026-06-02T18:30:30.220582Z"}},"outputs":[],"execution_count":null},{"id":"5d796df1-e4d6-47e7-9cb1-1b4b28a075b3","cell_type":"markdown","source":"**Prédictions test** : sans ground truth disponible, on évalue visuellement la cohérence physique des prédictions. Les 4 exemples montrent des structures géologiques réalistes et variées — couches horizontales, interfaces nettes, variations latérales. Les plages de vitesse prédites (1619-4212 m/s) sont cohérentes avec la distribution du dataset d'entraînement.","metadata":{}},{"id":"73d494fa-b21b-48fb-ba81-25861d78e999","cell_type":"code","source":"# Section 10.3 — Submission CSV\npredictions_sub = predictions[:,:,1::2]\nprint(f'Shape soumission : {predictions_sub.shape}')\nrows=[]\nmax_samples=min(len(all_test_names), predictions_sub.shape[0])\nfor i in range(max_samples):\n    base_name=os.path.basename(all_test_names[i]).replace('.npy','')\n    pred=predictions_sub[i]\n    for y_pos in range(pred.shape[0]):\n        rows.append([f'{base_name}_y_{y_pos}'] + pred[y_pos].tolist())\nx_cols=[f'x_{x}' for x in range(1,70,2)]\nsubmission_df=pd.DataFrame(rows, columns=['oid_ypos']+x_cols)\nassert submission_df.isnull().sum().sum()==0\nassert len(submission_df.columns)==36\nsubmission_df.to_csv('submission.csv', index=False)\nprint(f'submission.csv genere !')\nprint(f'  Lignes   : {len(submission_df):,}')\nprint(f'  Colonnes : {len(submission_df.columns)}')\nprint(f'  Plage    : [{submission_df[x_cols].values.min():.1f}, {submission_df[x_cols].values.max():.1f}] m/s')\nprint(submission_df.head(3))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:30:30.222768Z","iopub.execute_input":"2026-06-02T18:30:30.223123Z","iopub.status.idle":"2026-06-02T18:35:37.609808Z","shell.execute_reply.started":"2026-06-02T18:30:30.223086Z","shell.execute_reply":"2026-06-02T18:35:37.609015Z"}},"outputs":[],"execution_count":null},{"id":"7c184abc-39a5-47f5-b80a-6e41889e67ba","cell_type":"code","source":"# ============================================================\n#  RECAPITULATIF FINAL\n# ============================================================\nprint('='*65)\nprint('RECAPITULATIF FINAL')\nprint('='*65)\nprint(f'Modele         : ResUNet + CBAM Attention (float32)')\nprint(f'Parametres     : {sum(p.numel() for p in model.parameters()):,}')\nprint(f'Epochs         : {len(train_losses)} / {config[\"epochs\"]}')\nprint(f'Temps total    : {sum(epoch_times)/60:.1f} minutes')\nprint(f'Best val loss  : {min(val_losses):.5f}')\nprint(f'RMSE val       : {rmse_all:.2f} m/s')\nprint(f'MAE val        : {mae_all:.2f} m/s')\nprint(f'Correlation    : {corr_all:.4f}')\nprint(f'Soumission     : {len(submission_df):,} lignes | submission.csv')\nprint('='*65)\nprint('Notebook termine. Bonne chance sur le leaderboard !')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-02T18:35:37.610774Z","iopub.execute_input":"2026-06-02T18:35:37.611188Z","iopub.status.idle":"2026-06-02T18:35:37.620425Z","shell.execute_reply.started":"2026-06-02T18:35:37.611160Z","shell.execute_reply":"2026-06-02T18:35:37.619550Z"}},"outputs":[],"execution_count":null},{"id":"e2d667cc-3e15-41b5-b1ff-1195047f46de","cell_type":"markdown","source":"---\n## Récapitulatif Final\n\n### Résultats obtenus\n\n| Métrique | Valeur |\n|---|---|\n| Architecture | ResUNet (33.5M params) + CBAM |\n| Entraînement | 99 epochs × 60s = ~100 min |\n| Meilleure val loss | 0.19989 |\n| **RMSE validation** | **292 m/s** |\n| **MAE validation** | **160 m/s** |\n| **Corrélation** | **0.9262** |\n| Fichiers test traités | 65 818 |\n| Soumission | submission.csv |\n\n### Points forts de l'approche\n\n1. **Stabilité d'entraînement** : float32 pur + warmup + cosine annealing — convergence sur 99 epochs sans instabilité\n2. **Architecture adaptée** : ResUNet 33M params avec CBAM — capacité suffisante pour les structures latérales\n3. **Loss composite** : MSE + L1 + SSIM + Gradient — optimise simultanément précision, robustesse et structure\n4. **Pipeline robuste** : split avant normalisation, chargement par chunks, dénormalisation correcte\n5. **Analyse complète** : EDA, feature maps, attention maps, MC Dropout, profils géophysiciens\n\n### Limites et perspectives\n\n- Les structures géologiques à **ondulations latérales complexes** restent difficiles — limite inhérente à la résolution 70×70\n- Une résolution plus élevée ou un entraînement sur plus de données améliorerait ces cas\n- L'approche **ensemble** (moyenne de plusieurs modèles entraînés avec différents seeds) permettrait de réduire la variance","metadata":{}},{"id":"23030be8-7810-4e2d-aa9e-5e0c63ba54f5","cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}