{"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":"I muted the direct wave components since their high amplitude would obscure the relatively low-amplitude reflected waves. This should make the previously hidden reflections below more visible. Additionally, we can estimate the sound velocity of the upper layer using geometric calculations based on the direct waves.","metadata":{}},{"cell_type":"code","source":"import os\nimport sys\nimport json\nimport glob\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport torch\nimport zipfile\nimport io","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-09T23:29:11.615157Z","iopub.execute_input":"2025-04-09T23:29:11.615498Z","iopub.status.idle":"2025-04-09T23:29:11.621368Z","shell.execute_reply.started":"2025-04-09T23:29:11.615475Z","shell.execute_reply":"2025-04-09T23:29:11.620118Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"TEST_DATA_DIR = '/kaggle/input/waveform-inversion/test'\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-09T23:29:11.623824Z","iopub.execute_input":"2025-04-09T23:29:11.624189Z","iopub.status.idle":"2025-04-09T23:29:11.644924Z","shell.execute_reply.started":"2025-04-09T23:29:11.624150Z","shell.execute_reply":"2025-04-09T23:29:11.643765Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"DEBUG = True\nGAIN = 1.0 \n\nENABLE_FIRST_BREAK_MUTE = True\nDT = 0.001\n\nRECEIVER_INTERVAL = 10 \nDISTANCE_DECAY_POWER = 1.5 \nDISPLAY_SOURCE_INDEX = 1 \nMUTE_OFFSET_TIME = 0.08 \nPEAK_SEARCH_REC0_START = 0.07 \nPEAK_SEARCH_REC0_END = 0.09 \nPEAK_SEARCH_REC35_START = 0.20\nPEAK_SEARCH_REC35_END = 0.40 \nZIP_FILENAME = \"seismic_images_fb_mute.zip\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-09T23:29:11.646303Z","iopub.execute_input":"2025-04-09T23:29:11.646667Z","iopub.status.idle":"2025-04-09T23:29:11.669162Z","shell.execute_reply.started":"2025-04-09T23:29:11.646638Z","shell.execute_reply":"2025-04-09T23:29:11.667987Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"npy_files = sorted(glob.glob(os.path.join(TEST_DATA_DIR, '*.npy')))\ntotal_files = len(npy_files)\n\nprint(f\"Found {total_files} NPY files.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-09T23:29:11.671064Z","iopub.execute_input":"2025-04-09T23:29:11.671384Z","iopub.status.idle":"2025-04-09T23:29:11.847794Z","shell.execute_reply.started":"2025-04-09T23:29:11.671356Z","shell.execute_reply":"2025-04-09T23:29:11.846783Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if DEBUG:\n    npy_files = npy_files[:10]\n    print(f\"DEBUG mode enabled: Processing only the first {len(npy_files)} files.\")\n\n\nprint(f\"Processing files and creating {ZIP_FILENAME}...\")\n\nprocessed_count = 0\nerror_count = 0\n\n\ntry:\n    with zipfile.ZipFile(ZIP_FILENAME, 'w', zipfile.ZIP_DEFLATED) as zip_file:\n        for i, npy_path in enumerate(npy_files):\n            fig = None \n            velocity_est = None \n            mute_error_occurred = False \n\n            try:\n                file_name = os.path.basename(npy_path)\n                file_id = os.path.splitext(file_name)[0]\n\n                seismic_data_original = np.load(npy_path)\n                seismic_data = seismic_data_original.copy() \n\n                num_sources, time_steps, num_receivers = seismic_data.shape\n                print(f\"Processing {file_id}: shape=({num_sources}, {time_steps}, {num_receivers})\")\n\n                if ENABLE_FIRST_BREAK_MUTE:\n                    print(f\"  Applying FirstBreak Mute for {file_id}...\")\n                    try:\n                        if num_sources == 0:\n                             raise ValueError(\"FirstBreak Mute requires at least one source.\")\n\n                        data_src0 = seismic_data[0] # shape: (time_steps, num_receivers)\n\n                        search_start_idx_rec0 = int(PEAK_SEARCH_REC0_START / DT)\n                        search_end_idx_rec0 = int(PEAK_SEARCH_REC0_END / DT)\n                        if search_end_idx_rec0 > time_steps: search_end_idx_rec0 = time_steps\n                        if search_start_idx_rec0 >= search_end_idx_rec0:\n                            raise ValueError(f\"Invalid time range for receiver 0 peak detection ({PEAK_SEARCH_REC0_START}-{PEAK_SEARCH_REC0_END}s).\")\n                        if 0 >= num_receivers:\n                             raise ValueError(\"Receiver index 0 is out of bounds.\")\n\n                        peak_idx_rec0_rel = np.argmax(np.abs(data_src0[search_start_idx_rec0:search_end_idx_rec0, 0]))\n                        peak_idx_rec0 = peak_idx_rec0_rel + search_start_idx_rec0\n                        peak_time_rec0 = peak_idx_rec0 * DT\n                        print(f\"    Peak time at receiver 0 (source 0): {peak_time_rec0:.4f} s (index {peak_idx_rec0})\")\n\n\n                        target_receiver_idx = 35\n                        if target_receiver_idx >= num_receivers:\n                            raise ValueError(f\"Receiver index {target_receiver_idx} is out of bounds (num_receivers={num_receivers}). Cannot estimate velocity.\")\n\n                        search_start_idx_rec35 = int(PEAK_SEARCH_REC35_START / DT)\n                        search_end_idx_rec35 = int(PEAK_SEARCH_REC35_END / DT)\n                        if search_end_idx_rec35 > time_steps: search_end_idx_rec35 = time_steps\n                        if search_start_idx_rec35 >= search_end_idx_rec35:\n                            raise ValueError(f\"Invalid time range for receiver {target_receiver_idx} peak detection ({PEAK_SEARCH_REC35_START}-{PEAK_SEARCH_REC35_END}s).\")\n\n                        peak_idx_rec35_rel = np.argmax(np.abs(data_src0[search_start_idx_rec35:search_end_idx_rec35, target_receiver_idx]))\n                        peak_idx_rec35 = peak_idx_rec35_rel + search_start_idx_rec35\n                        peak_time_rec35 = peak_idx_rec35 * DT\n                        print(f\"    Peak time at receiver {target_receiver_idx} (source 0): {peak_time_rec35:.4f} s (index {peak_idx_rec35})\")\n\n\n                        delta_t = peak_time_rec35 - peak_time_rec0\n                        delta_x = target_receiver_idx * RECEIVER_INTERVAL\n                        if delta_t <= 1e-6: \n                            print(f\"    Warning: Calculated delta_t is too small ({delta_t:.4f}s). Skipping velocity estimation and mute.\")\n                            velocity_est = None\n                            mute_error_occurred = True\n                        else:\n                            velocity_est = delta_x / delta_t\n                            print(f\"    Estimated velocity: {velocity_est:.2f} m/s\")\n\n                            print(f\"    Applying mute...\")\n                            \n                            if num_sources > 1:\n                                source_receiver_scale = (num_receivers - 1) / (num_sources - 1)\n                            else:\n                                source_receiver_scale = 0 \n\n                            for s in range(num_sources):\n                                source_location_approx = s * source_receiver_scale \n                                for r in range(num_receivers):\n                                    receiver_location = r \n                                    horizontal_receiver_index_diff = abs(receiver_location - source_location_approx)\n                                    horizontal_distance = horizontal_receiver_index_diff * RECEIVER_INTERVAL\n\n                                    estimated_arrival_time = peak_time_rec0 + horizontal_distance / velocity_est\n\n                                    mute_end_time = estimated_arrival_time + MUTE_OFFSET_TIME\n                                    mute_end_idx = min(int(mute_end_time / DT), time_steps)\n\n                                    if mute_end_idx > 0:\n                                        seismic_data[s, :mute_end_idx, r] = 0.0\n                            print(f\"    Mute applied.\")\n\n\n                    except Exception as mute_error:\n                        print(f\"  Error during FirstBreak Mute for file {file_id}: {mute_error}\")\n                        velocity_est = None \n                        mute_error_occurred = True \n                        seismic_data = seismic_data_original.copy() \n                        print(f\"    Reverted to original data for {file_id} due to mute error.\")\n                else:\n                    print(f\"  FirstBreak Mute disabled for {file_id}.\")\n\n\n\n                if DISPLAY_SOURCE_INDEX < 0 or DISPLAY_SOURCE_INDEX >= num_sources:\n                    print(f\"  Warning: DISPLAY_SOURCE_INDEX ({DISPLAY_SOURCE_INDEX}) is out of range [0, {num_sources-1}]. Using central source index {num_sources // 2} instead.\")\n                    display_source_idx = num_sources // 2\n                else:\n                    display_source_idx = DISPLAY_SOURCE_INDEX\n\n                original_display_data = seismic_data_original[display_source_idx]\n                muted_display_data = seismic_data[display_source_idx] # ミュートされていなければオリジナルと同じ\n\n                # --- AGC ---\n                print(f\"  Applying AGC for {file_id} (source {display_source_idx})...\")\n                epsilon = 1e-8\n                times = np.arange(time_steps) * DT\n                gain_time = (times + epsilon)**DISTANCE_DECAY_POWER\n\n                # 1. 時間減衰補正 AGC\n                agc_data_time_original = original_display_data * gain_time[:, np.newaxis]\n                rms_original = np.sqrt(np.mean(agc_data_time_original**2, axis=1))\n                gain_rms_original = GAIN / (rms_original + epsilon)\n                agc_data_original = agc_data_time_original * gain_rms_original[:, np.newaxis]\n\n                # 2. RMSベース AGC \n                agc_data_time_muted = muted_display_data * gain_time[:, np.newaxis]\n                rms_muted = np.sqrt(np.mean(agc_data_time_muted**2, axis=1))\n                gain_rms_muted = GAIN / (rms_muted + epsilon)\n                agc_data_muted = agc_data_time_muted * gain_rms_muted[:, np.newaxis]\n\n                print(f\"  AGC applied.\")\n                # --- AGC  ---\n\n\n                fig, axes = plt.subplots(1, 2, figsize=(9.6, 4.8), sharey=True) \n\n                im_orig = axes[0].imshow(agc_data_original, aspect='auto', cmap='gray')\n                axes[0].set_title(\"Original (AGC Only)\")\n                axes[0].set_xlabel(\"Receivers\")\n                axes[0].set_ylabel(f\"Timesteps (dt={DT}s)\")\n\n                im_muted = axes[1].imshow(agc_data_muted, aspect='auto', cmap='gray')\n                mute_status_for_title = \"Muted\"\n                if not ENABLE_FIRST_BREAK_MUTE:\n                    mute_status_for_title = \"Mute Disabled\"\n                elif mute_error_occurred:\n                    mute_status_for_title = \"Mute Failed\"\n                axes[1].set_title(f\"{mute_status_for_title} (AGC Applied)\")\n                axes[1].set_xlabel(\"Receivers\")\n\n                title_lines = [\n                    f\"Seismic Data Comparison - ID {file_id}, Source Index {display_source_idx}\",\n                    f\"AGC: Time (n={DISTANCE_DECAY_POWER}), RMS (G={GAIN})\"\n                ]\n                if ENABLE_FIRST_BREAK_MUTE:\n                    if mute_error_occurred:\n                        mute_status = \"Failed\"\n                        title_lines.append(f\"FB Mute Status: {mute_status}\")\n                    elif velocity_est is not None:\n                        mute_status = \"Enabled\"\n                        title_lines.append(f\"FB Mute Status: {mute_status} (Est. Vel: {velocity_est:.2f} m/s)\")\n                    else:\n                        mute_status = \"Enabled (Vel. Est. Skipped)\"\n                        title_lines.append(f\"FB Mute Status: {mute_status}\")\n                else:\n                     mute_status = \"Disabled\"\n                     title_lines.append(f\"FB Mute Status: {mute_status}\")\n\n                fig.suptitle(\"\\n\".join(title_lines), y=1.02) \n\n                cbar_orig = fig.colorbar(im_orig, ax=axes[0], fraction=0.046, pad=0.04)\n                cbar_orig.set_label(\"Amplitude (AGC Only)\")\n                cbar_muted = fig.colorbar(im_muted, ax=axes[1], fraction=0.046, pad=0.04)\n                cbar_muted.set_label(\"Amplitude (Muted + AGC)\")\n\n                plt.tight_layout(rect=[0, 0.03, 1, 0.98]) \n\n                buffer = io.BytesIO()\n                plt.savefig(buffer, format='png', dpi=100)\n                buffer.seek(0)\n\n                zip_file.writestr(f\"{file_id}.png\", buffer.getvalue())\n                buffer.close() \n\n                if i < 10:\n                    print(f\"Displaying plot for {file_id}...\")\n                    plt.show() \n\n                plt.close(fig) \n                processed_count += 1\n                print(f\"  Successfully processed and added {file_id}.png to ZIP.\")\n\n\n            except Exception as e:\n                print(f\"Error processing file {npy_path}: {e}\")\n                error_count += 1\n                if fig is not None and plt.fignum_exists(fig.number):\n                    plt.close(fig)\n\nexcept Exception as e:\n    print(f\"An error occurred during ZIP file creation or processing loop: {e}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-09T23:29:11.848970Z","iopub.execute_input":"2025-04-09T23:29:11.849312Z","iopub.status.idle":"2025-04-09T23:29:22.746974Z","shell.execute_reply.started":"2025-04-09T23:29:11.849281Z","shell.execute_reply":"2025-04-09T23:29:22.745817Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 最終的な結果を表示\nprint(f\"\\nProcessing finished.\")\nif DEBUG:\n    print(\"(Ran in DEBUG mode)\")\nprint(f\"Successfully processed and added to ZIP: {processed_count} files.\")\nprint(f\"Failed to process: {error_count} files.\")\nif processed_count > 0:\n    print(f\"Output saved to {ZIP_FILENAME}\")\nelse:\n    print(f\"{ZIP_FILENAME} was not created as no files were processed successfully.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-09T23:29:22.748432Z","iopub.execute_input":"2025-04-09T23:29:22.749282Z","iopub.status.idle":"2025-04-09T23:29:22.755360Z","shell.execute_reply.started":"2025-04-09T23:29:22.749255Z","shell.execute_reply":"2025-04-09T23:29:22.754167Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}