{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","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":30919,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T06:16:50.643920Z","iopub.execute_input":"2025-04-17T06:16:50.644420Z","iopub.status.idle":"2025-04-17T06:16:51.084609Z","shell.execute_reply.started":"2025-04-17T06:16:50.644373Z","shell.execute_reply":"2025-04-17T06:16:51.083394Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"------\n## **MATLAB code → Python code**\nWe can find the contents of \"Semitic Forward Modeling Details\" in Appendix C of the paper.\n- discussion : https://www.kaggle.com/competitions/waveform-inversion/discussion/572329#3176072\n  \nThe code used by Forward Modeling can be found [here](https://csim.kaust.edu.sa/files/SeismicInversion/Chapter.FD/lab.FD2.8/lab.html), but it's a matlab code. So I used **chatGPT** to convert it into python code.\n- 1. Please convert the given MATLAB code (in units of functions) into python code.\n- 2. Does the given MATLAB code (in functional units) and python code act exactly the same?\n\nI have completed the python code below through these two steps.  \nThe input information of the main function (bel_to_seis function below) was also estimated based on the paper. However, there was some trial and error in the value of the coord, which is estimated to be the location of the receiver.","metadata":{"execution":{"iopub.status.busy":"2025-04-16T13:47:55.989917Z","iopub.execute_input":"2025-04-16T13:47:55.990275Z","iopub.status.idle":"2025-04-16T13:47:55.996695Z","shell.execute_reply.started":"2025-04-16T13:47:55.990250Z","shell.execute_reply":"2025-04-16T13:47:55.995305Z"}}},{"cell_type":"code","source":"# https://arxiv.org/pdf/2111.02926\n# https://csim.kaust.edu.sa/files/SeismicInversion/Chapter.FD/lab.FD2.8/lab.html\n\ndef ricker(f, dt, nt=None):\n    nw = int(2.2 / f / dt)\n    nw = 2 * (nw // 2) + 1\n    nc = nw // 2 + 1  # 중심 인덱스를 1-based 기준으로 설정\n\n    k = np.arange(1, nw + 1)  # 1-based index\n    alpha = (nc - k) * f * dt * np.pi\n    beta = alpha ** 2\n    w0 = (1.0 - 2.0 * beta) * np.exp(-beta)\n\n    # 1-based wavelet 생성\n    if nt is not None:\n        if nt < len(w0):\n            raise ValueError(\"nt is smaller than condition!\")\n        w = np.zeros(nt + 1)  # dummy 포함\n        w[1:len(w0) + 1] = w0\n    else:\n        w = np.zeros(len(w0) + 1)\n        w[1:] = w0\n\n    # 1-based time axis 생성\n    if nt is not None:\n        tw = np.arange(1, len(w)) * dt\n    else:\n        tw = np.arange(1, len(w)) * dt\n\n    return w, tw\n\ndef padvel(v0, nbc):\n    v_padded = np.pad(v0, ((nbc, nbc), (nbc, nbc)), mode='edge')\n    nz, nx = v_padded.shape\n    v = np.zeros((nz + 1, nx + 1))\n    v[1:, 1:] = v_padded\n    return v\n\t\ndef expand_source(s0, nt):\n    s0 = np.asarray(s0).flatten()\n    s = np.zeros(nt + 1)\n    s[1:len(s0) + 1] = s0\n    return s\n\t\ndef adjust_sr(coord, dx, nbc):\n    isx = int(round(coord['sx'] / dx)) + 1 + nbc\n    isz = int(round(coord['sz'] / dx)) + 1 + nbc\n    igx = (np.round(np.array(coord['gx']) / dx) + 1 + nbc).astype(int)\n    igz = (np.round(np.array(coord['gz']) / dx) + 1 + nbc).astype(int)\n\n    if abs(coord['sz']) < 0.5:\n        isz += 1\n    igz = igz + (np.abs(np.array(coord['gz'])) < 0.5).astype(int)\n    return isx, isz, igx, igz\n\t\ndef AbcCoef2D(vel, nbc, dx):\n    nzbc, nxbc = vel.shape[1] - 1, vel.shape[0] - 1  # 실제 사이즈\n    velmin = np.min(vel[1:, 1:])\n    nz = nzbc - 2 * nbc\n    nx = nxbc - 2 * nbc\n\n    a = (nbc - 1) * dx\n    kappa = 3.0 * velmin * np.log(1e7) / (2.0 * a)\n\n    damp1d = kappa * (((np.arange(1, nbc + 1) - 1) * dx / a) ** 2)\n    damp = np.zeros((nzbc + 1, nxbc + 1))\n\n    for iz in range(1, nzbc + 1):\n        damp[iz, 1:nbc + 1] = damp1d[::-1]\n        damp[iz, nx + nbc + 1 : nx + 2 * nbc + 1] = damp1d\n\n    for ix in range(nbc + 1, nbc + nx + 1):\n        damp[1:nbc + 1, ix] = damp1d[::-1]\n        damp[nz + nbc + 1 : nz + 2 * nbc + 1, ix] = damp1d\n\n    return damp\n\t\ndef a2d_mod_abc24(v, nbc, dx, nt, dt, s, coord, isFS):\n    ng = len(coord['gx'])\n    seis = np.zeros((nt + 1, ng))  # 1-based time axis\n\n    c1 = -2.5\n    c2 = 4.0 / 3.0\n    c3 = -1.0 / 12.0\n\n    v = padvel(v, nbc)\n    abc = AbcCoef2D(v, nbc, dx)\n\n    alpha = (v * dt / dx) ** 2\n    kappa = abc * dt\n    temp1 = 2 + 2 * c1 * alpha - kappa\n    temp2 = 1 - kappa\n    beta_dt = (v * dt) ** 2\n    s = expand_source(s, nt)\n    isx, isz, igx, igz = adjust_sr(coord, dx, nbc)\n\n    p0 = np.zeros_like(v)\n    p1 = np.zeros_like(v)\n\n    # Time Loop (1-based)\n    for it in range(1, nt + 1):\n        p = (temp1 * p1 - temp2 * p0 +\n             alpha * (\n                 c2 * (np.roll(p1, 1, axis=1) + np.roll(p1, -1, axis=1) +\n                       np.roll(p1, 1, axis=0) + np.roll(p1, -1, axis=0)) +\n                 c3 * (np.roll(p1, 2, axis=1) + np.roll(p1, -2, axis=1) +\n                       np.roll(p1, 2, axis=0) + np.roll(p1, -2, axis=0))\n             ))\n\n        # Source\n        p[isz, isx] += beta_dt[isz, isx] * s[it]\n\n        # Free Surface\n        if isFS:\n            p[nbc, :] = 0.0\n            p[nbc - 1 : nbc + 1, :] = -p[nbc + 1 : nbc + 3, :]\n\n        for ig in range(ng):\n            seis[it, ig] = p[igz[ig], igx[ig]]\n\n        p0, p1 = p1, p\n\n    return seis\n    \ndef vel_to_seis(vel):\n    \"\"\"\n    vel : (70, 70)\n    output : (5, 1001, 70)\n    \"\"\"\n    # 1. 모델 및 파라미터 설정\n    nz = 70\n    nx = 70\n    dx = 10\n    nbc = 120\n    nt = 1001\n    dt = 1e-3\n    freq = 15\n    isFS = False  # 자유표면 사용 여부\n    \n    # 2. Ricker 파형 생성\n    s, _ = ricker(freq, dt)\n    \n    # 3. 소스 및 수신기 설정\n    coord = {}\n    coord['sz'] = 1 * dx\n    coord['gx'] = np.arange(1, nx + 1) * dx\n    coord['gz'] = np.ones_like(coord['gx']) * dx\n    \n    # 4. 파동장 시뮬레이션 수행 : 소스 위치만 바꿔서\n    seis_data = []\n    for source_x in [1, 18, 35, 53, 70]:\n        coord['sx'] = source_x * dx        \n        \n        # 시뮬레이션\n        seis = a2d_mod_abc24(vel, nbc, dx, nt, dt, s, coord, isFS)\n\n        seis_data += [seis]\n        \n    return np.stack(seis_data, axis=0)\n    \n    \ndef plot_seis(seis):\n    \"\"\"\n    seis : (5, 1000, 70)\n    \"\"\"\n    fig,ax=plt.subplots(1,5,figsize=(20,5))\n    ax[0].imshow(seis[0, :, :],extent=[0,70,1000,0],aspect='auto',cmap='gray',vmin=-0.5,vmax=0.5)\n    ax[1].imshow(seis[1, :, :],extent=[0,70,1000,0],aspect='auto',cmap='gray',vmin=-0.5,vmax=0.5)\n    ax[2].imshow(seis[2, :, :],extent=[0,70,1000,0],aspect='auto',cmap='gray',vmin=-0.5,vmax=0.5)\n    ax[3].imshow(seis[3, :, :],extent=[0,70,1000,0],aspect='auto',cmap='gray',vmin=-0.5,vmax=0.5)\n    ax[4].imshow(seis[4, :, :],extent=[0,70,1000,0],aspect='auto',cmap='gray',vmin=-0.5,vmax=0.5)\n    for axis in ax:\n        axis.set_xticks(range(0, 70, 10))\n        axis.set_xticklabels(range(0, 700, 100))\n        axis.set_yticks(range(0, 2000, 1000))\n        axis.set_yticklabels(range(0, 2,1))\n        axis.set_ylabel('Time (s)', fontsize=12)\n        axis.set_xlabel('Offset (m)', fontsize=12)\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T06:16:52.432541Z","iopub.execute_input":"2025-04-17T06:16:52.433108Z","iopub.status.idle":"2025-04-17T06:16:52.460851Z","shell.execute_reply.started":"2025-04-17T06:16:52.433072Z","shell.execute_reply":"2025-04-17T06:16:52.459414Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"------\n## **Comparison**\nWhen I enter the velocity map provided by OpenFWI into the converted python code, I will check if Seismic data is generated well.  \nThe generated Seismic data and the Seismic data provided by OpenFWI should be almost identical.","metadata":{}},{"cell_type":"code","source":"seis_files = [\n    \"/kaggle/input/waveform-inversion/train_samples/FlatVel_A/data/data1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/FlatVel_B/data/data1.npy\",\n    \n    \"/kaggle/input/waveform-inversion/train_samples/CurveVel_A/data/data1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/CurveVel_B/data/data1.npy\",\n\n    \"/kaggle/input/waveform-inversion/train_samples/FlatFault_A/seis2_1_0.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/FlatFault_B/seis6_1_0.npy\",\n    \n    \"/kaggle/input/waveform-inversion/train_samples/CurveFault_A/seis2_1_0.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/CurveFault_B/seis6_1_0.npy\",\n\n    \"/kaggle/input/waveform-inversion/train_samples/Style_A/data/data1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/Style_B/data/data1.npy\",\n    \n]\nvel_files = [\n    \"/kaggle/input/waveform-inversion/train_samples/FlatVel_A/model/model1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/FlatVel_B/model/model1.npy\",\n    \n    \"/kaggle/input/waveform-inversion/train_samples/CurveVel_A/model/model1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/CurveVel_B/model/model1.npy\",\n\n    \"/kaggle/input/waveform-inversion/train_samples/FlatFault_A/vel2_1_0.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/FlatFault_B/vel6_1_0.npy\",\n    \n    \"/kaggle/input/waveform-inversion/train_samples/CurveFault_A/vel2_1_0.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/CurveFault_B/vel6_1_0.npy\",\n\n    \"/kaggle/input/waveform-inversion/train_samples/Style_A/model/model1.npy\",\n    \"/kaggle/input/waveform-inversion/train_samples/Style_B/model/model1.npy\",\n    \n]\ntypes = [\n    'FlatVel_A', 'FlatVel_B', 'CurveVel_A', 'CurveVel_B', 'FlatFault_A', 'FlatFault_B', 'CurveFault_A', 'CurveFault_B', 'Style_A', 'Style_B'\n]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T06:18:44.116382Z","iopub.execute_input":"2025-04-17T06:18:44.116829Z","iopub.status.idle":"2025-04-17T06:18:44.122891Z","shell.execute_reply.started":"2025-04-17T06:18:44.116764Z","shell.execute_reply":"2025-04-17T06:18:44.121579Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nfor i, tp in enumerate(types):\n    seis_data = np.load(seis_files[i])\n    vel_data = np.load(vel_files[i])\n    seis_data_sim = vel_to_seis(vel_data[0][0])\n    print(f\"TYPE : {tp}\")\n    print(\"\\n< VELOCITY MAP >\")\n    fig, ax = plt.subplots(1, 1, figsize=(5, 2.5))\n    ax.imshow(vel_data[0, 0])\n    plt.show()\n\n    print(f\"\\n< FORWARD SIMULATION (from velocity map) : shape={seis_data_sim.shape} >\")\n    plot_seis(seis_data_sim[:, 2:, :])\n    \n    print(\"\\n< INPUT (origin) >\")\n    plot_seis(seis_data[0])\n    \n    print(\"\\n< ERROR >\")\n    errors = []\n    for j in range(5):\n        error = np.mean(np.abs(seis_data[0, j, :, :] - seis_data_sim[j, 2:, :])) # MAE\n        errors.append(error)\n        print(f\"Receiver{j+1} Error : {error:.6f}\")\n    print(f\"Mean Error : {np.mean(errors):.6f}\")\n    print(\"#########################################\")\n    print()\n\n    # if i==2:\n    #     break\n    # break","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-16T21:24:35.256211Z","iopub.execute_input":"2025-04-16T21:24:35.256630Z","iopub.status.idle":"2025-04-16T21:27:16.602626Z","shell.execute_reply.started":"2025-04-16T21:24:35.256602Z","shell.execute_reply":"2025-04-16T21:27:16.601147Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"------\n## **If there is noise..**\nAssuming that there is noise in the velocity map, it will be similar to the situation where our predictive model is mispredicting the velocity map.  \nOf course, you can see that the error in seis_data has increased when comparing the same types in the previous results.  \nHowever, it is somewhat unreasonable to estimate the predictive quality of the Velocity map by the error level in seis_data.  ","metadata":{}},{"cell_type":"markdown","source":"#### **INPUT (origin) vs FORWARD SIMULATION (from noise velocity map)**","metadata":{}},{"cell_type":"code","source":"%%time\nfor i, tp in enumerate(types):\n    seis_data = np.load(seis_files[i])\n    vel_data = np.load(vel_files[i])    \n\n    noise = np.random.normal(loc=(np.min(vel_data) + np.max(vel_data))/500, scale=np.std(vel_data)/15, size=vel_data[0][0].shape)\n    vel_data_noise = vel_data[0, 0] + noise\n\n    seis_data_sim = vel_to_seis(vel_data_noise)   \n\n    print(f\"TYPE : {tp}\")\n    print(f\"\\n< VELOCITY MAP (origin vs noise) : error={np.mean(np.abs(vel_data[0, 0] - vel_data_noise)):.2f}>\")\n    fig, ax = plt.subplots(1, 2, figsize=(5, 2.5))\n    ax[0].imshow(vel_data[0, 0])\n    ax[1].imshow(vel_data_noise)\n    plt.show()\n\n    print(f\"\\n< FORWARD SIMULATION (from noise velocity map) : shape={seis_data_sim.shape} >\")\n    plot_seis(seis_data_sim)\n    \n    print(\"\\n< INPUT (origin) >\")\n    plot_seis(seis_data[0])\n    \n    print(\"\\n< ERROR >\")\n    errors = []\n    for j in range(5):\n        error = np.mean(np.abs(seis_data[0, j, :, :] - seis_data_sim[j, 2:, :])) # MAE\n        errors.append(error)\n        print(f\"Receiver{j+1} Error : {error:.6f}\")\n    print(f\"Mean Error : {np.mean(errors):.6f}\")\n    print(\"#########################################\")\n    print()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-16T14:49:20.537655Z","iopub.execute_input":"2025-04-16T14:49:20.537983Z","iopub.status.idle":"2025-04-16T14:49:54.743669Z","shell.execute_reply.started":"2025-04-16T14:49:20.537955Z","shell.execute_reply":"2025-04-16T14:49:54.742579Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#### **FORWARD SIMULATION (from velocity map) vs FORWARD SIMULATION (from noise velocity map)**","metadata":{}},{"cell_type":"code","source":"%%time\nfor i, tp in enumerate(types):\n    seis_data = np.load(seis_files[i])\n    vel_data = np.load(vel_files[i])   \n    seis_data_sim = vel_to_seis(vel_data[0][0])\n\n    noise = np.random.normal(loc=(np.min(vel_data) + np.max(vel_data))/500, scale=np.std(vel_data)/15, size=vel_data[0][0].shape)\n    vel_data_noise = vel_data[0, 0] + noise\n\n    seis_data_sim_noise = vel_to_seis(vel_data_noise)   \n\n    print(f\"TYPE : {tp}\")\n    print(f\"\\n< VELOCITY MAP (origin vs noise) : error={np.mean(np.abs(vel_data[0, 0] - vel_data_noise)):.2f}>\")\n    fig, ax = plt.subplots(1, 2, figsize=(5, 2.5))\n    ax[0].imshow(vel_data[0, 0])\n    ax[1].imshow(vel_data_noise)\n    plt.show()\n\n    print(f\"\\n< FORWARD SIMULATION (from noise velocity map) : shape={seis_data_sim.shape} >\")\n    plot_seis(seis_data_sim_noise)\n    \n    print(\"\\n< FORWARD SIMULATION (from noise velocity map) >\")\n    plot_seis(seis_data_sim)\n    \n    print(\"\\n< ERROR >\")\n    errors = []\n    for j in range(5):\n        error = np.mean(np.abs(seis_data_sim[j, 2:, :] - seis_data_sim_noise[j, 2:, :])) # MAE\n        errors.append(error)\n        print(f\"Receiver{j+1} Error : {error:.6f}\")\n    print(f\"Mean Error : {np.mean(errors):.6f}\")\n    print(\"#########################################\")\n    print()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T06:19:05.965003Z","iopub.execute_input":"2025-04-17T06:19:05.965439Z","iopub.status.idle":"2025-04-17T06:19:06.473845Z","shell.execute_reply.started":"2025-04-17T06:19:05.965408Z","shell.execute_reply":"2025-04-17T06:19:06.472738Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## **Conclusion**\n1. We can convert the simulation code to a Python version to recover the seis data to a small level of error.\n2. Visually, the simulated seis data is very similar to the seis data of openFWI.\n3. However, due to the nature of the physical field, this error can be very large.\n4. It should be noted that the level of error between the original seis data and the restored seis data is different depending on the characteristics of the fault.","metadata":{}}]}