{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Seismic Analysis Pipeline: SOM + Vision LLM + RAG\n\n**Author:** Yosmely Bermúdez  \n**Dataset:** [OpenFWI — Kaggle Waveform Inversion Competition](https://www.kaggle.com/competitions/waveform-inversion)  \n**Last updated:** 2025\n\n---\n\n## Overview\n\nThis notebook implements a **3-layer universal seismic analysis pipeline** applied to the OpenFWI synthetic dataset:\n\n| Layer | Component | Purpose |\n|-------|-----------|----------|\n| 1 | **Self-Organizing Map (SOM)** | Unsupervised seismic facies classification |\n| 2 | **Vision LLM (GPT-4o mini)** | Automatic geological interpretation of seismic images |\n| 3 | **RAG (arXiv API)** | Scientific report enrichment with relevant literature |\n\n### Datasets used\n- `FlatVel_A` — flat horizontal layers, velocity gradient with depth only\n- `FlatFault_A` — horizontal layers with vertical fault discontinuities\n- `CurveVel_A` — curved layers (future work)\n\n### Cost summary\n- Vision LLM: ~$0.003 per sample (GPT-4o mini, `detail:high`)\n- RAG: ~$0.0003 per sample (text only)\n- **Total for 6 samples analyzed: ~$0.021**\n\n---","metadata":{}},{"cell_type":"markdown","source":"## 0. Environment Setup\n\nStandard Kaggle environment check. Lists all available input files.\nThe `kagglehub.dataset_download` call below is the Kaggle default template — not used in this notebook since data is already available via the competition input path.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\n\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:32:19.920066Z","iopub.execute_input":"2026-06-15T19:32:19.923162Z","iopub.status.idle":"2026-06-15T19:34:08.688979Z","shell.execute_reply.started":"2026-06-15T19:32:19.922976Z","shell.execute_reply":"2026-06-15T19:34:08.687882Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 1. Dataset Exploration\n\n### 1.1 Directory Structure\n\nThe OpenFWI competition data lives under `/kaggle/input/competitions/waveform-inversion/train_samples/`.\nWe explore the top-level structure to understand the file layout before loading anything.","metadata":{}},{"cell_type":"code","source":"import os\n\nbase = '/kaggle/input/competitions/waveform-inversion'\nfor item in os.listdir(base):\n    full_path = os.path.join(base, item)\n    if os.path.isdir(full_path):\n        files = os.listdir(full_path)\n        print(f\"📁 {item}/ — {len(files)} items\")\n        for f in files[:3]:\n            print(f\"     {f}\")\n    else:\n        print(f\"📄 {item}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:08.691509Z","iopub.execute_input":"2026-06-15T19:34:08.691868Z","iopub.status.idle":"2026-06-15T19:34:08.730385Z","shell.execute_reply.started":"2026-06-15T19:34:08.691836Z","shell.execute_reply":"2026-06-15T19:34:08.729174Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 1.2 Shape and Range Inspection\n\nWe iterate over the three dataset families and inspect shapes and value ranges.\n\n**Key observations:**\n- `FlatVel_A` and `CurveVel_A`: organized as `data/` and `model/` subdirectories\n- `FlatFault_A`: files named directly as `seis_*.npy` and `vel_*.npy`\n- All seismic arrays: `(500, 5, 1000, 70)` — 500 samples, 5 sources, 1000 time steps, 70 receivers\n- All velocity maps: `(500, 1, 70, 70)` — 500 samples, 70×70 pixel grid, range 1500–4500 m/s","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport os\n\nbase = '/kaggle/input/competitions/waveform-inversion/train_samples'\n\nfor familia in ['FlatVel_A', 'FlatFault_A', 'CurveVel_A']:\n    path = os.path.join(base, familia)\n    items = os.listdir(path)\n    print(f\"\\n📁 {familia}/\")\n\n    if os.path.isdir(os.path.join(path, items[0])):\n        for sub in items:\n            subpath = os.path.join(path, sub)\n            files = os.listdir(subpath)\n            arr = np.load(os.path.join(subpath, files[0]))\n            print(f\"   📁 {sub}/ — {len(files)} files\")\n            print(f\"        Shape: {arr.shape} | dtype: {arr.dtype}\")\n            print(f\"        Min: {arr.min():.4f} | Max: {arr.max():.4f}\")\n    else:\n        seis = [f for f in items if f.startswith('seis')]\n        vels = [f for f in items if f.startswith('vel')]\n        print(f\"   📄 seis: {len(seis)} files | vel: {len(vels)} files\")\n        if seis:\n            arr = np.load(os.path.join(path, seis[0]))\n            print(f\"   seis shape: {arr.shape} | dtype: {arr.dtype}\")\n            print(f\"   Min: {arr.min():.4f} | Max: {arr.max():.4f}\")\n        if vels:\n            arr = np.load(os.path.join(path, vels[0]))\n            print(f\"   vel  shape: {arr.shape} | dtype: {arr.dtype}\")\n            print(f\"   Min: {arr.min():.4f} | Max: {arr.max():.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:08.732142Z","iopub.execute_input":"2026-06-15T19:34:08.732912Z","iopub.status.idle":"2026-06-15T19:34:13.090922Z","shell.execute_reply.started":"2026-06-15T19:34:08.732872Z","shell.execute_reply":"2026-06-15T19:34:13.089876Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 1.3 Visual EDA — FlatVel_A\n\nWe load `FlatVel_A/data2.npy` and `FlatVel_A/model2.npy` for initial visual exploration.\nThe 7-panel figure shows all 5 seismic sources, the ground truth velocity map, and the velocity distribution histogram.\n\n> **Note:** The variable names `data` and `model` are used here for the initial EDA only. \n> Later in the notebook they are superseded by `data_v`/`model_v` (FlatVel_A) and `data_f`/`model_f` (FlatFault_A) \n> to avoid accidental overwriting between datasets.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\n\nbase = '/kaggle/input/competitions/waveform-inversion/train_samples'\n\ndata  = np.load(f'{base}/FlatVel_A/data/data2.npy')   # (500, 5, 1000, 70)\nmodel = np.load(f'{base}/FlatVel_A/model/model2.npy') # (500, 1, 70, 70)\n\nsample_idx = 0\nfig, axes = plt.subplots(1, 7, figsize=(20, 5))\n\nfor i in range(5):\n    axes[i].imshow(data[sample_idx, i], aspect='auto', cmap='seismic', vmin=-5, vmax=5)\n    axes[i].set_title(f'Source {i+1}')\n    axes[i].set_xlabel('Receivers')\n    axes[i].set_ylabel('Time')\n\naxes[5].imshow(model[sample_idx, 0], aspect='auto', cmap='jet', vmin=1500, vmax=4500)\naxes[5].set_title('Vel Map (GT)')\naxes[5].set_xlabel('X')\naxes[5].set_ylabel('Depth')\n\naxes[6].hist(model[sample_idx, 0].flatten(), bins=50, color='steelblue')\naxes[6].set_title('Velocity distribution')\naxes[6].set_xlabel('m/s')\n\nplt.suptitle('FlatVel_A — Sample 0 | EDA', fontsize=14)\nplt.tight_layout()\nplt.savefig('eda_flatvel.png', dpi=150)\nplt.show()\nplt.close()\n\nprint(f\"Total samples available: {data.shape[0] * 2} per family\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:13.092195Z","iopub.execute_input":"2026-06-15T19:34:13.092478Z","iopub.status.idle":"2026-06-15T19:34:15.11981Z","shell.execute_reply.started":"2026-06-15T19:34:13.092452Z","shell.execute_reply":"2026-06-15T19:34:15.118831Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 2. Layer 1 — Self-Organizing Map (SOM) for FlatVel_A\n\nWe train an unsupervised SOM to classify seismic facies from acoustic attributes.\nThis is analogous to the SEG-Y horizon attribute workflow used in real exploration.\n\n### Feature engineering strategy\n\nFor `FlatVel_A`, velocity only varies with **depth** (horizontal layers). \nWe extract features **per depth row** — one feature vector per depth level per sample:\n\n```\nFor each depth d (0–69):\n    map depth → time index (depth_to_time)\n    extract snapshot at that time across all 70 receivers × 5 sources\n    compute: mean, max, std, energy, zero-crossings  →  25 features\n```\n\nThis gives a `(n_samples × 70, 25)` feature matrix — one row per depth horizon.","metadata":{}},{"cell_type":"code","source":"!pip install minisom -q\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom minisom import MiniSom\nfrom sklearn.preprocessing import StandardScaler","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:15.122725Z","iopub.execute_input":"2026-06-15T19:34:15.12306Z","iopub.status.idle":"2026-06-15T19:34:19.985739Z","shell.execute_reply.started":"2026-06-15T19:34:15.123032Z","shell.execute_reply":"2026-06-15T19:34:19.984009Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 2.1 First attempt: features per receiver column (v1)\n\nInitial approach: extract features per **receiver column X** (vertical traces).\nThis yielded vertical stripe patterns in the facies map — not geologically meaningful,\nsince FlatVel_A has horizontal layering, not vertical.\n\n> ⚠️ **This version was superseded by v2.** Kept here to document the iteration process.","metadata":{}},{"cell_type":"code","source":"base = '/kaggle/input/competitions/waveform-inversion/train_samples/FlatVel_A'\ndata  = np.load(f'{base}/data/data2.npy')\nmodel = np.load(f'{base}/model/model2.npy')\n\nprint(f\"Seismic : {data.shape}\")\nprint(f\"Vel map : {model.shape}\")\n\nN_SAMPLES = 100\n\ndef extract_features_v1(data, model, n_samples=100):\n    \"\"\"Features per receiver column X. Produces vertical stripes — not ideal for FlatVel_A.\"\"\"\n    features_list, vel_list = [], []\n    for s in range(n_samples):\n        seis = data[s]\n        vmap = model[s, 0]\n        for x in range(70):\n            traces = seis[:, :, x]  # (5, 1000)\n            feat = np.concatenate([\n                traces.mean(axis=1), traces.max(axis=1),\n                traces.std(axis=1), (traces**2).mean(axis=1),\n                (np.diff(np.sign(traces)) != 0).sum(axis=1).astype(float)\n            ])\n            features_list.append(feat)\n            vel_list.append(vmap[:, x].mean())\n    return np.array(features_list), np.array(vel_list)\n\nX, y_vel = extract_features_v1(data, model, n_samples=N_SAMPLES)\nprint(f\"Features shape: {X.shape}\")\nprint(f\"Velocity range: {y_vel.min():.0f} – {y_vel.max():.0f} m/s\")\n\nscaler = StandardScaler()\nX_scaled = scaler.fit_transform(X)\n\nSOM_X, SOM_Y = 5, 5\nsom = MiniSom(SOM_X, SOM_Y, input_len=25, sigma=1.5, learning_rate=0.5, random_seed=42)\nsom.random_weights_init(X_scaled)\nsom.train_random(X_scaled, num_iteration=5000, verbose=True)\n\nwinners = np.array([som.winner(x) for x in X_scaled])\nfacies_ids = winners[:, 0] * SOM_Y + winners[:, 1]\nprint(f\"\\nUnique facies found: {len(np.unique(facies_ids))}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:19.987936Z","iopub.execute_input":"2026-06-15T19:34:19.989008Z","iopub.status.idle":"2026-06-15T19:34:22.058215Z","shell.execute_reply.started":"2026-06-15T19:34:19.988927Z","shell.execute_reply":"2026-06-15T19:34:22.057165Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 2.2 Corrected approach: features per depth row (v2)\n\nWe switch to extracting features **per depth row** (horizontal horizons).\nThis aligns with the physical structure of FlatVel_A and produces horizontal band patterns\nthat correctly reflect the layered velocity model.\n\nWe test three SOM grid sizes:\n- 5×5 (25 neurons) — too granular\n- 3×3 (9 neurons) — too coarse\n- **4×4 (16 neurons) — best balance** ✅\n\nThe final model `som2` / `scaler2` (QE ≈ 0.88) is used throughout the rest of the notebook.","metadata":{}},{"cell_type":"code","source":"def extract_features_v2(data, model, n_samples=100):\n    \"\"\"\n    Features per depth row — one vector per depth level per sample.\n    Maps depth index (0-69) to time index (0-999) linearly.\n    Returns X (n_samples*70, 25) and y_vel (n_samples*70,).\n    \"\"\"\n    features_list, vel_list = [], []\n    depth_to_time = np.linspace(0, 999, 70).astype(int)\n\n    for s in range(n_samples):\n        seis = data[s]       # (5, 1000, 70)\n        vmap = model[s, 0]   # (70, 70)\n\n        for d in range(70):\n            t = depth_to_time[d]\n            traces_at_t = seis[:, t, :]  # (5, 70)\n            feat = np.concatenate([\n                traces_at_t.mean(axis=1),\n                traces_at_t.max(axis=1),\n                traces_at_t.std(axis=1),\n                (traces_at_t**2).mean(axis=1),\n                np.abs(np.diff(np.sign(traces_at_t))).sum(axis=1).astype(float)\n            ])\n            features_list.append(feat)\n            vel_list.append(vmap[d, :].mean())\n\n    return np.array(features_list), np.array(vel_list)\n\nX2, y_vel2 = extract_features_v2(data, model, n_samples=100)\nprint(f\"Features shape: {X2.shape}\")\n\nscaler2 = StandardScaler()\nX2_scaled = scaler2.fit_transform(X2)\n\n# Final SOM: 4x4 grid — best balance between granularity and geological meaning\nsom2 = MiniSom(4, 4, input_len=25, sigma=1.2, learning_rate=0.5, random_seed=42)\nsom2.random_weights_init(X2_scaled)\nsom2.train_random(X2_scaled, num_iteration=5000, verbose=True)\n\nfeats_s0 = X2_scaled[0*70 : 1*70]\nwinners_s0 = np.array([som2.winner(x) for x in feats_s0])\nfacies_s0 = winners_s0[:, 0] * 4 + winners_s0[:, 1]\nfacies_map2 = np.tile(facies_s0.reshape(-1, 1), (1, 70))\n\nprint(f\"Unique facies active: {len(np.unique(facies_s0))} of 16 possible\")\n\ncmap = plt.colormaps['tab20'].resampled(16)\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\n\naxes[0].imshow(data[0, 0], aspect='auto', cmap='seismic', vmin=-5, vmax=5)\naxes[0].set_title('Seismic (Source 1)')\n\nim1 = axes[1].imshow(model[0, 0], aspect='auto', cmap='jet', vmin=1500, vmax=4500)\naxes[1].set_title('Vel Map Ground Truth (m/s)')\nplt.colorbar(im1, ax=axes[1])\n\nim2 = axes[2].imshow(facies_map2, aspect='auto', cmap=cmap, vmin=0, vmax=15)\naxes[2].set_title(f'SOM Facies 4x4 ({len(np.unique(facies_s0))} active)')\nplt.colorbar(im2, ax=axes[2])\n\nplt.suptitle('FlatVel_A — Sample 0 | SOM 4x4 (final)', fontsize=14)\nplt.tight_layout()\nplt.savefig('som_facies_4x4.png', dpi=150)\nplt.show()\nplt.close()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:22.060323Z","iopub.execute_input":"2026-06-15T19:34:22.061366Z","iopub.status.idle":"2026-06-15T19:34:24.275814Z","shell.execute_reply.started":"2026-06-15T19:34:22.061324Z","shell.execute_reply":"2026-06-15T19:34:24.274838Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 2.3 Sample selection by velocity variability\n\nTo validate the SOM across diverse geological scenarios, we select 3 samples\nspanning the full range of subsurface complexity:\n\n- **Max variability** (idx=194, std≈1219 m/s) — many distinct velocity layers\n- **Median variability** (idx=69, std≈648 m/s) — typical sample\n- **Min variability** (idx=347, std≈92 m/s) — nearly homogeneous\n\nThis stratified selection tests whether the SOM correctly collapses to fewer facies\nwhen the subsurface is geologically simple.","metadata":{}},{"cell_type":"code","source":"vel_variability = []\nfor i in range(500):\n    vmap = model[i, 0]\n    variability = np.std(vmap.mean(axis=1))\n    vel_variability.append(variability)\n\nvel_variability = np.array(vel_variability)\n\nidx_max = vel_variability.argmax()\nidx_min = vel_variability.argmin()\nidx_med = np.argsort(vel_variability)[len(vel_variability)//2]\n\nprint(f\"Highest variability : idx={idx_max} | std={vel_variability[idx_max]:.1f} m/s\")\nprint(f\"Median variability  : idx={idx_med} | std={vel_variability[idx_med]:.1f} m/s\")\nprint(f\"Lowest variability  : idx={idx_min} | std={vel_variability[idx_min]:.1f} m/s\")\n\nfig, axes = plt.subplots(1, 3, figsize=(15, 5))\nfor ax, idx, label in zip(axes,\n                           [idx_max, idx_med, idx_min],\n                           ['High variability', 'Median', 'Low variability']):\n    im = ax.imshow(model[idx, 0], aspect='auto', cmap='jet', vmin=1500, vmax=4500)\n    ax.set_title(f'{label}\\nidx={idx} | std={vel_variability[idx]:.1f} m/s')\n    ax.set_xlabel('X'); ax.set_ylabel('Depth')\n    plt.colorbar(im, ax=ax)\n\nplt.suptitle('Sample selection for SOM validation', fontsize=13)\nplt.tight_layout()\nplt.show()\nplt.close()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:24.276995Z","iopub.execute_input":"2026-06-15T19:34:24.277652Z","iopub.status.idle":"2026-06-15T19:34:26.398662Z","shell.execute_reply.started":"2026-06-15T19:34:24.277614Z","shell.execute_reply":"2026-06-15T19:34:26.397602Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 2.4 SOM validation on 3 selected samples","metadata":{}},{"cell_type":"code","source":"test_samples = [\n    (194, 'High variability',  1218.9),\n    (69,  'Median',             648.0),\n    (347, 'Low variability',     92.1),\n]\n\nfig, axes = plt.subplots(3, 3, figsize=(18, 15))\ncmap_facies = plt.colormaps['tab20'].resampled(16)\ndepth_to_time = np.linspace(0, 999, 70).astype(int)\n\nfor row, (idx, label, std) in enumerate(test_samples):\n    seis = data[idx]\n    vmap = model[idx, 0]\n\n    feats = []\n    for d in range(70):\n        t = depth_to_time[d]\n        traces_at_t = seis[:, t, :]\n        feat = np.concatenate([\n            traces_at_t.mean(axis=1), traces_at_t.max(axis=1),\n            traces_at_t.std(axis=1), (traces_at_t**2).mean(axis=1),\n            np.abs(np.diff(np.sign(traces_at_t))).sum(axis=1).astype(float)\n        ])\n        feats.append(feat)\n\n    feats_scaled = scaler2.transform(np.array(feats))\n    winners = np.array([som2.winner(x) for x in feats_scaled])\n    facies  = winners[:, 0] * 4 + winners[:, 1]\n    facies_map = np.tile(facies.reshape(-1, 1), (1, 70))\n    n_unique = len(np.unique(facies))\n\n    axes[row, 0].imshow(seis[0], aspect='auto', cmap='seismic', vmin=-5, vmax=5)\n    axes[row, 0].set_title(f'Seismic — {label} (idx={idx})')\n\n    im1 = axes[row, 1].imshow(vmap, aspect='auto', cmap='jet', vmin=1500, vmax=4500)\n    axes[row, 1].set_title(f'Vel Map | std={std:.0f} m/s')\n    plt.colorbar(im1, ax=axes[row, 1])\n\n    im2 = axes[row, 2].imshow(facies_map, aspect='auto', cmap=cmap_facies, vmin=0, vmax=15)\n    axes[row, 2].set_title(f'SOM Facies | {n_unique} active of 16')\n    plt.colorbar(im2, ax=axes[row, 2])\n\n    print(f\"idx={idx:3d} | {label:20s} | active facies: {n_unique} | {np.unique(facies)}\")\n\nplt.suptitle('SOM Validation — 3 extreme samples | FlatVel_A', fontsize=14)\nplt.tight_layout()\nplt.savefig('som_validation_3samples.png', dpi=150)\nplt.show()\nplt.close()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:26.402393Z","iopub.execute_input":"2026-06-15T19:34:26.403403Z","iopub.status.idle":"2026-06-15T19:34:30.039487Z","shell.execute_reply.started":"2026-06-15T19:34:26.403362Z","shell.execute_reply":"2026-06-15T19:34:30.038485Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 2.5 Homogeneity threshold\n\nFor nearly-homogeneous samples (velocity std < 200 m/s), the SOM may assign\nmultiple facies to what is effectively one geological unit.\nWe add a threshold: if `std(vel_profile) < threshold`, assign a single facies.","metadata":{}},{"cell_type":"code","source":"def predict_facies(seis, vmap, som, scaler, threshold_std=200):\n    \"\"\"Predict SOM facies with homogeneity threshold.\"\"\"\n    depth_to_time = np.linspace(0, 999, 70).astype(int)\n    sample_std = np.std(vmap.mean(axis=1))\n\n    feats = []\n    for d in range(70):\n        t = depth_to_time[d]\n        traces_at_t = seis[:, t, :]\n        feat = np.concatenate([\n            traces_at_t.mean(axis=1), traces_at_t.max(axis=1),\n            traces_at_t.std(axis=1), (traces_at_t**2).mean(axis=1),\n            np.abs(np.diff(np.sign(traces_at_t))).sum(axis=1).astype(float)\n        ])\n        feats.append(feat)\n\n    feats_scaled = scaler.transform(np.array(feats))\n\n    if sample_std < threshold_std:\n        facies = np.zeros(70, dtype=int)\n        print(f\"  → Homogeneous (std={sample_std:.1f}) — single facies assigned\")\n    else:\n        winners = np.array([som.winner(x) for x in feats_scaled])\n        facies  = winners[:, 0] * 4 + winners[:, 1]\n        print(f\"  → std={sample_std:.1f} m/s | {len(np.unique(facies))} active facies\")\n\n    return np.tile(facies.reshape(-1, 1), (1, 70))\n\nfor idx, label, std in test_samples:\n    print(f\"\\n{label} (idx={idx})\")\n    fm = predict_facies(data[idx], model[idx, 0], som2, scaler2)\n    print(f\"  Final unique facies: {len(np.unique(fm))}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:30.040905Z","iopub.execute_input":"2026-06-15T19:34:30.041309Z","iopub.status.idle":"2026-06-15T19:34:30.083241Z","shell.execute_reply.started":"2026-06-15T19:34:30.041268Z","shell.execute_reply":"2026-06-15T19:34:30.082157Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 2.6 Lithological interpretation\n\nWe map SOM facies to lithological classes using the velocity ranges defined in `VEL_RANGES`.\nEach facies is assigned to a lithology class based on its mean velocity.\n\nThe `vel_to_litho` function implements this mapping.","metadata":{}},{"cell_type":"code","source":"# Geological velocity ranges — physically meaningful classification\nVEL_RANGES = [\n    (1500, 2000, 'Soft sediments',      '#4575b4'),\n    (2000, 2500, 'Compacted sediments', '#74add1'),\n    (2500, 3000, 'Sedimentary rock',    '#abd9e9'),\n    (3000, 3500, 'Hard rock',           '#fdae61'),\n    (3500, 4500, 'Very hard rock',      '#d73027'),\n]\n\ndef vel_to_litho(vel):\n    \"\"\"Map velocity (m/s) to lithology class index, label, and color.\"\"\"\n    for i, (vmin, vmax, label, color) in enumerate(VEL_RANGES):\n        if vmin <= vel < vmax:\n            return i, label, color\n    return len(VEL_RANGES)-1, VEL_RANGES[-1][2], VEL_RANGES[-1][3]\n\nprint(\"vel_to_litho() defined ✅\")\nprint(\"VEL_RANGES:\")\nfor vmin, vmax, label, color in VEL_RANGES:\n    print(f\"  {vmin}–{vmax} m/s → {label} ({color})\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:30.084532Z","iopub.execute_input":"2026-06-15T19:34:30.084827Z","iopub.status.idle":"2026-06-15T19:34:30.094847Z","shell.execute_reply.started":"2026-06-15T19:34:30.0848Z","shell.execute_reply":"2026-06-15T19:34:30.092696Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 3. Layer 1 — SOM for FlatFault_A\n\nFlatFault_A introduces **vertical fault discontinuities** on top of horizontal layering.\nThe velocity model is no longer purely 1D — faults create lateral variations.\n\n### Key difference from FlatVel_A\n\n| | FlatVel_A | FlatFault_A |\n|--|-----------|-------------|\n| SOM type | 1D (per depth row) | **2D (per pixel)** |\n| Features | 25 (5 sources × 5 attrs) | **7 (point-local attrs)** |\n| Grid | 4×4 | 4×4 |\n| Input len | 25 | **7** |\n| QE | 0.88 | **0.44** |\n\nThe 2D pixel-wise approach captures lateral discontinuities that the 1D row-wise approach would miss.","metadata":{}},{"cell_type":"code","source":"import os\n\nbase_fault = '/kaggle/input/competitions/waveform-inversion/train_samples/FlatFault_A'\nfiles = os.listdir(base_fault)\nseis_files = sorted([f for f in files if f.startswith('seis')])\nvel_files  = sorted([f for f in files if f.startswith('vel')])\n\nprint(\"Seismic files:\", seis_files)\nprint(\"Vel files:    \", vel_files)\n\ndata_f  = np.load(f'{base_fault}/{seis_files[0]}')\nmodel_f = np.load(f'{base_fault}/{vel_files[0]}')\n\nprint(f\"\\nSeismic shape : {data_f.shape}\")\nprint(f\"Vel map shape : {model_f.shape}\")\n\n# Quick visual — 3 samples to see fault variation\nfig, axes = plt.subplots(2, 3, figsize=(18, 8))\nfor col, idx in enumerate([0, 100, 400]):\n    axes[0, col].imshow(data_f[idx, 0], aspect='auto', cmap='seismic', vmin=-5, vmax=5)\n    axes[0, col].set_title(f'Seismic idx={idx}')\n\n    im1 = axes[1, col].imshow(model_f[idx, 0], aspect='auto', cmap='jet', vmin=1500, vmax=4500)\n    axes[1, col].set_title(f'Vel Map idx={idx}')\n    plt.colorbar(im1, ax=axes[1, col], label='m/s')\n\nplt.suptitle('FlatFault_A — Initial exploration (idx 0, 100, 400)', fontsize=13)\nplt.tight_layout()\nplt.show()\nplt.close()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:30.096332Z","iopub.execute_input":"2026-06-15T19:34:30.096675Z","iopub.status.idle":"2026-06-15T19:34:31.699622Z","shell.execute_reply.started":"2026-06-15T19:34:30.096637Z","shell.execute_reply":"2026-06-15T19:34:31.69863Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 3.1 Feature extraction — 2D pixel-wise\n\nFor each pixel `(depth d, lateral x)` in the 70×70 velocity map, we extract\n7 local features from the 5-source seismic snapshot at the corresponding time:\n\n```\nvals = seis[:, t, x]  # 5 values (one per source)\nfeatures = [mean, max, min, std, energy, range, zero_crossings]\n```\n\nThis produces `n_samples × 70 × 70 = 245,000` feature vectors for 50 training samples.","metadata":{}},{"cell_type":"code","source":"def extract_features_2d(data_arr, model_arr, n_samples=100):\n    \"\"\"\n    Pixel-wise feature extraction for FlatFault_A.\n    Returns X (n*70*70, 7), y (n*70*70,), coords list.\n    \"\"\"\n    depth_to_time = np.linspace(0, 999, 70).astype(int)\n    features_list, vel_list, coords = [], [], []\n\n    for s in range(n_samples):\n        seis = data_arr[s]\n        vmap = model_arr[s, 0]\n        for d in range(70):\n            t = depth_to_time[d]\n            for x in range(70):\n                vals = seis[:, t, x]  # (5,)\n                feat = np.array([\n                    vals.mean(), vals.max(), vals.min(), vals.std(),\n                    (vals**2).mean(),\n                    vals.max() - vals.min(),\n                    np.abs(np.diff(np.sign(vals))).sum()\n                ])\n                features_list.append(feat)\n                vel_list.append(vmap[d, x])\n                coords.append((s, d, x))\n\n    return np.array(features_list), np.array(vel_list), coords\n\nprint(\"Extracting 2D features... (1-2 min)\")\nX_f, y_f, coords_f = extract_features_2d(data_f, model_f, n_samples=50)\nprint(f\"Features shape : {X_f.shape}\")\nprint(f\"Vel range      : {y_f.min():.0f} – {y_f.max():.0f} m/s\")\n\nscaler_f = StandardScaler()\nX_f_scaled = scaler_f.fit_transform(X_f)\n\nsom_f = MiniSom(4, 4, input_len=7, sigma=1.2, learning_rate=0.5, random_seed=42)\nsom_f.random_weights_init(X_f_scaled)\nprint(\"Training FlatFault_A SOM...\")\nsom_f.train_random(X_f_scaled, num_iteration=10000, verbose=True)\nprint(f\"Quantization error: {som_f.quantization_error(X_f_scaled):.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:31.700939Z","iopub.execute_input":"2026-06-15T19:34:31.701938Z","iopub.status.idle":"2026-06-15T19:34:48.839761Z","shell.execute_reply.started":"2026-06-15T19:34:31.701893Z","shell.execute_reply":"2026-06-15T19:34:48.838637Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 3.2 Visualization — FlatFault_A lithology maps\n\nThree representative samples:\n- **idx=0**: No fault (flat horizontal layers only)\n- **idx=100**: Central fault (clear vertical discontinuity)\n- **idx=400**: Complex fault geometry\n\nThe `plot_fault_sample` function generates a 4-panel figure:\nseismic | vel map | SOM lithology 2D | geological legend","metadata":{}},{"cell_type":"code","source":"from matplotlib.colors import ListedColormap, BoundaryNorm\n\ndef plot_fault_sample(idx, label, data_arr, model_arr, som, scaler):\n    seis = data_arr[idx]\n    vmap = model_arr[idx, 0]\n    depth_to_time = np.linspace(0, 999, 70).astype(int)\n\n    feats = []\n    for d in range(70):\n        t = depth_to_time[d]\n        for x in range(70):\n            vals = seis[:, t, x]\n            feat = np.array([\n                vals.mean(), vals.max(), vals.min(), vals.std(),\n                (vals**2).mean(), vals.max()-vals.min(),\n                np.abs(np.diff(np.sign(vals))).sum()\n            ])\n            feats.append(feat)\n\n    feats_scaled = scaler.transform(np.array(feats))\n    winners      = np.array([som.winner(x) for x in feats_scaled])\n    facies_flat  = winners[:, 0] * 4 + winners[:, 1]\n    facies_map   = facies_flat.reshape(70, 70)\n\n    facies_vel = {}\n    for d in range(70):\n        for x in range(70):\n            f = facies_map[d, x]\n            if f not in facies_vel:\n                facies_vel[f] = []\n            facies_vel[f].append(vmap[d, x])\n    facies_vel = {f: np.mean(v) for f, v in facies_vel.items()}\n\n    litho_map = np.zeros((70, 70), dtype=int)\n    for f, vel in facies_vel.items():\n        clase, _, _ = vel_to_litho(vel)\n        litho_map[facies_map == f] = clase\n\n    colors     = [r[3] for r in VEL_RANGES]\n    cmap_litho = ListedColormap(colors)\n    norm_litho = BoundaryNorm(boundaries=range(len(VEL_RANGES)+1), ncolors=len(VEL_RANGES))\n\n    fig = plt.figure(figsize=(22, 6))\n    gs  = fig.add_gridspec(1, 4, width_ratios=[3, 3, 3, 2.5], wspace=0.38)\n    ax0, ax1, ax2, ax3 = [fig.add_subplot(gs[i]) for i in range(4)]\n\n    ax0.imshow(seis[0], aspect='auto', cmap='seismic', vmin=-5, vmax=5)\n    ax0.set_title('Seismic (Source 1)', fontsize=11, pad=8)\n    ax0.set_xlabel('Receivers'); ax0.set_ylabel('Time')\n\n    im1 = ax1.imshow(vmap, aspect='auto', cmap='jet', vmin=1500, vmax=4500)\n    ax1.set_title('Vel Map Ground Truth', fontsize=11, pad=8)\n    ax1.set_xlabel('X'); ax1.set_ylabel('Depth')\n    plt.colorbar(im1, ax=ax1, label='m/s', fraction=0.046, pad=0.04)\n\n    clases_unicas = np.unique(litho_map)\n    im2 = ax2.imshow(litho_map, aspect='auto', cmap=cmap_litho, norm=norm_litho)\n    ax2.set_title(f'SOM Lithology 2D\\n{len(clases_unicas)} class(es) detected', fontsize=11, pad=8)\n    ax2.set_xlabel('X'); ax2.set_ylabel('Depth')\n\n    ax3.axis('off')\n    ax3.set_title('Lithology → Velocity', fontsize=11, pad=10)\n    present = set(litho_map.flatten())\n    y_pos   = np.linspace(0.88, 0.08, len(VEL_RANGES))\n\n    for i, (vmin, vmax, desc, color) in enumerate(VEL_RANGES):\n        y = y_pos[i]\n        alpha = 1.0 if i in present else 0.25\n        ax3.add_patch(plt.Rectangle((0.02, y-0.045), 0.14, 0.08,\n                                     color=color, alpha=alpha, transform=ax3.transAxes, clip_on=False))\n        ax3.text(0.20, y+0.015, f'{vmin}–{vmax} m/s', transform=ax3.transAxes,\n                 fontsize=9.5, va='center', fontfamily='monospace', alpha=alpha)\n        ax3.text(0.20, y-0.022, desc, transform=ax3.transAxes,\n                 fontsize=8.5, va='center', color='#444444', style='italic', alpha=alpha)\n\n    ax3.text(0.02, 0.0, '* faded = absent in this sample',\n             transform=ax3.transAxes, fontsize=7.5, color='gray', style='italic')\n\n    std_val = np.std(vmap.mean(axis=1))\n    fig.suptitle(f'{label}   |   idx={idx}   |   vel std = {std_val:.0f} m/s', fontsize=13, y=1.02)\n    plt.savefig(f'fault_litho_{idx}.png', dpi=150, bbox_inches='tight')\n    plt.show()\n    plt.close()\n    print(f\"✅ Saved: fault_litho_{idx}.png\")\n\nfor idx, label in [(0,   'FlatFault_A — no fault'),\n                   (100, 'FlatFault_A — central fault'),\n                   (400, 'FlatFault_A — complex fault')]:\n    plot_fault_sample(idx, label, data_f, model_f, som_f, scaler_f)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:48.844432Z","iopub.execute_input":"2026-06-15T19:34:48.844811Z","iopub.status.idle":"2026-06-15T19:34:53.451064Z","shell.execute_reply.started":"2026-06-15T19:34:48.84478Z","shell.execute_reply":"2026-06-15T19:34:53.449809Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 4. Layer 2 — Vision LLM: Connectivity Tests\n\nBefore committing to a Vision LLM provider, we tested several free options from Kaggle.\nThis section documents the failed attempts — included for transparency and reproducibility.\n\n### HuggingFace token setup","metadata":{}},{"cell_type":"code","source":"!pip install huggingface_hub -q\n\nfrom kaggle_secrets import UserSecretsClient\nsecrets = UserSecretsClient()\nHF_TOKEN = secrets.get_secret(\"HF_TOKEN\")\n\nfrom huggingface_hub import whoami\ninfo = whoami(token=HF_TOKEN)\nprint(f\"✅ Connected as: {info['name']}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:53.45317Z","iopub.execute_input":"2026-06-15T19:34:53.453479Z","iopub.status.idle":"2026-06-15T19:34:58.322828Z","shell.execute_reply.started":"2026-06-15T19:34:53.453452Z","shell.execute_reply":"2026-06-15T19:34:58.321644Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 4.1 Connectivity test — which endpoints are reachable from Kaggle?\n\nKaggle's network blocks `api-inference.huggingface.co` but allows `huggingface.co` (Spaces).\nWe probe each endpoint before attempting any API calls.","metadata":{}},{"cell_type":"code","source":"import requests\n\nurls = [\n    'https://api-inference.huggingface.co',  # ❌ blocked by Kaggle\n    'https://huggingface.co',               # ✅ accessible\n    'https://www.google.com',               # ✅ accessible\n    'https://api.openai.com',               # ✅ accessible\n    'https://export.arxiv.org',             # ✅ accessible\n]\n\nfor url in urls:\n    try:\n        r = requests.get(url, timeout=5)\n        print(f\"✅ {url} → {r.status_code}\")\n    except Exception as e:\n        print(f\"❌ {url} → {str(e)[:60]}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:58.324748Z","iopub.execute_input":"2026-06-15T19:34:58.325116Z","iopub.status.idle":"2026-06-15T19:34:58.501403Z","shell.execute_reply.started":"2026-06-15T19:34:58.325067Z","shell.execute_reply":"2026-06-15T19:34:58.500233Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 4.2 Failed attempts — HuggingFace Gradio Spaces\n\nWe tried four Vision models via Gradio Spaces. All failed for different reasons:\n\n| Model | Error | Reason |\n|-------|-------|--------|\n| `vikhyatk/moondream2` | upstream error | Space overloaded |\n| `microsoft/Florence-2` | 401 | Auth required |\n| `qnguyen3/nanollava` | RUNTIME_ERROR | Space broken |\n| `Salesforce/BLIP2` | 429 | Rate limited |\n\n> ⚠️ The cells below are kept for documentation. They will fail — this is expected.","metadata":{}},{"cell_type":"code","source":"from gradio_client import Client, handle_file\n\n# Attempt 1 — Moondream2\n# Result: upstream error (Space overloaded or endpoint changed)\ntry:\n    c = Client(\"vikhyatk/moondream2\")\n    print(\"Connected. Endpoints:\", c.view_api())\nexcept Exception as e:\n    print(f\"❌ Moondream2: {str(e)[:120]}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:58.503278Z","iopub.execute_input":"2026-06-15T19:34:58.50367Z","iopub.status.idle":"2026-06-15T19:34:58.895861Z","shell.execute_reply.started":"2026-06-15T19:34:58.50363Z","shell.execute_reply":"2026-06-15T19:34:58.894788Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Attempt 2 — Florence-2, nanollava, BLIP2 via Gradio\n# Results: 401, RUNTIME_ERROR, 429 respectively\nspaces_to_try = [\n    (\"microsoft/Florence-2-large\", \"/process_image\"),\n    (\"qnguyen3/nanollava\",          \"/chat\"),\n    (\"Salesforce/BLIP2\",            \"/predict\"),\n]\n\nfor space, api in spaces_to_try:\n    try:\n        c = Client(space)\n        print(f\"✅ Connected to {space}\")\n        break\n    except Exception as e:\n        print(f\"❌ {space}: {str(e)[:80]}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:58.897287Z","iopub.execute_input":"2026-06-15T19:34:58.897653Z","iopub.status.idle":"2026-06-15T19:34:59.202289Z","shell.execute_reply.started":"2026-06-15T19:34:58.897615Z","shell.execute_reply":"2026-06-15T19:34:59.201176Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Attempt 3 — BLIP2 via direct HTTP (bypassing gradio_client websocket bug)\n# Result: 429 Too Many Requests — public Space rate limited\nimport requests, json\n\n# NOTE: img_b64 must be defined from a previous cell\nurl = \"https://salesforce-blip2.hf.space/run/predict\"\npayload = {\n    \"fn_index\": 2,\n    \"data\": [\n        f\"data:image/png;base64,{img_b64}\",\n        \"What geological structures are visible in this seismic velocity map?\"\n    ]\n}\n\nr = requests.post(url, json=payload, timeout=90)\nprint(f\"Status: {r.status_code}\")\nif r.status_code == 200:\n    print(r.json()['data'])\nelse:\n    print(r.text[:300])  # Expected: 429 or HTML error page","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:59.203826Z","iopub.execute_input":"2026-06-15T19:34:59.204216Z","iopub.status.idle":"2026-06-15T19:34:59.257583Z","shell.execute_reply.started":"2026-06-15T19:34:59.204176Z","shell.execute_reply":"2026-06-15T19:34:59.256629Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 4.3 Solution — OpenAI GPT-4o mini Vision\n\n`api.openai.com` is accessible from Kaggle. GPT-4o mini with `detail:high` provides\ngenuine image understanding at low cost (~$0.003/call).\n\n**Cost analysis:**\n- `detail:low` → ~2,988 tokens → $0.00072 — too generic, doesn't actually see the image\n- `detail:high` + large PNG → ~48,409 tokens → $0.00737 — accurate but expensive\n- `detail:high` + JPEG 1200px → **~20,000 tokens → $0.003** ✅ — best balance\n\n**Key optimization:** convert matplotlib figure → PNG → PIL resize → JPEG quality=82\nbefore encoding to base64. This reduces tokens 2.4× vs. sending raw PNG.","metadata":{}},{"cell_type":"code","source":"# Connectivity test\nimport requests\ntry:\n    r = requests.get(\"https://api.openai.com\", timeout=5)\n    print(f\"✅ OpenAI accessible → {r.status_code}\")\nexcept Exception as e:\n    print(f\"❌ {str(e)[:80]}\")\n\nOPENAI_KEY = secrets.get_secret(\"OPENAI_API_KEY\")\nprint(\"OPENAI_KEY loaded ✅\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:59.2589Z","iopub.execute_input":"2026-06-15T19:34:59.259376Z","iopub.status.idle":"2026-06-15T19:34:59.353758Z","shell.execute_reply.started":"2026-06-15T19:34:59.259334Z","shell.execute_reply":"2026-06-15T19:34:59.35274Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 5. Layer 2 — Vision LLM: Production Functions\n\n### 5.1 make_seismic_plot\n\nGenerates the 3-panel figure sent to the Vision LLM:\n- LEFT: seismic waveform (source 1)\n- CENTER: velocity map (ground truth)\n- RIGHT: SOM lithology classification","metadata":{}},{"cell_type":"code","source":"from matplotlib.colors import ListedColormap, BoundaryNorm\n\ndef make_seismic_plot(idx, data_arr, model_arr, litho_map):\n    \"\"\"\n    Generates 3-panel seismic figure for Vision LLM input.\n    Returns matplotlib figure (not shown, not saved).\n    \"\"\"\n    colors     = [r[3] for r in VEL_RANGES]\n    cmap_litho = ListedColormap(colors)\n    norm_litho = BoundaryNorm(boundaries=range(len(VEL_RANGES)+1), ncolors=len(VEL_RANGES))\n\n    fig = plt.figure(figsize=(16, 5))\n    gs  = fig.add_gridspec(1, 3, wspace=0.35)\n    ax0, ax1, ax2 = [fig.add_subplot(gs[i]) for i in range(3)]\n\n    ax0.imshow(data_arr[idx, 0], aspect='auto', cmap='seismic', vmin=-5, vmax=5)\n    ax0.set_title('Seismic Data (Source 1)', fontsize=11)\n    ax0.set_xlabel('Receivers'); ax0.set_ylabel('Time')\n\n    im1 = ax1.imshow(model_arr[idx, 0], aspect='auto', cmap='jet', vmin=1500, vmax=4500)\n    ax1.set_title('Velocity Map (m/s)', fontsize=11)\n    ax1.set_xlabel('X'); ax1.set_ylabel('Depth')\n    plt.colorbar(im1, ax=ax1, label='m/s', fraction=0.046)\n\n    im2 = ax2.imshow(litho_map, aspect='auto', cmap=cmap_litho, norm=norm_litho)\n    ax2.set_title('SOM Lithology Classification', fontsize=11)\n    ax2.set_xlabel('X'); ax2.set_ylabel('Depth')\n\n    legend_labels = [f'{vmin}–{vmax} m/s: {desc}' for vmin, vmax, desc, _ in VEL_RANGES]\n    handles = [plt.Rectangle((0,0),1,1, color=c) for _,_,_,c in VEL_RANGES]\n    ax2.legend(handles, legend_labels, loc='lower right', fontsize=7, framealpha=0.85)\n\n    fig.suptitle(f'Seismic Analysis — idx={idx}', fontsize=13)\n    return fig\n\nprint(\"make_seismic_plot() defined ✅\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:59.355344Z","iopub.execute_input":"2026-06-15T19:34:59.355734Z","iopub.status.idle":"2026-06-15T19:34:59.37288Z","shell.execute_reply.started":"2026-06-15T19:34:59.355697Z","shell.execute_reply":"2026-06-15T19:34:59.371879Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 5.2 compute_litho_map — FlatFault_A\n\nComputes the SOM lithology map for any FlatFault_A sample.\nUses the 2D pixel-wise SOM (`som_f`, `scaler_f`).","metadata":{}},{"cell_type":"code","source":"def compute_litho_map(idx, data_arr, model_arr, som, scaler):\n    \"\"\"\n    Computes 2D lithology map for FlatFault_A using pixel-wise SOM.\n    Returns: litho_map (70,70), facies_map (70,70), vmap (70,70).\n    \"\"\"\n    seis = data_arr[idx]\n    vmap = model_arr[idx, 0]\n    depth_to_time = np.linspace(0, 999, 70).astype(int)\n\n    feats = []\n    for d in range(70):\n        t = depth_to_time[d]\n        for x in range(70):\n            vals = seis[:, t, x]\n            feats.append(np.array([\n                vals.mean(), vals.max(), vals.min(), vals.std(),\n                (vals**2).mean(), vals.max()-vals.min(),\n                np.abs(np.diff(np.sign(vals))).sum()\n            ]))\n\n    feats_scaled = scaler.transform(np.array(feats))\n    winners      = np.array([som.winner(f) for f in feats_scaled])\n    facies_map   = (winners[:, 0] * 4 + winners[:, 1]).reshape(70, 70)\n\n    facies_vel = {}\n    for d in range(70):\n        for x in range(70):\n            f = facies_map[d, x]\n            if f not in facies_vel:\n                facies_vel[f] = []\n            facies_vel[f].append(vmap[d, x])\n    facies_vel = {f: np.mean(v) for f, v in facies_vel.items()}\n\n    litho_map = np.zeros((70, 70), dtype=int)\n    for f, vel in facies_vel.items():\n        clase, _, _ = vel_to_litho(vel)\n        litho_map[facies_map == f] = clase\n\n    return litho_map, facies_map, vmap\n\nprint(\"compute_litho_map() defined ✅\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:59.374372Z","iopub.execute_input":"2026-06-15T19:34:59.374748Z","iopub.status.idle":"2026-06-15T19:34:59.399792Z","shell.execute_reply.started":"2026-06-15T19:34:59.37471Z","shell.execute_reply":"2026-06-15T19:34:59.398362Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 5.3 detect_fault_from_som\n\nNumerically detects fault presence from the SOM lithology map.\nComputes mean lateral difference across the lithology map — \nhigh scores indicate abrupt horizontal changes consistent with faulting.\n\nThis result is injected into the Vision LLM prompt as prior context,\nimproving fault detection accuracy significantly (especially for no-fault samples).","metadata":{}},{"cell_type":"code","source":"def detect_fault_from_som(litho_map):\n    \"\"\"\n    Detects lateral discontinuity in SOM lithology map.\n    Returns: has_fault (bool), fault_score (float), fault_x (int).\n    \"\"\"\n    lateral_diff = np.diff(litho_map, axis=1)  # (70, 69)\n    fault_score  = np.abs(lateral_diff).mean()\n    col_changes  = np.abs(lateral_diff).sum(axis=0)\n    fault_x      = col_changes.argmax()\n    has_fault    = fault_score > 0.3  # empirical threshold\n\n    return {\n        'has_fault':   has_fault,\n        'fault_score': round(float(fault_score), 3),\n        'fault_x':     int(fault_x)\n    }\n\nprint(\"detect_fault_from_som() defined ✅\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:59.401018Z","iopub.execute_input":"2026-06-15T19:34:59.40139Z","iopub.status.idle":"2026-06-15T19:34:59.424524Z","shell.execute_reply.started":"2026-06-15T19:34:59.401351Z","shell.execute_reply":"2026-06-15T19:34:59.423428Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 5.4 analyze_seismic_vision — main Vision LLM function\n\nFull pipeline:\n1. matplotlib figure → PNG → PIL resize (max 1200px) → JPEG q=82 → base64\n2. Build prompt with SOM fault context injected\n3. Call GPT-4o mini with `detail:high`\n4. Return interpretation text + token usage + cost","metadata":{}},{"cell_type":"code","source":"from PIL import Image\nimport io, base64, requests\n\ndef analyze_seismic_vision(fig, openai_key, litho_map=None, prompt=None, max_width=1200):\n    \"\"\"\n    Converts matplotlib figure to optimized JPEG and calls GPT-4o mini Vision.\n    \n    Args:\n        fig: matplotlib Figure object\n        openai_key: OpenAI API key string\n        litho_map: optional (70,70) array — enables SOM fault context injection\n        prompt: optional custom prompt (uses default if None)\n        max_width: maximum image width in pixels before downsampling\n    \n    Returns:\n        dict with keys: text, prompt_tokens, completion_tokens, cost_usd\n    \"\"\"\n    # Step 1 — figure → optimized JPEG base64\n    buf_raw = io.BytesIO()\n    fig.savefig(buf_raw, format='png', dpi=120, bbox_inches='tight')\n    buf_raw.seek(0)\n\n    img_pil = Image.open(buf_raw)\n    w, h = img_pil.size\n    if w > max_width:\n        img_pil = img_pil.resize((max_width, int(h * max_width / w)), Image.LANCZOS)\n\n    buf = io.BytesIO()\n    img_pil.convert('RGB').save(buf, format='JPEG', quality=82)\n    buf.seek(0)\n    img_b64 = base64.b64encode(buf.read()).decode('utf-8')\n    plt.close(fig)  # avoid memory warning for >20 open figures\n\n    # Step 2 — SOM context injection\n    som_context = \"\"\n    if litho_map is not None:\n        fault_info = detect_fault_from_som(litho_map)\n        som_context = (\n            f\"CONTEXT FROM SOM ANALYSIS: \"\n            f\"{'Fault likely present' if fault_info['has_fault'] else 'NO fault detected'} \"\n            f\"(discontinuity score={fault_info['fault_score']}, \"\n            f\"most likely at lateral index {fault_info['fault_x']}).\\n\\n\"\n        )\n\n    # Step 3 — prompt\n    if prompt is None:\n        prompt = f\"\"\"{som_context}You are an expert geophysicist. Analyze this image carefully.\nTHREE panels: LEFT=seismic waveform, CENTER=velocity map (m/s), RIGHT=SOM lithology.\n\nAnswer based ONLY on what you actually see. Be conservative about faults —\nonly report one if there is a CLEAR offset or discontinuity, not just noise.\nThis is a synthetic dataset (OpenFWI FlatFault_A): some samples have faults, others do not.\n\nNOTE: Axes show pixel indices (0–70), NOT metric units. Use 'depth index' / 'lateral index'.\n\n1. FAULT: Clear discontinuity in CENTER panel? YES/NO. If yes: lateral index and depth index.\n2. VELOCITY RANGE: Min and max from colorbar.\n3. LAYERS: How many distinct horizontal layers? Approximate depth index ranges.\n4. SOM MAP: Dominant colors in top half vs bottom half of RIGHT panel.\n5. SEISMIC: Reflection hyperbolas or offsets in LEFT panel?\"\"\"\n\n    # Step 4 — API call\n    r = requests.post(\n        \"https://api.openai.com/v1/chat/completions\",\n        headers={\"Authorization\": f\"Bearer {openai_key}\", \"Content-Type\": \"application/json\"},\n        json={\n            \"model\": \"gpt-4o-mini\",\n            \"max_tokens\": 500,\n            \"messages\": [{\n                \"role\": \"user\",\n                \"content\": [\n                    {\"type\": \"image_url\",\n                     \"image_url\": {\"url\": f\"data:image/jpeg;base64,{img_b64}\", \"detail\": \"high\"}},\n                    {\"type\": \"text\", \"text\": prompt}\n                ]\n            }]\n        },\n        timeout=60\n    )\n\n    if r.status_code != 200:\n        return {\"error\": r.text[:200]}\n\n    resp  = r.json()\n    usage = resp['usage']\n    cost  = (usage['prompt_tokens'] * 0.00015 + usage['completion_tokens'] * 0.0006) / 1000\n\n    return {\n        \"text\":              resp['choices'][0]['message']['content'],\n        \"prompt_tokens\":     usage['prompt_tokens'],\n        \"completion_tokens\": usage['completion_tokens'],\n        \"cost_usd\":          round(cost, 5)\n    }\n\nprint(\"analyze_seismic_vision() ready ✅\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:59.425944Z","iopub.execute_input":"2026-06-15T19:34:59.426342Z","iopub.status.idle":"2026-06-15T19:34:59.448399Z","shell.execute_reply.started":"2026-06-15T19:34:59.426301Z","shell.execute_reply":"2026-06-15T19:34:59.447412Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 6. Layer 3 — RAG with arXiv\n\nThe third layer enriches the Vision LLM interpretation with scientific literature.\n\n### Architecture\n\n```\nVision LLM output text\n    → extract_geo_keywords()  (keyword detection from text)\n    → search_arxiv()          (arXiv Atom API, no key required)\n    → build context string    (title + abstract, 400 chars each)\n    → GPT-4o mini (text only) (synthesize report with citations)\n    → geological report\n```\n\n**Why arXiv over a local FAISS index?**\n- No installation, no embedding model, no disk space\n- Always up-to-date papers\n- Free, no API key\n- Sufficient for this prototype — FAISS would be Layer 3 v2","metadata":{}},{"cell_type":"code","source":"# Connectivity test\nimport requests\ntry:\n    r = requests.get(\"https://export.arxiv.org/api/query?search_query=seismic&max_results=1\", timeout=10)\n    print(f\"✅ arXiv accessible → {r.status_code}\")\n    print(r.text[:150])\nexcept Exception as e:\n    print(f\"❌ {str(e)[:80]}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:59.449823Z","iopub.execute_input":"2026-06-15T19:34:59.450218Z","iopub.status.idle":"2026-06-15T19:34:59.504607Z","shell.execute_reply.started":"2026-06-15T19:34:59.450178Z","shell.execute_reply":"2026-06-15T19:34:59.503659Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import xml.etree.ElementTree as ET\n\ndef extract_geo_keywords(interpretation_text):\n    \"\"\"\n    Extracts geophysical keywords from Vision LLM output text.\n    Keywords are tuned to retrieve FWI and deep-learning seismic papers.\n    \"\"\"\n    keywords = []\n    text_lower = interpretation_text.lower()\n\n    if 'fault' in text_lower:\n        keywords.append('fault detection full waveform inversion')\n    if 'layer' in text_lower:\n        keywords.append('seismic velocity model building synthetic')\n    if 'velocity' in text_lower:\n        keywords.append('FWI velocity estimation neural network')\n\n    keywords.append('OpenFWI seismic benchmark deep learning')\n    return keywords[:3]\n\n\ndef search_arxiv(keywords, max_results=4):\n    \"\"\"Searches arXiv Atom API. Returns list of dicts: title, authors, summary, link.\"\"\"\n    query = ' AND '.join(f'all:{kw}' for kw in keywords)\n    params = {'search_query': query, 'max_results': max_results,\n              'sortBy': 'relevance', 'sortOrder': 'descending'}\n\n    r = requests.get(\"https://export.arxiv.org/api/query\", params=params, timeout=15)\n    if r.status_code != 200:\n        return []\n\n    ns   = {'atom': 'http://www.w3.org/2005/Atom'}\n    root = ET.fromstring(r.text)\n    papers = []\n\n    for entry in root.findall('atom:entry', ns):\n        papers.append({\n            'title':   entry.find('atom:title',   ns).text.strip().replace('\\n', ' '),\n            'summary': entry.find('atom:summary', ns).text.strip().replace('\\n', ' ')[:400],\n            'link':    entry.find('atom:id',      ns).text.strip(),\n            'authors': [a.find('atom:name', ns).text\n                        for a in entry.findall('atom:author', ns)][:3]\n        })\n    return papers\n\n\ndef rag_enrich_interpretation(interpretation, openai_key, verbose=True):\n    \"\"\"\n    Full RAG pipeline:\n    1. Extract keywords from Vision LLM output\n    2. Search arXiv for relevant papers\n    3. Build context and call GPT-4o mini (text only)\n    4. Return geological report with citations\n    \n    Cost: ~$0.0003 per call (text only, no image).\n    \"\"\"\n    keywords = extract_geo_keywords(interpretation)\n    if verbose:\n        print(f\"🔍 Keywords: {keywords}\")\n\n    papers = search_arxiv(keywords, max_results=4)\n    if verbose:\n        print(f\"📄 Papers found: {len(papers)}\")\n        for p in papers:\n            print(f\"   • {p['title'][:70]}...\")\n\n    if not papers:\n        return {\"error\": \"No papers found\", \"interpretation\": interpretation}\n\n    context = \"\\n\\n\".join([\n        f\"[{i+1}] {p['title']}\\nAuthors: {', '.join(p['authors'])}\\nAbstract: {p['summary']}\"\n        for i, p in enumerate(papers)\n    ])\n\n    rag_prompt = f\"\"\"You are an expert geophysicist writing a scientific interpretation report.\n\nSEISMIC IMAGE INTERPRETATION (from Vision AI):\n{interpretation}\n\nRELEVANT LITERATURE (from arXiv):\n{context}\n\nWrite a concise geological report (200-250 words) that:\n1. Confirms or refines the structural interpretation with scientific context\n2. Cites at least 2 papers above as [1], [2], etc.\n3. Suggests the best FWI approach for this subsurface\n4. Identifies the most likely geological scenario\n\nIMPORTANT: Use only 'depth index' and 'lateral index' — no metric units.\nThis is a synthetic dataset (OpenFWI) with pixel coordinates, not real meters.\"\"\"\n\n    r = requests.post(\n        \"https://api.openai.com/v1/chat/completions\",\n        headers={\"Authorization\": f\"Bearer {openai_key}\", \"Content-Type\": \"application/json\"},\n        json={\"model\": \"gpt-4o-mini\", \"max_tokens\": 400,\n              \"messages\": [{\"role\": \"user\", \"content\": rag_prompt}]},\n        timeout=60\n    )\n\n    resp  = r.json()\n    usage = resp['usage']\n    cost  = (usage['prompt_tokens'] * 0.00015 + usage['completion_tokens'] * 0.0006) / 1000\n\n    if verbose:\n        print(f\"\\n💰 RAG cost: ${cost:.5f} | Tokens: {usage['total_tokens']}\")\n\n    return {\n        \"report\":   resp['choices'][0]['message']['content'],\n        \"papers\":   papers,\n        \"keywords\": keywords,\n        \"cost_usd\": round(cost, 5)\n    }\n\nprint(\"✅ RAG pipeline ready\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:59.506011Z","iopub.execute_input":"2026-06-15T19:34:59.506413Z","iopub.status.idle":"2026-06-15T19:34:59.528058Z","shell.execute_reply.started":"2026-06-15T19:34:59.506383Z","shell.execute_reply":"2026-06-15T19:34:59.526419Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 7. Full Pipeline — FlatFault_A (3 samples)\n\nWe run the complete SOM → Vision LLM → RAG pipeline on three FlatFault_A samples:\n\n| Sample | Expected | SOM fault score |\n|--------|----------|-----------------|\n| idx=0  | No fault | low score → NO fault context |\n| idx=100 | Central fault | high score → fault context |\n| idx=400 | Complex fault | medium score |\n\n**Result:** idx=0 correctly reported \"absence of significant faulting\" after adding SOM context.\nWithout context, all three samples reported faults — demonstrating the value of the SOM bridge layer.","metadata":{}},{"cell_type":"code","source":"samples = [\n    (0,   'FlatFault_A — no fault'),\n    (100, 'FlatFault_A — central fault'),\n    (400, 'FlatFault_A — complex fault'),\n]\n\nall_results = {}\n\nfor idx, label in samples:\n    print(f\"\\n{'='*60}\")\n    print(f\"📍 {label} (idx={idx})\")\n    print('='*60)\n\n    # Layer 1 — SOM lithology map\n    litho_map, _, _ = compute_litho_map(idx, data_f, model_f, som_f, scaler_f)\n    fig = make_seismic_plot(idx, data_f, model_f, litho_map)\n\n    # Layer 2 — Vision LLM with SOM context\n    print(\"🔭 Vision LLM...\")\n    vision = analyze_seismic_vision(fig, OPENAI_KEY, litho_map=litho_map)\n    print(f\"   ${vision['cost_usd']} | {vision['prompt_tokens']} tokens\")\n\n    # Layer 3 — RAG\n    print(\"📚 RAG arXiv...\")\n    rag = rag_enrich_interpretation(vision['text'], OPENAI_KEY, verbose=False)\n    print(f\"   ${rag['cost_usd']} | {len(rag['papers'])} papers\")\n\n    total = vision['cost_usd'] + rag['cost_usd']\n    all_results[idx] = {'label': label, 'vision': vision['text'],\n                        'report': rag['report'], 'papers': rag['papers'], 'cost': total}\n\n    print(f\"\\n{rag['report']}\")\n    print(\"\\n📎 References:\")\n    for i, p in enumerate(rag['papers']):\n        print(f\"   [{i+1}] {p['title'][:65]}...\")\n    print(f\"\\n💰 Sample cost: ${total:.5f}\")\n\nprint(f\"\\n{'='*60}\")\ntotal_fault = sum(r['cost'] for r in all_results.values())\nprint(f\"FlatFault_A TOTAL: ${total_fault:.5f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:34:59.529705Z","iopub.execute_input":"2026-06-15T19:34:59.530143Z","iopub.status.idle":"2026-06-15T19:35:26.188358Z","shell.execute_reply.started":"2026-06-15T19:34:59.5301Z","shell.execute_reply":"2026-06-15T19:35:26.187131Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 8. FlatVel_A — Adapting the Pipeline\n\nFlatVel_A has no faults. We adapt the pipeline with:\n1. `compute_litho_map_vel()` — uses 1D SOM (`som2`/`scaler2`), features per depth row\n2. `detect_gradient_from_som()` — replaces fault detection with velocity gradient analysis\n3. `analyze_seismic_vision_vel()` — prompt adapted for layered compaction scenario\n\n### Variable naming note\n> ⚠️ During the notebook session, `model` was accidentally overwritten as a string \n> by a cell that used `model` as a generic variable name. \n> We reload FlatVel_A data as `data_v` / `model_v` to avoid collisions.","metadata":{}},{"cell_type":"code","source":"# Reload FlatVel_A with unambiguous variable names\nbase = '/kaggle/input/competitions/waveform-inversion/train_samples/FlatVel_A'\ndata_v  = np.load(f'{base}/data/data2.npy')   # (500, 5, 1000, 70)\nmodel_v = np.load(f'{base}/model/model2.npy') # (500, 1, 70, 70)\n\nprint(f\"data_v  : {data_v.shape}  | {type(data_v)}\")\nprint(f\"model_v : {model_v.shape} | {type(model_v)}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:35:26.189637Z","iopub.execute_input":"2026-06-15T19:35:26.189916Z","iopub.status.idle":"2026-06-15T19:35:26.731322Z","shell.execute_reply.started":"2026-06-15T19:35:26.189888Z","shell.execute_reply":"2026-06-15T19:35:26.730113Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def compute_litho_map_vel(idx, data_arr, model_arr, som, scaler):\n    \"\"\"\n    Computes lithology map for FlatVel_A using the 1D row-wise SOM (som2/scaler2).\n    Returns: litho_map (70,70), facies_1d (70,), vmap (70,70).\n    \"\"\"\n    seis = data_arr[idx]\n    vmap = model_arr[idx, 0]\n    depth_to_time = np.linspace(0, 999, 70).astype(int)\n\n    feats = []\n    for d in range(70):\n        t = depth_to_time[d]\n        traces_at_t = seis[:, t, :]\n        feat = np.concatenate([\n            traces_at_t.mean(axis=1), traces_at_t.max(axis=1),\n            traces_at_t.std(axis=1), (traces_at_t**2).mean(axis=1),\n            np.abs(np.diff(np.sign(traces_at_t))).sum(axis=1).astype(float)\n        ])\n        feats.append(feat)\n\n    feats_scaled = scaler.transform(np.array(feats))\n    winners      = np.array([som.winner(f) for f in feats_scaled])\n    facies_1d    = winners[:, 0] * 4 + winners[:, 1]\n    facies_map   = np.tile(facies_1d.reshape(-1, 1), (1, 70))\n\n    facies_vel = {}\n    for d in range(70):\n        f = facies_1d[d]\n        if f not in facies_vel:\n            facies_vel[f] = []\n        facies_vel[f].append(vmap[d, :].mean())\n    facies_vel = {f: np.mean(v) for f, v in facies_vel.items()}\n\n    litho_map = np.zeros((70, 70), dtype=int)\n    for f, vel in facies_vel.items():\n        clase, _, _ = vel_to_litho(vel)\n        litho_map[facies_map == f] = clase\n\n    return litho_map, facies_1d, vmap\n\n\ndef detect_gradient_from_som(litho_map, vmap):\n    \"\"\"\n    For FlatVel_A: no faults. Analyzes vertical velocity gradient and layer count.\n    Returns: n_layers, gradient (m/s per depth index), vel_min, vel_max, vel_std, homogeneous.\n    \"\"\"\n    vel_profile = vmap.mean(axis=1)\n    gradient    = np.diff(vel_profile).mean()\n    litho_col   = litho_map[:, 35]\n    transitions = np.sum(np.diff(litho_col) != 0)\n\n    return {\n        'n_layers':    int(transitions + 1),\n        'gradient':    round(float(gradient), 2),\n        'vel_min':     round(float(vel_profile.min()), 1),\n        'vel_max':     round(float(vel_profile.max()), 1),\n        'vel_std':     round(float(vel_profile.std()), 1),\n        'homogeneous': vel_profile.std() < 200\n    }\n\n\ndef analyze_seismic_vision_vel(fig, openai_key, litho_map, vmap, max_width=1200):\n    \"\"\"\n    Vision LLM call adapted for FlatVel_A.\n    Uses gradient analysis context instead of fault detection.\n    \"\"\"\n    buf_raw = io.BytesIO()\n    fig.savefig(buf_raw, format='png', dpi=120, bbox_inches='tight')\n    buf_raw.seek(0)\n\n    img_pil = Image.open(buf_raw)\n    w, h = img_pil.size\n    if w > max_width:\n        img_pil = img_pil.resize((max_width, int(h * max_width / w)), Image.LANCZOS)\n\n    buf = io.BytesIO()\n    img_pil.convert('RGB').save(buf, format='JPEG', quality=82)\n    buf.seek(0)\n    img_b64 = base64.b64encode(buf.read()).decode('utf-8')\n    plt.close(fig)\n\n    grad_info = detect_gradient_from_som(litho_map, vmap)\n    som_context = (\n        f\"CONTEXT FROM SOM ANALYSIS: {grad_info['n_layers']} lithological layer(s) detected. \"\n        f\"Mean velocity gradient = {grad_info['gradient']:.1f} m/s per depth index \"\n        f\"({'increasing' if grad_info['gradient'] > 0 else 'decreasing'} with depth). \"\n        f\"Velocity range: {grad_info['vel_min']:.0f}–{grad_info['vel_max']:.0f} m/s. \"\n        f\"{'Nearly homogeneous sample.' if grad_info['homogeneous'] else 'Significant velocity variation.'}\"\n    )\n\n    prompt = f\"\"\"{som_context}\n\nYou are an expert geophysicist. Analyze this seismic image.\nTHREE panels: LEFT=seismic waveform, CENTER=velocity map (m/s), RIGHT=SOM lithology.\n\nIMPORTANT: This is OpenFWI FlatVel_A — flat horizontal layers, NO faults.\nFocus on sedimentary layering and compaction trends.\nAxes show pixel indices (0–70), NOT metric units.\n\n1. LAYERS: How many distinct horizontal layers in CENTER? Depth index ranges.\n2. VELOCITY TREND: Increases, decreases, or constant with depth?\n3. LITHOLOGY: What rock types in RIGHT panel? Top vs bottom.\n4. SEISMIC PATTERN: Reflection patterns in LEFT panel.\n5. GEOLOGICAL SCENARIO: What depositional environment fits this velocity profile?\"\"\"\n\n    r = requests.post(\n        \"https://api.openai.com/v1/chat/completions\",\n        headers={\"Authorization\": f\"Bearer {openai_key}\", \"Content-Type\": \"application/json\"},\n        json={\"model\": \"gpt-4o-mini\", \"max_tokens\": 500,\n              \"messages\": [{\"role\": \"user\", \"content\": [\n                  {\"type\": \"image_url\",\n                   \"image_url\": {\"url\": f\"data:image/jpeg;base64,{img_b64}\", \"detail\": \"high\"}},\n                  {\"type\": \"text\", \"text\": prompt}\n              ]}]},\n        timeout=60\n    )\n\n    if r.status_code != 200:\n        return {\"error\": r.text[:200]}\n\n    resp  = r.json()\n    usage = resp['usage']\n    cost  = (usage['prompt_tokens'] * 0.00015 + usage['completion_tokens'] * 0.0006) / 1000\n\n    return {\n        \"text\":          resp['choices'][0]['message']['content'],\n        \"prompt_tokens\": usage['prompt_tokens'],\n        \"cost_usd\":      round(cost, 5),\n        \"grad_info\":     grad_info\n    }\n\nprint(\"✅ FlatVel_A pipeline functions ready\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:35:26.733418Z","iopub.execute_input":"2026-06-15T19:35:26.733796Z","iopub.status.idle":"2026-06-15T19:35:26.792318Z","shell.execute_reply.started":"2026-06-15T19:35:26.733763Z","shell.execute_reply":"2026-06-15T19:35:26.790544Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 9. Full Pipeline — FlatVel_A (3 samples)\n\nSelected by velocity variability extremes (computed in Section 2.3):\n- **idx=194**: high variability → many distinct layers, active compaction\n- **idx=69**: median → typical layered sedimentary section  \n- **idx=347**: low variability → nearly homogeneous, shallow burial\n\n**Validation criteria:** Reports should describe horizontal layering and velocity increase \nwith depth — and should NOT mention faults.","metadata":{}},{"cell_type":"code","source":"samples_vel = [\n    (194, 'FlatVel_A — high variability (std≈1219 m/s)'),\n    (69,  'FlatVel_A — median variability (std≈648 m/s)'),\n    (347, 'FlatVel_A — low variability (std≈92 m/s)'),\n]\n\nall_results_vel = {}\n\nfor idx, label in samples_vel:\n    print(f\"\\n{'='*60}\")\n    print(f\"📍 {label} (idx={idx})\")\n    print('='*60)\n\n    # Layer 1 — SOM\n    litho_map, facies_1d, vmap = compute_litho_map_vel(idx, data_v, model_v, som2, scaler2)\n    fig = make_seismic_plot(idx, data_v, model_v, litho_map)\n\n    # Layer 2 — Vision LLM (gradient context)\n    print(\"🔭 Vision LLM...\")\n    vision = analyze_seismic_vision_vel(fig, OPENAI_KEY, litho_map, vmap)\n    grad   = vision['grad_info']\n    print(f\"   ${vision['cost_usd']} | {vision['prompt_tokens']} tokens\")\n    print(f\"   SOM: {grad['n_layers']} layers | gradient={grad['gradient']} m/s/idx | \"\n          f\"vel {grad['vel_min']:.0f}–{grad['vel_max']:.0f} m/s\")\n\n    # Layer 3 — RAG\n    print(\"📚 RAG arXiv...\")\n    rag = rag_enrich_interpretation(vision['text'], OPENAI_KEY, verbose=False)\n    print(f\"   ${rag['cost_usd']} | {len(rag['papers'])} papers\")\n\n    total = vision['cost_usd'] + rag['cost_usd']\n    all_results_vel[idx] = {\n        'label': label, 'vision': vision['text'], 'report': rag['report'],\n        'papers': rag['papers'], 'grad_info': grad, 'cost': total\n    }\n\n    print(f\"\\n{rag['report']}\")\n    print(\"\\n📎 References:\")\n    for i, p in enumerate(rag['papers']):\n        print(f\"   [{i+1}] {p['title'][:65]}...\")\n    print(f\"\\n💰 Sample cost: ${total:.5f}\")\n\nprint(f\"\\n{'='*60}\")\ntotal_vel   = sum(r['cost'] for r in all_results_vel.values())\ntotal_fault = sum(r['cost'] for r in all_results.values())\ngrand_total = total_vel + total_fault\nprint(f\"FlatVel_A TOTAL  : ${total_vel:.5f}\")\nprint(f\"FlatFault_A TOTAL: ${total_fault:.5f}\")\nprint(f\"GRAND TOTAL      : ${grand_total:.5f}\")\nprint(f\"Estimated remaining (from $5.00): ~${5.00 - grand_total:.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-15T19:35:26.79399Z","iopub.execute_input":"2026-06-15T19:35:26.794426Z","iopub.status.idle":"2026-06-15T19:35:59.793292Z","shell.execute_reply.started":"2026-06-15T19:35:26.794385Z","shell.execute_reply":"2026-06-15T19:35:59.7921Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 10. Summary and Future Work\n\n### Results\n\n| Dataset | Samples | SOM QE | LLM fault detection | Total cost |\n|---------|---------|--------|---------------------|------------|\n| FlatVel_A | 3 | 0.88 | N/A (no faults) | ~$0.011 |\n| FlatFault_A | 3 | 0.44 | ✅ Correct (with SOM context) | ~$0.010 |\n| **Total** | **6** | — | — | **~$0.021** |\n\n### Key findings\n\n1. **SOM bridge layer is critical** — injecting SOM fault scores as context reduced false positives in GPT-4o mini from 3/3 to 1/3 for no-fault samples.\n2. **Image optimization matters** — resizing to 1200px JPEG cut token usage 2.4× vs. raw PNG with no quality loss for geological interpretation.\n3. **arXiv RAG provides scientifically grounded context** — papers on FWI with CNNs and multi-scale inversion are directly relevant and consistently retrieved.\n4. **Cost is negligible** — ~$0.003/sample makes large-scale batch analysis viable.\n\n### Known limitations\n\n- GPT-4o mini sometimes misidentifies fault location (lateral index can be off by ±10–15 indices)\n- The SOM fault threshold (0.3) is empirical — needs calibration on more samples\n- `VEL_RANGES` label names leaked into RAG reports in Spanish (now fixed to English)\n\n### Future work\n\n- **CurveVel_A**: curved velocity layers require the 2D pixel-wise SOM (like FlatFault_A), not the 1D row-wise approach used for FlatVel_A. The pipeline structure is ready — only `compute_litho_map` and the prompt context need adaptation.\n- **Layer 3 v2**: replace arXiv live search with FAISS local index over pre-downloaded FWI papers for faster, more controlled retrieval.\n- **Quantitative evaluation**: compare SOM lithology maps against ground truth velocity bins using IoU or accuracy metrics.\n- **GitHub + LinkedIn writeup**: share as open-source pipeline with reproducible Kaggle notebook.","metadata":{}}]}