{"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":"## GWI 2025 - Velocity Patterns and a Training Dataframe\n\nThe [Yale/UNC-CH - Geophysical Waveform Inversion](https://www.kaggle.com/competitions/waveform-inversion/overview) competition is a bit intimidating: abstractly, the task is to learn how to transform 5 given 1000x70 input images into one corresponding 70x35 output image. There are large, pretrained ML models already created which seem to do very well on the clean, synthetic benchmark data. So, here I'm taking a simple tabular approach to make crude predictions.\n\nTo get an idea of the data, images of the velocity maps and seismic data are made with simple [plots based on the tutorial notebook](https://www.kaggle.com/code/hanchenwang114/waveform-inversion-kaggle-competition-tutorial). One novelty is to scale the seismic data by `f(t) = 1 + (t/a)^b` (with, e.g., a=200 ms, b=1.5) to help equalize the amplitudes vs time; this is similar to applying AGC when visualizing.\n\nA Training Dataframe is created and used as a place to add properties of the velocity maps (the targets) and the seismic data (the features). So far the only seismic features extracted are the surface velocities in the left (0-34) and right (34-69) halves combined to give Average and R-L Difference velocities for each sample. (The velocities are accurately measured using times based on polyfits to points around peaks. Quickly-reflected waves can affect the measurements.) \n\nScatter plots of the Average vs R-L Difference seismic velocities show an interesting pattern, similar for both the Train and Test data. The plots divide into samples with the surface velocity R-L difference relatively near 0 (blue points) and ones with asymmetric surface velocities (red points). The differences are probably due to the different vmap types: some are symmetric, some are not.\n\nThe targets to predict for each sample are the median training vmap values in coarse horizontal row ranges: 0-9 L&R, 10-29, 30-49, and 50-69. Correctly predicting all 5 medians for each sample would give MAE ~ 250. Simple polynomial models are fit to the vmap median values (y) versus the seismic average surface velocity (x) separately for the blue and red points. To approximate an MAE metric, the points are converted to median points in quantile ranges (by `xy_medians()`). These models are then used to make test predictions based on the measured test surface velocities.\n","metadata":{}},{"cell_type":"code","source":"# Notes:\n#   The 5 sources are x-located closest to: 0, 17, 34, 52, 69\n\n# Looked through the public notebooks, noted these:\n#   _This notebook has nice image-making code for velocity and data:\n#     https://www.kaggle.com/code/hanchenwang114/waveform-inversion-kaggle-competition-tutorial \n#   _\"Muting\" surface pulses is done in: https://www.kaggle.com/code/ozhiro/topmute-analysis \n#     Looks like \"AGC\" ~ adjusts gain within a (range of) row(s). Instead, I'll scale the data by the time.\n#   _This notebook mentions multiply-reflected waves, I wonder if they are seen/important :\n#     https://www.kaggle.com/code/nikita7364777/u-net-lb-413 \n\n# Some things to do next:\n# . . .","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-05-24T12:53:10.383651Z","iopub.execute_input":"2025-05-24T12:53:10.383930Z","iopub.status.idle":"2025-05-24T12:53:10.389948Z","shell.execute_reply.started":"2025-05-24T12:53:10.383905Z","shell.execute_reply":"2025-05-24T12:53:10.389049Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Usual things to use","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-24T12:53:10.391469Z","iopub.execute_input":"2025-05-24T12:53:10.391841Z","iopub.status.idle":"2025-05-24T12:53:12.510076Z","shell.execute_reply.started":"2025-05-24T12:53:10.391812Z","shell.execute_reply":"2025-05-24T12:53:12.508870Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Useful functions","metadata":{}},{"cell_type":"code","source":"# To show and process the Velocity Maps\n\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-24T13:54:44.335492Z","iopub.execute_input":"2025-05-24T13:54:44.335797Z","iopub.status.idle":"2025-05-24T13:54:44.353777Z","shell.execute_reply.started":"2025-05-24T13:54:44.335777Z","shell.execute_reply":"2025-05-24T13:54:44.352391Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# To show and process the Seismic Data\n\n# 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-24T12:53:20.789034Z","iopub.execute_input":"2025-05-24T12:53:20.789332Z","iopub.status.idle":"2025-05-24T12:53:20.806151Z","shell.execute_reply.started":"2025-05-24T12:53:20.789307Z","shell.execute_reply":"2025-05-24T12:53:20.805268Z"}},"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-24T12:53:27.361361Z","iopub.execute_input":"2025-05-24T12:53:27.361702Z","iopub.status.idle":"2025-05-24T12:53:27.370124Z","shell.execute_reply.started":"2025-05-24T12:53:27.361676Z","shell.execute_reply":"2025-05-24T12:53:27.368983Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Routine to read in the training data files given the training dataframe index value\n\n# 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-24T12:53:35.874055Z","iopub.execute_input":"2025-05-24T12:53:35.874356Z","iopub.status.idle":"2025-05-24T12:53:35.882357Z","shell.execute_reply.started":"2025-05-24T12:53:35.874333Z","shell.execute_reply":"2025-05-24T12:53:35.881397Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Convert many x,y points into a quartile-based set of x_median,y_median points.\n# 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-24T12:56:43.627791Z","iopub.execute_input":"2025-05-24T12:56:43.628129Z","iopub.status.idle":"2025-05-24T12:56:43.635860Z","shell.execute_reply.started":"2025-05-24T12:56:43.628105Z","shell.execute_reply":"2025-05-24T12:56:43.634712Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Make an Initial Training Dataframe","metadata":{}},{"cell_type":"code","source":"# There are 10,000 training samples on kaggle, organized as: 10 x 2 x 500 data-vel pairs\n##!ls /kaggle/input/waveform-inversion/train_samples/*\n\n# 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-24T12:56:51.542136Z","iopub.execute_input":"2025-05-24T12:56:51.542428Z","iopub.status.idle":"2025-05-24T12:56:51.563968Z","shell.execute_reply.started":"2025-05-24T12:56:51.542406Z","shell.execute_reply":"2025-05-24T12:56:51.562880Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"traindf","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-24T12:56:58.622203Z","iopub.execute_input":"2025-05-24T12:56:58.622584Z","iopub.status.idle":"2025-05-24T12:56:58.656360Z","shell.execute_reply.started":"2025-05-24T12:56:58.622498Z","shell.execute_reply":"2025-05-24T12:56:58.655539Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Look at a Training Velocity-Data Pair","metadata":{}},{"cell_type":"code","source":"# Select a dataframe index to look at\ndfind = 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-24T13:56:58.921181Z","iopub.execute_input":"2025-05-24T13:56:58.921532Z","iopub.status.idle":"2025-05-24T13:57:00.521801Z","shell.execute_reply.started":"2025-05-24T13:56:58.921505Z","shell.execute_reply":"2025-05-24T13:57:00.520583Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Look at the velocity map for the training sample\n# isample defined above\nplot_velocity(velocity, isample)\ninfo_velocity(velocity, isample)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-24T13:57:03.677286Z","iopub.execute_input":"2025-05-24T13:57:03.677981Z","iopub.status.idle":"2025-05-24T13:57:04.194392Z","shell.execute_reply.started":"2025-05-24T13:57:03.677955Z","shell.execute_reply":"2025-05-24T13:57:04.193508Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Look at the seismic data for the sample\n# 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-24T13:57:22.913106Z","iopub.execute_input":"2025-05-24T13:57:22.913382Z","iopub.status.idle":"2025-05-24T13:57:23.972884Z","shell.execute_reply.started":"2025-05-24T13:57:22.913362Z","shell.execute_reply":"2025-05-24T13:57:23.971950Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Look at a Test Data Sample","metadata":{}},{"cell_type":"code","source":"# The test data consists of 65818 data sets to be predicted\n## !ls /kaggle/input/waveform-inversion/test | wc\n# 65818   65818  987270\n\n##!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-21T19:19:39.837912Z","iopub.execute_input":"2025-05-21T19:19:39.838214Z","iopub.status.idle":"2025-05-21T19:19:40.930337Z","shell.execute_reply.started":"2025-05-21T19:19:39.838193Z","shell.execute_reply":"2025-05-21T19:19:40.929379Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Fill the Training Dataframe","metadata":{}},{"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-24T14:09:51.097494Z","iopub.execute_input":"2025-05-24T14:09:51.097916Z","iopub.status.idle":"2025-05-24T14:10:43.298559Z","shell.execute_reply.started":"2025-05-24T14:09:51.097891Z","shell.execute_reply":"2025-05-24T14:10:43.297735Z"}},"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-24T14:41:40.928814Z","iopub.execute_input":"2025-05-24T14:41:40.929255Z","iopub.status.idle":"2025-05-24T14:41:40.975732Z","shell.execute_reply.started":"2025-05-24T14:41:40.929229Z","shell.execute_reply":"2025-05-24T14:41:40.974718Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Save the Training Dataframe -- Save it further along after predictions are added.\n##traindf.to_csv(\"traindf.csv\", header=True, index=False, float_format='%.2f')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-24T12:59:27.045929Z","iopub.execute_input":"2025-05-24T12:59:27.046161Z","iopub.status.idle":"2025-05-24T12:59:27.226800Z","shell.execute_reply.started":"2025-05-24T12:59:27.046141Z","shell.execute_reply":"2025-05-24T12:59:27.225853Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Velocity Plots from the Training Dataframe","metadata":{}},{"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-24T12:59:27.233826Z","iopub.execute_input":"2025-05-24T12:59:27.234054Z","iopub.status.idle":"2025-05-24T12:59:28.406217Z","shell.execute_reply.started":"2025-05-24T12:59:27.234036Z","shell.execute_reply":"2025-05-24T12:59:28.405201Z"}},"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-24T12:59:28.407419Z","iopub.execute_input":"2025-05-24T12:59:28.407935Z","iopub.status.idle":"2025-05-24T12:59:29.185044Z","shell.execute_reply.started":"2025-05-24T12:59:28.407911Z","shell.execute_reply":"2025-05-24T12:59:29.183962Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Look into the y=0 Average and Difference\n\n# 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-24T12:59:37.142530Z","iopub.execute_input":"2025-05-24T12:59:37.142853Z","iopub.status.idle":"2025-05-24T12:59:39.231899Z","shell.execute_reply.started":"2025-05-24T12:59:37.142829Z","shell.execute_reply":"2025-05-24T12:59:39.230815Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Compare measured surface velocity with the y=0 average.\n# 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-24T13:01:18.424063Z","iopub.execute_input":"2025-05-24T13:01:18.424421Z","iopub.status.idle":"2025-05-24T13:01:19.102747Z","shell.execute_reply.started":"2025-05-24T13:01:18.424397Z","shell.execute_reply":"2025-05-24T13:01:19.101735Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Model the Row Ranges' Median Velocities","metadata":{}},{"cell_type":"code","source":"# Look at the velocities in row ranges vs the surface velocity and velocity difference.\n# 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-24T13:02:07.739570Z","iopub.execute_input":"2025-05-24T13:02:07.739889Z","iopub.status.idle":"2025-05-24T13:02:16.378633Z","shell.execute_reply.started":"2025-05-24T13:02:07.739868Z","shell.execute_reply":"2025-05-24T13:02:16.377729Z"}},"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-24T13:05:08.671135Z","iopub.execute_input":"2025-05-24T13:05:08.671495Z","iopub.status.idle":"2025-05-24T13:05:08.855819Z","shell.execute_reply.started":"2025-05-24T13:05:08.671470Z","shell.execute_reply":"2025-05-24T13:05:08.854543Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Mystery blues lines are when the region median equals the y=0 velocity","metadata":{}},{"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-24T14:43:37.866679Z","iopub.execute_input":"2025-05-24T14:43:37.867121Z","iopub.status.idle":"2025-05-24T14:43:39.447514Z","shell.execute_reply.started":"2025-05-24T14:43:37.867094Z","shell.execute_reply":"2025-05-24T14:43:39.446223Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Evaluate MAE for the Training Data","metadata":{}},{"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-24T15:42:44.878412Z","iopub.execute_input":"2025-05-24T15:42:44.879690Z","iopub.status.idle":"2025-05-24T15:42:45.159189Z","shell.execute_reply.started":"2025-05-24T15:42:44.879649Z","shell.execute_reply":"2025-05-24T15:42:45.157828Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 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-24T14:51:52.130167Z","iopub.execute_input":"2025-05-24T14:51:52.130584Z","iopub.status.idle":"2025-05-24T14:51:52.185276Z","shell.execute_reply.started":"2025-05-24T14:51:52.130555Z","shell.execute_reply":"2025-05-24T14:51:52.184400Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Look at the errors in predicting the medians vs the surface average velocity\n# 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-24T15:45:15.509409Z","iopub.execute_input":"2025-05-24T15:45:15.509828Z","iopub.status.idle":"2025-05-24T15:45:17.103971Z","shell.execute_reply.started":"2025-05-24T15:45:15.509797Z","shell.execute_reply":"2025-05-24T15:45:17.102592Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Calculate the average MAE based on the 5 predicted medians.\n# 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-24T15:16:36.911646Z","iopub.execute_input":"2025-05-24T15:16:36.912058Z","iopub.status.idle":"2025-05-24T15:17:12.325318Z","shell.execute_reply.started":"2025-05-24T15:16:36.912028Z","shell.execute_reply":"2025-05-24T15:17:12.324298Z"}},"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-24T15:47:47.576363Z","iopub.execute_input":"2025-05-24T15:47:47.576884Z","iopub.status.idle":"2025-05-24T15:47:47.925154Z","shell.execute_reply.started":"2025-05-24T15:47:47.576856Z","shell.execute_reply":"2025-05-24T15:47:47.924011Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"\\nOverall MAE of the predictions is: {:.2f}\\n\".format(np.mean(traindf[\"MAE_pred\"])))\n\nplt.figure(figsize=(6,3))\nplt.hist(traindf[\"MAE_pred\"],bins=100)\nplt.xlabel(\"MAE of the sample\")\nplt.title(\"Histogram of the MAEs of the 5-median Predictions\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-24T15:44:28.245641Z","iopub.execute_input":"2025-05-24T15:44:28.245987Z","iopub.status.idle":"2025-05-24T15:44:28.548774Z","shell.execute_reply.started":"2025-05-24T15:44:28.245965Z","shell.execute_reply":"2025-05-24T15:44:28.547605Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Create and Fill the Test Dataframe","metadata":{}},{"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},"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},"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},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Save the Test Dataframe\noiddf.to_csv(\"oiddf.csv\", header=True, index=False, float_format='%.2f')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Velocity Plots from the Test Dataframe","metadata":{}},{"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},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Submit a prediction","metadata":{}},{"cell_type":"code","source":"# Enter the predictions into the submis dataframe\n# 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},"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},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Check it\n##!ls -s submission.csv\n##!tail -5 submission.csv","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Submissions, fyi\n#  v4  1980.7    Predicted 1000 for all (so average is 2980.7)\n#       710.5    Submit the sample submission\n#  v5   605.5    Predicted 2781.0 for y=0-34 and 3781.0 for y=35-69\n# v11   -----    Made the Image: Add all 5 sources with surface pulses zeroed - not useful!\n# v13   551.4    Predicted 2481.0 for y=0-34 and 3481.0 for y=35-69. +300 not a good idea :)\n# v14   543.6    Predict 1966 (median surface vel) in rows 0-9, 2481 in rows 10-39, 3481 in rows 40-69\n# v15   495.1    Enter the actual surface velocity into rows 0-9, rows 10-39 = 2481, rows 40-69 = 3481\n# v19   493.5    Separate the rows 0-9 predictions into L and R values.\n# v20   489.0    Add x1.03 for rows 09 L,R, and use train average medians for 1039 and 4069.\n# v21   488.7    Use x1.017 for rows 09 L,R, same otherwise\n# v22   489.5    Use the pred_09 model values for rows 09 L,R, same otherwise.\n# v24   491.3    Use pred_09 but same for L,R (pred_09 is based on all xs)\n# v25   468.2    Filled the submission with the predicted values for each sample in each row range.\n#                Changed the seismic velocity calculation to be more accurate. Put back y0 vs surfAve plot.\n# v27   475.3    Has y09 model linear, but strange fit.\n# v28   484.1    Changed the y09 model to quadratic - strange fit also. Otherwise v27 with some cleaning up. \n# v29   473.1    Set y09 model to a constant.\n# v30   478.7    *** Use means in row regions instead of the medians ***\n# v33            Back to medians. Added a third predicted layer.\n#       442.1    Adjustments to surface velocity meas.: ignore peaks that are too far in time.\n# v34   440.3    Use surfDiff to adjust the L and R halves of the y09 predictions.\n# v35            Added some words.\n# v36   441.6    Better to plot medians directly vs surface velocity (not the ratios.)\n#                More adjustments to the surface velocity determination - it's tricky  :)\n#                Try to understand the \"blue lines\" in the medians vs surface velocity plots.\n# v37   435.3    Make separate models for the blue (R-L near 0) and red samples, not so different.\n# v38   435.1    Add surfAve > 3000 to the 'blue' samples; use degree 3 for model fits.\n# v39   434.9    Both fits are degree=3 (v38 was blue degree=3 but red degree=2)\n# v40   434.5    Redefined/Improved red/blue determination after comparing to the y0 R-L=0 truth.\n# v41   434.1    Reduced fitting range to surfAve < 4000. Changed plot ranges.\n# v43   --4.2    Fit 09L and 09R separately and adjust surfDiff factors added on.\n# v44   --6.0    v43 wasn't what I intended, try again. yeesh, v44 had double adjustments.\n# v45   434.3    Hmmm... Made slight tweak to 09L,R models...\n# v46   434.4    Still messing with the 09 stuff\n# v47   --4.4    Fit combined 09L,09R so they use same models, equivalent to adding flipped data.\n# v48   434.3    Fixed tiny error in v47.\n# v50   427.4    Used xy_medians() to transform x,y before the model fitting ~ as if using MAE. \n# v51   427.8    Increased fitting degree to 5; surfAve clipped at 4100; other little changes.\n# v52            Calculated MAE values: if all 5 medians in each sample were known MAE = 249.07 .\n#                Using each sample's 5 predicted medians gives MAE = 432.84 <-- Similar to LB score.","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}