{"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":31040,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Yale/UNC-CH - Geophysical Waveform Inversion\n## Develop physics-guided machine learning models to solve full-waveform inversion problems","metadata":{}},{"cell_type":"markdown","source":"# Add Libraries","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nfrom matplotlib.colors import ListedColormap\nimport pandas as pd","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:27.690203Z","iopub.execute_input":"2025-05-30T14:56:27.691252Z","iopub.status.idle":"2025-05-30T14:56:28.081565Z","shell.execute_reply.started":"2025-05-30T14:56:27.691216Z","shell.execute_reply":"2025-05-30T14:56:28.080426Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# To show and process the Velocity Maps","metadata":{}},{"cell_type":"code","source":"\n# Plot the velocity map\ndef plot_velocity(velocity, sample):\n    fig, ax = plt.subplots(1, 1, figsize=(11, 5))\n    img=ax.imshow(velocity[sample,0,:,:],cmap='jet')\n    ax.set_xticks(range(0, 70, 10))\n    ax.set_xticklabels(range(0, 700, 100))\n    ax.set_yticks(range(0, 70, 10))\n    ax.set_yticklabels(range(0, 700, 100))\n    ax.set_ylabel('Depth (m)', fontsize=12)\n    ax.set_xlabel('Offset (m)', fontsize=12)\n    clb=plt.colorbar(img, ax=ax)\n    clb.ax.set_title('km/s',fontsize=8)\n    plt.show()\n    # And a simple ave velocity vs depth plot\n    plt.figure(figsize=(8, 2.5))\n    plt.plot(np.arange(5,700,10), np.mean(velocity[sample,0,:,:],axis=1))\n    plt.xlabel(\"Depth (m)\")\n    plt.ylabel(\"Ave Velocity (m/s)\")\n    plt.show()\n\n# Get information from the velocity map\ndef info_velocity(velocity, sample, for_show=True):\n    # When for_show=True display results and plots.\n    # When for_show=False work silently and return measured values.\n    # Indices are: sample, 0, depth, xloc\n    ave_vel = np.mean(velocity[sample,0,:,:])\n    std_vel = np.std(velocity[sample,0,:,:])\n    min_vel = np.min(velocity[sample,0,:,:])\n    max_vel = np.max(velocity[sample,0,:,:])\n    medi_vel = np.median(velocity[sample,0,:,:])\n    MAE_1medi = np.mean(np.abs(velocity[sample,0,:,:] - medi_vel))\n    # Number of unique velocities\n    num_vels = len(np.unique(velocity[isample,0,:,:]))\n    # Average velocities in first row halves ~ surface velocity\n    y0_velL = np.mean(velocity[sample,0, 0 , 0:35  ])\n    y0_velR = np.mean(velocity[sample,0, 0 , 35:  ])\n    # Median velocities in rows 0-9, 10-29, 30-49, 50-69\n    # and keep track of MAE wrt to these\n    MAE_5medi = 0.0\n    y09L_medi = np.median(velocity[sample,0, 0:10 , 0:34+1  ])\n    MAE_5medi += 5.0*np.mean(np.abs(velocity[sample,0, 0:10 , 0:34+1  ] - y09L_medi))\n    y09R_medi = np.median(velocity[sample,0, 0:10 , 35:  ])\n    MAE_5medi += 5.0*np.mean(np.abs(velocity[sample,0, 0:10 , 35:  ] - y09R_medi))\n    y1029_medi = np.median(velocity[sample,0, 10:29+1 , :  ])\n    MAE_5medi += 20.0*np.mean(np.abs(velocity[sample,0, 10:29+1 , :  ] - y1029_medi))\n    y3049_medi = np.median(velocity[sample,0, 30:49+1 , :  ])\n    MAE_5medi += 20.0*np.mean(np.abs(velocity[sample,0, 30:49+1 , :  ] - y3049_medi))\n    y5069_medi = np.median(velocity[sample,0, 50: , :  ])\n    MAE_5medi += 20.0*np.mean(np.abs(velocity[sample,0, 50: , :  ] - y5069_medi))\n    MAE_5medi = MAE_5medi / 70.0\n    # Means\n    y09L_mean = np.mean(velocity[sample,0, 0:10 , 0:34+1  ])\n    y09R_mean = np.mean(velocity[sample,0, 0:10 , 35:  ])\n    y1029_mean = np.mean(velocity[sample,0, 10:29+1 , :  ])\n    y3049_mean = np.mean(velocity[sample,0, 30:49+1 , :  ])\n    y5069_mean = np.mean(velocity[sample,0, 50: , :  ])\n    if for_show:\n        print(\"Number of distinct velocities: {}\".format(num_vels))\n        print(\"Average velocity: {:.2f} m/s  SD: {:.2f}\".format(ave_vel, std_vel))\n        print(\"Median velocity: {:.2f} m/s\".format(medi_vel),\n             \"   Min, Max: {:.2f}, {:.2f}\".format(min_vel, max_vel))\n        print(\"MAE from median: {:.2f}  \".format(MAE_1medi))\n        print(\"Ave y=0 velocities L,R: {:.2f}, {:.2f}\".format(y0_velL, y0_velR))\n        print(\"Median velocities in rows:  {:.2f}(0-9:L), {:.2f}(0-9:R),\".format(\n                y09L_medi, y09R_medi),\n                \"{:.2f}(10-29), {:.2f}(30-49), {:.2f}(50-69)\".format(\n                y1029_medi, y3049_medi, y5069_medi))\n        print(\"MAE from 5 medians: {:.2f}\".format(MAE_5medi))\n        print(\"  Mean velocities in rows:  {:.2f}(0-9:L), {:.2f}(0-9:R),\".format(\n                y09L_mean, y09R_mean),\n                \"{:.2f}(10-29), {:.2f}(30-49), {:.2f}(50-69)\".format(\n                y1029_mean, y3049_mean, y5069_mean))\n        \n    else:\n        return (num_vels, y0_velL, y0_velR, y09L_medi, y09R_medi,\n                    y1029_medi, y3049_medi, y5069_medi, MAE_1medi, MAE_5medi)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:28.083443Z","iopub.execute_input":"2025-05-30T14:56:28.084520Z","iopub.status.idle":"2025-05-30T14:56:28.101817Z","shell.execute_reply.started":"2025-05-30T14:56:28.084431Z","shell.execute_reply":"2025-05-30T14:56:28.100730Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# To show and process the Seismic Data","metadata":{}},{"cell_type":"code","source":"# Make a gray-scale image of the seismic data\ndef plot_data(data, sample=-1):\n    fig,ax=plt.subplots(1,5,figsize=(20,7))\n    # Is it a Train (multiple) or Test (single) data?\n    if len(data.shape) == 3: \n        thisdata = data[:,:,:]\n    else:\n        thisdata = data[sample,:,:,:]\n    # Scale the color range. Use symmetric values to have 0 in the middle.\n    # Use the values in the source columns to avoid source pulses.\n    maxabs = []\n    for srclocid, xloc in enumerate([0,17,34,52,69]):\n        maxabs.append(np.max(np.abs(thisdata[srclocid,180:,xloc])))\n    vrange = np.max(maxabs) * 0.5  # use less than max, some saturation is OK\n    for iax in range(5):\n        ax[iax].imshow( thisdata[iax,:,:], extent=[0,70,1000,0],\n                       aspect='auto', cmap='gray', vmin=-vrange, vmax=vrange)\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()\n\n\n# Get an accurate time of the max (usually first) peak from given source in given xloc.\ndef time_max_peak(isrc, xloc, thisdata):\n    ipeak = np.argmax(thisdata[isrc,:,xloc])\n    # fit 7 points with degree=2\n    peakvals = thisdata[isrc, ipeak-3:ipeak+3+1, xloc]\n    timevals = np.linspace(ipeak-3,ipeak+3, num=7, endpoint=True)\n    if len(peakvals) == len(timevals):\n        fitcoefs = np.poly1d(np.polyfit(timevals, peakvals, 2)).coef\n        # max is at -b/(2a)   :)\n        return -0.5*fitcoefs[1]/fitcoefs[0]\n    else:\n        print(\"mis-matched lengths:\\n\",timevals, \"\\n\", peakvals)\n        return ipeak\n\n\n# Get information from the seismic data\n# When for_show=True display results and plots.\n# When for_show=False work silently and return measured values.\ndef info_data(data, sample=-1, for_show=True):\n    # Train (multiple) or Test (single) data?\n    if len(data.shape) == 3:\n        thisdata = data[:,:,:]\n    else:\n        thisdata = data[sample,:,:,:]\n    # Calculate the surface velocity, don't use source columns\n    # Ignore peaks that are too distant to be the surface peak at 200 m away from source:\n    # Times less than: 225 ms. From: 100 ms (source peak) + 100 m / 800 m/s * 1000 s/ms.\n    # Still have issues with quick reflections messing up the velocities -\n    # use short baselines and average 4 on each side.\n    partdata = thisdata[ : , 0:225 ,: ]\n    vsurfaceL = 4*5*10*1000/( time_max_peak(0, 6, partdata) - time_max_peak(0, 1, partdata) +\n                            time_max_peak(1, 11, partdata) - time_max_peak(1, 16, partdata) +\n                            time_max_peak(1, 23, partdata) - time_max_peak(1, 18, partdata) +\n                            time_max_peak(2, 28, partdata) - time_max_peak(2, 33, partdata) )\n    vsurfaceR = 4*5*10*1000/( time_max_peak(2, 40, partdata) - time_max_peak(2, 35, partdata) +\n                            time_max_peak(3, 46, partdata) - time_max_peak(3, 51, partdata) +\n                            time_max_peak(3, 58, partdata) - time_max_peak(3, 53, partdata) +\n                            time_max_peak(4, 63, partdata) - time_max_peak(4, 68, partdata) )\n    # Clip the average velocity to 4100 - don't trust higher values are real.\n    vsurface = np.clip((vsurfaceL + vsurfaceR)/2, 1400.0, 4100.0)\n    if for_show:\n        # Make a plot of surface wave distance vs time\n        # Use the middle source location to avoid/reduce reflected wave interference\n        idists = np.arange(0,70)\n        dists = []\n        times = []\n        # Time is relative to src 2 peak time.\n        timeref = np.argmax(thisdata[2,:,34])\n        for idist in idists:\n            dists.append(10*idist) # 0 to 690\n            # signed time for before/after xloc=34\n            times.append(np.sign(idist-34)*(np.argmax(thisdata[2,:,idist])-timeref)/1000.0)\n        times = -1.0*(times - times[0])  # adjust orientation of time axis\n        plt.figure(figsize=(6,3))\n        plt.plot(dists, times, '.b', alpha=0.7)\n        plt.plot([dists[0],dists[-1]],[times[0],times[-1]],c='orange',alpha=0.6)\n        plt.ylabel(\"$-$ Time (s)\")\n        plt.xlabel(\"Surface Distance (m)\")\n        plt.title(\"Time vs Distance from Source 2\")\n        plt.show()\n        ##v_orange = -1.0*(dists[-1] - dists[0])/(times[-1] - times[0])\n        print(\"Surface velocities : {:.2f}-Left, {:.2f}-Average, {:.2f}-Right\".format(\n                        vsurfaceL, vsurface, vsurfaceR))\n    else:\n        return vsurfaceL, vsurface, vsurfaceR","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:28.102887Z","iopub.execute_input":"2025-05-30T14:56:28.103230Z","iopub.status.idle":"2025-05-30T14:56:28.140679Z","shell.execute_reply.started":"2025-05-30T14:56:28.103197Z","shell.execute_reply":"2025-05-30T14:56:28.139594Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Make a plot of the waveforms at each of the source locations when that source is active.\ndef sources_data(data, sample=-1, for_show=True):\n    # The 5 sources are located closest to: 0, 17, 34, 52, 69\n    # The peak amplitude ~ 40 for each.\n    # Train (multiple) or Test (single) data?\n    if len(data.shape) == 3:\n        thisdata = data[:,:,:]\n    else:\n        thisdata = data[sample,:,:,:]\n    # Get the max, min amplitudes for t > 180 for each source-location\n    maxamps = []\n    minamps = []\n    for srclocid, xloc in enumerate([0,17,34,52,69]):\n        ##print(\"Source peak at xloc={} is: {:.2f}\".format(\n        ##        xloc, np.max(thisdata[srclocid,:,xloc]) ))\n        # Max and min after the source peak\n        maxamps.append(np.max(thisdata[srclocid,180:,xloc]))\n        minamps.append(np.min(thisdata[srclocid,180:,xloc]))\n    max_amp = np.max(maxamps)\n    min_amp = np.min(minamps)\n    delta_amp = 0.05*(max_amp - min_amp)\n    plt.figure(figsize=(8,5))\n    for srclocid, xloc in enumerate([0,17,34,52,69]):\n        timeseries = thisdata[srclocid,:,xloc]  # srclocid, time, xloc\n        offset = delta_amp*(xloc - 34)/35.0\n        plt.plot(np.array(range(1000)) + 0*offset, timeseries + offset, alpha=0.7) \n    plt.plot([0,1000],[0.0,0.0],c='gray',alpha=0.5)\n    plt.ylim(1.10*min_amp - delta_amp, 1.10*max_amp + delta_amp)\n    plt.xlabel('Time (ms)')\n    plt.ylabel(\"Amplitude     Traces are offset.\")\n    plt.title(\"Waveforms at the 5 source locations\")\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:28.142856Z","iopub.execute_input":"2025-05-30T14:56:28.143213Z","iopub.status.idle":"2025-05-30T14:56:28.172827Z","shell.execute_reply.started":"2025-05-30T14:56:28.143180Z","shell.execute_reply":"2025-05-30T14:56:28.171631Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Routine to read in the training data files given the training dataframe index value","metadata":{}},{"cell_type":"code","source":"# There are two directory formats for getting the data-velocity pairs depending on the type:\n\n# FlatVel_[A,B], CurveVel_[A,B], Style_[A|B]\n# Each of these 6 dirs contain: /data/data[1,2].npy and /model/model[1,2].npy\n# Total # of velocity-meaurement pairs: 12 x 500\n##velocity = np.load('/kaggle/input/waveform-inversion/train_samples/FlatVel_A/model/model2.npy')\n##data = np.load('/kaggle/input/waveform-inversion/train_samples/FlatVel_A/data/data2.npy')\n##isample = 13 \n\n# [FlatFault,CurveFault]_A has files: seis[2,4]_1_0.npy, vel[2,4]_1_0.npy\n# [FlatFault,CurveFault]_B has files: seis[6,8]_1_0.npy, vel[6,8]_1_0.npy\n# Total # of velocity-meaurement pairs: 8 x 500\n##velocity = np.load('/kaggle/input/waveform-inversion/train_samples/CurveFault_A/vel4_1_0.npy')\n##data = np.load('/kaggle/input/waveform-inversion/train_samples/CurveFault_A/seis4_1_0.npy')\n##isample = 23\n\n# Keep track of the last train file read in to avoid re-reading when just isample changes\nlast_data_file = \"None\"\n\ndef get_train_sample(dfind, ftscale=True):\n    # Assumes traindf is defined.  And uses global values:\n    global velocity, data, last_data_file\n    train_dir = \"/kaggle/input/waveform-inversion/train_samples/\"\n    veltype, ifile, isample = traindf.loc[dfind, [\"veltype\",\"ifile\",\"isample\"]]\n    if (\"Vel\" in veltype) or (\"Style\" in veltype):\n        data_file = train_dir+veltype+\"/data/data\"+str(ifile)+\".npy\"\n        model_file = train_dir+veltype+\"/model/model\"+str(ifile)+\".npy\"\n        ##print(\"got Vel or Style:\\n   \", data_file, \"\\n   \", model_file)\n    else:  # it is a Fault type\n        fault_num = 2*ifile + 4*(\"_B\" in veltype)\n        data_file = train_dir+veltype+\"/seis\"+str(fault_num)+\"_1_0.npy\"\n        model_file = train_dir+veltype+\"/vel\"+str(fault_num)+\"_1_0.npy\"\n        ##print(\"got Fault:\\n   \", data_file, \"\\n   \", model_file)\n    # Read them in if not already available\n    if data_file != last_data_file:\n            data = np.load(data_file)\n            # Scale the seismic data as a function of time:\n            if ftscale:\n                for itime in range(1000):\n                    data[ : , : , itime, : ] = (1.0+(itime/200)**1.5)*data[ : , : , itime, : ]\n            velocity = np.load(model_file)\n            last_data_file = data_file\n    return velocity, data, isample\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:28.174210Z","iopub.execute_input":"2025-05-30T14:56:28.174530Z","iopub.status.idle":"2025-05-30T14:56:28.196388Z","shell.execute_reply.started":"2025-05-30T14:56:28.174443Z","shell.execute_reply":"2025-05-30T14:56:28.195222Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Convert many x,y points into a quartile-based set of x_median,y_median points.","metadata":{}},{"cell_type":"code","source":"# Fitting these median points is similar to fitting x,y with a MAE metric.\ndef xy_medians(xin, yin, nqs):\n    # Outputs x,y medians from about nqs quartiles\n    sortinds = np.argsort(xin)\n    xsort = xin[sortinds]\n    ysort = yin[sortinds]\n    lenxs = len(xsort)\n    nsample = int(lenxs/nqs)\n    # Have a first and last range of ~ nsample/3 points\n    nfirstlast = int(nsample/4)\n    indups = list(range(nfirstlast, lenxs - nfirstlast, nsample))\n    indups.insert(0,0) # start with 0\n    indups.append(lenxs - nfirstlast)\n    indups.append(lenxs)\n    ##print(indups)\n    xmeds = []; ymeds = []\n    for iup in range(0, len(indups)-1):\n        indlow = indups[iup]\n        indhi = indups[iup+1]\n        xmed = np.median(xsort[indlow:indhi])\n        ymed = np.median(ysort[indlow:indhi])\n        xmeds.append(xmed)\n        ymeds.append(ymed)\n    # Add a value at xmax: average of last value and linear trend by quantile\n    xmeds.append(xsort[-1])\n    ymeds.append(ymeds[-1] + 0.5*(ymeds[-1] - ymeds[-2]))\n    return xmeds, ymeds","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:28.197595Z","iopub.execute_input":"2025-05-30T14:56:28.197930Z","iopub.status.idle":"2025-05-30T14:56:28.227666Z","shell.execute_reply.started":"2025-05-30T14:56:28.197885Z","shell.execute_reply":"2025-05-30T14:56:28.226782Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# There are 10,000 training samples on kaggle, organized as: 10 x 2 x 500 data-vel pairs","metadata":{}},{"cell_type":"markdown","source":"##!ls /kaggle/input/waveform-inversion/train_samples/*","metadata":{}},{"cell_type":"code","source":"# Make a dataframe with 10,000 rows labeled by:\n#   type - 5 x 2 string values\n#   ifile - two numeric values: 0,1 or 1,2 or 2,4 or 6,8 depending on type\n#   isample - 0 to 499\nveltypes = [\"FlatVel\",\"FlatFault\", \"CurveVel\", \"CurveFault\", \"Style\"]\nveltype = []; ifile = []; isample = []\nfor this_type in veltypes:\n    for this_AB in [\"_A\",\"_B\"]:\n        for this_ifile in [1,2]:\n            for this_isample in range(500):  # **************************************\n                veltype.append(this_type+this_AB); ifile.append(this_ifile); isample.append(this_isample)\n# Make a dataframe from these\ntraindf = pd.DataFrame({\"veltype\":veltype, \"ifile\":ifile, \"isample\":isample})","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:28.228946Z","iopub.execute_input":"2025-05-30T14:56:28.229293Z","iopub.status.idle":"2025-05-30T14:56:28.268001Z","shell.execute_reply.started":"2025-05-30T14:56:28.229264Z","shell.execute_reply":"2025-05-30T14:56:28.266353Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"traindf","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:28.269930Z","iopub.execute_input":"2025-05-30T14:56:28.270209Z","iopub.status.idle":"2025-05-30T14:56:28.310632Z","shell.execute_reply.started":"2025-05-30T14:56:28.270189Z","shell.execute_reply":"2025-05-30T14:56:28.309698Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Select a dataframe index to look at","metadata":{}},{"cell_type":"code","source":"dfind = int(0.91*len(traindf))\n\n\nprint(list(traindf.loc[dfind,[\"veltype\",\"ifile\",\"isample\"]]))\nvelocity, data, isample = get_train_sample(dfind)\n\nprint('Velocity map size:', velocity.shape)\nprint('Seismic data size:', data.shape)\n# Velocity map size: (500, 1, 70, 70)   sample, 0, yloc(10m), xloc(10m)\n# Seismic data size: (500, 5, 1000, 70) sample, srclocid, time(ms), xloc(10m)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:28.311537Z","iopub.execute_input":"2025-05-30T14:56:28.311783Z","iopub.status.idle":"2025-05-30T14:56:35.953916Z","shell.execute_reply.started":"2025-05-30T14:56:28.311766Z","shell.execute_reply":"2025-05-30T14:56:35.952827Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Look at the velocity map for the training sample","metadata":{}},{"cell_type":"code","source":"# isample defined above\nplot_velocity(velocity, isample)\ninfo_velocity(velocity, isample)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:35.956948Z","iopub.execute_input":"2025-05-30T14:56:35.957231Z","iopub.status.idle":"2025-05-30T14:56:36.485233Z","shell.execute_reply.started":"2025-05-30T14:56:35.957212Z","shell.execute_reply":"2025-05-30T14:56:36.484396Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Look at the seismic data for the sample","metadata":{}},{"cell_type":"code","source":"# isample = same as for the velocity map above\nplot_data(data, isample)\ninfo_data(data, isample)\nsources_data(data, isample)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:36.486227Z","iopub.execute_input":"2025-05-30T14:56:36.486576Z","iopub.status.idle":"2025-05-30T14:56:37.768490Z","shell.execute_reply.started":"2025-05-30T14:56:36.486547Z","shell.execute_reply":"2025-05-30T14:56:37.767544Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# The test data consists of 65818 data sets to be predicted\n\n# 65818   65818  987270","metadata":{}},{"cell_type":"markdown","source":"##!ls /kaggle/input/waveform-inversion/test | wc","metadata":{}},{"cell_type":"code","source":"##!ls -s /kaggle/input/waveform-inversion/test/c*.npy | head -5\n# total 90039024\n# 1368 000039dca2.npy\n# 1368 0000fd8ec8.npy\n# 1368 0001026c8a.npy\n# 1368 00015b24d5.npy\n#      a00269f1eb.npy\n#      c001726adb.npy\n#      c0021521e5.npy\n\n# Look at one of them (they seem to be shuffled)\n##testdata = np.load('/kaggle/input/waveform-inversion/test/000039dca2.npy')  # messy\n##testdata = np.load('/kaggle/input/waveform-inversion/test/0001026c8a.npy')  # very simple\ntestdata = np.load('/kaggle/input/waveform-inversion/test/00015b24d5.npy')  # weird straight lines\n##testdata = np.load('/kaggle/input/waveform-inversion/test/800222ab0d.npy')  # messy\n##testdata = np.load('/kaggle/input/waveform-inversion/test/a00269f1eb.npy')  # messy\n##testdata = np.load('/kaggle/input/waveform-inversion/test/c0021521e5.npy')  # simple-ish\n\n# Scale the seismic data by ~ (1+(t/a)^b) to help equalize the amplitudes vs time.\n# (This is similar to applying AGC for visualization, but is included in the analysis too.)\nfor itime in range(1000):\n    testdata[ : , itime, : ] = (1.0+(itime/200)**1.5)*testdata[ : , itime, : ]\n\nprint('Test data size:', testdata.shape)\n\nplot_data(testdata)\ninfo_data(testdata)\nsources_data(testdata)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:37.769629Z","iopub.execute_input":"2025-05-30T14:56:37.769902Z","iopub.status.idle":"2025-05-30T14:56:38.930142Z","shell.execute_reply.started":"2025-05-30T14:56:37.769883Z","shell.execute_reply":"2025-05-30T14:56:38.929047Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# For each sample add:\n\n# The y_ targets: y0_velL, y0_velR, y09L_medi, y09R_medi, y1039_medi, y4069_medi\nnunique = []; y0_aves = []; y0_diffs = []\ny09L_medis = []; y09R_medis = []; y1029_medis = []; y3049_medis = []; y5069_medis = []\n\n# These properties of the target velocity map\nMAE_1medis = []; MAE_5medis = []\n\n# The x_ features: surface velocity average and R-L difference\nsurf_aves = []; surf_diffs = []\n\nfor dfind in traindf.index:\n    velocity, data, isample = get_train_sample(dfind, ftscale=False)\n    # velocity, target, values\n    (num_vels, y0_velL, y0_velR, y09L_medi, y09R_medi,\n         y1029_medi, y3049_medi, y5069_medi, MAE_1medi, MAE_5medi) = info_velocity(\n                                            velocity, isample, for_show=False)\n    nunique.append(num_vels)\n    y0_aves.append((y0_velL + y0_velR)/2); y0_diffs.append(y0_velR - y0_velL)\n    y09L_medis.append(y09L_medi); y09R_medis.append(y09R_medi); y1029_medis.append(y1029_medi)\n    y3049_medis.append(y3049_medi); y5069_medis.append(y5069_medi)\n    MAE_1medis.append(MAE_1medi); MAE_5medis.append(MAE_5medi)\n    \n    # features are seismic measured values\n    velL, velave, velR = info_data(data, isample, for_show=False)\n    surf_aves.append(velave)\n    surf_diffs.append(velR-velL)\n\ntraindf[\"y_numVels\"] = nunique\ntraindf[\"y_y0Ave\"] = y0_aves\ntraindf[\"y_y0Diff\"] = y0_diffs\ntraindf[\"y_09LMedi\"] = y09L_medis\ntraindf[\"y_09RMedi\"] = y09R_medis\ntraindf[\"y_1029Medi\"] = y1029_medis\ntraindf[\"y_3049Medi\"] = y3049_medis\ntraindf[\"y_5069Medi\"] = y5069_medis\ntraindf[\"MAE_1Medi\"] = MAE_1medis\ntraindf[\"MAE_5Medi\"] = MAE_5medis\ntraindf[\"x_surfAve\"] = surf_aves\ntraindf[\"x_surfDiff\"] = surf_diffs\n\n# Add color-coding based on the surfDiff and surfAve values\n# Red = R-L not zero; Blue = R-L near zero\ntraindf[\"diff_clr\"] = 'red'\n# Use measured R-L difference to set color\nseldiff = traindf[\"x_surfAve\"] > (1300.0 + 1200.0*np.log10(1+np.abs(traindf[\"x_surfDiff\"])))\n# Use known R-L difference from the target (can't do this for test)\n##seldiff = (np.abs(traindf[\"y_y0Diff\"]) < 0.1*diff_color_change)\ntraindf.loc[seldiff, \"diff_clr\"] = 'blue'\n\ntraindf","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:56:38.931389Z","iopub.execute_input":"2025-05-30T14:56:38.932168Z","iopub.status.idle":"2025-05-30T14:58:49.017161Z","shell.execute_reply.started":"2025-05-30T14:56:38.932132Z","shell.execute_reply":"2025-05-30T14:58:49.015825Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Summary values for the columns\ntraindf_means = traindf.describe().loc[\"mean\",]\n##traindf.describe()\n\n# Some summary information\n# Number of discrete velocities:\n# 6003 samples have 2 through 16 (except 9)\n# 3997 sampes have 41 and above.\nnp.clip(traindf[\"y_numVels\"],0,40).value_counts()\n\nprint(\"\\nAverage MAE wrt the median of each sample: {:.2f}\".format(\n            traindf_means[\"MAE_1Medi\"]))\nprint(\"Average MAE wrt 5 medians in each sample: {:.2f}\\n\".format(\n            traindf_means[\"MAE_5Medi\"]))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:58:49.019012Z","iopub.execute_input":"2025-05-30T14:58:49.019373Z","iopub.status.idle":"2025-05-30T14:58:49.075894Z","shell.execute_reply.started":"2025-05-30T14:58:49.019343Z","shell.execute_reply":"2025-05-30T14:58:49.074694Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Save the Training Dataframe -- Save it further along after predictions are added.","metadata":{}},{"cell_type":"code","source":"\ntraindf.to_csv(\"traindf.csv\", header=True, index=False, float_format='%.2f')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:58:49.077155Z","iopub.execute_input":"2025-05-30T14:58:49.077486Z","iopub.status.idle":"2025-05-30T14:58:49.306095Z","shell.execute_reply.started":"2025-05-30T14:58:49.077434Z","shell.execute_reply":"2025-05-30T14:58:49.305215Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# For the plots\nvelocity_range = (1400,4200)\n\n\nprint(\"\\nMedian Ave Surface Velocity: {:.2f}\".format(np.median(traindf[\"x_surfAve\"])))\nprint(\"Average Ave Surface Velocity: {:.2f}\\n\".format(np.mean(traindf[\"x_surfAve\"])))\n\ndiffs = traindf[\"x_surfDiff\"]\n\nplt.figure(figsize=(8,4))\nplt.hist(traindf[\"x_surfAve\"],bins=100)\nplt.title(\"Train: Histogram of the Average Surface Velocity\")\nplt.xlabel(\"Surface Velocity (m/s)\")\nplt.xlim(velocity_range)\nplt.savefig(\"train_hist_surface_velocity.png\")\nplt.show()\n\nplt.figure(figsize=(8,4))\nplt.hist(np.sign(diffs)*np.log10(np.abs(diffs) + 1.0), log=True, bins=100)\nplt.title(\"Train: Histogram of the R-L Velocity Difference\")\nplt.xlabel(\"Signed Log10[1+ R-L Surface Velocity Difference (m/s) ]\")\nplt.savefig(\"train_hist_velocity_difference.png\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:58:49.307219Z","iopub.execute_input":"2025-05-30T14:58:49.307551Z","iopub.status.idle":"2025-05-30T14:58:50.410958Z","shell.execute_reply.started":"2025-05-30T14:58:49.307520Z","shell.execute_reply":"2025-05-30T14:58:50.409823Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(8,5))\nplt.scatter( np.sign(diffs)*(np.log10(np.abs(diffs) + 1.0)), traindf[\"x_surfAve\"],\n                             color=traindf[\"diff_clr\"], s=2, alpha=0.25)\nlindiffs = np.linspace(-3.0,3.0,100)  # <-- This is log10(1+ abs(surfDiff) )\nplt.plot(lindiffs, 1300.0 + 1200.0*np.abs(lindiffs),c='gray',alpha=0.5)\nplt.ylabel(\"Average Surface Velocity (m/s)\")\nplt.xlabel(\"Signed Log10[1+ R-L Velocity Difference (m/s) ]\")\nplt.title(\"Train: Average Surface Velocity vs. R-L Velocity Difference\")\nplt.ylim(velocity_range)\nplt.savefig(\"train_scatter_velocity_vs_difference.png\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:58:50.412341Z","iopub.execute_input":"2025-05-30T14:58:50.412758Z","iopub.status.idle":"2025-05-30T14:58:51.175644Z","shell.execute_reply.started":"2025-05-30T14:58:50.412723Z","shell.execute_reply":"2025-05-30T14:58:51.174584Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Look into the y=0 Average and Difference","metadata":{}},{"cell_type":"code","source":"# Scatter plot of the y=0 row Average and y=0 R-L Difference velocities\nif True:\n    diffs = traindf[\"y_y0Diff\"]\n\n    plt.figure(figsize=(6,3))\n    plt.scatter( np.sign(diffs)*(np.log10(np.abs(diffs) + 1.0)), traindf[\"y_y0Ave\"]/1000,\n                             color=traindf[\"diff_clr\"], s=2, alpha=0.25)\n    plt.ylabel(\"Ave y=0 Velocity (km/s)\")\n    plt.xlabel(\"Signed Log10[1+ y=0 R-L Velocity Diff (m/s) ]\")\n    plt.title(\"Train: y=0 Average Velocity vs. y=0 R-L Velocity Difference\")\n    plt.savefig(\"train_y0_scatter_velocity_vs_diff.png\")\n    plt.ylim(1.4,4.6) # in km/s\n    plt.show()\n\n    # Histogra of the y=0 R-L Diff\n    plt.figure(figsize=(6,3))\n    plt.hist(np.sign(diffs)*np.log10(np.abs(diffs) + 1.0), log=True, bins=100)\n    plt.title(\"Train: Histogram of the y=0 R-L Velocity Difference\")\n    plt.xlabel(\"Signed Log10[1+ R-L y=0 Velocity Difference (m/s) ]\")\n    plt.savefig(\"train_hist_y0_difference.png\")\n    plt.show()\n\n    # Scatter plot of the Seismic R-L Diff vs the y=0 R-L Diff\n    diffs = traindf[\"x_surfDiff\"]\n    diffy0 = traindf[\"y_y0Diff\"]\n\n    plt.figure(figsize=(6,3))\n    plt.scatter( np.sign(diffy0)*(np.log10(np.abs(diffy0) + 1.0)),\n                    np.sign(diffs)*(np.log10(np.abs(diffs) + 1.0)),\n                             color=traindf[\"diff_clr\"], s=2, alpha=0.25)\n    plt.xlabel(\"y=0  Log10[1+ R-L Velocity Diff (m/s) ]\")\n    plt.ylabel(\"Seismic  Log10[1+ R-L Velocity Diff (m/s) ]\")\n    plt.title(\"Train: Seismic R-L Difference vs the y=0 R-L Difference\")\n    plt.savefig(\"train_scatter_diff_vs_diff.png\")\n    plt.show()\n\n# Find some with y=0 diff = 0 and yet seismic R-L is high\n##traindf[(traindf[\"y_y0Diff\"] == 0) & (traindf[\"x_surfDiff\"] > 200)]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:58:51.176736Z","iopub.execute_input":"2025-05-30T14:58:51.177046Z","iopub.status.idle":"2025-05-30T14:58:53.364922Z","shell.execute_reply.started":"2025-05-30T14:58:51.177021Z","shell.execute_reply":"2025-05-30T14:58:53.363682Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Compare measured surface velocity with the y=0 average.","metadata":{}},{"cell_type":"code","source":"# Include a simple degree 1 polynomial fit \nmodel = np.poly1d(np.polyfit(np.array(traindf[\"y_y0Ave\"]), \n                             np.array(traindf[\"x_surfAve\"]), 1))\n# for polynomial line visualization \npolyline = np.linspace(1400, 4500, 100)  \n\nplt.figure(figsize=(4,4))\nplt.scatter( traindf[\"y_y0Ave\"], traindf[\"x_surfAve\"],\n                color=traindf[\"diff_clr\"], s=2, alpha=0.25)\nplt.plot(polyline, model(polyline), c='orange',alpha=0.6)\nplt.xlabel(\"y=0 Average Velocity\")\nplt.ylabel(\"Seismic Average Surface Velocity (m/s)\")\nplt.title(\"Train: Seismic Surface Velocity vs. y=0 Velocity\")\nplt.xlim(1400,4600)\nplt.ylim(velocity_range)\nplt.savefig(\"train_surf_vs_y0.png\")\nplt.show()\n\nprint(\"   Fit coefs [slope, intercept]:\", model.coef,\"\\n\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:58:53.366271Z","iopub.execute_input":"2025-05-30T14:58:53.366613Z","iopub.status.idle":"2025-05-30T14:58:54.078883Z","shell.execute_reply.started":"2025-05-30T14:58:53.366589Z","shell.execute_reply":"2025-05-30T14:58:54.077843Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Look at the velocities in row ranges vs the surface velocity and velocity difference.","metadata":{}},{"cell_type":"code","source":"# Create simple model fits for each region.\n\nsurfAves = traindf[\"x_surfAve\"]\nsurfDiffs = traindf[\"x_surfDiff\"]\nlog_surfDiffs = np.sign(surfDiffs)*(np.log10(np.abs(surfDiffs) + 1.0))\n\n# Fit red, blue separately, limit the surfAve range used\nselblue = (traindf[\"diff_clr\"] == 'blue') & (traindf[\"x_surfAve\"] < 4100)\nselred = (traindf[\"diff_clr\"] == 'red') & (traindf[\"x_surfAve\"] < 4100)\n\n# for polynomial line visualization \npolyline = np.linspace(1400, 4200, 100)\n\n# Save the fit models\nrows_models = []\nfor y_rows in [\"09L\", \"09R\", \"1029\", \"3049\", \"5069\"]:\n    \n    rows_values = traindf[\"y_\"+y_rows+\"Medi\"]\n    surf_values = surfAves.copy()\n    vel_axis_label = \"Ave Surface Velocity (m/s)\"\n    degree = 5\n    # Modify surf_values for the 09L,R data\n    if \"09L\" in y_rows:\n        surf_values = surf_values - 0.5*surfDiffs\n        vel_axis_label = \"L Surface Velocity (m/s)\"\n        degree = 3\n    if \"09R\" in y_rows:\n        surf_values = surf_values + 0.5*surfDiffs\n        vel_axis_label = \"R Surface Velocity (m/s)\"\n        degree = 3\n\n    plt.figure(figsize=(7,4))\n    plt.scatter(surf_values, rows_values, color=traindf[\"diff_clr\"], s=2, alpha=0.25)\n\n    # Blue polynomial fit:\n    if \"09\" in y_rows:\n        # Use combined L and R data for the model, selblue:\n        surf_RLvalues = np.concatenate( ( (surfAves - 0.5*surfDiffs)[selblue], \n                                            (surfAves + 0.5*surfDiffs)[selblue] ) )\n        rows_RLvalues = np.concatenate( ( traindf.loc[selblue,\"y_09LMedi\"], \n                                            traindf.loc[selblue,\"y_09RMedi\"] ) )\n        xmeds, ymeds = xy_medians(surf_RLvalues, rows_RLvalues, 25)\n        plt.scatter(xmeds, ymeds, s=8, alpha=1.0, c='darkblue')\n        model = np.poly1d(np.polyfit(xmeds, ymeds, degree))\n    else:\n        xmeds, ymeds = xy_medians(np.array(surf_values[selblue]),\n                                    np.array(rows_values[selblue]), 25)\n        plt.scatter(xmeds, ymeds, s=8, alpha=1.0, c='darkblue')\n        model = np.poly1d(np.polyfit(xmeds, ymeds, degree))\n    \n    rows_models.append(model)\n    blue_resids = (-1.0*model(np.array(surf_values[selblue])) + \n                             np.array(rows_values[selblue]))\n    print(\"  Blue Fit coefs:\", model.coef)\n    plt.plot(polyline, model(polyline), c='blue',alpha=1.0)\n    #\n    # Red polynomial fit:\n    if \"09\" in y_rows:\n        # Use combined L and R data for the model, selred:\n        surf_RLvalues = np.concatenate( ( (surfAves - 0.5*surfDiffs)[selred], \n                                            (surfAves + 0.5*surfDiffs)[selred] ) )\n        rows_RLvalues = np.concatenate( ( traindf.loc[selred,\"y_09LMedi\"], \n                                            traindf.loc[selred,\"y_09RMedi\"] ) )\n        xmeds, ymeds = xy_medians(surf_RLvalues, rows_RLvalues, 25)\n        plt.scatter(xmeds, ymeds, s=8, alpha=1.0, c='darkred')\n        model = np.poly1d(np.polyfit(xmeds, ymeds, degree))\n    else:\n        xmeds, ymeds = xy_medians(np.array(surf_values[selred]), \n                                     np.array(rows_values[selred]), 25)\n        plt.scatter(xmeds, ymeds, s=8, alpha=1.0, c='darkred')\n        model = np.poly1d(np.polyfit(xmeds, ymeds, degree))\n    rows_models.append(model)\n    red_resids = (-1.0*model(np.array(surf_values[selred])) + \n                             np.array(rows_values[selred]))\n    print(\"  Red Fit coefs:\", model.coef)\n    plt.plot(polyline, model(polyline), c='purple',alpha=1.0)\n    \n    plt.xlabel(vel_axis_label)\n    plt.xlim(1400, 4200) # reduce because of fitting range\n    plt.ylabel(\"y_\"+y_rows+\" Median\")\n    plt.ylim(1400, 4600)\n    plt.title(\"Train: y_\"+y_rows+\" Median vs. Surface Velocity\")\n    plt.savefig(\"train_rows\"+y_rows+\"_vs_average.png\")\n    plt.show()\n\n\n    # Show the residuals vs surface difference for the 09L, 09R\n    if \"09\" in y_rows:\n        plt.figure(figsize=(7,2))\n        plt.scatter( log_surfDiffs[selblue], blue_resids,\n                             color=traindf.loc[selblue,\"diff_clr\"], s=2, alpha=0.25)\n        plt.scatter( log_surfDiffs[selred], red_resids,\n                             color=traindf.loc[selred,\"diff_clr\"], s=2, alpha=0.25)\n        plt.ylim(-1000,1000)\n        plt.xlabel(\"Signed Log10[1+ R-L Velocity Diff (m/s) ]\")\n        plt.ylabel(\"y_\"+y_rows+\" Residuals\")\n        plt.title(\"Train: y_\"+y_rows+\" * Residuals * vs. Surface Difference\")\n        plt.savefig(\"train_residuals\"+y_rows+\"_vs_difference.png\")\n        plt.show()\n        \n    \n    # Show the median values vs surface difference\n    plt.figure(figsize=(7,2))\n    plt.scatter(log_surfDiffs, rows_values,\n                             color=traindf[\"diff_clr\"], s=2, alpha=0.25)\n    plt.xlabel(\"Signed Log10[1+ R-L Velocity Diff (m/s) ]\")\n    plt.ylabel(\"y_\"+y_rows+\" Median\")\n    plt.ylim(1400, 4600)\n    plt.title(\"Train: y_\"+y_rows+\" Median vs. Surface Difference\")\n    plt.savefig(\"train_rows\"+y_rows+\"_vs_difference.png\")\n    plt.show()\n\n    print(\"\\n\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:58:54.080079Z","iopub.execute_input":"2025-05-30T14:58:54.080373Z","iopub.status.idle":"2025-05-30T14:59:02.792836Z","shell.execute_reply.started":"2025-05-30T14:58:54.080351Z","shell.execute_reply":"2025-05-30T14:59:02.791588Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"polyline = np.linspace(1400, 4100, 100)  \nplt.figure(figsize=(6,3))\nfor imod in range(5):\n    plt.plot(polyline, rows_models[2*imod](polyline), c='blue',alpha=0.6)\n    plt.plot(polyline, rows_models[2*imod+1](polyline), c='red',alpha=0.6)\nplt.xlabel(\"Average Surface Velocity (m/s)\")\nplt.ylabel(\"Median of Rows\")\nplt.title(\"Fits of Row-Ranges Medians vs. Surface Velocity\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:59:02.794069Z","iopub.execute_input":"2025-05-30T14:59:02.794382Z","iopub.status.idle":"2025-05-30T14:59:03.133173Z","shell.execute_reply.started":"2025-05-30T14:59:02.794355Z","shell.execute_reply":"2025-05-30T14:59:03.131663Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# What/why are the blue lines in the 1029 and 3049 median vs surface velocity plots?\n# Find the samples in these lines\ntrainblue = traindf[traindf[\"diff_clr\"] == 'blue']\n\nprint(\"\\n\\n  Look for 'blue' samples that have Rows Medians equal to the y=0 Average.\")\nprint(\"  - List the counts of Velocity-Map Types.\")\nprint(\"  - Check the y0Diff values: they are all 0, so vmaps are R-L symmetric.\\n\\n\")\n\nfor yrows in [\"1029\",\"3049\"]:\n    plt.figure(figsize=(6,2))\n    plt.hist(np.clip(trainblue[\"y_\"+yrows+\"Medi\"] - trainblue[\"y_y0Ave\"],-800,800),\n             log=True, bins=160)\n    plt.xlim(-500,500)\n    plt.xlabel(\"Rows \"+yrows+\" Median  -  y=0 Average\")\n    plt.show()\n\n    matchdf = trainblue[np.abs(trainblue[\"y_\"+yrows+\"Medi\"] - trainblue[\"y_y0Ave\"]) < 0.0001]\n    print(matchdf[\"veltype\"].value_counts())\n    print(matchdf[\"y_y0Diff\"].value_counts())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:59:03.134199Z","iopub.execute_input":"2025-05-30T14:59:03.134485Z","iopub.status.idle":"2025-05-30T14:59:04.129126Z","shell.execute_reply.started":"2025-05-30T14:59:03.134433Z","shell.execute_reply":"2025-05-30T14:59:04.128137Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Compare actual MAE with these values:\nprint(\"\\nMAE if predicted the median of each sample: {:.2f}\".format(\n            traindf_means[\"MAE_1Medi\"]))\nprint(\"MAE if predicted the 5 row-range medians in each sample: {:.2f}\\n\".format(\n            traindf_means[\"MAE_5Medi\"]))\n\nplt.figure(figsize=(6,3))\nplt.hist(traindf[\"MAE_5Medi\"],bins=100)\nplt.xlabel(\"MAE of the sample\")\nplt.title(\"Histogram of the MAEs using 5 known medians\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:59:04.130290Z","iopub.execute_input":"2025-05-30T14:59:04.130662Z","iopub.status.idle":"2025-05-30T14:59:04.391230Z","shell.execute_reply.started":"2025-05-30T14:59:04.130640Z","shell.execute_reply":"2025-05-30T14:59:04.390246Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Add model-predicted columns to the training dataframe based on x_surfAve","metadata":{}},{"cell_type":"code","source":"# Model order is blue then red for each rows range.\nsurf_values = traindf[\"x_surfAve\"]\nsurf_diffs = traindf[\"x_surfDiff\"]\nsurf_L_values = surf_values - 0.5*surf_diffs\nsurf_R_values = surf_values + 0.5*surf_diffs\nselblue = traindf[\"diff_clr\"] == 'blue'\n# Blue and Red model for each row range\nfor imod, y_rows in enumerate([\"09L\", \"09R\", \"1029\", \"3049\", \"5069\"]):\n    if y_rows == \"09L\":\n        traindf.loc[selblue,\"pred_\"+y_rows] = rows_models[2*imod](surf_L_values[selblue])\n        traindf.loc[-selblue,\"pred_\"+y_rows] = rows_models[2*imod+1](surf_L_values[-selblue])\n    elif y_rows == \"09R\":\n        traindf.loc[selblue,\"pred_\"+y_rows] = rows_models[2*imod](surf_R_values[selblue])\n        traindf.loc[-selblue,\"pred_\"+y_rows] = rows_models[2*imod+1](surf_R_values[-selblue])\n    else:\n        traindf.loc[selblue,\"pred_\"+y_rows] = rows_models[2*imod](surf_values[selblue])\n        traindf.loc[-selblue,\"pred_\"+y_rows] = rows_models[2*imod+1](surf_values[-selblue])\n# Add model-predicted columns to the training dataframe based on x_surfAve\n# Model order is blue then red for each rows range.\nsurf_values = traindf[\"x_surfAve\"]\nsurf_diffs = traindf[\"x_surfDiff\"]\nsurf_L_values = surf_values - 0.5*surf_diffs\nsurf_R_values = surf_values + 0.5*surf_diffs\nselblue = traindf[\"diff_clr\"] == 'blue'\n# Blue and Red model for each row range\nfor imod, y_rows in enumerate([\"09L\", \"09R\", \"1029\", \"3049\", \"5069\"]):\n    if y_rows == \"09L\":\n        traindf.loc[selblue,\"pred_\"+y_rows] = rows_models[2*imod](surf_L_values[selblue])\n        traindf.loc[-selblue,\"pred_\"+y_rows] = rows_models[2*imod+1](surf_L_values[-selblue])\n    elif y_rows == \"09R\":\n        traindf.loc[selblue,\"pred_\"+y_rows] = rows_models[2*imod](surf_R_values[selblue])\n        traindf.loc[-selblue,\"pred_\"+y_rows] = rows_models[2*imod+1](surf_R_values[-selblue])\n    else:\n        traindf.loc[selblue,\"pred_\"+y_rows] = rows_models[2*imod](surf_values[selblue])\n        traindf.loc[-selblue,\"pred_\"+y_rows] = rows_models[2*imod+1](surf_values[-selblue])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:59:04.392503Z","iopub.execute_input":"2025-05-30T14:59:04.392867Z","iopub.status.idle":"2025-05-30T14:59:04.428925Z","shell.execute_reply.started":"2025-05-30T14:59:04.392840Z","shell.execute_reply":"2025-05-30T14:59:04.427940Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Look at the errors in predicting the medians vs the surface average velocity","metadata":{}},{"cell_type":"code","source":"# These are just the deviations of the points from the model curves in the plots above.\nsurfAves = traindf[\"x_surfAve\"]\nsurfDiffs = traindf[\"x_surfDiff\"]\nsurf_L_values = surfAves - 0.5*surfDiffs\nsurf_R_values = surfAves + 0.5*surfDiffs\nfor y_rows in [\"09L\", \"09R\", \"1029\", \"3049\", \"5069\"]:\n    surf_values = surfAves\n    vel_axis_label = \"Ave Surface Velocity (m/s)\"\n    # Modify surf_values for the 09L,R data\n    if \"09L\" in y_rows:\n        surf_values = surf_L_values\n        vel_axis_label = \"L Surface Velocity (m/s)\"\n    if \"09R\" in y_rows:\n        surf_values = surf_R_values\n        vel_axis_label = \"R Surface Velocity (m/s)\"\n\n    plt.figure(figsize=(6,3))\n    plt.scatter(surf_values, traindf[\"y_\"+y_rows+\"Medi\"] - traindf[\"pred_\"+y_rows],\n                        c=traindf[\"diff_clr\"],s=2, alpha=0.15)\n    plt.ylim(-1500,1500)\n    plt.xlim(1400, 4200)\n    plt.title(\"Error in Predicted Medians for Rows \"+y_rows)\n    plt.xlabel(vel_axis_label)\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:59:04.430022Z","iopub.execute_input":"2025-05-30T14:59:04.430352Z","iopub.status.idle":"2025-05-30T14:59:06.127137Z","shell.execute_reply.started":"2025-05-30T14:59:04.430326Z","shell.execute_reply":"2025-05-30T14:59:06.125894Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Calculate the average MAE based on the 5 predicted medians.","metadata":{}},{"cell_type":"code","source":"# Add MAE_pred to the dataframe for each sample.\nMAE_preds = []\nfor dfind in traindf.index:\n    # Read in the data for this sample\n    velocity, data, isample = get_train_sample(dfind, ftscale=False)\n    # Go through the 5 row regions and calculate MAE wrt their predicted medians\n    MAE_5medi = 0.0\n    MAE_5medi += 5.0*np.mean(np.abs(velocity[isample,0, 0:10 , 0:34+1  ] - \n                                    traindf.loc[dfind,\"pred_09L\"]))\n    MAE_5medi += 5.0*np.mean(np.abs(velocity[isample,0, 0:10 , 35:  ] - \n                                    traindf.loc[dfind,\"pred_09R\"]))\n    MAE_5medi += 20.0*np.mean(np.abs(velocity[isample,0, 10:29+1 , :  ] - \n                                     traindf.loc[dfind,\"pred_1029\"]))\n    MAE_5medi += 20.0*np.mean(np.abs(velocity[isample,0, 30:49+1 , :  ] - \n                                     traindf.loc[dfind,\"pred_3049\"]))\n    MAE_5medi += 20.0*np.mean(np.abs(velocity[isample,0, 50: , :  ] - \n                                     traindf.loc[dfind,\"pred_5069\"]))\n    MAE_preds.append(MAE_5medi / 70.0)\n\ntraindf[\"MAE_pred\"] = MAE_preds","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:59:06.128246Z","iopub.execute_input":"2025-05-30T14:59:06.128571Z","iopub.status.idle":"2025-05-30T14:59:31.744682Z","shell.execute_reply.started":"2025-05-30T14:59:06.128549Z","shell.execute_reply":"2025-05-30T14:59:31.743549Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Save the training dataframe with predictions, etc.\ntraindf.to_csv(\"traindf.csv\", header=True, index=False, float_format='%.2f')\n\ntraindf","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:59:31.745634Z","iopub.execute_input":"2025-05-30T14:59:31.745961Z","iopub.status.idle":"2025-05-30T14:59:32.089761Z","shell.execute_reply.started":"2025-05-30T14:59:31.745918Z","shell.execute_reply":"2025-05-30T14:59:32.088593Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Use the sample submission to get the test ids\nsubmis = pd.read_csv(\"/kaggle/input/waveform-inversion/sample_submission.csv\")\n\n# Create a df of just the test ids (with _y_0)\noiddf = submis.loc[0:4607260:70,[\"oid_ypos\"]].copy()\noiddf = oiddf.reset_index(drop=True)\n##oiddf","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T14:59:32.094225Z","iopub.execute_input":"2025-05-30T14:59:32.094551Z","iopub.status.idle":"2025-05-30T15:00:05.585716Z","shell.execute_reply.started":"2025-05-30T14:59:32.094529Z","shell.execute_reply":"2025-05-30T15:00:05.584554Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# For each sample, add the measured surface velocity average and the R-L difference\nave_vels = []\ndiff_vels = []\nfor indoid in oiddf.index:\n    testdata = np.load('/kaggle/input/waveform-inversion/test/' + \n                   oiddf.loc[indoid,\"oid_ypos\"][0:10]+'.npy')\n    velL, velave, velR = info_data(testdata, for_show=False)\n    ave_vels.append(velave)\n    diff_vels.append(velR-velL)\n\noiddf[\"x_surfAve\"] = ave_vels\noiddf[\"x_surfDiff\"] = diff_vels\n##oiddf","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T15:00:05.586874Z","iopub.execute_input":"2025-05-30T15:00:05.587174Z","iopub.status.idle":"2025-05-30T15:22:49.326947Z","shell.execute_reply.started":"2025-05-30T15:00:05.587152Z","shell.execute_reply":"2025-05-30T15:22:49.323097Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Add color-coding based on the surfDiff and surfAve values\noiddf[\"diff_clr\"] = 'red'\n# Set blue, same criteria as for the training data\nseldiff = oiddf[\"x_surfAve\"] > (1300.0 + 1200.0*np.log10(1+np.abs(oiddf[\"x_surfDiff\"])))\noiddf.loc[seldiff, \"diff_clr\"] = 'blue'\n\n# Add model-predicted columns to the test dataframe based on x_surfAve\n# Order is blue then red model for each rows range.\nsurf_values = oiddf[\"x_surfAve\"]\nsurf_diffs = oiddf[\"x_surfDiff\"]\nsurf_L_values = surf_values - 0.5*surf_diffs\nsurf_R_values = surf_values + 0.5*surf_diffs\nselblue = oiddf[\"diff_clr\"] == 'blue'\n# Blue and Red model for each row range\nfor imod, y_rows in enumerate([\"09L\", \"09R\", \"1029\", \"3049\", \"5069\"]):\n    if y_rows == \"09L\":\n        oiddf.loc[selblue,\"pred_\"+y_rows] = rows_models[2*imod](surf_L_values[selblue])\n        oiddf.loc[-selblue,\"pred_\"+y_rows] = rows_models[2*imod+1](surf_L_values[-selblue])\n    elif y_rows == \"09R\":\n        oiddf.loc[selblue,\"pred_\"+y_rows] = rows_models[2*imod](surf_R_values[selblue])\n        oiddf.loc[-selblue,\"pred_\"+y_rows] = rows_models[2*imod+1](surf_R_values[-selblue])\n    else:\n        oiddf.loc[selblue,\"pred_\"+y_rows] = rows_models[2*imod](surf_values[selblue])\n        oiddf.loc[-selblue,\"pred_\"+y_rows] = rows_models[2*imod+1](surf_values[-selblue])\n\noiddf","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T15:22:49.329473Z","iopub.execute_input":"2025-05-30T15:22:49.329879Z","iopub.status.idle":"2025-05-30T15:22:49.450560Z","shell.execute_reply.started":"2025-05-30T15:22:49.329849Z","shell.execute_reply":"2025-05-30T15:22:49.449207Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Save the Test Dataframe\n","metadata":{}},{"cell_type":"code","source":"\noiddf.to_csv(\"oiddf.csv\", header=True, index=False, float_format='%.2f')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T15:22:49.451584Z","iopub.execute_input":"2025-05-30T15:22:49.451854Z","iopub.status.idle":"2025-05-30T15:22:50.556178Z","shell.execute_reply.started":"2025-05-30T15:22:49.451834Z","shell.execute_reply":"2025-05-30T15:22:50.554661Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"\\nMedian Ave Surface Velocity: {:.2f}\".format(np.median(oiddf[\"x_surfAve\"])))\nprint(\"Average Ave Surface Velocity: {:.2f}\\n\".format(np.mean(oiddf[\"x_surfAve\"])))\n\nplt.figure(figsize=(8,4))\nplt.hist(oiddf[\"x_surfAve\"],bins=100)\nplt.title(\"Test: Histogram of the Average Surface Velocity\")\nplt.xlabel(\"Surface Velocity (m/s)\")\nplt.xlim(velocity_range)\nplt.savefig(\"test_hist_surface_velocity.png\")\nplt.show()\n\nplt.figure(figsize=(8,4))\ndiffs = oiddf[\"x_surfDiff\"]\nplt.hist(np.sign(diffs)*np.log10(np.abs(diffs) + 1.0), log=True, bins=100)\nplt.title(\"Test: Histogram of the R-L Velocity Difference\")\nplt.xlabel(\"Signed Log10[1+ R-L Surface Velocity Difference (m/s) ]\")\nplt.savefig(\"test_hist_velocity_difference.png\")\nplt.show()\n\nplt.figure(figsize=(8,5))\nplt.scatter( np.sign(diffs)*(np.log10(np.abs(diffs) + 1.0)), oiddf[\"x_surfAve\"],\n                             color=oiddf[\"diff_clr\"], s=2, alpha=0.15)\nplt.ylabel(\"Ave Surface Velocity (m/s)\")\nplt.xlabel(\"Signed Log10[1+ R-L Velocity Diff (m/s) ]\")\nplt.title(\"Test: Average Surface Velocity vs R-L Velocity Difference\")\nplt.ylim(velocity_range)\nplt.savefig(\"test_scatter_velocity_vs_difference.png\")\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T15:22:50.557579Z","iopub.execute_input":"2025-05-30T15:22:50.557960Z","iopub.status.idle":"2025-05-30T15:22:54.995250Z","shell.execute_reply.started":"2025-05-30T15:22:50.557926Z","shell.execute_reply":"2025-05-30T15:22:54.994262Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Enter the predictions into the submis dataframe","metadata":{}},{"cell_type":"code","source":"# The predictions are 3 values for each of the 65818 test samples:\n# pred_09L value --> rows 0-9\n# pred_09R value --> rows 0-9\n# pred_1029 value --> rows 10-29\n# pred_3049 value --> rows 30-49\n# pred_5069 value --> rows 50-69\n\n# For each range of rows,\n# fill all 35 x_j values of the 65818 y_i values with the 65818 predicted values.\nall_xs = list(submis.columns[1:])\nleft_xs = list(submis.columns[1:17+1])\nright_xs = list(submis.columns[18:])\n\n# Loop over each set of y_i rows and set them equal to the corresponding predicted values\n\nlen_oiddf = len(oiddf)\n\n# Rows 0-9, with values adjusted for L (1,3,...,33) and R (35,37,...,69) halves.\nfill_values = (np.ones([17,len_oiddf]) * np.array(oiddf[\"pred_09L\"])).T\nfor iy in range(10):\n    rowsel = (submis.index % 70) == iy\n    submis.loc[rowsel, left_xs] = fill_values\nfill_values = (np.ones([18,len_oiddf]) * np.array(oiddf[\"pred_09R\"])).T\nfor iy in range(10):\n    rowsel = (submis.index % 70) == iy\n    submis.loc[rowsel, right_xs] = fill_values\n\n\n# Rows 10-29\nfill_values = (np.ones([35,len_oiddf]) * np.array(oiddf[\"pred_1029\"])).T\nfor iy in range(10,29+1):\n    rowsel = (submis.index % 70) == iy\n    submis.loc[rowsel, all_xs] = fill_values\n    \n# Rows 30-49\nfill_values = (np.ones([35,len_oiddf]) * np.array(oiddf[\"pred_3049\"])).T\nfor iy in range(30,49+1):\n    rowsel = (submis.index % 70) == iy\n    submis.loc[rowsel, all_xs] = fill_values\n\n# Rows 50-69\nfill_values = (np.ones([35,len_oiddf]) * np.array(oiddf[\"pred_5069\"])).T\nfor iy in range(50,69+1):\n    rowsel = (submis.index % 70) == iy\n    submis.loc[rowsel, all_xs] = fill_values\n\nsubmis","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T15:22:54.996429Z","iopub.execute_input":"2025-05-30T15:22:54.996823Z","iopub.status.idle":"2025-05-30T15:23:12.859398Z","shell.execute_reply.started":"2025-05-30T15:22:54.996801Z","shell.execute_reply":"2025-05-30T15:23:12.858297Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Generate the submission file\nsubmis.to_csv(\"submission.csv\", header=True, index=False, float_format='%.0f')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T15:23:12.861241Z","iopub.execute_input":"2025-05-30T15:23:12.861833Z","iopub.status.idle":"2025-05-30T15:27:18.410119Z","shell.execute_reply.started":"2025-05-30T15:23:12.861801Z","shell.execute_reply":"2025-05-30T15:27:18.409073Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Check it\n!ls -s submission.csv\n!tail -5 submission.csv","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-30T15:27:18.411138Z","iopub.execute_input":"2025-05-30T15:27:18.411466Z","iopub.status.idle":"2025-05-30T15:27:18.848919Z","shell.execute_reply.started":"2025-05-30T15:27:18.411420Z","shell.execute_reply":"2025-05-30T15:27:18.847858Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}