{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":39763,"databundleVersionId":11756775,"sourceType":"competition"}],"dockerImageVersionId":31012,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 🌋 Seismic Waveform Inversion - EDA \n\nThis notebook explores a dataset used in **seismic waveform inversion**, where the goal is to predict subsurface velocity structures from time-series waveforms recorded by sensors.\n\nWe'll:\n- Load and visualize the data\n- Plot multichannel time series and 2D sensor heatmaps\n- Extract statistical and frequency-domain features\n\nUseful for geophysics, signal processing, and inverse modeling tasks.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"# 📦 Imports","metadata":{}},{"cell_type":"code","source":"import torch\nimport numpy as np\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\nfrom scipy.stats import skew, kurtosis\nfrom scipy.signal import welch\nfrom torch.utils.data import Dataset","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-14T07:41:49.564656Z","iopub.execute_input":"2025-04-14T07:41:49.565071Z","iopub.status.idle":"2025-04-14T07:41:55.422323Z","shell.execute_reply.started":"2025-04-14T07:41:49.565038Z","shell.execute_reply":"2025-04-14T07:41:55.421356Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🔧 Configuration and Path Management\n","metadata":{}},{"cell_type":"code","source":"class Config:\n    BASE_PATH = '/kaggle/input/waveform-inversion/'\n    N_SAMPLES_PER_FILE = 500\n    SENSOR_GRID = (7, 10) \n\nclass PATHS:\n    train = f'{Config.BASE_PATH}/train_samples'\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-14T07:42:12.356871Z","iopub.execute_input":"2025-04-14T07:42:12.357382Z","iopub.status.idle":"2025-04-14T07:42:12.362648Z","shell.execute_reply.started":"2025-04-14T07:42:12.357356Z","shell.execute_reply":"2025-04-14T07:42:12.361515Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Dataset Manager","metadata":{"execution":{"iopub.status.busy":"2025-04-14T07:43:30.843693Z","iopub.execute_input":"2025-04-14T07:43:30.844786Z","iopub.status.idle":"2025-04-14T07:43:30.849464Z","shell.execute_reply.started":"2025-04-14T07:43:30.844753Z","shell.execute_reply":"2025-04-14T07:43:30.847864Z"}}},{"cell_type":"code","source":"class FileManager:\n    @staticmethod\n    def input_to_output(files):\n        return [\n            Path(str(f).replace('seis','vel').replace('data','model')) for f in files\n        ]\n\n    @staticmethod\n    def get_file_lists():\n        all_inputs = [f for f in Path(PATHS.train).rglob('*.npy') if ('seis' in f.stem) or ('data' in f.stem)]\n        train_inputs = [all_inputs[i] for i in range(0, len(all_inputs), 2)]\n        valid_inputs = [f for f in all_inputs if f not in train_inputs]\n\n        train_outputs = FileManager.input_to_output(train_inputs)\n        valid_outputs = FileManager.input_to_output(valid_inputs)\n\n        return train_inputs, train_outputs, valid_inputs, valid_outputs\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-14T07:44:03.420116Z","iopub.execute_input":"2025-04-14T07:44:03.420464Z","iopub.status.idle":"2025-04-14T07:44:03.427437Z","shell.execute_reply.started":"2025-04-14T07:44:03.420438Z","shell.execute_reply":"2025-04-14T07:44:03.426395Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 📦 Custom PyTorch Dataset\nEach waveform is a 3D array: `(channels, time steps, sensors)`  \nEach target is a 2D seismic velocity map.","metadata":{"execution":{"iopub.status.busy":"2025-04-14T07:44:13.839254Z","iopub.execute_input":"2025-04-14T07:44:13.839600Z","iopub.status.idle":"2025-04-14T07:44:13.845467Z","shell.execute_reply.started":"2025-04-14T07:44:13.839576Z","shell.execute_reply":"2025-04-14T07:44:13.844349Z"}}},{"cell_type":"code","source":"class SeismicDataset(Dataset):\n    def __init__(self, input_files, output_files, n_samples_per_file=500):\n        self.input_files = input_files\n        self.output_files = output_files\n        self.n_samples_per_file = n_samples_per_file\n\n    def __len__(self):\n        return len(self.input_files) * self.n_samples_per_file\n\n    def __getitem__(self, idx):\n        file_idx = idx // self.n_samples_per_file\n        sample_idx = idx % self.n_samples_per_file\n        \n        X = np.load(self.input_files[file_idx], mmap_mode='r')\n        y = np.load(self.output_files[file_idx], mmap_mode='r')\n\n        try:\n            return X[sample_idx].copy(), y[sample_idx].copy()\n        finally:\n            del X, y","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-14T07:44:29.381853Z","iopub.execute_input":"2025-04-14T07:44:29.383116Z","iopub.status.idle":"2025-04-14T07:44:29.389928Z","shell.execute_reply.started":"2025-04-14T07:44:29.383086Z","shell.execute_reply":"2025-04-14T07:44:29.388813Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 📥 Load One Sample","metadata":{"execution":{"iopub.status.busy":"2025-04-14T07:44:44.736202Z","iopub.execute_input":"2025-04-14T07:44:44.736595Z","iopub.status.idle":"2025-04-14T07:44:44.740988Z","shell.execute_reply.started":"2025-04-14T07:44:44.736568Z","shell.execute_reply":"2025-04-14T07:44:44.739985Z"}}},{"cell_type":"code","source":"train_inputs, train_outputs, valid_inputs, valid_outputs = FileManager.get_file_lists()\ndataset = SeismicDataset(train_inputs, train_outputs)\nx, y = dataset[1000]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-14T07:44:58.856192Z","iopub.execute_input":"2025-04-14T07:44:58.856520Z","iopub.status.idle":"2025-04-14T07:44:58.893425Z","shell.execute_reply.started":"2025-04-14T07:44:58.856496Z","shell.execute_reply":"2025-04-14T07:44:58.892689Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🗺️ Target Heatmap (Velocity Field)\nThis is the target velocity model. Each pixel indicates subsurface properties (e.g., velocity of wave propagation).","metadata":{"execution":{"iopub.status.busy":"2025-04-14T07:45:13.115160Z","iopub.execute_input":"2025-04-14T07:45:13.115478Z","iopub.status.idle":"2025-04-14T07:45:13.122170Z","shell.execute_reply.started":"2025-04-14T07:45:13.115458Z","shell.execute_reply":"2025-04-14T07:45:13.120837Z"}}},{"cell_type":"code","source":"plt.figure(figsize=(6, 6))\nplt.imshow(y[0], cmap='seismic', origin='lower')\nplt.colorbar(label='Amplitude')\nplt.title('Seismic Field (1 Channel, 70x70 Grid)')\nplt.xlabel('Receiver Index')\nplt.ylabel('Source Index')\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-14T07:45:25.908347Z","iopub.execute_input":"2025-04-14T07:45:25.908639Z","iopub.status.idle":"2025-04-14T07:45:26.343166Z","shell.execute_reply.started":"2025-04-14T07:45:25.908620Z","shell.execute_reply":"2025-04-14T07:45:26.342331Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 📉 Multichannel Time Series + Sensor Grid Heatmap\nWe'll plot time series for each channel and overlay the target heatmap for spatial context.\n","metadata":{"execution":{"iopub.status.busy":"2025-04-14T07:45:35.547554Z","iopub.execute_input":"2025-04-14T07:45:35.547871Z","iopub.status.idle":"2025-04-14T07:45:35.554192Z","shell.execute_reply.started":"2025-04-14T07:45:35.547848Z","shell.execute_reply":"2025-04-14T07:45:35.552835Z"}}},{"cell_type":"code","source":"heatmap = np.squeeze(y)\n\nfor ch in range(5):\n    fig, axs = plt.subplots(1, 2, figsize=(14, 5))\n\n    for sensor in range(70):\n        axs[0].plot(x[ch, :, sensor], alpha=0.4)\n    axs[0].set_title(f'All Sensors - Channel {ch}')\n    axs[0].set_xlabel('Time')\n    axs[0].set_ylabel('Amplitude')\n    axs[0].grid(True)\n\n    im = axs[1].imshow(heatmap, cmap='viridis', origin='lower')\n    axs[1].set_title('Sensor Grid Heatmap (Target)')\n    plt.colorbar(im, ax=axs[1], fraction=0.046, pad=0.04)\n\n    plt.suptitle(f'Channel {ch} - Time Series & Heatmap')\n    plt.tight_layout(rect=[0, 0.03, 1, 0.95])\n    plt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-14T07:45:46.801737Z","iopub.execute_input":"2025-04-14T07:45:46.802156Z","iopub.status.idle":"2025-04-14T07:45:50.090113Z","shell.execute_reply.started":"2025-04-14T07:45:46.802126Z","shell.execute_reply":"2025-04-14T07:45:50.088980Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🛰️ 2D Sensor Grid at a Single Time Slice\nEach subplot represents sensor amplitudes at a fixed time step.\n","metadata":{"execution":{"iopub.status.busy":"2025-04-14T07:46:05.700981Z","iopub.execute_input":"2025-04-14T07:46:05.701320Z","iopub.status.idle":"2025-04-14T07:46:05.714028Z","shell.execute_reply.started":"2025-04-14T07:46:05.701300Z","shell.execute_reply":"2025-04-14T07:46:05.712722Z"}}},{"cell_type":"code","source":"time_index = 800\nfor channel in range(5):\n    snapshot = x[channel, time_index, :].reshape(Config.SENSOR_GRID)\n    plt.figure(figsize=(6, 5))\n    plt.imshow(snapshot, cmap='seismic', origin='lower')\n    plt.colorbar(label='Amplitude')\n    plt.title(f'2D Sensor Grid at Time {time_index} - Channel {channel}')\n    plt.xlabel('Sensor X')\n    plt.ylabel('Sensor Y')\n    plt.tight_layout()\n    plt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-14T07:46:26.958528Z","iopub.execute_input":"2025-04-14T07:46:26.958844Z","iopub.status.idle":"2025-04-14T07:46:28.530179Z","shell.execute_reply.started":"2025-04-14T07:46:26.958820Z","shell.execute_reply":"2025-04-14T07:46:28.529300Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 📊 Statistical & Spectral Features Per Sensor\nWe'll extract descriptive statistics from each waveform including:\n- **Time domain**: mean, std, energy, RMS, skewness, kurtosis, zero-crossings\n- **Frequency domain**: spectral entropy, band power (0.5-3 Hz, 3-8 Hz, 8-20 Hz)\n","metadata":{}},{"cell_type":"code","source":"from pprint import pprint\nclass SensorFeatureExtractor:\n    @staticmethod\n    def compute_features(x, fs=100):\n        C, T, S = x.shape\n        stats = {}\n\n        for sensor in range(S):\n            stats[f\"Sensor_{sensor}\"] = {}\n            for ch in range(C):\n                sig = x[ch, :, sensor]\n\n                # Time domain\n                mean_val = np.mean(sig)\n                std_val = np.std(sig)\n                rms = np.sqrt(np.mean(sig ** 2))\n                energy = np.sum(sig ** 2)\n                peak = np.max(np.abs(sig))\n                zcr = ((sig[:-1] * sig[1:]) < 0).sum()\n                skewness = skew(sig)\n                kurt_val = kurtosis(sig)\n\n                # Frequency domain\n                freqs, psd = welch(sig, fs=fs)\n                psd_norm = psd / np.sum(psd)\n                spec_entropy = -np.sum(psd_norm * np.log2(psd_norm + 1e-8))\n\n                def band_power(fmin, fmax):\n                    idx = np.logical_and(freqs >= fmin, freqs <= fmax)\n                    return np.sum(psd[idx])\n\n                bp_0_5_3 = band_power(0.5, 3)\n                bp_3_8 = band_power(3, 8)\n                bp_8_20 = band_power(8, 20)\n\n                stats[f\"Sensor_{sensor}\"][f\"Channel_{ch}\"] = {\n                    \"mean\": mean_val,\n                    \"std\": std_val,\n                    \"rms\": rms,\n                    \"energy\": energy,\n                    \"peak_amplitude\": peak,\n                    \"zero_crossings\": zcr,\n                    \"skewness\": skewness,\n                    \"kurtosis\": kurt_val,\n                    \"spectral_entropy\": spec_entropy,\n                    \"bandpower_0.5-3Hz\": bp_0_5_3,\n                    \"bandpower_3-8Hz\": bp_3_8,\n                    \"bandpower_8-20Hz\": bp_8_20\n                }\n\n        return stats\n\nfeatures = SensorFeatureExtractor.compute_features(x)\npprint(features)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-14T07:47:45.280613Z","iopub.execute_input":"2025-04-14T07:47:45.280891Z","iopub.status.idle":"2025-04-14T07:47:46.042337Z","shell.execute_reply.started":"2025-04-14T07:47:45.280872Z","shell.execute_reply":"2025-04-14T07:47:46.041296Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## ✅ Summary\n\nWe’ve:\n- Loaded and sampled seismic waveform data\n- Visualized targets, multichannel waveforms, and 2D sensor heatmaps\n- Extracted statistical and frequency-based features for ML models\n","metadata":{"execution":{"iopub.status.busy":"2025-04-14T07:47:05.096154Z","iopub.execute_input":"2025-04-14T07:47:05.096460Z","iopub.status.idle":"2025-04-14T07:47:05.103127Z","shell.execute_reply.started":"2025-04-14T07:47:05.096439Z","shell.execute_reply":"2025-04-14T07:47:05.101951Z"}}}]}